劉立剛，周 凌，孫 輝

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

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

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

本文給出了低、中超音速導(dǎo)彈彈翼結(jié)構(gòu)的氣動熱彈性分析流程，并給出Ma=5時氣動加熱引起的彈翼結(jié)構(gòu)的溫度場、熱應(yīng)力場、熱應(yīng)變場和熱位移場的仿真結(jié)果。

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

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

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

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

1.1 層流區(qū)、紊流區(qū)氣動加熱計算

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

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

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

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

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

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

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

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

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

1.2 彈翼前緣氣動加熱計算

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

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

式(11)中，Λ為彈翼后掠角，下標(biāo)sl為前緣線上的值。θ'w的數(shù)值計算關(guān)聯(lián)式為

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

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

2 熱彈性分析的有限元法

在求得超音速導(dǎo)彈彈翼表面的溫度分布之后，將其作為彈翼結(jié)構(gòu)溫度場的邊界條件，采用有限元法計算彈翼結(jié)構(gòu)的溫度場以及熱應(yīng)力場、熱應(yīng)變場和熱位移場。

2.1 穩(wěn)態(tài)溫度場有限元法

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

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

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

2.2 熱應(yīng)力場、熱應(yīng)變場和熱位移場有限元法

引起熱應(yīng)力的根本原因是在約束作用下且有溫度的變化。用位移法分析結(jié)構(gòu)熱應(yīng)力，應(yīng)先按溫度場的改變計算熱變形，進(jìn)而計算熱應(yīng)力。一個三維彈翼受熱有溫度變化會發(fā)生形狀變化，空間內(nèi)各點都有一定位移，溫度變化引起的位移求解方程:

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

進(jìn)而得到熱應(yīng)力的表達(dá)式為

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

3 氣動熱彈性分析的流程

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

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

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

彈翼翼面結(jié)構(gòu)為夾層結(jié)構(gòu)翼面，其外形與有限元網(wǎng)格模型如圖2所示。彈翼表面采用蜂窩夾芯蒙皮，內(nèi)部由3根翼肋與前墻、后墻組成，與彈身的連接是由翼根處前墻和后墻截面上的接頭進(jìn)行連接［12］。其飛行狀態(tài)和大氣物理參數(shù)如表1所示。

圖2 彈翼結(jié)構(gòu)及有限元網(wǎng)格

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

表1 飛行狀態(tài)和大氣物理參數(shù)

表2 材料物理參數(shù)

表3 彈性模量和熱參數(shù)與溫度的關(guān)系

4.1 彈翼表面氣動加熱計算

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

4.1.1 彈翼前緣氣動加熱計算

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

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

4.1.2 層流彈翼表面壁溫計算

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

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

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

4.1.3 紊流彈翼表面壁溫計算

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

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

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

4.2 彈翼結(jié)構(gòu)的溫度場

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

圖5 彈翼結(jié)構(gòu)及內(nèi)部骨架溫度分布

圖6 蒙皮弦向溫度分布

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

圖8 翼肋弦向溫度分布

圖9 前墻展向溫度分布

圖10 后墻展向溫度分布

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

4.3 彈翼結(jié)構(gòu)的熱彈性分析

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

圖11 彈翼結(jié)構(gòu)熱應(yīng)力(Mises應(yīng)力)及內(nèi)部骨架結(jié)構(gòu)Mises應(yīng)力云圖

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

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

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

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

圖16 彈翼結(jié)構(gòu)熱應(yīng)變(Mises應(yīng)變)及內(nèi)部骨架結(jié)構(gòu)Mises應(yīng)變云圖

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

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

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

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

圖21 彈翼結(jié)構(gòu)熱位移及內(nèi)部骨架結(jié)構(gòu)合位移云圖

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

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

5 結(jié)論

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

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

［2］ 劉芹.彈體結(jié)構(gòu)熱振動分析［D］.西安:西北工業(yè)大學(xué)，2004.

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

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

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

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

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

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

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

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

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

［12］酈正能.飛行器結(jié)構(gòu)學(xué)［M］.北京:北京航空航天大學(xué)出版社，2005.

(責(zé)任編輯周江川)