適用于混合網格的約束最小二乘重構方法

康忠良， 閆超

北京航空航天大學 航空科學與工程學院， 北京 100191

適用于混合網格的約束最小二乘重構方法

康忠良， 閆超*

北京航空航天大學 航空科學與工程學院， 北京 100191

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

混合網格； 任意多面體； 重構方法； 約束最小二乘法； 消去法； 加權法

在計算流體力學(CFD)的應用中，采用混合網格可以充分利用非結構網格的網格生成優勢和類結構網格的計算優勢，因此是一種先進有效的方法。混合網格的單元類型可以是常用的四面體、金字塔、三棱柱或六面體，也可以是這些單元之間以及與其他任意多面體的混合。

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

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

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

1 格心有限體積法

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

1.1 控制方程

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

(1)

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

1.2 數值離散方法

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

|ARoe|(QR-QL)]

(2)

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

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

{

(3)

線性重構方法在規整網格上可以使空間離散達到二階精度[13]。假如梯度是準確值，該方法在任意網格上都可以準確重構，因此梯度計算是重構過程的關鍵。

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

2 約束最小二乘法

由1.2節可知，梯度計算是線性重構的關鍵，最小二乘法是目前常用的一種梯度計算方法，且非常適合于混合網格計算，下面將基于原始最小二乘法給出一種新型的約束最小二乘法。

2.1 算法構造

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

(4)

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

(5)

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

(6)

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

Ax=b

(7)

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

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

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

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

(8)

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

2.2 方程求解

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

(9)

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

(10)

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

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

(11)

式中：w為大于1的加權法權系數，在保證系數矩陣不病態的前提下，可對其適當加大來強化邊界對梯度計算的貢獻。

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

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

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

3 算例驗證及分析

3.1 計算精度驗證

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

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

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

圖2 平板計算網格G1Fig.2 Computational grid G1 of flat plate

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

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

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

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

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

3.2 復雜混合網格應用

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

圖5 ONERA-M6計算網格Fig.5 Computational grid of ONERA-M6

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

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

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

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

4 結 論

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

2) 針對層流平板的數值實驗，加權法的加權系數取到5基本已經足夠大。

3) 約束最小二乘法強化了邊界對梯度計算的影響，相比原始最小二乘法流場重構更加準確，特別對邊界質量較差網格的計算優勢更為突出。

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

[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)

朱揚明, 王志中. 病態矩陣判別的一種新方法. 上海交通大學學報, 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">網絡出版時間

時間：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. 康忠良，閆超．適用于混合網格的約束最小二乘重構方法．航空學報,2012,33(9):1598-1605.

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

1000-6893(2012)09-1598-08

V211.3

A

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

Tel: 010-82318071

E-mail: KZL929@163.com

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

Tel: 010-82317019

E-mail: chyan@vip.sina.com