基于遺傳算法的復(fù)合推進(jìn)劑粘彈性本構(gòu)參數(shù)研究

鄭 健，韓 波，周長省

（南京理工大學(xué) 機(jī)械工程學(xué)院，南京210094）

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

本文提出了一種基于遺傳算法的HTPB復(fù)合推進(jìn)劑粘彈性本構(gòu)參數(shù)獲取方法，為固體火箭發(fā)動機(jī)裝藥結(jié)構(gòu)完整性分析提供準(zhǔn)確的本構(gòu)參數(shù)。

1 本構(gòu)參數(shù)獲取

1.1 問題描述

HTPB推進(jìn)劑是一種典型的粘彈性材料，其力學(xué)特性呈現(xiàn)出彈性固體和粘性流體的性質(zhì)，在固體火箭發(fā)動機(jī)裝藥完整性分析中常使用線粘彈性本構(gòu)模型來表征推進(jìn)劑的力學(xué)行為。三維積分型的各向同性線粘彈性材料應(yīng)力－應(yīng)變關(guān)系可以表示為

應(yīng)力－應(yīng)變關(guān)系的應(yīng)力偏量部分sij和體變部分σkk分別為

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

將式（2）進(jìn)行Laplace變換，得到：

式中：s為Laplace變換的復(fù)變量。彈性力學(xué)中材料常數(shù)存在以下關(guān)系：

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

將式（4）中的參數(shù)進(jìn)行Laplace變換，得到：

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

因此，如何準(zhǔn)確地獲得推進(jìn)劑的松弛模量對結(jié)構(gòu)完整性分析十分重要。

1.2 傳統(tǒng)實(shí)驗(yàn)方法

現(xiàn)行的復(fù)合推進(jìn)劑松弛模量測定方法是使用拉伸試驗(yàn)機(jī)快速地將推進(jìn)劑試樣拉至某一恒定應(yīng)變，同時測量出推進(jìn)劑的松弛應(yīng)力，從而擬合出該推進(jìn)劑材料的靜態(tài)松弛模量［1－2］。

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

圖1 階躍位移加載示意圖

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

Sorvaris［1－2］在 Zapas－Phillips的研究基礎(chǔ)上，推導(dǎo)出了t≥0時刻的松弛模量數(shù)據(jù)表達(dá)式：

式中 ：σ為應(yīng)力為應(yīng)力變化率 ，ε為應(yīng)變?yōu)閼?yīng)變變化率。

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

許進(jìn)升等人［4］基于直接積分方法來獲取復(fù)合HTPB推進(jìn)劑的松弛模量數(shù)據(jù)，即通過理論推導(dǎo)出Prony級數(shù)形式的松弛函數(shù)在圖1（c）加載條件下松弛應(yīng)力的顯式表達(dá)式，再通過最小二乘法來獲得靜態(tài)松弛模量。

雖然這些方法在一定程度上提高了實(shí)驗(yàn)數(shù)據(jù)的準(zhǔn)確性，但是與實(shí)驗(yàn)過程中真實(shí)的加載過程仍存在一定的差異。

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

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

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

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

②產(chǎn)生初始種群。在編碼空間中隨機(jī)產(chǎn)生一定數(shù)量的碼形成初始種群，這些初始種群對應(yīng)著一組初始解。

③適應(yīng)度評估。使用適應(yīng)度函數(shù)來判斷種群中個體的優(yōu)劣。

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

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

使用MATLAB編制了擬合HTPB推進(jìn)劑松弛曲線的計算程序，以擬合函數(shù)和實(shí)驗(yàn)曲線的重合度作為遺傳算法中的目標(biāo)函數(shù)，目標(biāo)函數(shù)的表達(dá)式為

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

2 實(shí)驗(yàn)方法和結(jié)果處理

為了獲得準(zhǔn)確的推進(jìn)劑松弛模量數(shù)據(jù)，使用美國BOSS動態(tài)力學(xué)實(shí)驗(yàn)機(jī)進(jìn)行實(shí)驗(yàn)，該實(shí)驗(yàn)機(jī)具有加載速度快、采集系統(tǒng)精度高等特點(diǎn)。為了能夠滿足動態(tài)力學(xué)實(shí)驗(yàn)機(jī)的材料安裝要求，參照文獻(xiàn)［5］中的單軸拉伸實(shí)驗(yàn)方法設(shè)計了單軸應(yīng)力松弛試樣。

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

圖2 松弛試樣示意圖

圖3 DMA實(shí)驗(yàn)系統(tǒng)

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

圖4 遺傳算法計算過程中適應(yīng)度最優(yōu)值和平均值變化情況

表1 松弛實(shí)驗(yàn)數(shù)據(jù)

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

圖5 遺傳算法和傳統(tǒng)算法結(jié)果對比

3 結(jié)論

本文建立了基于遺傳算法的復(fù)合推進(jìn)劑松弛模量獲取方法，為該推進(jìn)劑的粘彈性本構(gòu)模型提供了較精確的本構(gòu)參數(shù)，獲得如下結(jié)論：①由于拉伸試驗(yàn)機(jī)系統(tǒng)存在加載歷程，傳統(tǒng)松弛模量獲取方法獲取的數(shù)據(jù)較真實(shí)值偏小；②基于遺傳算法的松弛模量擬合方法考慮了真實(shí)實(shí)驗(yàn)過程中的加載歷程，因此可以準(zhǔn)確地獲得復(fù)合推進(jìn)劑的松弛模量。

［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］孟紅磊.改性雙基推進(jìn)劑裝藥結(jié)構(gòu)完整性數(shù)值仿真方法研究［D］.南京：南京理工大學(xué)，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ìn)升，鞠玉濤，鄭健.復(fù)合固體推進(jìn)劑松弛模量的獲取方法［J］.火炸藥學(xué)報，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.