湍流模型在復雜流場數(shù)值模擬中的應(yīng)用

1、文章編號 1001 246X(2003 04 0335 06收稿日期2002-05-23; 修回日期2002-08-28基金項目總裝備部 十 五! 及民機預研資助項目作者簡介肖志祥(1974- , 男, 四川資陽, 博士生, 主要從事計算流體力學方面的研究.湍流模型在復雜流場數(shù)值模擬中的應(yīng)用肖志祥, 李鳳蔚, 鄂 秦(西北工業(yè)大學114#流體力學研究所, 陜西西安 710072摘 要 采用4種湍流模型:代數(shù)Bald win Lomax (B L 模型、半方程Johnson King (J K 模型的兩個版本(J K90A 和J K92 以及兩方程k g 模型, 分別數(shù)值模擬了導彈超音速流動、N

2、ASA TN D 712標模和民機翼身組合體(兩區(qū)C O 網(wǎng)格 跨音速流動. 采用中心有限體積和多步Runge Kutta 方法數(shù)值積分三維可壓縮雷諾平均Navier Stokes(N S 方程組. k g 湍流模型方程的求解采用類似于N S 方程組的方法進行. 所有湍流模型均能很好地模擬附體及小分離流動; 對于大攻角、分離劇烈的導彈流動, k g 和J K92模型與實驗吻合更好; B L 模型在模擬民機跨音速流動時, 它所捕捉的激波位置較其余3種模型靠后. 利用多塊網(wǎng)格模擬民機翼身組合體流場時, k g 模型的模擬能力強于其余3種模型.關(guān)鍵詞 多塊網(wǎng)格; 翼身組合體; 湍流模型中圖分類號 V

3、211 3文獻標識碼 A0 引言隨著計算機CPU 速度的提高, 硬盤及內(nèi)存容量的持續(xù)增大, 以及數(shù)值計算方法的不斷完善, 計算流體力學(CFD 已經(jīng)成為民用飛機氣動設(shè)計過程中的一個強有力的工具. 然而湍流模型仍然是阻礙人們應(yīng)用N S 方程組進行飛機設(shè)計的瓶頸之一. 近20年來, 人們對粘性流動進行了比較深入的研究, 提出了許多湍流模型, 具有代表性且應(yīng)用較廣的模型是:代數(shù)B L 模型, 半方程J K 模型, 一方程Spalart Allma ras 模型, 以及兩方程k 和k 模型等.評價一種湍流模型的優(yōu)劣, 不能片面地只看它在某一方面的表現(xiàn), 還應(yīng)綜合考慮精度、效率和魯棒性的影響. 幾乎所有

4、的湍流模型都能較好地模擬附體及小分離流動, 對于大攻角、大分離復雜流動如導彈繞流, 沒有任何一種湍流模型能夠表現(xiàn)出超群的模擬能力, 甚至于兩方程模型也是如此1.翼身組合體的求解是進行全機粘性繞流數(shù)值模擬的基礎(chǔ), 解決了翼身組合體的粘性網(wǎng)格生成、湍流模型長度尺度處理以及在此基礎(chǔ)上的其它N S 方程求解技術(shù), 便為數(shù)值分析全機粘性繞流打下了良好的基礎(chǔ). 準確而有效地數(shù)值模擬單獨機翼和機身的粘性流場, 則是求解翼身組合體流場的基礎(chǔ).1 計算方法采用求解橢圓型偏微分方程和代數(shù)方法相結(jié)合, 生成計算網(wǎng)格:物面附近的網(wǎng)格采用代數(shù)方法沿物面法向外推; 遠離物面的網(wǎng)格采用橢圓型方法生成; 將兩種網(wǎng)格合并成單塊

5、整體網(wǎng)格. 該方法生成的網(wǎng)格正交性良好, 網(wǎng)格尺度易于控制.采用中心差分格式加人工粘性的有限體積法離散N S 方程組, 用四步Runge Kutta 方法進行顯式時間推進. 在計算過程中采用當?shù)貢r間步長加速收斂; 對殘值進行光順處理, 能夠有效抑制奇偶波動, 增大依賴域, 增加CFL 數(shù), 加快收斂速度; 對人工粘性的差分格式和邊界處理進行了改進, 考慮到粘性網(wǎng)格大長寬比的特點, 減少了人工粘性對物理粘性的干擾, 使該方法更魯棒, 收斂性也得到了很大的改善.為更好地捕捉激波位置、模擬強激波誘導的分離流動以及增加數(shù)值穩(wěn)定性, 本文對Jameson 人工粘性格式中壓力感應(yīng)因子進行了改進, 以i 方

6、向為例,p , i =|P i -1-2P i +P i +1|(P i -1+2P i +P i +1 +(1- (|P i -P i -1|+|P i +1-P i |,第20卷第4期計 算 物 理C HI NESE JOURNAL OF C OMP UTATI ONAL PHYSICSVol. 20, No. 4其中0 1. 當 =1時, 上式為原始的Ja meson 壓力感應(yīng)因子, 改進后的壓力感應(yīng)因子取 為0 5. 修改后的人工粘性具有類TVD 的性質(zhì), 能有效地避免不真實的數(shù)值振蕩.跨音速流場采用一維Riemann 不變量處理遠場邊界條件; 超音速流場, 在入口處直接用自由來流值,

7、 在出口處, 流動變量由流場內(nèi)變量外插得到. 光滑物面處應(yīng)用無滑移、法向零壓力梯度和絕熱壁邊界條件. 三維流場中, 假定流場關(guān)于x y 平面對稱.2 湍流模型2 1 Baldwin Lomax 模型B L 模型2被廣泛應(yīng)用于C FD 領(lǐng)域, 其優(yōu)勢和缺陷都非常明顯. 即:B L 模型能夠較好地模擬附體流動; 對于較小的局部分離流動, 該模型也有一定的模擬能力; 當采用B L 模型計算跨音速分離流動時, 其激波總是較實際位置靠后.B L 模型是一個基于Prandtl 混合長度理論的雙層代數(shù)模型, 內(nèi)層粘性系數(shù)定義為! t, i(B L =(0 4D (B L # 2, (1其中D (B L =1

8、-exp (-#+/26 , #+=#(w %w 0 5/! w , 為密度, #是網(wǎng)格單元距物面或尾跡割線的法向距離, 為旋度的模, 下標w 代表物面值.外層粘性系數(shù)! t, o(B L =(0 0168 (1 6 F wake &(B L , (2其中F wake =min0 25#max (q max -q min 2/F max , #max F max ,&(B L =1/1+5 5(0 3#/#max 6, F (# =#D (B L ,F max 是函數(shù)F (# 以#為自變量的最大值, #max 為F max 時#的值, q 代表速度. 在尾跡區(qū)域, 令D (B L =1.采用D

9、egani 及Schiff 代數(shù)修正, 選取函數(shù)F (# 距物面第一個極大值作為F max . 同時采用與J K 模型類似的指數(shù)形式合成粘性系數(shù), 使內(nèi)外層粘性系數(shù)光滑過渡, 即! t(B L =! t, o(B L 1-exp(-! t, i(B L ! t, o(B L , (3計算實踐表明, 修正后的B L 模型魯棒性得以提高, 收斂性較原始B L 模型更好. 2 2 Johnson King 模型J K 模型也稱作半方程模型, 需要額外求解一. JK 模型的兩個版本, 其一為Abid 等3于1990年改進的J K90A 模型, 另一個為Johnson4于1992年改進的J K92模型.

10、J K 模型也是基于Prandtl 混合長度理論的兩層模型, 內(nèi)層粘性系數(shù)為! t, i(J K90A =0 4D 2(J K90A #u m ,(4其中D (J K 90A =1-exp-(w #u T /(! w 17 , u T =max(u m , (%m /w 0 5, u m =(-u #v #m 0 5, 設(shè)雷諾剪應(yīng)力為%=! t /, %m 通過求解常微分方程得到.外層粘性系數(shù)定義為! t, o(J K90A = (0 0168 (1 6 F w &(J K 90A ,(5其中&(J K90A =1/1+5 5(#/( 6, F w =#max F max , F (# =#D

11、 (J K90A . 為使內(nèi)外層粘性系數(shù)光滑過渡, 采用指數(shù)函數(shù)合成湍流粘性系數(shù),! t(J K 90A =! t, o(J K 90A 1-exp(-! t, i(J K90A /! t,o(J K90A .(6 模型參數(shù) 將湍流粘性系數(shù)分布同最大雷諾剪應(yīng)力變化率聯(lián)系起來. 最大雷諾剪應(yīng)力隨時間變化的關(guān)系D %m /Dt =%m (%m , eq 0 5-%0 5m a 1/L m -a 1D m , (7 下標m 表示%取最大時的值. 定義L m =min(0 4#m , 0 09( , 其中(為附面層厚度. %m, eq 表示平衡態(tài)下的最大雷諾剪應(yīng)力. 定義湍流耗散項D m =C D %

12、1 5|1- 0 5|/a 1(0 7(-#m , 其中常系數(shù)a 1=0 25, C D =0 5.利用變換h =%-0 5m 線化偏微分方程(7 后求解h , 得到的h 用以更新模型參數(shù)t + t= t m (! m h 2.(8J K92模型是Johnson 和Coakley 在1990年對原始J K 模型修正的基礎(chǔ)上發(fā)展而來. J K92模型的內(nèi)層粘性系數(shù)定義為! t, i(J K92 =(1- ! t, i(B L +! t, i(J K 90A ,(9其中=tanh#(%0 5w +%0 5m /(#m %0 5w , 對u m 進行壓縮性修正, u m =(m /0 5(-u #v

13、 #m 0 5.外層湍流粘性系數(shù)修改為! t, o(J K 92 = (0 0168 (F max (1+ /(1+1 82+ ,(10其中+=0 4F max -u T /|u T |/1 82, 令(=1 2#F /Fmax=1/2. 內(nèi)外層粘性系數(shù)采用雙曲正切函數(shù)合成粘性系數(shù), ! t(J K92 =! t, o(J K92 tanh! t, i(J K92 /! t, o(J K 92 .(11 336計 算 物 理第20卷兩方程湍流模型相對較少地使用代數(shù)或一方程模型的經(jīng)驗公式, 引入了較多的物理概念; 與此同時, 采用Boussinesq 假設(shè)有效避免使用6分量的雷諾應(yīng)力模型, 減少

14、了待求方程數(shù)目. 將湍流模型應(yīng)用到工程CFD 領(lǐng)域的過程中, 人們遇到了許多難題, 對于多塊結(jié)構(gòu)網(wǎng)格系統(tǒng), 具有工程實用價值的兩方程湍流模型應(yīng)滿足以下要求:1 無需使用到物面的垂直距離. 大多數(shù)實用外形包括了兩個甚至更多相交的物面邊界, 比如飛機的翼身結(jié)合部等. 準確定義空間網(wǎng)格單元到物面的距離非常困難, 應(yīng)用需計算到物面距離的湍流模型勢必會增加流場處理難度, 同時造成較大的計算誤差.2 形式簡單的源項. 在物面附近, 源項會變得非常大, 如果處理不當, 將導致數(shù)值計算困難.3 直接的邊界條件. 對于多塊網(wǎng)格尤其重要.Kalitzin 5等提出的k g 模型是從k 模型推導而得, 滿足上述要求

15、, 具有類似于k 模型的優(yōu)點, 同時降低了k 模型在物面附近數(shù)值計算的難度, 易于推廣到具有復雜外形的多塊網(wǎng)格流場.k g 模型方程包含湍動能方程和輸運方程, (u j k j =, j (! l +k ! t , k j +P k -*22R, (12, (u j g j =, j(! l +g ! t j -. *g32R k +22R -(! l +g ! t *R j j ,(13其中k 為湍動能, g 為輸運方程自變量, ! l 為層流粘性系數(shù), P k 為生成項. R =max(0 01! l , *kg 2是一個限制器, 可以有效減弱源項在湍流粘性系數(shù)小于1%層流粘性系數(shù)區(qū)域內(nèi)的

16、影響. 湍流粘性系數(shù)定義為! t =*kg 2, 各常系數(shù)分別為k =g =0 5, . =5/9, *=0 09, =0 075.采用與N S 方程類似的方法求解湍流模型方程. 模型方程邊界條件為:物面處, k =g =0; 入口處, 取來流值; 出口處, k 和g 由場內(nèi)外插得到.3 算例分析采用上述4種湍流模型分別數(shù)值模擬導彈超音速流場, NASA TN D 712標模及民機翼身組合體跨音速流場.3 導彈空間網(wǎng)格如圖1所示, 彈體直徑為D , 前端頭部長3D , 后部平直段長10D . 計算狀態(tài)為M =2 5, . =14%, 基于彈體直徑D 的Re =1 123&106. 該狀態(tài)流動非

17、常復雜, 要準確模擬該流場是一件非常困難的事情.圖1 導彈空間網(wǎng)格(局部 Fig 1 Space grid of the missile(local圖2給出4種模型模擬的彈體表面壓力分布與實驗1對比圖. 所有模型都能準確模擬迎風面附著流動, 其區(qū)別主要集中在背風面, 能較準確模擬背風面流動的湍流模型模擬能力較好. 隨著流動分離, 各湍流模型之間的差異也逐漸明顯起來.圖2 4種湍流模型計算的壓力分布對比Fig 2 Pressure di stribution comparison among the four models圖3給出x /D =6 5截面處各模型數(shù)值模擬的流線對比圖:k g 模型、

18、J K92模型的模擬能力較為出眾, 不僅模擬出強烈的分離主渦, 還模擬出主渦誘導的二次分離渦; k g 模型模擬的二次渦較小, 緊貼于彈體表面; 而J K92模型模擬的二次渦較大, 已發(fā)圖3 不同湍流模型模擬的流線對比Fig 3 Comparison of streamlines among the four models337第4期肖志祥等:湍流模型在復雜流場數(shù)值模擬中的應(yīng)用展成為一個旋渦對; 至于B L 和J K90A 模型, 它們未能模擬出二次渦, 且分離位置也更靠后. k g 模型模擬的主渦要小于其它模型, 渦也更扁平; 其它湍流模型模擬的主渦的強度和形狀則比較相似. 3 2 NASA

19、 TN D 712標模NASA TN D 712標模的機翼具有較小展弦比、大后掠角的外形特點, 以O(shè) H 網(wǎng)格作為計算網(wǎng)格, 見圖4. 機翼上展向位置網(wǎng)格數(shù)一致, 可較好地描述機翼前緣, 有利于捕捉前緣分離渦 .圖4 NASA TN D 712標模表面及空間網(wǎng)格Fig 4 Grids of NASA TN D 712standard model圖5 多湍流模型計算的機翼壓力分布對比Fig 5 Pressure dis tribution comparison on thewing among the four models采用4種湍流模型對NASA TN D 712翼身組合體標模進行數(shù)值模擬,

20、 其狀態(tài)為M =0 9, . =8 4%, Re =7 5&106(文中翼身組合體的Re 數(shù)均基于機翼展長 , 數(shù)值模擬結(jié)果同實驗測量6的壓力分布對比如圖5所示, 所有模型都能較好地描述標模在該狀態(tài)下的流動特點, 且各模型之間差異不大, 區(qū)別在于以下3點:翼根截面靠近后緣處的壓力分布、前緣低壓區(qū)域的大小和翼尖附近的分離形態(tài):B L 模型和J K 模型均能準確模擬翼身結(jié)合部附近的壓力分布, 而k g 模型則遜色一些. 在2z B =95%截面處, k g 模型模擬的壓力分布與實驗吻合較好, B L 模型次之, J K90A 和J K92模型由于在前緣附近模擬出了兩個壓力波動, 與實驗吻合度最小.

21、圖6給出4種湍流模型模擬標模上表面極限流線的比較圖:B L 和k g 模型模擬的流動形態(tài)近似; J K 模型與B L 和k g 模型的最大區(qū)別在于翼尖有一條指向翼尖后緣的斜向分離線, J K90A 和J K92模型則區(qū)別甚微. 通過極限流線圖可以看出:大攻角、大后掠角機翼上表面流動非常復雜, 且分離嚴重, 分離幾乎起始于機翼前緣, 在機翼后緣附近流動重新附體, 幾種湍流模型在模擬如此復雜流動時都有不足之處 .圖6 多湍流模型模擬的機翼上表面流線Fig 6 Streamlines over the upper wing surfaceamong the four models3 3 民機翼身組合

22、體民機機翼具有小后掠角, 大展弦比的特點, 如果仍然按TN D 712那樣生成O H 網(wǎng)格, 那么物面附近網(wǎng)格的正交性和尺度將難以控制. 為使網(wǎng)格線均勻且正交地貼近機翼和機身, 選擇C O 網(wǎng)格作為計算網(wǎng)格, 繞民機翼身組合體的表面及空間網(wǎng)格見圖7,圖中同時給出典型截面網(wǎng)格. 單塊整體網(wǎng)格難以描述整個流場, 為此生成包含上、下兩塊網(wǎng)格的多塊網(wǎng)338計 算 物 理第20卷似, 在塊邊界上以虛擬單元值實現(xiàn) . 圖7 民機表面、空間及典型截面網(wǎng)格圖Fig 7 Grids of surface, space and typical sectionfor the civil ai rplane對民機翼身

23、組合體的如下狀態(tài)進行了計算, M =0 78, . =3%, Re =2 5&106. 該狀態(tài)下民機的攻角較大, 激波較強, 數(shù)值計算和實驗7的壓力分布對比見圖8. 4種模型雖然都能很好地模擬這個狀態(tài), 但它們之間的差別還是比較明顯:B L 模型模擬的激波位置較其余3個模型靠后, 但靠后并不明顯; 雖然采用k g 模型模擬的激波位置在翼尖附近截面稍靠前, 但是k g 模型模擬該狀態(tài)的表現(xiàn)卻是最好的, 它可模擬出壓力分布的細節(jié)變化; J K90A 和J K92模型計算的結(jié)果幾無區(qū)別, 它們模擬的激波位置位于B L 和k g 模型之間, 但是更靠近B L 模型 .圖8 多湍流模型模擬的機翼表面壓力

24、分布對比Fig 8 Pressure distribution comparison on the wingamong the four models圖9給出民機機翼上表面流線對比圖, 采用B L 和k g 模型模擬的結(jié)果表明機翼上表面已經(jīng)分離:B L 模型模擬的分離區(qū)域很小, 位于機翼中間截面; 而k g 模型的分離區(qū)域較大, 已經(jīng)發(fā)展成為一個劇烈的展向渦; 但采用J K90A 和J K92模型模擬的機翼上表面都沒有發(fā)生分離, 且它們之間的差別甚微. 在翼身結(jié)合部的后部, 4種湍流模型均模擬出了一個分離氣泡, 緊貼在機翼機身壁面處.圖9 采用不同模型模擬的流線對比Fig 9 Streamli

25、nes over the civil airplane upper surfacewith the four models4 結(jié)論網(wǎng)格拓撲結(jié)構(gòu)的合理選擇和生成高質(zhì)量的計算網(wǎng)格是精確數(shù)值模擬復雜流場的前提和保證. 文中所采用的湍流模型在模擬附體流動時均可獲得與實驗吻合良好的結(jié)果:修正后的B L 模型在模擬民機強激波流動時, 其激波位置較其余3種模型稍微靠后, 但不明顯; J K90A 和J K92模型在模擬翼身組合體流場時差別很小, 但是在模擬導彈流場時, J K92模型較J K90A 模型模擬能力要好; 雖然在模擬民機時靠近翼尖截面激波位置靠前, 但是k g 模型卻是4種湍流模型中表現(xiàn)最好的,

26、 由于不必計算到物面的垂直距離及直接的邊界條件, 兩方程k g 模型可以非常容易地推廣到復雜的氣動外形和多塊網(wǎng)格系統(tǒng), 這是k g 模型的優(yōu)勢所在.參 考 文 獻1 Sturek W B, et al. The Application of CFD to the Prediction of Missile Body Vortices R. AIAA Paper 97 0637,USA. American Institute of Aeronautics and Astronautics, 1997.2 Bald win B S, Lomax H. Thin Layer Approximati

27、on and Algebraic Model for Separated Turbulent Flows R. AIAA Paper 78 257, USA. American Institute of Aeronau tics and 339第4期肖志祥等:湍流模型在復雜流場數(shù)值模擬中的應(yīng)用340 3 計 算 物 理 第 20 卷 Abid R, et al. Prediction of separated transonic wing flows with non equilibrium algebraic turbulence model J . AIAA Journal, 1990,

28、28: 1426 1431. Models in Aircraft Design R . AIAA Paper 96 0327, USA. American Institute of Aeronautics and Astronautics, 1996. 6 7 陳海昕. 復雜外形繞流 Euler/ N S 方 程數(shù)值分 析 博 士 學位論文 D . 西北工業(yè)大學, 1999. 李杰. 先 進民用飛機 復雜外 形跨音 速流數(shù) 值分析 博 士學位 論 文 D . 西 北 工 業(yè) 大 學 博 士 學 位 論 文, 1999. 4 Johnson D A. Non equilibrium Algeb

29、raic Turbulence Model ing Considerations for Transonic Airfoils and Wings R . AIAA Paper 92 0026, USA. American Institute of Aeronau tics and Astronautics, 1992. 5 Kalitzin G, et al. Aplication of Two equation Turbulence Applications of Turbulence Models in Simulation of Complex Flows XIAO Zhi xiang

30、 , LI Feng wei, E Qin ( I nstitute o Fluids Mechanics, 114# Northwestern Polytechnical University , Xi an 710072, China f Abstract: Four turbulence models: algebraic Baldwin Lomax model with Degani Schiff modification, two versions of Johnson King model ( J K90A and J K92 and two equation k g model,

31、 are described and evaluated for missile supersonic flow, NASA TN D 712 standard model and the civil airplane wing body configuration ( two block C O mesh transonic flows. The 3 D compressible Reynolds averaged Navier Stokes ( RANS equations are integrated numerically by central difference with artificial viscosity scheme, finite volume formulation and explicit multi step Runge Kutta algorithm. The turbulence equations of k g model are explicitly solved in the same way as the RANS. Results show that