超音速導彈彈翼結構的氣動熱彈性分析

劉立剛，周 凌，孫 輝

(中國科學院長春光學精密機械與物理研究所，長春 130033)

隨著有翼導彈的飛行速度越來越高，尤其是在超音速下進行機動飛行的有翼導彈，氣動加熱引起的高溫將導致彈翼結構部件內產生不均勻的熱應力和熱變形等現象，繼而導致飛行器氣動外形、結構強度及結構剛度的變化，影響導彈的飛行性能和戰斗性能，有時甚至能導致結構破壞［1，2］，這在工程實際中是經常遇到的，所以超音速導彈彈翼這類飛行器的典型結構的氣動熱彈性分析是一個非常重要的問題。

李建林等給出了復雜外形高超聲速飛行器氣動熱的快速工程估算方法［3］。楊愷等發展了一套高超聲速飛行器關鍵部位氣動熱的計算方法［4，5］。趙曉利對圓柱殼高超聲速氣動加熱風洞實驗進行了三維非定常數值模擬，實驗結果驗證了耦合計算方法的準確性［6］。李國曙等建立了考慮熱效應影響的高超聲速飛行器從氣動熱計算到靜氣動彈性分析的快速分析方法［7］。楊超等建立了氣動熱、氣動彈性雙向耦合高超聲速二維曲面壁板顫振分析方法［8］。

本文給出了低、中超音速導彈彈翼結構的氣動熱彈性分析流程，并給出Ma=5時氣動加熱引起的彈翼結構的溫度場、熱應力場、熱應變場和熱位移場的仿真結果。

1超音速彈翼表面氣動加熱工程計算

在計算彈翼表面氣動加熱分布時，將其分為高速可壓層流區、高速可壓紊流區與彈翼前緣3個區域分別計算，由式(1)先確定彈翼蒙皮變量的層流與紊流區的轉捩點坐標xtri

式(1)中:ρ為空氣密度;u為氣流速度;μ為空氣動力黏度;Retri為轉捩點處的雷諾數，下標帶e為邊界層外緣值。對于光滑壁，工程上考慮馬赫數影響的轉捩準則計算，Retri表達式為

式(2)中，Ma為馬赫數。

1.1 層流區、紊流區氣動加熱計算

將彈翼結構簡化成平板，由平板的熱流密度公式和Eckert參考溫度法，可以得到高速可壓層流的平板氣動加熱公式為［9］

式(3)中:q為熱流密度;Pr為普朗特數，這里取為0.7;ρ為空氣密度;u為氣流速度;μ為空氣動力黏度;Rex為當地雷諾數為空氣定壓比熱容;Taw為恢復溫度;Tw為彈翼表面壁溫;下標帶e為邊界層外緣值，式中 μ*，ρ*，的工程計算方法如下:

R=287.1(J˙K－1˙kg－1)為單位氣體常數，w=0.76。

同樣可得高速可壓紊流的平板加熱公式為［10］

在高溫邊界層中，傳熱的基本形式有4種:導熱、對流傳熱、擴散傳熱和輻射傳熱。這里只考慮邊界層對流傳熱和彈翼表面輻射傳熱兩種方式。彈翼表面向外輻射的熱流量q1表達式為

式(6)中:ε 稱為表面黑度因數;σ =5.67 ×10－8W˙m－2˙ K－4，是斯蒂芬—波爾茲曼常數。

高溫邊界層對彈翼表面加熱的熱流密度為q，在考慮定常、無內壁冷卻情況的熱流量平衡方程為

在層流狀態下，q為式(3)，將其代入式(7)，可以得到彈翼表面壁溫公式

在紊流狀態下，q為式(5)，代入到式(7)，同樣可以得到彈翼表面壁溫公式

1.2 彈翼前緣氣動加熱計算

若不考慮彈體對彈翼的影響，彈翼前緣的熱流密度可近似按后掠圓柱前緣的熱流公式進行計算，并代入到熱流量平衡方程式(7)得

式(10)中，α為熱交換系數

式(11)中，Λ為彈翼后掠角，下標sl為前緣線上的值。θ'w的數值計算關聯式為

式(12)中，T∞o和Tno分別為來流總溫和垂直于翼前緣方向的動能轉化為熱能相應的總溫

式(15)中，Man∞為來流馬赫數在前緣的法向分量;ρno為氣體密度;R0為翼前緣半徑。PwslΛ為相應于Man∞的正激波后總壓，表達式如下:

2 熱彈性分析的有限元法

在求得超音速導彈彈翼表面的溫度分布之后，將其作為彈翼結構溫度場的邊界條件，采用有限元法計算彈翼結構的溫度場以及熱應力場、熱應變場和熱位移場。

2.1 穩態溫度場有限元法

將求解域分為有限個單元，設單元結點為 i，j，m，…，p，結點溫度為 Ti，Tj，Tm，…，Tp，單元內任一點溫度 Te(x，y，z)用結點溫度表示為

式(19)中，［N］為形函數矩陣，是坐標 x，y，z的函數，{T}e為結點溫度矩陣。根據變分原理，對于第一類邊界條件，泛函I(T)取得極小值，即對于求解域的全部結點，得到

式(20)中，{T}=［T1，T2，T3，…，Tn］T為包含全部結點的結點溫度矩陣，［H］為熱傳導矩陣。另在第一類邊界上的結點溫度為已知邊界溫度，在此條件下求解方程組式(20)，即可得到全部結點溫度，各單元內部任一點溫度可由式(19)插值求得。

2.2 熱應力場、熱應變場和熱位移場有限元法

引起熱應力的根本原因是在約束作用下且有溫度的變化。用位移法分析結構熱應力，應先按溫度場的改變計算熱變形，進而計算熱應力。一個三維彈翼受熱有溫度變化會發生形狀變化，空間內各點都有一定位移，溫度變化引起的位移求解方程:

式(21)中，［K］為結構總剛度矩陣;RΔT為全部單元熱載荷的迭加，是結構由于受熱膨脹而引起的等效為結點的載荷列陣，{δ}為總結點位移。熱應變可由下式得到［11］:

進而得到熱應力的表達式為

式(23)中，［D］為空間問題的彈性系數矩陣。

3 氣動熱彈性分析的流程

將氣動加熱工程計算與有限元法相結合，圖1給出彈翼氣動熱彈性計算流程。

圖1 彈翼氣動熱彈性計算流程

4 氣動熱彈性分析工程算例

彈翼翼面結構為夾層結構翼面，其外形與有限元網格模型如圖2所示。彈翼表面采用蜂窩夾芯蒙皮，內部由3根翼肋與前墻、后墻組成，與彈身的連接是由翼根處前墻和后墻截面上的接頭進行連接［12］。其飛行狀態和大氣物理參數如表1所示。

圖2 彈翼結構及有限元網格

彈翼主要由兩種材料制造而成，表面蒙皮采用的是蜂窩夾層板材料，內部骨架結構的材料牌號為HasteloyX，是一種輕質合金。它們的物理參數和熱性能參數分別見表2和表3。

表1 飛行狀態和大氣物理參數

表2 材料物理參數

表3 彈性模量和熱參數與溫度的關系

4.1 彈翼表面氣動加熱計算

在進行氣動加熱估算時，彈翼前緣為x坐標原點。當Ma=3時，由式(1)、式(2)計算得到轉捩點xtri=0.04 m。

4.1.1 彈翼前緣氣動加熱計算

由式(7)及式(10)～式(18)可得方程

對式(24)進行求解，可得Tw=889.576 K。

4.1.2 層流彈翼表面壁溫計算

層流區為0 ＜x≤0.04，由式(8)、式(14)可得方程:

計算得到Tw沿x的變化如圖3所示，其中0＜x≤0.04。

圖3 Ma=3時層流Tw變化曲線

4.1.3 紊流彈翼表面壁溫計算

紊流區為 0.04 ＜x≤0.605，由式(9)、式(14)可得方程:

計算得到Tw沿x的變化如圖4所示，其中0.04＜x≤0.605。

圖4 Ma=3時紊流Tw變化曲線

4.2 彈翼結構的溫度場

在仿真計算彈翼結構溫度場時，將4.1節由氣動加熱工程估算公式計算得到的彈翼表面蒙皮的氣動熱分布作為外壁邊界條件，彈翼結構內壁的邊界條件是空氣對流邊界條件，空氣對流換熱系數為 12.5 w˙m－2˙°C－1，空氣的初始溫度為25℃。圖5是計算得到的彈翼結構溫度場分布云圖。圖6～圖10給出了蒙皮沿弦向、厚度方向、翼肋沿弦向、前墻與后墻沿展向溫度分布。

圖5 彈翼結構及內部骨架溫度分布

圖6 蒙皮弦向溫度分布

圖7 蒙皮厚度方向溫度分布

圖8 翼肋弦向溫度分布

圖9 前墻展向溫度分布

圖10 后墻展向溫度分布

從以上的溫度場分布圖可以看出:彈翼結構最高溫度為289.108℃，是蒙皮表面所加載的最高溫度，出現在層流與紊流的轉捩點處。沿蒙皮厚度方向，溫度呈梯度分布，從外壁到內壁遞減變化。沿弦長方向，翼肋在靠近前緣、后緣附近的部分溫度較高，中間的翼肋溫度較低;前墻、后墻在翼根部溫度較低，翼梢處溫度較高。

4.3 彈翼結構的熱彈性分析

將溫度場的結點溫度當作溫度載荷施加在應力場分析中。并且對于在翼根與彈身連接處的約束作如下處理，將約束加在翼根處前墻和后墻截面上與上下蒙皮截面連接的結點上，圖11～圖22分別是熱應力場、熱應變場和熱位移場分布云圖與變化趨勢圖。

圖11 彈翼結構熱應力(Mises應力)及內部骨架結構Mises應力云圖

圖12 蒙皮沿弦向熱應力(Mises應力)變化曲線

圖13 翼肋沿弦向熱應力(Mises應力)變化曲線

圖14 前墻沿展向熱應力(Mises應力)變化曲線

圖15 后墻沿展向熱應力(Mises應力)變化曲線

圖16 彈翼結構熱應變(Mises應變)及內部骨架結構Mises應變云圖

圖17 蒙皮沿弦向熱應變(Mises應變)變化曲線

圖18 翼肋沿弦向熱應變(Mises應變)變化曲線

圖19 前墻沿展向熱應變(Mises應變)變化曲線

圖20 后墻沿展向熱應變(Mises應變)變化曲線

圖21 彈翼結構熱位移及內部骨架結構合位移云圖

圖22 蒙皮沿弦向熱位移變化曲線

從彈翼結構的熱應力場、熱應變場和熱位移場分布與趨勢圖可以得出以下結論:熱應力最大值為7.93×108Pa，出現在彈翼與彈身的連接處;沿弦長方向，蒙皮表面前緣、后緣的熱應力較小，靠近連接處的熱應力較大。由于溫度載荷的對稱性，在蒙皮表面的應力場也是上下對稱的。翼肋沿弦向熱應力中間較大，兩端較小。前墻熱應力沿展向翼梢較大，后墻沿展向熱應力翼根與翼梢較大，中間較小。熱應變與熱應力分布基本一致，最大值為4.367×10－3。最大熱位移出現在翼梢的尾部，其值為8.52×10－4m;從計算結果來看，最大熱應力、最大熱應變、最大熱位移均在許可值范圍內，可見該彈翼的材料及結構滿足設計要求。

5 結論

從算例分析結果得到以下結論:最高溫度出現在層流與紊流的轉捩點處，其值為289.108℃;最大熱應力和熱應變均出現在彈翼與彈身的連接處。最大熱應力值為7.93×108Pa，最大熱應變值為4.367×10－3，可見熱應力數值是非常大的，所以對于高超聲速飛行的導彈彈翼進行熱應力校核是完全必要的;最大熱位移出現在翼梢的尾部，其值為8.52×10－4m;最大熱應力、最大熱應變、最大熱位移均在許可值范圍內，可見該彈翼的材料及結構滿足設計要求。

［1］ 姜貴慶，劉連元.高速氣流傳熱與燒蝕熱防護［M］.北京:國防工業出版社，2003.

［2］ 劉芹.彈體結構熱振動分析［D］.西安:西北工業大學，2004.

［3］ 李建林，唐乾剛，霍霖，等.復雜外形高超聲速飛行器氣動熱快速工程估算［J］.國防科技大學學報，2012，34(6):89－93.

［4］ 楊愷，高效偉.高超聲速飛行器關鍵部位氣動熱計算［J］.計算力學學報，2012，29(1):13 －18.

［5］ 趙曉利，孫振旭，安亦然，等.高超聲速氣動熱的耦合計算方法研究［J］.科學技術與工程，2010，10(22):5450－5455.

［6］ 李國曙，萬志強，楊 超.高超聲速翼面氣動熱與靜氣動彈性綜合分析［J］.北京航空航天大學學報，2012，38(1):53－58.

［7］ 楊超，李國曙，萬志強.氣動熱－氣動彈性雙向耦合的高超聲速曲面壁板顫振分析方法［J］.中國科學:技術科學，2012，42(4):369 －377.

［8］ 楊 愷，高效偉.高超聲速氣動熱環境工程算法［J］.導彈與航天運載技術，2010(4):19－23.

［9］ 錢翼稷.空氣動力學［M］.北京:北京航空航天大學出版社，2004.

［10］中國人民解放軍總裝備部軍事訓練教材編輯工作委員會.高超聲速氣動熱和熱防護［M］.北京:國防工業出版社，2003.

［11］朱加銘，歐貴寶，何蘊增.有限元與邊界元法［M］.哈爾濱:哈爾濱工程大學出版社，2002.

［12］酈正能.飛行器結構學［M］.北京:北京航空航天大學出版社，2005.

(責任編輯周江川)