跨流域跨聲速平板繞流計(jì)算研究

姚思齊

(北京航空航天大學(xué) 北京 100191)

高超聲速飛行器的整個(gè)飛行階段涉及連續(xù)流到自由分子流整個(gè)流域范圍，是典型的跨流域流動(dòng)問題，需要考慮局部的稀薄氣體效應(yīng)、可壓縮效應(yīng)、粘性效應(yīng)和非平衡效應(yīng)的共同作用。全面分析跨流域高超聲速流動(dòng)的特點(diǎn)，準(zhǔn)確預(yù)測流動(dòng)的氣動(dòng)力與氣動(dòng)熱，具有重要的理論價(jià)值和實(shí)際意義[1]。

平板繞流作為流體力學(xué)中最經(jīng)典的問題之一，經(jīng)常被作為理論分析的研究對象[2]。BIRD G[3]研究了Knudsen（Kn）數(shù)在0.1～10 的高超聲速平板的氣動(dòng)特性，結(jié)果表明隨著Kn數(shù)減小，與平板碰撞分子增多，平板表面的剪切應(yīng)力也隨之增大。胡遠(yuǎn)等人[4]結(jié)合了計(jì)算流體力學(xué)（CFD）方法和直接模擬蒙特卡洛（DSMC）方法進(jìn)行模擬，研究了超聲速平板繞流的全流域阻力特性。HUANG A等人[5]最早通過求解BGK模型方程，對高超聲速前緣尖端的平板繞流問題進(jìn)行了研究，得到了流場速度、溫度、壓力和密度等流場宏觀量。

本文采用基于SPARTA平臺(tái)開發(fā)的USP[6-7]算法的求解器SPARTACUS[8]，針對典型的高超聲速平板繞流問題進(jìn)行研究，分別以來流馬赫數(shù)（Ma）和Kn數(shù)為參數(shù)，研究來流變化對流動(dòng)傳熱的影響；并針對主要關(guān)注的平板阻力系數(shù)進(jìn)行分析。

1 計(jì)算參數(shù)設(shè)置

平板繞流示意圖如圖1所示，平板的長度為0.1 m，起始位置在歸一化坐標(biāo)的原點(diǎn)處。平板的厚度為0.005 m，溫度為Twall=300 K，完全漫反射邊界條件。由于平板繞流問題是關(guān)于水平方向軸對稱的流動(dòng)，采用對稱邊界條件，僅考慮平板上半部分的流場。流場左側(cè)和上方為入口，來流氣體分子從左側(cè)水平方向過來，右側(cè)為流場出口。來流的工質(zhì)為單原子氣體——?dú)鍤猓瑴囟葹門∞=300 K，分子碰撞模型選擇變徑軟球模型（VSS），粘性指數(shù)（ω）為0.81，參考直徑（dref）為4.11×10-10m。筆者以Knudsen 數(shù)和Mach 數(shù)作為分析對象，Kn為0.02、0.2、2、20，Ma為2、5、10、20，計(jì)算流域從滑移流域到自由分子流，流速從超聲速到高超聲速。本文針對粒子密度較高且宏觀量變化梯度較大的區(qū)域進(jìn)行網(wǎng)格加密，從而能夠?qū)崿F(xiàn)對流動(dòng)傳熱的準(zhǔn)確捕捉。

圖1 高超聲速平板繞流示意圖

2 流場結(jié)果分析

2.1 流場密度分布

圖2為流場無量綱密度數(shù)值模擬的結(jié)果。為了方便展示，這里給出了代表性的計(jì)算結(jié)果，Kn=0.02 和20，Ma=2 和20。在近連續(xù)流超聲速流動(dòng)中，邊界層的影響導(dǎo)致出現(xiàn)斜激波，從平板前緣附近向下游發(fā)展。在邊界層內(nèi)，由于靠近壁面的連續(xù)流無滑移或滑移流域速度滑移條件，從外部無粘流動(dòng)到近表面流動(dòng)會(huì)出現(xiàn)速度和溫度梯度。由于邊界層內(nèi)的粘性效應(yīng)，流體在外部無粘流動(dòng)和粘性區(qū)域之間進(jìn)行質(zhì)量、動(dòng)量和能量的交換。當(dāng)超聲速來流在粘性區(qū)域內(nèi)流動(dòng)時(shí)，流體受到壓縮，前緣產(chǎn)生壓縮波并收斂形成激波。當(dāng)超聲速氣流通過斜激波時(shí)，流線發(fā)生偏轉(zhuǎn)，流體密度增大。

圖2 不同參數(shù)下的平板繞流密度云圖

激波的厚度大致相當(dāng)于分子平均自由程的幾倍。連續(xù)流中的激波在流場中產(chǎn)生不連續(xù)結(jié)構(gòu)，激波前后物理量跳躍。事實(shí)上，在平均自由程的尺度內(nèi)，激波中會(huì)有劇烈的質(zhì)量、動(dòng)量和能量的輸運(yùn)，盡管這個(gè)尺度對于連續(xù)流來說太小。在激波內(nèi)，流體處于強(qiáng)非平衡狀態(tài)，流體的熵會(huì)通過激波增大。激波結(jié)構(gòu)內(nèi)部的現(xiàn)象需要從分子動(dòng)力學(xué)理論的角度來描述。

隨著來流Kn數(shù)的增加，激波厚度也隨之增加。在連續(xù)流中，非連續(xù)激波與邊界層之間應(yīng)存在無粘性區(qū)域。對于滑移流動(dòng)，激波開始與整個(gè)平板上的邊界層合并，變得難以區(qū)分。在自由分子流中，不存在典型的激波結(jié)構(gòu)，密度通過壁面擾動(dòng)（特別是來自前緣的擾動(dòng)）連續(xù)變化。在近連續(xù)統(tǒng)區(qū)，激波沿前緣向下游的發(fā)展大致呈線性，而在自由分子流區(qū)，激波沿前緣向下游發(fā)展的區(qū)域大致呈線性，擾動(dòng)影響逐漸減小，形成弧形區(qū)域。低馬赫數(shù)（Ma）時(shí)，擾動(dòng)集中在前緣，隨著馬赫數(shù)（Ma）的增加，擾動(dòng)逐漸向下游和靠近壁面移動(dòng)。激波強(qiáng)度隨Kn數(shù)的增加而逐漸減弱，波后的密度變化明顯減小。

馬赫數(shù)（Ma）可以反映流動(dòng)的壓縮效應(yīng)強(qiáng)弱。對于氣體的流動(dòng)，可壓縮性對氣體的密度、溫度和壓力分布有顯著影響。隨著馬赫數(shù)（Ma）的增加，激波角度減小，激波層變薄，這在滑移流區(qū)中可以明顯觀察到。對于自由分子流區(qū)，在馬赫數(shù)（Ma）為2時(shí)，擾動(dòng)主要集中在前緣，而在馬赫數(shù)（Ma）為20 時(shí)，擾動(dòng)逐漸向壁面靠近并覆蓋整個(gè)平板區(qū)域。激波強(qiáng)度隨馬赫數(shù)（Ma）增大而增大，激波密度隨馬赫數(shù)（Ma）增大而增大。

在馬赫數(shù)（Ma）為2時(shí)，流場的最大密度位于激波后和前緣附近，見圖2（a）。對于滑移流動(dòng)，隨著馬赫數(shù)（Ma）的增加，最大密度位置沿激波逐漸向下游移動(dòng)。這是因?yàn)殡S著馬赫數(shù)（Ma）的增加，激波角減小，激波與邊界層合并的區(qū)域變大。在合并層中，流體的密度和壓力會(huì)由于邊界層內(nèi)質(zhì)量和動(dòng)量的交換而降低。對于自由分子流，隨著馬赫數(shù)（Ma）的增加，最大密度的位置會(huì)向壁面靠近。原因是對于高超聲速流動(dòng)，來流的總溫度遠(yuǎn)高于300 K的壁溫。由于墻體的冷卻作用，密度會(huì)明顯增加，并在壁面附近形成高密度區(qū)域。

為明確投資量的影響，回歸分析多個(gè)港口的固定資產(chǎn)投資和泊位能力的數(shù)據(jù)，得到根據(jù)固定資產(chǎn)計(jì)算泊位能力的公式。然后，基于泊位能力與其他設(shè)施指標(biāo)(如裝卸設(shè)備、存儲(chǔ)能力等)間的關(guān)系，確定投資對港口競爭力的影響。設(shè)泊位長度與其他q個(gè)基礎(chǔ)設(shè)施指標(biāo)間的配比率為g，則單位固定資產(chǎn)投資與基礎(chǔ)設(shè)施指標(biāo)的關(guān)系為

2.2 流場馬赫數(shù)（Ma）分布

流場馬赫數(shù)（Ma）分布的數(shù)值模擬結(jié)果如圖3 所示。氣體流經(jīng)激波后會(huì)被壓縮和減速，近連續(xù)流Kn=0.02 中激波前后的馬赫數(shù)可用斜向理論描述，激波的后面，在遠(yuǎn)離表面的地方有一個(gè)粘性區(qū)域，在靠近平板的地方有一個(gè)低速區(qū)，這就是邊界層。邊界層內(nèi)的流體受外部無粘流動(dòng)的驅(qū)動(dòng)，受到壁面阻力的約束。因此邊界層內(nèi)存在劇烈的動(dòng)量輸運(yùn)，可以看出壁面附近存在較大的速度梯度。

圖3 不同參數(shù)下的平板繞流馬赫數(shù)云圖

在平板末端，邊界約束消失，但外部流動(dòng)與來自邊界層的低速之間仍然存在梯度，因此由于粘性效應(yīng)，產(chǎn)生自由剪切流動(dòng)。后緣附近有3層，即無粘性外流層、邊界層和粘性層。子層與附面層相遇后，出現(xiàn)尾流，這是一種受來自子層的低速無粘流驅(qū)動(dòng)的自由剪切流。邊界層和尾跡的厚度隨著馬赫數(shù)的增大而增大，馬赫數(shù)越大，雷諾數(shù)越大。壁面附近的速度明顯不為零，這是速度滑移現(xiàn)象?；扑俣入SKn數(shù)的增加而增大，速度梯度減小，剪切應(yīng)力減小。

2.3 流場溫度分布

圖4為流場無量綱溫度的數(shù)值模擬結(jié)果?？倻囟热Q于來流馬赫數(shù)，激波前后的總溫度保持不變，這是能量守恒的結(jié)果。因此，來流總溫度可以用來表征整個(gè)流場的總能量。由于激波的減速和邊界層動(dòng)量的耗散，流動(dòng)的動(dòng)能會(huì)轉(zhuǎn)化為焓。對于高超聲速流動(dòng)，壁面附近的靜態(tài)溫度可以高達(dá)數(shù)十倍，顯然，這將在壁面附近產(chǎn)生較大的溫度梯度，并產(chǎn)生巨大的對流熱流。

在滑移流Kn=0.02 時(shí)，仍然可以觀察到典型的激波結(jié)構(gòu)，激波結(jié)構(gòu)從前緣到下游大致呈線性發(fā)展。在Kn=20時(shí)，平板對氣體流動(dòng)的擾動(dòng)向上游傳播，在平板前緣上游相當(dāng)大的區(qū)域溫度分布與來流不同。特別是Ma=20，有兩個(gè)區(qū)域溫度非常高，一個(gè)在靠近壁面處，另一個(gè)在前緣正前方。前者是由于壁面的減速作用引起的，后者是由于前緣的阻礙作用造成的。對于連續(xù)介質(zhì)流動(dòng)，前緣擾動(dòng)在高超聲速流動(dòng)中直接向下游傳播，而對于稀薄流動(dòng)，前緣擾動(dòng)可以像在亞音速連續(xù)介質(zhì)流動(dòng)中一樣向上游傳播。隨著Kn數(shù)的增加，高溫區(qū)向下游延伸，甚至覆蓋整個(gè)平板。隨馬赫數(shù)（Ma）增加，高溫區(qū)向壁面移動(dòng)。

3 平板表面參數(shù)研究

圖5 為平板上表面摩阻系數(shù)的數(shù)值模擬結(jié)果，表面摩阻系數(shù)Cf為

圖5 平板上表面摩阻系數(shù)

隨著Ma1增大，激波角減小，流線偏轉(zhuǎn)更劇烈，從而Ma2增大，但壓縮性或激波強(qiáng)度增加會(huì)使Ma2減小更明顯。因此，流馬赫數(shù)越大，激波后速度梯度越小，摩擦力也就越小。

如圖5 所示，在近連續(xù)流Ma=2 時(shí)，摩阻系數(shù)沿板先減小后增大；而高馬赫數(shù)時(shí)，前緣附近會(huì)有一個(gè)局部最大值。這是前緣局部稀薄效應(yīng)作用的結(jié)果。隨著馬赫數(shù)增加，激波強(qiáng)度增大，密度梯度增大。由局部Kn數(shù)定義可知，前緣附近的局部稀薄效應(yīng)會(huì)使馬赫數(shù)增大，導(dǎo)致速度滑移增大，速度梯度變小。因此，離開稀薄效應(yīng)的非平衡區(qū)域后，摩擦在前緣附近減小，沿平板方向增大。

一般情況下，稀薄條件下的摩阻系數(shù)要大于近連續(xù)流條件下。隨著Kn數(shù)的增大，摩阻系數(shù)沿平板方向的變化更為平緩，這是由于稀薄氣體效應(yīng)的范圍擴(kuò)大到整個(gè)平板。對于自由分子流，平板對流動(dòng)的影響在壁面的任何位置都是相同的，因此摩阻沿平板不會(huì)發(fā)生變化。對于自由分子流動(dòng)，摩阻系數(shù)與來流馬赫數(shù)成反比，這與自由分子流動(dòng)理論給出的結(jié)果是一致的：

圖6為平板上表面的壓力分布，壓力系數(shù)（Cp）為

圖6 平板上表面壓力系數(shù)

一般情況下，壓力系數(shù)隨馬赫數(shù)的增加而減小。如圖6（a）所示，前緣附近存在局部最大值，這也是產(chǎn)生稀薄效應(yīng)的原因。一般來說，稀薄流動(dòng)的壓力系數(shù)要大于近連續(xù)流動(dòng)的壓力系數(shù)。隨著Kn數(shù)的增大，壓力系數(shù)沿平板變化的趨勢逐漸緩和。壓力系數(shù)也與來流馬赫數(shù)成反比，這與自由分子流動(dòng)理論給出的結(jié)果一致：

圖7為平板上表面的熱流系數(shù)，熱流系數(shù)為

圖7 平板上表面熱流系數(shù)

當(dāng)熱流系數(shù)約為摩阻系數(shù)的1/2時(shí)，這不僅適用于任何馬赫數(shù)的稀薄流動(dòng)，而且還適用于近連續(xù)流流動(dòng)（Ma>5）。

平板阻力是工程應(yīng)用中最關(guān)心的氣動(dòng)特性參數(shù)，通過對摩阻系數(shù)沿表面積分得到：

式（7）中，Cf,up和Cf,down分別為上下表面摩阻系數(shù)，對于零攻角平板則有Cf,up=Cf,down。在實(shí)際應(yīng)用中，一個(gè)簡單準(zhǔn)確的阻力計(jì)算公式有非常重要的價(jià)值[9]。對于跨流域跨聲速平板繞流問題，需要同時(shí)考慮馬赫數(shù)和努森數(shù)的影響。針對低速平板繞流，樊菁等人[10]提出了阻力計(jì)算的橋函數(shù)公式：

式（8）中，CD,C和CD,FM分別為摩阻系數(shù)的連續(xù)流和自由分子流解。橋函數(shù)是利用連續(xù)流和自由分子流的氣動(dòng)特性，對過渡流的氣動(dòng)特性進(jìn)行插值預(yù)測。對于連續(xù)流，胡遠(yuǎn)等人[4]根據(jù)分析得到了高超聲速平板繞流的摩阻表達(dá)式：

對于自由分子流，同樣給出了拓展到過渡流域的公式：

把式（9）和式（11）帶入橋函數(shù)式（8）中，就可以得到不同參數(shù)下的平板繞流阻力系數(shù)，具體如圖8所示。研究了馬赫數(shù)從2～20變化的幾組數(shù)據(jù)，與通過USP方法模擬得到的結(jié)果進(jìn)行對比發(fā)現(xiàn)，USP 方法計(jì)算得到的跨流域跨聲速平板阻力系數(shù)與橋函數(shù)的計(jì)算結(jié)果吻合的較好。

圖8 跨流域跨聲速平板阻力系數(shù)

4 結(jié)語

本文基于能夠處理跨流域復(fù)雜流動(dòng)傳熱問題的多尺度粒子USP方法，通過開發(fā)后的SPARTACUS求解器分析了典型的高超聲速流動(dòng)，對跨流域跨聲速條件下的平板繞流問題進(jìn)行了系統(tǒng)性的流動(dòng)特性分析；并結(jié)合固壁表面的摩阻、熱流及壓力系數(shù)進(jìn)行研究。結(jié)果表明：隨著馬赫數(shù)增大，壓縮效應(yīng)導(dǎo)致激波厚度減小，平板表面的摩阻系數(shù)熱流系數(shù)都相應(yīng)的增大；隨著Kn數(shù)增大，分子間碰撞不充分導(dǎo)致更多的氣體分子與壁面發(fā)生碰撞，形成了局部的高溫高壓區(qū)域。結(jié)合阻力系數(shù)計(jì)算的結(jié)果與橋函數(shù)公式進(jìn)行對比，結(jié)果吻合得較好。