鄧 斌，申志彬，段靜波，李 東

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

0 引言

熱載荷是固體發(fā)動(dòng)機(jī)全壽命周期的主要載荷之一［1］。某固體發(fā)動(dòng)機(jī)從生產(chǎn)到服役階段，要經(jīng)歷固化降溫、低溫試驗(yàn)及長期貯存過程，此期間與環(huán)境不斷發(fā)生熱交換，并產(chǎn)生變化溫度及不均勻溫度場，同時(shí)發(fā)動(dòng)機(jī)殼體與襯層及推進(jìn)劑等材料的熱膨脹系數(shù)至少相差一個(gè)量級，這將使發(fā)動(dòng)機(jī)承受熱載荷作用，導(dǎo)致藥柱及各粘接界面產(chǎn)生較大熱應(yīng)力，嚴(yán)重影響發(fā)動(dòng)機(jī)的使用可靠性。因此，準(zhǔn)確分析熱載荷的影響，對發(fā)動(dòng)機(jī)結(jié)構(gòu)完整性分析及壽命預(yù)估都具有重要意義。

針對熱載荷引起固體發(fā)動(dòng)機(jī)結(jié)構(gòu)響應(yīng)問題，人們普遍采用簡化的分析方法。采用均勻溫度場變化方式，雖能基本反映結(jié)構(gòu)整體溫度變化趨勢，但無法反映瞬態(tài)傳熱過程中藥柱內(nèi)部不均勻溫度場引起的結(jié)構(gòu)響應(yīng)［2］。為了分析瞬態(tài)溫度場對發(fā)動(dòng)機(jī)結(jié)構(gòu)響應(yīng)的影響，目前大多采用熱力耦合分析方法，但普遍忽略了對流換熱的影響［3－4］。由于推進(jìn)劑藥柱屬粘彈性材料，其結(jié)構(gòu)響應(yīng)與載荷歷程緊密相關(guān)，一般在室內(nèi)溫度環(huán)境，物體與環(huán)境主要以對流形式進(jìn)行換熱［5］，忽略對流換熱影響，將導(dǎo)致響應(yīng)結(jié)果的較大誤差。因此，準(zhǔn)確分析熱載荷對粘彈性藥柱的歷史作用，應(yīng)充分考慮這種影響。

目前，考慮對流換熱因素的熱載荷對固體發(fā)動(dòng)機(jī)結(jié)構(gòu)響應(yīng)的影響研究在國內(nèi)少有報(bào)道。本文利用熱力耦合有限元法，就對流換熱因素對發(fā)動(dòng)機(jī)結(jié)構(gòu)響應(yīng)的影響展開研究，旨在為今后進(jìn)一步準(zhǔn)確分析長期貯存過程中熱載荷對藥柱結(jié)構(gòu)累積損傷的影響提供指導(dǎo)。

1 熱力耦合計(jì)算模型

1.1 熱粘彈性本構(gòu)關(guān)系

假設(shè)固體發(fā)動(dòng)機(jī)藥柱為各向同性的熱流變簡單材料，積分型線性熱粘彈本構(gòu)關(guān)系可寫成［6］

式中 K與G分別為體積模量與剪切模量，都是等效時(shí)間 ξ的函數(shù);Sij、eij、σkk、εkk分別為偏應(yīng)力、偏應(yīng)變、球應(yīng)力和球應(yīng)變。

式中 aT為時(shí)溫轉(zhuǎn)換因子。

1.2 熱力耦合有限元方程

基于能量守恒定律和傅里葉定律建立的無內(nèi)熱源瞬態(tài)導(dǎo)熱微分方程為［6］

式中 T為溫度;λ為導(dǎo)熱系數(shù);ρ為材料密度;c為比熱容。

為獲得熱平衡方程的解，需定義邊界條件和初始條件，且要同時(shí)求解。熱問題的基本有限元方程可由熱平衡方程推導(dǎo)求得［7］

當(dāng)熱邊界為對流換熱方式時(shí)，式(5)中的{Q}=∫Γh{Tf}{N}dΓ。

熱力耦合計(jì)算的基本方程為［7］

2 有限元模型及分析方法

2.1 有限元模型

以某雙傘盤圓管固體發(fā)動(dòng)機(jī)為研究對象，考慮發(fā)動(dòng)機(jī)推進(jìn)劑藥柱、包覆層、絕熱層及殼體等主要部件，并根據(jù)其結(jié)構(gòu)對稱特點(diǎn)創(chuàng)建了1/12結(jié)構(gòu)三維有限元模型，共劃分24 744個(gè)六面體單元，31 010個(gè)節(jié)點(diǎn);其中，模型前后端包覆層與絕熱層之間設(shè)置了脫粘層，如圖1所示。其中，整體坐標(biāo)系O－XYZ以殼體前裙根部環(huán)所圍圓的圓心為原點(diǎn)O，Z軸自前向后過發(fā)動(dòng)機(jī)軸線，X軸位于水平面內(nèi)，并垂直于Z軸，Y軸垂直于XZ平面向上，按右手法則構(gòu)成O－XYZ整體直角坐標(biāo)系。

2.2 邊界條件與計(jì)算工況

(1)位移邊界條件。對發(fā)動(dòng)機(jī)模型的殼體一端進(jìn)行軸向固定，并在模型兩側(cè)對稱面上，施加環(huán)向?qū)ΨQ約束。

(2)熱邊界條件。假設(shè)發(fā)動(dòng)機(jī)藥柱固化后，初始零應(yīng)力溫度為331.15 K(58℃)，發(fā)動(dòng)機(jī)與外界空氣的熱交換過程，考慮其殼體外表面與周圍空氣進(jìn)行對流換熱。

計(jì)算時(shí)根據(jù)具體情況，考慮以下3種載荷歷程:

(1)從初始零應(yīng)力溫度331.15 K(58℃)置于293.15 K(20℃)環(huán)境，進(jìn)行48 h的固化降溫;

(2)置于223.15 K(－50℃)環(huán)境，進(jìn)行48 h的低溫試驗(yàn);

(3)置于223.15 K環(huán)境，進(jìn)行48 h的低溫恒溫貯存。

計(jì)算時(shí)，取溫度331.15 K(58℃)為參考溫度。

2.3 相關(guān)材料性能參數(shù)

該型發(fā)動(dòng)機(jī)的藥柱為丁羥推進(jìn)劑，其剪切松弛模量為

參考溫度為293.15 K時(shí)，時(shí)溫等效WLF方程為

其他發(fā)動(dòng)機(jī)材料參數(shù)如表1所示。

表1 發(fā)動(dòng)機(jī)材料性能參數(shù)Table 1 Material performance parameters of SRM

3 數(shù)值計(jì)算及其結(jié)果分析

3.1 2種熱邊界處理方式的計(jì)算結(jié)果對比分析

固體發(fā)動(dòng)機(jī)與環(huán)境的換熱過程包含:

(1)發(fā)動(dòng)機(jī)內(nèi)部固體介質(zhì)的熱傳導(dǎo);

(2)殼體外表面與空氣的對流換熱;

(3)殼體外表面與環(huán)境的輻射換熱。

在一般的室內(nèi)環(huán)境，發(fā)動(dòng)機(jī)與環(huán)境溫差不大，其輻射換熱量較小，可忽略不計(jì)。因此，本文傳熱計(jì)算只考慮對流換熱的主要熱交換方式。

工程上關(guān)于傳熱的計(jì)算一般采用簡化方式。2種熱邊界簡化處理方式如表2所示。Mode I為文獻(xiàn)［2］給出的處理方式，Mode II為本文擬采用的處理方式。發(fā)動(dòng)機(jī)殼體外形呈圓管狀，根據(jù)對流換熱系數(shù)的經(jīng)驗(yàn)計(jì)算方法［6］，得到上述降溫過程的自然對流換熱系數(shù)大致范圍為1.0 ～10.0 W/(m2·K)，為獲取上述2 種方式引起的計(jì)算結(jié)果差別，這里取最大的10.0 W/(m2·K)。利用2.1節(jié)的計(jì)算模型，分別對上述2種方式進(jìn)行固化降溫過程的熱力耦合計(jì)算。

從圖2所示的藥柱局部部位的溫度與應(yīng)變隨時(shí)間變化結(jié)果可知，上述2種處理方式下的計(jì)算結(jié)果變化趨勢一致，溫度及應(yīng)變都在48 h基本達(dá)到平衡，但瞬態(tài)傳熱過程中的溫度與應(yīng)變結(jié)果存在顯著差別。取藥柱內(nèi)表面的部位C為例，此2種情況在同一時(shí)刻的最大溫度差值高達(dá)9 K，最大Von－Mises應(yīng)變相對差值在某些時(shí)刻甚至超過了30%，且實(shí)際對流換熱系數(shù)一般比10.0 W/(m2·K)要小，這種差別也將會(huì)更大。

Mode I為直接在表面施加環(huán)境溫度，相當(dāng)于認(rèn)為對流換熱系數(shù)無窮大，即對流換熱熱阻為0，這種簡化假設(shè)方式是不準(zhǔn)確的。分析發(fā)動(dòng)機(jī)與外界的真實(shí)傳熱過程可知，殼體與空氣之間存在一定的換熱熱阻，發(fā)動(dòng)機(jī)外表面溫度不可能立即達(dá)到環(huán)境溫度，而是主要通過對流等換熱方式與空氣進(jìn)行換熱［6］。因此，考慮對流換熱作用能更精確地計(jì)算得到熱載荷作用下結(jié)構(gòu)響應(yīng)歷程，這對進(jìn)一步研究載荷歷史對藥柱結(jié)構(gòu)的累積損傷作用至關(guān)重要。因此，在換熱過程的結(jié)構(gòu)分析時(shí)，采用Mode II所示的對流換熱方式進(jìn)行計(jì)算。

表2 2種熱邊界簡化假設(shè)Table 2 Two modes for thermal boundary predigestion

3.2 對流換熱系數(shù)對發(fā)動(dòng)機(jī)結(jié)構(gòu)的影響

固體推進(jìn)劑由粘彈性材料組成，其力學(xué)性質(zhì)依賴于時(shí)間、溫度等因素，其結(jié)構(gòu)響應(yīng)不僅與當(dāng)前應(yīng)力狀態(tài)有關(guān)，而且與以往歷史相關(guān)。固體發(fā)動(dòng)機(jī)主要通過殼體與周圍空氣的自然對流換熱進(jìn)行熱交換，然而物體表面的自然對流換熱系數(shù)受到溫度、重力加速度、表面與空氣間的溫差、物體特征長度以及運(yùn)動(dòng)粘度系數(shù)等外界因素的影響［5］，其實(shí)際數(shù)值在一定范圍內(nèi)波動(dòng)。為獲得自然對流換熱系數(shù)對結(jié)構(gòu)響應(yīng)的影響程度，結(jié)合自然對流換熱系數(shù)的經(jīng)驗(yàn)取值范圍，取2.0、4.0、6.0、8.0、10.0、12.0 共 6 個(gè)不同系數(shù)，分別對發(fā)動(dòng)機(jī)進(jìn)行固化降溫過程的結(jié)構(gòu)響應(yīng)分析。計(jì)算時(shí)，仍采用2.1節(jié)的計(jì)算模型。

對計(jì)算結(jié)果進(jìn)行后處理，得到應(yīng)變與溫度均較大的部位C在不同換熱系數(shù)下的溫度及等效Von－Mises應(yīng)變隨時(shí)間變化曲線，如圖3所示。

由圖3所示的計(jì)算結(jié)果可知，對流換熱系數(shù)為2.0 W/(m2·K)時(shí)，降溫速率很小，以至于在降溫過程第48 h的藥柱內(nèi)外溫度還遠(yuǎn)未達(dá)到平衡，但當(dāng)對流換熱系數(shù)大于4.0 W/(m2·K)時(shí)，隨對流換熱系數(shù)的不斷增大，即換熱熱阻不斷減小，溫度下降速率有所增大，但由此引起的溫度變化速率的敏感度下降，對流系數(shù)在6.0～12.0 W/(m2·K)變動(dòng)，引起的降溫速率變化已非常小，即對流換熱系數(shù)不斷增大，對溫度變化率的影響將逐漸減弱。

同樣，應(yīng)變也有類似的變化規(guī)律，由于零應(yīng)力溫度為較高的331.15 K，降溫過程的應(yīng)變與降溫溫差變化表現(xiàn)一致，即降溫溫差越大，其應(yīng)變越大，降溫速率快時(shí)，應(yīng)變變化率也較大，反之也對應(yīng)成立。但隨對流換熱系數(shù)的不斷增大，引起應(yīng)變率變化的影響逐漸減弱。

綜上所述，對換熱系數(shù)變化對結(jié)構(gòu)歷史響應(yīng)的影響顯著，但隨其不斷增大，這種影響逐漸減弱，在6.0～12.0 W/(m2·K)范圍引起的結(jié)果變化已很小。

3.3 考慮對流換熱的發(fā)動(dòng)機(jī)熱力耦合分析

3.3.1 計(jì)算邊界條件及分析工況

根據(jù)給出的發(fā)動(dòng)機(jī)有限元模型，采用對流換熱的熱邊界處理方式，考慮依次經(jīng)歷2.2節(jié)中的(1)～(3)載荷歷程，并進(jìn)行熱力耦合響應(yīng)計(jì)算。假設(shè)初始溫度和零應(yīng)力溫度都為331.15 K，整個(gè)過程考慮殼體外表面與空氣的對流換熱。發(fā)動(dòng)機(jī)殼體外表面為圓管狀，根據(jù)橫向圓管表面自然對流換熱系數(shù)的經(jīng)驗(yàn)計(jì)算方法［5］，可得上述3種載荷歷程下的平均換熱系數(shù)分別為 5.0、3.5、3.5 W/(m2·K)。

根據(jù)計(jì)算結(jié)果，得到圖4所示的降溫過程第96 h的溫度場和等效Von－Mises應(yīng)變場。

由于降溫過程中發(fā)動(dòng)機(jī)的內(nèi)部導(dǎo)熱及其與外界換熱過程一直處于非平衡態(tài)，其內(nèi)部溫度由此呈現(xiàn)不均勻分布狀態(tài)。為表現(xiàn)所關(guān)注部位的結(jié)果變化趨勢，下面給出其計(jì)算結(jié)果隨時(shí)間響應(yīng)曲線。根據(jù)圖4所示的結(jié)果，關(guān)注部位取藥柱中間段外表面A處、前脫粘層尖端B處、內(nèi)孔中間段表面C處以及后傘盤溝槽D處等(具體位置見圖1)，給出它們相應(yīng)的溫度及等效Von Mises應(yīng)變隨時(shí)間變化曲線，如圖5所示。

3.3.2 計(jì)算結(jié)果分析

(1)在整個(gè)降溫過程中，發(fā)動(dòng)機(jī)與外界溫差較大時(shí)，其降溫速率較大，反之也對應(yīng)成立;由于降溫過程的發(fā)動(dòng)機(jī)殼體對外進(jìn)行換熱，藥柱內(nèi)表面溫度相比于外表面的降溫要相對緩慢;藥柱應(yīng)變與降溫溫差變化趨勢一致，在降溫速率較大時(shí)，應(yīng)變率較大;在降溫速率較小時(shí)，應(yīng)變率也相應(yīng)較小(圖5)。

(2)發(fā)動(dòng)機(jī)在降溫過程中，主要以對流換熱方式與外界環(huán)境進(jìn)行換熱，并由此產(chǎn)生不均勻分布場。在經(jīng)歷充分的換熱過程后，整體溫度分布趨于一致。從圖4與圖5可看出，發(fā)動(dòng)機(jī)在固化降溫過程的第48 h接近于環(huán)境溫度293.15 K，而第96 h即經(jīng)過48 h的低溫降溫過程，發(fā)動(dòng)機(jī)最大溫度約為226 K，與環(huán)境溫度223.15 K仍存在較小溫差，再經(jīng)歷48 h的低溫恒溫保存過程的進(jìn)一步換熱后，內(nèi)外溫度則基本達(dá)到223.15 K的平衡溫度。

(3)整個(gè)降溫過程發(fā)動(dòng)機(jī)與周圍環(huán)境不斷進(jìn)行瞬態(tài)傳熱，不斷變化的溫度及溫度的不均勻分布導(dǎo)致藥柱內(nèi)部產(chǎn)生較大的熱應(yīng)變。發(fā)動(dòng)機(jī)整個(gè)降溫過程，從第48 h的331.15 K降到223.15 K，降溫溫差達(dá)108 K，使得發(fā)動(dòng)機(jī)在低溫下發(fā)生收縮，不同材料的膨脹系數(shù)巨大差異及不均勻溫度分布場，導(dǎo)致了其內(nèi)部局部區(qū)域產(chǎn)生較大應(yīng)變。隨熱應(yīng)變的不斷累積，位于藥柱內(nèi)孔表面的部位C在第144 h的等效Von Mises應(yīng)變最大達(dá)到了17.6%。由于低溫下藥柱延伸率低，結(jié)合應(yīng)變能準(zhǔn)則可知，其當(dāng)前結(jié)構(gòu)安全系數(shù)僅為2左右，發(fā)動(dòng)機(jī)長期貯存過程的結(jié)構(gòu)完整性受到嚴(yán)重考驗(yàn)。

4 結(jié)論

(1)考慮殼體對流換熱的影響，可更準(zhǔn)確地反映發(fā)動(dòng)機(jī)與外界的換熱過程，進(jìn)一步提高了熱載荷下的結(jié)構(gòu)響應(yīng)計(jì)算精度，有助于準(zhǔn)確分析熱載荷歷史對裝藥的累積損傷影響。

(2)對流換熱系數(shù)大小影響降溫過程的發(fā)動(dòng)機(jī)結(jié)構(gòu)響應(yīng)，但隨對流換熱系數(shù)的不斷增大，這種影響逐漸減弱。

(3)經(jīng)歷較大溫差降溫過程的藥柱內(nèi)表面熱應(yīng)變較大，考驗(yàn)長期貯存過程中的藥柱結(jié)構(gòu)完整性，并影響著固體發(fā)動(dòng)機(jī)的使用壽命。

［1］邢耀國，董可海，沈偉，等.固體火箭發(fā)動(dòng)機(jī)使用工程［M］.北京:國防工業(yè)出版社，2010.

［2］蒙上陽.基于粘彈性有限元方法的固體火箭發(fā)動(dòng)機(jī)結(jié)構(gòu)完整性分析［D］.長沙:國防科技大學(xué)，2005.

［3］Chyuan S W.A study of loading history effect for thermoviscoelastic solid propellant grains［J］.Computers and Structures，2000，(77):735－745.

［4］于洋.溫度載荷下帶筋套管形裝藥結(jié)構(gòu)完整性分析［J］.推進(jìn)技術(shù)，2006，27(6):493－496.

［5］楊世銘，陶文銓.傳熱學(xué)(第四版)［M］.北京:高等教育出版社，2006.

［6］楊挺青.粘彈性力學(xué)［M］.華中理工大學(xué)出版社，1990.

［7］Hu Z，Zhu L H，Wang B Y，et al.Computer simulation of the deep extrusion of a thin walled cup using the thermo mechanically coupled elastic plastic FEM［J］.Journal of Materials Processing Technology，2000，102(1):128－137.