程興華，李建林，楊 濤，李 理

(國防科技大學(xué)航天與材料工程學(xué)院，長沙 410073)

0 引言

高超聲速飛行器在出入大氣層或持續(xù)在大氣內(nèi)飛行時，將承受巨大的氣動力和氣動熱。氣動熱載荷表現(xiàn)為表面與大氣的摩擦而達(dá)到較高的溫度，尤其是飛行器頭錐、翼前緣等部位。氣動熱不僅使材料的力學(xué)性能降低，而且在結(jié)構(gòu)中產(chǎn)生熱應(yīng)力，從而使變形加劇，高溫?zé)彷d荷引起的熱應(yīng)力不容忽視。以往對高溫區(qū)的溫度、熱應(yīng)力分析多為非線性、正交各向異性，而對復(fù)合材料因編織方法引起的物性非均勻性分析很少[1－5]。

本文以正交各向異性復(fù)合材料端頭帽為對象，研究非均勻材料物性引起的應(yīng)力分布變化。

1 計算模型

1.1 熱彈性本構(gòu)方程

在彈性力學(xué)中，應(yīng)力和應(yīng)變關(guān)系由廣義胡克定律確定。胡克定律是固體彈性材料在均勻常值溫度下的本構(gòu)方程，如考慮溫度變化的影響，而不考慮溫度與變形之間的耦合關(guān)系，則材料的本構(gòu)方程[6]可記為

式中cijkl為彈性模量，也稱剛度系數(shù);εkl為應(yīng)變張量;βij是在應(yīng)變?yōu)榱銜r測得的熱模量，是一個對稱張量;T為當(dāng)前溫度;T0為參考溫度，即cijkl測定時的溫度。

式(1)中的溫度變化ΔT一般用符號Θ來表示，這樣式(1)即可寫作:

式(2)稱為杜哈梅爾-牛曼(Duhamel-Neumann)關(guān)系式。

1.2 瞬態(tài)傳熱本構(gòu)方程

熱彈性理論離不開對傳熱規(guī)律的了解，傳熱有傳導(dǎo)、輻射和對流3種形式。其基本的能量守恒方程為

式中V為固體材料的體積;S為表面積;ρ為材料密度;U·為內(nèi)能的時間變化速率;q為單位表面積上流入控制體的熱流率;r為單位體積內(nèi)的其他熱增量。

假設(shè)U=U(T)，T為材料的溫度，q和r與材料的應(yīng)力和位移無關(guān)。除相變產(chǎn)生的潛熱外，傳熱分析的本構(gòu)方程為

其中，c(T)為材料的比熱容，當(dāng)給定潛熱時，認(rèn)為它是比熱容的一個附加量。多數(shù)情況下，認(rèn)為相變在一定的溫度范圍內(nèi)發(fā)生是合理的。

導(dǎo)熱遵守傅里葉定律:

式中 λi為熱導(dǎo)率。

物體中的溫度分布，無論是定常還是非定常，都需要給定邊界條件[7－8]，對于非定常情況還須給出時間的初始條件，這才能求解方程(3)內(nèi)的溫度分布。

2 數(shù)值方法

2.1 順序耦合算法

持續(xù)受熱的特征決定了飛行器在某一時刻的受熱情況不僅與這一時刻的飛行條件有關(guān)，還與這之前飛行器自身熱環(huán)境所經(jīng)歷的時間歷程緊密相關(guān)。在以往氣動加熱分析中，多數(shù)研究的是瞬態(tài)情況，即在一定的飛行(來流)條件和外表面壁溫條件下給出外表面的熱流分布。這種瞬態(tài)結(jié)果尚不足以用來判斷飛行器在持續(xù)受熱過程中的各種中間狀態(tài)，因為在各個中間狀態(tài)對應(yīng)的物體自身熱環(huán)境都是經(jīng)歷了一定時間歷程后形成的，考慮各傳熱環(huán)節(jié)之間的耦合性十分必要。

順序耦合算法又稱松耦合方法，其特點是流場與結(jié)構(gòu)計算域在接觸面上周期性交換數(shù)據(jù)，主要體現(xiàn)為流體與結(jié)構(gòu)接觸面上的邊界條件處理。

本文通過Fortran用戶子程序?qū)Υ笮头蔷€性有限元計算軟件Abaqus進(jìn)行二次開發(fā)。每個時間增量步開始時，用戶子程序讀取端頭帽外表面上的溫度，將計算得到的氣動熱流傳遞給Abaqus主程序;Abaqus主程序根據(jù)邊界熱流進(jìn)行熱應(yīng)力計算;如此順序耦合循環(huán)，直至彈道結(jié)束點，最終完成端頭帽熱應(yīng)力分析。

2.2 物理數(shù)學(xué)模型

端頭帽為球錐結(jié)構(gòu)，如圖1所示，球頭半徑為RN，錐身半錐角θ，總長L。由對稱性，考慮計算工況含有攻角、材料的正交各向異性，數(shù)學(xué)模型取周向180°。xy平面為端頭帽對稱面，α為攻角，αg為子午面與xy平面的夾角。

網(wǎng)格劃分如圖2所示。由于頭部熱流梯度比較大，對頭部網(wǎng)格進(jìn)行了適當(dāng)加密。網(wǎng)格單元選取耦合溫度——位移單元C3D8T和C3D6T。

圖1 模型三維視圖Fig.1 3D view of them odel

圖2 網(wǎng)格劃分Fig.2 M esh for discretization

2.3 材料模型

端頭帽材料以周向0°、90°正交編織，材料沿端頭帽軸向性能一致。由于x、z方向(與端頭帽軸線垂直的平面)交錯疊加，消除了x、z方向性能差別，但不能消除xz平面內(nèi)45°方向的性能差異，即材料性能在xz平面內(nèi)呈現(xiàn)非均勻性。根據(jù)測試結(jié)果，在45°方向主要是力學(xué)性能(彈性模量E)的差異，熱物理性能(表1)差異很小，在本文研究中可忽略。

根據(jù)以往雙向拉伸性能試驗，可判斷在與炭布x向或z向的夾角為45°時，材料的模量和強(qiáng)度減小到x向或z向的70%。

端頭帽材料模型如圖3所示。材料力學(xué)性能在xz平面內(nèi)繞y軸以90°周期分布。在第Ⅰ象限內(nèi)，以材料在x方向的力學(xué)性能為標(biāo)準(zhǔn)(100%);當(dāng)0°＜αg＜45°時，材料力學(xué)性能隨αg增大而減小;至αg=45°時，力學(xué)性能減小至x方向的70%;當(dāng)45°＜αg＜90°時，材料力學(xué)性能隨αg增大而增大;至αg=90°時，力學(xué)性能增大至100%，與x方向性能一致。本文采用線性插值方法獲取αg在不同角度時的端頭帽力學(xué)性能。端頭帽在主軸方向的力學(xué)性能如表2所示。

2.4 初邊值條件

端頭帽彈道曲線如圖4所示。設(shè)端頭帽初始溫度均勻分布，取300 K，計算總時間為977 s。

表1 端頭帽材料熱物理性能Table 1 Thermal physical properties of the nosecap material

圖3 力學(xué)性能插值示意圖Fig.3 Interpolation sketch for mechanics properties

表2 端頭帽材料力學(xué)性能Table 2 M echanics properties of the nosecap material

圖4 端頭帽彈道曲線Fig.4 Flight trajectory of the nosecap

彈頭在大氣環(huán)境下高速飛行，不僅受氣動力作用，還受氣動加熱作用。端頭帽的邊界條件有熱邊界、氣動壓力邊界、對稱邊界、連接支撐邊界。

(1)熱邊界為端頭帽外表面，受氣動加熱作用，同時向外部環(huán)境空間輻射熱量。

本文采用Lees模型[9]計算氣動加熱熱流通量:

式中Q為進(jìn)入表面的熱量;S為邊界表面;t為時間;ψ(x，y，z，t)為給定的熱流通量函數(shù)。

端頭帽為熱結(jié)構(gòu)，表面輻射到環(huán)境空間的熱力密度為

式中 εr為端頭帽表面熱輻射率，一般假設(shè)εr=0.8;σr為波爾茲曼常數(shù)，σr=5.67×10－8;Tb為端頭帽的表面溫度;Tair為環(huán)境大氣溫度。

(2)氣動壓力邊界同樣作用在端頭帽的外表面，壓力值由已知條件給出。

(3)對稱邊界為端頭帽縱對稱面，在對稱面上:

式中U3為z方向的位移;UR1為繞x軸的旋轉(zhuǎn);UR2為繞y軸的旋轉(zhuǎn)。

(4)連接支撐邊界為端頭帽尾部，與彈體連接，假設(shè)連接支撐邊界為固支。

3 結(jié)果分析

3.1 溫度分布

表面節(jié)點溫度變化主要受氣動加熱和熱輻射影響，當(dāng)氣動加熱與熱輻射平衡時，節(jié)點溫度達(dá)到極大值。圖5為0°攻角母線溫度變化曲線，可見駐點附近溫差較小，駐點最高溫度為2 356 K;在彈道前段，表面溫度在球頭沿母線急劇下降，而在錐身下降緩慢;在彈道末段，由于高溫的強(qiáng)輻射作用，端頭帽球頭溫度較錐身下降快;977 s時球頭溫度低于錐身溫度，端頭帽最低溫度在駐點(約1 001 K)，最高溫度在尾部軸線上(約1 049 K)。由此可見，在彈道末段，端頭帽內(nèi)部溫度較高，表面熱輻射大于氣動加熱，端頭帽向外放熱。

圖6為0°攻角典型時刻的溫度云圖。由于熱導(dǎo)率較小，在870 s以前，高溫區(qū)主要集中在端頭帽球頭，溫度梯度較大，而在錐身部溫度分布較為均勻、梯度較小。在870 s后，高溫區(qū)向錐身擴(kuò)展，但此時最高溫度已小于1 200 K，最大溫差小于50 K，整體分布較均勻，溫度梯度較小。另外，由于輻射換熱與表面溫度的四次方成正比，在260 s以前，輻射換熱小于氣動加熱，駐點溫度急劇升高;之后，輻射換熱大于氣動加熱，駐點溫度逐漸減小，駐點位置由向內(nèi)加熱轉(zhuǎn)為向外放熱。

圖5 0°攻角母線溫度變化Fig.5 Temperature profile along generating line under 0°angle of attack

圖6 0°攻角溫度云圖Fig.6 Temperature distributions under 0°angle of attack

3.2 應(yīng)變-應(yīng)力分布

在復(fù)合材料的強(qiáng)度準(zhǔn)則中以最大應(yīng)變準(zhǔn)則和最大應(yīng)力準(zhǔn)則、蔡-希爾準(zhǔn)則、蔡-吳準(zhǔn)則[10]應(yīng)用的較多。Mises等效應(yīng)力可以看作蔡-吳準(zhǔn)則的一個特例，其定義為

其中，S為偏應(yīng)力張量，S=σ+p I(σ為應(yīng)力，I為單位矩陣，p=－σii/3為等效壓應(yīng)力，也就是常見的p=(σx+σy+σz)/3)。因此，以下重點對端頭帽的應(yīng)變和Mises等效應(yīng)力進(jìn)行分析。

圖7為0°攻角典型時刻主彈性應(yīng)變云圖，圖8為0°攻角典型時刻x向熱應(yīng)變云圖。由圖7、圖8可知，彈道初段，由于端頭帽頭部高溫區(qū)集中，局部大熱膨脹拉伸頭部，使得頭部出現(xiàn)環(huán)形大彈性應(yīng)變區(qū)域;彈性應(yīng)變較熱應(yīng)變?yōu)樾×?，總?yīng)變分布與熱應(yīng)變分布相同，大小基本相當(dāng)。在彈道中段、末段，高溫區(qū)向后移動，端頭帽尾部固定端的彈性應(yīng)變凸顯出來，總應(yīng)變在端頭帽尾部以約束條件產(chǎn)生的彈性應(yīng)變?yōu)橹?，而其余主體部分仍以熱應(yīng)力為主。

圖7 0°攻角主彈性應(yīng)變云圖Fig.7 Principal strains distributions under 0°angle of attack

圖8 0°攻角x向熱應(yīng)變云圖(t=120 s)Fig.8 x direction thermal strain distributions under 0°angle of attack(t=120 s)

圖9為0°攻角典型時刻Mises應(yīng)力云圖。由圖9可知，在彈道初段，由于熱膨脹作用，端頭帽外表面形成環(huán)狀Mises應(yīng)力分布，內(nèi)部也存在較大應(yīng)力;在彈道末段，熱應(yīng)力較小，尾部由于固定約束而形成較大的應(yīng)力。在熱應(yīng)力集中段，外表面Mises應(yīng)力最大;在尾部約束處，內(nèi)部Mises應(yīng)力最大。在垂直于y軸的剖面上，Mises應(yīng)力分布以90°為周期(如圖9中的尾端面應(yīng)力云圖)，在0°、90°方向最大，逐漸過渡到最小的45°方向。

3.3 攻角對溫度、應(yīng)力的影響

受攻角影響，頭部駐點下移(最大溫度2 351 K)，低溫點上移，在尾部背風(fēng)位置附近。高溫區(qū)域仍集中在頭部，由于頭部為球形，其以駐點為中心的分布趨勢與0°攻角時基本相同。圖10為10°攻角母線溫度分布曲線。由圖10可知，10°攻角時，迎風(fēng)母線(αg=180°)溫度最高，側(cè)面水平母線(αg=90°)近似為迎風(fēng)和背風(fēng)母線(αg=0°)的平均值。

圖9 0°攻角M ises應(yīng)力云圖Fig.9 M ises equivalent stress distributions under 0°angle of attack

圖10 10°攻角母線溫度分布(t=260 s)Fig.10 Tem perature profile along generating line under 10°angle of attack(t=260 s)

圖11為10°攻角x向熱應(yīng)變云圖。由圖11可見，除約束端外，10°攻角端頭帽總應(yīng)變在彈道初段仍以熱應(yīng)變?yōu)橹?。熱?yīng)變集中在球頭部，以攻角α=10°的球頭徑向?qū)ΨQ分布;總應(yīng)變最大值在駐點附近;主彈性應(yīng)變在迎風(fēng)面較大，高彈性應(yīng)變區(qū)域在迎風(fēng)面較寬。

對應(yīng)于應(yīng)變分布，彈道初段Mises應(yīng)力在迎風(fēng)面上出現(xiàn)極大值，且呈現(xiàn)較寬區(qū)域的應(yīng)力集中(圖12)。彈道末段，10°攻角Mises應(yīng)力軸向分布與無攻角時基本相同:熱應(yīng)力較小，尾部由于固定約束而形成較大的應(yīng)力。受攻角和材料性能參數(shù)的雙重影響，在垂直于y軸的剖面上，Mises應(yīng)力在αg=180°子午面內(nèi)最大，αg=45°子午面內(nèi)最小(圖13)，這種現(xiàn)象在彈道初段最為明顯。分析可知，最大Mises應(yīng)力為111 MPa，較無攻角時(104 MPa)略大。

圖11 10°攻角x向熱應(yīng)變云圖(t=260 s)Fig.11 x direction thermal strains distributions under 10°angle of attack(t=260 s)

圖12 10°攻角M ises應(yīng)力云圖(t=110 s)Fig.12 M ises equivalent stress distributions under 10°angle of attack(t=110 s)

圖13 10°攻角尾部外表面周向M ises應(yīng)力分布Fig.13 Circum ferential M ises stress distributions on tail surface under 10°angle of attack

在端頭帽結(jié)構(gòu)設(shè)計時，應(yīng)將材料的xz平面45°方向設(shè)置在迎風(fēng)線上。這樣可減小端頭帽應(yīng)力極值，使端頭帽周向應(yīng)力分布趨于均勻，有利于端頭帽與彈體連接。

4 結(jié)論

(1)彈道初段，表面熱輻射小于氣動加熱，端頭帽溫度迅速升高;彈道末段，端頭帽內(nèi)部溫度較高，表面熱輻射大于氣動加熱，端頭帽向外放熱。有攻角時，側(cè)面水平母線溫度近似為迎風(fēng)和背風(fēng)母線平均值。

(2)彈道初段，局部大熱膨脹拉伸頭部，在頭部形成環(huán)狀彈性應(yīng)變區(qū)。在彈道中、末段，高溫區(qū)向后移動，總應(yīng)變在端頭帽尾部以固定約束條件產(chǎn)生的彈性應(yīng)變?yōu)橹?，而其余主體部分仍以熱應(yīng)力為主。因此，端頭帽尾端的連接方式十分重要。

(3)在垂直于端頭帽y軸的剖面上，Mises應(yīng)力分布以90°為周期，在0°、90°方向最大，逐漸過渡到較小的45°方向。

(4)受攻角影響，Mises應(yīng)力在迎風(fēng)子午面內(nèi)最大，水平子午面內(nèi)最小，這種現(xiàn)象在彈道初段最為明顯。在端頭帽結(jié)構(gòu)設(shè)計時，應(yīng)將材料的xz平面45°方向設(shè)置在迎風(fēng)線上。

[1] Grant Palmer.A heating analysis of the nosecap and leading edges of the X-34 vehicle[R].AIAA 98-0878.

[2] 尹蓮花，劉莉，張帥.燒蝕材料的近程高速彈頭熱應(yīng)力數(shù)值仿真研究[J].系統(tǒng)仿真學(xué)報，2008，20(15):3966-3968.

[3] 姜志杰，張擘毅，何浩.高超聲速飛行器鼻錐的熱環(huán)境和結(jié)構(gòu)熱分析研究[J].導(dǎo)彈與航天運載技術(shù)，2009，(4):14-17.

[4] 易龍，孫秦，彭云.復(fù)合材料頭錐結(jié)構(gòu)氣動熱應(yīng)力分析方法研究[J].機(jī)械強(qiáng)度，2007，29(4):686-692.

[5] Rosario Borrellia，Aniello Riccio，Domenico Tescione.Thermo-structural behaviour of an UHTCmade nosecap of a reentry vehicle[J].Acta Astronautica，2009，65:442-456.

[6] 范緒箕.高速飛行器熱結(jié)構(gòu)分析與應(yīng)用[M].北京:國防工業(yè)出版社，2009.

[7] 朱濱.彈性力學(xué)[M].合肥:中國科學(xué)技術(shù)大學(xué)出版社，2008.

[8] 趙鎮(zhèn)南.傳熱學(xué)[M].北京:高等教育出版社，2002.

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

[10] 王寶來，吳世平，梁軍.復(fù)合材料失效及其強(qiáng)度理論[J].失效分析與預(yù)防，2006，1(2):13-19.