等離子體環境下高超聲速飛行器的流-固耦合機制

王琛, 田振國, 沈振興

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

0 引言

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

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

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

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

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

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

表1 化學反應式Table 1 Chemical reaction

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

(1)

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

(2)

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

(3)

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

Fsum=[-pI+K]·n

(4)

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

建立固體力學瞬態平衡方程[24]為

(5)

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

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

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

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

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

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

(6)

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

(7)

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

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

(8)

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

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

表2 算例檢驗Table 2 Example test

圖2 結構化網格劃分Fig.2 Structured meshing

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

3 流固耦合機制和計算結果討論

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

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

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

Fn=(Fsum·n)n

(9)

Fτ=Fsum-Fn

(10)

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

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

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

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

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

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

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

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

3.2 再入流固耦合應力分析

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

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

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

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

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

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

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

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

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

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

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

3.4 真實氣體效應對流固耦合的影響

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

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

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

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

4 結論

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

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

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

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

4) 當飛行器速度超過馬赫數5時,氣體振動能被激發,等離子體在描述氣體熱力學和化學狀態上更準確,在計算時必須考慮到這一差異。