含初始缺陷的復合固體推進劑力學性能

封 濤, 許進升, 韓 龍, 陳 雄, 周長省

復合固體推進劑作為一種含能復合材料被廣泛用作運載火箭及各種戰略戰術導彈的動力能源，其力學性能很大程度上影響著導彈的生存能力及作戰能力。早期對復合固體推進劑的研究大多基于連續介質力學，通過實驗從唯象的角度得到推進劑宏觀力學響應的本構關系[1-4]，但這種方法無法對推進劑內在結構變化機理進行有效研究。近年來隨著計算機性能的提高，國內外眾多學者開始著手從細觀角度出發，建立了復合固體推進劑細觀顆粒的填充模型，對推進劑的力學性能進行了定量的分析，并取得了一定的成果。

Matous等[5]最先通過自主開發的Rocpack軟件生成了固體推進劑代表性體積單元，并在顆粒與基體之間的界面層設置了粘接單元來模擬界面脫粘損傷的產生與發展。常武軍等[6]通過實驗證實了顆粒/基體界面是推進劑的薄弱環節，其脫粘是造成宏觀應力-應變曲線非線性的重要原因之一，獲得界面確切的力學性能參數對推進劑的細觀研究工作至關重要。張興高等[7]采用掃描電鏡，并結合界面化學原理，準確測量了AP及HTPB黏合劑基體的接觸角和表面能參數，得到了界面性能與老化之間的有效關系。Zhi等[8-10]在研究推進劑的細觀損傷和非線性力學性能時，用雙線性內聚力模型近似表征推進劑顆粒/基體界面的脫粘損傷，同時還研究了顆粒的體積分數、粒徑大小、位置的隨機分布對推進劑細觀損傷及宏觀力學性能的影響。韓龍等[11]研究了復合固體推進劑細觀界面性能隨拉伸速率的變化規律，通過參數反演的分析算法得到了優化后的界面參數值，該值能準確預測不同加載速率下推進劑的宏觀力學行為。

以往的研究中多認為推進劑為黏合劑基體、增強顆粒和界面相組成的三相復合材料，認為其內部細觀結構完好，沒有考慮到真實情況下推進劑生產過程中存在的初始缺陷的影響。事實上，復合固體推進劑細觀組成與混凝土極為相似，二者經常被用來作比較研究，受混凝土研究中的啟發[12-13]，本研究認為，由于工藝限制，推進劑細觀顆粒/基體界面處存在著包裹不完全的初始缺陷，并對含缺陷的界面力學性能進行了定義，研究當界面存在缺陷時推進劑的宏觀力學行為，得到推進劑的初始模量及拉伸強度隨界面缺陷含量的變化規律。

1 細觀建模及參數獲取

復合固體推進劑是由分散相顆粒、中間相界面和連續相基體組成的高填充比三相復合材料。要從細觀層面上研究其力學性能，首要的是建立一個能符合推進劑真實細觀結構的顆粒填充模型，然后賦予模型相應的力學參數和邊界條件，模擬推進劑在載荷作用下的力學響應。

1.1 細觀建模

1.1.1 代表性體積單元生成

顆粒增強復合材料細觀建模主要有基于蒙特卡羅算法[14]的隨機模擬方法和基于分子動力學算法[15]的碰撞模擬方法，由于復合固體推進劑的顆粒填充體積分數通常較大，可以達到70%～80%，再綜合考慮到算法實現的效率、成熟度以及有限元前處理的工作量等相關因素，本研究選用分子動力學方法建立復合固體推進劑的顆粒填充模型。

表1為復合固體推進劑的典型配方。圖1為推進劑典型配方中顆粒粒徑分布。

表1 復合固體推進劑典型配方Table 1 Typical formula of composite solid propellant

結合表1所示推進劑配方及圖1所示的顆粒分布規律，應用分子動力學方法對復合固體推進劑顆粒填充模型建模，如圖2所示。由于顆粒的填充體積分數較高，顆粒間大小差異很大，在這種模型上直接生成網格進行有限元計算是非常困難的。為了簡化計算模型，參考文獻[16]通過Mori-Tanaka等效的方法將Al顆粒對推進劑力學性能的影響等效到復合基體中，在計算時只單獨考慮AP顆粒的影響。

1.1.2 邊界條件

基于周期性假設認為，推進劑在外載荷的作用下，其內部的應力場和應變場存在連續性和周期性，因此在利用有限元方法預測復合材料固有屬性參數時，需要對代表性體積單元（representative volume element，RVE）施加周期性邊界條件以滿足宏觀上的連續均質假設，并確保得到合理的細觀應力、應變分布。周期性邊界可以表示為：

式中：L為RVE單元的尺寸；和為施加在邊界上的位移載荷。

由于周期邊界條件須在周期性網格的基礎上方能實現，在實際應用中較為不便，為了能在有限元軟件ABAQUS中既準確又簡便地模擬推進劑的單軸拉伸力學行為，本研究選用均勻位移邊界條件，如圖2（a）所示。均勻位移邊界條件下，RVE在受載變形過程中各邊始終保持平直，圖2（a）中虛線為變形后的RVE邊界位置。上邊界所給的均勻位移載荷為0.1667 mm/s，對應的應變率為0.003333 s–1，同時給右邊界施加一個限制x方向位移一致的耦合約束。文獻[17]研究表明，使用均勻位移邊界與周期性邊界計算的結果不存在明顯的差異，并且隨著RVE尺寸的增加誤差漸漸趨近于零。

1.1.3 網格劃分

本研究是基于推進劑二維顆粒填充模型展開，并視為平面應變問題進行數值求解。在有限元計算中，AP顆粒被視為彈性體，其模量相對基體材料大得多，近似認為在受載時不變形，于是在網格劃分時對其采用四節點的平面應變單元CPE4，并在界面附近處將網格細化以提高計算的收斂性。

在有限元分析中，除平面應力問題外，需采用雜交單元來模擬不可壓縮材料（泊松比為0.5）或近似不可壓縮材料（泊松比大于0.475）的響應。HTPB基體材料為近似不可壓縮材料，在受載荷時體積保持恒定，不隨應力應變狀態而發生改變，在計算時無法利用節點處的位移插值函數來獲得單元的壓應力，而雜交單元毋須通過節點位移插值，通過自身設定的額外自由度便可直接獲取單元的壓應力，其偏應力及偏應變的獲取則仍然依賴于節點的位移場計算。因此對基體材料的網格劃分選用平面應變4節點四邊形線性積分的雜交單元CPE4H，模型網格劃分如圖2（b）所示。

在顆粒與基體界面處設置粘接單元（cohesive element）模擬顆粒與基體的粘接特性，并沿界面的徑向設置掃掠網格，為保證粘接單元的收斂，在計算時設置較小的初始載荷步并對粘接單元設定一定的黏性。通過在界面處定義摩擦罰函數來施加通用接觸算法以防止界面脫粘后顆粒與基體可能會發生重疊或相互滲透。

1.2 組分參數獲取

復合推進劑對應的組分參數為黏合劑基體的松弛模量和固體填充顆粒的彈性模量及泊松比。其中顆粒的相關參數依據文獻[5]選取，AP：彈性模量32447 MPa，泊松比 0.1433；Al：彈性模量 68300 MPa，泊松比0.33。

對推進劑基體材料的等速率拉伸應力響應，本研究從線黏彈性理論出發，結合對基體材料松弛實驗[18]測得的應力響應，將松弛模量以Prony級數形式進行擬合，具體擬合表達式為[19]：

式中：為平衡模量；和分別為第i個Maxwell單元的模量和松弛時間；t為時間。松弛曲線擬合結果見圖3。由此可得HTPB基體的黏彈性松弛行為描述，進而得到基體基于松弛模量的應力-應變本構關系[1]：

式中：為基體應力響應，為應變率。

2 界面缺陷單元定義

內聚力模型（cohesive zone model，CZM）建立在彈塑性斷裂力學的基礎上，最早由Dugdale[20]和Barenblatt[21]在研究脆性材料的斷裂時提出來的，其基本思想是將界面理想化為具有一定粘接強度等力學特性的無厚度面，并通過牽引力-位移法則（tractionseparation law，TSL）定義粘接單元的力學響應，從而表征整個界面處的損傷起始和演化形式。

2.1 雙線性內聚力模型

在有限元ABAQUS中，內聚區采用一層厚度近似為0的內聚力單元表示，將內聚力單元嵌入到傳統單元之間使其上下表面與相鄰單元連接，外載荷引起的材料損傷只發生在內聚力單元中，周圍單元不受影響。通過定義內聚力單元的TSL本構關系，以實現內聚力單元在外力作用下損傷的起始與演化。雙線型損傷內聚力模型由于結構簡單，在較復雜的工程環境下也能很好地表征內聚區損傷演化等力學行為，被廣泛地應用于各種研究之中，典型的雙線型損傷內聚力模型[22]示意圖如圖4所示。

其滿足的牽引力-位移法則如下：

式中：下標，t 分別代表界面的法向和切向；為內聚應力；為界面的分離位移；和分別表示法向和切向的內聚強度，且；是臨界位移，表示內聚應力達到內聚強度時對應界面的分離位移；為定義的損傷變量。

二維界面的內聚力單元中存在兩個應力分量和位移分量，載荷作用初期，界面應力和位移為線性關系，當應力達到界面損傷初始準則后，界面處產生損傷，其應力和位移關系不再線性變化，本研究采用式（6）中的二次應力作為損傷起始準則：

損傷變量的定義如式（6）：

式中：為失效位移，表示界面分離位移的最大值，此后界面破壞失效，失去承載能力。

采用式（7）斷裂能準則判斷界面是否破壞失效：

式中：，分別為界面的法向和切向的內聚能，為界面參數（取值為2）。

2.2 初始缺陷定義

復合固體推進劑是典型的多相混合非均勻材料，在其制備過程中由于材料特性及制備工藝限制，使得推進劑中會隨機分布著大量的初始缺陷，包括：在機械混合時無法嚴格控制真空環境，使得有少許空氣進入黏合劑基體而導致孔洞和氣泡；在加入固體填充顆粒后攪拌至混合均勻過程中，顆粒間因相互碰撞導致的顆粒破碎；更多的缺陷是在固化降溫過程中，由于顆粒和基體導熱性能差異使得二者變形不一致而導致的顆粒/基體界面缺陷。這些初始缺陷的存在降低了推進劑的結構強度和剛度，嚴重降低了其力學性能。因此需要對含初始缺陷的復合固體推進劑顆粒填充模型進行數值計算，從細觀尺度探究初始缺陷的存在對推進劑力學性能的影響。

以往有關推進劑的宏觀力學實驗證實[6]，顆粒/基體界面作為推進劑結構中最薄弱的一環，在載荷作用下界面的脫粘是導致推進劑失效破壞的根本原因。本研究只考慮推進劑顆粒/基體界面處的初始缺陷，并注意到推進劑在固化降溫過程中形成的應力集中與顆粒間距有關這一事實，對界面缺陷作如下假設：

（1）對于某一給定位置的顆粒，在其周邊兩倍最大顆粒粒徑的同心圓域內，計算其他顆粒與該顆粒之間的距離，如果該間距小于設定的最小距離，即時認為在該顆粒界面處存在初始缺陷，本研究中，的取值為最小填充顆粒的粒徑。

（2）初始缺陷在兩顆粒外公切線區域內的界面上均勻分布，如圖5所示，采用弱化的界面單元來表征初始缺陷，定義缺陷系數為，則含缺陷的界面力學性能為，，。

（3）設界面初始缺陷含量為，其定義為，其中為單個界面缺陷單元周長；為所有界面單元總周長。

2.3 界面力學參數確定

細觀界面力學性能參數的取值直接影響到數值計算結果的準確性，目前還不能通過實驗技術測得其確切值，對于參數的選取更多的是依據文獻和經驗來預估，但是這樣得到的模型參數往往并不精確，因此，眾多學者提出了基于Hooke-Jeeves的參數反演方法[23]，其基本原理是通過不斷地調整模型參數的取值來改變仿真結果，最終使仿真結果與實驗值的誤差小于規定值的過程。

為反演需要，對復合固體推進劑進行單軸拉伸實驗，不考慮推進劑力學性能的溫度相關性和應變率相關性，實驗在室溫下進行，采用的拉伸速率為10 mm/min，對應的應變率為 0.003333 s–1，得到推進劑的工程應力-應變曲線，結果如圖6所示。

界面力學模型中主要有三個模型參數，界面臨界位移、內聚強度和破壞位移。經研究發現，改變臨界位移只會影響仿真結果中應力-應變曲線的線性上升段的斜率，曲線的下降段不會受到影響。內聚強度的取值主要影響應力-應變曲線的峰值載荷區域，究其原因，內聚強度是彈性段與損傷段的分界點，當應力達到內聚強度之后，損傷起始，界面單元力學性能開始發生衰減。破壞位移表示的是界面單元失效時界面層上下端的分離位移，其對應力-應變曲線下降段有顯著的影響，取值越大曲線下降段上載荷值越大。針對這三種參數對應力-應變曲線的不同影響特性，確定了分步反演的思路，構建表征實驗曲線與仿真曲線重合度目標函數：

當目標函數的值小于選定容差時，認為當前界面參數仿真所得結果與實驗結果吻合，界面參數的初始值選取及反演優化后的值列于表2中。

表2 界面內聚參數Table 2 Cohesive parameters of interface

2.4 界面參數驗證

在2.3節已通過推進劑的單軸拉伸實驗反演出粘結完好的顆粒/基體界面的內聚參數，但還需要進一步驗證參數的準確性。考慮到固體火箭發動機裝藥實際工作過程中承載的復雜性及多樣性，采用多階段的加載方式進行驗證實驗。首先將推進劑試樣以10 mm/min的速率拉伸至10%應變，接著松弛1 min，再以10 mm/min的速率對試樣加載至20%應變。同時，運用本研究構建的細觀模型和相關界面參數，對該實驗過程進行有限元仿真分析，輸出驗證曲線，驗證實驗曲線與仿真曲線，如圖7所示。由圖7可以看出，實驗結果與仿真結果吻合度較高，曲線的大體走勢基本一致，雖然兩者仍存在一定的偏差，但在一定程度上仿真結果能夠反映真實實驗過程中的載荷變化情況，說明本研究通過反演優化算法獲得的界面參數是準確可靠的，能用以描述推進劑顆粒與基體界面的粘接力學性能。

3 算例分析

采用商業有限元軟件ABAQUS，建立了4種顆粒含量均為0.72具有不同分布的復合固體推進劑細觀顆粒填充模型（圖8），并對所建立的模型進行數值計算。為滿足對比分析的需要，模型（a）暫未考慮初始缺陷的影響，模型（b），（c），（d）均基于2.2節中的假設定義初始缺陷，并計算得到其對應的缺陷含量分別為10%，20%和30%。限于篇幅，這里僅取界面缺陷含量為20%的模型（c）為例，圖9給出了相應初始缺陷分布及載荷作用下應力-應變云圖。

從圖9中可以看出，當整體應變為5%時，由于推進劑組分材料屬性的不同以及顆粒之間的相互作用，導致推進劑內部應力和應變分布很不均勻。顆粒作為彈性體，其模量相對于基體和界面大得多，在拉伸過程中，顆粒的變形較小，基體的變形較大，這就使得顆粒/基體界面成為整個推進劑中的薄弱環節，在載荷作用下很容易就會產生界面脫粘現象。尤其當顆粒/基體界面處存在著初始缺陷時，即使很小的外力也能使這些弱界面處產生很明顯的脫粘形貌，從而產生微孔洞。

當應變達到15%，顆粒/基體界面處的變形進一步增加，在定義的缺陷界面處開始出現脫粘，形成明顯的微孔洞。究其原因，推進劑中顆粒與基體形變的差異以及顆粒的相互作用很大程度決定了推進劑內部的應力和應變分布，進而影響了顆粒與基體粘接界面處的脫粘過程。由于在定義的缺陷界面處其局部應變會先達到相應的失效位移，導致界面破壞失去承載能力，形成明顯的孔洞。在以往文獻中多提出大顆粒附近界面會較先產生脫粘，其根本原因是在界面性能相同的前提下，較大顆粒界面處的局部應變更大，最先達到界面的失效位移致使界面失效，這與本研究得出的弱界面最先產生脫粘失效雖看似不同，但兩者的本質是一樣的，都是局部應力超過了當前的界面強度引起的。

整體應變持續增加達到50%，這時推進劑中因界面脫粘產生的微孔洞處就開始出現明顯的裂紋直至破壞。因為隨著載荷的進一步增加，界面進一步擴展，顆粒與基體的分離程度也不斷增加，不僅定義的弱界面幾乎都已經發生脫粘失效，在較大顆粒處粘接完好的界面也因局部應變達到失效位移而失效。界面完全脫粘后，顆粒不再承受應力，從而導致在因脫粘形成的裂隙尖端附近存在明顯應力集中現象，慢慢地這些由于界面脫粘而導致基體應力集中處的微孔洞逐漸匯合，最終使基體撕裂導致裂紋產生，進而造成推進劑的宏觀破壞。由圖9中真實推進劑拉斷的電鏡掃描圖可見，紅框標記的脫粘形貌與仿真結果吻合較好，驗證了本研究方法的可靠性。

將仿真結果中模型上邊界各節點在加載方向上的約束反力求和即為上邊界所受的拉力，經處理可得推進劑仿真所得的工程應力-應變曲線，將其與前節的實驗曲線的平均值進行對比，如圖10所示。

比較圖9中推進劑拉伸應力云圖和圖10應力-應變曲線可以看出：在較小應變（）時，推進劑處于線彈性段，應力隨著應變的增加表現出線性增長，在該階段推進劑內部顆粒/基體界面處雖有較大應變，但仍未達到界面的失效位移，即使存在缺陷的顆粒界面處也無脫粘現象。繼續加載，當應變達到15%時，由圖9可知，此時推進劑內部缺陷界面處存在著密集的脫粘，界面承載能力下降，在宏觀應力-應變曲線上此時開始出現非線性，進入“應力平臺區”。載荷繼續增加到圖9所示的50%應變時，這時即便當初粘接完好的顆粒基體界面也產生了嚴重的脫粘，脫粘處的顆粒喪失承載能力，黏合劑基體完全承擔外載荷，而后推進劑基體被不斷拉長，最終因微孔洞的匯合造成基體撕裂，在很短時間內推進劑宏觀破壞。

比較圖10中不含缺陷的模型仿真曲線與推進劑的宏觀實驗曲線，可以看出：經過反演優化獲得的仿真曲線與實際實驗曲線吻合度較高，但同時也存在著一定的差距。這可能是因為：（1）本研究的推進劑細觀仿真是建立在二維平面應變基礎上的，這與真實推進劑三維結構應力狀態存在著一定的偏差；（2）所建立的細觀模型中所有固體顆粒形貌均是規則的圓形，實際推進劑中可能含有一些不規則的顆粒，相比于規整的球體這些顆粒附近更容易出現應力集中，在載荷作用下更容易產生脫粘現象，這也是為什么實驗應力-應變曲線始終在仿真曲線下方的原因。

比較圖10中不同界面缺陷含量的推進劑仿真結果曲線，不難發現：（1）界面缺陷的存在使得推進劑的初始模量降低了，并且缺陷含量越高其模量下降得越劇烈，這是因為含有缺陷的界面傳遞載荷的能力相對較低，在相同的小應變下，其對應的應力也較低，缺陷含量越高則表現的更明顯；（2）界面缺陷的存在降低了推進劑的拉伸強度，這主要是因為存在缺陷的顆粒/基體界面在拉伸過程中會更早出現脫粘，造成推進劑應力下降，缺陷含量愈高則下降愈多。圖11為初始模量、拉伸強度與缺陷含量關系。初步探索表3中推進劑的初始模量及拉伸強度隨缺陷含量的變化規律，根據圖11中散點分布，認為其可能為線性或指數關系，在嘗試選用線性擬合時發現當缺陷含量為1時，推進劑的初始模量和拉伸強度均為負，這違背了物理規律，因此本研究選用指數函數擬合其關系，相應的擬合結果見圖11。

表3 推進劑力學性能Table 3 Mechanical properties of composite solid propellant

4 結論

（1）構建了復合固體推進劑的細觀顆粒填充模型，并基于合理的假設定義了含缺陷界面處的力學性能參數，通過仿真計算發現缺陷界面處會更早達到損傷條件發生脫粘，缺陷的存在加速了推進劑的斷裂失效的進程。

（2）比較了不同缺陷含量對應的推進劑應力-應變曲線，發現界面缺陷的存在降低了推進劑的力學性能，缺陷含量越高則力學性能下降越明顯，推進劑的初始模量及拉伸強度隨缺陷含量的增加呈指數下降的趨勢。

［1］JUNG G D, YOUN S K. A nonlinear viscoelastic constitutive model of solid propellant[J]. International Journal of Solids & Structures, 1999, 36（25）: 3755-3777.

［2］XU F, ARAVAS N, SOFRONIS P. Constitutive modeling of solid propellant materials with evolving microstructural damage[J]. Journal of the Mechanics & Physics of Solids, 2008, 56（5）: 2050-2073.

［3］ 王哲君, 強洪夫, 王廣, 等. 固體推進劑力學性能和本構模型的研究進展[J]. 含能材料, 2016, 24（4）: 403-416.

（WANG Z J, QIANG H F, WANG G, et al. Review on the mechanical properties and constitutive models of solid propellants[J]. Chinese Journal of Energetic Materials,2016, 24（4）: 403-416.）

［4］BECKER E B. Constitutive modeling of high-elongation solid propellants[J]. Journal of Engineering Materials and Technology, 1992, 114: 111.

［5］MATOU? K, INGLIS H M, GU X, et al. Multiscale modeling of solid propellants: from particle packing to failure[J]. Composites Science & Technology, 2007, 67（7/8）: 1694-1708.

［6］ 常武軍, 鞠玉濤, 王蓬勃. HTPB推進劑脫濕與力學性能的相關性研究[J]. 兵工學報, 2012, 32（3）: 261-266.

（CHANG W J, JU Y T, WANG P B. Research on correlation between dewetting and mechanical property of HTPB propellant[J]. Acta Armamentarii, 2012, 32（3）:261-266.）

［7］ 張興高, 張煒, 蘆偉, 等. HTPB推進劑填料/基體界面粘結性能老化特性研究[J]. 含能材料, 2009, 17（3）: 269-273.

（ZHANG X G, ZHANG W, LU W, et al. Aging characteristics of interfacial adhesive property of filler/binder matrix for HTPB propellant[J]. Chinese Journal of Energetic Materials, 2009, 17（3）: 269-273.）

［8］ZHI S J, SUN B, ZHANG J W. Multiscale modeling of heterogeneous propellants from particle packing to grain failure using a surface-based cohesive approach[J]. Acta Mechanica Sinica, 2012, 28（3）: 746-759.

［9］ 職世君, 曹付齊, 申志彬, 等. 復合固體推進劑顆粒脫濕損傷參數反演[J]. 推進技術, 2016, 37（10）: 1977-1983.

（ZHI S J, CAO F Q, SHEN Z B, et al. Parameters inversion of particle dewetting damage of composite solid propellants[J]. Journal of Propulsion Technology, 2016,37（10）: 1977-1983.）

［10］職世君, 孫冰, 張建偉. 基于表面粘結損傷的復合固體推進劑細觀損傷數值模擬[J]. 推進技術, 2013, 34（2）:273-279.

（ZHI S J, SUN B, ZHANG J W. Numerical simulation of solid sropellant mesoscopic damage using surfacebased cohesive approach[J]. Journal of Propulsion Technology, 2013, 34（2）: 273-279.）

［11］韓龍, 許進升, 周長省. HTPB/IPDI復合固體推進劑細觀界面率相關參數的反演識別研究[J]. 含能材料, 2016,24（10）: 928-935.

（HAN L, XU J S, ZHOU C S. Inverse identification of the rate-dependent micro interface parameters of HTPB/IPDI composite propellant[J]. Chinese Journal of Energetic Materials, 2016, 24（10）: 928-935.）

［12］王娟, 管巧艷, 孔宇田, 等. 界面初始缺陷對混凝土單軸抗拉強度的影響[J]. 水利水電科技進展, 2016, 36（5）:46-49.

（WANG J, GUAN Q Y, KONG Y T, et al. Influences of initial interfacial defects on uniaxial tensile strength of concrete[J]. Advances in Science and Technology of Water Resources, 2016, 36（5）: 46-49.）

［13］杜修力, 揭鵬力, 金瀏. 考慮初始缺陷影響的混凝土梁動態彎拉破壞模式分析[J]. 工程力學, 2015（2）: 74-81.

（DU X L, JIE P L, JIN L. Dynamic flexural-tensile failure mode analysis of concrete beam with initial defect[J].Engineering Mechanics, 2015（2）: 74-81.）

［14］FISH J. Multiscale modeling and simulation of composite materials and structures[J]. Multiscale Methods in Computational Mechanics, 2011, 55: 215-231.

［15］HUMPHREYS D D, FRIESNER R A, BERNE B J. A multiple-time-step molecular dynamics algorithm for macromolecules[J]. Journal of Physical Chemistry, 2002,98（27）: 6885-6892.

［16］馬昌兵. 復合固體推進劑細觀結構建模及其力學行為數值模擬[D]. 西安: 第二炮兵工程學院, 2011.

（MA C B. Research on mesostructured modeling and mechanical behaviors numerical simulation of composite solid propellant[D]. Xi’an: The Second Artillery Engineering University, 2011.）

［17］PIERARD O. Micromechanics of inclusion-reinforced composites in elasto-plasticity and elasto-viscoplasticity:modeling and computation[D]. Leuven: University Catholique de Louvain, 2006.

［18］常武軍. 復合固體推進劑細觀損傷及其數值仿真研究[D]. 南京: 南京理工大學, 2013.

（CHANG W J. Research on microstructural damage and its numerical simulation method for composite solid propellant[D]. Nanjing: Nanjing University of Science and Technology, 2013.）

［19］XU J, HAN L, ZHENG J, et al. Finite element implementation of a thermo-damage-viscoelastic constitutive model for hydroxyl-terminated polybutadiene composite propellant[J]. Mechanics of Time-Dependent Materials,2017: 1-19.

［20］DUGDALE D S. Yielding of steel sheets containing slits[J]. Journal of the Mechanics & Physics of Solids,1960, 8（2）: 100-104.

［21］BARENBLATT G I. The mathematical theory of equilibrium cracks in brittle fracture[J]. Advances in Applied Mechanics, 1962, 7: 55-129.

［22］ELICES M, GUINEA G V, GóMEZ J, et al. The cohesive zone model: advantages, limitations and challenges[J]. Engineering Fracture Mechanics, 2002, 69（2）: 137-163.

［23］JEEVES T A. “Direct Search” solution of numerical and statistical problems[J]. Journal of the Acm, 1961, 8（2）:212-229.