唐國金,鄧 斌,申志彬
(國防科技大學 航天科學與工程學院,長沙 410073)
固體發(fā)動機藥柱在長期服役過程中會發(fā)生化學老化,并在長期載荷作用下可能導致力學損傷。此外,推進劑是典型的粘彈性材料,其泊松比為時間的函數(shù)[1]。因此,藥柱精細結(jié)構(gòu)完整性分析需考慮推進劑變泊松比、老化及損傷等效應(yīng)的影響。
老化和損傷是導致推進劑本構(gòu)非線性的兩大主要因素。關(guān)于推進劑含老化或損傷粘彈性本構(gòu)理論,目前已有較多研究成果[2-6]。這些研究工作多采用定泊松比模型,即將泊松比視為常數(shù),未考慮粘彈性材料泊松比的時間相關(guān)性,且這些本構(gòu)理論尚未很好地實現(xiàn)工程應(yīng)用。眾所周知,泊松比參數(shù)的選取對藥柱結(jié)構(gòu)響應(yīng)結(jié)果影響顯著。因此,為更準確地進行藥柱結(jié)構(gòu)分析,有必要進一步開展考慮變泊松比、老化和損傷效應(yīng)的粘彈性本構(gòu)模型的分析方法及工程應(yīng)用研究。
開展藥柱老化損傷結(jié)構(gòu)分析需采用含老化和損傷的本構(gòu)模型,但分析過程涉及大量的非線性方程組求解運算,往往需要采用有限元法進行求解[7]。國內(nèi)有學者[8-9]利用有限元法對含損傷粘彈性本構(gòu)模型進行了研究,實現(xiàn)了簡單平面問題的應(yīng)用分析,但對于三維復雜結(jié)構(gòu)的有限元分析問題,由于涉及到復雜的前后處理、單元構(gòu)建以及大量非線性求解運算等問題,完全自主編程實現(xiàn)代價巨大,且也難以滿足面向工程應(yīng)用的需求。
利用Abaqus、Marc等商用有限元軟件,可對藥柱進行線粘彈性分析,但無法考慮變泊松比、老化和損傷效應(yīng)等的影響。然而,它們提供了用戶材料模型二次開發(fā)功能。一些學者[10-12]通過商業(yè)有限元軟件的二次開發(fā)技術(shù),實現(xiàn)了非線性粘彈性本構(gòu)模型的有限元計算,并取得了良好的效果。本文通過建立一種考慮變泊松比、老化和損傷效應(yīng)的三維粘彈性本構(gòu)模型,采用增量有限元法對其進行數(shù)值離散,同時為實現(xiàn)該本構(gòu)模型有效應(yīng)用于工程實踐,基于Abaqus軟件的二次開發(fā)技術(shù)編寫了相應(yīng)的材料子程序,并在已有的“發(fā)動機結(jié)構(gòu)分析系統(tǒng)”[13]基礎(chǔ)上,開發(fā)可同時考慮推進劑變泊松比、老化和損傷效應(yīng)的藥柱結(jié)構(gòu)分析模塊。
對于丁羥推進劑等高聚物材料,可采用交聯(lián)度近似地表征其化學老化過程。交聯(lián)度υ的發(fā)展方程一般可采用如下形式[3]
(1)
其中
式中t′為老化時間;T為老化溫度;Ea為老化反應(yīng)的活化能;k為Boltzmann常數(shù);h為Planck常數(shù);α為化學配分函數(shù)待定系數(shù);A為化學反應(yīng)速率常數(shù);υ0為初始交聯(lián)度;υm為最大交聯(lián)度。
用ω表示推進劑損傷,其發(fā)展方程可取為[3]
(2)
式中σ*為當量應(yīng)力。
g(x)函數(shù)具有如下形式
g(x)=0 (x>1)
式中γ、K、β和σth為與損傷相關(guān)的材料常數(shù)。
在損傷各向同性假設(shè)下,當量應(yīng)力σ*可采用如下形式[14]
其中
根據(jù)由松弛模量E和泊松比ν表示的線粘彈性本構(gòu)模型[1],考慮推進劑老化和損傷的影響[3,14],并假設(shè)老化時間與加載時間起點相同(取式(1)中t′=t),綜合可得含老化和損傷效應(yīng)的推進劑三維熱粘彈性本構(gòu)方程為
(3)
(4)
(5)
溫度平移因子aT滿足WLF方程
(6)
式中C1和C2為材料常數(shù);T為當前溫度,T0為參考溫度。
式(4)中的泊松比ν和含老化松弛模量E的形式分別為
(7)
(8)
對于式(3)所示的本構(gòu)方程,如果考慮材料在加載前已老化了時間td,則根據(jù)式(1)可知,只需將式(1)中的系數(shù)B等效替換成Be-β(td,T),此時的本構(gòu)方程(3)即可考慮加載之前的老化效應(yīng)。
下面采用增量法對本構(gòu)模型進行數(shù)值離散。將分析時間[0,t]劃分為[ti-1,ti](i=1,2,…,M)共M個子區(qū)間。在各增量步內(nèi),假設(shè)有效應(yīng)力和應(yīng)變均隨折算時間ξ′呈線性變化,采用積分算法對式(3)~(5)進行數(shù)值離散,并結(jié)合有效應(yīng)力形式,最終可得到在時間步[tm,tm+1]內(nèi)的本構(gòu)增量形式為
(9)
(10)
(11)
(12)
式中
(13)
(14)
(15)
(16)
(17)
式中
(18)
其中
(19)
(20)
(21)
(22)
(23)
利用一致切線剛度(Jacobian),可保證Abaqus計算所采用的Newton法具有二階收斂速率。根據(jù)切線剛度的定義,對于小變形或者小體變的大變形問題,其一致切線剛度的形式為[7]
(24)
綜合式(9)和(24),可得到tm+1時刻的一致切線剛度,其各個分量為
式中
對于各向同性材料,具有如下的關(guān)系式
C2222(tm+1)=C3333(tm+1)=C1111(tm+1)
C2233(tm+1)=C1133(tm+1)=C1122(tm+1)
C2323(tm+1)=C1313(tm+1)=C1212(tm+1)
針對增量型本構(gòu)方程(9),將其進一步改寫成
(25)
令
則方程組(25)可簡寫成
σij(tm+1)+Dij(tm+1){ωm+1[σij(tm+1)]-1}=0
(26)
式(26)是以應(yīng)力為變量的非線性方程組,對此一般需采用迭代法進行求解。本文擬采用Newton法[7]進行求解,進而可實現(xiàn)當前應(yīng)力等的更新。
為實現(xiàn)新型材料本構(gòu)模型的Abaqus二次開發(fā),用戶需通過其UMAT(User-defined Mechanical Material Behavior)接口開發(fā)相應(yīng)的材料子程序。針對上文給出的考慮變泊松比的含老化和損傷粘彈性本構(gòu)數(shù)值模型,開發(fā)了相應(yīng)的材料子程序,其基本實現(xiàn)流程如圖1所示。
實現(xiàn)的主要思路:根據(jù)主程序(Abaqus)提供的初始應(yīng)力、應(yīng)變、應(yīng)變增量、溫度以及材料參數(shù)等相關(guān)信息,材料子程序先后完成含老化和損傷的非線性方程組數(shù)值求解,實現(xiàn)當前時刻的應(yīng)力和狀態(tài)變量等的更新,最后向主程序返回當前應(yīng)力、本構(gòu)的一致切線剛度以及狀態(tài)變量等變量。
圖1 材料子程序?qū)崿F(xiàn)流程圖
在“發(fā)動機結(jié)構(gòu)分析軟件系統(tǒng)”[13]的基礎(chǔ)上,進一步嵌入本文含老化和損傷的藥柱結(jié)構(gòu)分析模塊。該模塊可實現(xiàn)用戶材料參數(shù)的輸入、材料子程序與求解器的連接、計算提交及結(jié)果后處理等功能,其分析過程關(guān)鍵的圖形化交互操作界面如圖2~圖4所示。
圖2 藥柱結(jié)構(gòu)分析模塊主菜單
為滿足不同平臺下的分析需要,開發(fā)了Abaqus和Marc 2種不同環(huán)境下的求解模塊(圖3)。用戶可根據(jù)需要選擇不同的本構(gòu)模型和求解器來進行計算。通過將開發(fā)的材料子程序以及分析模塊集成到綜合軟件分析系統(tǒng)中,可充分利用商用有限元軟件的強大前后處理和求解能力,方便地實現(xiàn)實際型號固體火箭發(fā)動機的工程應(yīng)用分析。
(a)結(jié)構(gòu)分析主面板 (b)用戶材料參數(shù)設(shè)置
圖4 松弛模量參數(shù)輸入界面
下面利用上述所得藥柱結(jié)構(gòu)分析模塊應(yīng)用于實際型號發(fā)動機藥柱結(jié)構(gòu)分析。
以某固體發(fā)動機為研究對象,根據(jù)其載荷和結(jié)構(gòu)對稱性,取其1/12結(jié)構(gòu)建立圖5所示有限元模型。式(8)所示的含老化松弛模量各系數(shù)如表 1所示,式(7)所示的變泊松比取5項時的系數(shù)如表 2所示,其他材料參數(shù)見表 3。
圖5 藥柱有限元模型
n0123τEn/s—5.295 5252.955 2529.552E0n/MPa5.933 40.956 90.697 20.509 7E1n/MPa7.328 50.857 11.599 31.200 5
表2 泊松比v(t)的 Prony 級數(shù)系數(shù)
表3 其他材料參數(shù)
取式(1)中的交聯(lián)度發(fā)展方程系數(shù)A=6.95×10-12,α=-8 710,υ0=3.13×10-5mol/ml,υm=1.012×10-4mol/ml、Ea=7.438×10-20J;式(6)中的WLF方程參數(shù)C1=4.97,C2=156.1,T0=293.15 K;損傷參數(shù)σth=1.82 MPa(受壓),γ=1.046,K=0.1,β=0.5。
載荷工況:先后經(jīng)歷1 d的固化降溫、10 a的20 ℃恒溫貯存,以及0.5 s的地面點火增壓過程。其中,固化降溫過程為從零應(yīng)力的T1=58 ℃自然降溫到T0=20 ℃,溫度變化近似服從規(guī)律[15]
T(t)=TI-(TI-T0)(1-e-6.8×10-5t)
(27)
在點火增壓過程,對其內(nèi)表面施加均布壓力載荷:p(t)=6(1-e-20t)MPa。
計算時約束模型的環(huán)向位移,并先后采用表 4 所示的8種不同效應(yīng)的本構(gòu)模型,計算上述過程藥柱的結(jié)構(gòu)響應(yīng)。表 4中,標“√”為含有該效應(yīng);標“×”為不含有該效應(yīng)。以模型8的計算結(jié)果為例,給出點火增壓結(jié)束時的藥柱Von Mises應(yīng)力(MPa)和損傷分布圖,分別如圖6和圖7所示。
表4 采用的本構(gòu)模型
圖6 藥柱Von Mises應(yīng)力分布圖
圖7 藥柱損傷分布圖
由圖6和圖7可知,藥柱Von Mises應(yīng)力在內(nèi)表面的過渡段和圓管段部位較大,而損傷主要發(fā)生在Von Mises應(yīng)力較大的區(qū)域,且?guī)缀跫性谶^渡段,這與式(2)中的損傷為Von Mises應(yīng)力的單調(diào)非遞減函數(shù)是相符的。以藥柱內(nèi)表面結(jié)果為例,取圓管段P1和過渡段P2兩個關(guān)鍵部位為例(見圖6),其在增壓過程第0.5 s時的Von Mises應(yīng)力、Von Mises應(yīng)變及損傷對比結(jié)果見表 5。根據(jù)表 5,采用定泊松比模型的計算結(jié)果中,(以P2處的結(jié)果為例),在同等條件下考慮老化效應(yīng)模型比不考慮老化效應(yīng)模型所得Von Mises應(yīng)力要高出約14%,而Von Mises應(yīng)變卻減少了約34%,這主要是由于老化效應(yīng)導致推進劑模量上升引起的。在同等條件下,考慮損傷比不考慮損傷時得到的Von Mises應(yīng)力有所減小,而Von Mises應(yīng)變卻有所增加。這是因為損傷效應(yīng)體現(xiàn)了推進劑的“軟化”效應(yīng),使其有效模量減小引起的。對于采用變泊松模型所得相應(yīng)結(jié)果也有上述類似結(jié)論。
表5 不同模型下的計算結(jié)果對比
對比采用變泊松比和定泊松比模型所得結(jié)果差異發(fā)現(xiàn),在同等情況下,變泊松比模型下對應(yīng)的應(yīng)力、應(yīng)變結(jié)果均相對較大。這主要是由于采用傳統(tǒng)的定泊松比模型分析時,其泊松比一般取值在較大的0.490~0.499之間,而本文變泊松比的取值在0.482~0.494之間(0.5 s內(nèi)),取值相對較小,但其應(yīng)力、應(yīng)變結(jié)果均卻相對較大,所得結(jié)果滿足泊松比越小、其應(yīng)力和應(yīng)變結(jié)果越大的規(guī)律。
綜上可知,推進劑的力學性能受老化和損傷效應(yīng)等的影響顯著,老化效應(yīng)可導致材料“硬“化,而損傷效應(yīng)導致材料“軟”化。因此,相比于傳統(tǒng)線粘彈性本構(gòu)模型,本文含老化和損傷的本構(gòu)模型可更準確地體現(xiàn)推進劑的實際力學性能變化。
(1)相比于傳統(tǒng)線粘彈性本構(gòu)模型,本文的粘彈性本構(gòu)模型體現(xiàn)了粘彈性泊松比的時間相關(guān)性,同時考慮了老化和損傷效應(yīng),適合于發(fā)動機長期服役過程中不同載荷下的藥柱結(jié)構(gòu)分析,具有實用價值。
(2)采用有限元法對該新型本構(gòu)模型進行了數(shù)值分析,并基于商用有限元軟件的二次開發(fā)技術(shù),創(chuàng)建了材料子程序,實現(xiàn)了復雜粘彈性本構(gòu)模型的三維有限元分析。研究方法對粘彈性本構(gòu)一類問題的分析具有通用性。
(3)將藥柱結(jié)構(gòu)分析模塊集成到統(tǒng)一的發(fā)動機綜合分析系統(tǒng)中,結(jié)合實際載荷工況,可實現(xiàn)具體型號發(fā)動機結(jié)構(gòu)分析。研究工作為進一步分析固體發(fā)動機藥柱精細結(jié)構(gòu)完整性與壽命評估提供了技術(shù)支撐,具有工程應(yīng)用價值。
參考文獻:
[1] 趙伯華. 固體推進劑粘彈泊松比的研究[J].北京理工大學學報,1994,14(1):87-90.
[2] Schapery R A.A micromechanical model for non-linear viscoelastic behavior of particle-reinforced rubber with distributed damage[J].Engineering Fracture Mechanics,1986,25(5):845-867.
[3] 周建平. 化學不穩(wěn)定的工程材料的粘彈性損傷本構(gòu)模型 [D].長沙: 國防科技大學,1989.
[4] Zhou J P.A constitutive model of polymer materials including chemical ageing and mechanical damage and its experimental verification[J].Polymer,1993,34(20):4252-4256.
[5] 彭威,周建平,任鈞國.復合固體推進劑非線性粘彈本構(gòu)方程的微觀力學分析[J].固體火箭技術(shù),1999,22(4):23-25.
[6] 陽建紅,俞茂宏,候根良,等.HTPB復合固體推進劑含損傷和老化本構(gòu)研究[J].推進技術(shù),2002,23(6):509-512.
[7] Zienkiewicz O C,Taylor R L. The finite element method for solid and structural mechanics (Sixth edition)[M].Oxford: Elsevier Butterworth-Heinemann publications,2005.
[8] 馮志剛,周建平.粘彈性有限元法研究[J].上海力學,1995,16(1):20-26.
[9] 杜建科,朱祖念,張善祁,等.固體火箭發(fā)動機藥柱損傷粘彈性有限元分析[J].固體火箭技術(shù),2001,24(1):1-6
[10] Hinterhoelzl R M,Schapery R A. FEM implementation of a three-dimensional viscoelastic constitutive model for particulate composites with damage growth[J].Mechanics of time-dependent materials,2004,8(1):65-94.
[11] Masad E,Huang C W,Airey G,et al.Nonlinear viscoelastic analysis of unaged and aged asphalt binders[J].Construction and Building Materials,2008,22:2170-2179.
[12] 孟紅磊,鞠玉濤.含損傷非線性粘彈性本構(gòu)模型及其數(shù)值仿真應(yīng)用[J].固體火箭技術(shù),2012,35(6):764-772.
[13] 申志彬,李磊,段靜波,等.基于Patran二次開發(fā)的固體發(fā)動機結(jié)構(gòu)分析系統(tǒng)[J].固體火箭技術(shù),2011,34(2):171-175.
[14] Lemaitre J. A course on damage mechanics[M].Berlin:Springer,1992.
[15] 李磊. 基于結(jié)構(gòu)完整性分析的固體火箭發(fā)動機藥形改進與優(yōu)化設(shè)計[D].長沙:國防科技大學,2011.