譜體積方法單元分割與問題單元標(biāo)記方法的改進(jìn)研究

柯耀祖， 張 兵

(合肥工業(yè)大學(xué) 機(jī)械工程學(xué)院,合肥 230000)

1 引 言

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

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

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

2 譜體積方法

積分形式的Euler方程為

(1)

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

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

(2)

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

(3)

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

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

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

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

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

(5)

式中多項式Pi的系數(shù)由控制體的守恒變量均值表達(dá)。此高階重構(gòu)多項式用于計算SVi內(nèi)任意位置處的守恒變量值。

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

(6)

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

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

(7)

式中J=(k+1)/2為每個面上的積分點(diǎn)數(shù)目，?r q為高斯積分權(quán)重，(xr q,yr q)為積分點(diǎn)坐標(biāo)。

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

3 分區(qū)方法和限制器

3.1 譜體積單元的分區(qū)

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

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

(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內(nèi)部邊界面。

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

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單元分區(qū)和通量點(diǎn)的位置

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

圖2 二維S V單元分區(qū)

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

3.2 限制器函數(shù)

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

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

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

(11)

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

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

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

(12)

式中下標(biāo)rq表示積分點(diǎn)。

(3) 上述的標(biāo)記規(guī)則過于嚴(yán)格，可能得到很多非問題單元。因此，本文通過改進(jìn)后的minmod限制器函數(shù)判斷最終需要限制的問題單元。

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

l=lxi+lyj

(13)

則CVi ,j和鄰居控制體這個方向上的一階和二階導(dǎo)數(shù)定義為

Ql ,i=Qx ,ilx+Qy ,ily

Ql ,neigh=Qx,neighlx+Qy ,neighly

(14)

式中下標(biāo)neigh表示鄰居控制體，Qx,Qy,Qx x,Qx y和Qy y是CVi ,j的一階和二階偏導(dǎo)數(shù)，可通過求解二次方最小二乘重構(gòu)[22]得到。

(15)

式中r為SV的體心位置。

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

(16)

式中下標(biāo)k為控制體的面編號，β取1.5。

③ 循環(huán)CVi ,j所有面，計算CVi ,j的標(biāo)量限制器最小值，即

(17)

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

(18)

4 算例驗證

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

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

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

圖3 不同密度的計算網(wǎng)格

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

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

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

圖5 不同分區(qū)方法的壓力系數(shù)計算結(jié)果

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

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

4.2 NACA0012跨音速流

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

圖7 NACA0012機(jī)翼流場網(wǎng)格

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

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

圖9 NACA0012機(jī)翼的跨音速流

4.3 計算效率分析

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

表1 SVM和FVM計算時間對比

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

5 結(jié) 論

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

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

(3) 采用二維可壓縮Euler方程典型算例驗證了上述方法的有效性。結(jié)果表明，采用較少的網(wǎng)格即可獲得與數(shù)倍密網(wǎng)格下傳統(tǒng)有限體積法相當(dāng)?shù)挠嬎憔?，并具有較好的計算效率。