周 威,張 健,胡雪岑

(1.沈陽(yáng)理工大學(xué),遼寧 沈陽(yáng) 110159;2.沈陽(yáng)工學(xué)院,遼寧 沈陽(yáng) 113122;3.遼沈工業(yè)集團(tuán)有限公司,遼寧 沈陽(yáng) 110045)

應(yīng)用計(jì)算流體力學(xué)(CFD)中的有限體積法是目前進(jìn)行飛行器氣動(dòng)仿真研究的主要計(jì)算方法。模型中,關(guān)于模型計(jì)算數(shù)據(jù)的可靠性,一直是仿真計(jì)算的主要問(wèn)題之一。很多學(xué)者就湍流模型進(jìn)行了相關(guān)研究,推導(dǎo)出了很多方程來(lái)模擬特定湍流,目前氣體外流場(chǎng)仿真計(jì)算應(yīng)用較多的模型是k-ε兩方程模型和k-ω兩方程模型[1]。Launder和Spalding在1972年建立了k和ε的兩方程模型用來(lái)計(jì)算湍流流動(dòng)[2],shih等人對(duì)k-ε進(jìn)行了優(yōu)化[3],實(shí)際上這些方程都是一種假設(shè),由于k-ε兩方程模型中ε的耗散率并不是唯一可能的長(zhǎng)度尺度決定變量,因此,Wilcox提出了k-ω兩方程模型[4],Menter后續(xù)通過(guò)將ε加入ω方程,對(duì)k-ω兩方程模型進(jìn)行了優(yōu)化[5],使其在計(jì)算湍流的零壓力梯度和逆壓力梯度方面具有優(yōu)勢(shì),從而有了更廣泛的應(yīng)用。

但是,以上研究?jī)H是針對(duì)模型外形較簡(jiǎn)單且易進(jìn)行實(shí)驗(yàn)對(duì)比的理論,對(duì)于實(shí)際使用還需要進(jìn)行相關(guān)的適應(yīng)性研究。目前,彈藥的空氣動(dòng)力研究,多集中在導(dǎo)彈和火箭彈,關(guān)于炮射修正尾翼彈的相關(guān)研究較少。本文對(duì)湍流模型的邊界層理論進(jìn)行分析,并對(duì)某口徑炮射超聲速修正尾翼彈的氣動(dòng)力仿真數(shù)據(jù)研究,找出Realizable k-ε和SST k-ω兩種湍流方程對(duì)于炮射修正尾翼彈的適應(yīng)性,以下簡(jiǎn)稱Realizable模型和SST模型。

1 邊界層理論及模型網(wǎng)格劃分

由于Realizable和SST模型都是針對(duì)充分發(fā)展的湍流才有效,而對(duì)于充分發(fā)展的湍流,一定要先進(jìn)行邊界層的研究,這樣才能保證對(duì)模型適應(yīng)性研究的正確性。限于研究范圍,本文對(duì)于k-ε兩方程模型的邊界層壁面函數(shù)不作研究。

1.1 邊界層理論

邊界層是垂直于所要研究的壁面一定厚度的區(qū)域,它有別于完全發(fā)展的湍流核心區(qū)。整個(gè)區(qū)域分為黏性底層、過(guò)渡層和對(duì)數(shù)律層,在CFD中第一層的網(wǎng)格節(jié)點(diǎn)一定要在對(duì)數(shù)律層內(nèi),這樣可以應(yīng)用相對(duì)應(yīng)的模型計(jì)算充分發(fā)展的湍流[6-7]。

由于壁面無(wú)量綱速度u+和與壁面法向距離的無(wú)量綱參數(shù)y+有如下關(guān)系式

(1)

(2)

(3)

(4)

其中,u為流體速度,uτ為壁面摩擦速度,τw為彈體表面剪切應(yīng)力,y為垂直于壁面的距離,κ為卡爾曼常數(shù),對(duì)于光滑壁面κ=0.4,E=9.8,E為與壁面粗糙度有關(guān)的常數(shù),壁面粗糙度越高,E值越小。根據(jù)這種數(shù)值關(guān)系,相關(guān)文獻(xiàn)確定了y+值在對(duì)數(shù)律區(qū)域大致的取值范圍為30

1.2 網(wǎng)格高度劃分

1.2.1 理論模型

炮射修正尾翼彈的理論模型如圖1所示。周圍為長(zhǎng)方體流場(chǎng),流場(chǎng)長(zhǎng)為20 m,寬和高都是8 m。

圖1 氣動(dòng)模型簡(jiǎn)化后的炮射修正尾翼彈三維圖

計(jì)算模型關(guān)于第一層網(wǎng)格的選取高度采用如下公式進(jìn)行計(jì)算：

(5)

Cf=0.058Re-0.2

(6)

(7)

(8)

(9)

式中,Cf為邊界層摩擦系數(shù),τw為邊界層剪切應(yīng)力,μ為氣體黏度,Uτ為計(jì)算剪切速度。

經(jīng)過(guò)計(jì)算,第一層網(wǎng)格高度y取值0.02 mm。將以上公式合并可分析出,第一層網(wǎng)格的高度y與無(wú)量綱壁面距離y+、動(dòng)力黏度μ、特征長(zhǎng)度L成正比,與空氣密度ρ、流體速度u成反比。高度隨著壁面距離、動(dòng)力黏度和特征長(zhǎng)度的數(shù)值增加而增加,隨著空氣密度和來(lái)流速度的增加而降低。溫度、壓強(qiáng)、濕度、細(xì)微粒越多,黏度越大[7]。

1.2.2 網(wǎng)格分布理論

網(wǎng)格越致密,得到的解越精確。如果密度過(guò)大,往往導(dǎo)致增大計(jì)算量而結(jié)果精度不變。因此,通常在對(duì)結(jié)果有影響的部分進(jìn)行局部網(wǎng)格加密。該理論也同樣適用于邊界層網(wǎng)格劃分,如果都以0.02 mm來(lái)劃分,那么網(wǎng)格數(shù)量同樣也會(huì)過(guò)大,增大計(jì)算量。采用比率遞增的網(wǎng)格劃分方法,即在給出第一層網(wǎng)格距離后,根據(jù)比率計(jì)算公式進(jìn)行邊界層網(wǎng)格劃分,公式如下:

(10)

Si=S1·i·eR(i-1)

(11)

其中,R為比率,N為邊界層總高度,S1為第一層網(wǎng)格高度,Si為開(kāi)始節(jié)點(diǎn)到第i節(jié)點(diǎn)的距離,i為邊界層節(jié)點(diǎn)數(shù)。模型劃分網(wǎng)格后,邊界層的節(jié)點(diǎn)分布如圖2所示,總網(wǎng)格數(shù)587萬(wàn),并經(jīng)過(guò)網(wǎng)格無(wú)關(guān)性檢測(cè),整體網(wǎng)格劃分符合仿真計(jì)算要求。

圖2 彈丸邊界層劃分后網(wǎng)格圖

2 湍流模型

在時(shí)均連續(xù)方程和雷諾應(yīng)力方程基礎(chǔ)上,目前大多采用兩種湍流模型計(jì)算氣體動(dòng)力:k-ε兩方程模型和 k-ω兩方程模型。

本文主要對(duì)Realizable模型和SST模型進(jìn)行研究。這兩種湍流模型的方程大體相同,一般由五項(xiàng)組成,方程組成用文字表述為:流體單元參量的變化率(change)+對(duì)流輸運(yùn)(convection)=擴(kuò)散輸運(yùn)(diffusion)+產(chǎn)生項(xiàng)(production)-消耗項(xiàng)(destruction)。

2.1 Realizable模型

它是湍流動(dòng)能和湍流動(dòng)能耗散率的方程,用來(lái)描述完全發(fā)展的湍流流動(dòng),分子的黏性忽略不計(jì)。由于兩方程模型k-ε對(duì)各向異性的雷諾應(yīng)力和自由流動(dòng)等情況的計(jì)算存在缺陷,因此,Realizable模型修正了標(biāo)準(zhǔn)方程中可能導(dǎo)致的負(fù)的正應(yīng)力。

本文基于高雷諾數(shù)的大渦流湍動(dòng)力方程,推導(dǎo)出了耗散率修正方程。對(duì)于完全發(fā)展的湍流流動(dòng)有很好的適應(yīng)性[3]。湍動(dòng)能黏度μt計(jì)算式中的系數(shù)Cμ不再是常數(shù),Cμ見(jiàn)公式(16),式中A0為常數(shù)，AS為流場(chǎng)中與角速度有關(guān)的參數(shù),U*為與角速度有關(guān)的時(shí)均轉(zhuǎn)動(dòng)速率張量。相比標(biāo)準(zhǔn)k-ε兩方程,Realizable模型見(jiàn)式(12)～(16)。

(12)

(13)

(14)

(15)

(16)

2.2 SST模型

它是湍流動(dòng)能和湍流動(dòng)能比耗散率的方程,湍流動(dòng)能比耗散率以ω表示,見(jiàn)式(17)。

(17)

由于多年來(lái)Wilcox對(duì)模型進(jìn)行了修改,在兩方程中的ω方程中加入了生產(chǎn)項(xiàng)Pk,提高了模型預(yù)測(cè)自由剪切流的精度,見(jiàn)式(18)(19)。

(18)

(19)

Menter的SST模型還修正了一些常數(shù),加入了湍流黏度的限制條件,對(duì)產(chǎn)生項(xiàng)湍流動(dòng)能生成率P進(jìn)行了修正,保證了迭代的穩(wěn)定性,見(jiàn)式(20)(21)。在ω方程末尾加入了源項(xiàng)Sω,源項(xiàng)是一種k與ε交叉擴(kuò)散項(xiàng),是在ε方程擴(kuò)散項(xiàng)的變換過(guò)程中產(chǎn)生的,提高了計(jì)算的穩(wěn)定性,更加接近實(shí)驗(yàn)值[4],見(jiàn)式(22)(23)。

(20)

(21)

(22)

(23)

3 仿真分析

本文對(duì)彈丸的外彈道空氣的湍流進(jìn)行研究,相關(guān)的計(jì)算參數(shù)見(jiàn)表1。

表1 彈丸仿真參數(shù)

該模型的仿真結(jié)果已與彈丸實(shí)際參數(shù)進(jìn)行比對(duì),誤差范圍在0.1%之內(nèi)。

兩種模型的湍流動(dòng)能云圖如圖3所示,仿真的迭代次數(shù)曲線圖如圖4所示。從圖3中可以看出,湍流動(dòng)能SST的模型擴(kuò)散得更好。圖4殘差趨于一個(gè)穩(wěn)定的波動(dòng)之中。之所以是穩(wěn)定的波動(dòng),主要因?yàn)橐韵氯c(diǎn)判據(jù):1)計(jì)算項(xiàng)全部采用了二階迎風(fēng)格式,離散化后,精度提高;2)邊界復(fù)雜,湍流的流動(dòng)軌跡線互相產(chǎn)生影響,而且極不穩(wěn)定;3)求得的阻力系數(shù)是穩(wěn)定值。所以這種波動(dòng)現(xiàn)場(chǎng)是正?，F(xiàn)象。

圖3 湍流動(dòng)能k云圖

圖4 仿真的迭代次數(shù)曲線圖

由于基于k-ε模型的計(jì)算中,所有數(shù)值在1300步之后都發(fā)生了劇烈跳動(dòng)(這種跳動(dòng)說(shuō)明整體的模型在計(jì)算時(shí),數(shù)值穩(wěn)定性不好),且湍流動(dòng)能耗散率ε的殘差在250步時(shí)就發(fā)生了跳動(dòng),因此,取整個(gè)尾翼座長(zhǎng)度上方40 mm處高度進(jìn)行研究,兩種模型湍流動(dòng)能耗散率如圖5所示,可以看出，Realizable模型的耗散率變化更大,更不穩(wěn)定。與湍流動(dòng)能耗散率相關(guān)的云圖如圖6所示,Realizable模型的耗散率都集中在尾翼處,其余部分沒(méi)有數(shù)值,而真實(shí)的湍流動(dòng)能耗散應(yīng)該發(fā)生在全彈體表面及周圍,因此,不符合實(shí)際要求。由公式(18)～(23)可知,SST模型修正了湍流動(dòng)能生成率P,加入了源項(xiàng)Sω,并使該模型有更豐富的剪切結(jié)構(gòu)計(jì)算,因此,從云圖看其分布較合理。

圖5 在尾翼座上方兩種模型的湍動(dòng)能耗散率對(duì)比圖

圖6 湍流動(dòng)能耗散率云圖

4 結(jié)束語(yǔ)

本文進(jìn)行了兩種湍流模型的仿真計(jì)算,并對(duì)殘差曲線進(jìn)行了分析,特別是不同模型的湍流動(dòng)能和湍流動(dòng)能耗散率進(jìn)行了對(duì)比,發(fā)現(xiàn)SST模型更適合進(jìn)行超聲速炮射修正尾翼彈的氣動(dòng)仿真。

1)從Realizable模型和SST模型的公式對(duì)比中可以發(fā)現(xiàn),湍流動(dòng)能和湍流動(dòng)能耗散率兩個(gè)方程中的湍流動(dòng)能生成率即產(chǎn)生項(xiàng)有區(qū)別,這就造成了在仿真過(guò)程中湍流過(guò)大時(shí),Realizable模型不能完好地進(jìn)行氣流仿真。

2)由仿真計(jì)算結(jié)果可知,SST模型對(duì)于氣動(dòng)仿真的計(jì)算是適合的。而Realizable模型計(jì)算過(guò)程殘差跳動(dòng)較大,殘差云圖結(jié)果不符合實(shí)際。因此,SST模型的計(jì)算結(jié)果可信度更高。

綜合以上結(jié)果,可以認(rèn)為SST模型更適合于本文所研究的超聲速?gòu)椡柰庑螝鈩?dòng)模擬。但是兩種模型都經(jīng)過(guò)了20多年的不斷改進(jìn),其對(duì)于仿真的適應(yīng)性都很強(qiáng),因此,今后需要針對(duì)這兩種模型所涉及的常數(shù)和其他參量進(jìn)行深入研究,找出所要研究問(wèn)題的最好計(jì)算方案。