1、15.1 RNG湍流模型在22000雷諾數(shù)二維方柱繞流中的應(yīng)用【算例內(nèi)容】：方柱尺寸為100mm×100mm，流場(chǎng)區(qū)域?yàn)?400mm×3500mm，方柱距離入口500mm。流體工質(zhì)采用水，其密度為998.2千克每立方米，動(dòng)力粘度為0.001003帕秒。由雷諾數(shù)為22000可以折算出入口流體的速度應(yīng)為0.2210579米每秒。實(shí)驗(yàn)所用來(lái)流湍流強(qiáng)度為2%，取入口湍流長(zhǎng)度尺度為0.007（入口湍流長(zhǎng)度尺度取值參考文獻(xiàn)十六P215）?！痉抡婺繕?biāo)】：對(duì)于雷諾數(shù)為22000的亞臨界方柱繞流，合理選擇湍流模型將是關(guān)鍵。對(duì)于分離流動(dòng)，不適合采用標(biāo)準(zhǔn)k-湍流模型，在高雷諾數(shù)湍流模型族中，RN

2、G湍流模型是值得研究的選擇。而SST k-模型、雷諾應(yīng)力模型、分離渦模型和大渦模擬也是適合的，本節(jié)僅討論RNG湍流模型。要求用兩種近壁面處理方式處理方柱表面，第一種方式是采用尺度化壁面函數(shù)，第二種是采用加強(qiáng)型近壁面處理方法。計(jì)算完成后比較兩種近壁面處理方法的結(jié)果。文獻(xiàn)十七中給出了雷諾數(shù)為22000的方柱繞流的實(shí)驗(yàn)結(jié)果、大渦模擬結(jié)果和RNG計(jì)算結(jié)果。可以用來(lái)驗(yàn)證本節(jié)計(jì)算結(jié)果。關(guān)于雷諾數(shù)22000方柱繞流（關(guān)于壁面函數(shù)法、湍流粘度比的計(jì)算等概念可以參考文獻(xiàn)十六）15.1.1 RNG模型與尺度化壁面函數(shù)在方柱繞流中的應(yīng)用圖15-1-1 分析流程圖15-1-1給出了二維方柱繞流問(wèn)題的分析流程，在Wor

3、kbench中建立模型，然后使用Mesh模塊進(jìn)行分網(wǎng)，最后提交到Fluent中進(jìn)行計(jì)算。圖15-1-2 流場(chǎng)模型按照算例內(nèi)容建立流場(chǎng)模型如圖15-1-2所示，然后進(jìn)行分網(wǎng)，由于采用的是尺度化壁面函數(shù)法，所以近壁面網(wǎng)格沒(méi)必要?jiǎng)澐值靥?xì)。進(jìn)入到Mesh模塊后首先到mesh模型樹(shù)下將Physics Preference改為CFD與Fluent，然后設(shè)置Sizing下的網(wǎng)格相關(guān)度、網(wǎng)格光順等參數(shù)。圖15-1-3 網(wǎng)格劃分結(jié)果圖15-1-3給出了網(wǎng)格劃分結(jié)果，在劃分網(wǎng)格時(shí)，約束整個(gè)面的網(wǎng)格尺寸為10mm，并且定義入口、出口、上下邊界、方柱表面的名稱集（Named Selections），以方便在Flue

4、nt中定義邊界條件。圖15-1-4 局部網(wǎng)格示意圖15-1-4給出了局部網(wǎng)格示意圖，從圖中可以看出，在方柱表面的網(wǎng)格高度是較高的，使用尺度化近壁面函數(shù)法時(shí)并沒(méi)有劃分很細(xì)的邊界層網(wǎng)格。最后共得到48900個(gè)網(wǎng)格，網(wǎng)格最大長(zhǎng)寬比為1.0069，滿足仿真計(jì)算的網(wǎng)格要求。圖15-1-5 采用4核并行計(jì)算在打開(kāi)Fluent時(shí)，按照?qǐng)D15-1-5所示，采用單精度、4核并行計(jì)算。雙精度一般用在以下幾個(gè)方面：（1）存在明顯的細(xì)長(zhǎng)比很大的幾何特征，如細(xì)長(zhǎng)的管道；（2）流體是幾個(gè)區(qū)域，區(qū)域之間用小尺寸的管道連接，其中一個(gè)區(qū)域的流體壓力大大高于整個(gè)流域的平均水平；（3）對(duì)于存在較高導(dǎo)熱率的共軛傳熱問(wèn)題；（4）對(duì)于流

5、場(chǎng)網(wǎng)格長(zhǎng)寬比較大的問(wèn)題。在本算例中，雖然不存在以上問(wèn)題，但仍采用雙精度4核并行計(jì)算求解。在進(jìn)行計(jì)算時(shí)，先進(jìn)行穩(wěn)態(tài)計(jì)算，然后將穩(wěn)態(tài)計(jì)算的結(jié)果作為瞬態(tài)計(jì)算的初場(chǎng)，再進(jìn)行瞬態(tài)計(jì)算。求解器采用壓力基求解器，并進(jìn)行湍流模型的設(shè)置。圖15-1-6 湍流模型的選擇如圖15-1-6所示，選擇RNG湍流模型與尺度化壁面函數(shù)處理方法，并且不啟用Differential Viscosity Model選項(xiàng)，在ANSYS Fluent 14.5的幫助文件中對(duì)該模式有以下論述：The RNG turbulence model in ANSYS FLUENT has an option of using a d

6、ifferential formula for the effective viscosity to account for the low-Reynolds-number effects.可見(jiàn)在考慮到邊界層的低雷諾數(shù)影響時(shí)應(yīng)該打開(kāi)該選項(xiàng)，另外，有人在做空氣的渦街仿真時(shí)遇到過(guò)不勾選該選項(xiàng)就得不到渦街的情況。（但是在實(shí)際計(jì)算中發(fā)現(xiàn)，啟用該選項(xiàng)后，計(jì)算容易變得不穩(wěn)定，出現(xiàn)了數(shù)值波動(dòng)甚至是不收斂。因此，在本算例中，不啟用該選項(xiàng)。流線曲率修正是因?yàn)槭褂酶飨蛲詼u粘假設(shè)的湍流模型對(duì)流線曲率不敏感，因此在流線彎曲時(shí)需要進(jìn)行流線曲率修正，啟用該選項(xiàng)。選擇Fluent中的water-liquid作為流體材料，

7、然后進(jìn)行邊界條件的設(shè)置：設(shè)置方柱表面為無(wú)滑移壁面條件wall；流場(chǎng)上下邊界為對(duì)稱邊界條件symmetry（在數(shù)學(xué)上即自由滑移邊界條件）；當(dāng)流場(chǎng)發(fā)展到出口位置時(shí)，認(rèn)為流動(dòng)已經(jīng)充分發(fā)展，因此設(shè)置出口為自由流outflow；入口在用速度入口邊界條件，入口速度、湍流強(qiáng)度、湍流長(zhǎng)度尺寸均采用算例中所給設(shè)定。在設(shè)置參考值時(shí)首先選定compute from->流場(chǎng)名，這樣就正確的給出了參考值中的物理參數(shù)（密度、粘度等），然后指定Area的值為0.1。關(guān)于參考值的設(shè)定是非常有講究的，根據(jù)文獻(xiàn)十八P173的升阻力系數(shù)定義式：Cd=Fd12v2SCl=Fl12v2S以上兩式中，Cd代表阻力系數(shù)，Cl代表升力

8、系數(shù)，F(xiàn)d為鈍體所受到的阻力，F(xiàn)l為鈍體所受到的升力，為流體的密度，v為前方來(lái)流的速度，S代表鈍體的投影面積，在二維問(wèn)題中S的數(shù)值與弦長(zhǎng)相等（弦長(zhǎng)乘以一，二維問(wèn)題中展向長(zhǎng)度為單位長(zhǎng)度一）。有關(guān)翼型的概念：翼型：飛機(jī)機(jī)翼具有獨(dú)特的剖面，其橫斷面（橫向剖面）的形狀稱為翼型，稱為翼型前緣：翼型最前面的一點(diǎn)。 后緣：翼型最后面的一點(diǎn)。 翼弦：前緣與后緣的連線。 弦長(zhǎng)：前后緣的距離稱為弦長(zhǎng)。如果機(jī)翼平面形狀不是長(zhǎng)方形，一般在參數(shù)計(jì)算時(shí)采用制造商指定位置的弦長(zhǎng)或平均弦長(zhǎng)攻角：也稱迎角。對(duì)于翼形來(lái)說(shuō)，攻角定義為翼弦與來(lái)流速度之間的夾角，抬頭為正，低頭為負(fù)，常用符號(hào)表示。其示意圖見(jiàn)圖15-1-7：圖15-1-

9、7 機(jī)翼攻角示意圖升力系數(shù)：物體（如飛機(jī)、導(dǎo)彈）所受到的升力與氣流動(dòng)壓和參考面積的乘積之比，是一個(gè)無(wú)量綱量。阻力系數(shù)：對(duì)于飛行器來(lái)說(shuō)，阻力系數(shù)定義為物體（如飛機(jī)、導(dǎo)彈）所受到的阻力與氣流動(dòng)壓和參考面積之比，是一個(gè)無(wú)量綱量。關(guān)于升力系數(shù)與阻力系數(shù)的計(jì)算參考圖15-1-8：圖15-1-8 計(jì)算升阻力系數(shù)所用的參考系需要注意的是，在計(jì)算升力與阻力時(shí)，升力方向垂直于來(lái)流方向（不是垂直于機(jī)翼！），阻力方向平行于來(lái)流方向（不是平行于機(jī)翼?。?。僅當(dāng)攻角為0°或180°時(shí)，升力/阻力才垂直/平行于來(lái)流方向。在Fluent二維計(jì)算中，機(jī)翼長(zhǎng)度方向（展向）垂直于幾何平面，軟件在處理時(shí)機(jī)翼長(zhǎng)度置

10、為1（在流體力學(xué)中處理二維問(wèn)題，展向長(zhǎng)度全部取單位長(zhǎng)度）。其余參考值需要參考升力系數(shù)公式進(jìn)行計(jì)算，并且一定要注意升力/阻力的方向。對(duì)于本二維算例，S的取值與弦長(zhǎng)數(shù)值相等，因此本算例S取值為0.1，即設(shè)定參考值時(shí)，與尺寸相關(guān)的進(jìn)設(shè)置Area為0.1平方米即可，其余的保持默認(rèn)不變；密度等材料的參考值已經(jīng)通過(guò)Compute from功能實(shí)現(xiàn)了。之后設(shè)置求解器選項(xiàng)，選擇耦合求解器，對(duì)流項(xiàng)的處理采用二階迎風(fēng)格式（對(duì)于之前劃分的網(wǎng)格，也可采用更高階的QUICK格式）。保持默認(rèn)的收斂殘差判定標(biāo)準(zhǔn)0.001不變，監(jiān)測(cè)了升阻力系數(shù)隨迭代步數(shù)的變化曲線。以入口條件初始化整個(gè)流場(chǎng)區(qū)域，設(shè)定最大迭代步數(shù)為1000步。

11、圖15-1-9 穩(wěn)態(tài)計(jì)算的殘差曲線圖15-1-9給出了穩(wěn)態(tài)計(jì)算殘差曲線，從曲線中可以看出，殘差曲線基本上隨著計(jì)算會(huì)繼續(xù)趨于收斂，其中連續(xù)性方程的殘差最大。在計(jì)算此類繞流問(wèn)題的實(shí)際應(yīng)用過(guò)程中，還會(huì)出現(xiàn)當(dāng)穩(wěn)態(tài)計(jì)算的殘差位于高位穩(wěn)定的情況，出現(xiàn)這種情況時(shí)，應(yīng)當(dāng)考慮的一種可能存在的情況就是物理模型中存在瞬態(tài)變化的物理效應(yīng)，導(dǎo)致穩(wěn)態(tài)迭代時(shí)前后迭代的范數(shù)殘差不能夠收斂而位于高位穩(wěn)定，這種現(xiàn)象在用穩(wěn)態(tài)雷諾平均方法求解瞬態(tài)問(wèn)題時(shí)經(jīng)常出現(xiàn)，此時(shí)可以將穩(wěn)態(tài)計(jì)算結(jié)果作為瞬態(tài)計(jì)算的初場(chǎng)，繼續(xù)計(jì)算。所以單看殘差收斂情況不一定有太大意義，要結(jié)合實(shí)際情況來(lái)看。另外發(fā)現(xiàn)當(dāng)啟用Differential Viscosity Mo

12、del選項(xiàng)后，本來(lái)收斂至非常小的殘差容易轉(zhuǎn)變?yōu)楦呶环€(wěn)定甚至發(fā)散。本例中升力系數(shù)與阻力系數(shù)的迭代歷史曲線如圖15-1-10所示：圖15-1-10 方柱繞流的穩(wěn)態(tài)升阻力迭代歷史從圖15-1-10中可以看出，升阻力系數(shù)在穩(wěn)態(tài)計(jì)算中最終歸于穩(wěn)定，這與殘差收斂至極小值是相互呼應(yīng)的。而在有些繞流問(wèn)題的穩(wěn)態(tài)計(jì)算中，也會(huì)出現(xiàn)瞬態(tài)特征。如圖15-1-11所示： 圖15-1-11 某圓柱繞流的穩(wěn)態(tài)計(jì)算圖15-1-11中給出了某圓柱繞流的穩(wěn)態(tài)計(jì)算結(jié)果，在穩(wěn)態(tài)計(jì)算結(jié)果中出現(xiàn)了瞬態(tài)特征。在穩(wěn)態(tài)計(jì)算中出現(xiàn)瞬態(tài)特征的現(xiàn)象，在方柱繞流與圓柱繞流的現(xiàn)象中都可能出現(xiàn)，當(dāng)出現(xiàn)這一現(xiàn)象時(shí)，殘差曲線將會(huì)居于高位，并且升力與阻力系數(shù)的值

13、不一定準(zhǔn)確，比如上圖中升力系數(shù)平均值部位0。還有一種現(xiàn)象是升阻力系數(shù)趨于平穩(wěn)，但殘差曲線局于高位，造成這種現(xiàn)象的一種原因是啟用了Differential Viscosity Model選項(xiàng)，這也說(shuō)明了Differential Viscosity Model選項(xiàng)對(duì)計(jì)算所造成的不穩(wěn)定性。根據(jù)文獻(xiàn)十八P119的關(guān)于斯特勞哈爾數(shù)的公式：St=fDv其中，St為斯特勞哈爾數(shù)，f為漩渦脫落頻率，v為來(lái)流速度。根據(jù)文獻(xiàn)十七中的實(shí)驗(yàn)結(jié)果（St=0.13），可以反推出漩渦脫落頻率應(yīng)為3.5s，所以在選擇瞬態(tài)步長(zhǎng)時(shí)就有了選擇的依據(jù)，瞬態(tài)求解采用二階隱式求解方式，針對(duì)方柱繞流問(wèn)題，隱式方法對(duì)時(shí)間步長(zhǎng)的要求不高，很多

14、人取漩渦脫落周期的1/50作為瞬態(tài)計(jì)算的時(shí)間步長(zhǎng)，本文取為0.02s，可以充分地描述出漩渦脫落的瞬態(tài)細(xì)節(jié)。對(duì)于收斂殘差的取值，很多人是取穩(wěn)態(tài)殘差穩(wěn)定值的1/50，本文取為1/100，即0.00001。指定每一時(shí)間步最多計(jì)算25個(gè)子步，共計(jì)算5000個(gè)時(shí)間步，共計(jì)100s。圖15-1-12 方柱繞流流場(chǎng)的靜壓云圖與動(dòng)壓運(yùn)圖根據(jù)圖15-1-12中可以看出，在方柱后方出現(xiàn)了明顯的周期性漩渦脫落的現(xiàn)象。漩渦脫落的現(xiàn)象在總壓云圖與速度云圖中也可以很明顯地看出來(lái)。文獻(xiàn)十七中認(rèn)為RNG方法在預(yù)測(cè)平均阻力系數(shù)和回流區(qū)再附著距離時(shí)與實(shí)驗(yàn)及大渦結(jié)果比較接近，但是在預(yù)測(cè)升阻力系數(shù)的均方根值時(shí)仍具有一定的差距。采用雷

15、諾平均的RNG模型，從模型精度上來(lái)說(shuō)要比大渦模擬至少低一個(gè)數(shù)量級(jí)。圖15-1-13 方柱繞流流場(chǎng)的總壓云圖與速度云圖圖15-1-14給出了方柱升阻力系數(shù)的歷史曲線。圖15-1-14 方柱繞流的升力與阻力系數(shù)的歷史曲線從圖中可以明顯看出，方柱繞流的升阻力系數(shù)隨著時(shí)間的推移呈周期性的波動(dòng)。經(jīng)過(guò)對(duì)數(shù)據(jù)的統(tǒng)計(jì)，可以知道，方柱阻力系數(shù)的平均值為2.3（取80-90s之間的每個(gè)時(shí)間步的阻力系數(shù)做平均），實(shí)驗(yàn)值為2.1，誤差為9.5%，計(jì)算出的漩渦脫落周期為3s，由此所得的斯特勞哈爾數(shù)為0.15，實(shí)驗(yàn)值為0.13，誤差為15.4%。圖15-1-15 方柱表面的Y+值分布范圍仿真后可以知道，近壁面Y+值的分布

16、范圍為7.07-68.6，在Fluent的官方文獻(xiàn)中推薦壁面函數(shù)法的Y+值適用范圍為30-500，對(duì)于尺度化壁面函數(shù)，由于在近壁面使用了預(yù)防Y+值過(guò)小導(dǎo)致結(jié)果惡化的技術(shù)，其Y+值范圍下限可以適當(dāng)放寬，但仍應(yīng)盡量取在11.25以上。過(guò)大或過(guò)小的Y+值都可能引起結(jié)果的惡化。15.1.2 RNG模型與加強(qiáng)型近壁面處理在方柱繞流中的應(yīng)用圖15-1-16 加強(qiáng)型近壁面處理所用到的流場(chǎng)網(wǎng)格返回到MESH模塊中，將方柱邊的尺寸約束為1mm，在方柱的邊上添加邊界層網(wǎng)格控制，指定第一層網(wǎng)格高度為0.25mm，網(wǎng)格增長(zhǎng)率為1.05，共劃分了15層邊界層網(wǎng)格，得到了135500個(gè)節(jié)點(diǎn)，共計(jì)68053個(gè)單元。所有單元

17、的最大長(zhǎng)寬比為4.53，完全滿足加強(qiáng)型近壁面處理的網(wǎng)格計(jì)算要求。圖15-1-17 加強(qiáng)型近壁面處理后的方柱升阻力系數(shù)歷史歷程通過(guò)圖15-1-17可以看出，使用了加強(qiáng)型近壁面處理方法之后，所得到的升阻力系數(shù)波動(dòng)更加穩(wěn)定直觀。將升阻力系數(shù)的保存文件使用Matlab讀入并處理后，可以知道，加強(qiáng)型近壁面處理計(jì)算出的波動(dòng)穩(wěn)定后的阻力系數(shù)平均值為2.15，與實(shí)驗(yàn)值2.10的誤差為2.38%，準(zhǔn)確度明顯高于尺度化壁面函數(shù)法；計(jì)算所得的斯特勞哈爾數(shù)為0.14，與實(shí)驗(yàn)值0.13的誤差為7.69%，其準(zhǔn)確度也明顯高于壁面函數(shù)法。計(jì)算所得的阻力系數(shù)波動(dòng)的均方根值為0.125，而W.Rodi等人通過(guò)大渦模擬得到的阻力

18、系數(shù)波動(dòng)的均方根值為0.14，誤差為-11.4%（文獻(xiàn)十七使用RNG模型計(jì)算的結(jié)果為0.11）；計(jì)算所得的升力系數(shù)波動(dòng)的均方根值為1.41，W.Rodi使用大渦計(jì)算的結(jié)果為1.01，誤差為39.6%（文獻(xiàn)十七使用RNG模型計(jì)算的結(jié)果為1.26）。通過(guò)以上結(jié)果可以看出，RNG模型與加強(qiáng)型近壁面處理方法的結(jié)合使用，能夠較好地預(yù)測(cè)阻力系數(shù)平均值、漩渦的脫落周期與阻力系數(shù)波動(dòng)的均方根值，但是在預(yù)測(cè)升力系數(shù)波動(dòng)的均方根值時(shí)，還有較大的誤差，與大渦結(jié)果存在一定的差距。（注意，升阻力系數(shù)波動(dòng)的結(jié)果均是與大渦模擬比較，在查閱資料時(shí)并未發(fā)現(xiàn)相關(guān)實(shí)驗(yàn)數(shù)據(jù)，因此在評(píng)判RNG計(jì)算所得的升阻力系數(shù)波動(dòng)時(shí)，均是與大渦模擬

19、結(jié)果相比較）同時(shí)必須檢查Y+值是否處于合理的范圍之內(nèi)，通過(guò)Contour中的Turbulence->Y Plus來(lái)顯示方柱表面上的Y+值分布范圍，經(jīng)查看后可以知道，Y+值分布在0.197-0.547之間，穩(wěn)態(tài)計(jì)算結(jié)果中的Y+值比瞬態(tài)計(jì)算的Y+值略小，但是前后兩次計(jì)算的Y+對(duì)于RNG與加強(qiáng)型近壁面處理來(lái)說(shuō)都是可以接受的。對(duì)于靜壓、動(dòng)壓、總壓、速度云圖如下圖15-1-18所示：圖15-1-18（a） 靜壓與動(dòng)壓云圖圖15-1-18（b） 總壓與速度云圖附錄：簡(jiǎn)單Fluent升力、阻力、力矩系數(shù)后處理程序%簡(jiǎn)單Fluent升力、阻力、力矩系數(shù)后處理程序 徐太棟 2013clcclearformat longfid=fopen('S