基于ABAQUS 的固體火箭發動機殼體外壓屈曲分析 *

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

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

0 引言

隨著設計制造技術的提升以及超高強度新材料的涌現，固體火箭發動機金屬殼體的筒體壁厚越做越薄。在火箭的飛行過程中，會出現固體火箭發動機承受外壓的情況，這無疑給越做越薄的固體火箭發動機殼體帶來全新的考驗。

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

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

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

借鑒上述研究經驗，本文以ABAQUS 有限元軟件為分析平臺，從光殼體入手探索外壓屈曲的分析方法，再將此方法應用于固體火箭發動機殼體的外壓屈曲分析中，對固體火箭發動機殼體的抗屈曲性能進行了研究。

1 光殼體的外壓屈曲分析

1.1 模型概述

殼體筒段長徑比為12，壁厚為2 mm。考慮一般殼體兩端為封頭或法蘭等加強結構，約束殼體兩端的圓度。殼體材料采用D406A，其彈性模量為210 GPa，泊松比取0.3。

1.2 理論計算方法

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

Pamm 公式為

式中：Pcr為臨界失穩壓力；δe為筒體厚度；D0為筒體外徑；L為筒體長度；R為筒體外徑半徑；n為筒體失穩的波紋數，，代入參數求得n為2。

將殼體筒段參數代入，Mises 公式計算得筒體臨界失穩壓力為0.087 9 MPa，Pamm 公式計算得筒體臨界失穩壓力為0.098 3 MPa。

1.3 有限元計算方法

通過編寫ABAQUS 子程序，實現結構臨界失穩壓力和屈曲模態預測。一是線性屈曲分析，也叫做特征值分析；另一種是非線性屈曲分析，包括幾何非線性和材料非線性。本文中同時考慮了材料非線性和幾何非線性的影響。

（1） 線性屈曲分析

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

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

臨界失穩時，ΔP=0，但ΔU≠0，因此的特征值為0，可解出其特征值以及對應的屈曲模態。

有限元分析網格單元類型為四節點殼單元S4R，經過網格獨立性驗證，網格大小設置為10 mm可以滿足要求，在殼體外表面施加1 MPa 外壓，計算得殼體的特征值為0.090，則線性屈曲分析殼體臨界失穩壓力為0.090 MPa，線性失穩模態如圖1 所示，可知殼體的臨界失穩波數為2，與理論計算的失穩波紋數一致。

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

（2） 非線性屈曲分析

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

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

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

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

由圖2 可以看出，隨著載荷比例因子的增加，監測點的位移逐漸變大；當比例因子增加到一定量時出現轉折點，在轉折點之前，隨著載荷比例因子的增加，監測點的位移緩慢變大，在轉折點之后，載荷比例因子只要增加一個微小量，監測點的位移就會出現較大的變化。這與殼體的屈曲過程一致，當外壓小于臨界失穩壓力時，殼體變形量較小，一旦外壓達到臨界失穩壓力，殼體瞬間垮塌，出現大變形。轉折點對應的載荷即為殼體臨界失穩載荷，通過載荷比例因子和施加載荷相乘可得臨界失穩壓力為0.088 8 MPa。

通過線性分析方法和非線性分析方法的對比可以看出，非線性分析方法計算的臨界失穩壓力比線性分析方法的計算結果低了1.3%。

1.4 對比分析

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

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

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

通過不同長徑比的光殼體屈曲分析計算可以看出，以Mises 公式計算結果為基準，當長徑比不小于8 時有限元計算結果的誤差都在2.7%內，表明本文基于ABAQUS 的殼體外壓屈曲分析方法準確可靠，可以用于后續分析計算。

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

為了增加固體火箭發動機殼體的剛度，通常會在殼體上增加一定數量的定心部結構，定心部起到與加強圈類似的功能，對殼體的剛度有較大的提升作用。本節將采用第1 節中驗證過的有限元屈曲分析方法分析帶定心部殼體的屈曲性能變化。

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

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

圖4 帶一個定心部的殼體結構示意圖Fig. 4 Schematic diagram of shell structure with a centering part

帶不同軸向長度定心部殼體的徑向位移隨外壓的變化如圖5 所示，可以看出隨著外壓增加，殼體的徑向位移逐漸變大。當外壓增加到一定量時，殼體徑向位移變化出現拐點，當外壓超過拐點，殼體的徑向位移迅速變大。拐點對應的外壓即為殼體的臨界失穩外壓。

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

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

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

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

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

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

由圖7 可以看出，當定心部厚度小于8 mm 時，隨著定心部厚度增加，殼體的臨界失穩外壓逐漸變大，定心部厚度由3 mm 增加到8 mm，殼體的臨界失穩外壓變大了1.35 倍；當定心部厚度不小于8 mm時，殼體的臨界失穩外壓趨于穩定，即使定心部厚度繼續增加，殼體的臨界失穩外壓變化也不會太大。

分析不同厚度定心部殼體在外壓作用下的線性屈曲模態，當定心部厚度小于8 mm 時，殼體的屈曲模態類似于圖8a）的形態，殼體呈現出整體失穩，失穩波紋數為2；當定心部厚度不小于8 mm 時，殼體的屈曲模態類似于圖8b）的形態，殼體兩側局部失穩，局部失穩波紋數為3。當定心部厚度不小于8 mm 時，由于定心部剛度太大，此時相當于將殼體以定心部為界限分為了2 個短圓筒，再增加定心部的厚度，對兩側短圓筒的屈曲性能不會產生太大影響，因此圖7 中當定心部厚度不小于8 mm 時，呈現出殼體的臨界失穩外壓趨于穩定的狀態。

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

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

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

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

由圖9 可以看出，以當前尺寸的定心部對發動機殼體進行補強，定心部位于殼體的任意位置，殼體的屈曲模態都為整體失穩，且失穩波紋數為2。定心部位于殼體的不同位置，其屈曲模態存在差異。隨著定心部位置向中間移動，殼體的屈曲中心逐漸偏離殼體中心向另一側移動。以殼體中心截面為對稱平面，位于殼體兩側對稱位置處的屈曲模態一致，其臨界失穩壓力也相同。

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

2.4 1 個定心部和3 個定心部計算結果對比

保持殼體其他尺寸不變，固定定心部軸向長度為100 mm，定心部厚度為5 mm。以殼體左側端面為原點，在1 100，2 200，3 300 mm 處各設計一個定心部進行外壓屈曲分析，其線性屈曲模態如圖11 所示，其臨界失穩壓力計算結果如表1 所示。

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

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

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

3 結論

文中以ABAQUS 分析軟件為平臺，通過編寫子程序模擬殼體的屈曲過程。光殼體的屈曲分析采用有限元分析法和理論公式計算法對比分析，對計算方法進行驗證。將此方法應用于帶定心部的固體火箭發動機殼體的屈曲分析中，對固體火箭發動機殼體的屈曲性能進行分析。得到以下幾點結論：

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

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

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