適用于混合網(wǎng)格的約束最小二乘重構(gòu)方法

康忠良， 閆超

北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院， 北京 100191

適用于混合網(wǎng)格的約束最小二乘重構(gòu)方法

康忠良， 閆超*

北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院， 北京 100191

為強化邊界條件對計算流體力學(xué)(CFD)計算中流場重構(gòu)過程的影響，基于三維任意多面體混合網(wǎng)格，對一種約束最小二乘法進(jìn)行了研究。通過對不同邊界類型的透明化處理，統(tǒng)一了約束方程組的形式，從而簡化了該方法的應(yīng)用。為了充分吸收消去法和加權(quán)法各自求解約束方程組的優(yōu)點，提出了混合采用兩種方法的思路。針對層流平板研究了加權(quán)法中加權(quán)系數(shù)如何取值問題，數(shù)值實驗表明加權(quán)系數(shù)取到5基本已經(jīng)足夠大。最后通過亞聲速層流平板和跨聲速湍流ONERA-M6算例對比了混合法、加權(quán)法與原始最小二乘法，計算結(jié)果表明：約束最小二乘法相比原始最小二乘法的流場重構(gòu)更加準(zhǔn)確，特別是對邊界質(zhì)量較差網(wǎng)格的計算優(yōu)勢更為突出；混合消去加權(quán)法相比單獨采用加權(quán)法的計算結(jié)果也有所改善。

混合網(wǎng)格； 任意多面體； 重構(gòu)方法； 約束最小二乘法； 消去法； 加權(quán)法

在計算流體力學(xué)(CFD)的應(yīng)用中，采用混合網(wǎng)格可以充分利用非結(jié)構(gòu)網(wǎng)格的網(wǎng)格生成優(yōu)勢和類結(jié)構(gòu)網(wǎng)格的計算優(yōu)勢，因此是一種先進(jìn)有效的方法?；旌暇W(wǎng)格的單元類型可以是常用的四面體、金字塔、三棱柱或六面體，也可以是這些單元之間以及與其他任意多面體的混合。

由于混合網(wǎng)格單元類型的任意性及拓?fù)涞碾S機性，所以必須采用無序的非結(jié)構(gòu)方式來存儲，從而給計算帶來了很大困難。在空間離散方面，雖然關(guān)于非結(jié)構(gòu)高階格式的研究已有了一定的進(jìn)展[1]，但由于計算效率和穩(wěn)定性等多方面的原因，目前混合網(wǎng)格中廣泛使用的仍然是二階空間離散格式，甚至一般二階精度都難以保證，其主要困難在于流場重構(gòu)方面[2]。對于普遍使用的Barth and Jespersen[3]提出的線性重構(gòu)方法，其中梯度計算是問題的關(guān)鍵[4-6]。

目前常用的梯度計算方法主要有格林-高斯法和最小二乘法。其中格林-高斯法在不同單元類型銜接處的梯度計算誤差可能會很大，不太適合于混合網(wǎng)格應(yīng)用[4]。而最小二乘法對于任意形狀網(wǎng)格的梯度計算都能達(dá)到一階精度，可滿足線性重構(gòu)的要求，從而保證空間離散達(dá)到二階精度，因此非常適合于混合網(wǎng)格計算。

本文主要針對一種新型的約束最小二乘法展開研究。Ollivier-Gooch and van Altena[7]首先在二維網(wǎng)格上采用消去法研究了約束最小二乘問題，Haselbacher[8-9]后來基于三維混合網(wǎng)格對比分析了消去法和加權(quán)法的優(yōu)劣，但其算法構(gòu)造對不同邊界類型不具有通用性，且沒有給出加權(quán)系數(shù)的取值影響。為方便于實際應(yīng)用，本文基于三維任意多面體混合網(wǎng)格對約束最小二乘法做了簡化和改進(jìn)，混合采用加權(quán)法和消去法進(jìn)行求解，并選用有效的算例進(jìn)行了驗證和分析。

1 格心有限體積法

本文算法構(gòu)造基于任意多面體混合網(wǎng)格單元。算法不考慮具體的網(wǎng)格拓?fù)?，對不同的單元類型統(tǒng)一處理(即對網(wǎng)格是透明的)，因此對任意的封閉網(wǎng)格都是通用的?？臻g離散采用格心有限體積法，控制體取網(wǎng)格單元，控制面取單元表面，流場變量存儲在單元中心。

1.1 控制方程

舍去源項的三維非定常可壓縮Navier-Stokes方程組在直角坐標(biāo)系下的守恒積分形式可表示為

(1)

式中：Ω為控制體；?Ω為控制面；Q為守恒變矢量；Fc為對流通量；Fv為黏性通量;S為面積。各項的具體描述詳見文獻(xiàn)[10]。

1.2 數(shù)值離散方法

對流通量離散采用Roe[11]格式近似求解Riemann問題。假設(shè)控制面IJ的左右單元分別為I和J,則對流通量可表示為

|ARoe|(QR-QL)]

(2)

式中：|ARoe|為Roe平均矩陣；QL和QR分別為控制面的左右狀態(tài)。如果QL和QR分別取左右單元中心的平均值，則空間離散只有一階精度，為了達(dá)到二階精度，就需要對其進(jìn)行重構(gòu)。

Barth and Jespersen[3]提出的線性重構(gòu)方法假定解在控制體內(nèi)呈線性分布，面左右兩側(cè)的值可表示為

{

(3)

線性重構(gòu)方法在規(guī)整網(wǎng)格上可以使空間離散達(dá)到二階精度[13]。假如梯度是準(zhǔn)確值，該方法在任意網(wǎng)格上都可以準(zhǔn)確重構(gòu)，因此梯度計算是重構(gòu)過程的關(guān)鍵。

另外，考慮黏性通量的橢圓形特征，對其采用中心格式離散，控制面上的值和梯度根據(jù)左右單元求取。時間離散采用LU-SGS(Lower-Upper Symmetric Gauss-Seidel)隱式格式[10]。

2 約束最小二乘法

由1.2節(jié)可知，梯度計算是線性重構(gòu)的關(guān)鍵，最小二乘法是目前常用的一種梯度計算方法，且非常適合于混合網(wǎng)格計算，下面將基于原始最小二乘法給出一種新型的約束最小二乘法。

2.1 算法構(gòu)造

采用最小二乘法計算梯度最早是由Barth[14]提出的，它是基于變量在體心值之間的一階泰勒級數(shù)近似，對于單元I和J，可表示為

(4)

將式(4)應(yīng)用到單元I的所有模板(一般取其所有鄰居單元)，便會得到一個線性方程組：

(5)

式中:Δ(*)IJ=(*)J-(*)I；?(*)U=?U/?(*)；q為模板個數(shù)；θ為模板加權(quán)系數(shù)，一般取為1，也可以根據(jù)網(wǎng)格幾何量或流場解取值，比如取距離倒數(shù)加權(quán)可表示為

(6)

為便于表示，式(6)可簡寫為

Ax=b

(7)

式中：A為q×n維系數(shù)矩陣，這里n=3，顯然對任意網(wǎng)格均能滿足q≥n，于是方程式(5)可采用最小二乘法求解[15]。

針對上述最小二乘法，Mavriplis[6]研究指出：加權(quán)系數(shù)θ采用距離倒數(shù)對于在彎曲壁面拉伸網(wǎng)格中保證計算精度比較重要。但是距離倒數(shù)加權(quán)的一大弊病在于它會導(dǎo)致計算穩(wěn)定性變差，因此本文中θ取為1。

實際上，可以認(rèn)為θ是對各個模板影響梯度程度的一種約束，采用距離倒數(shù)加權(quán)即是強化距離較近的模板的影響。與此類似，如果選取的模板中包括邊界，就考慮采取強化邊界條件對流場重構(gòu)的影響，也就是對每個邊界控制面施加約束，保證方程的最小二乘解優(yōu)先滿足邊界，該方法稱為約束最小二乘法，顯然這是一個十分合理的思路。

下面構(gòu)造約束最小二乘方程。Haselbacher[8-9]將不同邊界類型區(qū)分對待來構(gòu)造約束方程，導(dǎo)致了算法實現(xiàn)非常繁瑣，不利于該方法的實際應(yīng)用。考慮到算法對邊界條件的透明性，本文在每一迭代步采取將不同邊界類型均轉(zhuǎn)化為依賴流場變量來定義(比如對于無滑移壁面，速度均取為0，壓力取為首層網(wǎng)格格心值，絕熱壁處密度取為首層網(wǎng)格格心值，等溫壁處密度根據(jù)給定溫度和首層網(wǎng)格壓力由狀態(tài)方程求取)，這一處理大大簡化了約束最小二乘法的實施。對于式(4)，如果認(rèn)為UI是未知量，在邊界施加線性等式約束，構(gòu)造的約束方程組可統(tǒng)一表示為

(8)

式中：ΔxF，I1、ΔyF，I1和ΔzF，I1分別為邊界面模板rI1的x、y和z分量；前p個方程根據(jù)邊界模板來確定，求解過程應(yīng)該精確滿足這p個方程；后q個方程根據(jù)其他模板來確定，可近似滿足。

2.2 方程求解

式(8)可采用消去法或加權(quán)法求解。消去法是通過某種方法將前p個方程消去然后求解。比如對于圖1，要計算單元0處的梯度，選取模板為a-1-2-3-4，則構(gòu)造的方程組可寫為

(9)

對式(9)采用消去法，可降維變形為

(10)

可見式(10)與式(5)類似(注意各下標(biāo)的區(qū)別)，它可以認(rèn)為是將重構(gòu)位置從0移到了a，即優(yōu)先滿足了邊界。但是該方法只有當(dāng)模板中存在一個邊界時，才能通過消去降維到類似于式(5)的形式。

加權(quán)法是對式(8)前p個方程分別施加一個相對較大的權(quán)系數(shù)然后求解。同樣對于圖1，要計算單元1處的梯度，選取模板為b-c-0-2-3，則采用加權(quán)法構(gòu)造的方程組最終可寫為

(11)

式中：w為大于1的加權(quán)法權(quán)系數(shù)，在保證系數(shù)矩陣不病態(tài)的前提下，可對其適當(dāng)加大來強化邊界對梯度計算的貢獻(xiàn)。

圖1 約束最小二乘法構(gòu)造舉例圖(a、b和c為邊界面)Fig.1 Illustration of constrained least-squares reconstruction (a, b and c indicate boundary faces)

消去法的優(yōu)點在于可以使約束方程組精確滿足邊界，缺點是它限制了約束邊界個數(shù)，且對多個邊界的處理繁瑣并會導(dǎo)致計算魯棒性變差；加權(quán)法雖然對約束個數(shù)沒有限制，但其計算的梯度只能近似滿足邊界，而且加權(quán)系數(shù)過大又容易導(dǎo)致系數(shù)矩陣呈現(xiàn)病態(tài)[9]。本文實現(xiàn)中，考慮充分吸收兩種方法的優(yōu)點，采用混合消去和加權(quán)法：當(dāng)模板中只存在一個邊界時，采用消去法求解；當(dāng)存在多個邊界時，采用加權(quán)法求解。

另外，由于系數(shù)矩陣在某些條件下可能會呈現(xiàn)病態(tài)，從而導(dǎo)致計算的不穩(wěn)定，因此本文針對病態(tài)情況考慮選取一個優(yōu)化的模板。具體的模板構(gòu)造方法為：首先選取計算單元的所有鄰居單元組成其初始模板，對于邊界曲率較大處及多個邊界銜接處附近，適當(dāng)擴(kuò)充邊界模板；然后判斷構(gòu)造的系數(shù)矩陣是否呈病態(tài)，病態(tài)判斷方法詳見文獻(xiàn)[16]，對于病態(tài)模板，利用已選入模板的鄰居單元逐層擴(kuò)大模板，直至系數(shù)矩陣轉(zhuǎn)為良態(tài)。

3 算例驗證及分析

3.1 計算精度驗證

因為層流平板存在Blasius精確解，可以嚴(yán)格地測試速度和溫度分布的計算精度，因此是驗證本文算法的一個很有效算例。

選用3套網(wǎng)格G1、G2和G3，其區(qū)別在于邊界層內(nèi)網(wǎng)格數(shù)分別取為10、20和40個。計算網(wǎng)格G1如圖2所示，板長取一個單位長度，遠(yuǎn)場取到2倍板長，平板前后均多設(shè)置一段對稱面用于消除遠(yuǎn)場的影響。邊界層附近采用四邊形，以外區(qū)域采用三角形，由于程序是三維算法，因此在展向拉伸出一層網(wǎng)格，這樣就相當(dāng)于壁面附近網(wǎng)格單元是六面體，其他區(qū)域是三棱柱。

來流馬赫數(shù)Ma∞=0.3，以板長為特征長度雷諾數(shù)Re=1.0×105，壁面取為絕熱壁。

圖2 平板計算網(wǎng)格G1Fig.2 Computational grid G1 of flat plate

首先基于網(wǎng)格G1，通過數(shù)值實驗研究加權(quán)法中加權(quán)系數(shù)w的取值問題。暫對所有邊界模板的w取相同值，即取w=1，3，5，10，100進(jìn)行計算。由于摩擦阻力系數(shù)Cf的計算主要依賴于壁面處的梯度值和溫度分布，因此它對檢驗約束最小二乘法的計算精度十分有效。圖3為采用加權(quán)法計算中不同的w值時的摩擦阻力系數(shù)的分布，圖中x為距平板前緣的距離。從圖3可以看出，對于網(wǎng)格G1，加權(quán)系數(shù)取3左右影響計算結(jié)果變化較快，取到5基本已經(jīng)足夠大，再增大到10和100對結(jié)果影響不大，因此本文以下計算，加權(quán)系數(shù)均取為5。

圖3 不同加權(quán)系數(shù)時的摩擦阻力系數(shù)分布(網(wǎng)格G1)Fig.3 Skin-friction factor distribution for different weighting coefficients (grid G1)

然后對3套網(wǎng)格分別采用混合消去加權(quán)法(Mixed)和加權(quán)法(Weighting)，并與原始最小二乘法(Least-squares)進(jìn)行了對比，如圖4所示。

圖4 平板摩擦阻力系數(shù)分布Fig.4 Skin-friction factor distribution along flat plate

由圖4可見，對于較稀疏網(wǎng)格G1，約束最小二乘法相對原始最小二乘法的計算精度提高十分明顯，且混合采用消去法和加權(quán)法相比單獨加權(quán)法的精度也有所改善。隨著網(wǎng)格的加密，3種方法的計算結(jié)果都趨于理論解，約束最小二乘法的優(yōu)勢也變得不太明顯，分析主要是因為在邊界層特別是線性底層內(nèi)的網(wǎng)格量足夠多時，強化邊界影響梯度計算的意義會弱化。也就是說，約束最小二乘法在質(zhì)量較差網(wǎng)格上的優(yōu)勢更為突出，這對于解決一般混合網(wǎng)格在彎曲壁面附近加密困難的問題十分有效。

3.2 復(fù)雜混合網(wǎng)格應(yīng)用

由于ONERA-M6機翼具有簡單的幾何外形和復(fù)雜的跨聲速流動特征，因此是一個經(jīng)典的CFD驗證算例。圖5給出了其計算網(wǎng)格，機翼前后緣采用六面體以方便于法向和流向加密，機翼表面其他部分采用三棱柱以方便于法向加密，外圍采用四面體填充，四面體與六面體和棱柱之間用金字塔連接。網(wǎng)格總數(shù)為187萬，第1層網(wǎng)格高度約為機翼平均氣動弦長的10-5。Ma∞=0.839 5，迎角α=3.06°，相對于平均氣動弦長Re=1.172×107，壁面為絕熱壁，湍流模型采用Spalart-Allmaras模型[17]。

圖5 ONERA-M6計算網(wǎng)格Fig.5 Computational grid of ONERA-M6

圖6給出了采用混合消去加權(quán)法的對稱面和機翼表面的等密度線r分布，由圖可以清晰地看到陡峭的λ激波結(jié)構(gòu)。

圖6 ONERA-M6等密度線(混合消去加權(quán)法)Fig.6 ONERA-M6 density contour (mixed method)

圖7為ONERA-M6各站位壁面壓力系數(shù)Cp分布，橫坐標(biāo)x/c為各站位與氣動弦長的相對位置。從圖7中幾個站位的壓力分布對比可見，在本算例中壁面附近網(wǎng)格布置較為稀疏的條件下，約束最小二乘法相比原始最小二乘法的計算結(jié)果較好，特別是對前緣吸力峰值和激波位置的捕捉更為準(zhǔn)確，這主要是因為梯度計算過程強化了邊界的影響從而使流場重構(gòu)更加準(zhǔn)確。另外，混合消去加權(quán)法較單獨加權(quán)法的計算結(jié)果也更接近于實驗值，其差別主要在激波位置附近，原因是混合方法在前緣由于邊界模板擴(kuò)大使加權(quán)法生效，發(fā)揮了加權(quán)法的穩(wěn)定性優(yōu)勢；在激波位置附近由于模板只包含一個邊界使消去法生效，從而發(fā)揮了消去法的精度優(yōu)勢。

圖7 ONERA-M6各站位壁面壓力系數(shù)分布Fig.7 ONERA-M6 wall pressure coefficient distribution

4 結(jié) 論

1) 對不同邊界類型的透明化處理統(tǒng)一了約束最小二乘方程組的形式，從而大大簡化了該方法的應(yīng)用。

2) 針對層流平板的數(shù)值實驗，加權(quán)法的加權(quán)系數(shù)取到5基本已經(jīng)足夠大。

3) 約束最小二乘法強化了邊界對梯度計算的影響，相比原始最小二乘法流場重構(gòu)更加準(zhǔn)確，特別對邊界質(zhì)量較差網(wǎng)格的計算優(yōu)勢更為突出。

4) 混合采用消去法和加權(quán)法求解約束最小二乘方程，相比單獨采用加權(quán)法的計算結(jié)果有所改善。

[1] Wang Z J. High-order methods for the Euler and Navier-Stokes equations on unstructured grids. Progress in Aerospace Sciences, 2007, 43(1-3): 1-41.

[2] Mavriplis D J. Unstructured-mesh discretizations and solvers for computational aerodynamics. AIAA Journal, 2008, 46(6): 1280-1296.

[3] Barth T J, Jespersen D C. The design and application of upwind schemes on unstructured meshes. AIAA-1989-366, 1989.

[4] Haselbacher A, Blazek J. Accurate and efficient discretization of Navier-Stokes equations on mixed grids. AIAA Journal, 2000, 38(11): 2094-2101.

[5] Smith T M, Barone M F, Bond R B, et al. Comparison of reconstruction techniques for unstructured mesh vertex centered finite volume schemes. AIAA-2007-3958, 2007.

[6] Mavriplis D J. Revisiting the least-squares procedure for gradient reconstruction on unstructured meshes. AIAA-2003-3986， 2003.

[7] Ollivier-Gooch C, van Altena M. A high-order accurate unstructured mesh finite-volume scheme for the advection-diffusion equation. Journal of Computational Physics, 2002, 181(2): 729-752.

[8] Haselbacher A. A WENO reconstruction algorithm for unstructured grids based on explicit stencil construction. AIAA-2005-879, 2005.

[9] Haselbacher A. On constrained reconstruction operators. AIAA-2006-1274， 2006.

[10] Blazek J. Computational fluid dynamics principles and application. Oxford, UK： Elsevier Science Ltd, 2001: 5-25.

[11] Roe P L. Approximate Riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics, 1981, 43: 357-372.

[12] Venkatakrishnan V. On the accuracy of limiters andconvergence to steady state solutions. AIAA-1993-880, 1993.

[13] Aftosmis M, Gaitonde D, Tavares T S. Behavior of linear reconstruction techniques on unstructured meshes. AIAA Journal, 1995, 33(5): 2038-2049.

[14] Barth T J. A 3D upwind Euler solver for unstructured meshes. AIAA-1991-1548, 1991.

[15] Haselbacher A, Vasilyev O V. Commutative discrete filtering on unstructured grids based on least-squares techniques. Journal of Computational Physics, 2003, 187: 197-211.

[16] Zhu Y M, Wang Z Z. A new method for estimate ill-conditioning matrix. Journal of Shanghai Jiaotong University, 1992, 26(3): 110-112. (in Chinese)

朱揚明, 王志中. 病態(tài)矩陣判別的一種新方法. 上海交通大學(xué)學(xué)報, 1992, 26(3): 110-112.

[17] Spalart P R, Allmaras S R. A one-equation turbulence model for aerodynamic flows. AlAA-1992-439, 1992.

ConstrainedLeast-squaresReconstructionMethodforMixedGrids

KANGZhongliang,YANChao*

SchoolofAeronauticScienceandEngineering,BeihangUniversity,Beijing100191,China

Inordertoenforcetheeffectofboundaryconditionsonthesolutionreconstructionofcomputationalfluiddynamics(CFD)simulation,aconstrainedleast-squaresmethodbasedonarbitrarypolyhedralmixedgridsisinvestigated.Tomakethemethodsimpler,auniformconstrainedsystemisconstructedthroughtreatingdifferentboundaryconditionsinthesameway.Tomakeuseoftheadvantagesofboththeeliminationapproachandtheweightingapproach,amixedmethodisprovidedtosolvetheconstrainedsystem.Thenumericalexperimentsbasedonlaminarplatesshowthatsettingtheweightingcoefficientto5islargeenoughintheweightingapproach.ThenumericalresultsofasubsoniclaminarflowoveraplateandatransonicturbulentflowoverONERA-M6demonstratethattheconstrainedleast-squaresmethodismoreaccuratethantheoriginalleast-squaresmethod,especiallywhenbadqualitygridsareusedneartheboundary,andthatthemixedmethodshowsimprovedresultsascomparedwiththeweightingapproach.

mixedgrid;arbitarypolyhedral;reconstructionmethod;constrainedleast-squaresmethod;elimination;wei-

ghting

2011-10-23；Revised2011-12-22；Accepted2012-02-01；Publishedonline2012-02-291037

URL：www.cnki.net/kcms/detail/11.1929.V.20120229.1037.007.html

NationalBasicResearchProgramofChina(2009CB724104)

.Tel.：010-82317019E-mailchyan@vip.sina.com

2011-10-23;退修日期2011-12-22;錄用日期2012-02-01; < class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2012-02-291037

www.cnki.net/kcms/detail/11.1929.V.20120229.1037.007.html

國家“973”計劃(2009CB724104)

.Tel.：010-82317019E-mailchyan@vip.sina.com

KangZL,YanC.Constrainedleast-squaresreconstructionmethodformixedgrids.ActaAeronauticaetAstronauticaSinica,2012,33(9):1598-1605. 康忠良，閆超．適用于混合網(wǎng)格的約束最小二乘重構(gòu)方法．航空學(xué)報,2012,33(9):1598-1605.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

1000-6893(2012)09-1598-08

V211.3

A

康忠良男， 博士研究生。主要研究方向： 混合網(wǎng)格法、動網(wǎng)格法和大型CFD軟件平臺的研究及應(yīng)用。

Tel: 010-82318071

E-mail: KZL929@163.com

閆超男， 博士， 教授， 博士生導(dǎo)師。主要研究方向： 計算流體力學(xué)和高超聲速空氣動力學(xué)。

Tel: 010-82317019

E-mail: chyan@vip.sina.com