本發(fā)明屬于高超聲速飛行器熱防護(hù)系統(tǒng)設(shè)計(jì)領(lǐng)域，具體涉及一種可用于高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)的湍流模型。

背景技術(shù)：

臨近空間高超聲速飛行器已成為世界強(qiáng)國(guó)航空航天技術(shù)發(fā)展計(jì)劃的主要目標(biāo)。此類(lèi)飛行器由于需要進(jìn)行長(zhǎng)時(shí)間的高超聲速飛行將產(chǎn)生嚴(yán)重的氣動(dòng)加熱問(wèn)題，因此，飛行器熱防護(hù)系統(tǒng)的設(shè)計(jì)成為高超聲速飛行器發(fā)展的一個(gè)關(guān)鍵技術(shù)難點(diǎn)。通過(guò)準(zhǔn)確預(yù)測(cè)飛行器的氣動(dòng)熱環(huán)境不僅可以顯著提高熱防護(hù)系統(tǒng)的有效性，還能夠在一定程度上增加飛行器的有效載荷。

當(dāng)高超聲速飛行器在較稠密大氣層內(nèi)飛行時(shí)，繞飛行器的高超聲速流動(dòng)一般不再為層流狀態(tài)，此時(shí)飛行器的氣動(dòng)熱環(huán)境預(yù)測(cè)必須考慮湍流效應(yīng)。當(dāng)前計(jì)算流體力學(xué)(cfd)對(duì)高超聲速流動(dòng)中湍流的模擬大多是基于rans(reynolds-averagednavier-stoke)方法，而當(dāng)流場(chǎng)中存在強(qiáng)激波間斷時(shí)，現(xiàn)有湍流模型的模擬將會(huì)受到嚴(yán)重影響。例如，將當(dāng)前應(yīng)用最為廣泛的k-ωsst兩方程湍流模型應(yīng)用于美國(guó)nasa火星實(shí)驗(yàn)室所設(shè)計(jì)的橢圓鈍化再入飛行器(圖1)氣動(dòng)熱環(huán)境預(yù)測(cè)模擬時(shí)，氣動(dòng)熱模擬結(jié)果與試驗(yàn)數(shù)據(jù)相比產(chǎn)生了極大的誤差，頭部駐點(diǎn)區(qū)域誤差甚至超過(guò)了100％(圖2)。這是由于湍流模型中相關(guān)變量輸運(yùn)方程中的生成項(xiàng)為非守恒形式，強(qiáng)間斷附近的高速度梯度值會(huì)引起生成項(xiàng)的過(guò)度增加，從而引起湍流變量如湍動(dòng)能和耗散率等經(jīng)過(guò)激波時(shí)會(huì)發(fā)生非物理的變化，進(jìn)而導(dǎo)致氣動(dòng)熱環(huán)境的預(yù)測(cè)結(jié)果出現(xiàn)嚴(yán)重的誤差。sinha等人(見(jiàn)參考文獻(xiàn)[1]：a.a.pasha,k.sinha,simulationofhypersonicshock/turbulentboundary-layerinteractionsusingshock-unsteadinessmodel,journalofpropulsionandpower28(1)(2012)46-60)通過(guò)研究發(fā)現(xiàn)，若能對(duì)激波后湍流變量非物理增長(zhǎng)的現(xiàn)象進(jìn)行修正，可有效改善氣動(dòng)熱預(yù)測(cè)模擬的精度。然而，sinha等人提出的模型由于引入了需要進(jìn)行積分的邊界層厚度等非當(dāng)?shù)刈兞?，使得其所提出的修正湍流模型僅適用于簡(jiǎn)單外形流場(chǎng)的預(yù)測(cè)模擬，難以應(yīng)用于現(xiàn)代高超聲速?gòu)?fù)雜外形的飛行器。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明提出一種用于高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)的湍流模型，首先確定流場(chǎng)中存在強(qiáng)間斷的區(qū)域，隨后構(gòu)造了一個(gè)湍流變量生成項(xiàng)衰減函數(shù)以消除激波誘導(dǎo)湍流變量非物理的增長(zhǎng)現(xiàn)象，并與k-ωsst湍流模型結(jié)合獲得了一種適用于高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)的湍流模型，實(shí)現(xiàn)高馬赫數(shù)存在強(qiáng)激波間斷流場(chǎng)下提高氣動(dòng)熱預(yù)測(cè)精度的目的。

具體的，本發(fā)明提供的適用于高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)的湍流模型通過(guò)如下步驟得到：

第一步，計(jì)算網(wǎng)格節(jié)點(diǎn)(i,j,k)的光滑因子。

第二步，計(jì)算該網(wǎng)格節(jié)點(diǎn)(i,j,k)的衰減函數(shù)值，確定流場(chǎng)中強(qiáng)間斷的區(qū)域。

第三步，衰減函數(shù)耦合k-ωsst湍流模型，建立用于高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)的湍流模型。

本發(fā)明的優(yōu)點(diǎn)在于：

(1)計(jì)算精度高?，F(xiàn)有的湍流模型(如經(jīng)典的k-ωsst湍流模型)在應(yīng)用于存在強(qiáng)激波流場(chǎng)的氣動(dòng)熱環(huán)境預(yù)測(cè)時(shí)，駐點(diǎn)區(qū)域計(jì)算誤差甚至能達(dá)到100％以上，而在相同條件下本發(fā)明所提出的湍流模型計(jì)算精度顯著提高，誤差能夠降低至10％以內(nèi)。

(2)實(shí)用性強(qiáng)。相比于已有的湍流修正模型(如sinha的修正模型)僅能夠適用于平板等簡(jiǎn)單外形流動(dòng)的模擬，本發(fā)明所提出的基于光滑因子概念的激波檢測(cè)方法對(duì)于復(fù)雜外形鈍頭體飛行器仍能夠自動(dòng)檢測(cè)激波，并結(jié)合衰減函數(shù)實(shí)現(xiàn)高精度氣動(dòng)熱預(yù)測(cè)模擬，實(shí)用性強(qiáng)。

(3)易于融入現(xiàn)代并行化cfd計(jì)算程序。本發(fā)明所提出的用于高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)的湍流模型并未引入邊界層厚度等需要進(jìn)行積分的非當(dāng)?shù)刈兞浚Ｐ椭凶兞烤鶠楫?dāng)?shù)刈兞?，融入現(xiàn)代并行化cfd計(jì)算程序難度較低。

附圖說(shuō)明

圖1為現(xiàn)有技術(shù)中的橢圓鈍化飛行器幾何外形尺寸；

圖2為現(xiàn)有技術(shù)中的橢圓鈍化飛行器對(duì)稱面下壁面熱流分布與實(shí)驗(yàn)數(shù)據(jù)對(duì)比；

圖3為對(duì)于任意一維網(wǎng)格節(jié)點(diǎn)j光滑因子構(gòu)造中模板的選取方式；

圖4為本發(fā)明所提出的衰減函數(shù)在對(duì)稱面與馬赫數(shù)分布云圖對(duì)比；

圖5為本發(fā)明所提出的湍流模型和原始模型對(duì)稱面湍流渦粘性分布云圖對(duì)比；

圖6為本發(fā)明所提出的湍流模型和原始模型表面熱流預(yù)測(cè)結(jié)果對(duì)比。

具體實(shí)施方式

下面結(jié)合附圖和實(shí)施例對(duì)本發(fā)明進(jìn)行詳細(xì)說(shuō)明。

本發(fā)明提出一種用于高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)的湍流模型及其建立方法，所述的建立方法中提出一種適用于復(fù)雜外形的強(qiáng)激波間斷檢測(cè)方法，用于確定流場(chǎng)中存在強(qiáng)間斷的區(qū)域，隨后基于此構(gòu)造了一個(gè)湍流變量生成項(xiàng)衰減函數(shù)以消除激波誘導(dǎo)湍流變量非物理增長(zhǎng)的現(xiàn)象，并與k-ωsst湍流模型結(jié)合獲得了一種適用于高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)的湍流模型。該湍流模型可顯著提高高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)模擬的精度，對(duì)復(fù)雜鈍頭體外形實(shí)用性強(qiáng)并且易于融入現(xiàn)代并行化cfd計(jì)算程序。

本發(fā)明提供的用于高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)的湍流模型的建立方法，包括以下步驟：

步驟1：計(jì)算網(wǎng)格節(jié)點(diǎn)(i,j,k)光滑因子。

為了解決現(xiàn)有湍流模型在強(qiáng)激波附近產(chǎn)生的湍流變量非物理增長(zhǎng)現(xiàn)象，本發(fā)明首先引入了weno格式構(gòu)造中的光滑因子概念以實(shí)現(xiàn)流場(chǎng)中強(qiáng)間斷區(qū)域的自動(dòng)檢測(cè)。與weno格式不同的是，本發(fā)明采用了無(wú)量綱的壓強(qiáng)其中p為流場(chǎng)當(dāng)?shù)貕簭?qiáng)，u∞為自由來(lái)流速度，ρ∞為自由來(lái)流密度，上標(biāo)“-”代表無(wú)量綱)作為檢測(cè)變量，并在數(shù)值計(jì)算中針對(duì)一維空間任意網(wǎng)格節(jié)點(diǎn)j選擇了三個(gè)具有二階精度的插值模板獲得相應(yīng)的光滑因子，經(jīng)理論推導(dǎo)所獲得的三個(gè)插值模板對(duì)應(yīng)的光滑因子的表達(dá)式為：

其中，其中pj為網(wǎng)格節(jié)點(diǎn)j所在流場(chǎng)當(dāng)?shù)貕簭?qiáng)，is1j、is2j、is3j分別為網(wǎng)格節(jié)點(diǎn)j的三個(gè)插值模板1、插值模板2、插值模板3對(duì)應(yīng)的光滑因子。對(duì)于三維問(wèn)題，每一維上的光滑因子以同樣的方式構(gòu)造，將j替換為i和k，即得到相應(yīng)維度上的光滑因子。

在數(shù)值模擬過(guò)程中對(duì)于任意三維網(wǎng)格節(jié)點(diǎn)(i,j,k)，采用公式(1)中的表達(dá)式對(duì)每個(gè)維度計(jì)算三個(gè)插值模板各自對(duì)應(yīng)的光滑因子值，然后對(duì)每個(gè)插值模板內(nèi)的三個(gè)維度光滑因子按公式(2)的方法計(jì)算得到采用該插值模板計(jì)算得到的光滑因子。

每個(gè)插值模板m(m＝1,2,3)的光滑因子值可將三個(gè)維度上的對(duì)應(yīng)光滑因子值分別平方并求和后開(kāi)方得到：

然后，在三個(gè)插值模板的光滑因子中選取最大值is：

is＝max(is1,is2,is3)(3)

步驟2：計(jì)算該網(wǎng)格節(jié)點(diǎn)的衰減函數(shù)值。

構(gòu)造一個(gè)湍流變量輸運(yùn)方程生成項(xiàng)衰減函數(shù)fd：

fd＝1-tanh(camis)(4)

其中，cam為衰減函數(shù)系數(shù)，通過(guò)理論分析激波前后壓強(qiáng)變化規(guī)律確定該衰減函數(shù)系數(shù)為95。

將公式(3)中三個(gè)插值模板的光滑因子中的最大值is帶入公式(4)，得到該網(wǎng)格節(jié)點(diǎn)的衰減函數(shù)值。

若流場(chǎng)中存在強(qiáng)間斷，由本發(fā)明所建立的上述光滑因子三個(gè)維度中最大值is的取值在1左右，而流場(chǎng)越連續(xù)，光滑因子越接近于0。因此，可通過(guò)網(wǎng)格節(jié)點(diǎn)的三個(gè)插值模板中光滑因子的最大值判斷流場(chǎng)內(nèi)強(qiáng)間斷的區(qū)域，進(jìn)而實(shí)現(xiàn)對(duì)強(qiáng)間斷區(qū)域的自動(dòng)檢測(cè)，確定流場(chǎng)中需要對(duì)湍流變量非物理增長(zhǎng)現(xiàn)象抑制的區(qū)域。由于避免了邊界層厚度等非當(dāng)?shù)刈兞康囊?，使得這種基于光滑因子概念的激波檢測(cè)方法適用于任意一般形狀的激波。

步驟3：將所建立的衰減函數(shù)fd耦合原始k-ωsst湍流模型的生成項(xiàng)，以自動(dòng)在強(qiáng)激波間斷區(qū)域抑制湍流變量的非物理增長(zhǎng)，而在其他區(qū)域保持原始湍流模型計(jì)算方式，建立了一種用于高馬赫數(shù)強(qiáng)激波流場(chǎng)氣動(dòng)熱預(yù)測(cè)的湍流模型，輸運(yùn)方程具體表達(dá)式為：

其中，t代表時(shí)間，k為湍動(dòng)能，ω為定耗散率，ρ為當(dāng)?shù)亓鲌?chǎng)密度，u為當(dāng)?shù)亓鲌?chǎng)速度，為流場(chǎng)的剪應(yīng)力張量，下標(biāo)i,j為維度啞標(biāo)并遵循愛(ài)因斯坦求和約定，μ為分子粘性，μt為湍流渦粘性，β^*、σk、β、σω、σω2、γ和a1為k-ωsst湍流模型的系數(shù)，f1為混合函數(shù)。其模型系數(shù)及混合函數(shù)具體表達(dá)式與k-ωsst湍流模型相同(見(jiàn)參考文獻(xiàn)[2]：f.r.menter,two-equationeddy-viscosityturbulencemodelsforengineeringapplications,aiaajournal,32(8)(1994)1598-1605)。

實(shí)施例：選取了nasa火星實(shí)驗(yàn)室所設(shè)計(jì)的橢圓鈍化再入飛行器測(cè)試外形作為實(shí)施例，其幾何外形如圖1所示。其中空氣來(lái)流條件為：馬赫數(shù)6.03，溫度58.6k，壓強(qiáng)為2091pa，壁面溫度300k，攻角40°；圖4為對(duì)稱面馬赫數(shù)與本發(fā)明所構(gòu)造的衰減函數(shù)分布云圖，圖5給出了模型修正前后對(duì)稱面湍流渦粘性分布云圖對(duì)比。為了方便對(duì)比分析，圖4和圖5中虛線上半圖為鈍頭體下半部分流場(chǎng)繞x軸旋轉(zhuǎn)了180°。由圖4可以看出本發(fā)明構(gòu)造的衰減函數(shù)fd在三維模擬中實(shí)現(xiàn)了對(duì)任意形狀激波的準(zhǔn)確追蹤并在激波區(qū)域保持函數(shù)值為0而在其它區(qū)域?yàn)?。在圖5中，虛線上半圖和虛線下半圖分別為原始湍流模型和本發(fā)明所提出湍流模型的計(jì)算結(jié)果，可以看到采用衰減函數(shù)后湍流模型模擬的湍流渦粘性在激波處不再產(chǎn)生大幅度的增長(zhǎng)，而是集中在邊界層內(nèi)增長(zhǎng)，表明本發(fā)明所提出的衰減函數(shù)對(duì)激波處的湍流變量的非物理增長(zhǎng)現(xiàn)象起到了很好的抑制作用并且沒(méi)有影響到其在邊界層內(nèi)的正常發(fā)展過(guò)程。圖6給出了分別采用原始湍流模型、本發(fā)明提出的湍流模型以及層流計(jì)算所得對(duì)稱面下壁面的熱流分布與實(shí)驗(yàn)數(shù)據(jù)對(duì)比。相比于原始的湍流模型，本發(fā)明提出的衰減函數(shù)湍流模型由于消除了激波誘導(dǎo)的湍流變量的非物理增長(zhǎng)，從而大幅改善了湍流模型在頭部駐點(diǎn)區(qū)熱流的預(yù)測(cè)結(jié)果。另外，從層流計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比可以發(fā)現(xiàn)，實(shí)際流動(dòng)在x/l＝0.13(x為流向坐標(biāo)，l為該飛行器長(zhǎng)度位置開(kāi)始發(fā)生了轉(zhuǎn)捩。由于湍流模型本身無(wú)法準(zhǔn)確預(yù)測(cè)高超聲邊界層轉(zhuǎn)捩的過(guò)程，因此修正后的湍流模型所預(yù)測(cè)的轉(zhuǎn)捩位置提前，轉(zhuǎn)捩區(qū)域較短，導(dǎo)致下游湍流區(qū)域熱流預(yù)測(cè)結(jié)果相比于實(shí)驗(yàn)數(shù)據(jù)略微偏高。然而，與原始的湍流模型預(yù)測(cè)結(jié)果相比，本發(fā)明提出的湍流模型熱流預(yù)測(cè)結(jié)果基本位于實(shí)驗(yàn)數(shù)據(jù)誤差帶內(nèi)或接近誤差帶的上限，熱流預(yù)測(cè)結(jié)果顯著優(yōu)于原始湍流模型，顯著提高了氣動(dòng)熱環(huán)境預(yù)測(cè)模擬的精度。