王琛, 田振國, 沈振興

(燕山大學(xué) 河北省重型裝備與大型結(jié)構(gòu)力學(xué)可靠性重點實驗室, 河北 秦皇島 066004)

0 引言

飛行器做再入運動時,將以7～8 km/s的速度穿越地球大氣層[1],飛行器前部將形成弓形激波。激波內(nèi)的氣體被強烈壓縮,氣體壓力作功導(dǎo)致氣體內(nèi)能提升,在高能量環(huán)境下氣體發(fā)生熱電離、形成高溫等離子體氣流,飛行器頭部承受高溫等離子體氣流的劇烈作用。除流體荷載與高溫作用外,飛行器頭部防護罩在高溫下也與等離子體發(fā)生反應(yīng)[2],該防護罩承受了熱作用、荷載作用和化學(xué)作用,并在這三者的耦合作用下發(fā)生損毀。此外,在高速氣流中,等離子體化學(xué)反應(yīng)物存在尚未反應(yīng)即被輸運出流場的現(xiàn)象,稱為化學(xué)非平衡。上述現(xiàn)象表明考慮化學(xué)非平衡-流體-固體-熱耦合后,等離子體流體對于飛行器的流體與固體(簡稱流固)耦合作用相較于熱完全氣體將發(fā)生改變。

為了考慮這種作用并正確設(shè)計防護罩,研究人員采用了多種方法研究飛行器的等離子體流場分布。歐洲航天局的實驗數(shù)據(jù)證明了飛行器以高速在大氣層中飛行時[3]激波處能量被激發(fā),此時應(yīng)采用真實氣體描述熱力學(xué)規(guī)律。樂嘉陵等[1]對前人的經(jīng)典理論做了整理,給出了再入化學(xué)反應(yīng)繞流的基本方程組,其總結(jié)的多溫度熱化學(xué)不平衡模型成為近年來高超聲速等離子體流場研究所使用的基本理論。韋毅[4]基于單溫度和雙溫度熱化學(xué)不平衡方程組研究了再入飛行器的等離子體鞘套,計算了等離子體流場、熱場、物質(zhì)場與電子數(shù)分布,分析了熱化學(xué)非平衡對于物種分布和電子數(shù)的影響,計算結(jié)果表明在球頭區(qū)域熱化學(xué)不平衡效應(yīng)明顯。華彩成[5]基于鈍錐實驗?zāi)Ｐ瓦M行了數(shù)值模擬,并分析了化學(xué)反應(yīng)模型對等離子體流場的影響,發(fā)現(xiàn)不同的化學(xué)反應(yīng)模型對流場影響不大,是否考慮熱化學(xué)非平衡對流場影響較大。常用的化學(xué)反應(yīng)主要由Park[6]、Dunn等[7]和Gupta等[8]給出,由于等離子體反應(yīng)復(fù)雜,當(dāng)前對于化學(xué)反應(yīng)模型的選取沒有形成統(tǒng)一的標(biāo)準(zhǔn)。流固耦合作用計算對于高超聲速飛行器設(shè)計尤為重要,盡管模擬方法取得了重要進展,但考慮真實氣體效應(yīng)的多物理場耦合問題仍然是一個挑戰(zhàn)[9-12]。Sch?fer[13]在高焓條件下進行了飛行器流固耦合實驗,研究了通用碳化硅復(fù)合材料的熱變形,并與流固耦合仿真結(jié)果[13-15]進行了對比,計算了超聲速流固耦合問題,但在速度上仍然較低。在計算方法上,聶濤等[16]也對高超聲速飛行器前緣進行了流固耦合計算,并與圓管繞流實驗作對比,采用緊耦合的方法同時利用有限元法與有限體積法進行計算,目前該方法是飛行器流固耦合計算的常用辦法。緊耦合相對于松耦合[17]計算數(shù)值上更加精確,減小了異步迭代誤差。李佳偉[18]基于熱完全氣體研究了高超聲速流-固-熱多場耦合數(shù)值模擬的方法,認(rèn)為氣體物性參數(shù)會導(dǎo)致氣動力的改變。秦啟豪[19]建立了基于計算流體力學(xué)(CFD)/計算結(jié)構(gòu)力學(xué)(CSD)的多物理場求解策略,計算了熱完全氣體環(huán)境下進氣道的流固耦合作用,研究了流場在進氣道變形下的振蕩。Dennis等[20]進行了高超聲速流固耦合問題的薄板風(fēng)洞實驗,拓展了該領(lǐng)域流固耦合問題非常有限的實驗數(shù)據(jù)。Quan等[21]采用CFD方法研究了高超聲速熱完全氣體環(huán)境下板構(gòu)件的氣動熱彈性特性,認(rèn)為熱荷載與氣動荷載的組合可能會對飛行器結(jié)構(gòu)產(chǎn)生不利的作用。

截至目前,尚無一種可以綜合考慮等離子體真實氣體效應(yīng)與流固耦合計算的辦法,相關(guān)研究表明真實氣體效應(yīng)[9-12]與化學(xué)非平衡效應(yīng)[4-5]對飛行器流固耦合影響較大。

本文從考慮等離子體真實氣體效應(yīng)的角度,建立等離子體流固耦合方程組,在經(jīng)典流固耦合方程和瞬態(tài)固體力學(xué)平衡方程的基礎(chǔ)上,增加再入流場熱化學(xué)非平衡繞流方程組作為真實氣體方程;基于建立的模型計算等離子體對于再入飛行器的流固耦合效應(yīng),發(fā)現(xiàn)在較高的飛行速度下最大黏性力位置發(fā)生遷移,并給出了等離子體熱化學(xué)流固耦合的計算結(jié)果。

1 考慮等離子體的流固耦合模型

考慮到多物理場問題的非線性和復(fù)雜性,為利于問題的建模和求解,采用文獻[4]的假設(shè)。采用Gupta等[8]所提出的等離子體化學(xué)反應(yīng)模型,如表1所示。表1中第三碰撞體為催化反應(yīng)進行的物質(zhì),催化模型見文獻[1]。

表1 化學(xué)反應(yīng)式Table 1 Chemical reaction

由于飛行速度最高計算至馬赫數(shù)22,采用7組分18反應(yīng)模型,若要計算至更高馬赫數(shù),則必須考慮更多物種,本文采用的反應(yīng)參與物種為N2、O2、N、O、NO、NO+和e-。除了化學(xué)反應(yīng)的描述外,還需要各個組分的熱力學(xué)描述,采用文獻[4]給出的熱力學(xué)表達(dá)式。對于第i個物質(zhì),基于質(zhì)量守恒定律建立物質(zhì)輸運方程[1]:

(1)

流體黏度和熱傳導(dǎo)系數(shù)按照文獻[4]整理的計算方法得到。文獻[1]總結(jié)了等離子體混合物的性質(zhì)。根據(jù)參考文獻[22],采用二維黏性可壓縮流體動力學(xué)方程組描述混合氣體質(zhì)量輸運、動量輸運以及能量輸運。為便于分析熱源,采用熱力學(xué)溫度T描述能量方程:

(2)

式中:p為混合物壓力;I為單位矩陣;K為黏性張量;dφ為二維分析的周向厚度,dφ=1;Cp為恒壓熱容;q為熱傳導(dǎo)通量;Q為流體熱源。在模型中,流體熱源來源于等離子體黏性加熱,壓縮加熱和化學(xué)產(chǎn)熱,基于參考文獻[4],為便于計算,將文獻中的方程修改為統(tǒng)一的熱源項,為

(3)

式中:Qp為壓縮加熱源項;Qvd為黏性加熱源項;Qplasma為化學(xué)產(chǎn)熱源項;T為溫度;μ為混合物流體黏度;hi為物種i的質(zhì)量焓;Di為物種i的擴散系數(shù);ns為物種數(shù);r為半徑。飛行器流固耦合邊界承受高超聲速飛行的流體荷載,流固耦合邊界上的荷載為流體壓力和黏性力之和[23],即

Fsum=[-pI+K]·n

(4)

式中:Fsum為流固耦合邊界氣動荷載;n為流固耦合邊界法向量。

建立固體力學(xué)瞬態(tài)平衡方程[24]為

(5)

式中:S為應(yīng)力張量;F為體積力;us為固體位移場。在計算中,考慮流體對飛行器的作用,但不考慮飛行器的變形對流體的擾動。

2 考慮等離子體的流固耦合算例

采用美國國家航空航天局(NASA)的RAM-C Ⅱ再入實驗飛行器為計算模型,如圖1所示。在流體動力學(xué)方程和固體力學(xué)方程之間建立耦合關(guān)系,同時將等離子體化學(xué)反應(yīng)作為能量源項和物質(zhì)源項參與計算,應(yīng)用COMSOL軟件進行全耦合求解。

圖1 RAM C-Ⅱ的幾何尺寸Fig.1 Geometric dimensions of RAM C-Ⅱ

采用基于有限元法的瞬態(tài)求解方式,根據(jù)參考文獻[4]提供的61 km高空條件,有初始條件:

u=[0 m/s,0 m/s],p=27 Pa,T=254 K

(6)

入口邊界條件為馬赫數(shù)入口(馬赫數(shù)為Ma∞),來流馬赫數(shù)的大小隨著時間的增長而線性增大,出口邊界條件為壓力出口(pout為出口壓力),來流溫度T∞和來流壓力p∞與初始條件相同。

(7)

考慮飛行器壁面為無滑壁,則壁面速度uw為

uw=[0 m/s,0 m/s]

(8)

考慮來流空氣質(zhì)量分?jǐn)?shù)比例中N2占比0.767,O2占比0.233,其余組分為0。為驗證真實氣體模型的正確性,將等離子體流場計算結(jié)果與文獻[4]中的算例作對比。此外,對計算采用的映射網(wǎng)格進行逐次加密,檢驗數(shù)值計算的網(wǎng)格無關(guān)性,如表2所示,其中應(yīng)力為馬赫數(shù)Ma=22下等離子體環(huán)境中固體的最大應(yīng)力。

由表2可知,不同網(wǎng)格劃分程度得到的結(jié)論相差不大,因此采用40×(50+4)網(wǎng)格,即采用指數(shù)型加密的方式,邊界第1層網(wǎng)格厚度為2.5 μm,飛行器邊界網(wǎng)格數(shù)量為40,垂直于邊界向外擴展網(wǎng)格數(shù)量為50。對固體區(qū)域采用結(jié)構(gòu)化網(wǎng)格,為捕捉到應(yīng)力梯度,劃分4層,層厚為10 mm,如圖2所示。

表2 算例檢驗Table 2 Example test

圖2 結(jié)構(gòu)化網(wǎng)格劃分Fig.2 Structured meshing

本文算法計算結(jié)果與文獻[4]的計算結(jié)果相差不大,證明了計算模型考慮真實氣體效應(yīng)的可靠性。電子數(shù)密度因數(shù)值敏感,受算法影響較大,一般在數(shù)量級上沒有差異即可認(rèn)定數(shù)值上相差不大[4]。

3 流固耦合機制和計算結(jié)果討論

研究考慮等離子體真實氣體效應(yīng)的流固耦合方程,得到以RAM C-Ⅱ飛行器為算例的計算結(jié)果。采用熱完全氣體模型和單溫度等離子體模型計算結(jié)果對比,分析等離子體氣體和其組分對于流固耦合計算的影響,采用勻加速再入體速度,不考慮飛行器的運動姿態(tài),計算出等離子體流場、壓力場、溫度場和物質(zhì)場后,對等離子體關(guān)于流固耦合的作用機制進行了分析。

3.1 等離子體流固耦合氣動沖擊和氣動摩擦分析

高層大氣中,空氣的壓力性質(zhì)和黏性力性質(zhì)相對于熱完全氣體發(fā)生較大變化,二者是改變飛行器氣動荷載環(huán)境的關(guān)鍵因素。有邊界氣動壓力Fn和邊界氣動黏性力Fτ分別為

Fn=(Fsum·n)n

(9)

Fτ=Fsum-Fn

(10)

式中:Fsum為飛行器所受的氣動合力。

圖3的計算結(jié)果顯示,飛行器由靜止加速到Ma=22的加速過程中熱完全氣體和等離子氣體兩種氣體模型對于再入飛行器的氣動作用存在差異,其中當(dāng)速度Ma=22時,熱完全氣體和等離子氣體的流固耦合作用差別顯著。

圖3 等離子體與熱完全氣體的氣動力對比Fig.3 Comparison of aerodynamic forces between plasma and thermal perfect gas

對比發(fā)現(xiàn),在Ma=22的高速氣流作用下,等離子體氣動壓力由于考慮了真實氣體作用,相較于完全氣體提升了16.5%,而氣動黏性力較完全氣體而言提升了611.8%,且最大作用力位置不同,最大黏性力位置和最大壓力位置將承受高速高溫氣體最大的摩擦作用和沖擊作用,摩擦位置和沖擊位置對于考慮燒蝕作用的再入飛行器設(shè)計較為重要,可能是最快磨損燒蝕的位置。由圖3可見:在Ma=22的速度下,對于熱完全氣體,摩擦位置位于中部,沖擊位置位于頂部,呈環(huán)狀;反之,對于等離子氣體,摩擦位置位于中上部,沖擊位置位于頂部,且覆蓋整個飛行器前端,表明了飛行器運動時的實際氣體摩擦和沖擊狀態(tài)。以圖4(a)所示飛行器邊界弧長s為橫坐標(biāo),以頭部為起點,計算不同來流馬赫數(shù)下的邊界氣動壓力和邊界氣動黏性力。

由圖4(b)、圖4(c)、圖4(d)可見,由于氣動外形相對固定,且來流方向不變,在馬赫數(shù)為16～22的流速下,最大氣動力位置幾乎不改變。

圖4 等離子體和熱完全氣體的沖擊與摩擦對比Fig.4 Shock and friction comparison for plasma and thermal perfect gas

等離子氣體與熱完全氣體之間的沖擊差異主要與等離子體平均粒子動量有關(guān),等離子體采用了精確的組分計算空氣的平均摩爾質(zhì)量。由于等離子氣體模型中,計算了各個組分黏性的不同,圖4(e)、圖4(f)、圖4(g)中等離子體氣體與熱完全氣體的摩擦位置相差很大,熱完全氣體沒有考慮氣體組分的化學(xué)反應(yīng),僅考慮了流場的輸運作用,因此摩擦位置在馬赫數(shù)為16～22的速度區(qū)間內(nèi)幾乎不發(fā)生遷移。反之,等離子體氣體考慮了化學(xué)平衡的改變,等離子體熱力學(xué)性質(zhì)的改變以及流場的輸運,隨著流速的提高,不同組分的流體黏度逐漸發(fā)生改變,摩擦位置逐漸向鈍體中部遷移。

由圖4(f)可見,黏性力在222 mm處的變化較為顯著,由圖4(e)可見,等離子體模型的計算結(jié)果表明,摩擦環(huán)較平滑且數(shù)值增大了611.8%,表明了等離子體氣體模型所計算的混合物黏性力更大。

3.2 再入流固耦合應(yīng)力分析

飛行器在高速運動時,真實氣體效應(yīng)對其受荷狀態(tài)影響較大。研究真實氣體的流固耦合作用對于正確計算飛行器受荷狀態(tài)有重要的意義。流體荷載所引起的流固耦合應(yīng)力雖不如熱作用引起的應(yīng)力大,但其對飛行器受荷狀態(tài)有著多方面的耦合影響。飛行器外部材料為鋼,其彈性模量為200 GPa,泊松比為0.30,外殼厚度為40 mm。

通過計算得到了Ma=22下熱完全氣體和等離子體下的再入飛行器外部結(jié)構(gòu)鋼殼的應(yīng)力場,圖5(a)為飛行速度Ma=22下等離子體流固耦合應(yīng)力σMises的分布云圖,圖5(b)為飛行速度Ma=22下熱完全氣體流固耦合應(yīng)力σMises的分布云圖。

圖5 流固耦合作用下飛行器的Mises應(yīng)力Fig.5 Mises stress of aircraft under fluid structure coupling

由圖5可見,等離子體流固耦合應(yīng)力與熱完全氣體流固耦合應(yīng)力的分布相似,這與飛行器上的荷載主要是氣動壓力有關(guān),等離子體氣動壓力與熱完全氣體氣動壓力作用形式和作用大小相似。從數(shù)值大小上分析,等離子體流固耦合應(yīng)力比熱完全氣體流固耦合應(yīng)力在Ma=22下的最大值小9.5%,這是由于等離子體氣動壓力更為集中的作用于飛行器尖端,其沖擊環(huán)可近似是為一個沖擊點,而熱完全氣體氣動壓力更集中作用在飛行器頭部靠近尖端的沖擊位置上。計算結(jié)果證明,等離子體氣動作用雖然較大,但實際造成的Mises應(yīng)力更小。

3.3 不同等離子體組分的流固耦合作用

由于等離子體是由多種組分混合而成的混合物,流固耦合作用也可以分解為不同等離子體組分的作用,不同組分在壁面上的組分質(zhì)量分?jǐn)?shù)顯示了壁面承受組分作用的貢獻。對再入飛行器壁面的等離子體組分分析,如圖6所示。

隨著馬赫數(shù)的增加,圖6(a)中壁面氧分子質(zhì)量分?jǐn)?shù)逐漸減小,至Ma=22時大部分離解。圖6(b)中,壁面氧原子質(zhì)量分?jǐn)?shù)逐漸增加,至Ma=22基本覆蓋飛行器前端。圖6(c)中壁面氮分子質(zhì)量分?jǐn)?shù)逐漸減小,至Ma=22時,飛行器前端基本離解完畢,中后端質(zhì)量分?jǐn)?shù)較大,離解程度逐漸減弱。圖6(e)中NO作為N與O的結(jié)合生成物,同時也是電子與NO正離子的反應(yīng)物,在Ma=16時質(zhì)量分?jǐn)?shù)較高,在Ma=22時質(zhì)量分?jǐn)?shù)較低。由圖6(f)和圖6(g)可見,電子和NO正離子的質(zhì)量分?jǐn)?shù)逐漸增大。

圖7為不同來流速度下壁面氣體的熱力學(xué)溫度。由圖7可見,隨著來流馬赫數(shù)的增加,再入飛行器流固邊界熱力學(xué)溫度逐漸上升,導(dǎo)致圖6中流固邊界的空氣發(fā)生了多種等離子體反應(yīng)。再入飛行器承受的沖擊粒子由氧氣和氮氣轉(zhuǎn)為承受氧原子和氮原子。這意味著此時飛行器承受流固耦合作用的物質(zhì)不同。

以馬赫數(shù)22來流下的再入飛行器壁面物種而言,結(jié)合3.1節(jié)得到的馬赫數(shù)22來流速度下的再入飛行器壁面流固耦合作用,如圖8所示:在弧長0～200 mm的位置上,O2已充分反應(yīng)為O,氮氣也充分反應(yīng)為N,少部分O與N發(fā)生結(jié)合生成NO,飛行器尖端主要承受來自這3種物質(zhì)組成的原子氣體的沖擊和摩擦;在弧長200～400 mm的位置上,即飛行器中后部,隨著溫度的減小和流場對于O與N的輸運,發(fā)生了化學(xué)平衡的逆向移動,此時O與N又生成了O2與N2,這部分壁面上主要承受O2與N2組成的分子氣體的沖擊和摩擦。對于電子與NO正離子組成的等離子體,由于主要反應(yīng)物NO的生成量較少,該部分氣體對于壁面的沖擊和摩擦可以忽略,從圖8中可以看出二者質(zhì)量分?jǐn)?shù)分布的一致性,從側(cè)面檢驗了模型的計算正確性。

圖6 不同速度下等離子體組分壁面分布Fig.6 Distribution of plasma components on the wall at different velocities

圖7 不同速度下的壁面熱力學(xué)溫度Fig.7 Wall thermodynamic temperature at different velocities

3.4 真實氣體效應(yīng)對流固耦合的影響

并非所有飛行器都在高溫條件下工作,對于低速飛行器,來流空氣幾乎不發(fā)生化學(xué)反應(yīng),熱力學(xué)性質(zhì)和輸運性質(zhì)較為穩(wěn)定,使用熱完全氣體模型可以得到比較精確的結(jié)果。隨著速度的增大,真實氣體效應(yīng)對飛行器的影響將逐漸加強。

圖8 Ma=22等離子體壁面分布Fig.8 Plasma distribution on boundary in Ma=22

圖9(a)為飛行器尖端的物種變化,圖9(b)為飛行器尖端的氣動荷載變化與兩個氣體模型的差異。由圖9(a)可知,當(dāng)飛行器的速度升至馬赫數(shù)5時,振動能量被激發(fā),此時氧原子開始生成,馬赫數(shù)5后需要使用真實氣體模型才能準(zhǔn)確描述真實氣體中各個組分對飛行器的影響。由圖9(b)可知,等離子體與熱完全氣體間的相對差異在馬赫數(shù)5后趨于穩(wěn)定。

圖9 等離子體氣動荷載影響Fig.9 Influence of plasma aerodynamic load

4 結(jié)論

本文對等離子體關(guān)于再入飛行器的流固耦合作用做了分析,并與熱完全氣體流固耦合作用進行對比,探討了高超聲速飛行器在等離子體中穿行所承受的流固耦合作用。得出以下主要結(jié)論:

1) 等離子氣體相較于熱完全氣體,由于考慮等離子體真實氣體效應(yīng),產(chǎn)生的最大氣動壓力更大,在馬赫數(shù)22的高速來流下增大了16.5%,黏性力增大了611.8%。

2) 等離子氣體相較于熱完全氣體,最大沖擊和摩擦位置將會改變,熱完全氣體在來流速度逐漸增大的過程中,該位置幾乎不變,而等離子體的最大摩擦位置會向鈍體中部移動,表明了高超聲速飛行器考慮真實氣體效應(yīng)的流固耦合作用最大值位置將發(fā)生遷移。

3) 高超聲速飛行器不同位置承受的沖擊和摩擦物質(zhì)不同,隨著來流馬赫數(shù)的增大,等離子體組分分布動態(tài)改變,在小于馬赫數(shù)20的速度下,等離子體反應(yīng)溫度不夠,沖擊和摩擦物質(zhì)主要為O2與N2。在馬赫數(shù)為20～22的速度下,沖擊物質(zhì)由分子向原子轉(zhuǎn)變。在馬赫數(shù)22的來流速度下,沖擊物質(zhì)主要為O與N,還存在有少量NO,此時還有極少數(shù)的NO正離子和電子。以馬赫數(shù)22來流速度下的壁面為例,飛行器頂部承受原子氣體作用,中后部承受分子氣體作用,表明化學(xué)非平衡對流固耦合機制影響較大,將導(dǎo)致流固耦合作用物質(zhì)分布的不同。

4) 當(dāng)飛行器速度超過馬赫數(shù)5時,氣體振動能被激發(fā),等離子體在描述氣體熱力學(xué)和化學(xué)狀態(tài)上更準(zhǔn)確,在計算時必須考慮到這一差異。