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

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

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

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

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

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

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

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

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

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

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

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

0引言

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

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

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

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

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

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

(1)

(2)

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

(3)

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

(4)

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

(5)

從而使Abdol-Hamid的限制只在rPD小于等于1時的流動區域有效，而在其他流動區域自動成為Menter的限制。模型中各常數取值如下：ζ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數值計算方法及邊界設置

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

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

文中所有驗證算例，若無特殊說明，均采用式(5)的限制方法，并且結果均與理論或實驗數據進行對比分析。

3算例與分析

3.1網格收斂性

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

圖1　不同密度網格計算的摩擦系數比較

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

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

圖2　計算摩擦系數隨網格密度變化

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

3.2y+對計算結果的影響

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

圖3　不同y+計算的摩擦系數比較

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

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

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

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

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

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

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

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

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

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

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

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

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

(a) 摩擦系數

(b) 速度型

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

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

(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模型給出的壓力分布在上表面后緣分離區之前均略高于SST模型，在分離區內，比SST模型略低。

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

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

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

圖8　后臺階下壁面摩擦系數對比

4結論

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

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

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

參考文獻

[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-)，男，碩士，高級工程師。主要研究方向：計算流體力學、飛行器設計。

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

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

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

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

(編輯：趙毓梅)