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

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

0 引 言

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

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

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

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

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

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

1.1 超彈性框架

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

1.2 粘彈性響應(yīng)

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

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

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

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

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

1.3 細(xì)觀損傷

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

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

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

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

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

1.4 切線模量

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

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

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

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

1.5 模型的選擇

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

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

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

2.1 試驗方法

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

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

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

2.2 模型參數(shù)確定

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

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

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

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

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

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

3 結(jié)果與討論

3.1 單軸拉伸響應(yīng)預(yù)測

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

圖4 在單軸拉伸加載條件下試驗及本構(gòu)模型預(yù)測結(jié)果Fig.4 Model predictions and experimental results for uniaxial tensile tests

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

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

3.2 圍壓作用預(yù)測

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

圖6 圍壓對推進(jìn)劑力學(xué)響應(yīng)的影響Fig.6 Effects of superimposed pressure on the response of the propellant

3.3 雙軸拉伸響應(yīng)預(yù)測

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

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

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

3 結(jié) 論

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

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

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