譜體積方法單元分割與問題單元標記方法的改進研究

柯耀祖， 張 兵

(合肥工業大學 機械工程學院,合肥 230000)

1 引 言

基于TVD性質發展的傳統有限體積法是計算流體力學中的重要方法，廣泛應用于各大商業CFD軟件和研究代碼中，但其計算精度與網格劃分密切相關,往往需要極大的網格數量才能得到網格收斂解。盡管網格自適應方法一定程度上可解決上述問題，但針對激波和接觸間斷處于運動狀態的非定常計算問題，不但會導致計算量劇增，還存在難以解決的并行負載平衡問題[1]。

譜體積方法(SVM)是一種從本質上解決求解器對網格依賴性的高精度方法，其基本原理是將一個網格單元(簡稱譜體積元SV)按精度需要分割成多個內部控制體(簡稱CV)，并假設譜體積內變量的空間分布為控制體平均值的多項式函數，進而采用傳統有限體積方法對方程進行空間和時間離散，并求解得到具有高階精度的解。由于譜體積元內部精度可按需要提高至任意階，故該方法可實現在任意密度網格上獲得網格收斂解[2]。譜體積方法由Wang[3]首次提出并應用于一維Euler方程求解，隨后推廣到二維及三維問題[4-7]。為了進一步滿足不同問題求解需求，Liang等[8]對譜體積方法進行擴展，如N-S方程的多重網格方法。Breviglieri[9]通過一種基于雙曲守恒定律的空間離散方法，實現對流場問題的高精度求解。高精度格式是現在CFD方法的研究熱點，侯冶[10]研究了一維SVM方法，并與有限差分法對比驗證SVM的高效性。鄭華盛等[11]基于譜體積方法思想，通過篩選出適當的單元組成模板，利用WENO重構方法和有限體積法，構造出適用于非結構網格二階和三階精度的有限體積WENO格式。針對多物理耦合計算耗費過大，劉娜等[12]設計出適合多介質流動模擬的譜體積方法。值得說明的是，由于N-S方程的特征決定了其解的形式是非光滑的，特別是對于存在激波的問題，譜體積方法的多項式假設必然會導致解的震蕩，需要采用合適的控制體分割方法、限制器函數及問題單元標記方法進行處理，三者決定了SV方法的穩定性和精確性。van den Abeele等[13]提出通過傅里葉分析檢測分區質量和色散誤差，來確定非結構網格單元的穩定性；Lamouroux等[14]提出一種基于空間加權投影的高階限制器；Hadadian Nejad Yousefi等[15]則提出了一種基于WENO方案的切比雪夫SV方法，用于求解一維和二維的守恒定律方程。

由于切比雪夫多項式的正交性、奇偶性和有界性，其在連續函數逼近問題中得到廣泛應用，如能最大限度地降低龍格現象，可以提供多項式在連續函數的最佳一致逼近等，因此本文提出一種采用切比雪夫多項式確定非結構網格單元的內部控制體劃分方法。同時,為防止不連續區域在求解過程中數值發散，在minmod TVB限制器[16]基礎上，采用參數自由的廣義限制器，通過標記單元和限制標記問題單元這兩步驟，以保證平滑區域的精度和避免不連續區域的解振蕩。最后，通過二維Euler方程的典型算例進行驗證。

2 譜體積方法

積分形式的Euler方程為

(1)

式中t為時間，Ω為計算域，Q為守恒變量，F為無粘通量。

將求解域離散成N個不重疊譜體積單元，即

(2)

根據有限體積假設，對于任意譜體積單元，式(1)可改寫為

(3)

對二維問題，為了達到k階精度，需要將譜體積元SVi分成m=k(k+1)/2個控制體CV(3.1節)，任意控制體CVi ,j的守恒變量平均值可表示為

(i=1,…,N；j=1,…,m)(4)

式中Vi ,j為控制體CVi ,j的體積。

當SVi內所有控制體的平均守恒變量已知時，可以重構出一個(k-1)次多項式Pi(x,y)，以k階精度在SVi內逼近Q(x,y)。

Q(x,y)=Pi(x,y)+O(hk -1)

(5)

式中多項式Pi的系數由控制體的守恒變量均值表達。此高階重構多項式用于計算SVi內任意位置處的守恒變量值。

對控制體CVi ,j，方程(1)可寫為

(6)

采用k階高斯數值積分方法對式(3)進行面積分可得

[Q(xr q,yr q)]Ar+O(Arhk)

(7)

式中J=(k+1)/2為每個面上的積分點數目，?r q為高斯積分權重，(xr q,yr q)為積分點坐標。

由于在每個SVi內部守恒變量分布滿足多項式假設，故在SVi內的控制體交界面上的通量F直接采用解析表達式計算。而在譜體積之間的交界上，由于變量的重構結果并不相同，故需要采用迎風格式計算，本文采用基于Harten熵修正方法[17]的Roe格式。

3 分區方法和限制器

3.1 譜體積單元的分區

插值多項式Pi(x,y)的性能是由形函數Li ,j(x,y)決定的，而Li ,j(x,y)又是通過CVi ,j單元的劃分計算得到。一般來說，可以采用任意分區方法，但是不恰當的分區方法可能會導致結果振蕩，降低計算精度[18,19]。因此，一個良好的分區方法是保證高精度數值計算的前提。

為提高模擬過程的穩定性和精確性，van den Abeele[13]和Wang[5]詳細探討了二維譜體積單元的劃分方法。在此基礎上，本文提出使用基于切比雪夫多項式[20]的分區方法，以保證計算結果的收斂性。以一維SV方法為例，對任意一個譜體積單元，建立單元局部坐標系x∈[-1,1]，其分區過程如下。

(1) 將切比雪夫多項式的根作為控制體的形心,

(k=0,1,…,n-1)(8)

式中xk為切比雪夫多項式的根，n為計算精度。

(2) 通過xk確定所有CVi ,j的邊界，得到所有CVi ,j邊界面的位置為

f1=-1，fn +1=1,fk +1=f1+2|xk|

(k=1,…,n)(9)

式中f1和fn為SVi外部邊界面，fk為SVi內部邊界面。

(3) 如果階次n是偶數階，則先確定SVi靠近外側的內部邊界面，即

f1=-1，fn + 1=1,fk + 1=2xk-fk

fn + 1 - k=2xn + 1 - k-fn + 2 - k

(k=1,…,n/2)(10)

通過上述的三個步驟確定最終的通量面位置，如 圖1 給出的就是一維CVi ,j邊界面的分布。

圖1 一維S V單元分區和通量點的位置

對于二維非結構網格的任意譜體積元，采用面積坐標定義分割點位置，將切比雪夫多項式的根作為SVi三個邊界面上數值點xk，然后根據xk求出面通量位置，圖2為2階和3階SVi單元的分割結果，其中，xk為切比雪夫多項式的根，而fk則為控制體邊界在譜體積元面上的位置。

圖2 二維S V單元分區

對于四階及四階以上的偶數階精度，和一維類似，先確定外側控制體邊界的位置，再逐步確定內部控制體的邊界。

3.2 限制器函數

對于Euler方程，如果解存在間斷或不連續，則需要在通量積分點處進行限制，以保證數值穩定性和收斂性。本文將限制器計算過程分為兩部分，首先找到并標記出需要限制的問題單元，然后對該單元上的所有積分點處的變量進行限制。在后一過程中，采用泰勒級數展開重構譜體積單元平均化導數[21]，然后從最高階導數到最低階導數限制問題單元。如果最高階導數不受限制，則取消該單元的標記。這種限制器能夠在保證平滑區域精度的前提下，有效地抑制不連續區域附近的振蕩。

與已有研究方法不同，Wang等[5，6]是直接通過任意譜體積的相鄰譜體積單元來判斷該譜體積單元是否需要限制，而本文使用的單元標記都是在控制體單元內進行的，但是限制卻在譜體積單元中展開。這是由于一旦控制體單元標記為問題單元，那么當前控制體所屬的譜體積元的插值多項式不再繼續適用。本文使用的限制器是基于ENO和WENO方法改進的TVB標記方法，以三階精度為例，將限制和重構的過程分為以下幾個階段。

(1) 對于任意一個控制體單元，采用Q表示其任一主變量(ρ，u，v，p)，首先計算出該單元與其共面的所有相鄰單元中的最大值和最小值主變量。

(11)

(2) 如果當前控制體單元中給的任一積分點處的變量滿足以下條件，則初步認為這個單元為問題單元。

Qi(xr q,yr q)>1.001Qmax,i

Qi(xr q,yr q)<0.999Qmin,i

(12)

式中下標rq表示積分點。

(3) 上述的標記規則過于嚴格，可能得到很多非問題單元。因此，本文通過改進后的minmod限制器函數判斷最終需要限制的問題單元。

① 定義從當前CVi ,j形心到與其共面的相鄰控制體形心的方向矢量為

l=lxi+lyj

(13)

則CVi ,j和鄰居控制體這個方向上的一階和二階導數定義為

Ql ,i=Qx ,ilx+Qy ,ily

Ql ,neigh=Qx,neighlx+Qy ,neighly

(14)

式中下標neigh表示鄰居控制體，Qx,Qy,Qx x,Qx y和Qy y是CVi ,j的一階和二階偏導數，可通過求解二次方最小二乘重構[22]得到。

(15)

式中r為SV的體心位置。

② 對該面上積分點處的物理量，采用式(16)進行初步限制。

(16)

式中下標k為控制體的面編號，β取1.5。

③ 循環CVi ,j所有面，計算CVi ,j的標量限制器最小值，即

(17)

Qx y(ΔxΔy-Ix y)]

(18)

4 算例驗證

通過兩個算例結果驗證本文的分區方法和限制器性能，所用的時間積分全部為隱式LU-SGS方法。

4.1 15°壓縮拐角超聲速流動

第一個算例是NASA算例庫的15°壓縮拐角超聲速流動算例[23]，其中來流邊界位于左側x=-0.3左，右側出口邊界位于x=0.5，計算域高度取0.75，馬來流馬赫數為2.0,采用疏密不同的計算網格進行求解，以便對比不同方法的計算性能。其中，圖3(a)的網格用于譜體積法和常規有限體積法(簡稱FVM)求解,單元數量為1137。圖3(c,d)僅使用FVM求解，單元數量分別為4633和18731。圖3(b)展示的則是三階精度下,劃分控制體后的網格。

圖3 不同密度的計算網格

如圖4為物面壓力系數分布結果曲線，可以看出，粗網格下SVM和密網格下FVM的求解結果吻合較好，計算得到的激波比粗網格FVM結果更陡峭。值得一提的是，SVM僅使用1137個譜體積元(6822個控制體元)，而FVM則使用18731個網格。

圖4 SVM和FVM的壓力系數分布曲線對比

圖5展示的是在相同網格數條件下，采用改進后的控制體分割方法和使用Wang等[5]提出的分區方法計算結果與密網格下FVM計算結果的對比，可以看出，本文方法的激波更陡。

圖5 不同分區方法的壓力系數計算結果

圖6(a)為三階SV方法計算得到的楔形板壓力分布情況。采用改進后的限制器，只對激波區域的譜體積單元進行標記和限制，如圖6(b)所示，激波處的黑色單元為標記的問題單元，其出現的位置全部位于激波位置處，故可有效識別不連續區域。

圖6 15°楔形板的超音速流，3階

4.2 NACA0012跨音速流

第二個算例是Yang等[16]的NACA0012跨音速流算例，其中翼型展長為1，來流馬赫數為0.8，迎角1.25°，分別采用三種不同密度網格進行計算，網格單元數量分別為5128,9814和20880。采用如圖7(a)所示的粗網格用于SVM和FVM求解,圖7(c，d)僅使用FVM計算,圖7(b)為二階譜體積元內部的控制體單元劃分情況。圖8為翼型表面的壓力系數分布,可以明顯看出，粗網格下的FVM結果激波過于平滑，說明網格密度對FVM求解的精度影響較大。粗網格下SVM方法得到的壓力分布與密網格下FVM結果吻合較好，說明SVM方法大大降低了網格需求。另一方面，通過對相同網格數下不同分割方法計算結果的對比可以發現，相比于Wang[5]的分割方法，采用切比雪夫多項式劃分譜體積元的方法得到的激波更陡峭。

圖7 NACA0012機翼流場網格

圖8 SVM和FVM以及不同分割方法的的壓力系數分布

為清楚觀察標記的問題單元，只放大機翼周圍的單元，圖9(a)為使用三階SV方法求解獲得的馬赫分布圖，圖9(b)給出標記問題單元的位置在上下表面激波處。測試結果表明，本文采用的分區方法和限制器可以有效標記激波不連續區域，并對其進行合理限制。

圖9 NACA0012機翼的跨音速流

4.3 計算效率分析

盡管SV方法在一定程度上減少了求解器對網格的依賴性，降低了CFD分析時間，但其基本原理決定了SVM在求解過程中需占據更多的計算資源。表1為相同網格數下三階SVM和FVM的計算效率對比，其中第一行為FVM計算時間，而第二行為SVM的計算時間。所有程序由C++語言編寫，并在i7-8750H，2.20 GHz處理器上完成測試。

表1 SVM和FVM計算時間對比

由表1可知，在相同網格數條件下，SV方法需要的計算時間大于傳統的二階FVM；對本文算例而言，在滿足相同精度條件下，粗網格下的SVM計算時間大致與4倍網格密度下的FVM相當。

5 結 論

(1) 建立了基于切比雪夫多項式方法的譜體積分割方法，典型算例表明，該方法具有更好的計算精度。

(2) 發展了基于WENO的變量限制方法，以及綜合考慮譜體積和控制體的問題單元標記方法，以有效識別不連續區域。

(3) 采用二維可壓縮Euler方程典型算例驗證了上述方法的有效性。結果表明，采用較少的網格即可獲得與數倍密網格下傳統有限體積法相當的計算精度，并具有較好的計算效率。