孟紅磊，鞠玉濤

(1.中國航天科技集團(tuán)公司四院四十一所，西安 710025;2.南京理工大學(xué)機(jī)械工程學(xué)院，南京 210094)

0 引言

固體火箭發(fā)動機(jī)的裝藥結(jié)構(gòu)完整性問題一直是制約固體火箭發(fā)動機(jī)發(fā)展的關(guān)鍵因素，為了精確地進(jìn)行裝藥結(jié)構(gòu)完整性數(shù)值仿真，需要建立精確的本構(gòu)方程和數(shù)值仿真方法。固體推進(jìn)劑是典型的粘彈性材料，在一定的載荷作用下，表現(xiàn)出非線性粘彈性特性，現(xiàn)有的裝藥結(jié)構(gòu)完整性分析更多地采用線性粘彈性本構(gòu)方程［1－7］，而線性粘彈性本構(gòu)方程與實(shí)際誤差較大。近年來，非線性粘彈性本構(gòu)方程發(fā)展較快，應(yīng)用較廣的有Schapery非線性粘彈性本構(gòu)［8－9］、Leaderman 非線性粘彈性本構(gòu)［10］、微分型非線性粘彈性本構(gòu)［11］及含損傷的非線性粘彈性本構(gòu)方程［12－14］。其中，含損傷的非線性本構(gòu)方程認(rèn)為非線性是由于材料的內(nèi)部損傷引起的，能較好地描述材料的非線性粘彈性特性，在瀝青混凝土、固體推進(jìn)劑等方面應(yīng)用較多。陽建紅［12］針對復(fù)合固體推進(jìn)劑建立了含損傷的非線性粘彈性本構(gòu)方程，彭威［13］針對復(fù)合固體推進(jìn)劑建立了含細(xì)觀損傷的非線性粘彈性本構(gòu)方程，但都未針對本構(gòu)方程進(jìn)行二次開發(fā)和數(shù)值計(jì)算應(yīng)用;Hinterhoelzl R M和Schapery R A［14］針對瀝青混凝土材料，建立了含損傷的非線性粘彈性本構(gòu)方程，并進(jìn)行二次開發(fā)，但其本構(gòu)形式復(fù)雜，二次開發(fā)難度大，導(dǎo)致計(jì)算量大。

本文針對固體推進(jìn)劑材料，提出了一種新形式的含累積損傷的非線性粘彈性本構(gòu)方程，并擴(kuò)展到三維形式，形式簡單，便于二次開發(fā)，且能較好反映材料的非線性粘彈性本構(gòu)方程。對裝藥結(jié)構(gòu)完整性數(shù)值仿真提供了精確可行的本構(gòu)方程及數(shù)值方法。

1 本構(gòu)方程

1.1 彈性粘彈性對應(yīng)原理

Schapery最初針對彈性材料發(fā)展起來的損傷理論，可通過彈性-粘彈性對應(yīng)原理擴(kuò)展到粘彈性材料中，通過引入偽應(yīng)變的概念，將彈性本構(gòu)中的應(yīng)變值用偽應(yīng)變替代后，可變?yōu)檎硰椥员緲?gòu)方程，偽應(yīng)變定義為

式中ER為參考彈性模量，與彈性模量單位一致，引入?yún)⒖紡椥阅Ａ?，主要是為了能使偽?yīng)變同樣是無量綱的，ER可任意選取，通常選取1值或選用材料的初始彈性模量值作為ER。

偽應(yīng)變是一個卷積分的形式，是一個包含時間遺傳因素的量，將偽應(yīng)變值帶入線性粘彈性本構(gòu)方程后，積分形式的線性粘彈性本構(gòu)方程變?yōu)?/p>

式(2)表示的粘彈性本構(gòu)方程有著與彈性本構(gòu)方程相似的形式。對于線性粘彈性來說，ER是一個常數(shù)，不隨著加載過程應(yīng)變水平的變化而變化。

1.2 含累積損傷的非線性本構(gòu)方程

當(dāng)應(yīng)變值過大時，應(yīng)力-偽應(yīng)變曲線出現(xiàn)明顯的非線性特性，且不同的加載速率和不同的加載路徑時，應(yīng)力-偽應(yīng)變曲線并不重合，這說明ER值不僅和加載水平有關(guān)，還和加載過程有關(guān)。因此，為了能準(zhǔn)確描述固體推進(jìn)劑在較大變形下的力學(xué)響應(yīng)特性，必須建立起ER關(guān)于加載過程和加載水平的函數(shù)形式?？赏ㄟ^引入損傷因子來準(zhǔn)確描述非線性粘彈性，如式(3)所示:

式中 S為損傷內(nèi)變量，是表征材料內(nèi)部損傷程度的物理量;C(S)稱為軟化函數(shù)，是關(guān)于損傷內(nèi)變量S的函數(shù);C并不是常數(shù)，是一個跟加載過程相關(guān)的變量，C值是關(guān)于損傷內(nèi)變量S的函數(shù)。

本文結(jié)合累積損傷的概念，建立含累積損傷的非線性粘彈性本構(gòu)方程，并證明含累積損傷的非線性粘彈性本構(gòu)方程能較好地反映不同應(yīng)變率條件下及復(fù)雜應(yīng)變率條件下的應(yīng)力應(yīng)變響應(yīng)，含累積損傷的非線性粘彈性本構(gòu)方程為

由于改性雙基推進(jìn)劑力學(xué)特性不具有明顯的方向性，而且在三向應(yīng)力狀態(tài)下，材料在未出現(xiàn)明顯裂紋之前損傷不具有明顯的方向性，在工程中也經(jīng)常將推進(jìn)劑力學(xué)特性視為各向同性。所以，本文近似認(rèn)為三維應(yīng)力狀態(tài)下，損傷是各向同性的，即損傷內(nèi)變量S是標(biāo)量形式，損傷變量S用單一標(biāo)量表示，那么軟化函數(shù)C也為標(biāo)量。

則含損傷的三維粘彈性本構(gòu)方程表示為

Lai J和hai-Ali Rami M建立的三維非線性粘彈性本構(gòu)方程中，進(jìn)行了各向同性假設(shè)，其中的非線性系數(shù)假設(shè)為八面體剪應(yīng)力的函數(shù)［15］。本文基于該思想，假設(shè)三維應(yīng)力條件下，損傷內(nèi)變量是關(guān)于Mises等效應(yīng)力的函數(shù)形式。在三向應(yīng)力狀態(tài)下，Mises等效應(yīng)力常用來判斷材料包括固體推進(jìn)劑材料的屈服，Mises等效應(yīng)力的水平可衡量材料的損傷軟化程度。所以，三維應(yīng)力條件下，損傷內(nèi)變量S演化方程中的應(yīng)力值用Mises等效應(yīng)力代替。

式中 σeq為Mises等效應(yīng)力。

當(dāng)單軸拉伸時，σ1=σx，σ2=σ3=0。其中，σx指拉伸方向應(yīng)力，等效應(yīng)力可退化為單軸應(yīng)力形式:

1.3 本構(gòu)方程參數(shù)擬合

含累積損傷的非線性粘彈性本構(gòu)中，需擬合的有損傷參數(shù)λ值和η值，及軟化函數(shù)C(S)的形式。擬合過程需先通過松弛試驗(yàn)獲得松弛模量，然后根據(jù)5組等速拉伸曲線，前面損傷變量λ值和η值是先賦初值，λ值可任意選擇，一般選為1，η值初值的選擇可先選擇利用累積損傷模型預(yù)測材料損傷時得到的值，具體見文獻(xiàn)［16］;在上述給定的損傷變量的情況下，得到的不同速率下的C-S曲線不一定能很好地重合;經(jīng)過分析，λ值只會影響C(S)函數(shù)的擬合結(jié)果，不會影響不同速率下C-S曲線的重合性，而η值的選取會影響不同速率下C-S曲線的重合性;所以，λ值就取初始值為最終值，η值需不斷調(diào)整以獲取最優(yōu)值;根據(jù)最終的5組重合較好的C-S曲線，擬合出C(S)的函數(shù)形式。根據(jù)曲線的形式選擇合適的函數(shù)形式描述C-S關(guān)系，如C(S)=1－aSb的形式。

2 數(shù)值仿真應(yīng)用

2.1 本構(gòu)方程增量形式及切線模量求解

將通過對含累積損傷的非線性粘彈性本構(gòu)方程數(shù)值展開，來獲得應(yīng)力更新表達(dá)式和切線模量矩陣(Jacobian矩陣)。

由主應(yīng)力 σ1、σ2、σ3求解 Von Mises等效應(yīng)力，并計(jì)算損傷內(nèi)變量值S。

偽應(yīng)變由偏偽應(yīng)變部分和體積偽應(yīng)變部分組成

則應(yīng)力張量分量形式等于

對偏偽應(yīng)變進(jìn)行展開，可得到

將求和符號中的各項(xiàng)積分分別單獨(dú)表示。其中，qij，n(t)表示偏偽應(yīng)變中prony級數(shù)中的第n項(xiàng)的遺傳積分，qij，∞(t)則表示偏偽應(yīng)變的穩(wěn)態(tài)項(xiàng)。

將式(13)和式(14)代入式(12)中，偏偽應(yīng)變可表示為各項(xiàng)遺傳積分之和的形式。

由式(15)可得到t+Δt時刻偏偽應(yīng)變增量:

同理，體積偽應(yīng)變展開可得到

將體積偽應(yīng)變中求和符號內(nèi)的各個積分項(xiàng)也分別單獨(dú)表示。其中，qkk，n(t)表示體積偽應(yīng)變中prony級數(shù)中的第n項(xiàng)的遺傳積分，qkk，∞(t)則表示偏偽應(yīng)變的穩(wěn)態(tài)項(xiàng)。

將式(18)和式(19)代入式(17)中，體積偽應(yīng)變可表示為各項(xiàng)遺傳積分之和的形式:

由式(20)得到t+Δt時刻體積偽應(yīng)變增量:

對損傷內(nèi)變量S進(jìn)行離散，可得

對偏偽應(yīng)變中prony級數(shù)的第n項(xiàng)的遺傳積分進(jìn)行數(shù)值離散，將界限(t，t+Δt)內(nèi)的積分表示為兩段積分，第一段是極限(0，t)內(nèi)積分，第二段是(t，t+Δt)內(nèi)積分，并將t～t+Δt時間段內(nèi)的應(yīng)變率近似為恒值，用平均應(yīng)變率代替，整理后可得到:

由式(24)得到

由式(13)得到

同理，對體積偽應(yīng)變中prony級數(shù)的第n項(xiàng)的遺傳積分進(jìn)行數(shù)值離散，整理后可得到:

由式(18)得

由式(27)得

將式(25)和式(26)代入式(16)中，得到偏偽應(yīng)變增量:

將式(28)和式(29)代入式(21)中，得到體積偽應(yīng)變增量:

總的偽應(yīng)變增量可表示為偏偽應(yīng)變增量和體積偽應(yīng)變增量之和的形式:

進(jìn)而可得到

對于粘彈性材料來說，一致切線模量分量形式表示為［14］

可得出三維時切線模量表達(dá)式［14］:

式(36)中認(rèn)為，粘彈性和損傷對Jacobian矩陣的影響是解耦的。由于軟化函數(shù)C(S)和損傷內(nèi)變量是關(guān)于Mises等效應(yīng)力的復(fù)雜的函數(shù)形式，因此精確的一致切線模量的計(jì)算，將是相當(dāng)復(fù)雜并耗費(fèi)相當(dāng)多的計(jì)算時間，Hinterhoelz.l R M 和 Schapery R A［14］在對其建立的三維含損傷非線性粘彈性本構(gòu)方程求解一致切線模量時，同樣遇到這樣的問題，他們在對應(yīng)力關(guān)于偽應(yīng)變的求偏導(dǎo)時，認(rèn)為損傷內(nèi)變量是固定的，即可得到式(36)的結(jié)論。由此獲得的切線模量并不和應(yīng)力的更新相一致，相當(dāng)于割線模量，這樣會導(dǎo)致Newton-Raphson迭代過程收斂變慢，但在Hinterhoelzl R M和Schapery R A及本文進(jìn)行的實(shí)際計(jì)算中，并沒有出現(xiàn)收斂問題，每個增量步一般通過2～4次迭代過程即可收斂。所以，式(36)既能順利實(shí)現(xiàn)Newton-Raphson迭代收斂，又便于應(yīng)用。Simo和Taylor［17］在1985年針對率相關(guān)材料提出了一致切線算子(CTO)的概念，一致切線算子為應(yīng)力關(guān)于應(yīng)變增量的偏導(dǎo)。Hinterhoelzl R M基于其一致切線算子獲得切線模量。

將式(30)和式(31)代入到式(32)中，得到

令

將式(39)和式(40)代入式(38)中，可得

由式(41)可得

結(jié)合式(37)和式(42)可得

2.2 UMAT子程序介紹

UMAT在配合ABAQUS求解非線性問題時的主要任務(wù):

(1)提供Jacobian矩陣，以便組裝切線剛度矩陣，用于迭代求解平衡方程。其中，Jacobian矩陣不影響計(jì)算結(jié)果的準(zhǔn)確性，只影響收斂速度。

(2)提供應(yīng)力增量和應(yīng)變增量之間的關(guān)系，用于更新應(yīng)力。

首先在每個增量步開始時，程序給節(jié)點(diǎn)外力一個增量，平衡方程不再平衡，ABAQUS會根據(jù)上一增量步提供的切線剛度矩陣計(jì)算出節(jié)點(diǎn)位移變量，然后形成新的剛度矩陣并進(jìn)行迭代計(jì)算，每一次迭代步都需調(diào)用UMAT提供的Jacobian矩陣來計(jì)算節(jié)點(diǎn)位移變量，進(jìn)而得到節(jié)點(diǎn)應(yīng)變增量，然后再調(diào)用UMAT，進(jìn)而得到應(yīng)力增量，進(jìn)行應(yīng)力更新，積分得到新的內(nèi)力，如果節(jié)點(diǎn)內(nèi)力和節(jié)點(diǎn)外力仍未平衡，繼續(xù)迭代，直至收斂，獲得收斂后的應(yīng)變增量和應(yīng)力增量，再通過差值函數(shù)得到單元應(yīng)變和單元應(yīng)力。至此，該增量步已經(jīng)完成，進(jìn)入下一個增量步進(jìn)行求解計(jì)算，UMAT中定義的狀態(tài)變量也將隨之傳遞給下一個增量步。

2.3 數(shù)值仿真及驗(yàn)證

基于改性雙基推進(jìn)劑，獲取本構(gòu)方程參數(shù)。利用獲取的應(yīng)力增量表達(dá)式和Jacabian矩陣，利用FORTRAN語言編制UMAT子程序，將UMAT子程序代入到ABAQUS有限元軟件中，數(shù)值模擬等速拉伸過程推進(jìn)劑材料應(yīng)力-應(yīng)變關(guān)系，結(jié)果如圖1所示。理論預(yù)測值能較好地預(yù)測材料的應(yīng)力-應(yīng)變關(guān)系。所以，本文建立的含損傷的非線性粘彈性本構(gòu)方程能較好地用于裝藥結(jié)構(gòu)完整性分析。

圖1 等速拉伸過程應(yīng)力-應(yīng)變理論預(yù)測與實(shí)驗(yàn)對比Fig.1 Comparison between prediction and experimental measurement at constant strain rate tension

3 結(jié)論

提出了一種含累積損傷的非線性粘彈性本構(gòu)方程，并獲得其增量形式和切線模量矩陣，結(jié)合ABAQUS有限元軟件，利用FORTRAN語言編制了UMAT子程序，將本構(gòu)方程應(yīng)用到數(shù)值仿真中去，并結(jié)合試驗(yàn)和數(shù)值仿真結(jié)果進(jìn)行對比，驗(yàn)證了該本構(gòu)方程和數(shù)值仿真方法的準(zhǔn)確性，為精確的裝藥結(jié)構(gòu)完整性數(shù)值分析提供了本構(gòu)和方法。本文建立的本構(gòu)方程沒有考慮溫度的影響，可在本文建立的本構(gòu)基礎(chǔ)上，通過時溫等效引入溫度，這樣可建立更為精確的本構(gòu)，其數(shù)值離散方法與本文方法近似。

［1］于洋，王寧飛，張平.一種自由裝填式組合藥柱的低溫三維結(jié)構(gòu)完整性分析［J］.固體火箭技術(shù)，2007，30(1):34-38.

［2］鐘濤，張為華.點(diǎn)火瞬態(tài)過程對復(fù)合推進(jìn)劑力學(xué)響應(yīng)特性的影響［J］.國防科技大學(xué)學(xué)報，2004，26(3):2004.

［3］于勝春，趙汝巖.固體火箭發(fā)動機(jī)快速升壓過程的流固耦合分析［J］.固體火箭技術(shù)，2008，31(3):232-235.

［4］隋欣，魏志軍.炮射導(dǎo)彈發(fā)射過程發(fā)動機(jī)裝藥強(qiáng)度分析［J］.彈道學(xué)報，2009，21(2):19-22.

［5］Fiedler R，Namazifard A，Campbell M，et al.Detailed simulation of propellant slumping in the TitanIV SRMU PQM-1［R］.AIAA 2006-4592.

［6］Chyuan S W.Dynamic analysis of solid propellant grains subjected to ignition pressurization loading［J］.Journal of Sound and Vibration，2003，268(3):.465-483.

［7］Chyuan S W.Studies of poisson＇s ratio variation for solid propellant grains under ignition pressure loading［J］.Pressure Vessels and Piping，2003，80:871-877.

［8］Schapery R A.On the characterization of nonlinear viscoelastic materials［J］.Polymer Engineering ＆ Science，1969，9:295-310.

［9］Huang Chien-wei，Eyad Masad，Anastasia H，et al.Nonlinearly viscoelastic analysis of asphalt mixes subjected to shear loading［J］.Mech.Time-Depend Mater.，2007，11:91-110.

［10］Leaderman H.Large longitudinal retarded elastic deformation of rubberlike network polymers［J］.Journal of Polymer Science，1962:361-382.

［11］Himanshu Shekhar，Sahasrabudhe A D.Viscoelastic modeling of solid rocket propellants using maxwell fluid moedel［J］.Defence Science Journal，2010，60:423-427.

［12］陽建紅，俞茂宏，侯根良，等.HTPB復(fù)合推進(jìn)劑含損傷和老化本構(gòu)研究［J］.推進(jìn)技術(shù)，2002，23(6):509-512.

［13］彭威，鄭堅(jiān)，白鴻柏，等.復(fù)合推進(jìn)劑微裂紋損傷本構(gòu)模型研究［J］.固體火箭技術(shù)，2003，26(2):33-37.

［14］Hinterhoelzl R M，Schapery R A.FEM implementation of a three-demensional viscoelastic constitutive model for particulate composites with damage growth ［J］.Mechanics of Time-Dependent Materials，2004，8:65-94.

［15］Haj-Ali R M，Muliana A H.Numerical finite element formulation of the schapery non-linear viscoelastic material model［J］.International Journal for Numerical Methods in Engineering，2004，59(1):25-45.

［16］孟紅磊，趙秀超，鞠玉濤，等.基于累積損傷的雙基推進(jìn)劑強(qiáng)度準(zhǔn)則及實(shí)驗(yàn)［J］.推進(jìn)技術(shù)，2011，32(1):109-112.

［17］Simo J C，Taylor R L.Consistent tangent operators for rateindependent elastoplasticity［J］.Computer Methods in Applied Mechanics and Engineering，1985，48(1):101-118.