基于遺傳算法的復合推進劑粘彈性本構參數研究

鄭 健，韓 波，周長省

（南京理工大學 機械工程學院，南京210094）

端羥基聚丁二烯（hydroxyl terminated polybutadiene，HTPB）推進劑是一種常用的復合推進劑，其力學行為十分復雜，其性能依賴于溫度和應變率，為了能夠準確地描述出該種推進劑裝藥在受力狀態下的力學行為，需要使用粘彈性本構模型。粘彈性本構參數獲取的準確性影響仿真分析的準確性，而常規粘彈性實驗參數的獲取方法均存在較大的實驗誤差。為了解決這一問題，Zapas－Phillips和Sorvaris等人［1－2］在研究了由試驗機加載產生的誤差問題之后，分別使用近似積分處理的方法來獲取較為準確的松弛模量數據。Zapas－Phillips法使用積分中值定理來獲得修正的松弛模量數據，但是僅能得到t≥t1／2時刻的松弛模量數據，其中，t為松弛時間，t1為達到恒應變的時間。Sorvaris法在此基礎上進一步研究、推導出了t≥0時刻的松弛模量數據。孟紅磊等人［3］在Sorvari方法的基礎上提出了改進型的Sorvari方法，即引入修正因子α來更加精確地近似積分過程。許進升等人［4］基于直接積分方法來獲取復合HTPB推進劑的松弛模量數據。

本文提出了一種基于遺傳算法的HTPB復合推進劑粘彈性本構參數獲取方法，為固體火箭發動機裝藥結構完整性分析提供準確的本構參數。

1 本構參數獲取

1.1 問題描述

HTPB推進劑是一種典型的粘彈性材料，其力學特性呈現出彈性固體和粘性流體的性質，在固體火箭發動機裝藥完整性分析中常使用線粘彈性本構模型來表征推進劑的力學行為。三維積分型的各向同性線粘彈性材料應力－應變關系可以表示為

應力－應變關系的應力偏量部分sij和體變部分σkk分別為

式中：λ（t），G（t），K（t）分別為拉梅常數對應的松弛函數、剪切松弛函數和體積松弛函數；δij為位移分量，eij為偏應變張量，εkk為應變偏量。對于各向同性線粘彈性材料，必須具備2種材料函數才能準確描述其應力－應變特性。但是通過實驗直接獲得粘彈性材料的材料函數G（t）和K（t）比較困難，而材料的松弛模量E（t）較容易獲得。因此，如果獲得了材料函數之間的轉化關系，就可以通過較容易測量的材料函數來表征線粘彈性材料的力學行為。

將式（2）進行Laplace變換，得到：

式中：s為Laplace變換的復變量。彈性力學中材料常數存在以下關系：

式中：E為松弛模量，ν為泊松比。

將式（4）中的參數進行Laplace變換，得到：

對式（5）、式（6）進行Laplace逆變換可以得到線粘彈性材料中E（t），ν（t）和G（t），K（t）之間的關系。實際工程應用中，可以根據實際情況簡化材料函數之間的變化關系。在固體推進劑有限元計算中經常采用下列假設：由于泊松比ν（t）隨時間變化很小，故假設泊松比為常數，則根據松弛模量E（t）可以得到剪切松弛模量G（t）和體積松弛模量K（t）：

因此，如何準確地獲得推進劑的松弛模量對結構完整性分析十分重要。

1.2 傳統實驗方法

現行的復合推進劑松弛模量測定方法是使用拉伸試驗機快速地將推進劑試樣拉至某一恒定應變，同時測量出推進劑的松弛應力，從而擬合出該推進劑材料的靜態松弛模量［1－2］。

圖1給出了階躍位移和加載示意圖。

圖1 階躍位移加載示意圖

圖1 （a）表示試驗機系統的理想加載過程，即不考慮加載歷程。但真實情況下試驗機通過調節控制系統的PID參數，實際的上升階段可能存在超調、緩慢爬坡等現象，即存在加載歷程，如圖1（b）所示。且HTPB推進劑的松弛響應較快，因此采用圖1（a）這種數據處理方法會造成實驗獲得的靜態松弛模量數據偏小。為了解決這一問題，很多研究學者提出了各種數據處理方法，但這些方法均假設試驗機的加載階段近似為一個直線上升的階段，即如圖1（c）中所示的加載歷程。國外Zapas－Phillips采用積分中值定理獲得修正的t≥t1／2時刻的松弛模量數據，其計算式如下：

Sorvaris［1－2］在 Zapas－Phillips的研究基礎上，推導出了t≥0時刻的松弛模量數據表達式：

式中 ：σ為應力為應力變化率 ，ε為應變為應變變化率。

國內孟紅磊［3］在Sorvari方法的基礎上提出了改進型的Sorvari方法，即引入修正因子α來更加精確地近似積分過程：

許進升等人［4］基于直接積分方法來獲取復合HTPB推進劑的松弛模量數據，即通過理論推導出Prony級數形式的松弛函數在圖1（c）加載條件下松弛應力的顯式表達式，再通過最小二乘法來獲得靜態松弛模量。

雖然這些方法在一定程度上提高了實驗數據的準確性，但是與實驗過程中真實的加載過程仍存在一定的差異。

1.3 基于遺傳算法的松弛模量獲取

使用基于遺傳算法的全階段松弛模量數據擬合方法解決實驗系統和數據處理帶來的松弛數據不準確問題。遺傳算法（genetic algorithm，GA）是由美國Michgan大學的Holland J教授于1975年最早提出，它是一種借鑒了生物學進化規律的算法，該算法模擬了生物學的遺傳、突變、自然選擇和雜交等現象。對于一個優化問題，首先從解空間中隨機抽取一定數量的候選解形成一個初始種群，使用一種編碼方法對初始種群進行編碼。上一代種群內部通過遺傳、變異、雜交等過程形成新一代種群，通過算法的適應度函數來判斷新種群的優劣，通過這一過程來模擬自然界的優勝劣汰。經過若干代的優勝劣汰獲得最優的種群，即最優解。

遺傳算法的基本步驟如下：

①編碼。使用某種編碼方式將解空間映射到遺傳算法的編碼空間中。

②產生初始種群。在編碼空間中隨機產生一定數量的碼形成初始種群，這些初始種群對應著一組初始解。

③適應度評估。使用適應度函數來判斷種群中個體的優劣。

④遺傳算子的確定。遺傳算子中確定了選擇、交叉、變異的整個過程。

⑤算法的終止條件。確定遺傳算法的終止條件。

使用MATLAB編制了擬合HTPB推進劑松弛曲線的計算程序，以擬合函數和實驗曲線的重合度作為遺傳算法中的目標函數，目標函數的表達式為

式中：f為適應度，n為計算過程中的迭代次數，即種群繁衍的代數；σl（X）代表優化過程中解空間內某一個體在給定載荷作用下得到的應力值；σe為實驗獲得的應力值。使目標函數表達式獲得最小值的X即為最優解，即為所需要的松弛模量擬合參數。為避免目標函數在尋優過程中數值過小影響計算精度，將式（11）乘以一個較大的系數k可以將目標函數適當放大（此處k取10 000）。

2 實驗方法和結果處理

為了獲得準確的推進劑松弛模量數據，使用美國BOSS動態力學實驗機進行實驗，該實驗機具有加載速度快、采集系統精度高等特點。為了能夠滿足動態力學實驗機的材料安裝要求，參照文獻［5］中的單軸拉伸實驗方法設計了單軸應力松弛試樣。

圖2為松弛實驗的試樣示意圖。使用刀片將HTPB推進劑切割成長條形試樣，試樣尺寸B＝5mm，W＝5mm，H＝30mm。使用改性丙烯酸膠粘劑將HTPB推進劑試樣上下兩端和木質接頭粘接，試樣制作完成后放入密封干燥箱中在室溫下靜置24h進行固化。實驗設備如圖3所示，整個實驗系統由控制柜、動態力學分析儀DMA、保溫箱、液氮罐和非接觸式應變測量系統組成。HTPB推進劑模量較低且容易造成損傷，不宜直接使用引伸計來測量拉伸應變。直接使用實驗機的加載位移來反推試樣應變存在較大的誤差，因此實驗中采用基于數字圖像相關性的非接觸式應變測量系統來測量拉伸過程中推進劑的應變。

圖2 松弛試樣示意圖

圖3 DMA實驗系統

以松弛實驗過程中的實際位移變化曲線作為輸入條件開始遺傳算法的計算，圖4給出了優化計算過程中的優化歷程，圖中，n為計算過程中種群繁衍的代數，f為種群在優化過程中的適應度。圖中給出了適應度的最優值和平均值。從圖中可以看出，在種群的繁衍過程中種群的適應度不斷下降，這表明種群的質量在逐步提高，經過292代之后種群的最佳適應度之間達到一個穩定值，并認為第292代中的最優個體即為最優解。計算結果如表1所示。表中，E0＝8MPa（瞬態松弛模量），N為Prony級數的階數，EN，tN分別為各階松弛模量、松弛時間。

圖4 遺傳算法計算過程中適應度最優值和平均值變化情況

表1 松弛實驗數據

分別使用傳統計算方法和遺傳算法獲得了HTPB推進劑的松弛模量數據，將該結果應用于真實松弛實驗的加載過程中，所得的應力σ隨時間變化結果如圖5所示。從圖中可以看出，遺傳算法結果和實驗結果吻合良好，而傳統方法所得到的結果明顯小于實驗結果；隨著時間的增加，傳統方法和遺傳算法的差距逐漸減小。由于傳統方法忽略了實驗過程中加載階段對實驗結果的影響，因此其獲得的結果偏低。但是加載階段對于HTPB推進劑的平衡模量的影響隨著時間逐漸減小，因此傳統方法和遺傳算法的差距會隨著時間的推移逐漸減小，直至重合。

圖5 遺傳算法和傳統算法結果對比

3 結論

本文建立了基于遺傳算法的復合推進劑松弛模量獲取方法，為該推進劑的粘彈性本構模型提供了較精確的本構參數，獲得如下結論：①由于拉伸試驗機系統存在加載歷程，傳統松弛模量獲取方法獲取的數據較真實值偏??；②基于遺傳算法的松弛模量擬合方法考慮了真實實驗過程中的加載歷程，因此可以準確地獲得復合推進劑的松弛模量。

［1］SORVARI J，MALINEN M.On the direct estimation of creep and relaxation functions［J］.Mechanics of Time－dependent Materials，2007（11）：143－157.

［2］SORVARI J，MALINEN M.Numerical interconversion between linear viscoelastic material function with regularization［J］.Solids and Structures，2007，44：1 291－1 303.

［3］孟紅磊.改性雙基推進劑裝藥結構完整性數值仿真方法研究［D］.南京：南京理工大學，2011.MENG Hong－lei.Research on numerical simulation method of structure integrity analysis for modified double base propellant grain［D］.Nanjing：Nanjing University of Science and Technology，2011.（in Chinese）

［4］許進升，鞠玉濤，鄭健.復合固體推進劑松弛模量的獲取方法［J］.火炸藥學報，2011，34（5）：58－62.XU Jin－sheng，JU Yu－tao，ZHENG Jian.Acquisition of the relaxation modulus of composite propellant［J］.Chinese Journal of Explosive ＆ Propellants，2011，34（5）：58－62.（in Chinese）［5］JUNG G D，YOUN S K.A nonlinear viscoelastic constitutive model of solid propellant［J］.International Journal of Solids and Structures，1999，36（25）：3 755－3 777.