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

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

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

1 四面體質量衡量準則

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

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

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

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

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

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

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

表1 各種單元質量衡量準則Table1 Different quality metrics

2 網格獨立解

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

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

圖1 90°彎管結構示意圖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 獲得網格獨立解的方法

2.2.1 算法流程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 網格劃分結果Fig.2 Mesh generation results

圖3 各網格等級質量Fig.3 Quality for different levels of mesh

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

3 分析與討論

3.1 驗證方法

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

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

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

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

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

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

3.2 結果統(tǒng)計

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

圖5 網格質量衡量準則與計算精度之間的關系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ū)間的趨勢，結果如圖6所示。由圖6可知：兩者的變化規(guī)律與圖5(c)和5(d)變化規(guī)律基本一致，都是在誤差較大的區(qū)間內，QNEW1和QNEW2 2種衡量準則與計算精度有較好的關系，且在相同誤差的情況下，QNEW2的度量值大于QNEW1的度量值。

4 結論

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

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

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

[2]黃曉東, 杜群貴.二維有限元網格的局部加密方法[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]董亮, 劉厚林, 談明高, 等.離心泵四面體網格質量衡量準則及優(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]聶春戈, 劉劍飛, 孫樹立.四面體網格質量度量準則的研究[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.