考慮細觀損傷的推進劑粘彈性多尺度本構模型研究

烏布力艾散·麥麥提圖爾蓀，周 濤，吳艷青，侯 曉，周程哲

（1.北京理工大學爆炸科學與技術國家重點實驗室， 北京 100081； 2.中國航天科技集團第四研究院西安航天化學動力有限公司， 陜西 西安 710025； 3.北京理工大學宇航學院， 北京 100081）

0 引 言

固體推進劑作為固體火箭發(fā)動機的關鍵組成部分，其力學性能被視為主要性能指標之一，對發(fā)動機的綜合性能有著重要影響［1］。準確預測推進劑在不同工況下的力學響應對于發(fā)動機的結構完整性及壽命評估具有重要的理論和工程價值。作為一種高填充顆粒復合材料，推進劑通常由高分子黏合劑和含能填料組成，在細觀層面呈現(xiàn)出明顯的非均勻性，導致其力學性能不僅與外部載荷（如溫度、應變率、圍壓和變形歷史）相關［2-3］，而且與內(nèi)部因素（如細觀結構特征、組分參數(shù)、相界面性質(zhì)）密切相關［4-5］。為了揭示推進劑變形過程中的應力響應及損傷演化過程，通常采用含顆粒-相界面-黏合劑“三相”細觀模型［6-8］。細觀模型能夠直觀地描述組份含量、顆粒級配及初始缺陷等對推進劑力學性能影響，從而為深入理解推進劑損傷演化的物理機制提供重要參考。然而，確定推進劑相界面力學參數(shù)的程序復雜、成本高，計算效率較低［9-10］；而且由于各組分之間的剛度不匹配導致計算收斂性困難，因此目前的研究僅針對相對較低的固體含量（0.1～0.4）或較小的變形（如應變1%～12%）進行，與實際工況存在較大的差距［11］。此外，過多的微觀細節(jié)會限制細觀模型在發(fā)動機裝藥的應力分析中的適用性。

為了兼顧推進劑細觀損傷演化及計算效率，Xu等［12］基于細觀力學均勻化方法，采用“自下而上”的思路提出了一種以孔隙率為損傷變量的粘彈性本構方程；而Hur 等［13］、王貴軍等［14］在此基礎上，引入循環(huán)加載損傷函數(shù)和溫度項進一步提高了模型適用范圍，然而損傷函數(shù)表達式和參數(shù)需要通過宏觀試驗反復擬合，適用的變形范圍較小，對于復雜應力狀態(tài)的適用性有限。Lei 等［15-16］采用單顆粒幾何參數(shù)“歸一化裂紋長度”作為損傷變量表征了顆粒/黏合劑相界面脫粘程度，并將溫度/應變率相關失效判據(jù)引入模型中，實現(xiàn)了非線性粘彈性本構與失效斷裂統(tǒng)一，但尚未討論模型在循環(huán)載荷及圍壓作用下的實用性。Wubuliaisan等［17］從細觀模型中量化提取“脫粘比”和體積膨脹率并引入到應變能密度函數(shù)中，表征了脫濕及黏合劑自身損傷對推進劑力學性能的影響；模型在一定程度上把細觀結構變化與推進劑的宏觀力學行為聯(lián)系起來，實現(xiàn)了推進劑力學性能的多尺度分析，然而并未對模型在復雜載荷（如不同溫度）下的預測能力進行討論。

考慮到推進劑大變形和實際工程需求，Tun? 等［18］在連續(xù)介質(zhì)理論框架中引入細觀損傷函數(shù)，基于剛度折減法建立了推進劑超彈性-粘彈性本構模型。該模型將應變能分解為體量與偏量部分推導應力，并引入基于空隙率的經(jīng)驗公式以描述推進劑應力軟化、循環(huán)損傷等特性。Yun 等［19］采用與Tun? 類似的思路，基于粘彈性脫濕準則建立了推進劑含損傷本構模型，并成功捕捉了雙軸板條加載和變溫下的體積膨脹和應力軟化現(xiàn)象。Kumar 等［20］開發(fā)了一種超彈粘彈性本構模型，揭示了推進劑從不可壓縮轉變?yōu)榭蓧嚎s狀態(tài)過程。Kantor 等［21］通過將應力三軸度相關損傷準則引入超彈粘彈性本構中，分析了推進劑損傷演化與應力狀態(tài)的相關性。孟紅磊等［22］、Xu 等［23］、Wang 等［24］及Talebi 等［25］基于Schapery 理論［26］，通過引入了損傷內(nèi)變量和軟化函數(shù)開發(fā)了推進劑含損傷本構模型，并采用ABAQUS 軟件UMAT 子程序進行了計算。上述模型大多采用“自上而下”的方法，主要關注推進劑的宏觀力學特性，雖然在一定程度上能呈現(xiàn)溫度、應變率、圍壓及循環(huán)加載對推進劑力學性能及損傷演化的影響，然而，由于經(jīng)驗性模型參數(shù)缺乏明確物理意義，難以刻畫推進劑真實細觀損傷演化過程。

本研究基于推進劑細觀結構及粘彈性特性，在有限變形下提出一種考慮細觀損傷的粘彈性多尺度本構模型。該模型考慮推進劑力學性能的溫度、應變率、壓力相關性及循環(huán)加載過程中的應力軟化。首先基于商用軟件ABAQUS 編寫了UMAT 子程序，并通過試驗數(shù)據(jù)確定了模型參數(shù)。之后采用不同載荷下的高能推進劑試驗數(shù)據(jù)對模型的預測能力及可靠性進行了驗證。該模型所需參數(shù)較少且實施難度小，可為發(fā)動機推裝藥結構完整性的多尺度分析提供一定的理論指導。

1 含損傷粘彈性多尺度模型

本研究約定A、B、C分別代表標量、二階張量和四階張量；上標T 表示張量的轉置，上標點表示對時間的導數(shù)。在有限變形下，變形梯度張量F的偏量部分可表示為［27］：

1.1 超彈性框架

推進劑往往呈現(xiàn)出大變形且卸載后形變基本恢復的力學特性，因此可通過應變能密度函數(shù)推導其（超）彈性響應。首先將應變能密度函數(shù)分解為體量φ和偏量φ部分之和［17］：

1.2 粘彈性響應

根據(jù)對應原理［26］，粘彈性響應可通過對上述彈性PK2 應力進行積分獲得：

式中，k(t)和g(t)是體積和剪切松弛函數(shù)；一般情況下取k(t) =g(t)。

當考慮溫度效應時，基于時-溫等效原理可將式（11）和（12）中真實時間t和ξ替換為：

式中，T為溫度，單位是K；aT為平移因子。

至此，含損傷粘彈性PK2 及對應的Cauchy 應力可表示為：

1.3 細觀損傷

當推進劑在外載荷下發(fā)生損傷時，相界面脫濕、黏合劑撕裂或顆粒斷裂導致微空洞的形成，從而導致其剪切模量和體積模量降低。假定固體推進劑的代表性體積單元中微空洞含量為c、顆粒含量Vp，那么基于細觀力學Hashin-Shtrikman 上下界均質(zhì)化方法可得推進劑當前的剪切模量［12］：

式中，f(0)表示推進劑初始狀態(tài)。

推進劑體積模量也可以采用類似的方法獲得，然而由于推進劑幾乎不可壓縮特性，即使發(fā)生損傷后其體積模量比剪切模量大幾個數(shù)量級，因此采用唯象方法表示推進劑體積模量的變化［28］：

式中，γ為材料常數(shù)，表示推進劑可壓縮程度。

根據(jù)試驗觀測及理論分析，假設微空洞演化隨著形變增大而隨著靜水壓Svol減?。?，18］，則：

1.4 切線模量

隱式求解中需要提供切線模量，即雅可比矩陣?？紤]到帶損傷粘彈性模型切線模量推導過程復雜，本節(jié)采用類似于應力推導的步驟，即先考慮彈性切線模量Ce并分解為體量部分和偏量部分［27］：

考慮損傷時的彈性切線模量Ce d為：

結合式（24）～（26）可得帶損傷粘彈性材料切線模量C為：

從式（19）～（27）可以看出，有限變形下含損傷粘彈性模型切線模量計算相當復雜，為了提高模型實用性，可采用ABAQUS 內(nèi)置UHYPER 子程序令軟件自動計算上述切線模量。

1.5 模型的選擇

式（4）～（18）呈現(xiàn)了含細觀損傷粘彈性模型的通用形式，實際應用過程中可根據(jù)材料應力應變特性選擇合適的函數(shù)簡化計算。為了簡單起見，采用Neo-Hookean 模型［5］作為應變能密度函數(shù)的偏量部分，即：

式中，c1，c2為常數(shù)，Tref是參考溫度。

2 試驗方法及參數(shù)確定

2.1 試驗方法

選用硝酸酯增塑聚醚（NEPE）推進劑作為研究對象進行試驗，其顆粒含量為69.5%。 首先根據(jù)GB770B-2005 所規(guī)定的方法制取啞鈴型推進劑試件，其標距及橫截面分別為70 mm 及10 mm×10 mm，如圖1a 所示。之后采用Instron5967 通用材料試驗機進行測試，每種條件進行3 次重復試驗，取其平均作為最終結果。試驗測試溫度為223～333 K，其中低溫試驗使用液氮進行冷卻；為了確保溫度均勻性，被測試樣在溫控箱中至少保存40 min，如圖1b 所示。單軸拉伸、應力松弛及循環(huán)加載試驗加載速率為100 mm·min-1。為了在應力松弛試驗中避免出現(xiàn)初始損傷，根據(jù)行業(yè)標準QJ 2487-1993 將初始恒定應變設定為5%。單軸拉伸試驗過程中利用HIKVISION MV-CH120 相機進行同步記錄，并采用VIC-2D 系統(tǒng)進行位移場處理，最后采用數(shù)字圖像相關技術（Digital image correlation， DIC）計算拉伸過程中試樣體積變化。單軸拉伸條件下體積變化J可表示為［29］：

圖1 單軸拉伸試驗Fig.1 Experimental setup used in the uniaxial loading test

式中，V0和V是初始及當前體積，ε||為沿拉伸方向的應變，ε⊥為垂直于拉伸方向的應變。

2.2 模型參數(shù)確定

由于應力松弛函數(shù)僅在5%的應變下進行，認為推進劑尚未出現(xiàn)損傷。因此，首先參考行業(yè)標準QJ 2487-1993，通過對不同溫度下的應力松弛曲線進行平移確定推進劑松弛模量主曲線，并通過曲線擬合確定松弛函數(shù)式（30）的系數(shù)，如圖2a 所示。此外，根據(jù)不同溫度下松弛曲線的平移距離，獲得式（31）中WLF 方程的系數(shù)，如圖2b 所示。

圖2 應力松弛主曲線及WLF 方程擬合Fig.2 The master relaxation curve and fitted WLF form

對于式（18）中微空洞演化，可以通過原位CT 等試驗直接獲?。?0］，然而由于缺乏相關數(shù)據(jù)，本研究將通過擬合298 K 及100 mm·min-1條件下的單軸拉伸試驗應力應變曲線及體積變化曲線獲得損傷相關參數(shù)，如圖3 所示。

圖3 損傷參數(shù)的確定Fig.3 Determination of the progressive damage parameters

如前文所述，由于松弛函數(shù)參數(shù)在無損傷情況下獲得，通過應力應變曲線擬合的參數(shù)僅有3 個，即本模型可通過較少的參數(shù)實現(xiàn)不同加載下的應力響應表征。模型參數(shù)如表1 所示。

表1 NEPE 推進劑模型參數(shù)Table 1 Model parameters for NEPE propellant

3 結果與討論

3.1 單軸拉伸響應預測

NEPE 推進劑在不同加載速率和溫度下的試驗應力應變曲線及模型預測結果如圖4 所示。由圖4a 可見，NEPE 推進劑力學行為具有明顯的應變率相關性，即初始模量及拉伸強度隨著加載速率增大而增大。盡管在斷裂應變附近存在一定的偏差，模型預測結果與整體上試驗吻合；考慮到預測結果涵蓋的加載速率范圍較寬，模型在分析推進劑大變形方面仍展現(xiàn)出了良好的預測能力。圖4b 呈現(xiàn)了推進劑在寬溫度條件下的應力應變曲線，隨著溫度降低，推進劑初始模量及拉伸強度增加；在不同溫度下模型預測結果與試驗數(shù)據(jù)吻合，進一步說明該模型具有良好的預測能力。

圖4 在單軸拉伸加載條件下試驗及本構模型預測結果Fig.4 Model predictions and experimental results for uniaxial tensile tests

為了更進一步展示模型的預測能力，對NEPE 推進劑循環(huán)加載響應進行了預測，加載歷史與試驗一致，如圖5a 所示。圖5b 對試驗及仿真結果進行了對比，模型能夠較好地預測循環(huán)加載過程中的應力相應。圖5c 呈現(xiàn)了循環(huán)加載過程中的體積變化，在卸載過程中體積有所減小，這是由于推進劑粘彈性恢復特性導致，然而如前文所述，雖然微空洞含量可以減小，由于已脫濕相界面失去承載能力，式（14）中的損傷并不減小，從而在循環(huán)加載卸載過程中會出現(xiàn)應力軟化現(xiàn)象。此外，圖5d 則呈現(xiàn)了在第三、第五循環(huán)（對應圖5c 虛線處）時刻DIC 及仿真應變場分布，模型預測結果的合理性進一步說明該模型也能復現(xiàn)推進劑循環(huán)加載損傷特性。

圖5 循環(huán)加載試驗與模型預測結果對比Fig.5 Comparison between cyclic loading test results and model predictions

3.2 圍壓作用預測

在發(fā)射程中，發(fā)動機裝藥經(jīng)歷點火建壓，導致燃燒室內(nèi)壓強快速攀升。為了分析圍壓對推進劑力學響應的影響，采用Li 等［3］冷增壓試驗結果對本構模型適用性及預測能力進行了考察。首先采用圍壓為5.4 MPa應力應變曲線擬合了模型參數(shù)，然后對圍壓分別為0，0.5 MPa 及2.0 MPa 工況進行了仿真，如圖6a 所示。從圖6a 可見，由于模型考慮了靜水壓效應（式（18）），模型能夠較好地預測不同圍壓下的力學響應。圖6b呈現(xiàn)了在不同圍壓下推進劑微空洞含量的變化，說明圍壓的存在一定程度上抑制了細觀損傷，即微空洞含量越小，推進劑斷裂應變越大。這種現(xiàn)象可以解釋為什么火箭發(fā)動機在點火過程中更容易發(fā)生故障，因為推進劑的老化會增加微空洞含量，導致斷裂應變較?。?8］。

圖6 圍壓對推進劑力學響應的影響Fig.6 Effects of superimposed pressure on the response of the propellant

3.3 雙軸拉伸響應預測

為了驗證本構模型在復雜應力狀態(tài)下的適用性，對Ranjan 等［31］雙軸試驗進行了仿真。試驗中采用的十字形雙軸拉伸試樣的尺寸如圖7a 所示。首先采用1 mm·min-1的測試數(shù)據(jù)（圖7b）以及Ranjan 等［32］應力松弛測試數(shù)據(jù)確定了模型參數(shù)。然后在不同加載速率下進行了仿真，由圖7b 可見模型成功預測了雙軸加載下推進劑應力響應。在雙軸拉伸過程中，試樣表面上形成了一些微空洞，這些微空洞的區(qū)域呈現(xiàn)出較亮的顏色（圖8 中的黃色虛線），圖8 呈現(xiàn)了（X∶Y）為1∶1 和0.5∶1 的比例加載下（加載速率1000 mm·min-1）的微空洞分布預測，結果顯示該模型能夠較好地預測損傷路徑，進一步說明該模型在復雜應力狀態(tài)下的預測能力。

圖7 雙軸拉伸試驗結果與模型預測對比Fig.7 Comparison between biaxial tensile test results and model prediction

圖8 不同比例加載下試樣變形圖Fig.8 Deformation of the sample at different loading rate ratio

3 結 論

（1）在寬溫和加載速率下對NEPE 推進劑進行了單軸拉伸及應力松弛試驗。試驗結果表明，隨著加載速率增加和溫度降低，NEPE 推進劑初始模量及拉伸強度增加。此外，圍壓能夠在一定程度上抑制微空洞演化，從而使推進劑呈現(xiàn)更高的拉伸強度及延伸率。

（2）基于推進劑細觀特征，在有限變形下開發(fā)了考慮細觀損傷的超彈粘彈性多尺度本構模型，該本構模型考慮了溫度、應變率、圍壓及循環(huán)加載應力軟化等因素；基于微空洞演化和細觀均質(zhì)化理論，實現(xiàn)了細觀損傷向宏觀本構的傳遞。

（3）將模型應用到ABAQUS 中對NEPE 推進劑力學性能進行了預測，結果表明該模型能夠準確預測溫度、應變率和圍壓對NEPE 推進劑力學性能的影響，且能夠很好地描述推進劑在循環(huán)加載損傷及雙軸拉伸下的力學響應。