基于推進劑復雜本構(gòu)模型的藥柱結(jié)構(gòu)分析模塊開發(fā)①

唐國金,鄧 斌,申志彬

(國防科技大學 航天科學與工程學院,長沙 410073)

0 引言

固體發(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)分析模塊。

1 一種新型粘彈性本構(gòu)模型

1.1 推進劑老化發(fā)展方程

對于丁羥推進劑等高聚物材料，可采用交聯(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)度。

1.2 推進劑損傷發(fā)展方程

用ω表示推進劑損傷，其發(fā)展方程可取為[3]

(2)

式中σ*為當量應(yīng)力。

g(x)函數(shù)具有如下形式

g(x)=0 (x>1)

式中γ、K、β和σth為與損傷相關(guān)的材料常數(shù)。

在損傷各向同性假設(shè)下，當量應(yīng)力σ*可采用如下形式[14]

其中

1.3 含老化和損傷的熱粘彈性本構(gòu)方程

根據(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)。

2 數(shù)值計算模型

2.1 本構(gòu)方程的增量形式

下面采用增量法對本構(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)

2.2 本構(gòu)方程的切線剛度

利用一致切線剛度(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)

2.3 應(yīng)力更新方法

針對增量型本構(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)力等的更新。

3 新型本構(gòu)模型的二次開發(fā)

3.1 Abaqus材料子程序開發(fā)

為實現(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)流程圖

3.2 藥柱結(jié)構(gòu)分析模塊開發(fā)

在“發(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ù)輸入界面

4 分析模塊在藥柱結(jié)構(gòu)分析中的應(yīng)用

下面利用上述所得藥柱結(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)推進劑的實際力學性能變化。

5 結(jié)論

(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.