耿發(fā)貴 ，安自朝 ，何坤 ，康健 ，張鵬坤 ，劉騰躍 ，肖航

（1. 四川航天系統(tǒng)工程研究所，四川 成都 610000；2. 某部駐成都地區(qū)第三軍事代表室，四川 成都 610041）

0 引言

隨著設(shè)計(jì)制造技術(shù)的提升以及超高強(qiáng)度新材料的涌現(xiàn)，固體火箭發(fā)動(dòng)機(jī)金屬殼體的筒體壁厚越做越薄。在火箭的飛行過程中，會(huì)出現(xiàn)固體火箭發(fā)動(dòng)機(jī)承受外壓的情況，這無疑給越做越薄的固體火箭發(fā)動(dòng)機(jī)殼體帶來全新的考驗(yàn)。

在以往固體火箭發(fā)動(dòng)機(jī)殼體外壓屈曲的研究中，研究者的關(guān)注點(diǎn)主要在復(fù)合材料殼體以及殼體的軸壓屈曲上。文獻(xiàn)［1-2］利用有限元結(jié)構(gòu)分析法，分別對某纖維纏繞復(fù)合材料固體火箭發(fā)動(dòng)機(jī)殼體在均布外壓作用下以及軸壓作用下的穩(wěn)定性進(jìn)行了分析；文獻(xiàn)［3］介紹了靜力軸壓下圓柱殼體屈曲應(yīng)力的工程計(jì)算方法，對某固體火箭發(fā)動(dòng)機(jī)殼體在復(fù)雜受載條件下殼體的屈曲應(yīng)力進(jìn)行了分析，并對該火箭發(fā)動(dòng)機(jī)殼體屈曲應(yīng)力如何估計(jì)提出了建議。

目前國內(nèi)對光滑圓柱殼體外壓屈曲的研究較多，研究方向主要在于通過特征值計(jì)算方法獲得圓柱殼體的屈曲模態(tài)，再基于特征值屈曲模態(tài)引入初始缺陷對圓柱殼體進(jìn)行非線性屈曲分析。例如，文獻(xiàn)［4］利用ANSYS 對外壓圓柱殼進(jìn)行了特征值、幾何非線性和幾何/材料雙非線性屈曲分析；文獻(xiàn)［5］采用ANSYS 分析軟件引入初始缺陷對某壓力容器筒體進(jìn)行非線性屈曲分析；文獻(xiàn)［6］基于單一模態(tài)缺陷法和組合模態(tài)缺陷法對圓柱殼體外壓屈曲進(jìn)行了仿真分析；文獻(xiàn)［7］利用ANSYS 的幾何/材料雙非線性分析技術(shù)對外壓圓筒臨界壓力進(jìn)行計(jì)算；文獻(xiàn)［8］將線性屈曲分析的屈曲模態(tài)作為初始缺陷的參考值，對受外壓海底管道進(jìn)行非線性有限元分析。

此外，部分研究者采用理論公式對光滑圓柱殼體的外壓穩(wěn)定性能進(jìn)行了分析。文獻(xiàn)［9］采用經(jīng)驗(yàn)公式和有限元分析方法對外壓圓柱殼體的屈曲進(jìn)行了對比分析；文獻(xiàn)［10］對鋁合金殼體的外壓屈曲計(jì)算進(jìn)行了探討，分析影響圓筒失穩(wěn)臨界壓力的因素；文獻(xiàn)［11］介紹了國內(nèi)科研人員對于圓柱殼體屈曲問題的研究現(xiàn)狀以及結(jié)果應(yīng)用。

借鑒上述研究經(jīng)驗(yàn)，本文以ABAQUS 有限元軟件為分析平臺(tái)，從光殼體入手探索外壓屈曲的分析方法，再將此方法應(yīng)用于固體火箭發(fā)動(dòng)機(jī)殼體的外壓屈曲分析中，對固體火箭發(fā)動(dòng)機(jī)殼體的抗屈曲性能進(jìn)行了研究。

1 光殼體的外壓屈曲分析

1.1 模型概述

殼體筒段長徑比為12，壁厚為2 mm?？紤]一般殼體兩端為封頭或法蘭等加強(qiáng)結(jié)構(gòu)，約束殼體兩端的圓度。殼體材料采用D406A，其彈性模量為210 GPa，泊松比取0.3。

1.2 理論計(jì)算方法

圓柱光殼體的外壓屈曲理論計(jì)算公式眾多，本文分別采用Mises 公式和Pamm 公式［12］計(jì)算殼體筒段的臨界失穩(wěn)壓力。Mises 公式為

Pamm 公式為

式中：Pcr為臨界失穩(wěn)壓力；δe為筒體厚度；D0為筒體外徑；L為筒體長度；R為筒體外徑半徑；n為筒體失穩(wěn)的波紋數(shù)，，代入?yún)?shù)求得n為2。

將殼體筒段參數(shù)代入，Mises 公式計(jì)算得筒體臨界失穩(wěn)壓力為0.087 9 MPa，Pamm 公式計(jì)算得筒體臨界失穩(wěn)壓力為0.098 3 MPa。

1.3 有限元計(jì)算方法

通過編寫ABAQUS 子程序，實(shí)現(xiàn)結(jié)構(gòu)臨界失穩(wěn)壓力和屈曲模態(tài)預(yù)測。一是線性屈曲分析，也叫做特征值分析；另一種是非線性屈曲分析，包括幾何非線性和材料非線性。本文中同時(shí)考慮了材料非線性和幾何非線性的影響。

（1） 線性屈曲分析

線性屈曲平衡方程增量形式為

式中：Ke為彈性剛度矩陣；Kσ(σ0)為初始應(yīng)力σ0下計(jì)算得到的初始應(yīng)力矩陣；ΔP為外壓載荷增量；ΔU為結(jié)構(gòu)位移增量；λ為特征值，也稱為應(yīng)力剛度矩陣的比例因子或載荷因子，結(jié)構(gòu)臨界失穩(wěn)載荷P=λ·P0，P0為施加在結(jié)構(gòu)上的初始載荷。

臨界失穩(wěn)時(shí)，ΔP=0，但ΔU≠0，因此的特征值為0，可解出其特征值以及對應(yīng)的屈曲模態(tài)。

有限元分析網(wǎng)格單元類型為四節(jié)點(diǎn)殼單元S4R，經(jīng)過網(wǎng)格獨(dú)立性驗(yàn)證，網(wǎng)格大小設(shè)置為10 mm可以滿足要求，在殼體外表面施加1 MPa 外壓，計(jì)算得殼體的特征值為0.090，則線性屈曲分析殼體臨界失穩(wěn)壓力為0.090 MPa，線性失穩(wěn)模態(tài)如圖1 所示，可知?dú)んw的臨界失穩(wěn)波數(shù)為2，與理論計(jì)算的失穩(wěn)波紋數(shù)一致。

圖1 筒體外壓作用下線性屈曲模態(tài)（放大100 倍）Fig. 1 Linear buckling mode of cylinder under external pressure (magnification 100 times)

（2） 非線性屈曲分析

線性分析方法未考慮結(jié)構(gòu)表面缺陷以及加工精度（圓度、直線度等），計(jì)算結(jié)果過于理想化；此外線性屈曲僅能計(jì)算出殼體的屈曲模態(tài)以及臨界屈曲載荷，并不能反映出殼體屈曲之后的后屈曲狀態(tài)，因此要增加非線性分析。

本文的非線性分析同時(shí)考慮材料/幾何雙非線性。材料非線性的實(shí)現(xiàn)通過改變材料性能參數(shù)，輸入材料應(yīng)力應(yīng)變曲線，殼體材料屈服強(qiáng)度為1 320 MPa，抗拉強(qiáng)度為1 620 MPa。幾何非線性的引入有2 種方法：①直接法，在筒體上直接加擾動(dòng)或設(shè)置缺陷；②間接法，將線性屈曲求得的變形乘以一個(gè)比例系數(shù)作為缺陷引入，比例系數(shù)越大，則引入的缺陷越大。在本文中采用間接法，以線性計(jì)算結(jié)果為基礎(chǔ)，通過編寫ABAQUS 關(guān)鍵字輸出線性屈曲模態(tài)，在非線性計(jì)算時(shí)通過編寫程序引入線性屈曲模態(tài)并乘以比例系數(shù)作為殼體缺陷，考慮加工精度，引入比例系數(shù)0.01。

非線性屈曲分析采用弧長法計(jì)算，監(jiān)測筒體中部一點(diǎn)的位移，分析殼體屈曲過程，加載過程中監(jiān)測點(diǎn)位移U隨載荷比例因子（load proportional factor，LPF）的變化曲線如圖2所示。

圖2 載荷比例因子-位移曲線圖Fig. 2 Load scale factor-curve

由圖2 可以看出，隨著載荷比例因子的增加，監(jiān)測點(diǎn)的位移逐漸變大；當(dāng)比例因子增加到一定量時(shí)出現(xiàn)轉(zhuǎn)折點(diǎn)，在轉(zhuǎn)折點(diǎn)之前，隨著載荷比例因子的增加，監(jiān)測點(diǎn)的位移緩慢變大，在轉(zhuǎn)折點(diǎn)之后，載荷比例因子只要增加一個(gè)微小量，監(jiān)測點(diǎn)的位移就會(huì)出現(xiàn)較大的變化。這與殼體的屈曲過程一致，當(dāng)外壓小于臨界失穩(wěn)壓力時(shí)，殼體變形量較小，一旦外壓達(dá)到臨界失穩(wěn)壓力，殼體瞬間垮塌，出現(xiàn)大變形。轉(zhuǎn)折點(diǎn)對應(yīng)的載荷即為殼體臨界失穩(wěn)載荷，通過載荷比例因子和施加載荷相乘可得臨界失穩(wěn)壓力為0.088 8 MPa。

通過線性分析方法和非線性分析方法的對比可以看出，非線性分析方法計(jì)算的臨界失穩(wěn)壓力比線性分析方法的計(jì)算結(jié)果低了1.3%。

1.4 對比分析

通過以上計(jì)算可以得出，有限元計(jì)算結(jié)果和理論公式計(jì)算結(jié)果十分接近，為避免計(jì)算結(jié)果的單一性和偶然性，增加本文計(jì)算方法的說服力，改變殼體的長徑比，通過多種長徑比的光殼體對比分析驗(yàn)證。上述殼體的長徑比為12，保持殼體的其他尺寸不變，通過改變殼體長度，研究長徑比，分別為14，10，8，6 的光殼體的臨界失穩(wěn)壓力的變化。同樣采用上述的4 種計(jì)算方法，計(jì)算結(jié)果如圖3 所示。

圖3 長徑比-臨界失穩(wěn)壓力曲線圖Fig. 3 Length-diameter ratio-critical instability pressure curve

通過圖3 中的計(jì)算結(jié)果可以看出，隨著長徑比增大，殼體的臨界失穩(wěn)壓力逐漸減小，殼體的抗屈曲性能逐漸變差；通過4 種方法的對比計(jì)算結(jié)果可以看出，當(dāng)長徑比小于8 時(shí)，4 種方法的計(jì)算結(jié)果較為分散，誤差較大，但有限元的2 種計(jì)算結(jié)果介于2種理論公式計(jì)算結(jié)果之間，相對可靠；當(dāng)長徑比不小于8 時(shí)，4 種方法的計(jì)算結(jié)果相差不大，Mises 公式計(jì)算結(jié)果和有限元線性計(jì)算結(jié)果、非線性計(jì)算結(jié)果十分接近，Mises 公式計(jì)算結(jié)果最小，線性計(jì)算結(jié)果最大，而非線性計(jì)算結(jié)果居中；以Mises 公式計(jì)算結(jié)果為基準(zhǔn)，本文所研究的殼體（長徑比為12）線性計(jì)算結(jié)果誤差為2.4%，非線性計(jì)算結(jié)果誤差為1.0%，同樣也表明了非線性計(jì)算結(jié)果更接近于殼體實(shí)際屈曲情況。

通過不同長徑比的光殼體屈曲分析計(jì)算可以看出，以Mises 公式計(jì)算結(jié)果為基準(zhǔn)，當(dāng)長徑比不小于8 時(shí)有限元計(jì)算結(jié)果的誤差都在2.7%內(nèi)，表明本文基于ABAQUS 的殼體外壓屈曲分析方法準(zhǔn)確可靠，可以用于后續(xù)分析計(jì)算。

2 帶定心部殼體的外壓屈曲分析

為了增加固體火箭發(fā)動(dòng)機(jī)殼體的剛度，通常會(huì)在殼體上增加一定數(shù)量的定心部結(jié)構(gòu)，定心部起到與加強(qiáng)圈類似的功能，對殼體的剛度有較大的提升作用。本節(jié)將采用第1 節(jié)中驗(yàn)證過的有限元屈曲分析方法分析帶定心部殼體的屈曲性能變化。

2.1 定心部軸向長度對殼體屈曲性能的影響

保持殼體其他尺寸不變，在殼體中心位置增加一個(gè)定心部，其結(jié)構(gòu)如圖4 所示。固定定心部厚度為5 mm，分別取定心部軸向長度為50，100，150，200，250 mm，非線性計(jì)算同時(shí)考慮材料/幾何雙非線性。

圖4 帶一個(gè)定心部的殼體結(jié)構(gòu)示意圖Fig. 4 Schematic diagram of shell structure with a centering part

帶不同軸向長度定心部殼體的徑向位移隨外壓的變化如圖5 所示，可以看出隨著外壓增加，殼體的徑向位移逐漸變大。當(dāng)外壓增加到一定量時(shí)，殼體徑向位移變化出現(xiàn)拐點(diǎn)，當(dāng)外壓超過拐點(diǎn)，殼體的徑向位移迅速變大。拐點(diǎn)對應(yīng)的外壓即為殼體的臨界失穩(wěn)外壓。

圖5 外壓-位移曲線圖Fig. 5 Curve of external pressure-displacement

不同軸向長度定心部殼體的臨界失穩(wěn)外壓如圖6 所示，可以看出隨著定心部軸向長度增加，殼體的臨界失穩(wěn)外壓逐漸增大，且臨界失穩(wěn)外壓呈線性增長。定心部軸向長度由50 mm 增加到250 mm，殼體的臨界失穩(wěn)外壓提高了73%，表明定心部軸向長度的變化對殼體抗失穩(wěn)性能有較大的提升作用?？梢钥闯?，考慮了材料和幾何雙非線性計(jì)算的結(jié)果低于線性計(jì)算的結(jié)果，但二者相差不大，誤差都在1.6%內(nèi)。

圖6 定心部軸向長度-臨界失穩(wěn)壓力曲線圖Fig. 6 Curve of axial length of centering part-critical instability pressure

2.2 定心部厚度對殼體屈曲性能的影響

保持殼體其他尺寸不變，殼體中心位置設(shè)置一個(gè)定心部，固定定心部軸向長度為100 mm，分別取定心部厚度為3，4，5，6，7，8，9，10 mm進(jìn)行外壓屈曲分析，計(jì)算結(jié)果如圖7所示。

圖7 定心部厚度-臨界失穩(wěn)壓力曲線圖Fig. 7 Curve of centering thickness-critical instability pressure

由圖7 可以看出，當(dāng)定心部厚度小于8 mm 時(shí)，隨著定心部厚度增加，殼體的臨界失穩(wěn)外壓逐漸變大，定心部厚度由3 mm 增加到8 mm，殼體的臨界失穩(wěn)外壓變大了1.35 倍；當(dāng)定心部厚度不小于8 mm時(shí)，殼體的臨界失穩(wěn)外壓趨于穩(wěn)定，即使定心部厚度繼續(xù)增加，殼體的臨界失穩(wěn)外壓變化也不會(huì)太大。

分析不同厚度定心部殼體在外壓作用下的線性屈曲模態(tài)，當(dāng)定心部厚度小于8 mm 時(shí)，殼體的屈曲模態(tài)類似于圖8a）的形態(tài)，殼體呈現(xiàn)出整體失穩(wěn)，失穩(wěn)波紋數(shù)為2；當(dāng)定心部厚度不小于8 mm 時(shí)，殼體的屈曲模態(tài)類似于圖8b）的形態(tài)，殼體兩側(cè)局部失穩(wěn)，局部失穩(wěn)波紋數(shù)為3。當(dāng)定心部厚度不小于8 mm 時(shí)，由于定心部剛度太大，此時(shí)相當(dāng)于將殼體以定心部為界限分為了2 個(gè)短圓筒，再增加定心部的厚度，對兩側(cè)短圓筒的屈曲性能不會(huì)產(chǎn)生太大影響，因此圖7 中當(dāng)定心部厚度不小于8 mm 時(shí)，呈現(xiàn)出殼體的臨界失穩(wěn)外壓趨于穩(wěn)定的狀態(tài)。

圖8 殼體外壓作用下線性屈曲模態(tài)（放大100 倍）Fig. 8 Linear buckling mode of shell under external pressure (100 times magnification)

2.3 定心部位置對殼體屈曲性能的影響

保持殼體其它尺寸不變，固定定心部軸向長度為100 mm，定心部厚度為5 mm。改變定心部的位置，以殼體左側(cè)端面為原點(diǎn)，分別取定心部位于550，1 100，1 650，2 200，2 750，3 300，3 850 mm 進(jìn)行外壓屈曲分析，其屈曲模態(tài)如圖9 所示。

圖9 殼體外壓作用下線性屈曲模態(tài)（放大100 倍）Fig. 9 Linear buckling mode of shell under external pressure （100 times magnification）

由圖9 可以看出，以當(dāng)前尺寸的定心部對發(fā)動(dòng)機(jī)殼體進(jìn)行補(bǔ)強(qiáng)，定心部位于殼體的任意位置，殼體的屈曲模態(tài)都為整體失穩(wěn)，且失穩(wěn)波紋數(shù)為2。定心部位于殼體的不同位置，其屈曲模態(tài)存在差異。隨著定心部位置向中間移動(dòng)，殼體的屈曲中心逐漸偏離殼體中心向另一側(cè)移動(dòng)。以殼體中心截面為對稱平面，位于殼體兩側(cè)對稱位置處的屈曲模態(tài)一致，其臨界失穩(wěn)壓力也相同。

由圖10 可以看出，當(dāng)定心部位于殼體中心時(shí)，殼體的臨界失穩(wěn)外壓最大，隨著定心部的位置逐漸向兩側(cè)移動(dòng)，殼體的臨界失穩(wěn)外壓逐漸減小。將定心部由殼體中心移到距離兩側(cè)端面550 mm 處，殼體的臨界失穩(wěn)外壓降低了27.9%。

2.4 1 個(gè)定心部和3 個(gè)定心部計(jì)算結(jié)果對比

保持殼體其他尺寸不變，固定定心部軸向長度為100 mm，定心部厚度為5 mm。以殼體左側(cè)端面為原點(diǎn)，在1 100，2 200，3 300 mm 處各設(shè)計(jì)一個(gè)定心部進(jìn)行外壓屈曲分析，其線性屈曲模態(tài)如圖11 所示，其臨界失穩(wěn)壓力計(jì)算結(jié)果如表1 所示。

表1 殼體外壓作用下臨界失穩(wěn)壓力Table 1 Critical instability pressure under external pressure of shell MPa

圖11 殼體外壓作用下線性屈曲模態(tài)（放大100 倍）Fig. 11 Linear buckling mode of shell under external pressure (100 times magnification)

由圖11 可以看出，在上述尺寸約束下，3 個(gè)定心部和1 個(gè)定心部的失穩(wěn)模態(tài)都為整體失穩(wěn)，失穩(wěn)波紋數(shù)為2。表明3 個(gè)定心部仍然不足以改變殼體的屈曲模態(tài)。通過表1 的計(jì)算結(jié)果對比可以看出，在上述定心部尺寸約束下，由中心1 個(gè)定心部增加到3個(gè)定心部均布，臨界失穩(wěn)壓力增加了36.1%。

3 結(jié)論

文中以ABAQUS 分析軟件為平臺(tái)，通過編寫子程序模擬殼體的屈曲過程。光殼體的屈曲分析采用有限元分析法和理論公式計(jì)算法對比分析，對計(jì)算方法進(jìn)行驗(yàn)證。將此方法應(yīng)用于帶定心部的固體火箭發(fā)動(dòng)機(jī)殼體的屈曲分析中，對固體火箭發(fā)動(dòng)機(jī)殼體的屈曲性能進(jìn)行分析。得到以下幾點(diǎn)結(jié)論：

（1） 考慮幾何/材料雙非線性的有限元分析方法計(jì)算的臨界失穩(wěn)壓力比線性分析方法的計(jì)算結(jié)果小。當(dāng)殼體長徑比不小于8 時(shí)有限元計(jì)算方法和Mises 公式、Pamm 公式計(jì)算的結(jié)果十分接近，計(jì)算結(jié)果適應(yīng)性較好，當(dāng)殼體長徑比小于8，Mises 公式和Pamm 公式的計(jì)算結(jié)果較為分散，有限元計(jì)算結(jié)果介于二者之間。

（2） 相較于增加定心部的軸向長度，增加定心部的厚度對殼體抗屈曲性能的提升更為顯著。定心部軸向長度由50 mm 增加到250 mm，殼體的臨界失穩(wěn)外壓提高了73%，定心部厚度由3 mm 增加到8 mm，殼體的臨界失穩(wěn)外壓變大了1.35 倍。

（3） 定心部位置越靠近殼體中間，殼體的抗屈曲性能越好。在文中定心部尺寸的約束下，由中心1 個(gè)定心部增加到3 個(gè)定心部均布，臨界失穩(wěn)壓力增加了36.1%。