高超聲速飛行器流-熱-固耦合數值模擬

陳梅潔，程珙，田楓林，何玉榮

(哈爾濱工業大學能源科學與工程學院，黑龍江 哈爾濱 150001)

引言

高超聲速飛行器主要包括高超聲速巡航導彈、高超聲速飛機以及航天飛機等。高超聲速巡航導彈具有反應迅速、突防能力強、精確打擊和機動作戰的優勢[1-3]；高超聲速飛機可以實現全球軍事偵察、遠程軍力快速部署和遠距離攻擊；而航天飛機是未來進入太空軌道、開發太空資源最經濟、最安全的運載工具。因此，突破高超聲速關鍵技術，將對一個國家科學技術和國民經濟的發展、綜合國力的提升產生重大影響。與傳統飛行器相比，高超聲速飛行器可以大量縮短防御響應時間，增強突防和反防御能力，提高飛行器的生存能力[4]。

分析高超聲速飛行器內部的溫度場及應力場的分布是氣動熱防護設計的關鍵，提高內部結構溫度場及應力場的計算準確性對于最終的飛行器的熱結構設計十分重要[5-7]。高超聲速飛行器的氣動加熱計算是一個包括外部流場、溫度場和內部結構應力場的多物理場耦合計算過程，即氣-熱-固耦合計算過程[8]。

目前較常用的流-熱-固耦合計算方法是順序耦合計算[9]，即先將飛行器壁面考慮成溫度保持不變，通過外流場計算出壁面的熱通量分布情況，再將此熱流作為不變的邊界條件，進行結構傳熱計算。本文即采用這種方法對高超聲速飛行器流-熱-固耦合場進行模擬。對于內部溫度場和結構場，采用一種耦合單元的方法進行計算，分析高超聲速飛行器在氣動加熱環境下的瞬態熱響應。

1 數學模型

1.1 控制方程

高超聲速流動需滿足自然界三大守恒定律，也就是動量守恒定律、質量守恒定律及能量守恒定律。流體力學中把三大守恒定律描述為流體動力學所應遵循的基本方程式[10]。

（1）連續性方程（質量守恒方程）：基于連續介質假設，在笛卡爾坐標系下，描述流動與傳熱問題的控制方程如下

式中，ρ為密度，kg·m-3；t為時間，s；v為速度矢量，m·s-1。

（2）運動方程（動量守恒方程）：流體動力學中的動量守恒方程，也就是常說的N-S方程。在笛卡爾坐標系下，描述動量守恒的控制方程如下

式中，p為壓力，Pa；g為重力加速度，m·s-2，其值為9.80665 m·s-2；F為體積力矢量，N；τ為應力張量。

（3）能量方程（能量守恒方程）：當計算的流動問題包含熱量交換時，應當引入能量守恒方程。用溫度表示的能量方程為

式中，cp為比定壓熱容，J·kg-1·K-1；T為溫度，K；λ為熱導率，W·m-1·K-1。

為了保證方程的可解，必須滿足方程組的封閉，需要引入狀態方程和本構方程。

（4）氣體狀態方程和本構方程：一般可以用溫度、密度和壓強等參數來描述氣體宏觀狀態，通過熱力學狀態方程，可以得到ρ、p、T三者之間有以下關系式

三者的關系為f(p,ρ,T) =0，此即狀態方程。當氣體的溫度不是太低，氣體壓強不太高的情況下，在計算過程中，氣體分子占用的體積以及氣體各分子間的吸引力可以都不用考慮。那也就是說，對于完全氣體來說，可以將狀態方程寫為

式中，R為摩爾氣體常數。

二維應力分量的本構方程表示為

1.2 湍流模型

由于以上方程并不能使平均N-S方程不滿足方程組的封閉性，人們才引入了湍流模型來完成方程組的封閉，所以湍流模型的準確度很大程度上影響了模擬計算結果的好壞。本文的模擬將選用k-kl-ω模型。

對于高超聲速飛行時出現的邊界層轉捩現象，二方程模型并不能很好地預測轉捩的出現和轉捩位置，因此本文又采用了另一種三方程湍流模型——Transitionk-kl-ω模型，該模型能夠很好地預測邊界層發展和轉捩位置，它引入一個“分裂機制”來描述層流與湍流脈動之間的相互作用，對強激波后的溫度計算相比常用的間歇因子轉捩模型與實驗值更吻合[11-12]。

Transitionk-kl-ω模型包含湍流動能kT、層流動能kL和尺度變量ω三個輸運方程。

式中，kT為湍流動能；kL為層流動能；PkL為平均變形速度下的層流動能；RBP為旁路轉捩生成項；RNAT為自然轉捩生成項；ω為逆湍流時間尺度；DT為湍流動能在近壁面各向異性耗散項；DL為層流動能在近壁面各向異性耗散項；xj為矢量位置；v為運動黏度，m2·s-1；aT為湍流有效擴散系數；Cω為湍流黏性系數；fW為非黏性近壁面阻尼函數；fω為邊界層項阻尼函數；d為壁面距離，m；

2 數值模擬參數設定

飛行器初步建模由前面章節研究分析得知，外流場氣動加熱計算是一個非常復雜的計算過程，而流-熱-固耦合計算需要反復迭代直到得到穩態解。因此要進行流-熱-固耦合運算[13]，需針對特定部件，簡化幾何模型，來達到節省計算資源與時間的目的。本文選擇高超聲速飛行器HiFire-1受氣動加熱最為嚴重的前端錐體部分進行研究，所取得幾何模型總長為1 m，如圖1所示。

圖1 高超聲速飛行器HiFire-1幾何模型Fig. 1 Geometric model of hypersonic vehicle HiFire-1

2.1 外流場建模與參數設定

利用 ICEM 14.0軟件對簡化模型重新劃分網格，網格劃分情況如圖2所示。網格總數為8萬個，本模型的網格全部為結構化網格，并盡可能保證了單元成長方形形狀，結構化網格的計算效率高并且計算精確，這將對流場計算收斂很有幫助。對流固交界面處邊界層進行網格加密，在轉捩出現的前端錐面沿著壁面方向進行網格加密處理。離壁面最近一層的網格高度為0.001 mm，其附近網格的增長比為1.1，保證了y+值在0.2左右。

耦合計算前，外流場的模擬工況：來流 Mach數為7.16，Reynolds數為10.2×106，密度為0.072 kg·m-3，壓力為4.62 kPa，溫度為231.7 K。選擇耦合隱式求解器，雙精度，定常求解，單元體中心處結果變量的梯度選擇Green-Gauss Node-Based來提高計算精度。湍流模型選擇Transitionk-kl-ω模型，外邊界選擇壓力遠場邊界條件與壓力出口邊界條件。壁面設置一個初始溫度，此處為300 K，在這個壁溫下，可以得到一個穩態的流場情況，進而得到壁面的熱通量。本處計算時用的是耦合隱式求解器。

計算得到一個穩態的流場后，再以此流場的壁面處近壁面熱通量作為耦合參數導入 ABAQUS 6.12中，用于計算壁面熱流。讀入此時的數據文件作為耦合初始狀態，來進行流-熱-固耦合計算。

2.2 結構場建模與參數設定

將幾何模型導入 ABAQUS 6.12軟件進行網格化分，由于錐頭部位尺寸較小，而熱通量較大，需要進行網格加密，以保證計算精度。錐形整體采用漸變網格，單元類型選擇軸對稱耦合單元CAX4RT，以便實現傳熱與結構變形的直接耦合計算。

在ABAQUS 6.12中進行如下參數設定。

圖2 外流場網格化分示意圖Fig. 2 Schematic diagram of external flow field grid

表1 30CrMnSiA合金鋼參數表Table 1 Steel parameters

（1）材料屬性：為了初步研究飛行器內部溫度場與應力場分布情況，首先選用工程中常用的30CrMnSiA合金鋼材料作為機體材料，具體參數見表1。

（2）分析步設定：建立分析步，選擇瞬態計算類型，分析類型為直接耦合，計算時間步長取0.0001 s，總時間為1 s。

（3）建立約束：約束對稱軸的縱向位移和錐尾處的橫向位移。

（4）定義耦合單元集：定義飛行器外壁面為耦合單元集，固體壁面形狀坐標要與流場中的壁面坐標一致，以便形成耦合區域進行變量交換傳遞。

（5）建立工作任務：建立工作任務時，輸出計算.inp文件，此文件記錄了ABAQUS6.12中的坐標、約束、分析步等信息，可以自行計算或導入到MPCCI4.2.1中計算，可以實現與FLUENT14.0的雙向數據對接。

3 模擬結果分析

3.1 外流場模擬結果分析

由于本文將壁面處的熱通量作為耦合變量，所以在對外流場的模擬結果中，只需考察其壁面熱通量與實驗的對比情況即可，對比后發現使用Transitionk-kl-ω湍流模型可以準確計算轉捩位置和熱流數值。

3.2 溫度場模擬結果分析

由熱通量分布情況得知，熱通量最大部位在錐頭，此處的氣動加熱情況最為嚴重。取0.03 m錐頭對稱截面進行分析，圖3為1 s時頭錐的截面溫度分布情況，由圖可以看出，氣動加熱1 s時，溫度最高點出現在頭錐頂端，為1000 K，在軸線方向上呈現環狀的分布，沿軸線逐漸降低，且變化梯度逐漸降低。這是因為，錐頭部位為高熱流區，開始加熱時，流場與壁面溫差較大，傳熱非常劇烈，溫升較快，進而形成高溫區，由于錐體內傳熱系數的限制，熱量不能及時被傳導至尾端，導致錐頭形成較大的熱梯度。

圖3 1 s時錐頭結構溫度分布情況Fig.3 Temperature distribution of conehead at 1 s

取錐頭頂端節點（錐頂）為研究對象。圖4為1 s內，頂端節點的溫度隨時間的變化情況，由圖可以看出，由于氣動加熱的原因，頂端節點的溫度會隨時間推移逐漸上升，且變化趨勢逐漸減緩。這是因為，錐頭區域溫度隨著時間的推移逐漸升高，30CrMnSiA合金鋼材料的熱容會逐漸升高，而導致頂端節點處壁面溫度升高的速率逐漸降低。

圖4 錐頭結構溫度隨時間變化曲線Fig.4 Curve of conehead temperature versus time

3.3 應力場模擬結果分析

圖5 1 s時錐頭結構變形分布情況Fig.5 Deformation distribution of conehead at 1 s

圖5為1 s時頭錐截面彈性變形分布情況，由圖可以看出，變形最大點出現在頭錐頂端，為5×10-5，且在軸線方向上呈現環狀的分布，沿軸線逐漸降低，且變化梯度逐漸降低。整體結構的彈性變形和溫度正相關，且主要變形區域集中在前端高溫區。

圖6為1 s內不同時刻頭錐的截面應力分布情況，由圖可以看出，在0.1 s時，錐頂處為應力最大值點，從ABAQUS中讀取數值為1.8×108Pa。但隨著事件的推移，應力最大值點會沿著壁面向后移動，且應力值會逐漸減小。最終1 s時應力最大值出現在距離頂部 4 mm左右的壁面處，數值為1.15×108Pa，應力分布較為復雜，由壁面向內部應力值先降低后升高，而錐頭內部會出現一個應力集中區域。

以錐頭頂端節點為研究對象。圖7為頂端節點的應力隨時間的變化情況，由圖7可以看出，由于氣動加熱的原因，頂端節點的應力會迅速上升，在0.008 s左右達到最大值1.82×108Pa。之后隨時間推移，應力值會逐漸下降，達到1 s時，應力值降為 11.0×107Pa。這是因為剛開始加熱時，錐頂節點升溫速率極快，錐頂和附近各節點的溫度及變形都存在很大的梯度，導致錐頂的熱應力迅速上升。之后因為錐身溫度不斷升高，錐頭部位溫度及變形梯度降低，而使得錐頂節點的熱應力值降低。

圖6 1 s內不同時刻錐頭結構應力分布情況Fig.6 Stress distribution of conehead in 1 s

圖7 1 s內頭錐結構應力分布情況Fig.7 Stress distribution of conehead in 1 s

4 結 論

本文選用了HiFire-1錐體部分作為研究對象，采用順序耦合方法，對飛行器錐體的流-熱-固耦合場進行了模擬分析。通過分析模擬結果，得出如下結論。

（1）對于外流場的數值模擬，Transitionk-kl-ω湍流模型可以準確計算轉捩位置和熱流數值。

（2）對于溫度場的數值模擬，氣動加熱1 s時，溫度最高點出現在頭錐頂端，在軸線方向上呈現環狀的分布，沿軸線逐漸降低，且變化梯度逐漸降低。且頂端節點的溫度會隨時間推移逐漸上升，且變化趨勢逐漸減緩。

（3）對于應力場的數值模擬，氣動加熱1 s時，變形最大點出現在頭錐頂端，整體結構的彈性變形和溫度正相關，且主要變形區域集中在前端高溫區。應力分布較為復雜，應力值呈現出由壁面向結構內部先降低后升高的趨勢，而錐頭內部會出現一個應力集中區域。由于升溫速率的變化，錐頂的熱應力會先迅速上升后緩慢降低。

符號說明

aT——湍流有效擴散系數

Cω——湍流黏性系數

cp——比定壓熱容，J·kg-1·K-1

DL——層流動能在近壁面各向異性耗散項

DT——湍流動能在近壁面各向異性耗散項

d——壁面距離，m

F——體積力矢量，N

fW——非黏性近壁面阻尼函數

fω——邊界層項阻尼函數

g——重力加速度，m·s-2，其值為9.80665 m·s-2

kT——湍流動能

PkL——平均變形速度下的層流動能

p——壓力，Pa

R——摩爾氣體常數

RBP——旁路轉捩生成項

RNAT——自然轉捩生成項

T——溫度，K

t——時間，s

v——速度矢量，m·s-1

xj——矢量位置

λ——熱導率，W·m-1·K-1

ν——運動黏度，m2·s-1

ρ——密度，kg·m-3

τ——應力張量

ω——逆湍流時間尺度

[1] Chen Yushu(陳予恕), Guo Hulun(郭虎倫), Zhong Shun(鐘順).Review of several problems about hypersonic aircraft [J].Winged Missiles Journal(飛航導彈), 2009, (8)：26-33.

[2] Yang Yazheng(楊亞政), Yang Jialing(楊嘉陵), Fang Daining(方岱寧). Research progress of hypersonic flight vehicle thermal protection materials and the structure [J].Applied Mathematics and Mechanics(應用數學和力學), 2008, 29(1)：47-56.

[3] Wang Xinzhi(汪新智), Ma Junjun(馬軍軍), Peng Wengen(彭穩根),et al. Optimal design for active cooling system of hypersonic vehicle[J].Acta Aeronautica et Astronautica Sinica(航空學報), 2014, 35(3)：624-633.

[4] Yentsch R J, Gaitonde D V. Numerical investigation of hypersonic phenomena encountered in hifire flight 1//50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition [C]. Nashville, TN：AIAA, 2012：1-12.

[5] Yentsch R J, Gaitonde D V, Kimmel R. Performance of turbulence modeling in simulation of the HIFiRE-1 flight test [J].J. Spacecr.Rockets, 2013, 51(1)：117-127.

[6] Henline W D, Palmer G E, Chen Y K,et al. Aerothermodynamic heating analysis and heatshield design of an SSTO rocket vehicle for Access-to-Space//AIAA, Thermophysics Conference[C]. San Diego,CA, 1995.

[7] Olynck D R, Henline W D. Navier-Stokes heating calculations for benchmark thermal protection system sizing [J].J. Spacecr. Rockets,1996, 33(6)：807-814.

[8] Chen Y K, Milos F S. Solution strategy for thermal response of nonblating thermal protection systems at hypersonic speeds [J].AIAA Paper, 1996：96-0615.

[9] Olynic D, Tam T. Trajectory-based validation of the shuttle heating environment [J].J. Spacecr. Rockets, 1997, 34(2)：172-181.

[10] Molis F S, Squire T H. Thermostructural analysis of SIRCA tile for X-34 wing leading edge TPS[J].AIAA Paper, 1998：98-0883.

[11] Hollis B R, Horvath T J, Berry S A,et al. X-33 computational aeroheating predictions and comparisons with experimental data [J].J.Spacecr. Rockets, 2001, 38(5)：658-669.

[12] Hassan B, Kuntz D W, Salguero D E,et al. A coupled fluid/thermal/flight dynamics approach for predicting hypersonic vehicle performance [J].AIAA Paper, 2001：2001-2903.

[13] Holden M S, Wadhams T P, Smolinski G J,et al. Experimental and numerical studies on hypersonic vehicle performance in the LENS shock and expansion tunnels [J].AIAA Paper, 2006, 125：2006.

[14] Walters D K, Cokljat D. A there-equation eddy-viscosity model for Reynolds-averaged Navier-Stokes simulations of transitional flows [J].J. Fluids Eng.-Trans. ASME, 2008, 130(1)：1-14.