固體火箭發(fā)動(dòng)機(jī)含裂紋藥柱應(yīng)力場有限元模擬的新方法①

任海峰，高 鳴

(海軍航空工程學(xué)院， 煙臺(tái) 264001)

?

固體火箭發(fā)動(dòng)機(jī)含裂紋藥柱應(yīng)力場有限元模擬的新方法①

任海峰，高 鳴

(海軍航空工程學(xué)院， 煙臺(tái) 264001)

針對(duì)研究固體火箭發(fā)動(dòng)機(jī)藥柱出現(xiàn)裂紋前、后藥柱內(nèi)應(yīng)力/應(yīng)變場的需要，提出利用奇異單元和生死單元技術(shù)模擬含三維裂紋藥柱的新方法，并利用該方法對(duì)固體火箭發(fā)動(dòng)機(jī)三維非貫穿裂紋進(jìn)行模擬，分析藥柱裂紋附近區(qū)域應(yīng)力分布的規(guī)律。結(jié)果表明，該方法便捷有效，尤其適用于對(duì)比研究裂紋、脫粘等藥柱缺陷引起的應(yīng)力釋放和應(yīng)力分布的變化。

固體火箭發(fā)動(dòng)機(jī)；藥柱；奇異單元；生死單元；裂紋

0 引言

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

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

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

1 裂紋奇異單元和生死單元技術(shù)

1.1 裂紋奇異單元技術(shù)

按照斷裂力學(xué)理論，裂紋尖端的應(yīng)力、位移場可表示為[10-11]

(1)

(2)

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

圖1 1/4節(jié)點(diǎn)二次等參奇異單元Fig.1 1/4 node isoparametric quadratic singular element

xA、xB、xC為A、B、C節(jié)點(diǎn)的坐標(biāo)，uA、uB、uC分別為三節(jié)點(diǎn)的水平位移。

裂紋線上任意一點(diǎn)的坐標(biāo)x和位移u都可用形函數(shù)插值：

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)

當(dāng)xA=0，xB=L/4，xC=L，uA=0時(shí)，式(3)簡化為

x=0.25(1+η)2L

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

(4)

位移u對(duì)坐標(biāo)進(jìn)行微分：

(5)

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

(6)

當(dāng)xB=L/2時(shí)：

x=0.5(1+η)L

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

(7)

(8)

由式(8)可知，應(yīng)力、應(yīng)變不具有奇異性。因此，裂紋尖端采用1/4節(jié)點(diǎn)等參二次奇異單元可較好地實(shí)現(xiàn)應(yīng)力、應(yīng)變場的模擬。

1.2 生死單元技術(shù)

在模擬藥柱是否存在裂紋時(shí)，本文采用了生死單元技術(shù)。有限元軟件進(jìn)行分析時(shí)，單元的“生死”是通過修改單元?jiǎng)偠鹊姆绞綄?shí)現(xiàn)的[13]。單元被“殺死”時(shí)，單元并非將其從實(shí)際模型中刪除，只是將其單元?jiǎng)偠染仃嚦松弦粋€(gè)極小的縮減系數(shù)(該參數(shù)的缺省值為1.0×10-6，該參數(shù)取數(shù)值足夠小即可，并非常量)，為了防止矩陣奇異，該參數(shù)不設(shè)置為0。與被“殺死”單元有關(guān)的單元載荷矢量(如應(yīng)力、溫度等)為0輸出。對(duì)于“殺死”單元質(zhì)量、阻尼和應(yīng)力剛度矩陣設(shè)置為0。因?yàn)椤皻⑺馈钡膯卧獩]有被刪除，不包括在載荷列陣中，仍然可以通過列表顯示它們。所以“殺死”前后剛度矩陣尺寸總是保持不變的，也就是說“殺死”前后模型單元一致，無須重構(gòu)網(wǎng)絡(luò)，這對(duì)求解分析帶來極大便利。

單元的“激活”則是將被“殺死”的單元重新“激活”，在被“殺死”單元被重新“激活”時(shí)，其剛度、質(zhì)量、單元載荷等均可以恢復(fù)至真實(shí)取值。但被重新“激活”的單元應(yīng)變?yōu)?，若存在初應(yīng)變，則可以通過實(shí)常數(shù)方式輸入，并不受單元生死操作的影響。

2 數(shù)值模擬過程

以約束邊界的推進(jìn)劑方坯出現(xiàn)裂紋為例，利用奇異單元和生死單元技術(shù)模擬裂紋應(yīng)力、應(yīng)變場，涉及如下技術(shù)。

2.1 三維奇異裂紋單元建模

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

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

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

(3)采用結(jié)構(gòu)化和掃略2種策略劃分模型網(wǎng)格，采用C3D8T單元，整體和局部網(wǎng)格如圖2。

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

2.2 生死單元建模

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

2.3 模擬結(jié)果

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

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

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

3 發(fā)動(dòng)機(jī)藥柱三維非貫穿裂紋模擬

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

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

(a) 全景應(yīng)力云圖

(b) 截面應(yīng)力云圖

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

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

圖6 模擬發(fā)動(dòng)機(jī)粘接界面中點(diǎn)圓周方向最大主應(yīng)力曲線Fig.6 Max main stress along circumferential direction by a bore crack of SRM

圖7 裂紋產(chǎn)生前后模擬發(fā)動(dòng)機(jī)粘接界面 應(yīng)力分布相對(duì)變化著色圖Fig.7 Relative bond stress 3D shaded surface and contour by a bore crack of SRM

4 結(jié)論

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

(2)裂紋的出現(xiàn)使藥柱粘接界面應(yīng)力分布呈一定的方向特性，使粘接界面應(yīng)力監(jiān)測系統(tǒng)的設(shè)計(jì)和優(yōu)化成為可能。

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

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

[2] 劉著卿,顏世東,丁彪. 藥柱裂紋對(duì)固體火箭發(fā)動(dòng)機(jī)工作過程的影響[J]. 海軍航空工程學(xué)院學(xué)報(bào),2007,22(4): 443-446.

[3] 劉承武，陽建紅，鄧凱,等.基于Mori-Tanaka 有限元法的粘彈復(fù)合推進(jìn)劑非線性界面脫粘[J]. 推進(jìn)技術(shù),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ìn)劑拉伸斷裂行為的掃描電鏡研究[J].固體火箭技術(shù)，1999，22(4):69-72.

[6] 彭威，周建平，趙維昌，等． 復(fù)合固體推進(jìn)劑應(yīng)力分布的數(shù)值模擬及損傷萌生分[J]. 固體火箭技術(shù)，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] 高鳴,任海峰.固體儀器發(fā)動(dòng)機(jī)健康檢測技術(shù)評(píng)述與研究[J].固體火箭技術(shù)，2013,36 (2)：278-284.

[9] 高鳴,任海峰. 固體儀器發(fā)動(dòng)機(jī)粘接界面應(yīng)力監(jiān)測系統(tǒng)設(shè)計(jì)[J].火炸藥學(xué)報(bào)，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] 唐立強(qiáng),譚英杰,蔡艷紅. 剛性-粘彈性材料界面I 型動(dòng)態(tài)擴(kuò)展裂紋的尖端場[J].力學(xué)季刊，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] 徐學(xué)文，孫建國，牟俊林. 某固體火箭發(fā)動(dòng)機(jī)藥柱上三維裂紋擴(kuò)展的判定[J].固體火箭技術(shù)，2008，31(4)：331-335.

[15] 王勖成. 有限單元法[M].北京：清華大學(xué)出版社，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。

國防預(yù)研項(xiàng)目(51328050101)。

任海峰(1978—)，男，博士，研究方向?yàn)榛鸺l(fā)動(dòng)機(jī)健康監(jiān)測技術(shù)。E-mail:haifeng_ren@163.com

高鳴(1957—)，男，博士，博士生導(dǎo)師,研究方向?yàn)榛鸺l(fā)動(dòng)機(jī)技術(shù)。

V438

A

1006-2793(2015)01-0050-05

10.7673/j.issn.1006-2793.2015.01.09