網格對CFD模擬結果的影響分析

劉明亮，張思青，李勝男

（1.昆明理工大學冶能學院，云南省昆明市 650093;2.云南電網有限責任公司電力科學研究院，云南省昆明市 650000）

網格對CFD模擬結果的影響分析

劉明亮1，張思青2，李勝男3

（1.昆明理工大學冶能學院，云南省昆明市 650093;2.云南電網有限責任公司電力科學研究院，云南省昆明市 650000）

在CFD模擬的過程中，高質量的網格對計算結果的影響至關重要。本文通過一個繼電器內部流場的實例，描述了網格獨立性檢驗的詳細過程。通過Fluent雙精度求解器，應用k-omega模型對比添加與不添加邊界層時計算誤差，進而探討邊界層網格對計算精度的影響。

網格無關解；網格獨立性檢驗；k-omega模型；邊界層網格劃分

0 引言

作為流體動力學研究方法之一的CFD模擬法，隨著理論研究和實際應用的深入以及商用軟件的日漸成熟，越來越受到工程應用者的重視，甚至產生一定的依賴。這就對計算可信度的要求也越來越大，國內已有大量的文獻對不同方法給出誤差的理論分析[1]-[5]。

但不管理論如何完善，實際計算中的大量誤差還是避免不了的，如不確定的邊界條件、近似的計算模型、簡化的幾何模型、低質量的網格等，以上幾個因素的綜合，可能會導致計算結果毫無意義。通常的解決辦法是配以模型試驗的數據驗證。然而，并不是所有的CFD計算都有可行的實驗來驗證，數值模擬結果的誤差本身是可以接受的，不過這需要CFD軟件應用者能深刻的理解計算模型的原理，每個模型的適用條件，從而可以最大程度的減少人為誤差。

實際計算中，我們可以用商用軟件畫出質量很好的網格，但是網格質量高，并不代表網格數量足夠滿足計算精度的要求。找到合適計算的網格數量需要網格獨立性檢驗，進而得到的結果叫網格無關解。本文就是在假定模型選擇合理的條件下，研究網格對計算結果的影響。

1 誤差理論

CFD的模擬就是用計算機數值方法，求解控制流體運動的包含初邊值問題的偏微分方程。若不考慮模型誤差，數值計算中代入的誤差主要包括：用有窮表示無窮而產生的截斷誤差和由計算機的位寬有限因估值而產生的舍入誤差。

CFD計算實質上是將連續的初邊值問題離散化，轉化為離散的初邊值問題。常見的離散方法有有限元、有限體積、有限差分等方法，理論上都是網格數越多，結果越接近真值。當網格無限多且無限小時，差分方程將變成精確的理論解。

實際上：網格不可能無限制的加密。主要的問題有：首先，網格越密，則計算量越大，計算周期也越長，通常我們的計算資源總是有限的。其次，隨著網格的加密，為了計算的比斂時間步長也必須減小。實際的經驗告訴我們當網格尺寸及時間步長減小時，截斷誤差確實減小，而舍入誤差卻反而增加（如圖1所示）。持續減小步長尺寸并不意味著能獲得更加精確的結果。事實恰好相反，在很小的時間步長處，因為舍入誤差的增加，使計算精度反而下降。計算機浮點運算造成的舍入誤差也會增大。因此在實際應用中，我們尋找的就是計算精度與計算開銷間一個平稀點，也就是達到網格無關的閾值[6]。

所謂網格無關性檢驗，就是驗證網格數獨立于計算精度。通常是找到一個滿足CFD計算的最少網格數。執行時，讓網格密度逐漸增加，并選擇一個判斷參數，觀察其隨網格數量增加的變化情況，當網格的增加幾乎不再影響該參數的變化時，網格數量達到計算要求。

圖1 網格尺寸或時間步長與離散誤差、舍入誤差、總誤差的函數關系Fig.1 Function relation between mesh size or time step and discrete error，rounding error，total error

2 計算實例

2.1 幾何模型

本算例所選計算模型是國產沈陽四興QJ-80型瓦斯繼電器，瓦斯繼電器是變壓器的非電量保護重要的元器件，用三維繪圖軟件ProE建立的幾何模型，用ANSYS中的幾何處理模塊Design Modeler提取幾何流道生成流體計算域如圖2所示，相應的部件部分經過簡化，標記在該圖上。

圖2 瓦斯繼電器三維模型（上）與流體計算域三維模型（下）Fig.2 Buchholz relay three-dimensional model（on）and computational fluid domain three-dimensional model（down）

繼電器的工作原理見示意簡圖（如圖3所示），當變壓器內部發生嚴重故障（絕緣擊穿、匝間短路、鐵芯事故等）時，產生強烈的瓦斯氣體，變壓器油箱內壓力瞬時突增，產生很大的流向油枕方向的油流沖擊，沖擊繼電器中擋板，產生擋板動作，進而觸發報警信號。

圖3 瓦斯繼電器工作原理示意圖Fig.3 Schematic diagram of the Buchholz relay working principle

本算例的目的是為精確計算擋板上油流沖擊產生的阻力力矩而準備高質量網格。下面先從沒有邊界層的網格起，到添加邊界層的網格，用ICEMCFD軟件來詳細說明準備計算網格的過程。

2.2 無邊界層網格的生成

進行網格獨立性檢驗首先要確定判斷參數。通常不可壓計算，選用進出口壓差來作為判斷收斂的參數。當進口的壓強不變時，壓差只取決于出口值，因此本算例用出口壓強面積加權平均值來作為判斷收斂參數。

其次要繪制初始網格，初始的網格不需要太細致。模型的入口管徑為80mm，因為我們關心擋板處的結果，因此將擋板處最大網格邊長設為5mm，其余部分最大網格邊長定為8mm。生成的網格記為“網格1”。

下面對網格進行加密。ICEM-CFD中，功能“Scale Factor”用來改變全局的網格尺寸（體、面、線），通過將它乘以設置的參數來得到實際網格參數。即可以通過調整Scale Factor值的大小來改變網格密度。

我們選用加倍的網格加密原則，即加密后的網格數量是加密前的兩倍。為達到該目的，只需另加密后的邊長變為原來的，這個值約等于0.7937。第二次加密就使Scale Factor變為0.79372，以此類推，第n次加密，Scale Factor 的值就為0.7937n。

以網格1為基礎，依據以上原則，生成網格2，網格3，一直到網格5，每次畫好的網格再經過光順，保證網格質量在0.4以上。用 fluent計算每個網格，觀察出口相對靜壓值，并將相應的數據匯總于表1中，并以圖示的形式展示在圖4中。

表1 出口相對靜壓與網格數的關系Tab.1 Relationship between export relatively static and number of grid

圖4 出口相對靜壓隨網格數的變化Fig.4 Exports relatively static pressure changes with the number of grid

從圖4可以看出當網格達到1337948時，出口壓強平均值已經逐漸穩定。從網格4到網格5，出口壓強變化很小，也就是說，當網格數量達到1337948時，再增加網格數量已經不會顯著提高計算精度，進而可以判斷，網格4在一定程度上已經滿足計算要求。

2.3 邊界層網格的生成

邊界層網格對定量的計算問題，如表面摩擦，升力系數阻力系數的求取影響比較大，對換熱問題尤為顯著。

盡管在工程精度允許的范圍內，有時可以對復雜問題不畫邊界層，取而代之的是加密邊界層網格，但由于本問題屬于大雷諾數流動，且我們關注的就是阻力系數，因此有必要添加邊界層網格，來提高計算精度。

Fluent中k-epsilon模型是針對充分發展的湍流模型，即高Re數的湍流模型。對近壁區內的流動，Re數比較低，湍流的發展并不充分，湍流的脈動影響不如分子黏性的影響大，這樣，在這個區域內就不能直接使用該模型，通常采用壁面函數法或近壁模型法。

壁面函數法本身并不對貼壁層進行求解，而是用一組半經驗公式來將壁面與充分發展的湍流區聯系起來。其缺點在于，當流動情況偏離了壁面函數的理想條件時，如強體積力流動，近壁區三維性很強的流動問題，計算誤差就會明顯增大。其次壁面函數的缺點還在于：沿壁面法向細化網格時，會導致數值結果惡化。當y+小于15時，將會在壁面剪切力及熱傳遞方面逐漸導致產生無界錯誤。

近壁模型法。修改湍流模型以使其能夠求解近壁黏性影響區域，包括黏性底層。本文使用的方法即近壁模型法。近壁模型不需要使用壁面函數，但需要壁面網格的精細。 fluent中的k-omega湍流模型就是一種典型的近壁湍流模型。但要保證y+＜=1，以1為最佳。

根據給定的y+值，計算第一層網格與壁面的距離方法如下：

（1）計算雷諾數Re：

（2）計算表面摩擦系數：

（3）計算壁面切應力：

（4）計算摩擦速度：

（5）計算第一層網格高度：

繼電器所使用的工質為25號變壓器油，其密度ρ=895kg/m3，運動黏度ν=9.6×10-6m2/s。計算出動力黏度μ=νρ=8.592×10-3。計算中y+取1，特征速度U取2m/s。特征長度L取方形擋板的邊長0.05m。代入式（1）～式（5）得第一層網格高度ΔS為8.1517×10-5m。

使用近壁模型法時，覆蓋邊界層的最小網格數量在 10層左右，最好能達到20層。還有一點需要注意的是，提高邊界層分辨率，常常可以取得穩健的數值計算結果，因為只需要細化壁面法向方向網格。對于非結構網格，建議劃分10～20層棱柱層網格以提高壁面邊界層的預測精度。棱柱層厚度應當被設計為保證有15層或更多網格節點。另外，棱柱層大于邊界層厚度是必要的，否則棱柱層會限制邊界層的增長。這可以在獲得計算結果后，通過查看邊界層中心的最大湍流黏度，該值提供了邊界層的厚度（最大值的兩倍位置即邊界層的邊）。

平板湍流邊界層厚度計算公式[7]為：

運動黏度ν=9.6×10-6，U=2m/s，L=0.05，代入式（6）得δ=3.4997452×10-3m。邊界層網格厚度S應大于平板邊界層取3.5×10-3。

假設網格按指數律分布，設總網格層數為n（=15），公比為r。

由式（7）可反解出公比r（=1.1374）。

至此，邊界層參數首層網格高度h，公比r，層數n均得到，依據這些參數生成邊界層棱柱網格。畫好的網格如圖5所示。

圖5 計算網格剖面及邊界層網格Fig.5 Compute Grid profile and boundary layer mesh

本文劃分邊界層后還有一個好處，由于本算例中繼電器內部擋板是可以轉動的，當擋板轉動后，fluent的局部重構和彈性光順算法動態生成的新網格，難以保證計算y+值的要求。所以定義邊界層網格，并在 fluent將邊界層指定為隨同擋板一同運動的流體域。這就可以保證邊界層不會因擋板的運動而破壞[8]。

2.4 邊界條件

以入口管徑80mm為特征長度，由式（1）可求出入口處雷諾數Re為10416.6，大于內部流動臨界雷諾數Recr=2320，屬于湍流，選用低雷諾數模型中的標準k-omega模型。

圖中相應的邊界均用文字加以說明。在 fluent中邊界in表示進口。將其定義為 fluent中的速度入口條件velocity-inlet，取入口測試速度1.3m/s。Flunet對于亞聲速流動，速度入口邊界條件直接忽略壓強的大小，故模擬中無壓強條件，選用默認值即可。以下公式均來自ANSYS Fluent 14.5幫助手冊。

（1）湍流強度I的計算：

（2）湍流尺度I的計算：

（3）湍動能k的計算：

（4）特性耗散率w的計算：

上式中為模型中一個經驗常數其值近似為0.09。

依次求得式（9）～式（11）可得入口邊界條件湍動能k=0.0152040。特性耗散率ω=40.2003909/m。

由于是不可壓流體，出口速度和壓強未知，出口邊界“out”只能選用out flow條件。這實際上假定流場已充分發展，各個變量（壓力除外）沿流動方向的梯度值為零。其他邊界均定義為壁面。邊界“wall”定義為默認的速度為零的壁面條件。

流場的計算采用SIMPLEC算法，梯度項、壓力及動量方程均用二階格式以提高精度。設置收斂殘差為1×10-6。為提高計算精度，選用雙精度求解器。

2.5 結果分析

首先看y+值分布。影響精度的只是擋板壁面的y+值，在這里首先觀察擋板的y+值云圖，如圖6所示，左邊是正面迎接沖擊油流側，雷諾數大，導致相應的y+值也大，右邊是擋板背向油流側，雷諾數小，相應的的y+值也小。

圖6 y+值云圖Fig.6 cloud of y+ value

為進一步衡量最大y+值是否滿足低雷諾數模型k-omega模型的網格要求，即y+≤1。在擋板的正面取豎直線段up-down和水平線段left-right。分別繪制兩條線段上y+值的散點圖如圖7所示。由圖7可知，擋板上y+值除少數點大于1外，絕大多數均滿足條件。所以擋板壁面處網格質量滿足計算要求。

圖7 線段left-right和up-down上y+值散點圖Fig.7 y+ value on the scatterplot in line left-right and updown

在工程計算中對精度要求高時，網格和模型的選擇尤為重要。本文參考文獻[9]中給出第一層網格高度對敏感參數產生決定性的影響。本文結合理論計算來對比有無邊界層網格的差異。

由于前邊繪制的網格參考y+值是在特征速度為2m/s的條件進行的。由公式（1）～（5）可得，在一定范圍內，y+值會隨特征速度的增加而增加，為保證本計算y+值不會顯著增高，滿足y+＜=1的要求，我們分別在特征速度小于2m/s，取為1m/s、1.25m/s、1.5m/s、1.75m/s不同工況下，在 fluent中監控擋板表面摩擦系數。并將由式（2）計算的表面摩擦系數匯總于表2。

表2 表面摩擦因數計算結果

由圖8數據整理可見，有邊界層的數值計算結果和理論值吻合良好，沒有邊界層則使結果產生較大偏差。由此可見邊界層網格對計算結果的影響很大。合理選擇首層網格高度繪制的邊界層網格，會大大提高計算的精度。

如圖9所示，詳細給出有邊界層時，計算流場的細節。

圖8 表面摩擦系數對比

圖9 速度大小云圖及流場的流線圖Fig.9 Streamline velocity magnitude contours and streamline of flow field

3 結論

本文討論了網格數量對CFD模擬精度的影響，用進出口壓差作為對比參數，用倍增網格數的方法，對比不同網格密度對壓差的影響，找到適合計算的最佳網格數。并嚴格按照計算模型的要求，繪制有助于提高計算精度的邊界層網格，使壁面網格y+值滿足模型準則要求。

在以表面摩擦系數理論計算值為參照，對比有無邊界層網格時得到如下結論：邊界層網格對計算結果的影響很大。合理選擇首層網格高度繪制的邊界層網格，會大大提高計算的精度。本文的意義還在于生成了分析計算所需的高質量的網格。為后續進一步計算（動網格計算）做好鋪墊。

[1]康順，石磊，戴麗萍，范忠瑤.CFD模擬的誤差分析及網格收斂性研究[J].工程熱物理學報，2010，12：2009-2013.

[2]劉重陽，于芳，徐讓書.CFD計算網格誤差分析的一個算例[J].沈陽航空工業學院學報，2006，04：21-24.

[3]楊永.基于網格計算的CFD模擬可信度分析[A].中國空氣動力學會、中國宇航學會、中國航空學會、中國力學學會.計算流體力學研究進展——第十二屆全國計算流體力學會議論文集[C].中國空氣動力學會、中國宇航學會、中國航空學會、中國力學學會：2004：6.

[4]康順，劉強，祁明旭.一個高壓比離心葉輪的CFD結果確認[J].工程熱物理學報，2005，03：400-404.

[5]晏麗琴，徐宏彤.基于CFD軟件模擬精度與可信度影響因素分析[J].自動化與儀器儀表，2014，09：101-103.

[6]（澳）Jiyuan Tu Guan Heng Yeoh，（美）Chaoqun Liu 著，王曉冬譯.計算流體力學——從實踐中學習[M].第1版.沈陽：東北大學出版社，2009，10 ：147-148.

[7]陳卓如，金朝銘，王洪杰，王成敏.工程流體力學[M].第2版.北京：高等教育出版社，2004，1：360-361.

[8]隋洪濤，李鵬飛，馬世虎，馬富銀，胡穎.精通CFD動網格工程仿真與案例實戰[M].第1版.北京：人民郵電出版社，2013，05：92-94.

[9]紀兵兵，陳金瓶.ANSYS ICEM CFD網格劃分技術實例詳解[M].第1版.北京：中國水利水電出版社，2012，01：255-259.

[10]吳逸飛，鄒正平.考慮來流邊界條件不確定性的數值模擬方法[J].航空發動機，2015，02：35-39.

[11]康順.計算域對CFD模擬結果的影響[J].工程熱物理學報，2005，S1：57-60.

[12]劉智益.不確定性CFD模擬方法及其應用研究[D].華北電力大學，2014.

劉明亮（1990—），男，碩士，主要研究方向：流體機械內部流場分析，E-mail：1052746541@qq.com

李勝男（1971—），女，碩士，高級工程師，主要研究方向：電力系統繼電保護等，E-mail：2450564588@qq.com

張思青（1957—），男，教授，E-mail：363685119@qq.com

Analysis of Mesh on the CFD Simulation Results

LIU Mingliang1，ZHANG Siqing1，LI Shengnan2
（1.Kunming University of Science and Technology，Kunming 650093，China；2.Yunnan Electric Power Research Institute，Kunming 650000，China）

In the CFD simulation process，the impact of highquality mesh on the calculation results is essential.This article describes in detail the process grid independence test.By Fluent double precision solver and k-omega model，the paper compares the calculation error of add and without boundary layer mesh，and then，Investigates the effect of the boundary layer mesh on calculation accuracy.

grid independent solution; grid independence test;k-omega models; boundary layer meshing

國家自然科學基金資助（編號：51369012）。