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

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

引言

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

分析高超聲速飛行器內(nèi)部的溫度場及應(yīng)力場的分布是氣動熱防護(hù)設(shè)計的關(guān)鍵，提高內(nèi)部結(jié)構(gòu)溫度場及應(yīng)力場的計算準(zhǔn)確性對于最終的飛行器的熱結(jié)構(gòu)設(shè)計十分重要[5-7]。高超聲速飛行器的氣動加熱計算是一個包括外部流場、溫度場和內(nèi)部結(jié)構(gòu)應(yīng)力場的多物理場耦合計算過程，即氣-熱-固耦合計算過程[8]。

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

1 數(shù)學(xué)模型

1.1 控制方程

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

（1）連續(xù)性方程（質(zhì)量守恒方程）：基于連續(xù)介質(zhì)假設(shè)，在笛卡爾坐標(biāo)系下，描述流動與傳熱問題的控制方程如下

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

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

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

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

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

為了保證方程的可解，必須滿足方程組的封閉，需要引入狀態(tài)方程和本構(gòu)方程。

（4）氣體狀態(tài)方程和本構(gòu)方程：一般可以用溫度、密度和壓強(qiáng)等參數(shù)來描述氣體宏觀狀態(tài)，通過熱力學(xué)狀態(tài)方程，可以得到ρ、p、T三者之間有以下關(guān)系式

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

式中，R為摩爾氣體常數(shù)。

二維應(yīng)力分量的本構(gòu)方程表示為

1.2 湍流模型

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

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

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

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

2 數(shù)值模擬參數(shù)設(shè)定

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

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

2.1 外流場建模與參數(shù)設(shè)定

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

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

計算得到一個穩(wěn)態(tài)的流場后，再以此流場的壁面處近壁面熱通量作為耦合參數(shù)導(dǎo)入 ABAQUS 6.12中，用于計算壁面熱流。讀入此時的數(shù)據(jù)文件作為耦合初始狀態(tài)，來進(jìn)行流-熱-固耦合計算。

2.2 結(jié)構(gòu)場建模與參數(shù)設(shè)定

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

在ABAQUS 6.12中進(jìn)行如下參數(shù)設(shè)定。

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

表1 30CrMnSiA合金鋼參數(shù)表Table 1 Steel parameters

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

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

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

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

（5）建立工作任務(wù)：建立工作任務(wù)時，輸出計算.inp文件，此文件記錄了ABAQUS6.12中的坐標(biāo)、約束、分析步等信息，可以自行計算或?qū)氲組PCCI4.2.1中計算，可以實現(xiàn)與FLUENT14.0的雙向數(shù)據(jù)對接。

3 模擬結(jié)果分析

3.1 外流場模擬結(jié)果分析

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

3.2 溫度場模擬結(jié)果分析

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

圖3 1 s時錐頭結(jié)構(gòu)溫度分布情況Fig.3 Temperature distribution of conehead at 1 s

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

圖4 錐頭結(jié)構(gòu)溫度隨時間變化曲線Fig.4 Curve of conehead temperature versus time

3.3 應(yīng)力場模擬結(jié)果分析

圖5 1 s時錐頭結(jié)構(gòu)變形分布情況Fig.5 Deformation distribution of conehead at 1 s

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

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

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

圖6 1 s內(nèi)不同時刻錐頭結(jié)構(gòu)應(yīng)力分布情況Fig.6 Stress distribution of conehead in 1 s

圖7 1 s內(nèi)頭錐結(jié)構(gòu)應(yīng)力分布情況Fig.7 Stress distribution of conehead in 1 s

4 結(jié) 論

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

（1）對于外流場的數(shù)值模擬，Transitionk-kl-ω湍流模型可以準(zhǔn)確計算轉(zhuǎn)捩位置和熱流數(shù)值。

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

（3）對于應(yīng)力場的數(shù)值模擬，氣動加熱1 s時，變形最大點出現(xiàn)在頭錐頂端，整體結(jié)構(gòu)的彈性變形和溫度正相關(guān)，且主要變形區(qū)域集中在前端高溫區(qū)。應(yīng)力分布較為復(fù)雜，應(yīng)力值呈現(xiàn)出由壁面向結(jié)構(gòu)內(nèi)部先降低后升高的趨勢，而錐頭內(nèi)部會出現(xiàn)一個應(yīng)力集中區(qū)域。由于升溫速率的變化，錐頂?shù)臒釕?yīng)力會先迅速上升后緩慢降低。

符號說明

aT——湍流有效擴(kuò)散系數(shù)

Cω——湍流黏性系數(shù)

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

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

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

d——壁面距離，m

F——體積力矢量，N

fW——非黏性近壁面阻尼函數(shù)

fω——邊界層項阻尼函數(shù)

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

kT——湍流動能

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

p——壓力，Pa

R——摩爾氣體常數(shù)

RBP——旁路轉(zhuǎn)捩生成項

RNAT——自然轉(zhuǎn)捩生成項

T——溫度，K

t——時間，s

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

xj——矢量位置

λ——熱導(dǎo)率，W·m-1·K-1

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

ρ——密度，kg·m-3

τ——應(yīng)力張量

ω——逆湍流時間尺度

[1] Chen Yushu(陳予恕), Guo Hulun(郭虎倫), Zhong Shun(鐘順).Review of several problems about hypersonic aircraft [J].Winged Missiles Journal(飛航導(dǎo)彈), 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(應(yīng)用數(shù)學(xué)和力學(xué)), 2008, 29(1)：47-56.

[3] Wang Xinzhi(汪新智), Ma Junjun(馬軍軍), Peng Wengen(彭穩(wěn)根),et al. Optimal design for active cooling system of hypersonic vehicle[J].Acta Aeronautica et Astronautica Sinica(航空學(xué)報), 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.