舵翼結(jié)構(gòu)氣動熱力耦合數(shù)值模擬與實驗研究

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

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

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

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

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

1 實驗裝置與數(shù)值模擬

1.1 實驗裝置

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

圖1 集束射流氣動熱環(huán)境模擬實驗裝置圖

圖2 集束射流氣動熱環(huán)境模擬實驗原理示意圖

1.2 數(shù)值模擬

1) 流場控制方程

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

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

(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等九個量為流體微元粘性應(yīng)力張量的各個分量是由分子間的粘性作用引起的(Pa)。對于牛頓流體，N-S方程中的粘性應(yīng)力張量的表達式為：

τ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)

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

(5)

(6)

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

本文使用完全氣體狀態(tài)方程

p=Z(R*/M)ρfT

(7)

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

2) 結(jié)構(gòu)傳熱方程

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

(8)

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

3) 結(jié)構(gòu)熱響應(yīng)方程

結(jié)構(gòu)受到氣動力和氣動加熱的同時作用，氣動熱使結(jié)構(gòu)表面溫度升高，通過熱量的傳導使結(jié)構(gòu)出現(xiàn)溫度梯度并發(fā)生熱膨脹，結(jié)構(gòu)受到約束產(chǎn)生熱應(yīng)力，與氣動力共同作用，使結(jié)構(gòu)產(chǎn)生變形。結(jié)構(gòu)的響應(yīng)方程為：

[σ]=[D][ε]

(9)

式(9)中，[σ]為應(yīng)力矩陣，[D]為結(jié)構(gòu)材料彈性矩陣，[ε]為結(jié)構(gòu)的總應(yīng)變矩陣。

2 舵翼結(jié)構(gòu)氣動熱力耦合模擬

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

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

2.1 噴管流場模擬

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

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

圖4 噴管流體域網(wǎng)格示意圖

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

圖5 噴管流場馬赫數(shù)分布云圖

表1 噴管流場參數(shù)

2.2 舵翼結(jié)構(gòu)建模

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

圖6 舵翼結(jié)構(gòu)示意圖

圖7 舵翼結(jié)構(gòu)網(wǎng)格示意圖

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

圖8 舵翼結(jié)構(gòu)的邊界條件示意圖

表2 D6AC鋼材料參數(shù)

2.3 模擬結(jié)果

1) 流場模擬結(jié)果

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

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

圖10 舵翼結(jié)構(gòu)前緣流場激波溫度與壓力分布曲線

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

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

2) 舵翼結(jié)構(gòu)模擬結(jié)果

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

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

圖12 舵翼結(jié)構(gòu)溫度分布和應(yīng)變分布云圖

圖13 舵翼結(jié)構(gòu)溫度與應(yīng)變沿x方向分布曲線

為了分析氣動熱和氣動力對結(jié)構(gòu)響應(yīng)的影響大小，將氣動熱和氣動力分別作用在舵翼結(jié)構(gòu)上，變形分布如圖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。可見，氣動熱作用下結(jié)構(gòu)變形值與增大幅度均遠遠大于氣動力的作用，因此，氣動熱是造成舵翼結(jié)構(gòu)破壞的主要因素。

3 舵翼結(jié)構(gòu)風洞實驗

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

圖14 氣動力和氣動熱引起的舵翼結(jié)構(gòu)變形分布云圖

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

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

(13)

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

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

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

圖16 舵翼結(jié)構(gòu)高溫應(yīng)變與溫度測量實驗裝置示意圖

圖17 舵翼結(jié)構(gòu)高溫應(yīng)變與溫度測量實驗裝置

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

圖18 應(yīng)變、溫度實驗測量與模擬曲線

4 結(jié)論

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

2) 舵翼結(jié)構(gòu)溫度場與應(yīng)變場規(guī)律相似，前緣駐點處最大，沿尾部方向逐漸降低，隨著時間推移整體溫度和應(yīng)變逐漸升高，在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。可見，氣動熱對舵翼結(jié)構(gòu)變形的影響遠大于氣動力的作用。

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