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

唐國(guó)金,鄧 斌,申志彬

(國(guó)防科技大學(xué) 航天科學(xué)與工程學(xué)院,長(zhǎng)沙 410073)

0 引言

固體發(fā)動(dòng)機(jī)藥柱在長(zhǎng)期服役過(guò)程中會(huì)發(fā)生化學(xué)老化，并在長(zhǎng)期載荷作用下可能導(dǎo)致力學(xué)損傷。此外，推進(jìn)劑是典型的粘彈性材料，其泊松比為時(shí)間的函數(shù)[1]。因此，藥柱精細(xì)結(jié)構(gòu)完整性分析需考慮推進(jìn)劑變泊松比、老化及損傷等效應(yīng)的影響。

老化和損傷是導(dǎo)致推進(jìn)劑本構(gòu)非線性的兩大主要因素。關(guān)于推進(jìn)劑含老化或損傷粘彈性本構(gòu)理論，目前已有較多研究成果[2-6]。這些研究工作多采用定泊松比模型，即將泊松比視為常數(shù)，未考慮粘彈性材料泊松比的時(shí)間相關(guān)性，且這些本構(gòu)理論尚未很好地實(shí)現(xiàn)工程應(yīng)用。眾所周知，泊松比參數(shù)的選取對(duì)藥柱結(jié)構(gòu)響應(yīng)結(jié)果影響顯著。因此，為更準(zhǔn)確地進(jìn)行藥柱結(jié)構(gòu)分析，有必要進(jìn)一步開(kāi)展考慮變泊松比、老化和損傷效應(yīng)的粘彈性本構(gòu)模型的分析方法及工程應(yīng)用研究。

開(kāi)展藥柱老化損傷結(jié)構(gòu)分析需采用含老化和損傷的本構(gòu)模型，但分析過(guò)程涉及大量的非線性方程組求解運(yùn)算，往往需要采用有限元法進(jìn)行求解[7]。國(guó)內(nèi)有學(xué)者[8-9]利用有限元法對(duì)含損傷粘彈性本構(gòu)模型進(jìn)行了研究，實(shí)現(xiàn)了簡(jiǎn)單平面問(wèn)題的應(yīng)用分析，但對(duì)于三維復(fù)雜結(jié)構(gòu)的有限元分析問(wèn)題，由于涉及到復(fù)雜的前后處理、單元構(gòu)建以及大量非線性求解運(yùn)算等問(wèn)題，完全自主編程實(shí)現(xiàn)代價(jià)巨大，且也難以滿足面向工程應(yīng)用的需求。

利用Abaqus、Marc等商用有限元軟件，可對(duì)藥柱進(jìn)行線粘彈性分析，但無(wú)法考慮變泊松比、老化和損傷效應(yīng)等的影響。然而，它們提供了用戶材料模型二次開(kāi)發(fā)功能。一些學(xué)者[10-12]通過(guò)商業(yè)有限元軟件的二次開(kāi)發(fā)技術(shù)，實(shí)現(xiàn)了非線性粘彈性本構(gòu)模型的有限元計(jì)算，并取得了良好的效果。本文通過(guò)建立一種考慮變泊松比、老化和損傷效應(yīng)的三維粘彈性本構(gòu)模型，采用增量有限元法對(duì)其進(jìn)行數(shù)值離散，同時(shí)為實(shí)現(xiàn)該本構(gòu)模型有效應(yīng)用于工程實(shí)踐，基于Abaqus軟件的二次開(kāi)發(fā)技術(shù)編寫了相應(yīng)的材料子程序，并在已有的“發(fā)動(dòng)機(jī)結(jié)構(gòu)分析系統(tǒng)”[13]基礎(chǔ)上，開(kāi)發(fā)可同時(shí)考慮推進(jìn)劑變泊松比、老化和損傷效應(yīng)的藥柱結(jié)構(gòu)分析模塊。

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

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

對(duì)于丁羥推進(jìn)劑等高聚物材料，可采用交聯(lián)度近似地表征其化學(xué)老化過(guò)程。交聯(lián)度υ的發(fā)展方程一般可采用如下形式[3]

(1)

其中

式中t′為老化時(shí)間；T為老化溫度；Ea為老化反應(yīng)的活化能；k為Boltzmann常數(shù)；h為Planck常數(shù)；α為化學(xué)配分函數(shù)待定系數(shù)；A為化學(xué)反應(yīng)速率常數(shù)；υ0為初始交聯(lián)度；υm為最大交聯(lián)度。

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

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

(2)

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

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

g(x)=0 (x>1)

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

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

其中

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

根據(jù)由松弛模量E和泊松比ν表示的線粘彈性本構(gòu)模型[1]，考慮推進(jìn)劑老化和損傷的影響[3，14]，并假設(shè)老化時(shí)間與加載時(shí)間起點(diǎn)相同(取式(1)中t′=t)，綜合可得含老化和損傷效應(yīng)的推進(jìn)劑三維熱粘彈性本構(gòu)方程為

(3)

(4)

(5)

溫度平移因子aT滿足WLF方程

(6)

式中C1和C2為材料常數(shù)；T為當(dāng)前溫度，T0為參考溫度。

式(4)中的泊松比ν和含老化松弛模量E的形式分別為

(7)

(8)

對(duì)于式(3)所示的本構(gòu)方程，如果考慮材料在加載前已老化了時(shí)間td，則根據(jù)式(1)可知，只需將式(1)中的系數(shù)B等效替換成Be-β(td,T)，此時(shí)的本構(gòu)方程(3)即可考慮加載之前的老化效應(yīng)。

2 數(shù)值計(jì)算模型

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

下面采用增量法對(duì)本構(gòu)模型進(jìn)行數(shù)值離散。將分析時(shí)間[0,t]劃分為[ti-1,ti](i=1,2,…,M)共M個(gè)子區(qū)間。在各增量步內(nèi)，假設(shè)有效應(yīng)力和應(yīng)變均隨折算時(shí)間ξ′呈線性變化，采用積分算法對(duì)式(3)～(5)進(jìn)行數(shù)值離散，并結(jié)合有效應(yīng)力形式，最終可得到在時(shí)間步[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計(jì)算所采用的Newton法具有二階收斂速率。根據(jù)切線剛度的定義，對(duì)于小變形或者小體變的大變形問(wèn)題，其一致切線剛度的形式為[7]

(24)

綜合式(9)和(24)，可得到tm+1時(shí)刻的一致切線剛度，其各個(gè)分量為

式中

對(duì)于各向同性材料，具有如下的關(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)力更新方法

針對(duì)增量型本構(gòu)方程(9)，將其進(jìn)一步改寫成

(25)

令

則方程組(25)可簡(jiǎn)寫成

σij(tm+1)+Dij(tm+1){ωm+1[σij(tm+1)]-1}=0

(26)

式(26)是以應(yīng)力為變量的非線性方程組，對(duì)此一般需采用迭代法進(jìn)行求解。本文擬采用Newton法[7]進(jìn)行求解，進(jìn)而可實(shí)現(xiàn)當(dāng)前應(yīng)力等的更新。

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

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

為實(shí)現(xiàn)新型材料本構(gòu)模型的Abaqus二次開(kāi)發(fā)，用戶需通過(guò)其UMAT(User-defined Mechanical Material Behavior)接口開(kāi)發(fā)相應(yīng)的材料子程序。針對(duì)上文給出的考慮變泊松比的含老化和損傷粘彈性本構(gòu)數(shù)值模型，開(kāi)發(fā)了相應(yīng)的材料子程序，其基本實(shí)現(xiàn)流程如圖1所示。

實(shí)現(xiàn)的主要思路：根據(jù)主程序(Abaqus)提供的初始應(yīng)力、應(yīng)變、應(yīng)變?cè)隽俊囟纫约安牧蠀?shù)等相關(guān)信息，材料子程序先后完成含老化和損傷的非線性方程組數(shù)值求解，實(shí)現(xiàn)當(dāng)前時(shí)刻的應(yīng)力和狀態(tài)變量等的更新，最后向主程序返回當(dāng)前應(yīng)力、本構(gòu)的一致切線剛度以及狀態(tài)變量等變量。

圖1 材料子程序?qū)崿F(xiàn)流程圖

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

在“發(fā)動(dòng)機(jī)結(jié)構(gòu)分析軟件系統(tǒng)”[13]的基礎(chǔ)上，進(jìn)一步嵌入本文含老化和損傷的藥柱結(jié)構(gòu)分析模塊。該模塊可實(shí)現(xiàn)用戶材料參數(shù)的輸入、材料子程序與求解器的連接、計(jì)算提交及結(jié)果后處理等功能，其分析過(guò)程關(guān)鍵的圖形化交互操作界面如圖2～圖4所示。

圖2 藥柱結(jié)構(gòu)分析模塊主菜單

為滿足不同平臺(tái)下的分析需要，開(kāi)發(fā)了Abaqus和Marc 2種不同環(huán)境下的求解模塊(圖3)。用戶可根據(jù)需要選擇不同的本構(gòu)模型和求解器來(lái)進(jìn)行計(jì)算。通過(guò)將開(kāi)發(fā)的材料子程序以及分析模塊集成到綜合軟件分析系統(tǒng)中，可充分利用商用有限元軟件的強(qiáng)大前后處理和求解能力，方便地實(shí)現(xiàn)實(shí)際型號(hào)固體火箭發(fā)動(dòng)機(jī)的工程應(yīng)用分析。

(a)結(jié)構(gòu)分析主面板 (b)用戶材料參數(shù)設(shè)置

圖4 松弛模量參數(shù)輸入界面

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

下面利用上述所得藥柱結(jié)構(gòu)分析模塊應(yīng)用于實(shí)際型號(hào)發(fā)動(dòng)機(jī)藥柱結(jié)構(gòu)分析。

以某固體發(fā)動(dòng)機(jī)為研究對(duì)象，根據(jù)其載荷和結(jié)構(gòu)對(duì)稱性，取其1/12結(jié)構(gòu)建立圖5所示有限元模型。式(8)所示的含老化松弛模量各系數(shù)如表 1所示，式(7)所示的變泊松比取5項(xiàng)時(shí)的系數(shù)如表 2所示，其他材料參數(shù)見(jiàn)表 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 級(jí)數(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的地面點(diǎn)火增壓過(guò)程。其中，固化降溫過(guò)程為從零應(yīng)力的T1=58 ℃自然降溫到T0=20 ℃，溫度變化近似服從規(guī)律[15]

T(t)=TI-(TI-T0)(1-e-6.8×10-5t)

(27)

在點(diǎn)火增壓過(guò)程，對(duì)其內(nèi)表面施加均布?jí)毫d荷：p(t)=6(1-e-20t)MPa。

計(jì)算時(shí)約束模型的環(huán)向位移，并先后采用表 4 所示的8種不同效應(yīng)的本構(gòu)模型，計(jì)算上述過(guò)程藥柱的結(jié)構(gòu)響應(yīng)。表 4中，標(biāo)“√”為含有該效應(yīng)；標(biāo)“×”為不含有該效應(yīng)。以模型8的計(jì)算結(jié)果為例，給出點(diǎn)火增壓結(jié)束時(shí)的藥柱Von Mises應(yīng)力(MPa)和損傷分布圖，分別如圖6和圖7所示。

表4 采用的本構(gòu)模型

圖6 藥柱Von Mises應(yīng)力分布圖

圖7 藥柱損傷分布圖

由圖6和圖7可知，藥柱Von Mises應(yīng)力在內(nèi)表面的過(guò)渡段和圓管段部位較大，而損傷主要發(fā)生在Von Mises應(yīng)力較大的區(qū)域，且?guī)缀跫性谶^(guò)渡段，這與式(2)中的損傷為Von Mises應(yīng)力的單調(diào)非遞減函數(shù)是相符的。以藥柱內(nèi)表面結(jié)果為例，取圓管段P1和過(guò)渡段P2兩個(gè)關(guān)鍵部位為例(見(jiàn)圖6)，其在增壓過(guò)程第0.5 s時(shí)的Von Mises應(yīng)力、Von Mises應(yīng)變及損傷對(duì)比結(jié)果見(jiàn)表 5。根據(jù)表 5，采用定泊松比模型的計(jì)算結(jié)果中，(以P2處的結(jié)果為例)，在同等條件下考慮老化效應(yīng)模型比不考慮老化效應(yīng)模型所得Von Mises應(yīng)力要高出約14%，而Von Mises應(yīng)變卻減少了約34%，這主要是由于老化效應(yīng)導(dǎo)致推進(jìn)劑模量上升引起的。在同等條件下，考慮損傷比不考慮損傷時(shí)得到的Von Mises應(yīng)力有所減小，而Von Mises應(yīng)變卻有所增加。這是因?yàn)閾p傷效應(yīng)體現(xiàn)了推進(jìn)劑的“軟化”效應(yīng)，使其有效模量減小引起的。對(duì)于采用變泊松模型所得相應(yīng)結(jié)果也有上述類似結(jié)論。

表5 不同模型下的計(jì)算結(jié)果對(duì)比

對(duì)比采用變泊松比和定泊松比模型所得結(jié)果差異發(fā)現(xiàn)，在同等情況下，變泊松比模型下對(duì)應(yīng)的應(yīng)力、應(yīng)變結(jié)果均相對(duì)較大。這主要是由于采用傳統(tǒng)的定泊松比模型分析時(shí)，其泊松比一般取值在較大的0.490～0.499之間，而本文變泊松比的取值在0.482～0.494之間(0.5 s內(nèi))，取值相對(duì)較小，但其應(yīng)力、應(yīng)變結(jié)果均卻相對(duì)較大，所得結(jié)果滿足泊松比越小、其應(yīng)力和應(yīng)變結(jié)果越大的規(guī)律。

綜上可知，推進(jìn)劑的力學(xué)性能受老化和損傷效應(yīng)等的影響顯著，老化效應(yīng)可導(dǎo)致材料“硬“化，而損傷效應(yīng)導(dǎo)致材料“軟”化。因此，相比于傳統(tǒng)線粘彈性本構(gòu)模型，本文含老化和損傷的本構(gòu)模型可更準(zhǔn)確地體現(xiàn)推進(jìn)劑的實(shí)際力學(xué)性能變化。

5 結(jié)論

(1)相比于傳統(tǒng)線粘彈性本構(gòu)模型，本文的粘彈性本構(gòu)模型體現(xiàn)了粘彈性泊松比的時(shí)間相關(guān)性，同時(shí)考慮了老化和損傷效應(yīng)，適合于發(fā)動(dòng)機(jī)長(zhǎng)期服役過(guò)程中不同載荷下的藥柱結(jié)構(gòu)分析，具有實(shí)用價(jià)值。

(2)采用有限元法對(duì)該新型本構(gòu)模型進(jìn)行了數(shù)值分析，并基于商用有限元軟件的二次開(kāi)發(fā)技術(shù)，創(chuàng)建了材料子程序，實(shí)現(xiàn)了復(fù)雜粘彈性本構(gòu)模型的三維有限元分析。研究方法對(duì)粘彈性本構(gòu)一類問(wèn)題的分析具有通用性。

(3)將藥柱結(jié)構(gòu)分析模塊集成到統(tǒng)一的發(fā)動(dòng)機(jī)綜合分析系統(tǒng)中，結(jié)合實(shí)際載荷工況，可實(shí)現(xiàn)具體型號(hào)發(fā)動(dòng)機(jī)結(jié)構(gòu)分析。研究工作為進(jìn)一步分析固體發(fā)動(dòng)機(jī)藥柱精細(xì)結(jié)構(gòu)完整性與壽命評(píng)估提供了技術(shù)支撐，具有工程應(yīng)用價(jià)值。

參考文獻(xiàn)：

[1] 趙伯華. 固體推進(jìn)劑粘彈泊松比的研究[J].北京理工大學(xué)學(xué)報(bào),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] 周建平. 化學(xué)不穩(wěn)定的工程材料的粘彈性損傷本構(gòu)模型 [D].長(zhǎng)沙: 國(guó)防科技大學(xué),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] 彭威,周建平,任鈞國(guó).復(fù)合固體推進(jìn)劑非線性粘彈本構(gòu)方程的微觀力學(xué)分析[J].固體火箭技術(shù),1999,22(4):23-25.

[6] 陽(yáng)建紅,俞茂宏,候根良,等.HTPB復(fù)合固體推進(jìn)劑含損傷和老化本構(gòu)研究[J].推進(jìn)技術(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].上海力學(xué),1995,16(1):20-26.

[9] 杜建科,朱祖念,張善祁,等.固體火箭發(fā)動(dòng)機(jī)藥柱損傷粘彈性有限元分析[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二次開(kāi)發(fā)的固體發(fā)動(dòng)機(jī)結(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ā)動(dòng)機(jī)藥形改進(jìn)與優(yōu)化設(shè)計(jì)[D].長(zhǎng)沙:國(guó)防科技大學(xué),2011.