帶有剪切應(yīng)力輸運(yùn)性質(zhì)的k-ε兩方程湍流模型構(gòu)造與應(yīng)用

王 偉,張 揚(yáng),陳利麗

(1.航空工業(yè) 第一飛機(jī)設(shè)計(jì)研究院,西安 710089;2.西安交通大學(xué) 機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,西安 710049)

0 引 言

在空氣動(dòng)力學(xué)所涉及的流動(dòng)中，逆壓梯度引起的分離流動(dòng)是極為常見(jiàn)的流動(dòng)現(xiàn)象之一，尤其是在航空器的氣動(dòng)設(shè)計(jì)領(lǐng)域，這一流動(dòng)現(xiàn)象多帶來(lái)氣動(dòng)性能損失、結(jié)構(gòu)疲勞和操作性穩(wěn)定性變差,會(huì)嚴(yán)重威脅航空器的飛行性能和安全。雖然直接數(shù)值模擬方法(DNS)和大渦數(shù)值模擬方法(LES)能夠描述湍流場(chǎng)的結(jié)構(gòu)和時(shí)間歷程，但是對(duì)于工業(yè)界中廣泛的氣動(dòng)力評(píng)估，仍然需要依賴于高效魯棒的雷諾平均方法(RANS)。然而，普遍存在的非平衡效應(yīng)是RANS方法所需要面對(duì)的最大求解障礙之一。RANS方法中對(duì)于湍流脈動(dòng)的求解仍然基于雷諾分解假設(shè)，即關(guān)于雷諾應(yīng)力的合理建模。早期的純代數(shù)封閉模型，如Baldwin-Lomax(BL)模型[1]，由于缺乏對(duì)于湍流變量的歷史效應(yīng)反饋，在求解帶有分離現(xiàn)象的流動(dòng)時(shí)精度較差。但是，BL模型提供了一種雷諾應(yīng)力的簡(jiǎn)化封閉思路——即完成湍流速度尺度、長(zhǎng)度尺度和時(shí)間尺度中的至少兩個(gè)量,即可完成在線性渦黏假設(shè)下的雷諾應(yīng)力封閉。為了求解這三個(gè)變量，學(xué)者們發(fā)展了“N”方程模型，其中N代表求解基本變量所需的偏微分方程數(shù)目。例如，由Baldwin等發(fā)展的一方程模型(BB模型)[2]，以沒(méi)有經(jīng)過(guò)壁面衰減的渦黏系數(shù)為輸運(yùn)變量構(gòu)造了一個(gè)輸運(yùn)方程，雖然表面上看似乎在效率與精度間達(dá)到了最好的平衡，并且引入了時(shí)間效應(yīng)，然而經(jīng)典的BB模型是由k-ε兩方程模型推導(dǎo)而來(lái)，數(shù)項(xiàng)擴(kuò)散相關(guān)項(xiàng)被省略，這必然引起湍流求解精度的下降。Menter等[3]的研究表明BB模型在剪切層邊緣的模擬精度較差，Rahman等[4]認(rèn)為BB模型的缺陷是由于破壞項(xiàng)的病態(tài)與對(duì)來(lái)流的極為敏感造成。Spalart等[5]基于經(jīng)驗(yàn)和類(lèi)渦黏性系數(shù)提出的Spalart-Allmaras模型(SA模型)是目前最受歡迎的一方程模型，其構(gòu)造理念采用“堆積木”的形式完成，通過(guò)簡(jiǎn)單流動(dòng)一步步進(jìn)行標(biāo)定，最終完成對(duì)于復(fù)雜流動(dòng)的求解，壁面效應(yīng)被融入了等價(jià)渦黏系數(shù)的比例系數(shù)中。SA模型在飛行器流動(dòng)中的應(yīng)用極為廣泛。但是，Menter[6]認(rèn)為SA模型在預(yù)測(cè)逆壓梯度流動(dòng)中精度欠佳，尤其是在回流問(wèn)題中容易將回流速度過(guò)高估計(jì)。Coakley提出了渦黏系數(shù)的混合求解模型[7]，Menter[3,6]將這個(gè)理念發(fā)展為基于k-ω兩方程的剪切應(yīng)力輸運(yùn)湍流模型，以Bradshaw假設(shè)[8]為限制，改善了在逆壓梯度和分離流動(dòng)出現(xiàn)時(shí)，傳統(tǒng)湍流模型將產(chǎn)生項(xiàng)做過(guò)大估計(jì)的缺陷，并且在該假設(shè)起效時(shí)，整個(gè)湍流模型實(shí)質(zhì)上變成了剪切應(yīng)力的輸運(yùn)方程。在此之后，有學(xué)者[9-11]將這種理念發(fā)展到一方程湍流模型，并在復(fù)雜流動(dòng)的求解中得到了很好的應(yīng)用。然而，在k-ε兩方程模型中，還未有剪切應(yīng)力概念的引入，本文將原始的Bradshaw假設(shè)采用直接數(shù)值模擬(DNS)數(shù)據(jù)重新標(biāo)定，基于結(jié)構(gòu)系綜動(dòng)力學(xué)[12]理論給出具有物理意義的標(biāo)定方程，以此消除原始SST方程中的混合函數(shù)項(xiàng)。

另一方面，如圖1所示，基于DNS數(shù)據(jù)的反求以及標(biāo)準(zhǔn)k-ε湍流模型的對(duì)比發(fā)現(xiàn)可知，近壁面附近的渦黏系數(shù)與DNS對(duì)應(yīng)值差距較大。其原因在于通過(guò)攝動(dòng)法分析后，采用k-ε兩方程模型中的渦黏系數(shù)求解方法(Cμk2/ε)與渦黏性系數(shù)(μt)漸進(jìn)解不一致，需要乘以一定的系數(shù)來(lái)滿足固壁面的影響來(lái)得到正確的邊界解。因此衍生了一種新的湍流模型子集即低雷諾數(shù)湍流模型(Low Reynolds Number,LRN)。該種模型的設(shè)計(jì)初衷并非針對(duì)特征雷諾數(shù)很低的流動(dòng)，而是為了更正如圖1所示的現(xiàn)象，修正當(dāng)?shù)乩字Z數(shù)較低時(shí)無(wú)法產(chǎn)生正確的邊界層漸進(jìn)行為而提出。在近壁面區(qū)域通過(guò)逐步衰減系數(shù)來(lái)完成均衡態(tài)假設(shè)的近壁效應(yīng)修正，該修正函數(shù)被稱(chēng)為壁面衰減函數(shù)(WDF)。因此低雷諾數(shù)湍流模型中的“雷諾數(shù)”指代的并非是模擬流動(dòng)的特征雷諾數(shù)。LRN湍流模型中，Jones等[13]提出了第一種WDF函數(shù)并得以應(yīng)用，在近十年來(lái)，幾十種湍流模型相繼被提出，隨之帶來(lái)的便是數(shù)種WDF函數(shù)以及額外的源項(xiàng)引入，但最終目的都是為了使湍流模型本身產(chǎn)生合理的近壁漸進(jìn)行為。Patel[14]評(píng)估了9種LRN模型后，通過(guò)湍流邊界層算例展示這些模型的精度，雖然所有模型都能夠完美滿足湍流對(duì)數(shù)率，然而，在湍動(dòng)能的預(yù)測(cè)上表現(xiàn)的差距非常大，其根本原因來(lái)自于WDF函數(shù)的構(gòu)造方法。常用來(lái)構(gòu)造的基本元素是兩種湍流雷諾數(shù)(Reτ和Rek)和無(wú)量綱的壁面距離y+，然而，y+引入的壁面摩擦速度在流動(dòng)的分離點(diǎn)和再附點(diǎn)處為0，容易引起WDF函數(shù)的病態(tài)。另外，如圖1所示，WDF的合理作用范圍應(yīng)該不高于對(duì)數(shù)層，但是，多數(shù)LRN模型在構(gòu)造WDF函數(shù)時(shí)，過(guò)多關(guān)注WDF的衰減而忽略了WDF函數(shù)的峰值變化。Rodi等人[15]的研究清晰地表現(xiàn)出WDF的峰值是隨著雷諾數(shù)的變化而變化，因此，本文另一個(gè)工作即為提出一個(gè)峰值可變的WDF函數(shù)形式。

新模型通過(guò)充分發(fā)展的湍流平板、NACA4412翼型、二維鼓包以及DLR-F6翼身組合體的計(jì)算來(lái)驗(yàn)證計(jì)算精度。

圖1 基于DNS和標(biāo)準(zhǔn)k-ε模型求得的渦黏性系數(shù)μtFig.1 Eddy-viscosity coefficient got by DNS and standard k-ε turbulence model

1 湍流模型演化思路

1.1 k-ε兩方程湍流模型

忽略掉浮力項(xiàng)的影響，標(biāo)準(zhǔn)不可壓縮k-ε湍流模型的形式為：

(1)

(2)

其中:ρ為密度，Uj為速度矢量，k為湍動(dòng)能，μ為層流黏性系數(shù)，μt為湍流黏性系數(shù)，Pk為湍動(dòng)能產(chǎn)生項(xiàng)，ε為湍動(dòng)能耗散率，t為時(shí)間變量，xj為空間張量，Cε1、Cε2為常數(shù)。

Pk=μtΩ2

(3)

1.2 WDF函數(shù)的構(gòu)造

通常，對(duì)于二維邊界層流動(dòng)，將近壁面附近速度脈動(dòng)量(即u′、v′和w′ )采用小擾動(dòng)方式泰勒展開(kāi)，可得：

(4)

由于受無(wú)滑移邊界條件的約束，式(4)中的速度脈動(dòng)右端均無(wú)常數(shù)項(xiàng)，同時(shí)，聯(lián)合連續(xù)方程可得b1=0。通過(guò)式(4)，可以得到近壁附近湍動(dòng)能和湍動(dòng)能耗散率的階數(shù)為：

(5)

式(5)表明，在近壁面附近保持Cμk2/ε的湍流渦黏性系數(shù)顯然是不合理的，需要一定的降階函數(shù)使μt保持正確的邊界漸進(jìn)性。此類(lèi)函數(shù)即為WDF函數(shù)，WDF的精確表達(dá)式為：

(6)

Zhang等人[16]結(jié)合DNS數(shù)據(jù)，將三種常見(jiàn)湍流模型(Abidk-ε[17]、Yang-Shihk-ε[18]和Launder-Sharmak-ε[19])中的WDF函數(shù)通過(guò)式(4)求解出來(lái)。如圖2，隨著雷諾數(shù)的變化，三種模型的WDF函數(shù)并未隨著一同改變，而是被機(jī)械地限制在了值為1附近。通過(guò)表1可以看出，三種模型的WDF函數(shù)在自變量趨于無(wú)窮時(shí)，峰值均為1，然而從DNS數(shù)據(jù)反求出的WDF峰值是隨著雷諾數(shù)不同而變化的。

參考表1給出三種對(duì)比模型的WDF形式，從數(shù)學(xué)性質(zhì)上可以看出，三個(gè)WDF的極值隨著自變量的增加均為1。而圖2的DNS數(shù)據(jù)反推表明，這樣的極值是不合理的。

根據(jù)觀察，Zhang等人[16]以因子Rek/Reτ來(lái)表征雷諾數(shù)的變化，通過(guò)構(gòu)造兩個(gè)不同的函數(shù)fμ1和fμ2，基于此形成新的WDF函數(shù)：

(5)

式(5)的構(gòu)造理念在于——fμ1是將現(xiàn)有的WDF進(jìn)行聯(lián)合，以此達(dá)到近壁面附近更精確的衰減作用，通過(guò)fμ2控制峰值的增減。圖3為不同雷諾數(shù)槽道中Rek/Reτ的表達(dá)。如圖4所示，新構(gòu)造的WDF函數(shù)的峰值與雷諾數(shù)變化具有一種自適應(yīng)關(guān)聯(lián),其中fnc代表非均衡湍流識(shí)別函數(shù)。

(a)Change factor with the change of y+

(b)New WDF’s structure function fμ2圖3 新構(gòu)造的WDF帶入DNS數(shù)據(jù)中的變化Fig.3 Changes of newly constructed WDF into DNS data

圖4 新WDF函數(shù)在不同雷諾數(shù)下的變化Fig.4 Variation of new WDF function at different Reynolds numbers

1.3 改進(jìn)的湍渦黏性系數(shù)計(jì)算方式

一般情況下，在k-ε兩方程模型中，湍渦黏性系數(shù)的確定方式為：

μt=Cμfμk2/ε

(6)

式(6)是從均勻各向同性平衡湍流中推導(dǎo)而來(lái)，因此對(duì)于帶有分離和逆壓梯度的流場(chǎng)適應(yīng)性不強(qiáng)。Menter[6]的研究表明，在逆壓梯度的流動(dòng)中，湍動(dòng)能產(chǎn)生項(xiàng)可以數(shù)倍大于其破壞項(xiàng)，導(dǎo)致雷諾應(yīng)力的過(guò)大估計(jì)，因此引入了Bradshaw假設(shè)進(jìn)行平衡。然而，Brashaw將湍動(dòng)能與雷諾應(yīng)力間的比值設(shè)定為常值(0.31)，而Nagano[11]和文曉慶等[20]的研究均表明該常數(shù)是一個(gè)可變值。張揚(yáng)[21]等根據(jù)已有的DNS數(shù)據(jù)，反推出Bradshaw常數(shù)，并以佘振蘇[12]等人提出的結(jié)構(gòu)系綜動(dòng)力學(xué)理論對(duì)新的常數(shù)進(jìn)行標(biāo)定。本文采用張揚(yáng)等[21]的研究成果，將這個(gè)常數(shù)定義為：

(7)

基于式(7)，湍渦黏性系數(shù)計(jì)算公式為：

(8)

式(8)考慮了湍動(dòng)能產(chǎn)生于耗散的平衡系數(shù)NP2D(式中fne代表非均衡湍流識(shí)別函數(shù))。該系數(shù)的引入主要為了在產(chǎn)生項(xiàng)過(guò)大時(shí)將渦黏性系數(shù)的確定方式指向剪切應(yīng)力輸運(yùn)模型，即將該系數(shù)與湍動(dòng)能直接聯(lián)系起來(lái)。

2 算例驗(yàn)證

2.1 湍流平板

采用湍流平板對(duì)湍流模型進(jìn)行校驗(yàn)是發(fā)展新模型的第一步。本文采用Wieghardt等[22]發(fā)布的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比對(duì)。計(jì)算的工況為：Ma=0.2，平板板長(zhǎng)設(shè)定為1，因此其特征雷諾數(shù)Re=6×106，采用3套不同粒度的網(wǎng)格進(jìn)行模擬，網(wǎng)格規(guī)模如表2所示。從圖5給出了Rex=2.6×106處的速度型?？梢钥闯觯履Ｐ驮诓煌木W(wǎng)格下均能給出準(zhǔn)確的湍流壁面率，證明新模型能夠給出湍流邊界層正確的速度結(jié)構(gòu)。

表2 平板計(jì)算使用網(wǎng)格Table 2 Computational grids

2.2 DLR-F15三段翼型

二維多段翼型 DLR-F15是德國(guó)宇航院在 2005 年試驗(yàn)的一個(gè)三段翼型，該翼型是低噪聲增升裝置項(xiàng)目[23](Low Noise Exposing Integrated Design for Start and Approach ，LEISA)的研究對(duì)象。DLR 在荷蘭 DNW-KKK 風(fēng)洞對(duì) DLR-F15 多段翼型進(jìn)行了多個(gè)馬赫數(shù)和雷諾數(shù)的測(cè)試，并且公開(kāi)了部分?jǐn)?shù)據(jù)。試驗(yàn)狀態(tài)為：Ma=0.15，Re=2.094×106(參考弦長(zhǎng)為c=0.6 m)，T=525R(R為蘭氏溫度單位)，α=6°。

(a)Coarse

(b)Medium

(c)Fine

圖5 湍流平板的速度型
Fig.5 Velocity profiles of turbulent plates at two positions

從圖6可以看出，新模型對(duì)于縫翼、主翼以及襟翼的峰值模擬精確，而SST模型過(guò)高地預(yù)測(cè)了前緣縫翼和主翼的峰值。

圖6 DLF-F15多段翼壓力分布Fig.6 Pressure distribution of DLR-F15

2.3 NACA4412翼型

NACA4412翼型[24]在特定的迎角下，后緣會(huì)出現(xiàn)一個(gè)穩(wěn)定的分離泡，對(duì)于分離泡的起始點(diǎn)、分離區(qū)域的速度型以及再附點(diǎn)的預(yù)測(cè)一直都是一個(gè)具有挑戰(zhàn)性的工作。計(jì)算的工況為：Ma=0.09，基于單位弦長(zhǎng)的特征雷諾數(shù)Re=1.52×106，迎角α=13.79°。為了簡(jiǎn)化起見(jiàn)，采用一套y+小于1的網(wǎng)格進(jìn)行計(jì)算，網(wǎng)格總量為57 921。遠(yuǎn)場(chǎng)邊界約為翼型弦長(zhǎng)的100倍。挑選6個(gè)站位進(jìn)行速度型取樣，分別為x/c=0.6753、x/c=0.7308、x/c=0.7863、x/c=0.8418、x/c=0.8973和x/c=0.9528。

從速度型上(圖7)可以看出，新模型由于引入了對(duì)于分離流動(dòng)這種非平衡效應(yīng)的感知公式，計(jì)算所得速度型與試驗(yàn)值極為接近，而其他模型不同程度地低估了分離區(qū)的回流速度。

圖7 NACA4412翼型不同站位的速度型Fig.7 Velocity type of NACA4412 airfoil at different locations

2.4 二維鼓包的計(jì)算

基于NASA官網(wǎng)發(fā)布的分離預(yù)測(cè)標(biāo)模(二維鼓包)[25]，計(jì)算了其表面壓力分布和表面摩擦力系數(shù)分布。

從壓力分布上可見(jiàn)(圖8)，Abid模型計(jì)算的壓力最大值與新模型的計(jì)算結(jié)果近似，SST模型對(duì)于壓力最大值的模擬偏低。對(duì)于負(fù)壓峰值的模擬，新模型和Abid模型結(jié)果類(lèi)似，SST模型模擬的負(fù)壓峰值偏大。

表面摩擦力曲線(圖9)顯示新模型對(duì)于摩擦力峰值的預(yù)測(cè)好于其他兩個(gè)模型，但是對(duì)于再附點(diǎn)的預(yù)測(cè)偏后。其他兩個(gè)模型也未能很好地貼合。

2.5 DLR-F6翼身組合體

DPW會(huì)議考量的是從空客某構(gòu)型中簡(jiǎn)化而來(lái)的一個(gè)機(jī)翼-機(jī)身-發(fā)動(dòng)機(jī)短艙構(gòu)型(DLR-F6)。風(fēng)洞試驗(yàn)由法國(guó) ONERA S2MA 風(fēng)洞承擔(dān)[26]。本文采用其中的翼身組合體(不含發(fā)動(dòng)機(jī)短艙)模型進(jìn)行計(jì)算，以考察湍流模型對(duì)于復(fù)雜構(gòu)型的計(jì)算精度。計(jì)算狀態(tài)和特征參數(shù)為：Ma=0.75，Re=3×106，T=520R，平均氣動(dòng)弦長(zhǎng)c=0.1412 m，參考面積S=0.07265 m2，半展長(zhǎng)B=0.5877 m，計(jì)算時(shí)定升力系數(shù)CL=0.5。圖10為不同站位的壓力分布計(jì)算結(jié)果，可以看出，新模型在三維構(gòu)型上與其它三個(gè)模型結(jié)果基本相近。

圖8 二維鼓包壓力分布對(duì)比Fig.8 Comparison of pressure distribution of 2D hump

圖9 二維鼓包表面摩擦力系數(shù)對(duì)比Fig.9 Comparison of surface friction coefficient

(a)x/B=0.15

(b)x/B=0.24

(c)x/B=0.33

(d)x/B=0.51

(e)x/B=0.64

(f)x/B=0.84

圖10 DLR-F6機(jī)翼不同站位壓力分布對(duì)比
Fig.10 Pressure distribution comparison for DLR-F6 wing at different stations

3 結(jié) 論

本文以改造的WDF函數(shù)和重新標(biāo)定的Bradshaw常數(shù)為基礎(chǔ)，構(gòu)造一種低雷諾數(shù)兩方程k-ε湍流模型。主要完成了以下工作：

1)采用湍流雷諾數(shù)構(gòu)造了雷諾數(shù)自適應(yīng)因子，修正WDF函數(shù)不能隨著雷諾數(shù)變化的缺陷；

2)基于平板DNS數(shù)據(jù)對(duì)Bradshaw假設(shè)進(jìn)行改進(jìn)，擬合而成了一種自適應(yīng)Bradshaw函數(shù)，使得雷諾應(yīng)力與湍動(dòng)能之間的比例隨流動(dòng)情況不同而進(jìn)行自適應(yīng)調(diào)節(jié)；

3)通過(guò)平板、翼型、二維鼓包及翼身組合體算例的驗(yàn)證，表明新模型的計(jì)算結(jié)果尚可，在帶有分離和激波的流動(dòng)中能夠準(zhǔn)確反映流動(dòng)的非平衡效應(yīng)，具有一定的工程實(shí)用價(jià)值。