李廣佳，李喜樂，張強，郝海兵，李典

(1.中國航天空氣動力技術(shù)研究院 第十一總體設(shè)計部，北京　100074)

(2.西北工業(yè)大學 航空學院，西安　710072)

(3.中國航空計算技術(shù)研究所 第七研究室，西安　710068)

k-kL兩方程湍流模型的改進及驗證研究

李廣佳1，李喜樂1，張強2，郝海兵3，李典3

(1.中國航天空氣動力技術(shù)研究院 第十一總體設(shè)計部，北京100074)

(2.西北工業(yè)大學 航空學院，西安710072)

(3.中國航空計算技術(shù)研究所 第七研究室，西安710068)

摘要：提高湍流數(shù)值模擬的準確性，從而明確湍流模型對數(shù)值模擬結(jié)果的影響具有重要的意義，應(yīng)用K.S.Abdol-Hamid給出的尺度自適應(yīng)k-kL兩方程模型封閉RANS方程，并修改von Karman長度尺度的限制方法，通過平板、翼型、后臺階等流動的模擬，考察k-kL模型在湍流模擬中的準確性，及其反映主要流動特征的能力和網(wǎng)格收斂性，并對影響流動模擬準確性的因素進行討論。結(jié)果表明：改進長度尺度限制之后的k-kL兩方程模型無論是對附著流動還是分離流動都可以給出比較準確的結(jié)果。

關(guān)鍵詞：k-kL兩方程湍流模型；RANS；von Karman長度尺度

0引言

A.N.Kolmogorov于1942年提出的兩方程湍流模型是研究所有基于統(tǒng)計平均的湍流模型的基礎(chǔ)。兩方程湍流模型的概念反映了湍流研究的基本思想，即在平均流動中模擬湍流效應(yīng)需要兩個相互獨立的湍流尺度變量，這兩個湍流尺度可以通過各自的輸運方程求解給出。此外，兩方程湍流模型也是其他更高階湍流模型的基礎(chǔ)和核心，例如雷諾應(yīng)力模型、顯式代數(shù)應(yīng)力模型或者其他基于非線性應(yīng)力-應(yīng)變關(guān)系的湍流模型。即使只采用渦粘系數(shù)作為唯一變量的一方程湍流模型(例如SA湍流模型[1])，也可以在平衡假設(shè)(Equilibrium Assumptions)的基礎(chǔ)上，從兩方程模型推導(dǎo)而來。在一方程湍流模型中，平衡假設(shè)認為湍流時間尺度與剪切應(yīng)變率成反比[2]，因此不需要第二個尺度方程。

在目前所有的兩方程湍流模型中，第一個湍流尺度方程幾乎都采用湍流脈動動能k方程。只需要對湍流擴散項(Turbulent Diffusion)進行建模，從而可以降低建模的難度。至于第二個湍流尺度方程，目前發(fā)展了很多[3]，最常見的是ε和ω方程。盡管第二個湍流尺度方程(例如ε方程)可以通過嚴格推導(dǎo)得到準確的形式，但鮮有實際應(yīng)用。準確的形式中包含更復(fù)雜和更高階的關(guān)聯(lián)項，對于這些項目前只能給出量級大小估計[4]，難以采用逐項建模實現(xiàn)，因此實際應(yīng)用中的ε和ω方程通過量綱分析或者直觀物理分析得到，并采用與k方程類似的形式。

針對第二個湍流尺度，J.C.Rotta[5]采用逐項建模的方法給出了一種(kL)方程，其中的L指的是湍流的積分長度尺度，k是湍流脈動動能。J.C.Rotta提出的kL方程源項中，包含速度的三階導(dǎo)數(shù)項，這在實際應(yīng)用存在諸多問題，限制了其應(yīng)用[6]。F.R.Menter等[7-9]采用二階速度導(dǎo)數(shù)的合理形式取代三階速度導(dǎo)數(shù)項，表明采用速度二階導(dǎo)數(shù)項后，k-kL方程可以自動滿足壁面對數(shù)律，并且使模型可以自動調(diào)整湍流模擬尺度以滿足復(fù)雜多尺度湍流結(jié)構(gòu)模擬的需求，而無需像傳統(tǒng)RANS方程耗散掉。在F.R.Menter修改的基礎(chǔ)上，K.S.Abdol-Hamid等[10-11]給出了完整形式的k-kL兩方程模型，并對模型所用的常數(shù)進行了標定。

目前，國內(nèi)對于基于湍流積分長度尺度的湍流模型，未見報道，因此本文對尺度自適應(yīng)k-kL兩方程模型開展初步的研究。為了提高k-kL模型的數(shù)值穩(wěn)定性和準確性，修改von Karman長度尺度的限制方法；通過對典型流動的模擬，考察改進后的k-kL模型在定常湍流模擬中的準確性、捕捉主要流動特征的能力及其網(wǎng)格收斂性，并對影響流動模擬的因素進行討論。

1k-kL湍流模型控制方程

守恒形式的k-kL兩方程湍流模型控制方程可以寫為

(1)

(2)

Lvk為von Karman長度尺度，為了避免長度尺度之比(L/Lvk)過大或者過小，對Lvk取值限制如下：

(3)

不同的rPD對應(yīng)不同的von Karman長度尺度限制方法。F.R.Menter根據(jù)附著流動湍流生成和耗散之比約等于1，取rPD=1。但是，當流動有分離時rPD卻不等于1，于是K.S.Abdol-Hamid定義

(4)

采用K.S.Abdol-Hamid的定義可以極大地改善k-kL模型在分離區(qū)的模擬精度，但是對全流場都采用此限制也不妥，因此把rPD修改為

(5)

從而使Abdol-Hamid的限制只在rPD小于等于1時的流動區(qū)域有效，而在其他流動區(qū)域自動成為Menter的限制。模型中各常數(shù)取值如下：ζ1=1.20，ζ2=0.97，ζ3=0.13，σk=1.00，σφ=1.00，κ=0.41，Cμ=0.09，Cl1=10.00，Cl2=1.30，Cd1=4.70。

2數(shù)值計算方法及邊界設(shè)置

流動控制方程基于多塊結(jié)構(gòu)網(wǎng)格采用基于格心的二階有限體積法進行空間離散，其中無粘項采用Roe-FDS格式離散，粘性項采用中心差分近似。但在湍流模型無粘項的離散中僅采用一階迎風格式。RANS方程和湍流模型方程采用主對角占優(yōu)DDADI格式以松耦合的方式進行隱式時間推進。為了更加快速地得到定常解，采用當?shù)貢r間步長和多重網(wǎng)格加速收斂技術(shù)。

湍流模型變量k和kL在物面均取0，自由來流條件只需保證流動能自動轉(zhuǎn)捩為全湍流。

文中所有驗證算例，若無特殊說明，均采用式(5)的限制方法，并且結(jié)果均與理論或?qū)嶒灁?shù)據(jù)進行對比分析。

3算例與分析

3.1網(wǎng)格收斂性

為了考察k-kL模型的網(wǎng)格收斂性，分別在545×385(極密網(wǎng)格，Super-Fine)，273×197(密網(wǎng)格，F(xiàn)ine)，137×97(中等密度網(wǎng)格，Medium)，69×49(稀疏網(wǎng)格，Coarse)，35×25(極稀疏，Tiny-Coarse)從密到疏的五套網(wǎng)格上，計算Ma為0.2、單位平板長度Re為5×106的平板附面層(各網(wǎng)格中的平板長度取為x/c=2.00)的摩擦系數(shù)。不同密度網(wǎng)格上沿平板的摩擦系數(shù)與理論值的比較如圖1所示。

圖1　不同密度網(wǎng)格計算的摩擦系數(shù)比較

從圖1可以看出：密網(wǎng)格和極密網(wǎng)格上的計算結(jié)果差別較小、幾乎重合，在中等密度網(wǎng)格上的計算結(jié)果與密網(wǎng)格和極密網(wǎng)格上的差別已經(jīng)很小，而在稀疏網(wǎng)格和極稀疏網(wǎng)格上的計算結(jié)果與網(wǎng)格和極密網(wǎng)格上的差別較大，并且網(wǎng)格越稀疏，計算結(jié)果越偏離理論解；受附面層粘性底層建模項的影響(k和kL方程右端第三項)，在平板前緣的轉(zhuǎn)捩間斷比其他兩方程模型更為明顯，并且轉(zhuǎn)捩點的位置隨著網(wǎng)格的加密前移，更接近平板前緣。

不同密度網(wǎng)格上計算的x/c=0.97處摩擦系數(shù)的比較如圖2所示。

圖2　計算摩擦系數(shù)隨網(wǎng)格密度變化

從圖2可以看出：隨著網(wǎng)格的加密，摩擦系數(shù)變化趨勢與CFL3D和FUN3D在同系列網(wǎng)格上用SSTk-ω模型所得變化趨勢一致[12]，并且在最密網(wǎng)格上的摩擦系數(shù)基本相同。

3.2y+對計算結(jié)果的影響

為了考察y+對計算結(jié)果的影響，首先在x/c=1.00長的平板上生成65×97的網(wǎng)格(流向網(wǎng)格節(jié)點總數(shù)65，物面法向網(wǎng)格節(jié)點總數(shù)97)，然后通過調(diào)整物面第一層網(wǎng)格到物面的距離和物面法線方向的網(wǎng)格增長率，得出六套網(wǎng)格。在Ma為0.2、單位平板長度Re為6×106時，k-kL模型在不同網(wǎng)格上計算的平均y+從小到大依次約為0.020，0.220，0.505，1.150，2.290，4.600。在不同y+時計算的摩擦系數(shù)與理論值的比較如圖3所示。可以看出：當y+=0.220，0.505，1.150，2.290時計算的摩擦系數(shù)沿平板變化曲線接近；當y+進一步增加或者減小時，平板摩擦系數(shù)都會變小，進一步偏離理論解。

圖3　不同y+計算的摩擦系數(shù)比較

為了進一步說明y+對計算結(jié)果的影響，給出x/c=0.79處不同y+時計算摩擦系數(shù)的比較，如圖4所示。

圖4　x/c=0.79處不同y+計算的摩擦系數(shù)

從圖4可以看出：當y+分別為0.220，0.505，1.150時計算所得摩擦系數(shù)幾乎相等；y+=2.290時的結(jié)果相比前三種計算值略小，但差別小于0.1個阻力單位(1×104為1阻力單位)；y+=0.020和4.600時的結(jié)果與其他結(jié)果差別較大，接近甚至超過0.25個阻力單位。結(jié)果表明：k-kL模型與其他兩方程模型相同，為了保證附面層模擬的精度，y+約等于1最優(yōu)。

3.3不同湍流初始條件的影響

分別計算兩種湍流初始條件下的平板附面層，如表1所示。

表1　Ma=0.2時不同的湍流初始條件

在兩種湍流初始條件下，不同密度網(wǎng)格上計算x/c=0.97處的摩擦系數(shù)，如圖5(a)所示，可以看出：各網(wǎng)格上計算結(jié)果幾乎相同，差別可以忽略不計。在兩種湍流初始條件下，中等密度網(wǎng)格上的摩擦系數(shù)如圖5 (b)所示，摩擦系數(shù)曲線也幾乎處處重合。

(a) x/c=0.97摩擦系數(shù)

(b) 不同湍流條件下的摩擦系數(shù)

綜上所述，只要給定的自由來流湍流度Tu和μt/μ∞能夠保證流動轉(zhuǎn)捩為全湍流，k-kL模型的計算結(jié)果與自由來流湍流條件無關(guān)，即k-kL模型對自由來流的湍流條件并不敏感。

3.4不同長度尺度限制的影響

本節(jié)驗證算例考慮無分離平板流動、NACA4412翼型粘性繞流和后臺階流動。

在無分離的平板流動中，采用式(5)的限制和采用Menter的限制計算的摩擦阻力系數(shù)和速度型幾乎完全相同，如圖6所示。

(a) 摩擦系數(shù)

(b) 速度型

從圖6可以看出：在無分離或遠離物面的區(qū)域，式(5)采用的是Menter的限制；與采用其他兩種限制方法相比，采用Abdol-Hamid的限制會使計算的摩擦系數(shù)略微增加、在離開物面相同位置處計算的速度稍偏大，但差別不大。

為了進一步說明式(5)的修改在小分離流動中的效果，對NACA4412翼型粘性繞流進行模擬。計算的翼型表面壓力分布及上表面六處站位的速度型與實驗值的比較如圖7所示，前三處站位在分離之前，后三處站位處于分離區(qū)。

(a) 翼型表面壓力分布

(b) 在x/c=0.620處的速度型

(c) 在x/c=0.675處的速度型

(d) 在x/c=0.731處的速度型

(e) 在x/c=0.786處的速度型

(f) 在x/c=0.842處的速度型

(g) 在x/c=0.897處的速度型

(h) 在x/c=0.953處的速度型

從圖7(a)可以看出：k-kL模型給出的壓力分布在上表面后緣分離區(qū)之前均略高于SST模型，在分離區(qū)內(nèi)，比SST模型略低。

從圖7(b)～圖7(h)可以看出：采用式(5)限制方法的k-kL模型所給出的速度型，在前三個站位以及后三個站位的附面層外部區(qū)域，與采用Menter的限制給出的結(jié)果接近，而在分離區(qū)內(nèi)靠近物面區(qū)域，與采用Abdol-Hamid的限制給出的結(jié)果類似；采用式(5)限制方法給出的結(jié)果在無分離區(qū)與SST模型相當，但在分離區(qū)的結(jié)果要優(yōu)于SST模型，與實驗值更為吻合。

綜上所述，不管有無分離，采用式(5)的限制方法可以給出最為準確的數(shù)值模擬結(jié)果。

針對后臺階流動算例，生成四塊網(wǎng)格129×129，49×129，193×225，65×225，離開物面第一層網(wǎng)格的平均y+約為0.887。在該網(wǎng)格上，采用不同von Karman長度尺度限制方法計算后臺階下壁面的摩擦系數(shù)，如圖8所示?？梢钥闯觯翰捎肕enter的限制方法，計算結(jié)果不收斂，流動是非定常的，而采用Abdol-Hamid和式(5)進行計算，均可以得到準定常的收斂結(jié)果，并且給出的下壁面摩擦系數(shù)差別不大，與實驗值吻合良好，k-kL模型的結(jié)果要略好于SSTk-w模型。

圖8　后臺階下壁面摩擦系數(shù)對比

4結(jié)論

(1)k-kL模型具有優(yōu)異的網(wǎng)格收斂性；為了準確求解附面層，物面第一層網(wǎng)格的平均y+約為1的量級最優(yōu)。

(2)k-kL模型對湍流自由來流條件不敏感，只需保證流動能轉(zhuǎn)捩為全湍流即可。

(3) 采用式(5)的von Karman長度尺度限制方法，不管流動是否有分離，k-kL模型都可以給出較為準確的結(jié)果。

參考文獻

[1] Spalart P R, Allmaras S R. A one-equation turbulence model for aerodynamic flows[J]. La Recherche Aerospatialen, 1994, 1: 5-21.

[2] Menter F R. Eddy viscosity transport equations and their relation to thek-εmodel[J]. Journal of Fluids Engineering, 1997, 119: 876-884.

[3] Kantha L H. The length scale equation in turbulence models[J]. Nonlinear Processes in Geophysics, 2004, 11(1): 83-97.

[4] Temmerman L, Leschziner M A. Large eddy simulation of separated flow in a streamwise periodic channel constriction[C]. Stockholm: 2nd Symposium on Turbulence and Shear-Flow Phenomena, 2001.

[5] Rotta J C. Statistische theorie nichthomogener turbulenz[J]. Zeitschrift für Physik, 1951, 129: 547-572.

[6] Rodi W. Turbulence modelling for boundary layer calculations[C]. G?ttingen: Springer, IUTAM Symposium One Hundred Years of Boundary Layer Research, 2004.

[7] Menter F R, Egorov Y. The scale-adaptive simulation me-thod for unsteady turbulent flow predictions(Part 1): theory and model description[J]. Flow Turbulence Combust, 2010, 85: 113-138.

[8] Menter F R, Egorov Y. The scale-adaptive simulation me-thod for unsteady turbulent flow predictions(Part 2): application to complex flows[J]. Flow Turbulence Combust, 2010, 85: 139-165.

[10] Abdol-Hamid K S. Assessmants ofk-kLturbulence model based on Menter’s modification to Rotta’s two-equation model[J]. International Journal of Aerospace Engineering, 2015(1): 1-18.

[11] Abdol-Hamid K S. Assessments of a turbulence model based on Menter’s modification to Rotta’s two-equation model[C]. AIAA-2013-0341, 2013.

[12] Chris Rumsey. Turbulence modeling resource[EB/OL].(2015-11-05)[2016-04-20].http:∥turbmodels.larc.nasa.gov.

Improvement and Validation ofk-kLTwo-equation Turbulence Model

Li Guangjia1, Li Xile1, Zhang Qiang2, Hao Haibing3, Li Dian3

(1.The Eleventh Design Department, China Academy of Aerospace Aerodynamics， Beijing 100074， China)

(2.School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China)

(3.No. 7 Department, Aeronautical Computing Technique Research Institute, Xi’an 710068, China)

Abstract:It is of important significance to improve the accuracy of turbulent flow simulations, and identify the effects of turbulence model on simulation results. The scale adaptive k-kL two-equation turbulence model constructed by K.S.Abdol-Hamid is applied to close RANS equations and assessed. A modification to the upper limit in von Karman length scale is proposed. Some factors effecting simulation accuracy are also discussed. A number of test cases including flat plate case, airfoil case and backward facing step case are run to demonstrate the accuracy of this model and its capability to resolve important flow features and grid convergence. Test results show that this turbulence model is superior or at least comparable to other two-equation turbulence models in simulating either attached or separated turbulent flows.

Key words:k-kL two-equation turbulence model; RANS; von Karman length scale

收稿日期：2016-04-20；修回日期：2016-05-10

通信作者：李喜樂,lxl1027@163.com

文章編號:1674-8190(2016)02-209-07

中圖分類號：V211

文獻標識碼：A

DOI：10.16615/j.cnki.1674-8190.2016.02.011

作者簡介：

李廣佳(1979-)，男，碩士，高級工程師。主要研究方向：計算流體力學、飛行器設(shè)計。

李喜樂(1985-)，男，博士，工程師。主要研究方向：計算流體力學、設(shè)計空氣動力學。

張強(1979-)，男，博士，副教授。主要研究方向：理論與計算流體力學。

郝海兵(1981-)，男，博士，高級工程師。主要研究方向：理論與計算流體力學、氣動優(yōu)化設(shè)計。

李典(1986-)，男，博士，工程師。主要研究方向：計算流體力學。

(編輯：趙毓梅)