固體火箭發動機含裂紋藥柱應力場有限元模擬的新方法①

任海峰，高 鳴

(海軍航空工程學院， 煙臺 264001)

?

固體火箭發動機含裂紋藥柱應力場有限元模擬的新方法①

任海峰，高 鳴

(海軍航空工程學院， 煙臺 264001)

針對研究固體火箭發動機藥柱出現裂紋前、后藥柱內應力/應變場的需要，提出利用奇異單元和生死單元技術模擬含三維裂紋藥柱的新方法，并利用該方法對固體火箭發動機三維非貫穿裂紋進行模擬，分析藥柱裂紋附近區域應力分布的規律。結果表明，該方法便捷有效，尤其適用于對比研究裂紋、脫粘等藥柱缺陷引起的應力釋放和應力分布的變化。

固體火箭發動機；藥柱；奇異單元；生死單元；裂紋

0 引言

固體火箭發動機在制造、運輸、貯存和使用等過程可能產生裂紋、脫粘等典型藥柱缺陷。發動機點火增壓時，燃氣壓力、熱量和振動等因素可能導致裂紋擴展，致使內彈道改變，甚至引發災難性事故。因此，藥柱完整性方面的研究是固體火箭發動機領域關注的重點問題之一[1]。但作為藥柱主體的推進劑是以高聚物粘合劑為基體，并在粘合劑中填加固體氧化劑、金屬燃燒劑和少量其他成分的多相混合物。材料本身的復雜性對裂紋的萌生和擴展產生的復雜影響很難描述[2-5]，一些裂紋開裂與止裂準則是隨溫度和時間變動的，試驗測定環境與實際使用環境有較大差異[6]。發動機工作過程中裂紋的擴展過程還是一個復雜的熱量-流體-結構耦合過程[7]。此外，斷裂力學的相關理論如應力強度因子斷裂準則、COD斷裂準則、最大能量釋放率準則和J積分準則等判斷裂紋是否發生擴展仍需深入研究。利用有限元方法分析裂紋的萌生和擴展涉及預設裂紋擴展路徑、重新劃分化網格等操作，建模過程紛繁復雜[7]。因此，固體火箭發動機裂紋擴展機理復雜，現有裂紋模擬方法紛繁復雜，且效果欠佳。

固體火箭發動機實際服役過程中，往往只關心藥柱是否存在裂紋，發生藥柱裂紋缺陷的發動機一般直接判廢或進行針對性的處理。近年隨著固體儀器發動機技術的發展[8-9]，業界開始嘗試使用傳感器監測裂紋。該研究方向更多地關心裂紋產生、擴展對藥柱應力場的影響，尋求粘接界面應力敏感點和缺陷特征點，并不關心裂紋萌生擴展等復雜過程。

本文以Abaqus軟件為平臺，研究利用奇異單元和生死單元技術模擬藥柱有無三維裂紋情況下藥柱應力場的新方法，并利用該方法對發動機三維非貫穿裂紋進行模擬，分析裂紋對發動機藥柱應力分布的影響。此外，該方法同樣適用于藥柱脫粘的分析，也可利用J積分準則等裂紋斷裂準則判斷裂紋是否擴展，并預測裂紋擴展方向，具有較好的應用前景。

1 裂紋奇異單元和生死單元技術

1.1 裂紋奇異單元技術

按照斷裂力學理論，裂紋尖端的應力、位移場可表示為[10-11]

(1)

(2)

式中σij(N)(i，j=1，2，3)為應力分量；ui(N)(i=1，2，3)為位移分量；N=Ⅰ～Ⅲ為裂紋類型；fij(N)、gij(N)(θ)為極角θ的函數。

圖1 1/4節點二次等參奇異單元Fig.1 1/4 node isoparametric quadratic singular element

xA、xB、xC為A、B、C節點的坐標，uA、uB、uC分別為三節點的水平位移。

裂紋線上任意一點的坐標x和位移u都可用形函數插值：

x=-0.5η(1-η)xA+(1+η)(1-η)xB+0.5η(1+η)xC

u=-0.5η(1-η)uA+(1+η)(1-η)uB+0.5η(1+η)uC

-1≤η≤1

(3)

當xA=0，xB=L/4，xC=L，uA=0時，式(3)簡化為

x=0.25(1+η)2L

u=(1+η)[(1-η)uB+0.5ηuC]

(4)

位移u對坐標進行微分：

(5)

在極坐標中x=r，θ=0，故式(5)可寫為

(6)

當xB=L/2時：

x=0.5(1+η)L

u=(1+η)[(1-η)uB+0.5ηuC]

(7)

(8)

由式(8)可知，應力、應變不具有奇異性。因此，裂紋尖端采用1/4節點等參二次奇異單元可較好地實現應力、應變場的模擬。

1.2 生死單元技術

在模擬藥柱是否存在裂紋時，本文采用了生死單元技術。有限元軟件進行分析時，單元的“生死”是通過修改單元剛度的方式實現的[13]。單元被“殺死”時，單元并非將其從實際模型中刪除，只是將其單元剛度矩陣乘上一個極小的縮減系數(該參數的缺省值為1.0×10-6，該參數取數值足夠小即可，并非常量)，為了防止矩陣奇異，該參數不設置為0。與被“殺死”單元有關的單元載荷矢量(如應力、溫度等)為0輸出。對于“殺死”單元質量、阻尼和應力剛度矩陣設置為0。因為“殺死”的單元沒有被刪除，不包括在載荷列陣中，仍然可以通過列表顯示它們。所以“殺死”前后剛度矩陣尺寸總是保持不變的，也就是說“殺死”前后模型單元一致，無須重構網絡，這對求解分析帶來極大便利。

單元的“激活”則是將被“殺死”的單元重新“激活”，在被“殺死”單元被重新“激活”時，其剛度、質量、單元載荷等均可以恢復至真實取值。但被重新“激活”的單元應變為0，若存在初應變，則可以通過實常數方式輸入，并不受單元生死操作的影響。

2 數值模擬過程

以約束邊界的推進劑方坯出現裂紋為例，利用奇異單元和生死單元技術模擬裂紋應力、應變場，涉及如下技術。

2.1 三維奇異裂紋單元建模

方坯裂紋尖端采用1/4節點等參二次奇異單元處理，如圖2所示。本算例中Abaqus三維奇異單元建模分3步：

(1)用草圖分割面，并指定掃略路徑；

(2)采用掃略的方法分割部件；

(3)采用結構化和掃略2種策略劃分模型網格，采用C3D8T單元，整體和局部網格如圖2。

圖2 奇異裂紋單元建模過程Fig.2 3D singular crack element modeling process with Abaqus

2.2 生死單元建模

生死單元建模過程如圖3所示。在軟件中建立分析步：Step→Create step→Step-2；Step-2：Dynamic，Implicit(用動態過程分析裂紋產生過程)，時間1 s，增量步0.000 1。二定義接觸：在Interaction→Create interaction，選Step-2 ；如圖3所示，選Model change→Continue，對話框，在視圖中拾取裂紋，視圖紅色陰影部分，點選Deactivated in this step。

2.3 模擬結果

仿真結果如圖4所示，裂紋尖端應力最大，截面的應力場分布模式符合I型裂紋尖端彈塑性區應力場分布特征。由于邊界為固支邊界，在粘接界面裂紋根部形成產生了應力集中，這與物理過程相符，方法簡捷有效。

圖3 生死單元建模過程Fig.3 Birth-death element modeling process with Abaqus

圖4 三維曲線裂紋最大主應力云圖Fig.4 Max principal stress of 3D curvilinear crack

3 發動機藥柱三維非貫穿裂紋模擬

某試驗用模擬發動機采用貼壁澆注的中心孔型藥柱，外徑150 mm，內徑50 mm。降溫固化引起的初始應力分布，可通過單元“生”狀態下的有限元模擬獲得其計算結果。此后，在中心孔軸向中點，人為設置長10 mm、深15 mm、寬0.5 mm的楔形裂紋，可通過單元“死”狀態下的有限元模擬獲得其計算結果。研究裂紋對藥柱和包覆層粘接表面應力的分布的影響，以便為發動機裂紋監測系統及其傳感器的設計布局優化提供參考。具體仿真過程在此不再詳述，仿真結果如圖5所示。

圖5(a)給出了模擬發動機中心孔藥三維非貫穿裂紋應力云圖全景。如圖5(b)所示，該方法很好的模擬了裂紋尖端應力場，裂紋尖端應力最大，完全符合I型裂紋尖端彈塑性區應力場分布特征，與文獻[2，14]一致。如圖5(a)、(c)所示，該方法模擬了非貫穿裂紋兩端面應力集中情況，體現了結構對應力分布的影響,比如由于材料的破壞引起的應力集中,裂紋尖端兩側端面結構對應力應變的限制，以及靠近中心孔處裂紋應力的釋放。

(a) 全景應力云圖

(b) 截面應力云圖

(c) 1/4截面應力云圖

根據有限元的近似解性質[15]，應力和應變近似解一定在精確解上下震蕩，但某些點上解正好和精確解相等，即最佳應力點。根據以應力為自變量的最小位能原理，采用高斯數值積分，由高斯積分的性質，可知積分點處應力的精度最高。因此，Abaqus中的應力存儲在高斯積分點。一般來說，裂紋仿真都要重新劃分裂紋網格，導致網格節點和單元的重新排列，需要將單元積分點應力向外插值到節點，獲得節點應力，差值方法十分復雜。采用本文方法裂紋開裂前后無需重構網格，節點和單元一致，結果對比分析簡單可靠。圖6是模擬發動機粘接界面中點圓周方向最大主應力曲線，顯示了中心孔裂紋對固體火箭發動機粘接界面應力分布的影響，對比分析不難看出裂紋開裂后，0°應力最大，40°～60°應力最小，180°應力恢復大約為未有裂紋的90%以上，與文獻[16-17]一致。圖7是裂紋產生前后粘接界面最大主應力分布相對變化著色面，著色面圖代表應裂紋發生前后粘接界面應力的差值，該圖非常直觀地反映了粘接界面應力變化情況。裂紋引起的粘接界面應力相對于無裂紋狀態的藥柱粘接界面應力呈現出規律性和方向性變化。粘接界面應力變化整體保持在0 之下，較好地體現了中心孔裂紋引起的“附加位移”，產生的應力釋放，導致藥柱的整體“變軟”，承載能力下降。裂紋尖端局部應力集中產生了應力激增，0°附近位置粘接界面應力局部激增達最大值。局部應力集中引起應力增加和整體應力釋放的較量在0°～20°之間，應力集中占據上峰，在20°～30°左右達到平衡，在30°～60°之間附件位移產生的應力釋放占據上峰，40°～60°粘接界面應力達到最小值，最大值與最小值之間差值為3 kPa。此后，應力集中和應力釋放的影響隨著遠離裂紋逐漸衰減。由于裂紋的出現引起的應力存在以上的方向特征，因此發動機藥柱缺陷的監測系統可以通過合理地選擇和布置傳感器監測裂紋的萌生并判斷其位置。對于脫粘等其他會引起藥柱內應力分布變化的缺陷，其分析方法和監測技術也同樣有效。

圖6 模擬發動機粘接界面中點圓周方向最大主應力曲線Fig.6 Max main stress along circumferential direction by a bore crack of SRM

圖7 裂紋產生前后模擬發動機粘接界面 應力分布相對變化著色圖Fig.7 Relative bond stress 3D shaded surface and contour by a bore crack of SRM

4 結論

(1)基于奇異單元和生死單元技術的含三維裂紋藥柱完整性模擬的技術，便捷有效，且可為結果對比分析提供了便利。

(2)裂紋的出現使藥柱粘接界面應力分布呈一定的方向特性，使粘接界面應力監測系統的設計和優化成為可能。

(3)裂紋引起的藥柱粘接界面應力變化幅值較小，監(檢)測有一定難度，要求監測傳感器具有很高的靈敏度和穩定性，或采取降溫等措施使應力變化幅值增大。

[1] 胡松啟，周宴星，劉迎吉,等.固體推進劑裂紋擴展研究綜述[J].火箭推進，2012，38 (5):86-92.

[2] 劉著卿,顏世東,丁彪. 藥柱裂紋對固體火箭發動機工作過程的影響[J]. 海軍航空工程學院學報,2007,22(4): 443-446.

[3] 劉承武，陽建紅，鄧凱,等.基于Mori-Tanaka 有限元法的粘彈復合推進劑非線性界面脫粘[J]. 推進技術,2011,32(2):225-229.

[4] Bencher C D． Microstructural damage and fracture processes in a composite solid rocket propellant[J]． Journal of Spacecraft and Rockets，1995，32(2):328-334.

[5] 曾甲牙．丁羥推進劑拉伸斷裂行為的掃描電鏡研究[J].固體火箭技術，1999，22(4):69-72.

[6] 彭威，周建平，趙維昌，等． 復合固體推進劑應力分布的數值模擬及損傷萌生分[J]. 固體火箭技術，2002，25(1):12-15.

[7] Dick W A,Health M T,Fiedler R A. Integrated 3-D simulation of solid propellant rockets[R]. AIAA 2001-3949.

[8] 高鳴,任海峰.固體儀器發動機健康檢測技術評述與研究[J].固體火箭技術，2013,36 (2)：278-284.

[9] 高鳴,任海峰. 固體儀器發動機粘接界面應力監測系統設計[J].火炸藥學報，2013,36 (4)：69-73.

[10] Rice J R,Rosengren G F. Plane strain deformation near crack tip in a power law hardening material [J].Journal of Mechanics and Physics of Solids,1968,16(1):1-12.

[11] 唐立強,譚英杰,蔡艷紅. 剛性-粘彈性材料界面I 型動態擴展裂紋的尖端場[J].力學季刊，2004,25(2):188-194.

[12] Homas J R Hughes,Akin J E.Techniques for developing special finite shape function with particular reference to singularities[J]. International Journal for Numerical Methods in Engineering，1980，15：733-751.

[13] ABAQUS Inc. ABAQUS Analysis User's Manual/6.12.1 ABAQUS/Aqua analysis [K].2010.

[14] 徐學文，孫建國，牟俊林. 某固體火箭發動機藥柱上三維裂紋擴展的判定[J].固體火箭技術，2008，31(4)：331-335.

[15] 王勖成. 有限單元法[M].北京：清華大學出版社，2003.

[16] Tussiwand G,et al. Application of embedded sensor technology to a full-scale experimental nozzleless rocket motor[R].AIAA 2007-5790.

[17] Reeling Brouwer G,et al. Evaluation of ageing in composite propellant grains part 2[R].AIAA 2005-3803.

(編輯：呂耀輝)

A new simulation method of 3D crack of solid rocket motor grain

REN Hai-feng，GAO Ming

(Naval Aeronautical and Astronautical University，Yantai 264001,China)

To study the variations in stress/strain field when the crack appears,a new technique to simulate 3D crack of solid rocket motor grain based on singular element and birth-death element was proposed. A part-through crack for solid rocket motor grain simulated by utilizing this method and stress distribution caused by crack were analyzed. Simulation results indicate that the method is correct,effective,and convenient for further research,especially suitable for comparative analysis on the stress distribution changes and the stress reliefs which are caused by crack or debond.

solid rocket motor；grain；singular element;birth-death element;crack

2014-11-06；

：2014-12-01。

國防預研項目(51328050101)。

任海峰(1978—)，男，博士，研究方向為火箭發動機健康監測技術。E-mail:haifeng_ren@163.com

高鳴(1957—)，男，博士，博士生導師,研究方向為火箭發動機技術。

V438

A

1006-2793(2015)01-0050-05

10.7673/j.issn.1006-2793.2015.01.09