非均勻正交各向異性端頭帽熱應(yīng)力分析①

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

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

0 引言

高超聲速飛行器在出入大氣層或持續(xù)在大氣內(nèi)飛行時(shí)，將承受巨大的氣動力和氣動熱。氣動熱載荷表現(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 計(jì)算模型

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(shí)測得的熱模量，是一個對稱張量;T為當(dāng)前溫度;T0為參考溫度，即cijkl測定時(shí)的溫度。

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

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

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

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

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

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

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

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

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

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

2 數(shù)值方法

2.1 順序耦合算法

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

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

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

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

端頭帽為球錐結(jié)構(gòu)，如圖1所示，球頭半徑為RN，錐身半錐角θ，總長L。由對稱性，考慮計(jì)算工況含有攻角、材料的正交各向異性，數(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ù)以往雙向拉伸性能試驗(yàn)，可判斷在與炭布x向或z向的夾角為45°時(shí)，材料的模量和強(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°時(shí)，材料力學(xué)性能隨αg增大而減小;至αg=45°時(shí)，力學(xué)性能減小至x方向的70%;當(dāng)45°＜αg＜90°時(shí)，材料力學(xué)性能隨αg增大而增大;至αg=90°時(shí)，力學(xué)性能增大至100%，與x方向性能一致。本文采用線性插值方法獲取αg在不同角度時(shí)的端頭帽力學(xué)性能。端頭帽在主軸方向的力學(xué)性能如表2所示。

2.4 初邊值條件

端頭帽彈道曲線如圖4所示。設(shè)端頭帽初始溫度均勻分布，取300 K，計(jì)算總時(shí)間為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)熱邊界為端頭帽外表面，受氣動加熱作用，同時(shí)向外部環(huán)境空間輻射熱量。

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

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

圖6為0°攻角典型時(shí)刻的溫度云圖。由于熱導(dǎo)率較小，在870 s以前，高溫區(qū)主要集中在端頭帽球頭，溫度梯度較大，而在錐身部溫度分布較為均勻、梯度較小。在870 s后，高溫區(qū)向錐身擴(kuò)展，但此時(shí)最高溫度已小于1 200 K，最大溫差小于50 K，整體分布較均勻，溫度梯度較小。另外，由于輻射換熱與表面溫度的四次方成正比，在260 s以前，輻射換熱小于氣動加熱，駐點(diǎn)溫度急劇升高;之后，輻射換熱大于氣動加熱，駐點(diǎn)溫度逐漸減小，駐點(diǎn)位置由向內(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)。因此，以下重點(diǎn)對端頭帽的應(yīng)變和Mises等效應(yīng)力進(jìn)行分析。

圖7為0°攻角典型時(shí)刻主彈性應(yīng)變云圖，圖8為0°攻角典型時(shí)刻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)橹鳎溆嘀黧w部分仍以熱應(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°攻角典型時(shí)刻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)力的影響

受攻角影響，頭部駐點(diǎn)下移(最大溫度2 351 K)，低溫點(diǎn)上移，在尾部背風(fēng)位置附近。高溫區(qū)域仍集中在頭部，由于頭部為球形，其以駐點(diǎn)為中心的分布趨勢與0°攻角時(shí)基本相同。圖10為10°攻角母線溫度分布曲線。由圖10可知，10°攻角時(shí)，迎風(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)變最大值在駐點(diǎn)附近;主彈性應(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)力軸向分布與無攻角時(shí)基本相同:熱應(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，較無攻角時(shí)(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è)計(jì)時(shí)，應(yīng)將材料的xz平面45°方向設(shè)置在迎風(fēng)線上。這樣可減小端頭帽應(yīng)力極值，使端頭帽周向應(yīng)力分布趨于均勻，有利于端頭帽與彈體連接。

4 結(jié)論

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

(2)彈道初段，局部大熱膨脹拉伸頭部，在頭部形成環(huán)狀彈性應(yīng)變區(qū)。在彈道中、末段，高溫區(qū)向后移動，總應(yīng)變在端頭帽尾部以固定約束條件產(chǎn)生的彈性應(yīng)變?yōu)橹鳎溆嘀黧w部分仍以熱應(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è)計(jì)時(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é)報(bào)，2008，20(15):3966-3968.

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