考慮防熱層以及燒蝕的高超聲速飛行器動特性計算模型研究

王亮 張妍 蔡毅鵬 蔡文杰

摘 要：針對全大氣層滑翔的高超聲速飛行器動特性預(yù)示需求， 本文提出了考慮防熱層以及燒蝕的高超聲速飛行器截面剛度計算模型。 首先， 提出了三層組合的高超聲速飛行器截面剛度計算模型， 分別為最內(nèi)層承力金屬殼體、 中間防熱層原始層和最外防熱層燒蝕層。 其次， 根據(jù)剛度線性疊加原理， 給出了各層以及總的截面剛度參數(shù)的計算方法。 最后， 針對實際使用， 給出了高超聲速飛行器結(jié)構(gòu)動特性具體的計算分析程序和步驟。

關(guān)鍵詞：???? 高超聲速;? 動特性;? 熱防護(hù)系統(tǒng);? 燒蝕效應(yīng);? 抗彎剛度

中圖分類號：???? TJ760; V212.11 文獻(xiàn)標(biāo)識碼：??? A? 文章編號：???? 1673-5048（2021）06-0072-04

0 引? 言

高超聲速飛行器是結(jié)構(gòu)組成十分復(fù)雜的彈性體， 其動力學(xué)穩(wěn)定性設(shè)計是飛行器穩(wěn)定飛行的前提， 而結(jié)構(gòu)動力學(xué)特性是飛行器姿態(tài)穩(wěn)定控制系統(tǒng)設(shè)計的關(guān)鍵參數(shù)之一， 參數(shù)的準(zhǔn)確性將直接影響飛行成敗。 結(jié)構(gòu)動力學(xué)特性預(yù)示是理論性和實踐性較強的一項工作， 因此其重要性不言而喻。

在分析細(xì)長體的高超聲速飛行器時， 常采用梁單元進(jìn)行模態(tài)分析。 學(xué)者們在梁單元精確分析方面進(jìn)行了大量的研究[1-9]， 提出了較多梁單元剛度精細(xì)化建模方法。 如多種截面剪切系數(shù)計算方法， 對梁單元中的剪切效應(yīng)進(jìn)行有效修正;? 包含自重項的細(xì)長梁單元剛度矩陣， 其中考慮了軸力對彎曲剛度的影響及彎矩引起的軸向剛度修正系數(shù);? 針對管狀復(fù)合材料梁提出一種復(fù)合梁單元的有限元建模新方法。

典型的高超聲速飛行器， 如傳統(tǒng)的戰(zhàn)術(shù)彈道導(dǎo)彈， 其為再入式飛行， 雖然再入速度較高， 遭受嚴(yán)酷的氣動加熱影響， 但由于再入飛行時間較短， 且再入體外部防熱層較厚， 因此氣動加熱效應(yīng)不明顯， 殼體溫升和防熱層燒蝕嚴(yán)酷度不高。 再入體的結(jié)構(gòu)動特性預(yù)示難度較低， 采用常溫材料參數(shù)進(jìn)行動特性預(yù)示精度較高， 偏差帶可以包絡(luò)由于短時防熱層燒蝕和殼體溫升帶來的動特性變化。 在此前提下， 學(xué)者們在研究高超聲速飛行器動特性預(yù)示時， 一般采用考慮常溫下飛行器的模態(tài)特性預(yù)示方法， 或考慮高溫對承力殼體的影響下的模態(tài)預(yù)示方法。 何成[10]基于結(jié)構(gòu)動力學(xué)理論、 熱傳導(dǎo)理論、 模型修正與參數(shù)辨識方法， 以理論分析與仿真試驗驗證相結(jié)合的方式開展熱環(huán)境下的高超聲速飛行器結(jié)構(gòu)動力學(xué)的建模方法研究。 王亮等[11-13]以高超聲速戰(zhàn)術(shù)導(dǎo)彈為應(yīng)用背景， 研究了各類全彈的梁-質(zhì)量塊建模方法， 給出了常溫下多種導(dǎo)彈全彈有限元質(zhì)量模型、 扭轉(zhuǎn)模態(tài)模型和橫向模態(tài)模型， 對比了各類模型對計算結(jié)果的影響， 并與試驗結(jié)果進(jìn)行了對比， 給出了全彈有限元建模時的推薦模型。 張琪等[14]針對以往方法0維點單元和3維實體單元等效非結(jié)構(gòu)質(zhì)量中存在的問題， 提出等效密度法， 即根據(jù)密度和剛度等效的原則， 將非結(jié)構(gòu)質(zhì)量等效均布到殼體單元， 修改殼體材料密度以保證其質(zhì)量不變。 應(yīng)用3種方法分別建立了導(dǎo)彈有限元模型， 并進(jìn)行了導(dǎo)彈模態(tài)試驗， 將3種方法得到的模態(tài)頻率與試驗值比較， 發(fā)現(xiàn)等效密度法建立的模型更為合理， 應(yīng)用簡便。

綜上所述， 學(xué)者們在研究梁單元精細(xì)建模和高超聲速飛行器常規(guī)建模方面做出了很多工作。 但是考慮高超聲速飛行工況下， 氣動加熱下全飛行器動特性模型建模的研究偏少， 尤其是考慮防熱層燒蝕情況下的動特性模型尚無公開研究。

隨著全大氣層高超滑翔技術(shù)的廣泛應(yīng)用， 與再入飛行模式相比， 氣動加熱同樣惡劣， 而且一方面氣動加熱時間較長， 另一方面為了追求方案最優(yōu)， 材料的承載能力也進(jìn)一步挖潛， 使用溫度的上限也逐步提高。 由此出現(xiàn)的熱模態(tài)問題顯著， 對飛行器結(jié)構(gòu)動特性預(yù)示帶來了難題， 尤其是在對飛行過程中的結(jié)構(gòu)動特性變化趨勢的精確預(yù)示， 涉及到殼體熱模態(tài)和考慮部分燒蝕的外部防熱層的熱模態(tài)預(yù)示， 其中關(guān)鍵是缺少飛行過程中的熱模態(tài)預(yù)示的模型， 這一模型的研究顯得意義重大。 因此， 本文針對全大氣層滑翔的高超聲速飛行器動特性預(yù)示需求， 提出了考慮防熱層以及燒蝕的高超聲速飛行器截面剛度計算模型， 并進(jìn)行了詳細(xì)的推導(dǎo)和分析。

1 模態(tài)計算理論

王 亮， 等： 考慮防熱層以及燒蝕的高超聲速飛行器動特性計算模型研究

飛行器結(jié)構(gòu)動特性求解方程可以表示為

KΦ=ΛMΦ（1）

式中： K， M為飛行器結(jié)構(gòu)的總剛度、 總質(zhì)量矩陣;? Λ， Φ為飛行器結(jié)構(gòu)的特征值與特征向量矩陣。

因此， 在計算飛行器結(jié)構(gòu)動特性時， 首先要確定飛行器結(jié)構(gòu)動力學(xué)模型的總剛度與總質(zhì)量矩陣。 而結(jié)構(gòu)動力學(xué)模型的總剛度矩陣由選取的梁模型特性所決定， 因此截面剛度模型對梁模型剛度特性影響是顯而易見的。

2 動特性模型

高超聲速飛行器結(jié)構(gòu)動力學(xué)建模時， 一般將飛行器結(jié)構(gòu)簡化為梁-質(zhì)量塊模型， 并將飛行器結(jié)構(gòu)分為若干個分站（這些分站一般為艙段連接點、 儀器放置點， 或者某些大型部組件的重心）， 給出高超聲速飛行器各分站處的集中質(zhì)量以及各分站之間梁的材料和截面尺寸。

2.1 質(zhì)量模型

每個節(jié)點為一分站， 分站質(zhì)量與轉(zhuǎn)動慣量添加在對應(yīng)的自由度上， 對應(yīng)分站的轉(zhuǎn)動慣量計算公式為

J=me2（R2+46h2），對于端點me2（R2+16h2），對于中間點 （2）

式中： me為該站的集中質(zhì)量;? R為該站的彈體半徑;? h處于端點時， 取與相鄰點間距之半， 處于中間點時， 取相鄰前后兩點間距之半。

2.2 剛度模型

高超聲速飛行器總剛度矩陣由各梁單元剛度矩陣組裝得到。 一方面由于高超聲速飛行器為連續(xù)氣動外形， 各節(jié)點處截面半徑不相等， 需采用變截面梁單元技術(shù);? 另一方面， 由于高超聲速飛行器長時間氣動加熱效應(yīng)的影響， 內(nèi)部金屬殼體溫度隨時間逐漸升高， 外部防熱層沿周向和厚度方向溫度存在差異， 而且迎風(fēng)面會由于氣動加熱產(chǎn)生燒蝕和碳化， 厚度較原始狀態(tài)會發(fā)生改變。 因此， 計算梁單元剛度矩陣時， 使用型函數(shù)導(dǎo)數(shù)與材料參數(shù)矩陣乘積在單元長度上積分的方法， 得到各單元的剛度矩陣。 梁單元的剛度計算公式為

[Ke]=∫l0（EIBTB+kGAε2NTN）dx（3）

式中： G為截面剪切模量;? A為截面面積;? k為剪應(yīng)力修正因子， 這里取4/3[3];? ε=6EIl/12EIl+GAlk;? [B]=[B]-6εl2xl-1[N]， B為梁單元的型函數(shù)對坐標(biāo)的二次導(dǎo)數(shù);? [N]=[2/l1-2/l1]T， l為梁單元長度;? EI（x）為梁單元的抗彎剛度， 有

N=1-3x2/l2+2x3/l3x-2x2/l+x3/l23x2/l2-2x3/l3-x2/l+x3/l2;

B=dN2dx2=-6/l2+12x/l3-4/l+6x/l26/l2-12x/l3-2/l+6x/l2。

針對不同截面的抗彎剛度， 本文提出了考慮帶燒蝕的防熱層的飛行器截面抗彎剛度計算模型。

首先以旋成體高超聲速飛行器為例， 給出了飛行器典型截面的示意圖， 其中包含金屬殼體、 防熱層原始層和防熱層燒蝕層， 如圖1所示。

從圖中可以看出， 截面共包括三層：

金屬殼體位于最內(nèi)層， 厚度沿周向相同， 其抗彎剛度計算公式中彈性模量為E1， 隨溫度t變化， 而溫度隨周向子午角θ變化， 截面慣性矩I1（x）是外徑R1和內(nèi)外徑比例α的函數(shù)， I1（x） = πR414（1-α4）， 面積A1為外徑R1和內(nèi)外徑比例α的函數(shù)， A1（x） = πR21（1-α2）。

防熱層原始層是中間層， 位于金屬殼體外部， 是防熱層除燒蝕層以外的部分， 厚度沿周向子午角θ不同， 其中由于迎風(fēng)面燒蝕嚴(yán)重， 厚度較背風(fēng)面小， 因此其外徑R2（θ）抗彎剛度計算公式中彈性模量為E2， 隨溫度t變化， 與E1類似。 截面慣性矩I2（x）和面積A2均為外徑R2（θ）的函數(shù)， 即θ的函數(shù)。

防熱層燒蝕層是最外層， 是防熱層燒蝕碳化后的部分， 厚度沿周向子午角θ不同， 其中由于迎風(fēng)面燒蝕嚴(yán)重， 厚度較背風(fēng)面大， 因此其外徑R3（θ）抗彎剛度計算公式中彈性模量為E3， 隨溫度t變化， 與E1和E2類似。 截面慣性矩I3 （x）和面積A3均為外徑R3（θ）的函數(shù)， 即θ的函數(shù)。

因此， 截面的抗彎剛度計算公式為

EI（x）= E1I1（x）+ E2I2（x）+ E3I3（x）（4）

式中：

E1I1（x）=∫AE1[t（θ）]dθ×πR414（1-α4）;

E2I2（x）=∫AE2[t（θ）]R22（θ）sin2（θ）R2（θ）dR2（θ）dθ;

E3I3（x）=∫AE3[t（θ）]R23（θ）sin2（θ）R3（θ）dR3（θ）dθ。

截面的抗剪剛度計算公式為

GA（x）= G1A1（x）+ G2A2（x）+ G3A3（x）? （5）

式中：

G1A1（x）=∫AE1[t（θ）]2（1+μ）dθ×πR21（1-α2）;

G2A2（x）=∫AE2[t（θ）]2（1+μ）R2（θ）dR2（θ）dθ;

G3A3（x）=∫AE3[t（θ）]2（1+μ）R3（θ）dR3（θ）dθ。

從而求得各截面的抗彎剛度和抗剪剛度， 代入式（3）中， 通過積分可求得梁單元的剛度矩陣， 再通過對相應(yīng)節(jié)點的剛度矩陣進(jìn)行組裝可以獲得飛行器結(jié)構(gòu)的總剛度矩陣。 實際使用時， 需要已知截面各層的溫度隨周向子午角的分布， 以及各層材料隨溫度變化的力學(xué)性能。

3 應(yīng)用程序

根據(jù)上述單元剛度和質(zhì)量矩陣的計算模型和計算公式， 具體使用的計算程序和流程如圖2所示。

4 算? 例

由于燒蝕對飛行器整體質(zhì)量的變化影響較小， 因此模態(tài)的變化主要取決非于飛行器整體的剛度分布。 本節(jié)采用給出的新模型計算截面抗彎剛度， 并與傳統(tǒng)方法進(jìn)行對比。

假設(shè)截面的參數(shù)如表1所示。 其中， 由于金屬殼體的傳熱系數(shù)較大， 算例中假設(shè)整個殼體溫度空間分布是一致的， 另外傳統(tǒng)模型僅為原始層模型， 因此傳統(tǒng)模型中防熱層剛度為考慮溫度和燒蝕對剛度的影響。

新模型采用本文的計算方法， 與傳統(tǒng)方法相比， 金屬殼體剛度計算方法一致。 考慮溫度效應(yīng)的影響， 防熱層將原始層和燒蝕層合并， 計算參數(shù)采用原始層材料參數(shù)。 兩種方法計算得到的截面抗彎剛度和抗剪剛度的結(jié)果如表2所示。

從計算結(jié)果可以看出， 新模型計算的截面抗彎剛度和抗剪剛度較傳統(tǒng)模型偏小， 相對差值分別為3.38%和2.69%。 一方面， 這是由于新模型對防熱層原始層和燒蝕層尺寸建模細(xì)化， 這兩層的物理參數(shù)存在一定的差異; 另一方面， 防熱層細(xì)化模型考慮了溫度的影響， 這說明了有必要對防熱層的剛度進(jìn)行細(xì)化建模。

5 結(jié)? 論

針對全大氣層滑翔的高超聲速飛行器動特性預(yù)示需求， 本文提出了考慮防熱層以及燒蝕的高超聲速飛行器截面剛度計算模型， 通過分析得出以下結(jié)論：

（1） 提出了3層組合的高超聲速飛行器截面剛度計算模型， 分別為最內(nèi)層承力金屬殼體、 中間防熱層原始層和最外層的防熱層燒蝕層。

（2） 根據(jù)剛度線性疊加原理， 給出了各層以及總的截面剛度參數(shù)的計算方法和具體計算公式。

（3） 針對實際使用， 給出了高超聲速飛行器結(jié)構(gòu)動特性具體的計算分析程序和步驟， 從材料參數(shù)獲取直到飛行器結(jié)構(gòu)動特性計算。

（4） 由于新模型對防熱層進(jìn)行細(xì)化建模， 而且考慮了尺寸和溫度對防熱層剛度的影響， 新模型計算的截面抗彎剛度和抗剪剛度較傳統(tǒng)模型偏小。

通過本研究， 細(xì)化了全大氣層滑翔的高超聲速飛行器動特性預(yù)示模型。 不僅可以獲得高超聲速飛行器在各飛行時間點下的動特性數(shù)值， 而且可以建立隨飛行彈道和軌跡的動特性模型， 為姿態(tài)控制提供了精細(xì)的飛行器彈性模型。

參考文獻(xiàn)：

[1] Cleghorn W L， Tabarrok B. Finite Element Formulation of a Tapered Timoshenko Beam for Free Lateral Vibration Analysis[J]. Journal of Sound and Vibration， 1992， 152（3）： 461-470.

[2] Al-Gahtani H J. Exact Stiffnesses for Tapered Members[J]. Journal of Structural Engineering， 1996， 122（10）： 1234-1239.

[3] Berman A. Mass Matrix Correction Using an Incomplete Set of Mea-sured Modes[J]. AIAA Journal， 1979， 17（10）： 1147-1148.

[4] 陸念力， 張宏生. 計及二階效應(yīng)的一種變截面梁精確單元剛度陣[J]. 工程力學(xué)， 2008， 25（12）： 60-64.

Lu Nianli， Zhang Hongsheng. Exact Elemental Stiffness Matrix of a Tapered Beam Considering Second-Order Effects[J]. Engineering Mechanics， 2008， 25（12）： 60-64.（in Chinese）

[5] 周凌遠(yuǎn)， 李彤梅， 李喬. 基于網(wǎng)格截面法的非線性梁單元[J]. 西南交通大學(xué)學(xué)報， 2009， 44（5）： 726-731.

Zhou Lingyuan， Li Tongmei， Li Qiao. Nonlinear Beam Element with Meshed Sections[J]. Journal of Southwest Jiaotong University， 2009， 44（5）： 726-731.（in Chinese）

[6] 周世軍， 朱晞. 一組新的Timoshenko梁單元一致矩陣公式[J]. 蘭州鐵道學(xué)院學(xué)報， 1994， 13（2）： 1-7.

Zhou Shijun， Zhu Xi. A Set of New Consistent Matrix for Formulations of Timoshenko Beam Element[J]. Journal of Lanzhou Railway Institute， 1994， 13（2）： 1-7.（in Chinese）

[7] 李華， 武蘭河， 黃羚. 一種新的考慮剪切變形影響的梁單元[J]. 石家莊鐵道學(xué)院學(xué)報， 1999， 12（4）： 58-61.

Li Hua， Wu Lanhe， Huang Ling. A New Element in Consideration? of the Effect of the Shear Deformation on Beam Deflection[J]. Journal of Shijiazhuang Railway Institute， 1999， 12（4）： 58-61.（in Chinese）

[8] 劉文煉. 一種新的Timoshenko梁單元[J]. 廣東土木與建筑， 2007， 14（12）： 13-14.

Liu Wenlian. A New Timoshenko Beam Element [J]. Guangdong Architecture Civil Engineering， 2007， 14（12）： 13-14.（in Chinese）

[9] 王樂， 王亮. 一種新的計算Timoshenko梁截面剪切系數(shù)的方法[J]. 應(yīng)用數(shù)學(xué)和力學(xué)， 2013， 34（7）： 756-763.

Wang Le， Wang Liang. A New Method of Obtaining Timoshenko’s Shear Coefficients[J]. Applied Mathematics and Mechanics， 2013， 34（7）： 756-763.（in Chinese）

[10] 何成. 高溫環(huán)境下結(jié)構(gòu)動力學(xué)建模關(guān)鍵技術(shù)研究[D]. 南京： 南京航空航天大學(xué)， 2014.

He Cheng. Key Technology of Dynamic Modeling in High-Tempera-ture Environment[D]. Nanjing： Nanjing University of Aeronautics and Astronautics， 2014. （in Chinese）

[11] 王亮， 商霖， 牛智玲， 等. 集中質(zhì)量法與一致質(zhì)量法對導(dǎo)彈模態(tài)計算的影響研究[J]. 強度與環(huán)境， 2013， 40（2）： 28-32.

Wang Liang， Shang Lin， Niu Zhiling， et al. Study on the Influence of the Consistent Mass Method and Lumped Mass Method in the Mode of the Rocket[J]. Structure & Environment Engineering， 2013， 40（2）： 28-32.（in Chinese）

[12] 王亮， 商霖， 王樂. 戰(zhàn)術(shù)導(dǎo)彈結(jié)構(gòu)動力學(xué)建模研究[J]. 戰(zhàn)術(shù)導(dǎo)彈技術(shù)， 2013（1）： 24-27.

Wang Liang， Shang Lin， Wang Le. Study on the Structure Dynamic Modeling of Tactical Missile[J]. Tactical Missile Technology， 2013（1）： 24-27.（in Chinese）

[13] 王亮， 周曉麗， 廖選平， 等. 飛行器扭轉(zhuǎn)模態(tài)建模分析研究[J]. 戰(zhàn)術(shù)導(dǎo)彈技術(shù)， 2015（6）： 53-56.

Wang Liang， Zhou Xiaoli， Liao Xuanping， et al. Study of the Structure Dynamic Modeling of Aircraft’s Torsion Mode[J]. Tactical Missile Technology， 2015（6）： 53-56.（in Chinese）

[14] 張琪， 劉莉. 導(dǎo)彈等效密度有限元模型的模態(tài)分析[J]. 兵工學(xué)報， 2008， 29（11）： 1395-1399.

Zhang Qi， Liu Li. Modal Analysis of Missile’s Equivalent Density Finite Element Model[J]. Acta Armamentarii， 2008， 29（11）： 1395-1399.（in Chinese）

Study on ?Dynamic Model of? Hypersonic

Aircraft Considering TPS with Ablation Effects

Wang Liang*， Zhang Yan， Cai Yipeng， Cai Wenjie

（China Academy of Launch Vehicle Technology， Beijing 100076， China）

Abstract： To indicate the dynamic characteristics of? hypersonic aircraft gliding in the whole atmosphere， a calculation model of section stiffness considering thermal protection system（TPS） and ablation effects is proposed. Firstly， the section stiffness model with three-layer combination is put forward， which include the inside main backbone， the in-between original layer of TPS， and the outmost ablative layer of TPS. Secondly，

based on the linear superposition principle，

the calculation methods of each layer and total section stiffness parameters are given.? At last， the analysis strategy of the hypersonic aircraft’s dynamic characteristics is given for actual use.

Key words： hypersonic; dynamic characteristics; thermal protection system（TPS）; ablation effect; flexural stiffness