一種驗證網(wǎng)格質(zhì)量與CFD計算精度關系的方法

董亮，劉厚林，談明高，王勇，王凱

(江蘇大學 流體機械工程技術研究中心，江蘇 鎮(zhèn)江，212013)

網(wǎng)格質(zhì)量對計算的精度和效率有顯著的影響。劣質(zhì)單元的出現(xiàn)將導致在CFD計算時出現(xiàn)剛度矩陣病態(tài)，嚴重影響計算結(jié)果的精度，甚至可能導致無法計算[1-3]。因此，在CFD計算之前，需要對網(wǎng)格質(zhì)量進行評估，以判斷其是否能夠滿足求解器要求。四面體網(wǎng)格因其對復雜邊界有較強的適應性，已成為CFD數(shù)值計算中最常用的網(wǎng)格單元。雖然公認正四面體是四面體單元的理想形狀，但長期以來，對四面體單元質(zhì)量的度量或評判并沒有公認的或絕對的標準[4-6]。在計算幾何和計算科學等領域，所謂“高質(zhì)量”單元并沒有一個明確的定義，而是將與規(guī)則形狀偏離較大的單元認為是“劣質(zhì)”單元，而將接近于規(guī)則形狀的單元認為是“高質(zhì)量”單元。目前，大多數(shù)針對網(wǎng)格質(zhì)量衡量準則的研究工作，主要集中在各衡量準則對于劣質(zhì)單元的評判效果以及它們是否等價等方面[7-10]。如Parthasarathy等[11]研究了7種常用的網(wǎng)格質(zhì)量衡量準則判別劣質(zhì)單元的敏感性問題。Liu等[12]對3種網(wǎng)格質(zhì)量衡量準則之間的關系進行了研究，認為這些網(wǎng)格質(zhì)量衡量準則在某種意義上是相互等價的，即對形狀較差單元，若一種度量值趨于0，則其他2種度量值也趨于0。文獻[13]也認為這幾種衡量準則是相互等價的，并通過算例表明，采用任何一種衡量準則優(yōu)化，都可以獲得更高質(zhì)量的網(wǎng)格。聶春戈等[14]通過數(shù)值試驗發(fā)現(xiàn)，采用不同的準則評判單元形狀的變化,有可能產(chǎn)生矛盾的評判效果，從而影響網(wǎng)格質(zhì)量優(yōu)化的過程和結(jié)果。然而，網(wǎng)格單元質(zhì)量的提高與CFD計算精度提高的關系研究較少。本文作者以90°彎管為研究對象，探討構建網(wǎng)格獨立解的方法，并對比4種網(wǎng)格質(zhì)量衡量準則與計算精度的關系。

1 四面體質(zhì)量衡量準則

判斷單元質(zhì)量最直接的標準是數(shù)值計算模塊利用該單元進行求解時的精度和速度。這個衡量標準并不適合直接作為網(wǎng)格生成和網(wǎng)格優(yōu)化調(diào)整的準則，因此，從不同角度考慮的四面體單元質(zhì)量衡量準則被提出。一般而言，一個合理的四面體網(wǎng)格質(zhì)量衡量準則應滿足以下原則：

(1)單元的平移、旋轉(zhuǎn)、均勻縮放都應不改變其度量值；

(2)退化單元(單元體積為0)的質(zhì)量度量值為0，正四面體的度量值最大為1；

(3)反轉(zhuǎn)單元(單元體積為負)的度量值為負；

(4)衡量準則應為四面體各個節(jié)點坐標的函數(shù)，并且能夠反映單元形狀的變化；

(5)能夠檢測出易導致計算發(fā)散的劣質(zhì)單元。

各種質(zhì)量衡量準則如表1所示。其中，QEVS和QEAS 2種網(wǎng)格質(zhì)量衡量準則為商業(yè)網(wǎng)格劃分軟件GAMBIT[15]中EquiSize Skew和EquiAngle Skew準則的變形，它們分別用來驗證單元等尺寸和等角度的失真程度。同時，本文提出了QNEW1和QNEW2 2種質(zhì)量衡量準則，用來驗證單元的形狀和尺寸特性，大量數(shù)值實驗表明它們滿足作為一個合理質(zhì)量衡量準則的所有特性。衡量準則的度量值越大，說明單元質(zhì)量越好；且當單元為正四面體時度量值達到最大值1，單元為退化單元時度量值達到最小值0。

表1 各種單元質(zhì)量衡量準則Table1 Different quality metrics

2 網(wǎng)格獨立解

2.1 計算模型及數(shù)值方法

為了驗證網(wǎng)格質(zhì)量衡量準則與計算精度之間的關系，以90°彎管作為研究對象，其幾何尺寸如圖1所示[16]。為了便于分析，將彎管分成上游直線段、彎曲段和下游直線段3部分。并定義彎曲段的主流入口截面處θ=0°，彎曲段出口截面處θ=90°，坐標系原點O位于圓管的圓心，有關90°彎管具體參數(shù)如表2所示。

圖1 90°彎管結(jié)構示意圖Fig.1 Sketch of 90°bending duct and definition of coordinates

應用商業(yè)軟件Fluent，計算采用RNGk-ε湍流模型，非耦合隱式方案進行求解。通過有限體積法進行空間離散，對控制方程中的源項和擴散項應用二階中心差分格式，對流項應用二階迎風格式，采用SIMPLE算法來實現(xiàn)速度和壓力的耦合求解。具體邊界條件如下。

進口邊界條件：流體均勻地進入管道中，初始速度為1 m/s。

出口邊界條件：壓力出口。

壁面條件：固體壁面采用邊壁無滑移條件；靠近固體邊界區(qū)域，采用標準壁面函數(shù)處理。

表2 90°彎管參數(shù)Table2 90° bending duct parameters

2.2 獲得網(wǎng)格獨立解的方法

2.2.1 算法流程

本文網(wǎng)格獨立解的獲得是基于各套網(wǎng)格數(shù)值解計算得到的相應速度誤差。通過調(diào)整網(wǎng)格步長使網(wǎng)格加密，然后比較相鄰網(wǎng)格中各單元的速度，進而得到相應的誤差。當所有單元的平均相對誤差小于某一規(guī)定值ε時，即可認為計算結(jié)果是網(wǎng)格獨立解。具體過程如下。

步驟1 生成初始四面體(作為第k套網(wǎng)格)，進行數(shù)值計算并記錄網(wǎng)格中各節(jié)點的速度(作為第k套網(wǎng)格數(shù)值解)；

步驟2 按照第k套網(wǎng)格單元數(shù)的1.3倍增加網(wǎng)格單元數(shù)作為第k+1套網(wǎng)格單元數(shù)，進行數(shù)值計算并記錄網(wǎng)格中各節(jié)點的速度(作為第k+1套網(wǎng)格數(shù)值解)；

步驟3 采用三維線性插值算法將k+1套數(shù)值解映射到k套的數(shù)值解中，獲得第k+1套網(wǎng)格在第k套網(wǎng)格中相應的數(shù)值解；

步驟4 計算第k+1套與第k套的速度相對誤差：

其中：vk，i為第k套網(wǎng)格節(jié)點i的速度；vK+1，i為第k+1套網(wǎng)格映射后在節(jié)點i處的速度；N為k套網(wǎng)格節(jié)點總數(shù)。

步驟5 如果小于給定的誤差ε(本文取ε=2%)，那么就認為第k套網(wǎng)格計算結(jié)果是網(wǎng)格獨立解，否則重復步驟1～4。

2.2.2 映射獲得相鄰數(shù)值解的算法

將k+1套網(wǎng)格的數(shù)值解映射到k套數(shù)值解中，獲得第k+1套網(wǎng)格在第k套網(wǎng)格中相應數(shù)值解的具體過程如下。

步驟1 將k套網(wǎng)格中的所有節(jié)點的坐標存儲到Φ1中；

步驟2 將k+1套網(wǎng)格中的所有節(jié)點的坐標和速度存儲到Φ2中；

步驟3 依次讀取Φ1中的節(jié)點P，并以該節(jié)點為球心，R為半徑，找出Φ2中所有位于該球內(nèi)的節(jié)點，如果Φ2中存在于P點重合的節(jié)點，則將其速度存儲到Ω中，并接著執(zhí)行步驟7；

步驟4 判斷位于球內(nèi)節(jié)點的個數(shù)，若小于給定值N，則繼續(xù)增加區(qū)間，直到位于球內(nèi)的節(jié)點數(shù)大于N為止；

步驟5 若位于球內(nèi)的節(jié)點數(shù)大于N，則只取最近的N個節(jié)點；

步驟6 采用三維線性插值算法計算出節(jié)點P在k+1套網(wǎng)格相同位置的速度，并將其存儲到Ω中；

步驟7 重復步驟3～6，直到Φ1中沒有節(jié)點為止。

其中R和N為預先給定值，從模型尺寸以及計算精度考慮，本文取R=0.015m，N=20。

2.3 結(jié)果統(tǒng)計及分析

圖2所示為90°彎管在不同網(wǎng)格等級下的網(wǎng)格單元以及節(jié)點數(shù)，其中相鄰網(wǎng)格的單元數(shù)比值大約為1.3。各套網(wǎng)格的最差單元質(zhì)量和平均網(wǎng)格單元質(zhì)量如圖3所示。從圖3可以看出，隨著單元數(shù)的增加，網(wǎng)格最差單元質(zhì)量以及平均網(wǎng)格質(zhì)量都隨之增加。

圖4所示為相鄰網(wǎng)格之間計算得到的速度相對誤差，從圖4可以看出，網(wǎng)格數(shù)較小時，相鄰網(wǎng)格之間速度存在較大的誤差，如第2套網(wǎng)格與第1套網(wǎng)格的速度相對誤差達到24.59%，然而，隨著網(wǎng)格數(shù)的增加，相鄰網(wǎng)格之間的速度相對誤差逐漸減少，第10套網(wǎng)格與第9套網(wǎng)格之間的速度相對誤差為1.64%，小于給定的誤差2%，因此，可以將第9套網(wǎng)格計算得到的數(shù)值解作為90°彎管的網(wǎng)格獨立解。

圖2 網(wǎng)格劃分結(jié)果Fig.2 Mesh generation results

圖3 各網(wǎng)格等級質(zhì)量Fig.3 Quality for different levels of mesh

圖4 相鄰網(wǎng)格之間計算得到的速度相對誤差Fig.4 Average relative error between adjacent meshes

3 分析與討論

3.1 驗證方法

將單元質(zhì)量的度量值以及相對于獨立解的速度相對誤差組成一個數(shù)據(jù)庫，通過研究該數(shù)據(jù)庫來獲得網(wǎng)格質(zhì)量衡量準則與計算精度之間的關系。具體步驟如下。

步驟1 獲得不完善網(wǎng)格(前8套網(wǎng)格)中各單元的網(wǎng)格質(zhì)量度量值(采用表1中的4種衡量準則分別進行計算)；

步驟2 獲得不完善網(wǎng)格中各單元相對于獨立解(第9套網(wǎng)格)的速度相對誤差(與2.2節(jié)方法一致)；

步驟3 將算例中所有單元質(zhì)量的度量值和相對于獨立解的誤差建立一個數(shù)據(jù)庫；

步驟4 按照單元的速度相對誤差從小到大順序?qū)?shù)據(jù)庫中的數(shù)據(jù)重新排列；

步驟5 將該數(shù)據(jù)庫分為若干區(qū)間(本文取20)，并計算各個區(qū)間的統(tǒng)計數(shù)據(jù)值(本文采用該區(qū)間中數(shù)據(jù)的平均值或最小值)。

3.2 結(jié)果統(tǒng)計

圖5所示為將前8套網(wǎng)格中所有單元質(zhì)量度量值和計算誤差建立的1個數(shù)據(jù)庫，獲得4種質(zhì)量衡量準則與計算精度之間的關系。從圖5可以看出：當采用區(qū)間最小值表示該區(qū)間的趨勢時，這4種衡量準則與計算精度沒有任何規(guī)律。而當以區(qū)間平均值表示該區(qū)間的趨勢時，QEVS和QEAS這2種準則在誤差較小的區(qū)間內(nèi)(＜0.1)，與計算精度有較好的關系(隨著誤差的增加，單元的度量值降低)；然而在其他區(qū)間內(nèi)，隨著誤差的增加，它們的度量值出現(xiàn)了波動，也就是說在網(wǎng)格生成或優(yōu)化的過程中，使用這2種準則可能會使網(wǎng)格朝著誤差較大的方向進行優(yōu)化。對于QNEW1和QNEW2，除了某些區(qū)間外(0.1＜＜0.12)，在其他區(qū)間內(nèi)，它們與計算精度具有較好的關系，尤其在計算誤差較大的區(qū)間內(nèi)(＞0.15)，變化更加明顯。因此，QNEW1和QNEW2這2種衡量準則能夠較好地區(qū)分計算誤差較大的單元，更加準確地指導網(wǎng)格的生成。

圖5 網(wǎng)格質(zhì)量衡量準則與計算精度之間的關系Fig.5 Relationship between four metrics and errors

圖6 QNEW1與QNEW2計算精度之間的關系Fig.6 Relationship between two metrics (QNEW1 and QNEW2)and errors

為了進一步驗證QNEW1和QNEW2這2種衡量準則與計算精度之間的關系，以第4套與第7套算例作為研究對象，且采用區(qū)間的平均值表示該區(qū)間的趨勢，結(jié)果如圖6所示。由圖6可知：兩者的變化規(guī)律與圖5(c)和5(d)變化規(guī)律基本一致，都是在誤差較大的區(qū)間內(nèi)，QNEW1和QNEW2 2種衡量準則與計算精度有較好的關系，且在相同誤差的情況下，QNEW2的度量值大于QNEW1的度量值。

4 結(jié)論

(1)提出了2種四面體網(wǎng)格質(zhì)量衡量準則QNEW1和QNEW2，研究了它們及衡量準則QEAS，QEVS與CFD計算精度之間的關系。

(2)為了準確地驗證網(wǎng)格質(zhì)量衡量準則與計算精度之間的關系，提出了一種通過數(shù)值計算構建網(wǎng)格獨立解的方法，由于該方法考慮了每個單元的計算誤差，因而獲得的獨立解更加接近真實解。

[1]Munson T.Mesh shape-quality optimization using the inverse mean-ratio metric[J].Mathematical Programming, 2007, 110(3):561-590.

[2]黃曉東, 杜群貴.二維有限元網(wǎng)格的局部加密方法[J].華南理工大學學報: 自然科學版, 2004, 32(12): 56-60.HUANG Xiao-dong, DU Qun-gui.Local refinement of 2D finite element meshes[J].Journal of South China University of Technology: Natural Science Edition, 2004, 32(12): 56-60.

[3]Park J, Shontz S M.Two derivative-free optimization algorithms for mesh quality improvement[J].International Conference on Computational Science, 2010, 1(1): 387-396.

[4]Hetmaniuk U, Knupp P.A mesh optimization algorithm to decrease the maximum interpolation error of linear triangular finite elements[J].Engineering with Computers, 2010, 27(1):3-15.

[5]Knupp P M.Algebraic mesh quality metric[J].SIAM Journal on Scientific Computing, 2001, 23(1): 193-218.

[6]Babuska I, Flaherty J E, Henshaw W D, et al.Modeling, mesh generation, and adaptive numerical methods for partial differential equations[M].New York: Springer, 1995: 97-127.

[7]Gosselin S, Ollivier-Gooch C.Tetrahedral mesh generation using Delaunay refinement with non-standard quality measures[J].International Journal for Numerical Methods in Engineering,2011, 87(8): 795-820.

[8]Dompierre J, Vallet M G, Labbé P, et al.An analysis of simplex shape measures for anisotropic meshes[J].Computer Methods in Applied Mechanics and Engineering, 2005, 194(48/49):4895-4914.

[9]Freitag L A, Knupp P M.Tetrahedral element shape optimization via the Jacobian determinant and condition number[C]//Proceedings of the 8th International Meshing.California, 1999:247-258.

[10]董亮, 劉厚林, 談明高, 等.離心泵四面體網(wǎng)格質(zhì)量衡量準則及優(yōu)化算法[J].西安交通大學學報, 2011, 45(11): 31-36.DONG Liang, LIU Hou-lin, TAN Ming-gao, et al.Quality measurement criteria and optimization algorithm of tetrahedral mesh for centrifugal pumps[J].Journal of Xi’an Jiao Tong University, 2011, 45(11): 31-36.

[11]Parthasarathy V N, Graichen C M, Hathaway A F.A comparison of tetrahedron quality measures[J].Finite Elements in Analysis and Design, 1993, 15(3): 255-261.

[12]Liu A, Joe B.Relationship between tetrahedron shape measures[J].BIT Numerical Mathematics, 1994, 34(2):268-287.

[13]Lo S H.Optimization of tetrahedral meshes based on element shape measures[J].Computers& Structures, 1997, 63(5):951-961.

[14]聶春戈, 劉劍飛, 孫樹立.四面體網(wǎng)格質(zhì)量度量準則的研究[J].計算力學學報, 2003, 20(5): 579-582.NIE Chun-ge, LIU Jian-fei, SUN Shu-li.Study on quality measures for tetrahedral mesh[J].Chinese Journal of Computational Mechanics, 2003, 20(5): 579-582.

[15]Fluent Inc, Gambit User’s Guide[EB/OL].[2003-12-30].http://combust.hit.edu.cn:8080/fluent/Gambit13_heop/users_guide.

[16]Taylor A, Whitelaw J H.Curved ducts with strong secondary motion: Velocity measurements of developing laminar and turbulent flow[J].Journal of Fluids Engineering, 1982, 104:350-359.