非線性粘彈性本構(gòu)模型在固體火箭發(fā)動(dòng)機(jī)裝藥結(jié)構(gòu)完整性分析中應(yīng)用

李 毅,田維平,董新剛,姚 東

(中國航天科技集團(tuán)公司四院四十一所，西安 710025)

0 引言

復(fù)合固體推進(jìn)劑是一種顆粒增強(qiáng)的高聚型含能材料，廣泛用于運(yùn)載火箭、導(dǎo)彈等裝備的動(dòng)力系統(tǒng)，其性能直接影響裝備的運(yùn)載能力、射程、貯存性等關(guān)鍵性能。復(fù)合固體推進(jìn)劑力學(xué)特性的粘彈屬性與其自身的配方特點(diǎn)有關(guān)，也與服役環(huán)境如溫度、加載應(yīng)變率、應(yīng)力應(yīng)變狀態(tài)等因素高度相關(guān)。在變形較小時(shí)，可認(rèn)為其本構(gòu)模型為線粘彈性；隨著變形量增大，推進(jìn)劑內(nèi)部微裂紋、微孔洞等損傷逐漸演化，出現(xiàn)明顯的脫濕現(xiàn)象，此時(shí)線性理論不再適用，應(yīng)選取合適的非線性粘彈性本構(gòu)模型。

早期基于復(fù)合推進(jìn)劑的線粘彈性本構(gòu)模型開展了工程領(lǐng)域的結(jié)構(gòu)分析，結(jié)合產(chǎn)品研制的經(jīng)驗(yàn)修正取得了良好的應(yīng)用效果[1]。Scahpery[2]及其合作者一直致力于含損傷非線性粘彈性本構(gòu)模型的研究。Ulrich Mandel等[3]對纖維增強(qiáng)復(fù)合材料進(jìn)行研究，提出一種三維的非線性粘彈性本構(gòu)模型，并指出此模型可適用于任意載荷情況下。Muliana[4]從細(xì)觀角度入手，提出一種非線性粘彈性本構(gòu)模型。張曉等[5]從應(yīng)變能出發(fā)推導(dǎo)出HTPB推進(jìn)劑非線性粘彈性本構(gòu)模型。常新龍等[6]研究了HTPB高應(yīng)變率下的非線性粘彈性本構(gòu)模型。王哲君等[7]對HTPB壓應(yīng)力狀態(tài)下的本構(gòu)模型進(jìn)行了研究。姚東等[8]從損傷效應(yīng)的擬合出發(fā)，提出了考慮應(yīng)力狀態(tài)的非線性本構(gòu)模型。唐志平、朱兆祥等[9]提出了被后來的學(xué)者稱為朱-王-唐非線性本構(gòu)模型的成果，均為復(fù)合固體推進(jìn)劑非線性粘彈性本構(gòu)模型的研究作出了貢獻(xiàn)。

本文基于朱-王-唐非線性本構(gòu)模型，針對其積分項(xiàng)繁多、本構(gòu)參數(shù)獲取難度較大的特點(diǎn)，提出一種形式簡單、便于工程應(yīng)用的本構(gòu)模型。

1 推進(jìn)劑本構(gòu)模型的提出

朱-王-唐非線性本構(gòu)模型如式(1)所示：

σ=E0ε+αε2+βε3+

(1)

式中E0、E1、E2、α、β、θ1、θ2均為材料常數(shù)，前3項(xiàng)表示非線性彈性部分，后2項(xiàng)積分項(xiàng)表示粘性部分。

式(1)表明，非線性粘彈性本構(gòu)模型可寫成彈性部分和粘性部分之和的形式。參考該結(jié)論，本文將推進(jìn)劑本構(gòu)模型表述為如下形式：

σ(t,ε)=fE(t,ε)+fv(t,ε)

(2)

式中fE(t,ε)為彈性項(xiàng)；fv(t,ε)為粘性項(xiàng)。

Staverman和Schwarzl[10]在大量實(shí)驗(yàn)觀測基礎(chǔ)上指出，采用式(3)：

σ(t)=E(t)f(ε)

(3)

代替σ(t)=E(t)ε可較好地描述實(shí)驗(yàn)結(jié)果。趙榮國[11]由此假設(shè)過去不同時(shí)刻的應(yīng)變對現(xiàn)時(shí)應(yīng)力的影響可忽略不計(jì)，并對其進(jìn)行了驗(yàn)證。聯(lián)想到應(yīng)力松弛實(shí)驗(yàn)，本文認(rèn)為推進(jìn)劑的受力分為應(yīng)變變化引起的應(yīng)力及松弛現(xiàn)象導(dǎo)致的應(yīng)力耗散。對于應(yīng)變從0～ε的變化，都可認(rèn)為是以上2種應(yīng)力疊加產(chǎn)生的效果，這樣應(yīng)力的表達(dá)形式就可簡單地寫成：

σ(t,ε)=fE(ε)-σ松弛

(4)

式中fE(ε)、σ松弛分別為彈性、粘性應(yīng)力項(xiàng)，其函數(shù)形式及其中的材料參數(shù)還有待通過實(shí)驗(yàn)數(shù)據(jù)進(jìn)行確定。

2 基于松弛及拉伸實(shí)驗(yàn)的本構(gòu)模型建模

推進(jìn)劑在變形較小時(shí)，服從線粘彈性規(guī)律，當(dāng)變形達(dá)到一定值后，損傷的逐漸演化導(dǎo)致其非線性顯著增加，本文假設(shè)應(yīng)力耗散項(xiàng)不受損傷變量的影響。因此，式(4)中的σ松弛就可用式(5)表示：

σ松弛=[E(0)-E(t)]ε

(5)

式中E(0)為初始時(shí)刻的松弛模量值；E(t)為t時(shí)刻松弛模量值；ε為t時(shí)刻應(yīng)變。

粘性應(yīng)力耗散項(xiàng)是狀態(tài)量，僅與當(dāng)前狀態(tài)有關(guān)，如應(yīng)變、時(shí)間等，不受以往應(yīng)力-應(yīng)變狀態(tài)的影響。

彈性項(xiàng)的應(yīng)力-應(yīng)變關(guān)系包括大應(yīng)變狀態(tài)的數(shù)據(jù)，因此，還需進(jìn)行一定的拉伸實(shí)驗(yàn)。

2.1 松弛實(shí)驗(yàn)及其分析

推進(jìn)劑松弛實(shí)驗(yàn)采用啞鈴模型，如圖1所示。

試驗(yàn)件拉伸到應(yīng)變?yōu)?%，隨后測其應(yīng)力值，并將數(shù)據(jù)轉(zhuǎn)化為松弛模量。為得到完備的結(jié)果，需進(jìn)行多組溫度下的松弛實(shí)驗(yàn)。

根據(jù)時(shí)溫等效原理E(t,T)=E(t/αT,T0)，對不同溫度下的松弛模量進(jìn)行處理，使所有數(shù)據(jù)基本重合為一條曲線，即松弛模量主曲線。隨后，將松弛模量主曲線擬合為6階Prony級數(shù)形式：

(6)

式中各參數(shù)的值如表1所示。在對不同溫度下松弛模量進(jìn)行處理時(shí)，得到不同溫度對應(yīng)的時(shí)溫等效因子的值，時(shí)溫等效因子lgαT僅和溫度有關(guān)。根據(jù)WLF方程(7)，擬合得到的C1和C2的值分別為8.719 0、144.695 8。擬合后的松弛模量主曲線與實(shí)驗(yàn)數(shù)據(jù)的對比如圖2所示。

(7)

2.2 拉伸實(shí)驗(yàn)及其分析

拉伸試件同圖1所示模型，實(shí)驗(yàn)在Instron 4210型試驗(yàn)機(jī)上完成，測試溫度 25 ℃，應(yīng)變率 0.023 8 s-1，測得數(shù)據(jù)結(jié)果如圖3所示。

式(4)中的應(yīng)力彈性項(xiàng)不包括粘性應(yīng)力。因此，根據(jù)拉伸實(shí)驗(yàn)得到的松弛模量及式(5)，對拉伸數(shù)據(jù)進(jìn)行修正，得到彈性應(yīng)力項(xiàng)的表達(dá)式：

(8)

fE(ε)=σ實(shí)驗(yàn)+σ松弛

修正后的應(yīng)力用式(9)形式的函數(shù)進(jìn)行擬合效果如圖4所示。其中，E0和Es均為材料參數(shù)，不同溫度下的E0和Es值如表2所示。

(9)

將E0和Es表示為溫度T的函數(shù)關(guān)系式，擬合的函數(shù)關(guān)系式如式(10)、式(11)所示，擬合曲線如圖5所示。

(10)

(11)

3 本構(gòu)模型及其參數(shù)的數(shù)值驗(yàn)證

3.1 本構(gòu)模型一維MATLAB驗(yàn)證

本文采用MATLAB編寫代碼對本構(gòu)模型進(jìn)行一維數(shù)值驗(yàn)證，MATLAB中僅模擬邊界條件的變化，無需定義單元類型。在應(yīng)變率0.023 8 s-1、25 ℃的情況下，計(jì)算出的應(yīng)力-應(yīng)變走勢與實(shí)驗(yàn)數(shù)據(jù)的對比如圖6所示。可看出，應(yīng)力隨應(yīng)變變化的趨勢基本和實(shí)驗(yàn)數(shù)據(jù)吻合，脫濕點(diǎn)對應(yīng)的應(yīng)變約為20%，脫濕后的誤差明顯大于線性段，相對偏差達(dá)8%。本文假設(shè)應(yīng)力耗散與損傷無關(guān)，損傷的增加有可能會導(dǎo)致粘性應(yīng)力耗散加大，具體規(guī)律需進(jìn)一步研究。

粘彈性材料的力學(xué)性能具有溫度相關(guān)性和載荷率相關(guān)性，需進(jìn)行多組不同邊界條件的數(shù)值仿真，現(xiàn)將不同條件下的數(shù)值仿真結(jié)果同實(shí)驗(yàn)的對比匯總?cè)鐖D7所示。從圖7可看出，初始階段應(yīng)力應(yīng)變曲線差別不大，符合線粘彈性規(guī)律；在相同拉伸速率條件下，隨著溫度的升高，應(yīng)力值有一個(gè)減小的趨勢，溫度越低，應(yīng)力水平越高；當(dāng)溫度不變時(shí)，拉伸速率越大，應(yīng)力值越大。數(shù)值擬合結(jié)果與實(shí)驗(yàn)、理論都具有相同的規(guī)律，相對偏差不超過10%。

用MATLAB對松弛實(shí)驗(yàn)進(jìn)行數(shù)值仿真，結(jié)果如圖8所示。結(jié)果表明，從拉伸實(shí)驗(yàn)數(shù)據(jù)推導(dǎo)出的本構(gòu)模型同樣適用于非拉伸情況，且相對偏差不超過10%。

3.2 本構(gòu)模型模型三維ABAQUS驗(yàn)證

利用ABAQUS用戶材料子程序UMAT對本構(gòu)模型進(jìn)行二次開發(fā)，將本文提出的本構(gòu)模型編寫成UMAT子程序，在ABAQUS中進(jìn)行本構(gòu)模型的三維算例驗(yàn)證。

將本構(gòu)模型從一維擴(kuò)展到三維[12-13]如式(12)。

Aijkl(E(0)-E(t))εkl

(12)

其中，Aijkl為一個(gè)四階張量，此處可表示為矩陣形式：

(13)

(14)

粘性項(xiàng)的應(yīng)力計(jì)算可按照線彈性部分進(jìn)行三維計(jì)算，需注意一點(diǎn)：粘性項(xiàng)是狀態(tài)量，不可按照增量形式進(jìn)行計(jì)算。

應(yīng)力更新過程分為2步：(1)計(jì)算增量步結(jié)束時(shí)的彈性應(yīng)力。t時(shí)刻傳入的應(yīng)力σ(t)為真實(shí)應(yīng)力，欲得到增量步結(jié)束時(shí)的彈性應(yīng)力，需先得到增量步開始時(shí)的彈性應(yīng)力值，加上上一步減去的粘性應(yīng)力項(xiàng)，才可得到增量步結(jié)束時(shí)的彈性應(yīng)力。(2)計(jì)算增量步結(jié)束時(shí)的真實(shí)應(yīng)力。用步驟(1)中得到的彈性應(yīng)力減去增量步結(jié)束時(shí)的粘性應(yīng)力，便可得到增量結(jié)束時(shí)的真實(shí)應(yīng)力。增量步結(jié)束時(shí)真實(shí)應(yīng)力為

(15)

則Δt時(shí)間段內(nèi)應(yīng)力增量為

Δσij(t+Δt)=σij(t+Δt)-σij(t)

=[Dij-E(0)+E(t+Δt)]AijklΔεkl(t+Δt)+

[E(t+Δt)-E(t)]Aijklεkl(t)

(16)

令：

(17)

依據(jù)上述有限元增量形式，將其編寫為UMAT在ABAQ中進(jìn)行數(shù)值仿真計(jì)算。

數(shù)值仿真模型尺寸與圖1所示試件尺寸相同，在ABAQUS中建立模型如圖9所示。模型包括試驗(yàn)件和夾具兩部分。試驗(yàn)件網(wǎng)格尺寸1 mm，網(wǎng)格數(shù)量23 230，單元類型為C3D20H。夾具網(wǎng)格尺寸2 mm，與圖9所示垂直紙面方向尺寸為1 mm，網(wǎng)格數(shù)量5410，材質(zhì)為鋼。計(jì)算過程中固定一個(gè)夾具，另一夾具勻速向另一方向移動(dòng)。

推進(jìn)劑泊松比取為0.496，單元類型選擇雜交單元(Hybrid Formulation)以消除沙漏效應(yīng)，應(yīng)力更新過程有UMAT實(shí)現(xiàn)。

對不同溫度、不同加載速率條件進(jìn)行數(shù)值仿真，取試件中心節(jié)點(diǎn)為研究對象，查看其應(yīng)力-應(yīng)變曲線并與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，結(jié)果如圖10所示。三維數(shù)值仿真結(jié)果符合推進(jìn)劑拉伸實(shí)驗(yàn)結(jié)果，不同溫度、拉伸速率下的變化規(guī)律也符合理論預(yù)測。

從圖10可看出，三維結(jié)果較一維結(jié)果誤差有所增大，相對偏差最大為14.8%，這是因?yàn)槿S計(jì)算在ABAQUS下進(jìn)行，所用應(yīng)力應(yīng)變?yōu)檎鎸?shí)應(yīng)力應(yīng)變，較名義應(yīng)力應(yīng)變誤差增大。

誤差原因分析：

(1)單向拉伸試驗(yàn)得到的應(yīng)力-應(yīng)變均為名義值，轉(zhuǎn)換為真實(shí)值，進(jìn)而擬合本構(gòu)參數(shù)時(shí)存在一定偏差；

(2)泊松比按經(jīng)驗(yàn)取為0.496，與真實(shí)值存在一定偏差；

(3)粘性應(yīng)力耗散項(xiàng)的計(jì)算由松弛模量確定，文中假設(shè)這一項(xiàng)與損傷無關(guān)，而損傷對推進(jìn)劑本構(gòu)模型影響極大。因此，在損傷較大的情況下，松弛模量的形式不變必然導(dǎo)致誤差的存在。

4 模擬發(fā)動(dòng)機(jī)熱-結(jié)構(gòu)分析

根據(jù)本文提出的非線性粘彈性本構(gòu)模型，對貼壁澆注固體火箭發(fā)動(dòng)機(jī)硫化降溫過程進(jìn)行分析。計(jì)算模型如圖11所示，主要結(jié)構(gòu)參數(shù)采取無量綱表示。

在ABAQUS/CAE中建立模型。藥柱、絕熱層、殼體之間采取布爾加運(yùn)算進(jìn)行合并。為控制網(wǎng)格對分析結(jié)果的影響，對關(guān)鍵部位采取了不同的種子控制策略(尺寸及過渡要求)，并對藥柱、絕熱結(jié)構(gòu)等對應(yīng)材料泊松比較大的部位采用雜交單元技術(shù)以消除沙漏效應(yīng)，如圖12所示。在藥柱兩側(cè)施加對稱邊界條件，在后裙端面施加固支條件以消除剛體位移。

推進(jìn)劑、絕熱層、殼體各部件的材料屬性如表3所示。推進(jìn)劑模量部分本文提出的本構(gòu)模型在UMAT進(jìn)行計(jì)算。

熱邊界：藥柱硫化降溫過程是放在一定的溫度環(huán)境中，該環(huán)境溫度為10 ℃；藥柱、絕熱層、殼體初始溫度均為60 ℃，推進(jìn)劑和絕熱層與空氣的對流換熱系數(shù)為1 K/(m2·W)，殼體與空氣的對流換熱系數(shù)為10 K/(m2·W)，時(shí)間歷程設(shè)置為30 d。

表3 推進(jìn)劑、絕熱層、殼體的材料屬性

藥漿澆注后的硫化過程涉及粘合劑基體高分子網(wǎng)絡(luò)的交聯(lián)、粘合劑基體與固體顆粒組分的粘接等一系列化學(xué)-熱力學(xué)耦合過程，反應(yīng)機(jī)理復(fù)雜、影響因素眾多。本文所述硫化降溫過程，以藥漿硫化反應(yīng)達(dá)到“零應(yīng)力溫度”為起點(diǎn)，并假設(shè)：

(1)“零應(yīng)力溫度”直至目標(biāo)環(huán)境溫度的過程中，藥柱熱物理屬性及力學(xué)性能均可由推進(jìn)劑試件的測試結(jié)果進(jìn)行有效表征；

(2)在“零應(yīng)力溫度”直至目標(biāo)環(huán)境溫度范圍內(nèi)，推進(jìn)劑熱物理屬性不隨溫度變化。

受限于實(shí)驗(yàn)測試條件，在進(jìn)行硫化降溫實(shí)驗(yàn)中，僅可測試藥柱中孔處的位移變化，因此選取此作為本構(gòu)模型的驗(yàn)證標(biāo)準(zhǔn)。

取模型上3個(gè)典型位置如圖13所示，其輸出溫度隨時(shí)間變化的數(shù)據(jù)如圖14所示。外壁面與目標(biāo)環(huán)境的對流換熱系數(shù)最大，故溫度下降快；藥柱內(nèi)孔處與目標(biāo)環(huán)境的對流換熱系數(shù)較小，因此溫度下降較為緩慢；肉厚中間處熱損失最少，故溫度下降最為緩慢。

硫化降溫過程結(jié)構(gòu)的變化主要是因?yàn)闇囟葓龅淖兓鸬?，初始階段溫度變化劇烈，所以中孔位移變化也會較快，當(dāng)溫度場平衡后，位移也不會發(fā)生變化。從圖14可看出，約200 h后，溫度場趨于平衡。因此，中孔位移必然也有類似的規(guī)律。將中孔位移隨時(shí)間變化的數(shù)據(jù)繪制在圖15中，數(shù)值仿真結(jié)果與分析一致。

數(shù)值仿真計(jì)算得到的中孔位移最大為0.014 48R(R為殼體半徑)。采用相同結(jié)構(gòu)尺寸的發(fā)動(dòng)機(jī)進(jìn)行了10 ℃環(huán)境下的溫度平衡試驗(yàn)，8 d后測得藥柱內(nèi)孔處溫度與環(huán)境溫度一致，表明藥柱溫度達(dá)到了平衡；此時(shí)進(jìn)行藥柱內(nèi)徑測量，并與設(shè)計(jì)值對比，發(fā)現(xiàn)藥柱內(nèi)孔徑向位移約為0.013 8R。數(shù)值仿真與實(shí)驗(yàn)數(shù)據(jù)的相對偏差僅為4.9%，驗(yàn)證了本構(gòu)模型的準(zhǔn)確性。

5 結(jié)論

(1)本文構(gòu)建的考慮損傷、溫度的非線性粘彈性本構(gòu)模型，由彈性項(xiàng)和粘性項(xiàng)組成；彈性項(xiàng)僅為應(yīng)變的函數(shù)，結(jié)合松弛實(shí)驗(yàn)和拉伸實(shí)驗(yàn)對參數(shù)進(jìn)行了擬合；粘性項(xiàng)為狀態(tài)量，且不受損傷的影響，可由松弛模量和應(yīng)變計(jì)算得到。

(2)通過MATLAB和QBAQUS對本構(gòu)模型進(jìn)行了數(shù)值仿真驗(yàn)證，仿真結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合較好，一維仿真結(jié)果與實(shí)驗(yàn)數(shù)據(jù)相對偏差不超過10%，三維仿真結(jié)果相對偏差略有增大，最大為14.8%。

(3)基于本文所建立的本構(gòu)模型對發(fā)動(dòng)機(jī)硫化降溫過程進(jìn)行了數(shù)值仿真，得到中孔徑向位移為0.144 8R，與實(shí)驗(yàn)測得的數(shù)據(jù)(0.013 8R)相對偏差僅為4.9%，說明本構(gòu)模型在發(fā)動(dòng)機(jī)熱-結(jié)構(gòu)分析具有較好的精度。

[1] 劉中兵,利鳳祥,李越森,等.高過載條件下固體推進(jìn)劑藥柱結(jié)構(gòu)完整性分析計(jì)算[J].固體火箭技術(shù),2003,26(2):12-16.

LIU Zhongbing,LI Fengxiang,LI Yuesen,et al.A calculation analysis for structural integrity of solid propellant grains under high overload[J].Journal of Solid Rocket Technology,2003,26(2):12-16.

[2] Ha K,Schapery R A.A three dimensional viscoelastic constitutive model for particulate composites with growing damage and its experimental validation[J].International Journal of Solid Structures,1998,35(26-27):3497-3517.

[3] Ulrich Mandel,Robin Taubert,Roland Hinterh?lzl.Three-dimensional nonlinear constitutive model for composites[J].Composite Structures,2016,142:78-86.

[4] Muliana A,Rajagopal K R,scharnuterD T,et al.A nonlinear viscoelastic constitutive model for polymeric solids based on multiple natural configuration theory[J].International Journal of Solid Structures,2016,100-101:95-100.

[5] 張曉,鄭堅(jiān),彭威,等.HTPB復(fù)合固體推進(jìn)劑粘彈性應(yīng)變能及非線性本構(gòu)模型[J].固體火箭技術(shù),2015,38(6):827-832.

ZHANG Xiao,ZHENG Jian,PENG Wei,et al.Viscoelastic strain energy and nonlinear constitutive model for HTPB composite solid propellant[J].Journal of Solid Rocket Technology,2015,38(6):827-832.

[6] 常新龍,賴建偉,張曉軍,等.HTPB推進(jìn)劑高應(yīng)變率粘彈性本構(gòu)模型研究[J].推進(jìn)技術(shù),2014,35(1):123-127.

CHANG Xinlong,LAI Jianwei,ZHANG Xiaojun,et al.High strain-rate viscoelastic constitutive model for HTPB propellant[J].Journal of Propulsion Technology,2014,35(1):123-127

[7] 王哲君,強(qiáng)洪夫,王廣,等.中應(yīng)變率下HTPB推進(jìn)劑壓縮力學(xué)性能和本構(gòu)模型研究[J].推進(jìn)技術(shù),2016,37(4):776-782.

WANG Zhejun,QIANG Hongfu,WNAG Guang,et al.Mechanical properties and constitutive model for HTPB propellant under intermediate strain rate compression[J].Journal of Propulsion Technology,2016,37(4):776-782.

[8] 姚東,張光喜,高波.考慮應(yīng)力狀態(tài)的HTPB/AP推進(jìn)劑含損傷熱-粘彈性本構(gòu)方程[J].固體火箭技術(shù),2014,34(4):496-499.

YAO Dong,ZHANG Guangxi,GAO Bo.Constitutive equations involving damage for HTPB/AP propellant considering stress state[J].Journal of Solid Rocket Technology,2014,34(4):496-499.

[9] 唐志平,田蘭橋,朱兆祥.高應(yīng)變率下環(huán)氧樹脂的力學(xué)性能[C]//全國第二屆爆炸力學(xué)會議文集.揚(yáng)州,1981.

TANG Zhiping,TIAN Lanqiao,ZHU Zhaoxiang.Mechanical properties of epoxy resin under high strain rate[C]//Anthology of the Second National Conference on Mechanics of Explosion.Yangzhou,1981.

[10] Staverman A J,Schwarzl F.Nonlinear deformation behavior of high polymers[M].Stuart H A,Ed.Die Phyaile der Hochpolymeren,1956.

[11] 趙榮國,張淳源.一個(gè)考慮損傷的非線性粘彈性本構(gòu)關(guān)系[J].湘潭大學(xué)自然學(xué)科學(xué)報(bào),2001,23(3):28-33.

ZHAO Rongguo,ZHANG Chunyuan.A nonlinear viscoelastic constitutive relation with damage[J].National Science Journal of Xiangtan University,2001,23(3):28-33.

[12] 馮震宙,王新軍,王富生,等.朱-王-唐非線性粘彈性本構(gòu)模型在有限元分析中的實(shí)現(xiàn)及其應(yīng)用[J].材料科學(xué)與工程學(xué)報(bào),2007,25(2):269-272.

FENG Zhenzhou,WANG Xinjun,WANG Fusheng,et al.Implementation and its application in finite element analysis of constitutive model for ZWT nonlinear viscoelastic material[J].Joural of Materials Sciense and Engineering,2007,25(2):269-272.

[13] 彭云.基于 ABAQUS 的非線性粘彈性本構(gòu)模型二次開發(fā)[J].南昌航空大學(xué)學(xué)報(bào),2011,25(2):38-41.

PENG Yun.Developing of nonlinear viscoelastic constitutive model based on ABAQUS[J].Joural of Nanchang Hangkong University (Natural Sciences),2011,25(2):38-41.