一種固體推進劑藥柱結構完整性的快速評估方法

張 路，余 瑞，鄧康清，龐愛民，楊 玲

(湖北航天化學技術研究所，湖北 襄陽 441003)

引 言

固體火箭發(fā)動機在生產和工作期間，會受到各種載荷作用，這些載荷作用產生的應力和應變可能會導致發(fā)動機失效，因此需要對藥柱結構完整性進行分析[1]。其中，模量是衡量固體推進劑藥柱變形能力的重要指標，準確獲得模量變化對藥柱結構完整性評估具有非常重要的意義[2]。目前，通常將推進劑的模量分為靜態(tài)特性中的應力松弛模量和動態(tài)特性中的儲能模量進行研究，應力松弛模量一般通過推進劑試樣的靜態(tài)試驗得到[3-6]；儲能模量一般通過使用動態(tài)力學分析(DMA)儀的試驗得到[7-8]。針對等應變率和等應力率載荷條件下藥柱等效模量變化的情況，M. R. Lajczok[9]通過將推進劑松弛模量的Prony級數形式代入線黏彈性本構方程中，得到了等效模量隨時間變化的函數。但對于更一般的非等應力率、應變率點火增壓和固化降溫等動態(tài)載荷下藥柱等效模量的理論分析還需要進一步的研究。

對于固化降溫和點火增壓過程中藥柱的結構完整性分析，大都是基于推進劑的各種本構模型應用有限元軟件進行數值模擬分析。如于洋[10]、鄧斌[11]基于熱黏彈性本構模型，得到了變溫和增壓過程中藥柱應力、應變隨時間的變化情況; Shiang-Woei Chyuan[12]分別基于非線性熱黏彈性本構方程和線黏彈本構方程，分析了溫度載荷和點火增壓載荷下藥柱的瞬態(tài)應力、應變；Tejas Nikam[13]基于線彈性模型應用ANSYS有限元軟件分析了點火增壓條件下不同藥型對藥柱應力、應變的影響。通過有限元數值模擬分析，得到藥柱的應力、應變分布，進而可以得到藥柱模量的變化情況，不同載荷條件下藥柱模量的變化也不同。有限元數值模擬方法雖然可以評估藥柱的模量變化和結構完整性，但對特定發(fā)動機在不同載荷工況(不同固化降溫曲線、不同點火增壓曲線)下藥柱的結構完整性進行分析時，需多次模擬計算且每次計算耗時長，時效性不高。

本研究基于對藥柱的黏彈性分析，建立動態(tài)載荷作用下藥柱等效模量的實時評估方法，并基于溫度和壓力載荷下藥柱應變解析的理論分析，建立藥柱等效應變和等效應力的計算模型，可實現藥柱結構完整性的評估。然后應用該評估方法借助于數值計算軟件和有限元仿真軟件對固化降溫方程和點火增壓下藥柱的結構完整性進行分析，并與基于線黏彈性理論的有限元模擬計算結果進行對比。

1 計算方法

1.1 固化降溫過程中藥柱等效模量的計算方法

固化降溫過程中藥柱的溫度變化方程為T(t)，將溫度作用時間t劃分為n份，假設每一份時間內為一個階躍溫度，第i份的作用時間為ti-ti-1，溫度為Ti。這樣，藥柱的溫度變化方程就等效為一個階躍函數，溫度隨時間變化曲線如圖1所示。當時間步長ti-ti-1足夠小時，在誤差允許范圍內，認為兩種曲線重合。溫度變化為ΔTi時，時間ti-ti-1內藥柱產生的應變?yōu)棣う舏。時間t內所有溫度變化對藥柱在t時刻產生的總應變?yōu)椋?/p>

ε(t)=Δε1+Δε2+Δε3+…+Δεn

(1)

圖1 固化降溫過程中溫度隨時間的變化曲線Fig.1 The changing curve of temperature with time during cooling process after curing

已知，推進劑在參考溫度Ts下的松弛模量為：

(2)

WLF方程為：

(3)

式中：αT為溫度-時間轉換因子；C1、C2為常數。

根據時溫等效原理，恒定溫度T下的時間t在基準溫度Ts下的等效時間為：

(4)

因此，利用基準溫度Ts下的主曲線和溫度-時間轉換因子，可得到推進劑在任何恒定溫度T下的松弛模量ET(t)，表達式為：

(5)

而變溫過程中的時間t在參考溫度Ts下的等效時間可表示為：

(6)

固化降溫過程中，時間步ti-ti-1內的溫度變化ΔTi在t時刻對藥柱所產生的應力為E(ξn-ξi-1)·Δεi進行玻爾茲曼疊加就得到了t時刻藥柱產生的總應力σ(t)，可表示為：

(7)

定義藥柱在任一時刻的總應力與總應變的比值為等效模量，由式(1)、式(7)得到t時刻藥柱的等效模量，具體表達式為：

(8)

由于在溫度載荷作用下，藥柱的應變僅與溫差和推進劑的物理特性有關[14]，對于已知的推進劑藥柱，其物理特性不變，藥柱的應變與溫差成線性關系，可表示為：

ε(t)=St1·T+St2

(9)

式中：St1、St2為溫度應變系數，它們與發(fā)動機材料的幾何參數和物理特性有關，與藥柱的模量無關，對于給定的發(fā)動機，St1、St2為常數(注：不同應變形式的溫度應變系數不同)。

對于簡單圓管形發(fā)動機藥柱，溫度應變系數可通過理論計算公式求得；對于復雜藥形發(fā)動機藥柱，可通過有限元軟件加載簡單的穩(wěn)態(tài)溫度載荷激勵，再對求解的應變進行擬合得到溫度應變系數。

得到固化降溫過程中藥柱等效模量和溫度等效應變方程以后，根據線彈性理論，將對應的等效模量和等效應變相乘就得到了藥柱的等效應力，從而可實現固化降溫下藥柱結構完整性的分析。

當每個時間步內溫度變化相同時，即：

ΔT1=ΔT2=ΔT3=…=ΔTn

(10)

從而，每個時間步內溫度變化使藥柱產生的應變都是恒定且相等的，即：

Δε1=Δε2=Δε3=…=Δεn

(11)

此時，式(8)可化簡為：

(12)

綜上所述，由式(12)便可求出固化降溫過程中藥柱在不同時刻的等效模量。本研究將以上公式通過數值計算軟件Matlab編制程序，把溫度變化方程輸入該程序中，就可實現固化降溫過程中推進劑藥柱模量變化的快速評估。

1.2 點火增壓過程中藥柱等效模量的計算方法

發(fā)動機在點火過程中，藥柱會承受燃燒室內燃氣所帶來的增壓載荷。壓強會在很短時間內上升到最大，并保持不變，壓力作用的整個過程中藥柱的模量都會發(fā)生變化。這里主要討論點火增壓過程中藥柱的模量變化，不考慮溫度帶來的影響。點火增壓過程中增壓方程為p(t)，將增壓時間t劃分n份，認為每一份時間內為階躍增壓，第i份的作用時間為ti-ti-1，壓力為pi。這樣，點火增壓就等效為一系列的階躍增壓，壓強隨時間變化曲線如圖2所示，當時間步足夠小時，在誤差允許范圍內認為兩條曲線重合。

圖2 點火增壓過程中壓強隨時間變化曲線Fig.2 The changing curve of pressure with time during ignition pressurization process

通過實驗測量得到藥柱在參考溫度下的松弛模量主曲線，其Prony級數表達式見式(2)。已知，藥柱在階躍壓力載荷pi作用時的應變與壓強和松弛模量有關，可表示為：

ε(t)/pi+Sp1/E(t)+Sp2

(13)

式中：Sp1、Sp2為藥柱的壓力應變系數，它們與發(fā)動機材料的幾何參數和物理特性有關，與藥柱的力學性能無關，對于給定的發(fā)動機，Sp1、Sp2為常數(注：不同應變形式的壓力應變系數不同)。

文獻[15]中給出了理想發(fā)動機的Sp1、Sp2的計算公式，但對于真實的發(fā)動機應用公式計算會帶來較大的誤差，可通過數值模擬軟件計算得到藥柱在階躍壓力載荷下的應變分布，再對應變進行擬合得到。

第i份的壓力變化Δpi在t時刻所產生的應變?yōu)棣う舏，可表示為：

Δεi=Δpi·(Sp1/E(tn-ti-1)+Sp2) (t0=0)

(14)

值得注意的是，點火增壓載荷與溫度載荷下的Δεi是不同的，這里每個時間步內壓強變化產生的應變隨著作用總時間t是變化的。

壓強變化Δpi在t時刻產生的應力為：

Δσi=E(tn-ti-1)·Δεi(t0=0)

(15)

基于波爾茲曼疊加將每份時間內階躍壓力變化對藥柱產生的應力、應變相加，即可得到任意時刻藥柱總的應力、應變，然后即可求出各時刻藥柱的等效模量。t時刻藥柱的等效模量為：

(16)

為計算方便，在實際的計算中，可以將時間t等分為n份，每個時間步為t/n。同樣地，采用數值計算軟件Matlab編制程序，將時間步進行迭代，即可求出在點火增壓過程中藥柱在任意時刻的應力、應變，從而可以得到藥柱等效模量的變化情況。

2 數值模擬

2.1 有限元模型的建立

本研究選擇貼壁澆注式固體圓筒發(fā)動機為研究對象，該圓筒發(fā)動機包含殼體、絕熱層和推進劑藥柱。由于結構的對稱性，對圓筒發(fā)動機的1/12進行三維有限元建模，有限元模型如圖3所示，共劃分5120個單元，27931個節(jié)點。

圓筒發(fā)動機殼體、絕熱層和藥柱的基本材料參數見表1。

圖3 固體火箭發(fā)動機3D有限元模型Fig.3 3D Finite element model of the solid rocket motor

表1 發(fā)動機基本材料參數Table 1 Parameters of basic materials used for motor

注：ρ為密度；E為彈性模量；ν為泊松比；α為熱膨脹系數;λ為熱導率;C為比熱容。

2.2 推進劑力學性能參數

固體推進劑固化溫度為52℃，零應力溫度為60℃。60℃時，推進劑的松弛模量擬合成的Prony級數可表示為：

(17)

式中各個系數見表2。其中，推進劑的初始模量E0為3.861MPa，平衡模量Ee為0.803MPa。

表2 Prony級數的各個系數Table 2 The coefficients of Prony series

溫度—時間轉換因子函數式WLF方程為：

(18)

式中：Ts=333.15K；C1=13.025；C2=201.31。

固化降溫過程中，在發(fā)動機殼體端面的一條邊上對發(fā)動機施加3個方向的位移約束，對稱面滿足對稱邊界條件，不考慮應力過大產生的脫粘。

點火增壓過程中，藥柱內表面為自由邊界，藥柱與絕熱層之間、絕熱層與殼體之間均為粘接邊界條件，剖面上約束環(huán)向位移，藥柱、絕熱層和殼體兩端約束軸向位移，殼體外表面固定約束。

3 結果與討論

3.1 溫度載荷下藥柱等效模量變化的評估

3.1.1 穩(wěn)態(tài)溫度載荷下藥柱有效模量變化的評估

將藥柱放置在0℃環(huán)境溫度下進行保溫，藥柱溫度與環(huán)境溫度一致，保溫過程中假設任意時刻整個藥柱內的溫度相等。根據時溫等效原理，可得到藥柱在0℃下的保溫時間t在參考溫度60℃下的等效時間ξ，如圖4所示。

圖4 藥柱在0℃下的保溫時間t及其等效時間ξ的關系曲線Fig.4 Relation curve of time(t) vs. equivalent time(ξ) of grain at 0℃

將等效時間ξ代入式(17)中可得到0℃下藥柱的松弛模量，其隨保溫時間t的變化如圖5所示。基于線黏彈性本構模型應用有限元軟件ANSYS WORKBENCH對恒定溫度0℃下推進劑藥柱的結構響應進行分析，得到了藥柱內部應力、應變的分布。現取藥柱內表面最大等效應力、等效應變位置A點進行研究，在0℃溫度載荷下數值計算得到的A點等效模量隨時間分布如圖5所示。由圖5可知，藥柱的松弛模量和等效模量基本重合，所以，穩(wěn)態(tài)溫度載荷作用下可認為藥柱的等效模量就是其真實的松弛模量。

圖5 0℃下A點的松弛模量和等效模量分布曲線Fig.5 Distribution curve of relaxation modulus and equivalent modulus at point A of grain at 0℃

3.1.2 固化降溫過程中不同降溫條件下藥柱模量變化的評估

固化降溫過程中，假設發(fā)動機溫度場在任意時刻都是均勻的，藥柱溫度從60℃降至-40℃，降溫所用時間為100h。以下分析3種降溫條件，按溫度與時間的函數關系分為直線降溫，平方降溫和三次方降溫，其溫度變化函數為：

條件1： 直線降溫T1(t)=60-1/3600·t

(19)

條件2： 平方降溫T2(t)=-40+1/1.29×109·(t-3.6×106)2

(20)

條件3： 三次方降溫T3(t)=60-1/4.6656×1014·t3

(21)

3種降溫曲線如圖6所示。降溫過程中3種降溫曲線的等效時間ξ隨降溫時間t的分布情況如圖7所示，圖7中左側縱坐標為直線降溫和平方降溫過程中等效時間大小變化的參考坐標，右側縱坐標為三次方降溫過程中等效時間大小變化的參考坐標。

圖6 3種降溫條件下溫度隨時間分布曲線Fig.6 Distribution curves of temperature with time under three cooling conditions

圖7 3種降溫條件的降溫時間及其等效時間分布曲線Fig.7 Distribution curves of cooling time(t)and equivalent time(ξ) under three cooling conditions

取藥柱內的A點為研究對象。將3種降溫曲線分別輸入數值計算軟件Matlab編寫的程序中，并選取合適的時間步進行運算，可得到相對應時間點藥柱的等效模量。如直線降溫中，取每變化1℃的時間為時間步，即每一步的時間為1h，溫降為1℃，從而可以得到降溫時間段內100個時刻點藥柱的等效模量，現選取其中的一些點繪制在圖8中。當階躍溫度取的越小，時間步就越短，得到的等效模量也就越精確。同樣地，可求出其他兩種降溫條件下藥柱在不同時間點的等效模量，也選取一些點繪制在圖8中。然后，基于線黏彈性本構模型應用有限元軟件ANSYS WORKBENCH分別對3種降溫條件下藥柱的結構完整性進行數值模擬計算，得到不同降溫條件下不同時刻A點的等效應力和等效應變，從而求出各個時刻藥柱的等效模量，其分布情況如圖8所示。

圖8 3種降溫條件下兩種評估方法的等效模量分布Fig.8 Distribution of equivalent modulus obtained by two evaluation methods under three cooling conditions

將利用基于本研究提出的評估方法計算得到的等效模量和利用有限元模擬計算得到的等效模量進行對比，由圖8可知，兩種計算方法得到的等效模量基本重合，從而證明了固化降溫過程中藥柱等效模量快速評估方法的準確性。

3.1.3 固化降溫過程中藥柱結構完整性的快速評估

基于線彈性本構模型應用有限元分析軟件ANSYS WORKBENCH，通過輸入各溫度下相應的等效模量，計算任意穩(wěn)態(tài)溫度下藥柱的等效應力和等效應變。分別對恒定溫度10℃和 -40℃下藥柱的應力、應變進行穩(wěn)態(tài)分析，得到了兩種溫度載荷下藥柱的等效應變分別為0.096252和0.17801。結合初始溫度下藥柱的等效應變?yōu)榱悖瑪M合出該藥柱的溫度應變系數St1為-0.0017800，St2為0.10922，因此藥柱溫度載荷下的等效應變方程為：

εeq(t)=-0.0017800·T+0.10922

(22)

利用該方程就可以得到不同溫度下藥柱的等效應變，然后根據應用Matlab計算得到的各溫度下藥柱的等效模量，可求出不同溫度下藥柱的等效應力，從而達到藥柱結構完整性的評估。3種降溫條件下，由該評估方法得到的藥柱A點的等效應變、等效應力隨時間分布如圖9所示。然后，將基于線黏彈性模型應用有限元分析軟件ANSYS WORKBENCH對藥柱進行數值模擬計算得到的結果同樣見圖9。

圖9 3種降溫條件下藥柱A點的應力、應變曲線Fig.9 Curves of von Mises strain vs. stress at point A of grain under three cooling conditions

由圖9可知，在相同降溫條件下，基于該快速評估方法和基于線黏彈性理論的有限元模擬計算方法得到的藥柱等效應力和等效應變曲線基本重合。在計算精度方面，直線降溫條件下藥柱的等效應變和等效應力的最大差值分別為3.68%和3.83%；平方降溫條件下藥柱的等效應變和等效應力的最大差值分別為4.60%和4.01%；三次方降溫條件下藥柱的等效應變和等效應力最大差值分別為3.68%和4.74%。綜合以上，3種降溫條件下對藥柱進行完整性分析時，應用快速評估方法和有限元模擬計算方法得到的藥柱等效應變和等效應力的差值都在5%以內。在計算效率方面，根據本研究提出的理論評估方法，得到藥柱在溫度載荷下的等效應變方程之后，只要將藥柱的降溫方程輸入編好的程序中，就可以得到任意降溫條件下藥柱等效應力和等效應變的分布；而基于線黏彈性本構模型有限元模擬計算方法，需要針對降溫條件的數量，進行多次的模擬計算，計算所需時間較長。結果表明，此藥柱結構完整性快速評估方法的計算精度和時效性很高。

3.2 壓力載荷下藥柱等效模量變化的評估

3.2.1 穩(wěn)態(tài)壓力載荷下藥柱等效模量動態(tài)變化的評估

點火時，發(fā)動機燃燒室內部壓力在很短時間內上升到發(fā)動機的工作壓強為2MPa，并作用一段時間，該過程中不考慮溫度對藥柱的影響，忽略增壓的時間，認為藥柱遭受階躍內壓力載荷。已知給出了藥柱在零應力溫度下的松弛模量曲線，如圖10所示。基于線黏彈性本構模型應用有限元軟件ANSYS WORKBENCH對恒定壓強2MPa下推進劑藥柱的結構響應進行分析，得到不同時刻藥柱A點的等效應力和等效應變，然后求得A點的等效模量，如圖10所示。

圖10 2MPa下A點的松弛模量和等效模量分布曲線Fig.10 Distribution curve of relaxation modulus and equivalent modulus at point A of grain at 2MPa

由圖10可知，藥柱的松弛模量和等效模量基本重合，壓力載荷作用下可用等效模量來表示藥柱的實際松弛模量。

3.2.2 點火增壓過程中不同增壓條件下藥柱模量變化的評估

固體發(fā)動機在點火增壓過程中，發(fā)動機內部壓力經過約0.1s的建壓時間達到峰值2MPa，然后維持平衡壓強繼續(xù)工作。由于建壓時間短，因此，計算時不考慮燒蝕效應，且忽略建壓過程中溫度的變化，這里只討論增壓過程中藥柱模量隨時間的變化。以下分析3種增壓條件，按壓強與時間的函數分為直線增壓、平方增壓和指數增壓，增壓過程中壓強變化函數為：

條件1： 直線增壓p1(t)=20·t

(23)

條件2： 平方增壓p2(t)=200·t2

(24)

條件3： 指數增壓p3(t)=2·(1-e-60t)

(25)

3種增壓條件下壓強隨時間分布曲線如圖11所示。

圖11 3種增壓條件下壓強隨時間分布曲線Fig.11 Distribution curves of pressure with time under three pressurized conditions

基于線彈性本構模型應用有限元分析軟件ANSYS WORKBENCH，通過改變輸入的藥柱彈性模量，對藥柱進行穩(wěn)態(tài)熱結構分析，計算恒定壓強1MPa下藥柱的等效應變，分別取藥柱的彈性模量為2和3MPa，得到不同彈性模量下藥柱A點的等效應變分別為0.24732和0.16519。然后對其進行擬合，得到藥柱的壓力應變系數Sp1為0.24641，Sp2為0.00045502，因此藥柱在壓力載荷下的等效應變方程為：

εeq(t)/pi=0.24641/E(t)+0.00045502

(26)

利用該方程就可以得到不同階躍壓力下藥柱的等效應變，根據公式(14)、公式(15)可求出每個時間步內壓力變化產生的藥柱等效應力和等效應變，然后得到各時刻藥柱總的等效應力和等效應變，進而獲得了3種增壓條件下藥柱的等效模量，其隨時間分布如圖12所示。為了驗證評估結果的準確性，基于線彈性模型應用有限元分析軟件ANSYS WORKBENCH分別對3種增壓條件下藥柱進行數值模擬計算，分別得到各時刻藥柱A點的等效應力和等效應變，獲得各時刻藥柱的等效模量，其分布情況也如圖12所示。

圖12 3種增壓條件下兩種評估方法的等效模量分布曲線Fig.12 Distribution curves of equivalent modulus obtained by two evaluation methods under three pressurized conditions

將快速評估方法計算得到的等效模量和模擬計算得到的等效模量進行對比，由圖12可知，兩種計算方法得到的等效模量基本重合，從而證明了點火增壓過程中藥柱等效模量快速評估方法的準確性。

3.2.3 點火增壓過程中藥柱結構完整性的快速評估

將以上由藥柱結構完整性快速評估方法和有限元模擬計算方法求得的各個時刻藥柱的等效應力和等效應變分別繪制在圖13中。

圖13 3種增壓條件下藥柱A點的應力、應變曲線Fig.13 Curves of von Mises strain and stress at point A of grain under three pressurized conditions

由圖13可知，在相同增壓條件下，兩種分析方法得到的等效應力和等效應變曲線基本重合。在精度方面，直線增壓條件下藥柱的等效應變和等效應力最大差值分別為1.48%和1.06%；平方增壓條件下藥柱的等效應變和等效應力的最大差值分別為1.14%和1.01%；指數增壓條件下藥柱的等效應變和等效應力的最大差值分別為1.93%和1.23%。綜合以上，3種增壓條件下對藥柱結構完整性進行分析，應用快速評估方法和有限元模擬計算方法得到的藥柱等效應力和等效應變之間的差值都在2%以內。在計算效率方面，快速評估方法在得到藥柱的等效應變方程之后，只需將藥柱的增壓方程輸入編好的程序中，就可以得到任意增壓條件下藥柱的等效應力和等效應變。因此，此藥柱結構完整性快速評估方法具有很高的計算精度和時效性。

4 結 論

(1)建立了固化降溫和點火增壓動態(tài)載荷作用下推進劑藥柱等效模量的理論評估方法。具體根據溫度和壓力載荷加載歷程進行分段，基于波爾茲曼疊加原理計算分段載荷對藥柱應力、應變產生的疊加響應，推導出了固化降溫和點火增壓過程中推進劑藥柱等效模量的計算公式。

(2)提出了溫度應變系數和壓力應變系數的概念，建立了動態(tài)載荷工況下藥柱等效應力和等效應變的計算模型。溫度應變系數和壓力應變系數可通過對有限元模擬計算的結果進行擬合得到，其可作為藥柱本身的材料參數，不管在何種溫度和壓力載荷工況作用下，其溫度和壓力應變系數始終不變。

(3)驗證了動態(tài)模量評估方法和藥柱等效應變、等效應力計算方法的準確性。通過與基于線黏彈性理論的模擬計算方法得到的藥柱內部危險點位置的等效應力和等效應變進行對比，得到藥柱等效應力和等效應變之間的差值都在5%以內，證明了該快速評估方法的可行性。

參考文獻:

[1] Salita M. Modern SRM ignition transient modeling.I. introduction and physical models[C]∥ 37th AIAA Joint Propulsion Conference & Exhibit. Salt Lake City: American Institute of Aeronautics and Astronautics, 2001: 2001-3443.

[2] Bohwi S, Jaehoon K. Estimation of master curves of relaxation modulus and tensile properties for solid propellant [J]. Advanced Materials Research, 2014, 871: 247-252.

[3] Lees S, Knauss W G. A note on the determination of the relaxation and creep data from ramp tests[J]. Mechanics of Time-dependent Materials, 2000(4): 1-7.

[4] Sorvari J,Malinen M. On the direct estimation of creep and relaxation functions[J]. Mechanics of Time-dependent Materials, 2007(11): 143-157.

[5] 許進升，鞠玉濤，鄭健，等．復合固體推進劑松弛模量的獲取方法[J]．火炸藥學報，2011，34(5)：58-62．

XU Jin-sheng, JU Yu-tao, ZHENG Jian, et al. Acquisition of the relaxation modulus of composite solid propellant[J]. Chinese Journal of Explosives & Propellants (Huozhayao Xuebao), 2011, 34(5): 58-62.

[6] Adel W M, Liang Guo zhu. Different methods for developing relaxation modulus master curves of AP-HTPB solid propellant[J].Chinese Journal of Energetic Materials,2017,25(10): 810-816.

[7] Husband D M. Use of dynamic mechanical measurements to determine the behavior of solid propellant[J]. Propellants, Explosives, Pyrotech, 1992, 17: 196-201.

[8] 張昊，龐愛民，彭松．固體推進劑貯存壽命非破壞性評估方法(II)—動態(tài)力學性能主曲線監(jiān)測法[J]．固體火箭技術，2006，9(3)：190-194．

ZHANG Hao, PANG Ai-min, PENG Song. Nondestructive assessment approaches to storage life of solid propellants(II)—master curve of dynamic mechanical property surveillance method[J]. Journal of Solid Rocket Technology, 2006, 9(3): 190-194.

[9] Lajczok M R. Effective propellant modulus approach for solid rocket motor ignition structural analysis[J]. Computers & Structures,1995,56(1):101-110.

[10] 于洋，王寧飛，張平．溫度載荷下帶筋套管形裝藥結構完整性分析[J]．推進技術，2006，27(6)：492-496．

YU Yang, WANG Ning-fei, ZHANG Ping. Structural integrity analysis for the canular solid propellant grains subjected to temperature loading[J]. Journal of Propulsion Technology, 2006, 27(6): 492-496.

[11] Bin Deng, Yan Xie, Guo Jin tang. Three-dimensional structural analysis approach for aging composite solid propellant grains[J]. Propellants, Explosives, Pyrotech, 2014, 39: 117-124.

[12] Shiang-Woei Chyuan. Nonlinear thermoviscoelastic analysis of solis propellant grains subjected to temperature loading[J]. Finite Elements in Analysis and Design, 2002, 38(7): 613-630.

[13] Nikam T, Pardeshi M,Patil A, et al. Structural integrity analysis of propellant in solid rocket motor[J]. International Conference on Ideas, Impact and Innovation in Mechanical Engineering, 2017, 5(6): 896-902.

[14] 鮑福廷，侯曉．固體火箭發(fā)動機設計[M]．北京：中國宇航出版社，2016，1：400-403．

[15] 王元有，胡克嫻，蔡湘芬，等．固體火箭發(fā)動機設計[M]．北京：國防工業(yè)出版社，1984．