基于35°Ahmed模型的外流場(chǎng)仿真精確性與 被動(dòng)控制減阻研究

田 思， 譚文林， 張英朝， 吳 敏， 邵書鑫

(1. 奇瑞汽車股份有限公司， 安徽 蕪湖 241000； 2.中國(guó)汽車工程研究院股份有限公司， 重慶 401120； 3.吉林大學(xué) 汽車仿真與控制國(guó)家重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)春 130022)

計(jì)算流體動(dòng)力學(xué)也被稱作CFD[1]，它是用計(jì)算機(jī)數(shù)值計(jì)算去研究流體力學(xué)的一門學(xué)科，也是研究汽車空氣動(dòng)力學(xué)的一種手段。計(jì)算機(jī)技術(shù)的進(jìn)步使得數(shù)值仿真的效率和仿真結(jié)果的可靠度都有了很大的提升[2]，盡管如此，要想得到地面車輛外流場(chǎng)的精確模擬依然不容易。車輛外流場(chǎng)的流動(dòng)結(jié)構(gòu)是很復(fù)雜的，而且氣體的流動(dòng)特征是不規(guī)則的，因此外流場(chǎng)仿真的精確度受很多因素的影響，主要有幾何模型的前處理及建模、網(wǎng)格劃分的類型、物理模型的選取、計(jì)算資源的供給和仿真策略的制定等[3]。在汽車CFD仿真分析中，我們主要研究汽車外流場(chǎng)，主要研究的對(duì)象為：應(yīng)用快捷的CFD仿真方法來(lái)優(yōu)化汽車的氣動(dòng)性能；為能夠達(dá)到汽車氣動(dòng)性能的驗(yàn)證性要求，需要開(kāi)發(fā)出準(zhǔn)確度高的CFD仿真策略[4]。汽車的氣動(dòng)阻力主要由摩擦阻力和壓差阻力組成[5]。湍流模型[6]、網(wǎng)格、求解器這3個(gè)影響因素對(duì)汽車外流場(chǎng)仿真精確性的影響至關(guān)重要。汽車尾部流場(chǎng)因氣流分離嚴(yán)重而產(chǎn)生較復(fù)雜的湍流運(yùn)動(dòng)[7]，對(duì)這部分流場(chǎng)的仿真有一定的難度，也占用較大的計(jì)算資源。目前，應(yīng)用較多的仿真方法有直接數(shù)值模擬、大渦模擬、雷諾時(shí)均法。不同的數(shù)值仿真方法對(duì)計(jì)算機(jī)硬件的要求也是不同的，在滿足工程需要的前提下，本文選用常用的雷諾時(shí)均法進(jìn)行仿真計(jì)算[8]，探索出一種質(zhì)量高、精度高的CFD[9]仿真策略。雷諾時(shí)均法所用的模型包括Reynolds應(yīng)力模型、標(biāo)準(zhǔn)k-e模型、Realizablek-e模型和SSTk-w模型等[9]。不同的仿真模型采用的求解方程不一樣，也影響了計(jì)算結(jié)果的精確度，但基本都能滿足工程需要。當(dāng)采用數(shù)值方法求解控制方程時(shí)，都是想辦法將控制方程在空間區(qū)域上進(jìn)行離散，然后求解得到的離散方程組，要想在空間域上離散控制方程，必須使用網(wǎng)格。網(wǎng)格是離散的基礎(chǔ)，網(wǎng)格節(jié)點(diǎn)是離散化后物理量的存儲(chǔ)位置，網(wǎng)格在離散過(guò)程中起著關(guān)鍵的作用。網(wǎng)格類型[10]和布置對(duì)網(wǎng)格質(zhì)量來(lái)說(shuō)是至關(guān)重要的，常用的網(wǎng)格主要有4類，包括常用的四面體非結(jié)構(gòu)網(wǎng)格和六面體剪切體網(wǎng)格，以及不常用的多面體非結(jié)構(gòu)網(wǎng)格和十二面體剪切體網(wǎng)格。仿真的精確度受網(wǎng)格類型、網(wǎng)格尺寸和網(wǎng)格數(shù)量的影響很大[11]。本文要找到仿真精確度高的仿真策略，就要同時(shí)考慮湍流模型和網(wǎng)格對(duì)計(jì)算結(jié)果的影響，通過(guò)對(duì)不同組合的研究，得到最佳組合方案。Ahmed模型是研究外流場(chǎng)結(jié)構(gòu)常用的一種類車體模型，它是由德國(guó)人Ahmed發(fā)明并提出的。通過(guò)試驗(yàn)，發(fā)現(xiàn)模型的阻力系數(shù)隨著它的尾部?jī)A角的變化而改變，且對(duì)應(yīng)著不同的尾流結(jié)構(gòu)。Ahmed等探究了該模型空氣阻力的來(lái)源，這為后人研究減阻技術(shù)提供了極大的幫助。后人對(duì)Ahmed模型進(jìn)行了大量的仿真與風(fēng)洞試驗(yàn)，有豐富的經(jīng)驗(yàn)數(shù)據(jù)可以參考。本文對(duì)汽車外流場(chǎng)仿真精確性的研究，采用35°Ahmed模型[12]，當(dāng)尾部?jī)A角為35°時(shí)，Ahmed模型阻力系數(shù)較低，且尾流結(jié)構(gòu)穩(wěn)定，適合進(jìn)一步研究它的減阻機(jī)理。因先忽略網(wǎng)格尺寸的影響，根據(jù)經(jīng)驗(yàn)選擇適中的網(wǎng)格尺寸，對(duì)湍流模型和網(wǎng)格類型這兩種影響因素進(jìn)行探討，選擇4種湍流模型和4種網(wǎng)格類型，通過(guò)組合，共有16個(gè)Case。對(duì)這些Case進(jìn)行初步的評(píng)價(jià)，找出效果較好的4個(gè)Case，選定這4種仿真方案，進(jìn)一步去進(jìn)行網(wǎng)格尺寸的研究。在這一階段，會(huì)對(duì)網(wǎng)格類型、網(wǎng)格尺寸與湍流模型這3種影響因素進(jìn)行最終評(píng)價(jià)，得到最優(yōu)仿真結(jié)果。本部分仿真研究共涉及32個(gè)Case，選定最優(yōu)仿真方案后，在35°Ahmed模型尾部斜面上端添加凸起結(jié)構(gòu)，通過(guò)改變凸起與模型背部的夾角去研究被動(dòng)控制技術(shù)對(duì)減阻的影響，并實(shí)現(xiàn)一定程度的減阻。

1 數(shù)值仿真方案的設(shè)定

由于空氣相對(duì)于汽車的流速低于0.3個(gè)馬赫數(shù)，在汽車外流場(chǎng)的仿真中，空氣可以被看作是理想的氣體[13]，具有不可壓縮性，我們采用離散控制方程，應(yīng)用有限體積法去求解。我們從控制方程的空間離散格式與松弛因子這兩個(gè)方面去設(shè)置求解器?？臻g離散格式對(duì)計(jì)算的穩(wěn)定性有較大影響，松弛因子對(duì)方程的求解速度有較大影響，這兩個(gè)因素都影響著計(jì)算結(jié)果的精確性。本文采用滿足二階精度的離散格式。

汽車外流場(chǎng)仿真精確性研究的具體方案為：第一階段，先忽略網(wǎng)格尺寸的影響，根據(jù)經(jīng)驗(yàn)選擇適中的網(wǎng)格尺寸，選擇4種湍流模型和4種網(wǎng)格類型作為研究對(duì)象，組合出16個(gè)Case，如表1所示，對(duì)這些Case進(jìn)行初步的評(píng)價(jià)，找出效果較好的4個(gè)Case；第二階段，在第一階段選定的4種方案的基礎(chǔ)上，探討網(wǎng)格尺寸對(duì)仿真精確度的影響，這一階段也會(huì)有16個(gè)Case，如表2所示。最后，會(huì)得出這3種影響因素的最終評(píng)估結(jié)果，得出最優(yōu)仿真方案。研究中為方便描述，將四面體網(wǎng)格用符號(hào)TE表示，六面體網(wǎng)格用符號(hào)T表示，多面體網(wǎng)格用符號(hào)P表示，十二面體網(wǎng)格用符號(hào)TD表示。將標(biāo)準(zhǔn)k-e湍流模型用符號(hào)KS表示，Realizablek-e模型用符號(hào)KR表示，SSTk-w模型用符號(hào)KW表示，Reynolds應(yīng)力模型用符號(hào)TL表示。將最小加密體的網(wǎng)格尺寸用數(shù)字表示，4表示最大值，為20 mm，3表示網(wǎng)格尺寸為15 mm,2表示網(wǎng)格尺寸為10 mm，1表示網(wǎng)格尺寸為5 mm。鑒于湍流模型的選取受壁面函數(shù)的影響較大，本文根據(jù)經(jīng)驗(yàn)將壁面函數(shù)的值設(shè)置在允許的范圍內(nèi)。

表1 第一階段仿真方案設(shè)定

表2 第二階段仿真方案設(shè)定

本仿真采用具有35°后窗傾角的Ahmed模型，總的風(fēng)阻系數(shù)Cw由壓差阻力系數(shù)Cp和摩擦阻力系數(shù)Cr組成，其中，壓差阻力系數(shù)包含了前端的壓差阻力系數(shù)Ck、尾部斜面的壓差阻力系數(shù)Cs和尾部垂直面的壓差阻力系數(shù)Cb。本文采用的數(shù)字風(fēng)洞為15 m×3 m×1.5 m的長(zhǎng)方體計(jì)算域，模型的阻塞比不超過(guò)1%，基本滿足了要求的阻塞比。模型到入口的距離約為模型長(zhǎng)度的4倍，模型到出口的距離約為模型長(zhǎng)度的8倍。這樣的設(shè)置使模型周圍有足夠大的空間，使它附近特別是尾部的空氣流動(dòng)能夠充分發(fā)展。為保證與風(fēng)洞試驗(yàn)的數(shù)據(jù)保持一致，將模型底部與地面的距離設(shè)置為50 mm[5]，設(shè)定60 m/s的來(lái)流速度，雷諾數(shù)Re為4.29×106，湍流[14]強(qiáng)度小于0.5%。采用分離式的求解器，適中的松弛因子和二階精度的空間離散格式。

如圖1模型設(shè)定兩重加密體，為進(jìn)行網(wǎng)格尺寸的研究[15]，把第一重加密體的網(wǎng)格尺寸設(shè)置成可變化的；采用5層棱柱層網(wǎng)格對(duì)邊界層進(jìn)行模擬。圖2為邊界層網(wǎng)格示意圖，圖3為模型表面的Y+云圖，可見(jiàn)模型表面的Y+值基本控制在20～100之間，這保證了對(duì)邊界層流動(dòng)的仿真更準(zhǔn)確[16]。

圖1 幾何模型與網(wǎng)格空間整體布置圖

圖2 邊界層網(wǎng)格示意圖

圖3 模型表面Y+云圖

計(jì)算模型的物理?xiàng)l件及求解參數(shù)設(shè)置如表3所示。

表3 仿真相關(guān)參數(shù)設(shè)置

2 仿真精確性研究

工程上的誤差一般允許在5%左右，為得到更為精確的仿真策略，使仿真結(jié)果的誤差不超過(guò)1%，本文依次從整體阻力系數(shù)、局部阻力系數(shù)、模型尾部流場(chǎng)的物理量具體數(shù)據(jù)這3個(gè)層級(jí)進(jìn)行分析評(píng)價(jià)，使結(jié)論更有說(shuō)服力。

2.1 第一步仿真分析

此階段仿真采用適中的網(wǎng)格尺寸，這樣可以節(jié)約時(shí)間，快速地獲得計(jì)算結(jié)果，探討湍流模型和網(wǎng)格類型對(duì)仿真精確度的影響。在這一步分析中，主要考慮整體風(fēng)阻系數(shù)和局部風(fēng)阻系數(shù)這兩個(gè)評(píng)價(jià)指標(biāo)，且以整體風(fēng)阻系數(shù)為主。結(jié)果如表4所示。其中的試驗(yàn)值參考文章Ahmed et al.1984[5]；表4中Cp是總的壓差阻力系數(shù)，仿真的結(jié)果均精確到小數(shù)點(diǎn)后第3位即千分位；由于T+TL和TE+TL兩種情況下的計(jì)算結(jié)果是發(fā)散的，不予考慮。

表4 初步仿真結(jié)果

從表4中的仿真結(jié)果可知，各個(gè)工況下Cr的值基本一樣，即模型總摩擦阻力系數(shù)基本不變，因此，該模型的主要阻力為它的壓差阻力。在此，忽略模型摩擦阻力系數(shù)對(duì)阻力的影響，主要研究它的壓差阻力系數(shù)。

圖4是總的壓差阻力系數(shù)Cp的折線圖，試驗(yàn)數(shù)據(jù)在圖中用紅色水平直線表示，4種網(wǎng)格和4種湍流模型的仿真數(shù)據(jù)如圖所示。

圖4 第一步仿真的總壓差阻力系數(shù)

圖5是在第一步仿真中，采用適中的同一尺寸的網(wǎng)格時(shí)，4種類型體網(wǎng)格模型所生成的網(wǎng)格數(shù)量對(duì)比圖。對(duì)網(wǎng)格數(shù)量進(jìn)行歸一化處理，其中以1表示多面體網(wǎng)格P的數(shù)量，歸一化處理后，其余類型網(wǎng)格的數(shù)量如圖5所示。

圖5 第一步仿真體網(wǎng)格數(shù)量對(duì)比

經(jīng)過(guò)分析研究，可以排除P和TD這兩種網(wǎng)格，以及KS和TL這兩種湍流模型。下面將對(duì)由SSTk-w、Realizablek-e這兩種湍流與四面體網(wǎng)格、六面體網(wǎng)格這兩種網(wǎng)格類型組合的工況做進(jìn)一步的研究。

2.2 第二步仿真分析

在第一步仿真分析中，已經(jīng)確定了幾種有效的仿真方案，接下來(lái)還需對(duì)這些方案的網(wǎng)格進(jìn)行精細(xì)處理。通常局部較復(fù)雜的流場(chǎng)對(duì)網(wǎng)格的精細(xì)度要求比較高[17-18]，網(wǎng)格過(guò)大不易捕捉到流場(chǎng)的細(xì)節(jié)信息，因此我們選擇局部阻力系數(shù)作為第二級(jí)評(píng)價(jià)指標(biāo)，從而對(duì)仿真的精確度做進(jìn)一步的研究。

表5為第二步仿真的結(jié)果，表中的所有數(shù)據(jù)精確到千分位。從表5中，可以看出一個(gè)很明顯的特征，即仿真中的摩擦阻力值和試驗(yàn)值也相差10個(gè)count左右，這與之前第一步仿真中得到的值一樣都相對(duì)一致。因此，在這步的研究中，我們同樣忽略摩擦阻力的影響，將壓差阻力作為主要研究對(duì)象。

表5 第二步仿真結(jié)果

圖6縱坐標(biāo)Cp表示模型總的壓差阻力系數(shù)的值，橫坐標(biāo)代表網(wǎng)格的大小。網(wǎng)格尺寸經(jīng)過(guò)歸一化處理后如圖所示，紅線表示模型在風(fēng)洞中測(cè)得的試驗(yàn)值。從圖6中可以明顯地看出，總體上Cp值隨著網(wǎng)格尺寸的減小而逐漸減小，從大于試驗(yàn)值降到等于試驗(yàn)值，再到小于試驗(yàn)值。對(duì)比分析可以看出，仿真計(jì)算的精確度與網(wǎng)格設(shè)置的大小有很大關(guān)系，網(wǎng)格大小的設(shè)置應(yīng)該限制在一定量級(jí)，采用適中的網(wǎng)格尺寸。從圖6可知，這些仿真工況中精確度最高的工況為Case28，即T+KR+1組合，其次為Case19和Case23,即T+KW+2和TE+KW+2這兩個(gè)組合，下面將對(duì)這幾個(gè)工況的局部阻力系數(shù)作深入的分析。

圖6 第二步仿真總壓差阻力系數(shù)

圖7為選定3個(gè)工況與試驗(yàn)值對(duì)比的Cp值條形圖?？梢郧逦乜吹礁鱾€(gè)工況的總壓差阻力系數(shù)與試驗(yàn)值之間的誤差關(guān)系，由該圖所得的結(jié)論與前面所作的分析結(jié)果是一致的。

圖7 3種組合Cp仿真值與試驗(yàn)值比較

圖8為3種組合工況中Ahmed模型的前端、尾部斜面、尾部垂直面的局部風(fēng)阻系數(shù)與試驗(yàn)值對(duì)比。從圖8可以清晰地看出，在局部部位的風(fēng)阻系數(shù)中，方案TE+KW+2總是與試驗(yàn)值最為接近。而前面從圖6中得到的最優(yōu)組合T+KR+1的誤差之所以最小，由圖8可知這是在Ck、Cb值的過(guò)小與Cs值的過(guò)大雙重作用下形成的，這是一種假象，不是理想的結(jié)果。

因此，我們得到仿真精度最高的方案是TE+KW+2組合，這個(gè)方案才是最優(yōu)結(jié)果。其次較好的方案為T+KW+2組合，再次為T+KR+1組合。

2.3 局部流場(chǎng)物理量具體數(shù)據(jù)分析

為了進(jìn)一步驗(yàn)證第二步分析得到的3種組合的仿真精度，下面我們從局部流場(chǎng)的數(shù)據(jù)即分離渦剪切層附近的流速來(lái)進(jìn)行分析。

圖8 局部壓差阻力系數(shù)對(duì)比

如圖9所示，在模型尾部斜面的對(duì)稱面上選取4條豎直線，選定直線上特定位置的點(diǎn)，將仿真得到的這些點(diǎn)的速度值與已有的風(fēng)洞試驗(yàn)值作對(duì)比，進(jìn)一步驗(yàn)證前面所得的結(jié)果。

圖9 數(shù)據(jù)選取的位置

圖10為選定4條直線上各個(gè)方案仿真速度值與試驗(yàn)值的對(duì)比，從中我們可以明顯地看到：

1)選定的3個(gè)工況與試驗(yàn)數(shù)據(jù)相比，均沒(méi)有與試驗(yàn)值完全吻合，均存在一定的誤差。

2) 在這3個(gè)工況中，工況TE+KW+2與試驗(yàn)值吻合最好，誤差程度最小，可以確定為仿真精確度最好的方案，這與前面所得結(jié)論是一致的。

綜上可知，本文所得的最優(yōu)仿真方案為TE+KW+2這個(gè)工況，該方案應(yīng)用四面體網(wǎng)格。由于四面體網(wǎng)格生成的網(wǎng)格數(shù)量要比六面體網(wǎng)格多，相應(yīng)的計(jì)算量就會(huì)加大，在工程應(yīng)用中，需要考慮計(jì)算時(shí)間、成本等因素，因此可以根據(jù)具體情況選擇合理的方案進(jìn)行仿真。

圖10 剪切層處的氣流縱向速度變化對(duì)比圖

3 被動(dòng)控制減阻研究

3.1 減阻方案設(shè)定

減阻方案如圖11所示。

圖11 減阻方案示意圖

參考第2部分的研究結(jié)果，本文采用精確仿真的方法，在35°Ahmed模型尾部斜面上端添加厚度為1 mm，高度為20 mm的凸起結(jié)構(gòu)，如圖11所示，以凸起與尾部斜面的夾角為變量，從0°開(kāi)始以10°為間隔遞增，找到模型的最優(yōu)減阻工況。

3.2 仿真結(jié)果

由于Ahmed模型的風(fēng)阻主要來(lái)自于壓差阻力，下面的研究將以Ahmed模型的壓差阻力系數(shù)為評(píng)價(jià)指標(biāo)進(jìn)行研究。Ahmed模型被動(dòng)減阻仿真結(jié)果如表6所示。

表6 被動(dòng)減阻方案計(jì)算結(jié)果

從表6可以看出，隨著夾角的增大，模型的壓差阻力系數(shù)的變化趨勢(shì)是先減小再增大，然后再減小。當(dāng)附加板與斜面角度為30°時(shí)，模型的壓差阻力系數(shù)最小，此時(shí)的值為0.190，相對(duì)于原始數(shù)模的壓差阻力系數(shù)減小了11個(gè)count，有較好的減阻效果。

3.3 后處理分析

選出減阻效果最好的工況，即附加板與斜面角度為30°的工況作深入的分析，探索Ahmed模型的減阻機(jī)理。由之前的分析我們知道，Ahmed模型的阻力主要來(lái)源于模型的尾部，因此重點(diǎn)關(guān)注模型的尾部斜面和豎直面區(qū)域。表7為最優(yōu)減阻方案與原始方案的局部壓差阻力系數(shù)對(duì)比表。從中可以看出模型尾部斜面的壓差阻力系數(shù)變化更明顯。

表7 標(biāo)準(zhǔn)模型與最優(yōu)減阻方案結(jié)果對(duì)比

圖12為標(biāo)準(zhǔn)35°Ahmed模型與最優(yōu)減阻工況的尾部斜面壓力云圖[17]。最優(yōu)減阻工況斜面上的壓力相較于標(biāo)準(zhǔn)模型有所增大，這就有效地減小了模型的壓差阻力，而且C柱渦在斜面上的印記變的很弱。

圖12 尾部斜面壓力云圖

圖13為標(biāo)準(zhǔn)35°Ahmed模型與最優(yōu)減阻工況的縱向?qū)ΨQ面上的速度矢量圖，優(yōu)化模型的尾部分離渦相較于原始模型變化不明顯。

圖13 縱向?qū)ΨQ面上速度矢量圖

圖14為標(biāo)準(zhǔn)35°Ahmed模型與最優(yōu)減阻工況的尾部10 mm處橫截面上的速度云圖。從中可以看出減阻方案尾部斜面上的分離區(qū)域有所擴(kuò)展，與C柱渦融合得更強(qiáng)些，C柱渦的強(qiáng)度明顯減弱。

圖14 尾部10mm處橫截面上速度云圖

4 結(jié)論

本文采用數(shù)值模擬的方法，從網(wǎng)格類型、湍流模型、網(wǎng)格尺寸這3個(gè)影響因素出發(fā)，分3個(gè)評(píng)價(jià)層次對(duì)汽車外流場(chǎng)的仿真精確性進(jìn)行了研究，找到了精確度高的仿真策略，在保證仿真精確度的基礎(chǔ)上，對(duì)Ahmed模型被動(dòng)控制減阻技術(shù)作了研究，可以得到如下結(jié)論：

1)數(shù)值模擬的精確性受網(wǎng)格類型、湍流模型和網(wǎng)格尺寸等因素的影響比較大，在進(jìn)行仿真計(jì)算前，有必要對(duì)仿真精確度進(jìn)行研究，確保計(jì)算結(jié)果的可靠性。

2)對(duì)于35°Ahmed模型，它的阻力主要來(lái)源于模型尾部。本文研究了模型尾部添加凸起結(jié)構(gòu)對(duì)減阻的影響，當(dāng)附加板與斜面角度為30°時(shí)，有最優(yōu)減阻效果。

3)凸起結(jié)構(gòu)的作用主要是改善了35°Ahmed模型周圍流場(chǎng)的分布，從而有效地減小了模型的壓差阻力。