舵翼結構氣動熱力耦合數值模擬與實驗研究

張佳明，王文瑞，武宏宇，溫曉東

(1.北京科技大學 機械工程學院， 北京 100083; 2.流體與材料相互作用教育部重點實驗室， 北京 100083; 3.天津航天機電設備研究所 天津市宇航智能裝備技術企業重點實驗室， 天津 300301)

近空間飛行器擁有快速響應、突防能力強等特點，可在較短時間內對目標實施精準打擊，是各個軍事大國研究的重點，對國家安全、空間技術發展乃至科學技術的進步具有極大的戰略意義[1-2]。在近空間飛行器高速飛行過程中，周圍空氣受到激波的強烈壓縮，對飛行器表面產生劇烈的氣動加熱作用。舵翼結構作為飛行器關鍵部件，產生很高的溫度與溫度梯度，使其結構產生變形或破壞，直接影響飛行穩定性，嚴重時造成偏離飛行軌道甚至掉落，影響精確打擊效果[3-4]。因此，亟需開展舵翼結構在氣動熱力結構耦合作用下的響應研究。

在實驗方面，夏吝時[12]對金屬尖化前緣、吳大方[13]對飛行器耐高溫復合材料翼面進行了熱環境下的熱力性能實驗。歐朝等對錐-柱-裙構件在氣動環境下的結構響應開展了實驗，并分析了表面粗糙度對激波的影響[14]。蔡禮港等對舵翼結構在電弧風洞中的燒蝕速率開展了實驗[15]。

大多數研究者在模擬方面僅針對較低馬赫數工況開展，實驗方面基于單純熱環境進行強度分析，或對氣動環境進行流場驗證，而多物理場耦合模擬結果難以得到實驗驗證。本文根據集束射流風洞裝置，基于計算流體力學、結構力學、空氣動力學及傳熱學基本理論，利用ANSYS Workbench建立一體化的多物理場耦合數值模擬框架，得到舵翼結構在近空間飛行環境中的氣動加熱和結構響應過程，并在集束射流風洞裝置中開展實驗，驗證數值模擬結果。

1 實驗裝置與數值模擬

1.1 實驗裝置

本文實驗采用北京科技大學自主研制的集束射流氣動熱環境模擬實驗裝置如圖1所示，其工作原理如圖2所示。裝置由氣體及航空煤油供應系統、燃燒器、射流噴管、實驗艙、引射管道消音塔組成。利用引射系統在實驗艙內形成負壓環境，空氣、氧氣、航空煤油在燃燒器內反應產生的高溫高壓氣體，通過調節燃燒組分控制來流溫度，經過Laval噴管的膨脹和加速作用，在噴管出口形成速度4～6馬赫、總溫≤2 500 K的氣流，來模擬構件在近空間飛行狀態下多物理場耦合作用的結構響應。

圖1 集束射流氣動熱環境模擬實驗裝置圖

圖2 集束射流氣動熱環境模擬實驗原理示意圖

1.2 數值模擬

1) 流場控制方程

流場按照質量、動量和能量守恒方程來計算。

質量守恒方程的一般形式是由歐拉根據達朗貝爾的研究結果得到的，也被稱為連續性方程：

(1)

式(1)中，ρf為流體密度(kg/m3)；t為時間(s)； ▽代表拉普拉斯算子，V為速度矢量(u、v、w)，則

(2)

式(2)中，u、v、w分別指流體微元在x、y、z三個方向的速度分量(m/s)。

流體的動量守恒方程是指動量對時間的變化率與作用在該流體微元上的各種外力之和是相等的，其一般形式為：

(3)

式(3)中，p指作用在所劃分的流體微元上的壓強(Pa)；fx、fy、fz分別指的是流體微元體積力3個方向的分量(N)；τxx、τxy、τxz等九個量為流體微元粘性應力張量的各個分量是由分子間的粘性作用引起的(Pa)。對于牛頓流體，N-S方程中的粘性應力張量的表達式為：

τxx=λ▽·V+2μux

τyy=λ▽·V+2μvy

τzz=λ▽·V+2μwz

τxy=τyx=μ(uy+vx)

τxz=τzx=μ(uz+wx)

τyz=τzy=μ(vz+wy)

(4)

針對高速流動，在空氣與舵翼結構的相互作用過程中，伴隨著各種形式能量的相互轉換，但整個系統的總能量是保持不變的，這一規律即為能量守恒定律。能量守恒方程可表示為：

(5)

(6)

對于可壓縮粘性流動，將方程(1)、(3)和(5)合稱為納維-斯托克斯方程(N-S方程)。N-S方程中有5個方程式，而有ρf、P、V(u、v、w)、T共6個未知量，未知量個數大于方程個數，則方程組無法進行求解，需添加狀態方程，以使方程組封閉。

本文使用完全氣體狀態方程

p=Z(R*/M)ρfT

(7)

式(7)中，M表示氣體摩爾質量(kg/mol)；R*為摩爾氣體常數，R*=8.314 5 J/(mol·K)；Z為壓縮因子，若不計氣體分子之間的作用力及分子的體積，則壓縮因子Z=1。

2) 結構傳熱方程

根據傅里葉定律和能量守恒定律，可建立熱傳導問題的方程，即導熱微分方程：

(8)

式(8)中，T為溫度場(K)；t為時間(s)；ρ為材料密度(kg/m3)；λ為材料的導熱系數(W/(m·k))；C為比熱容(J/(kg·K))；Q為物體內部的熱源強度(W/kg)。

3) 結構熱響應方程

結構受到氣動力和氣動加熱的同時作用，氣動熱使結構表面溫度升高，通過熱量的傳導使結構出現溫度梯度并發生熱膨脹，結構受到約束產生熱應力，與氣動力共同作用，使結構產生變形。結構的響應方程為：

[σ]=[D][ε]

(9)

式(9)中，[σ]為應力矩陣，[D]為結構材料彈性矩陣，[ε]為結構的總應變矩陣。

2 舵翼結構氣動熱力耦合模擬

飛行器在近空間飛行環境中的氣動加熱過程是一個持續的、非穩態的過程，并且是流場、溫度場和結構場相互耦合作用的過程，本文采用流-熱-固耦合解法進行氣動熱力耦合問題的求解。求解模擬過程如圖3所示，同時建立流體域模型和固體域模型，利用基于有限體積法的Fluent求解器，通過求解連續方程、動量守恒方程和能量守恒方程，計算流體域的溫度、壓強、速度及耦合面上的溫度分布。然后通過System Coupling模塊將流體域網格節點上的熱流和壓力數據插值映射到固體域表面網格上，并作為結構場求解的邊界條件。利用Transient Structural求解器計算結構上的溫度、應力、應變與位移分布。然后再此通過System Coupling模塊將結構的溫度和位移場數據插值映射到流體域網格上，并以此作為邊界條件進行流體域的求解，直到達到所需的耦合計算時間。

圖3 流熱固耦合數值求解模擬過程框圖

2.1 噴管流場模擬

本文采用數值模擬和風洞實驗相結合的方法證明研究結果的準確性，因此首先根據裝置建立噴管流場模型，模擬得到噴管出口處的流場參數，再將該流場參數作為舵翼結構多物理場耦合計算的初始邊界條件，模擬舵翼周圍氣動環境和結構響應。

按照實驗裝置實際尺寸在SolidWorks軟件中建立噴管模型，使用ICEM軟件進行網格劃分。為保證計算結果準確高效，對模型進行結構化網格劃分，由于噴管為對稱結構，所以本文只建立噴管的1/2模型，網格數量為31 826，網格質量均大于0.98。噴管入口氣流為氧氣、航空煤油、空氣燃燒后的混合氣體，成分體積分數為O221%、CO211%、H2O 9%、N259%，噴管入口邊界條件為壓力入口，參數如表1所示，出口邊界條件為壓力出口，壁面為絕熱壁面邊界條件，對稱面為對稱邊界條件，采用剪切壓力傳輸k-w湍流模型，噴管網格如圖4所示。

圖4 噴管流體域網格示意圖

噴管流場速度分布云圖如圖5所示，在噴管出口處形成了速度、溫度、壓強穩定的菱形區域，作為實驗區域，并選擇此區域流場參數作為舵翼結構多物理場耦合計算的初始邊界條件，具體參數值如表1所示。

圖5 噴管流場馬赫數分布云圖

表1 噴管流場參數

2.2 舵翼結構建模

舵翼結構三維模型在SolidWorks軟件中建立，尺寸如圖6所示，取前緣處為原點建立坐標軸，下文相關位置信息以此圖作為參考。將三維模型導入ICEM軟件進行結構化網格劃分，得到其流體域和固體域網格，如圖7所示，由于在舵翼結構的外表面會存在邊界層，邊界層處的流場參數變化梯度較大，為能準確模擬舵翼結構與流場之間的相互作用，對邊界層網格進行加密處理。網格總數為831 246，其中舵翼結構網格數為94 826，流體域網格數為736 420。采用ICEM CFD中的Determinant(2×2×2)值判斷網格質量，得到網格質量均大于0.85，遠大于建議值0.1，因此網格質量滿足計算要求。

圖6 舵翼結構示意圖

圖7 舵翼結構網格示意圖

模型邊界條件如圖8所示，流體域壁面邊界選擇無滑移邊界條件，初始壁溫設置為300 K；流體域外邊界設置為壓力遠場邊界條件，選擇理想氣體為自由來流，來流參數選用噴管出口流場參數。為考慮舵翼結構內部傳熱及變形，將舵翼結構與外流場交界面設置為流-固耦合面。舵翼結構的對稱面添加對稱約束，尾部施加固定約束。舵翼結構材料為D6AC高強度合金結構鋼，其材料參數值如表2所示。

圖8 舵翼結構的邊界條件示意圖

表2 D6AC鋼材料參數

2.3 模擬結果

1) 流場模擬結果

在舵翼結構多物理場耦合數值模擬中，設置時間推進步長Δt=0.000 1 s，耦合計算時間為40 s。舵翼結構20 s時外流場溫度和壓力分布云圖如圖9，由于舵翼結構頭部圓弧半徑很小，對來流的阻滯作用較弱，在舵翼結構頭部軸線上x=-0.65 mm處形成了附體激波，壓力由3 495 Pa升至109 681 Pa，氣流速度急劇降低使動能轉換為熱能，使溫度由513 K升至2397 K，激波前后溫度與壓力數據如圖10所示。

圖9 舵翼外流場溫度分布和壓力分布云圖(20 s)

圖10 舵翼結構前緣流場激波溫度與壓力分布曲線

如圖11所示，由于舵翼前緣激波作用明顯，流場升溫劇烈，在1.8 s時間內由494 K升至2 300 K，最終穩定在2 397 K。舵翼結構尾部下端與夾持部分過渡不平滑，發生激波與膨脹波的摻混、引射，以及附面層的分離，貼近壁面區域的氣體低速流動，外部高速流動，形成渦流，在溫度監測區域形成一個局部高溫區。尾部溫度升溫相對緩慢，在10 s時間內由300 K升至1 542 K，最終穩定在1 723 K。

圖11 舵翼外流場溫度變化曲線

2) 舵翼結構模擬結果

如圖12(a)、圖9(a)溫度分布云圖所示，受流場氣動熱作用，舵翼結構前緣駐點溫度最高，沿尾部方向逐漸降低。隨著時間的推進，舵翼結構溫度逐漸升高，且高溫區域不斷向尾部移動，溫度梯度逐漸減小。當t=10 s時，前緣駐點溫度為1 467 K，當t=40 s時，升高到1 723 K。舵翼結構沿x方向溫度分布曲線如圖13(a)所示。

由圖12 (c)、圖9(b)所示，舵翼結構應變分布與溫度分布趨勢相同，前緣駐點應變最大，沿尾部方向逐漸降低。隨著時間的推進，應變隨溫度升高而逐漸升高，且大應變區域不斷向尾部移動，應變梯度逐漸減小。當t=10 s時，前緣駐點應變為17 864 με，當t=40 s時，升高到25 983 με。舵翼結構沿x方向應變分布曲線如圖13(b)所示。

圖12 舵翼結構溫度分布和應變分布云圖

圖13 舵翼結構溫度與應變沿x方向分布曲線

為了分析氣動熱和氣動力對結構響應的影響大小，將氣動熱和氣動力分別作用在舵翼結構上，變形分布如圖14所示。取模擬過程中t=5 s和t=10 s兩個時間點，在氣動熱作用下，舵翼前緣變形由0.418 69 mm增加至0.655 93 mm；在氣動力作用下，舵翼前緣變形由3.095 1×10-5mm 增加至3.149 6×10-5mm。可見，氣動熱作用下結構變形值與增大幅度均遠遠大于氣動力的作用，因此，氣動熱是造成舵翼結構破壞的主要因素。

3 舵翼結構風洞實驗

對集束射流氣動熱環境模擬實驗裝置噴管出口處的流場參數進行測量，將氣體壓力測試排架安裝在噴管出口，測量氣流壓力參數，如圖15所示。

圖14 氣動力和氣動熱引起的舵翼結構變形分布云圖

圖15 流場壓力測試實驗裝置

測量探頭使流場產生激波，測量得到的壓力為激波后的流場總壓，流場實際靜壓值為：

(13)

式(13)中，P為測點靜壓值(Pa)；PT為探頭測量的激波后總壓值(Pa)；M為流場馬赫數[16]。

實驗測得激波后的流場總壓為116 282.94 Pa，計算得到測點靜壓為3 388.46 Pa，相對模擬結果3 494.92 Pa誤差為3.5%。

舵翼結構高溫應變與溫度測量實驗如圖16、圖17所示，將高溫應變片布置在舵翼結構尾部背風面上部、將S型熱電偶布置在背風面下部，將舵翼試件安裝在支架上，通過鎧裝高溫電纜分別接入高溫應變信號采集系統與溫度信號采集系統。

圖16 舵翼結構高溫應變與溫度測量實驗裝置示意圖

圖17 舵翼結構高溫應變與溫度測量實驗裝置

舵翼結構應變、溫度實驗測量結果與模擬曲線如圖18所示，數值模擬與實驗測量得到的溫度數據變化曲線趨勢是一致的，在時間歷程內最大誤差為6.41%。由于實驗條件限制，在高溫應變測量實驗中只得到了前20 s的數據，從曲線中可以看出，數值模擬得到的舵翼結構尾部應變變化曲線與實驗測得的結果是相匹配的，在20 s的時間歷程內最大誤差為9.73%。

圖18 應變、溫度實驗測量與模擬曲線

4 結論

1) 流場在舵翼結構頭部軸線上x=-0.65 mm處形成了附體激波，壓力由3 495 Pa升至109 681 Pa；前緣在1.8 s時間內由494 K升至2 300 K，最終穩定在2 397 K；尾部形成高溫渦流，在10 s時間內由300 K升至1 542 K，最終穩定在1 723 K。

2) 舵翼結構溫度場與應變場規律相似，前緣駐點處最大，沿尾部方向逐漸降低，隨著時間推移整體溫度和應變逐漸升高，在40 s時分別達到1 723 K和25 983 με。

3) 在t=5 s和t=10 s時，氣動熱單獨作用下舵翼前緣變形為0.418 69 mm和0.655 93 mm；氣動力單獨作用下，則為3.095 1×10-5mm 和3.149 6×10-5mm。可見，氣動熱對舵翼結構變形的影響遠大于氣動力的作用。

4) 流場靜壓測量結果為3 388.46 Pa，較模擬結果相差3.5%；舵翼結構溫度與高溫應變測量結果與模擬基本一致，在時間歷程內最大誤差分別為6.41%和9.73%。證明本文使用的多物理場耦合數值模擬方法基本正確，建立的模型基本可行，得到的結果符合科學實驗規律，對未來采用模擬方法對近空間飛行器材料構件的結構響應分析有著重要的意義。