一種計算高空滑移流動的快速邊界條件

一種計算高空滑移流動的快速邊界條件

在空間高超速機場的飛行區(qū)域，整個層被覆蓋，飛機路徑被穿過連續(xù)流、滑動流、過渡流和自由分子流的區(qū)域。其高機動高精確打擊的特性及含較大的扁平升力面與尖化前緣等外形特征要求在氣動布局設(shè)計中能夠精確獲取飛行器在上述各個流域中的氣動特性及熱環(huán)境?；屏鲄^(qū)和過渡流區(qū)的流動研究比較困難,在基于地面試驗設(shè)備的環(huán)境中再現(xiàn)這樣的流動條件非常昂貴且在技術(shù)上具有很大的挑戰(zhàn),因此數(shù)值模擬已成為該領(lǐng)域的主要研究手段。目前國內(nèi)外滑移流領(lǐng)域已有的研究成果主要集中在低速微流動領(lǐng)域,高超聲速流動領(lǐng)域則相對較少。在低空時,由于空氣密度相對較高,稀薄氣體效應(yīng)對飛行器的影響不明顯,連續(xù)介質(zhì)模型假設(shè)仍然成立,可采用傳統(tǒng)的NavieR-Stokes(N-S)方程及物面無滑移邊界條件進行模擬;隨著高度的增加大氣密度逐漸降低,當(dāng)飛行器周圍的流動位于滑移流區(qū)域時,分子平均自由程的變大使得飛行器壁面附近的努森層內(nèi)連續(xù)介質(zhì)假設(shè)失效,此時速度滑移與溫度跳躍現(xiàn)象顯著。值得注意的是,從連續(xù)流到稀薄滑移流的變化是一個漸變的過程,飛行器在較低的高度亦可能出現(xiàn)局部滑移現(xiàn)象?；菩?yīng)對飛行器氣動力的影響主要反映在摩擦阻力上。對于低空高雷諾數(shù)下,飛行器的總阻力中壓差阻力占主導(dǎo)地位摩擦阻力相對較小;而高空滑移流雷諾數(shù)較小,飛行器所受的摩擦阻力占總阻力的很大部分,因此滑移效應(yīng)必須考慮。在稀薄流的數(shù)值模擬中DSMC(蒙特卡羅直接模擬)是常用的且比較精確的方法,但對于氣體密度相對較高的滑移流及連續(xù)流情況,DSMC計算耗費巨大,采用求解N-S方程加滑移邊界條件的方法可提高計算效率,同時可以給出相當(dāng)精度的模擬結(jié)果。1計算方法1.1t+fx+gy案的+表達直角坐標(biāo)系下,二維可壓縮N-S方程為?Q?t+?F?x+?G?y=?Fv?x+?Gv?y(1)?Q?t+?F?x+?G?y=?Fv?x+?Gv?y(1)方程組的具體形式見文獻。對上述控制方程組運用基于單元中心型的有限體積法求解,無粘通量采用M-AUSMPW+格式,通過MUSCL方法選用Vanleer平均限制器進行空間重構(gòu)達到二階精度,粘性通量采用二階中心格式離散,時間推進為改進的隱式LU-SGS方法。1.2up-gokcen邊界條件目前數(shù)值模擬中常用的滑移邊界條件主要有Maxwell滑移邊界條件及其改進形式,如Gokcen條件及Lockerby條件。文獻指出Gokcen條件所得到的結(jié)果優(yōu)于經(jīng)典的Maxwell條件,但在迭代求解過程中每個時間步都需要確定離開壁面一個分子平均自由程處的速度與溫度,其計算量的增加導(dǎo)致Gokcen條件難以在多區(qū)的復(fù)雜外形網(wǎng)格中的應(yīng)用;Lockerby條件雖然應(yīng)用了壁面函數(shù)方法提高了努森層內(nèi)的模擬精度但是其給出的氣動力/熱結(jié)果并不優(yōu)于Maxwell條件所得。因此,從適用范圍及精度來看經(jīng)典的Maxwell滑移邊界條件具有重要的工程應(yīng)用價值。Maxwell滑移邊界條件中速度滑移和溫度跳躍由下列表達式給出:Vs=A(2-σσ)λ?Vx?n+34μρΤ?Τ?x(2)Vs=A(2?σσ)λ?Vx?n+34μρT?T?x(2)Τ0-Τw=2-αα2γ(γ+1)prλ?Τ?n|0(3)T0?Tw=2?αα2γ(γ+1)prλ?T?n∣∣0(3)式中n代表垂直壁面的方向,x代表與壁面相切的方向,μ為氣體的粘性系數(shù),γ為氣體比熱比,pr為普朗特數(shù)。常系數(shù)A、動量調(diào)節(jié)系數(shù)σ和能量適應(yīng)系數(shù)α由壁面的物理特性及氣體的屬性決定,一般需要實驗來確定。分子平均自由程λ由式(4)給出:對于等溫壁情況,可以忽略式中的溫度梯度項得到滑移速度的簡化形式。2滑移邊界條件的求解應(yīng)用Maxwell滑移邊界條件時,往往直接按式(3)及式(5)的形式進行數(shù)值離散處理,即首先計算壁面法向的速度梯度及溫度梯度然后得到速度滑移和溫度跳躍量,為方便起見,將此數(shù)值處理方法定義為梯度法。計算中發(fā)現(xiàn),上述梯度法在網(wǎng)格尺度較密的情況下會出現(xiàn)計算發(fā)散的問題。下面以速度滑移邊界條件為例對梯度法進行具體分析:如圖1所示,已知某個n時間層的流場信息,求解控制方程可得到n+1時間層的物理量,此時壁面第一層網(wǎng)格中心處切向速度Vn+1xn+1x為已知量,但壁面處的物理量都為未知量。為了求得滑移速度Vn+1sn+1s,首先要確定式(5)中的λ,而壁面處n+1時間層的密度和溫度未知,考慮兩個相鄰時刻分子平均自由程變化不大,故可取n時間層計算。其次,因Vn+1sn+1s為待求量,在壁面法向網(wǎng)格垂直于壁面的假設(shè)下,n+1時間層的速度梯度通常可由下式得到?Vx?n|n+10=2(Vn+1x-Vns)d(6)?Vx?n∣∣n+10=2(Vn+1x?Vns)d(6)求得速度梯度之后,通過式(5)即可得到滑移速度Vn+1s。實際上,若將式(6)代入式(5),則上述梯度法求解過程等價于下面式子Vn+1s=2A(2-σσ)λd(Vn+1x-Vns)(7)式(7)對整個時間推進來說,其等價于雅克比迭代式,須滿足一定的條件才能保證解的存在性,其收斂性條件為對于溫度跳躍邊界條件,梯度法相應(yīng)的收斂性條件為因系數(shù)A、動量調(diào)節(jié)系數(shù)σ和能量適應(yīng)系數(shù)α均為常數(shù),故當(dāng)壁面附近的分子平均自由程λ與網(wǎng)格尺度d之比大于某個臨界值時收斂性條件不滿足。實際上,網(wǎng)格尺度d即第一層網(wǎng)格高度需根據(jù)網(wǎng)格無關(guān)性要求事先給定;分子平均自由程λ與當(dāng)?shù)氐南”〕潭扔嘘P(guān),在稀薄程度相差較大的流場中λ的值相差很大,網(wǎng)格尺度往往因局部網(wǎng)格分布過密而使上述收斂性條件不滿足。下面給出在任意網(wǎng)格密度下均收斂的邊界條件處理方法。實際上式(6)中的速度梯度應(yīng)取當(dāng)前時間層的滑移速度計算,將式(6)中Vns替換為Vn+1s再代入式(5),可以求解出滑移速度其中同理,相應(yīng)的溫度跳躍表達式為式中,T為第一層網(wǎng)格單元中心的溫度為方便起見,本文將此處理方法定義為求解法。分析式(10)及式(12)可知:若Vx與T隨著時間步的推進逐漸收斂,則滑移速度Vs及壁面處的溫度T0也必收斂;當(dāng)θ與ε趨于0時,滑移邊界條件趨于無滑移邊界條件,當(dāng)θ與ε趨于∞時,滑移邊界條件趨于無粘滑移邊界條件。應(yīng)用求解法后,Maxwell滑移邊界條件由第三類邊值問題轉(zhuǎn)換為第一類邊值問題。該方法避免了對梯度項的直接計算,保證了任意網(wǎng)格尺度下解的存在性。3網(wǎng)格結(jié)構(gòu)及結(jié)果分析本文選取了文獻中二維圓柱算例為驗證算例,以考察所給出的邊界條件處理方法的正確性及本文發(fā)展的計算方法的精度。文獻提供了在稀薄氣體領(lǐng)域中比較認(rèn)可的DSMC方法的結(jié)果。計算條件如下:圓柱半徑為0.1524m,氬氣,γ=4/3,pr=2/3,來流溫度200K,馬赫數(shù)10,壁面溫度500K,來流密度及努森數(shù)Kn見表1;粘性系數(shù)公式對應(yīng)DSMC的變徑硬球模型的冪次式,其冪次指數(shù)ω=0.734,參考溫度為1000K。常系數(shù)A、動量調(diào)節(jié)系數(shù)σ和能量適應(yīng)系數(shù)α都取定為1.0,其余參數(shù)見文獻。為了說明網(wǎng)格尺度對邊界條件處理方法的影響及確保計算結(jié)果的網(wǎng)格無關(guān)性,如表2所示,對表1中的兩個狀態(tài)生成了五套具有對比性的網(wǎng)格,各套網(wǎng)格對壁面法向進行了不同程度的加密,周向網(wǎng)格均勻分布。本文中圓柱角度Φ的定義與文獻一致,即以駐點為零點順時針旋轉(zhuǎn)為正。表中列出了各套網(wǎng)格的網(wǎng)格量(圓柱整體-180°至180°范圍)、壁面法向第一層網(wǎng)格高度d。表3給出了case1狀態(tài)從自由來流開始計算,五套網(wǎng)格下梯度法與求解法的收斂情況。對于整體網(wǎng)格情況,梯度法在所有網(wǎng)格下都計算發(fā)散;若僅考慮前體網(wǎng)格(圓柱-90°～90°),梯度法可以在Grid1和Grid2下得到收斂的結(jié)果;求解法在各套網(wǎng)格下均能得到收斂的結(jié)果。對于Grid1和Grid2情況梯度法在整體網(wǎng)格時計算發(fā)散是由于后體網(wǎng)格分布過密的原因。圓柱前體部分氣流受壓縮氣體密度相對較高,而后體部分為尾流區(qū)域氣流膨脹氣體密度相對稀薄,按前述分析,氣體越稀薄則分子平均自由程λ越大,因此圓柱后體需要較大的網(wǎng)格尺度才能滿足收斂性條件。表中的結(jié)果表明,梯度法的應(yīng)用受到了網(wǎng)格尺度的限制。圖2為Case2狀態(tài)下Maxwell滑移邊界條件應(yīng)用求解法得到的表面摩阻系數(shù)與熱流系數(shù)分布在不同網(wǎng)格下的對比。因?qū)ΨQ只給出了圓柱表面0°至180°區(qū)間的結(jié)果。由圖可知,隨著網(wǎng)格的加密,表面摩阻系數(shù)與熱流系數(shù)逐漸接近收斂解;圓柱表面后體即Φ=90°之后的部分各套網(wǎng)格所對應(yīng)的值都相差較小,可認(rèn)為Grid1為結(jié)果收斂所需的網(wǎng)格,而Φ=90°之前的前體部分到Grid3才為收斂網(wǎng)格。相對于后體,圓柱表面前體需要較密的網(wǎng)格才能達到收斂。圖中的結(jié)果表明,在滑移流區(qū)應(yīng)用求解法可以得到網(wǎng)格無關(guān)的計算結(jié)果。圖3給出了Case1及Case2狀態(tài)下滑移邊界條件應(yīng)用求解法及無滑移邊界條件所得到的摩阻系數(shù)及熱流系數(shù)分布與文獻滑移邊界條件及DSMC所得結(jié)果的對比。Case1狀態(tài)處于連續(xù)流區(qū),由圖可知:此狀態(tài)下本文滑移邊界條件的結(jié)果與無滑移邊界條件的結(jié)果相差較小,與文獻的滑移邊界條件及DSMC的結(jié)果都非常接近。Case2狀態(tài)位于滑移流區(qū),從圖可以看出:對于駐點附近之外的區(qū)域,無滑移邊界條件所得的摩阻系數(shù)分布明顯高于DSMC所得結(jié)果;整個圓柱表面無滑移邊界條件所得的熱流系數(shù)分布與DSMC的結(jié)果相比有較大的差異;滑移邊界條件所得的摩阻系數(shù)及熱流系數(shù)分布都與DSMC的結(jié)果比較接近。兩個狀態(tài)下,本文應(yīng)用滑移邊界條件所得的表面摩阻系數(shù)及熱流系數(shù)分布與文獻的結(jié)果都吻合較好。4天飛機模型中的滑移/無滑移邊界條件針對文獻所述的空天飛機模型,本文給出了馬赫數(shù)15,攻角20°,飛行高度60km-90km的情況下Maxwell滑移邊界條件應(yīng)用求解法所得到的模擬結(jié)果。由于實驗?zāi)Ｐ统叽巛^小與實際的飛行器大小有所差距,本文計算中將原始模型放大了10倍,此時空天飛機的長度L=2.9m。計算方法如前所述,按層流狀態(tài)考慮,壁溫1000K,粘性系數(shù)由Surtherland公式得到。圖4所示為空天飛機的表面及空間網(wǎng)格,網(wǎng)格總量為300萬左右,壁面第一層網(wǎng)格高度滿足網(wǎng)格無關(guān)性要求。計算坐標(biāo)系如圖4所示,x=0m的位置位于機身轉(zhuǎn)折處。圖5給出了空天飛機在90km處滑移/無滑移邊界條件的表面及對稱面上的壓力分布。從頭部激波的形狀與位置來看,兩種邊界條件的結(jié)果基本無差異;從機身表面的壓力分布及機翼上的壓力峰值來看,兩種邊界條件所得的結(jié)果差異較小。由此可見滑移效應(yīng)對壓力的影響不明顯。圖6給出了不同高度下空天飛機模型對稱面的摩擦阻力系數(shù)分布。由圖可知,60km與70km時,無論是處于迎風(fēng)的下表面(lower)還是處于背風(fēng)的上表面(lower)滑移邊界條件得到的結(jié)果與無滑移邊界條件的結(jié)果十分接近;隨著高度的增加兩者的差距越來越大,90km處背風(fēng)的上表面無滑移邊界條件所對應(yīng)的摩擦阻力系數(shù)明顯大于滑移邊界條件的值,其中局部位置的偏差高達400%。對比可知下表面摩擦阻力系數(shù)的差異隨高度的變化相差不大,相對上表面而言,由于下表面氣流受壓縮氣體的密度相對較高,因此滑移效應(yīng)對它的影響不顯著。圖7給出了不同高度下x=1m處的熱流系數(shù)分布。由圖可知,在相對較低的高度,滑移邊界條件與無滑移邊界條件所得結(jié)果較為一致,隨著高度的增加,兩者所得結(jié)果的差距逐漸增大。同時也可看出空天飛機沿翼展方向翼前緣的氣動加熱最為嚴(yán)重。5數(shù)值模擬結(jié)果(1)針對高空高超聲速滑移流動的數(shù)值模擬,本文給出了一種在任意網(wǎng)格密度下均收斂的邊界條件處理方法,解決了Maxwell滑移邊界條件在應(yīng)用時出現(xiàn)的迭代計算發(fā)散的問題。二維