類車體氣動性能的大渦模擬

朱 暉， 周永祥, 楊志剛, 史芳琳

(1. 同濟大學 上海地面交通工具風洞中心， 上海 201804； 2. 同濟大學 上海市地面交通工具空氣動力與熱環(huán)境模擬重點實驗室, 上海 201804； 3. 泛亞汽車技術中心有限公司， 上海 201201)

大渦模擬(LES)最早于1963年由Smagorinsky以亞格子應力模型的方式運用于氣象學研究中[1].而后LES被逐漸應用于平面射流、繞流、燃燒、槽道流和氣固兩相流等多個方面，并取得了一定的研究成果.

在針對二維流動的數值仿真中，馬金花等論證了LES方法對二維圓柱繞流計算的適用性[2];苑明順等采用LES方法研究了二維圓柱繞流中三維渦量的耗散效應[3]；Murakami等及Yu等的研究表明，LES對方柱繞流的預測結果更接近實際[4-5]，計算結果受網格尺度及亞格子模型構建方式的影響顯著[6]；蘇銘德的研究表明，LES對方直管內流中二次流現象的捕捉及統(tǒng)計平均量的解析與試驗符合良好[7];楊建明等采用LES方法對方彎管內流動的仿真結果與實驗數據亦符合良好[8]；劉沛清等利用LES成功預測了某翼型在不同迎角下的外流場特性，并得出相應的氣動力分量[9]；劉寧宇等采用動力學亞格子模型研究了圓柱繞流在穩(wěn)態(tài)和振蕩下的尾跡流態(tài)[10]； Cheng等的研究表明，與k-ε模型相比，LES方法針對方柱繞流中所涉及的分離流及回流的預測結果更為準確[11]；王兵等以自主開發(fā)的亞格子模型為基礎對后臺階流的再附著過程展開研究，其計算數據與試驗結果較為一致[12].

隨著計算機技術的發(fā)展，以及對LES濾波方法和亞格子模型研究的深入，LES被逐步應用于求解三維湍流流動.Constantinescu等的研究發(fā)現，雖然LES和時均湍流模型(RANS)皆能準確預測出鈍體繞流場的平均速度分布，但LES能更好地預測湍流動能[13]；Sinisa等利用LES方法對斜背傾角為25°的Ahmed類車體外部繞流場進行數值仿真，采用不同的網格方案對氣動阻力及類車體頭部、上方、斜背和尾跡的流場結構進行了詳細的對比分析[13-15];Eric等的研究表明大渦模擬可以得到與試驗數據總體符合的結果，但Smagorinsky亞格子湍流模型不能準確預測出斜背上的分離[16]；Aljure等以Ahmed類車體和Asmo 類車體為對象，采用4種亞格子湍流模型對其繞流場進行數值仿真，通過與試驗數據的對比，發(fā)現不同亞格子模型計算準確性存在差異[17].目前，針對繞流場結構更為復雜的三廂車體氣動性能LES計算方法的研究較為罕見.

本文以MIRA三廂車體為對象，對采用LES方法解算非定常特征顯著且具有大分離流動結構的近地鈍體外部繞流場所涉及的迭代步數、時間步長、網格方案等影響因素開展研究；同時分析不同亞格子湍流模型對MIRA車體所受氣動力系數時均值以及車身表面壓力系數時均值的計算準確性；最終提出適用于三廂車型的LES數值仿真策略.

1 大渦模擬

LES的基本思想：通過瞬時N-S方程求解大尺度渦，同時建立模型以量化小尺度渦對大尺度渦的影響，該模型即為亞格子模型(SGS).

(1)

式中：D為流動區(qū)域；x′為實際流動的空間坐標；x為過濾后大尺度空間坐標；G(x,x′)為濾波函數，決定了由N-S方程直接求解的渦尺度范圍.常用的濾波函數有3種：盒式濾波函數、富氏截斷濾波函數和高斯濾波函數.在本文所采用的FLUENT流體仿真平臺中，有限控制體在離散的過程中隱含了過濾運算，濾波函數的表達式為

(2)

(3)

應用式(3)的濾波函數，對質量守恒方程和瞬態(tài)N-S方程進行過濾，得到大渦模擬的控制方程

(4)

(5)

式中：τij為亞格子應力，該項表示了小尺度渦與大尺度渦之間的動量輸運，量化了小尺度脈動對整體流場的影響.τij定義為

(6)

為了使方程組封閉，必須構造相應的亞格子模型，建立關于亞格子應力的數學表達式.對于被過濾掉的小尺度渦團，LES思想認為這部分為各向同性的耗散項.

目前比較普遍的亞格子模型皆沿用RANS方法中Boussinesq提出的渦黏假設[18]，亞格子應力可以表示為

(7)

(8)

渦黏假設將對亞格子應力的求解轉變?yōu)閷喐褡映叨葴u黏系數的求解，并發(fā)展出多種模型.本文重點研究目前常用的3種亞格子模型：DSL(dynamic smagorinsky-lilly)模型[19]，WALE (wall-adapting local eddy-viscosity)模型[20]和DKE(dynamic kinetic energy subgrid-scale)模型[21].

Smagorinsky給出了最基本的亞格子模型，經Lilly改進后形成Smagorinsky-Lilly 模型.該模型對亞格子尺度渦黏系數μt的計算方法為

(9)

Ls=min(kd,CsΔ)

(10)

式中：k為卡門常數；d為節(jié)點到最近壁面的距離；Cs為Smagorinsky常數；Δ為過濾尺度.在研究中發(fā)現,當Cs取為常數時該模型具有一定的局限性，于是發(fā)展出根據求解尺度動態(tài)獲取Cs值的DSL模型.

在WALE亞格子尺度模型中，亞格子尺度渦黏系數的計算方法仍然沿用了混合長度的概念，計算式為

(11)

Ls=min(kd,CwV1/3)

(12)

(13)

DKE模型計入了亞格子尺度湍動能的輸運，其亞格子尺度湍動能定義為

(14)

亞格子尺度渦黏系數則通過ksgs進行計算：

(15)

式中：Ck為模型系數.

2 流場仿真相關信息

仿真對象為1∶3縮尺比MIRA三廂車型：長(L)1 388.3 mm、寬(W)541.7 mm、高(H)473.7 mm.具體構造如圖1所示.

圖1 仿真模型構造

仿真空間區(qū)域，長13.88 m、寬5.42 m、高2.37 m，阻塞比1.61%；x正向為從左到右的空氣流動方向,z正向垂直向上，y正向以右手螺旋定則確定，如圖2所示.

圖2 仿真計算區(qū)域

仿真采用混合網格結構，在車身周圍區(qū)域建立計算子域，采用適應性較好的四面體網格，在外圍區(qū)域采用六面體網格，具體見圖3.

將x、y、z3個方向速度記為u、v、w，計算域入口邊界設為速度入口，速度均勻分布：u=30 m·s-1、v=w=0 m·s-1.出口邊界設為壓力出口，表壓為0 Pa.車體及地面皆采用無滑移壁面邊界條件，計算域左、右兩側及頂部采用對稱邊界條件.按車長計算的雷諾數為Re≈2.81×106.

圖3 體網格布局

LES屬于非穩(wěn)態(tài)計算方法，在時間積分方案上保證二階精度；同時，采用二階精度的空間離散格式.根據文獻[22]的研究成果，以低雷諾數模型計算所得定常解作為LES計算的初始場.

MIRA車體外部繞流為充分發(fā)展的湍流流動，脈動特征顯著，具有很強的時間相關性，任一點的流動參數皆為時間的函數.因此引入物理量的時均值，以便將計算值與試驗值進行對比

(16)

為了對變量隨時間脈動的強弱程度進行評估，引入了物理量的標準差

(17)

式中：σ(φ′)表示變量的標準差.標準差是統(tǒng)計意義下的物理量，數值越大表示此處脈動越劇烈.

3 仿真參數的確定

依托現有計算資源，基于DSL亞格子模型，研究迭代步數、時間步長、初始場、網格數量對LES計算結果的影響規(guī)律.由于LES在近壁區(qū)域的網格策略與低雷諾數模型完全一致，故依據文獻[22]的研究成果，采用車身面網格尺寸1.5 mm，體網格6 745萬單元的網格方案.

3.1 時間步長內最大迭代步數

在MIRA車身上選擇一個監(jiān)測點A，在車體尾跡區(qū)選擇一個監(jiān)測點B，位置如圖4所示.

對于車體上的A點，監(jiān)測其壓力系數Cp隨迭代步數的變化；對于尾跡區(qū)內的B點，監(jiān)測其壓力系數及x、y、z3個方向的速度隨迭代步數的變化.A點、B點壓力系數的變化見圖5.由圖可知：對于2個監(jiān)測點，在任意一個時間步長內，第20次迭代后監(jiān)測點的壓力系數值均達到穩(wěn)定.

圖4 監(jiān)測點位置

a 點A

b 點B

圖6給出了B點處3個方向的速度值隨迭代步數的變化曲線，由圖可知：在每一個時間步長內，3個方向的速度值達到穩(wěn)定所需的迭代步數不同，但在20次迭代后皆達到穩(wěn)定狀態(tài).為確保所有監(jiān)測點都能夠在時間步長內達到穩(wěn)定，LES方法中每一時間步長內最大迭代步數建議為25步.

3.2 時間步長

在非定常流動計算中，當時間步長達到無窮小時，理論上數值仿真可以重現真實的流動狀態(tài)，然而卻不現實.由于時間步長的大小直接影響仿真精度及計算成本，故選取時間步長為5.0×10-4s、2.0×10-4s和1.0×10-4s 3種方案，分別對MIRA三廂車體繞流場進行非定常數值仿真.

在車身表面設置25個監(jiān)測點，如圖7所示，并獲取監(jiān)測點處的壓力系數隨計算時間的變化數據.

圖6 速度變化曲線

圖8給出了數值仿真的對比結果.由圖8可知：3種方案計算所得監(jiān)測點處壓力系數時均值相差較小，因此不同方案對時均結果的影響較小；對比脈動壓力系數的標準差時，發(fā)現在16號至22號監(jiān)測點處3種方案計算結果的標準差最大，并存在明顯差異；1.0×10-4s方案的標準差最大，2.0×10-4s方案與5.0×10-4s方案標準差的數值彼此接近，但2.0×10-4s方案的計算成本是5.0×10-4s方案的2.5倍；3種方案對其余18個監(jiān)測點壓力系數的標準差計算結果相差很小.綜合考慮3種方案的計算精度和計算成本，LES方法中推薦采用5.0×10-4s的時間步長.

圖7 監(jiān)測點位置

圖8 壓力系數平均值及標準差

3.3 網格數量及初始場

網格數量不僅直接決定了初始場的求解質量，而且對LES方法的計算精度及成本影響極大.本節(jié)對車身面網格尺寸1.5 mm、體網格總數6 745萬單元,以及面網格尺寸5 mm、體網格總數1 820萬單元兩個方案進行LES計算的對比分析.文中所涉及的測點位置見文獻[22].由文獻[22]可知，在MIRA三廂車體背部，5 mm網格方案的定常解不對稱性顯著，而1.5 mm網格方案則具有明顯的對稱分布特征.因此,本節(jié)研究同時明確了初始場的不對稱性對LES數值仿真結果的影響.

監(jiān)測氣動阻力系數CD及氣動升力系數CL值，兩種網格方案計算所得時均氣動力系數、標準差及相對試驗值的誤差見表1.由表可知，兩種網格方案的時均CD值相近，僅相差0.001 3，但5 mm網格方案的CD標準差略大；兩種網格方案的時均CL值相差較大，5 mm網格方案對應的CL相對誤差達到134.7%，且該方案的CL標準差也較大.由此表明：5 mm方案的準確性低于1.5 mm方案，且迭代過程的振幅較大，收斂性差于1.5 mm網格方案.

表1 CD 及CL 的對比

在后風窗、行李箱及尾部端面處各選擇一組監(jiān)測點進行詳細比較.圖9給出了壓力系數時均值分布結果.由圖可知，在后風窗區(qū)域，由于流動的非定常性顯著致使兩種網格方案計算結果的對稱性均不理想；在尾部端面及底部上翹面區(qū)域，1.5 mm網格方案的對稱性明顯優(yōu)于5.0 mm方案；1.5 mm方案計算結果的準確性明顯高于5.0 mm方案.

圖10給出了兩種網格方案在MIRA三廂車體表面254個測壓點處壓力系數時均值相對于風洞試驗值的誤差統(tǒng)計結果.由圖可知，1.5 mm方案的計算準確性明顯高于5.0 mm方案；在<10%的誤差范圍內，1.5 mm方案對應83個測點， 5.0 mm方案僅對應58個測點；在>70%的誤差范圍內，1.5 mm方案對應21個測點， 5.0 mm方案則對應49個測點.綜合表1、圖9及圖10的結果可知，網格數量和初始場對LES數值計算結果的準確性影響顯著，且初始場的不對稱性幾乎無法通過LES計算過程予以消除.

a 后風窗

b 行李箱

c 尾部端面

通過對兩種網格方案在氣動力預測準確性、車身表面壓力預測準確性及流場結構對稱性的綜合考察，面網格尺寸1.5 mm、體網格總數6 745萬單元的網格方案適用于LES仿真計算.

4 湍流模型的對比

基于針對LES方法的時間步長、步長內迭代步數、初始場和網格方案的研究成果，本節(jié)采用DSL模型、WALE模型和DKE模型對MIRA三廂車體氣動性能進行仿真計算，通過將氣動性能參數計算值與風洞試驗值進行對比分析，研究3種亞格子模型的計算準確性.文中涉及的測點位置見文獻[22].

圖10 壓力系數誤差統(tǒng)計

圖11為氣動阻力系數時均值的對比結果，倒T線為相對于風洞試驗值的誤差線.由圖可知： DKE模型計算的CD值為0.338，相對誤差13.83%，計算準確性最高；WALE模型計算的CD值為0.345，相對誤差16.32%，計算準確性最低.圖12為氣動升力系數時均值的對比結果.由圖可知， 3種亞格子模型計算所得CL值均小于試驗值； DSL模型的CL計算值為-0.119，誤差最小，WALE模型的CL計算值為-0.154，誤差最大.

圖11 3種模型氣動阻力系數對比

圖12 3種模型氣動升力系數對比

圖13顯示了車身縱向對稱面的壓力系數分布.3種模型針對車體后風窗及行李箱上部表面(x坐標處于[5.2 m，5.55 m]范圍)的壓力系數計算值之間差異顯著，與試驗值相比DKE模型的計算準確度最高；3種模型計算所得車身下表面壓力系數之間的差異較小，預測能力相當.

a 上表面

b 下表面

圖14顯示了在后風窗表面所選取的2組共10個測點的壓力系數時均值分布.由圖可知，3種亞格子模型在該區(qū)域的預測結果趨勢較為一致；與試驗值相比，DKE模型計算所得10個測點處的Cp值誤差皆較小， DSL模型次之， WALE模型的時均值計算誤差最大.

圖15給出了在行李箱表面所選取的2組共10個測點的壓力系數時均值分布.由圖可知， 3種亞格子模型計算所得時均結果與試驗值之間的趨勢一致性較差；與具體的試驗值相比，DKE模型的計算準確性最高，DSL模型次之，WALE模型在該區(qū)域的計算準確性最低.

車體尾部端面所選取的10個測點處壓力系數的分布如圖16所示.由圖可知，3種亞格子模型中，DKE模型的計算結果與試驗值之間的趨勢一致性最高，DSL與WALE模型的趨勢一致性相當；與具體的試驗值相比，3種模型的預測能力相當.

圖14 后風窗壓力系數分布

圖15 行李箱壓力系數分布

圖16 尾部端面壓力系數分布

車體底部上翹面所選取的10個測點處壓力系數的分布如圖17所示.由圖可知，3種亞格子模型中， DSL模型的計算結果與試驗值之間的趨勢一致性最高，DKE與WALE模型的趨勢一致性相當；與具體的試驗值相比，DSL模型的計算準確性最高，DKE與WALE模型的計算能力相當.

圖17 上翹面壓力系數分布

圖18給出了3種亞格子模型在MIRA三廂車體表面254個測壓點處壓力系數時均值相對于風洞試驗值的誤差統(tǒng)計結果.由圖可知，在<10%的誤差范圍內，DKE模型的監(jiān)測點數量最多；在>70%的誤差范圍內，DKE模型有20個監(jiān)測點，DSL模型和WALE模型分別有21個和26個監(jiān)測點，這些測點主要分布在行李箱邊緣和后風窗側面區(qū)域.

圖18 壓力系數誤差統(tǒng)計

5 結論

(1) 應用大渦模擬法計算MIRA三廂車體的氣動性能過程中，時間步長推薦值為5×10-4s，時間步長內最大迭代步數推薦值為25步，1∶3縮尺比模型的面網格推薦尺寸為1.5 mm，推薦采用低雷諾數湍流模型獲得定常初始場，定常初始場應具備明顯的左右對稱特征.

(2) 對三維性較弱的車體頭部及頂部區(qū)域流動，3種亞格子模型的時均結果預測能力相當，均與試驗值符合較好；對存在地面效應和局部分離的底部區(qū)域流動，DSL模型的預測能力最強.

(3) 對三維性較強、存在大分離結構且非定常特征顯著的車體后風窗、行李箱及尾部區(qū)域流動，DKE模型時均結果的計算準確性最高，WALE模型的計算準確性最低.

(4) DKE模型對CD值的計算結果最接近試驗值，DSL模型對CL值的計算結果誤差最?。灰攒嚿肀砻?54個測點處壓力系數時均值為準，DKE模型的計算準確性略優(yōu)于DSL模型，WALE模型的計算準確性最低.

參考文獻:

[1] SMAGORINSKY J. General circulation experiments with the primitive equations [J]. Monthly Weather Review, 1963, 91(3):99.

[2] 馬金花, 金生, 賀德馨, 等. 圓柱體繞流的數值模擬[J]. 山東建筑工程學院學報, 2001,2(6):45.

MA Jinhua, JIN Sheng， HE Dexin，etal. Numerical simulation of fluid flow around a single circular-cylinder [J]. Journal of Shandong Institute of Architecture and Engineering, 2001, 2(6):45.

[3] 苑明順. 高雷諾數圓柱繞流的二維大渦模擬[J]. 水動力學研究與進展(A輯), 1992(7):614.

YUAN Mingshun. Two-dimensional large eddy simulation of flow past a circular cylinder at high Reynold number [J]. Journal of Hydrodynamics (Series A), 1992(7):614.

[4] MURAKAMI S, MOCHIDA A. On turbulent vortex shedding flow past 2D square cylinder predicted by CFD[J]. Journal of Wind Engineering & Industrial Aerodynamics, 1995,54(94):191.

[5] YU D H, KAREEM A. Numerical simulation of flow around rectangular prism[J]. Journal of Wind Engineering & Industrial Aerodynamics, 1997, 67 (2):19.

[6] 童兵, 祝兵, 周本寬. 方柱繞流的數值模擬[J]. 力學季刊, 2002, 23(1):77.

TONG Bing, ZHU Bing, ZHOU Benkuan. Numerical simulation of flow around square cylinder[J]. Chinese Quarterly of Mechanice, 2002, 23(1):77.

[7] 蘇銘德. 直方管內充分發(fā)展湍流的大渦模擬第一部分[J]. 空氣動力學學報, 1995, 13(1):1.

SU Mingde. Larger eddy simulation of fully-developed turbulent flow in a straight duct: part I [J]. ACTA Aerodynamica Sinica, 1995, 13(1): 1.

[8] 楊建明, 吳玉林, 曹樹良. 流體機械中高雷諾數流動的大渦模擬[J]. 工程熱物理學報, 1998,19(2):184.

YANG Jianming, WU Yulin, CAO Shuliang. Large eddy simulation of high Reynolds number flow in fluid machinery [J]. Journal of Engineering Thermophysics, 1998,19(2):184.

[9] 劉沛清, 鄧學鎣. 繞翼型分離流結構的數值研究[J]. 航空學報, 1997,18(4):2.

LIU Peiqing, DENG Xueying. Numerical study of separated flows over an isolated airfoil [J]. Acta Aeronautica et Astronautica Sinica, 1997,18(4): 2.

[10] 劉寧宇, 陸夕云, 莊禮賢. 分層剪切湍流大渦模擬的一種動力學亞格子尺度模型[J].中國科學(A輯), 2000(2):145.

LIU Ningyu, LU Xiyun, ZHUANG LIxian. A subgrid scale model for large eddy simulation of shear turbulence [J]. Science in China (Series A), 2000(2):145.

[11] CHENG Y, LIEN F S, YEE E,etal. A comparison of large eddy simulations with a standardk-εReynolds-averaged Navier-Stokes model for the prediction of a fully developed turbulent flow over a matrix of cubes [J]. Journal of Wind Engineering & Industrial Aerodynamics, 2003, 91(11):1301.

[12] 王兵, 張會強, 王希麟, 等. 后臺階流動再附著過程的大渦模擬研究[J]. 應用力學學報, 2004, 21(3):17.

WANG Bing, ZHANG Huiqiang, WANG Xilin,etal. Investigation on the reattachment process of backward-facing step flow using large eddy simulation [J]. Chinese Journal of Applied Mechanics, 2004, 21(3):17.

[13] CONSTANTINESCU G S, KRAJEWSKI W F, OZDEMIR C E,etal. Simulation of airflow around rain gauges: comparison of LES with RANS models [J]. Advances in Water Resources, 2007,30(1):43.

[14] KRAJNOVIC S, DAVIDSON L. Flow around a simplified car, part 1: large eddy simulation [J]. Journal of Fluids Engineering, 2005, 127(5):907.

[15] KRAJNOVIC S, DAVIDSON L. Flow around a simplified car, part 2: understanding the flow[J]. Journal of Fluids Engineering Transactions of the Asme, 2005, 127(5):919.

[16] SERRE E, MINGUEZ M, PASQUETTI R,etal. On simulating the turbulent flow around the Ahmed body: a French-German collaborative evaluation of LES and DES [J]. Computers & Fluids, 2013, 78(12):10.

[17] ALJURE D E, LEHMKUHL O, RODRíGUEZ I,etal. Flow and turbulent structures around simplified car models [J]. Computers & Fluids, 2014, 96(96):122.

[18] HINZE J O. Turbulence [M]. New York: McGraw-Hill, 1975.

[19] LILLY D K. A proposed modification of the germano subgrid‐scale closure method[J]. Physics of Fluids A Fluid Dynamics, 1992, 4(4):633.

[20] NICOUD F, DUCROS F. Subgrid-scale stress modelling based on the square of the velocity gradient tensor [J]. Flow, Turbulence and Combustion, 1999, 62(3):183.

[21] KIM W W, MENON S. Application of the localized dynamic subgrid-scale model to turbulent wall-bounded flows[R].[S.l.]: AIAA Technical Paper, 1997.

[22] 楊志剛, 周永祥, 朱暉, 等. 低雷諾數k-ε模型鈍體繞流場預測能力研究[J].同濟大學學報(自然科學版), 2017,45(3):413.

YANG Zhigang, ZHOU Yongxiang, ZHU Hui,etal. Comparison of low Reynolds numberk-εmodels in predicting complicated flow field around a bluff body [J]. Journal of Tongji University (Natural Science), 2017,45(3):413.