帶有剪切應力輸運性質的k-ε兩方程湍流模型構造與應用

2019-08-21 10:30:46陳利麗

空氣動力學學報 2019年3期

關鍵詞：模型

王偉,張揚,陳利麗

(1.航空工業第一飛機設計研究院,西安 710089;2.西安交通大學機械結構強度與振動國家重點實驗室,西安 710049)

0 引言

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

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

新模型通過充分發展的湍流平板、NACA4412翼型、二維鼓包以及DLR-F6翼身組合體的計算來驗證計算精度。

圖1 基于DNS和標準k-ε模型求得的渦黏性系數μtFig.1 Eddy-viscosity coefficient got by DNS and standard k-ε turbulence model

1 湍流模型演化思路

1.1 k-ε兩方程湍流模型

忽略掉浮力項的影響，標準不可壓縮k-ε湍流模型的形式為：

(1)

(2)

其中:ρ為密度，Uj為速度矢量，k為湍動能，μ為層流黏性系數，μt為湍流黏性系數，Pk為湍動能產生項，ε為湍動能耗散率，t為時間變量，xj為空間張量，Cε1、Cε2為常數。

Pk=μtΩ2

(3)

1.2 WDF函數的構造

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

(4)

由于受無滑移邊界條件的約束，式(4)中的速度脈動右端均無常數項，同時，聯合連續方程可得b1=0。通過式(4)，可以得到近壁附近湍動能和湍動能耗散率的階數為：

(5)

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

(6)

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

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

根據觀察，Zhang等人[16]以因子Rek/Reτ來表征雷諾數的變化，通過構造兩個不同的函數fμ1和fμ2，基于此形成新的WDF函數：

(5)

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

(a)Change factor with the change of y+

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

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

1.3 改進的湍渦黏性系數計算方式

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

μt=Cμfμk2/ε

(6)

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

(7)

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

(8)

式(8)考慮了湍動能產生于耗散的平衡系數NP2D(式中fne代表非均衡湍流識別函數)。該系數的引入主要為了在產生項過大時將渦黏性系數的確定方式指向剪切應力輸運模型，即將該系數與湍動能直接聯系起來。

2 算例驗證

2.1 湍流平板

采用湍流平板對湍流模型進行校驗是發展新模型的第一步。本文采用Wieghardt等[22]發布的實驗數據進行比對。計算的工況為：Ma=0.2，平板板長設定為1，因此其特征雷諾數Re=6×106，采用3套不同粒度的網格進行模擬，網格規模如表2所示。從圖5給出了Rex=2.6×106處的速度型。可以看出，新模型在不同的網格下均能給出準確的湍流壁面率，證明新模型能夠給出湍流邊界層正確的速度結構。

表2 平板計算使用網格Table 2 Computational grids

2.2 DLR-F15三段翼型

二維多段翼型 DLR-F15是德國宇航院在 2005 年試驗的一個三段翼型，該翼型是低噪聲增升裝置項目[23](Low Noise Exposing Integrated Design for Start and Approach ，LEISA)的研究對象。DLR 在荷蘭 DNW-KKK 風洞對 DLR-F15 多段翼型進行了多個馬赫數和雷諾數的測試，并且公開了部分數據。試驗狀態為：Ma=0.15，Re=2.094×106(參考弦長為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可以看出，新模型對于縫翼、主翼以及襟翼的峰值模擬精確，而SST模型過高地預測了前緣縫翼和主翼的峰值。

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

2.3 NACA4412翼型

NACA4412翼型[24]在特定的迎角下，后緣會出現一個穩定的分離泡，對于分離泡的起始點、分離區域的速度型以及再附點的預測一直都是一個具有挑戰性的工作。計算的工況為：Ma=0.09，基于單位弦長的特征雷諾數Re=1.52×106，迎角α=13.79°。為了簡化起見，采用一套y+小于1的網格進行計算，網格總量為57 921。遠場邊界約為翼型弦長的100倍。挑選6個站位進行速度型取樣，分別為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)可以看出，新模型由于引入了對于分離流動這種非平衡效應的感知公式，計算所得速度型與試驗值極為接近，而其他模型不同程度地低估了分離區的回流速度。

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

2.4 二維鼓包的計算

基于NASA官網發布的分離預測標模(二維鼓包)[25]，計算了其表面壓力分布和表面摩擦力系數分布。

從壓力分布上可見(圖8)，Abid模型計算的壓力最大值與新模型的計算結果近似，SST模型對于壓力最大值的模擬偏低。對于負壓峰值的模擬，新模型和Abid模型結果類似，SST模型模擬的負壓峰值偏大。

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

2.5 DLR-F6翼身組合體

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