含損傷非線性粘彈性本構模型及數(shù)值仿真應用①

孟紅磊，鞠玉濤

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

0 引言

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

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

1 本構方程

1.1 彈性粘彈性對應原理

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

式中ER為參考彈性模量，與彈性模量單位一致，引入參考彈性模量，主要是為了能使偽應變同樣是無量綱的，ER可任意選取，通常選取1值或選用材料的初始彈性模量值作為ER。

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

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

1.2 含累積損傷的非線性本構方程

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

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

本文結合累積損傷的概念，建立含累積損傷的非線性粘彈性本構方程，并證明含累積損傷的非線性粘彈性本構方程能較好地反映不同應變率條件下及復雜應變率條件下的應力應變響應，含累積損傷的非線性粘彈性本構方程為

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

則含損傷的三維粘彈性本構方程表示為

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

式中 σeq為Mises等效應力。

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

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

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

2 數(shù)值仿真應用

2.1 本構方程增量形式及切線模量求解

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

由主應力 σ1、σ2、σ3求解 Von Mises等效應力，并計算損傷內變量值S。

偽應變由偏偽應變部分和體積偽應變部分組成

則應力張量分量形式等于

對偏偽應變進行展開，可得到

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

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

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

同理，體積偽應變展開可得到

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

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

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

對損傷內變量S進行離散，可得

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

由式(24)得到

由式(13)得到

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

由式(18)得

由式(27)得

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

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

總的偽應變增量可表示為偏偽應變增量和體積偽應變增量之和的形式:

進而可得到

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

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

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

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

令

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

由式(41)可得

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

2.2 UMAT子程序介紹

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

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

(2)提供應力增量和應變增量之間的關系，用于更新應力。

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

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

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

圖1 等速拉伸過程應力-應變理論預測與實驗對比Fig.1 Comparison between prediction and experimental measurement at constant strain rate tension

3 結論

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

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

［2］鐘濤，張為華.點火瞬態(tài)過程對復合推進劑力學響應特性的影響［J］.國防科技大學學報，2004，26(3):2004.

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

［4］隋欣，魏志軍.炮射導彈發(fā)射過程發(fā)動機裝藥強度分析［J］.彈道學報，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復合推進劑含損傷和老化本構研究［J］.推進技術，2002，23(6):509-512.

［13］彭威，鄭堅，白鴻柏，等.復合推進劑微裂紋損傷本構模型研究［J］.固體火箭技術，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］.推進技術，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.