有限元分段反演優(yōu)化方法測試固體推進劑的泊松比

張　曉，鄭　堅，彭　威，顧志旭

(軍械工程學院，石家莊　050003)

?

有限元分段反演優(yōu)化方法測試固體推進劑的泊松比

張曉，鄭堅，彭威，顧志旭

(軍械工程學院，石家莊050003)

提出了一種通過準確測量載荷響應，利用ABAQUS有限元軟件，對推進劑相應狀態(tài)下的泊松比進行分段反演的測試方法?；诜囱輧?yōu)化原理，建立了迭代模型和誤差分析模型；測試了某型固體推進劑應力松弛下的泊松比，給出了測試過程和結(jié)果，分析了測量誤差；利用時溫等效原理，獲得了泊松比變化主曲線；反演結(jié)果與光學非接觸測量結(jié)果進行了對比，得到了反演誤差，進一步驗證了測試方法的有效性。

固體推進劑；有限元；分段反演優(yōu)化；泊松比測試；應力松弛；主曲線

0　引言

泊松比研究的關(guān)鍵問題是泊松比的準確測定[1]。固體推進劑為近似不可壓材料，一般取為接近于0.5的定值，但由于藥柱結(jié)構(gòu)性能對泊松比的取值非常敏感[2-3]，因此采用確定性的分析方法，得出的結(jié)果可信性較差[4]。固體推進劑泊松比是溫度和載荷作用時間的函數(shù)，且很早就有人指出，泊松比的測量方法應能準確模擬推進劑的狀態(tài)[5]。

傳統(tǒng)引伸計法操作簡單方便，但由于引伸計造成試件的附加變形，測試結(jié)果很難使人信服；有人建議，采用雙激光杠桿引伸計[6]，測定試件變形，產(chǎn)生附加蠕變效應較小；非接觸式的激光全息照相法和云紋干涉法[7]提高測試精度，但不適用于非線性材料和初始泊松比測量[8]；趙伯華[9]推導了泊松比與其他參數(shù)之間的表達式，間接求解了泊松比，但存在多個參數(shù)不能通過一次試驗同時得到的問題。參數(shù)反演優(yōu)化方法作為獲取未知參量的有效方法，目前廣泛用于眾多工程領(lǐng)域，如地質(zhì)工程[10-11]、土木建筑[12]、復合材料[13-14]以及光譜分析[15-16]等。然而，國內(nèi)外至今尚無有關(guān)參數(shù)反演優(yōu)化用于固體推進劑泊松比測試的報道。

本文針對上述問題進行研究，結(jié)合參數(shù)反演優(yōu)化原理，提出一種通過準確測量載荷響應(位移、應力等)，利用ABAQUS有限元軟件，對推進劑相應狀態(tài)下的泊松比進行分段反演的測試方法；運用該方法測試某型固體推進劑應力松弛下的泊松比，揭示其變化規(guī)律，分析測量誤差；與試驗測量數(shù)據(jù)進行對比，結(jié)果基本一致。本方法實現(xiàn)了由載荷響應反演得到推進劑的材料參數(shù)，操作簡單，測量準確，能進一步推廣到蠕變或單向定速拉伸等條件下的泊松比測試，尤其適用于不易直接測量的場合，對進一步研究推進劑力學性能和藥柱結(jié)構(gòu)完整性具有重要意義。

1　測量原理

分段反演優(yōu)化的基本思想是建立載荷響應與材料參數(shù)之間的理論模型，對載荷響應進行準確測量，利用有限元軟件，將材料參數(shù)隨時間變化的整個過程劃分成多個時間段，分段運用參數(shù)反演優(yōu)化方法進行迭代尋優(yōu)，直到所取參數(shù)使得載荷響應的有限元計算結(jié)果與試驗測量結(jié)果誤差最小，所取參數(shù)值即為真值的最優(yōu)近似。根據(jù)以上分段反演優(yōu)化思想，本文將載荷作用過程進行劃分，假設所劃分段內(nèi)泊松比不變，從而實現(xiàn)對泊松比取值的離散化，運用FORTRAN有限元語言編制用戶子程序，分別進行滿足邊界和載荷條件的循環(huán)迭代運算，尋找能夠使得有限元計算結(jié)果與試驗測量結(jié)果誤差最小的泊松比。此時，所取值即為泊松比真值的最優(yōu)近似，視為此段實際泊松比，當段寬度趨近于零時，最優(yōu)近似值將趨近于瞬時值，整個分段反演優(yōu)化流程如圖1所示。

圖1　分段反演優(yōu)化流程圖

2　數(shù)學模型

2.1迭代模型

加載條件不變的情況下，載荷響應是材料參數(shù)的函數(shù)，即

(1)

松弛模量E(t)為試驗可得，因此有限元軟件計算的結(jié)果僅依賴于泊松比的取值，即

(2)

分析整個加載過程，將測量時間點對應的泊松比 寫成向量形式為

(3)

其中，N為時間段數(shù)，Xi為一個時間段內(nèi)的泊松比向量。

將式(3)代入式(2)，得

(4)

對于任意段，令第j個測點的載荷響應試驗測量值為Hje，有限元計算值為Hjs，m為測點數(shù)，建立段目標函數(shù)為

(5)

令νk+1=νk+Δν進行循環(huán)迭代，Δ2取得最小值時的ν即為泊松比反演的最優(yōu)值。

2.2誤差分析模型

分段反演方法的誤差來源主要有測量誤差、有限元模型誤差、反演誤差等，采用基于靈敏度計算的誤差分析方法[17]，可獲得反演參數(shù)的誤差，具體分析如下。

對于整個加載過程，式(5)改寫為

(6)

令Hjδ=Hje-Hjs為第j個測點總誤差(包括測量誤差、有限元模型誤差、反演誤差等)，將總誤差Hδ寫成向量形式為

(7)

假設νr為泊松比真值，νo為反演最優(yōu)值，令νjδ=νjr-νjo為第j個測點上的泊松比反演誤差，將Hj(νjr)在νjo附近作一階泰勒展開為

(8)

令Hj(νjr)-Hj(νjo)=Hjδ，由式(8)可得

(9)

?H(ν)/?v即為靈敏度矩陣，由式(9)可得

(10)

為了求得泊松比的反演誤差，關(guān)鍵工作是靈敏度矩陣的計算，下面推導有限元靈敏度矩陣計算過程。

有限元方程矩陣形式如下：

KU=F

(11)

由于F與v無關(guān)，方程兩邊對v求導可得

(12)

整理式(12)得

(13)

假設載荷響應函數(shù)H=SU，S為響應函數(shù)與位移的轉(zhuǎn)換矩陣，則由式(13)得

(14)

為了求得?H/?v，關(guān)鍵是要求出?K/?v，由有限元基本原理可知

K=?VBTDBdV

(15)

(16)

K可由有限元程序形成，文中S為應力矩陣，S=DB，因此，由式(14)～式(16)即可求得靈敏度矩陣，從而求出反演誤差。

3　測試實例

3.1試驗部分

試驗的主要內(nèi)容：(1)在70、50、20、0、-20、-40 ℃下，以500mm/min的拉伸速率，分別將6組試件(每組5個)單向等速拉伸至恒定應變5%；(2)記錄應力松弛時間2、4、8、20、40、80、200、600、1 000s時的松弛力值。

采用標準啞鈴型試件，標距70mm，截面10mm×10mm，用微機控制五頭電子式萬能試驗機進行單向應力松弛試驗。

3.2輸入?yún)?shù)計算

對于松弛特性較弱的材料，松弛模量可按式(17)計算：

(17)

式中E(t)為t時刻的應力松弛模量，MPa；F(t)為t時刻的松弛力，N；ε0為初始恒定應變；A0為試件初始橫截面積，mm2。

將松弛模量以Prony級數(shù)形式表示為[18]

(18)

其中，待定系數(shù)E∞、Ei和τi可通過最小二乘法擬合確定，取τi=τ110i-1。

對式(18)進行歸一化處理，可得

(19)

參數(shù)E0、gi、τi為粘彈本構(gòu)模型在ABAQUS中的輸入?yún)?shù)。

3.3有限元模型及仿真

利用ABAQUS建立推進劑試件的有限元模型如圖2所示。

圖2　有限元模型

輸入材料參數(shù)E0、gi、τi、ν進行仿真計算，其中泊松比ν即為待反演參數(shù)。

反演變量設置見表1。

表1　反演變量設置

3.4測試結(jié)果與分析

試驗測得的不同溫度下的松弛力見表2。

表2　不同溫度下的松弛力

試驗按照GJB 770B—2005嚴格執(zhí)行，采用的試驗機精度等級為0.5級，參數(shù)控制準確，系統(tǒng)誤差?。徊僮鬟^程基本由程序控制，且將多組測量結(jié)果取平均值，隨機誤差可忽略不計。因此，認為載荷響應測量準確。將不同溫度下的應力松弛模量擬合成Prony級數(shù)形式，并分別進行歸一化處理，可得有限元計算的輸入?yún)?shù)，如表3所示。

表3　有限元輸入?yún)?shù)

將參數(shù)輸入ABAQUS中，建立有限元模型，使得邊界以及載荷條件與試驗條件相同，利用分段反演優(yōu)化方法對推進劑的泊松比進行測試，結(jié)果如圖3所示。

圖3　泊松比反演結(jié)果

以上結(jié)果表明，隨著溫度的升高和載荷作用時間的增加，泊松比的值增大，且逐漸趨向于0.5。這是由于一方面隨著溫度的升高，推進劑中存在的主價力、次價力、界面之間的化學鍵力和物理吸附作用力等降低，推進劑的應變能力提高[19]；另一方面，隨著加載時間的增長，推進劑的應變表現(xiàn)出延遲性。運用誤差分析模型可得：泊松比反演誤差<5×10-4，主要的誤差來源：一是材料輸入?yún)?shù)的擬合誤差；二是有限元的模型誤差。

利用時溫等效原理，得到參考溫度為20 ℃時泊松比隨時間的變化關(guān)系，擬合成Prony級數(shù)形式為

103ν(t)=498.193 9-1.289 5e-0.0008019ξ-

2.142 3e-0.008019ξ-0.364 1e-0.08019ξ-0.188 2e-0.8019ξ

(20)

其中，ξ=t/aT為由溫度T變到T0時的折減時間，T0為參考溫度，aT為時溫轉(zhuǎn)換因子。

推進劑的時溫轉(zhuǎn)換因子可用WLF方程表示為[20]

(21)

式中T0為參考溫度，K；C1、C2為材料參數(shù)。

文中T0=20 ℃，C1=3，C2=99，轉(zhuǎn)換因子隨溫度的變化曲線，如圖4所示。

同時，采用了光學非接觸方法，根據(jù)粘彈性泊松比的定義，在相同松弛條件下，對推進劑泊松比進行了直接測量。

圖4　轉(zhuǎn)換因子變化曲線

粘彈性泊松比ν(t)表征在拉伸靜載作用下，橫向應變εy(t)對縱向應變εx(t)的響應[21]，可按式(22)計算：

(22)

縱向應變εx(t)可由上述試驗機準確測得；橫向應變εy(t)則由光學位移跟蹤器測定，測試精度可達1 μm，主要誤差來源：一是目標標定誤差；二是因外界光強變化和震動產(chǎn)生的誤差，以上均可通過多組測量取平均值的方法進行消除。

由于測試精度較高，可將多組直接測量結(jié)果的期望值作為泊松比真值，進行反演誤差分析，可得：

依據(jù)時溫等效原理，將試驗測得的泊松比，繪制成參考溫度為20 ℃的主曲線，與分段反演泊松比主曲線進行對比，如圖5所示。

由圖5可知，有限元分段反演測得的泊松比與試驗測量結(jié)果具有較好的一致性，進一步驗證了測試方法的有效性。

圖5　主曲線對比圖

4　結(jié)論

(1)有限元分段反演優(yōu)化方法能準確簡便地模擬固體推進劑的不同受力狀態(tài)，以載荷響應測試試驗為基礎，能測量固體推進劑不同狀態(tài)下的泊松比，結(jié)果具有實用性。

(2)測量了某型固體推進劑應力松弛下的泊松比，結(jié)果表明，隨著溫度的升高和載荷作用時間的增加，泊松比的值增大，且逐漸趨向于0.5，變化規(guī)律與理論分析相吻合。

(3)反演結(jié)果絕對誤差均值為6.4×10-4，與誤差分析模型所得結(jié)果一致，相對誤差均值為0.13%，最大相對誤差為0.2%，說明該方法是測量固體推進劑泊松比的有效途徑。

[1]單桂芳,楊偉,馮建民,等. 材料泊松比測試方法的研究進展[J]. 材料導報,2006,20(3):15-20.

[2]雷勇軍,袁端才,藺文峰. 固體發(fā)動機藥柱正交各向異性材料參數(shù)靈敏度分析[J]. 國防科技大學學報,2010,32(6):107-110.

[3]蒙上陽,唐國金,雷勇軍. 材料性能對固體發(fā)動機結(jié)構(gòu)完整性的影響[J]. 國防科技大學學報,2002,24(5):10-15.

[4]王君祺,職世君,張建偉. 非均布瞬態(tài)內(nèi)壓作用下固體推進劑藥柱泊松比隨機結(jié)構(gòu)分析[J]. 固體火箭技術(shù), 2012,35(6):747-750.

[5]Freudenthal M A. Solid propellant rocket research[M]. London: Academic Press Inc. ,1964.

[6]權(quán)鐵漢,郝松林. 雙激光杠桿/電測兩用橫向變形引伸計[J]. 儀表技術(shù)與傳感器,1999(4):36-38.

[7]O'brien D J,Sottos N R,White S R. Cure-dependent viscoelastic poisson's ratio of epoxy[J]. Experimental Mechanics,2007,47(2):237-249.

[8]郝松林. 高填充復合推進劑的初始粘彈性泊松比[J]. 固體火箭技術(shù),1996,19(4):46-50.

[9]趙伯華. 固體推進劑粘彈泊松比的研究[J]. 北京理工大學學報,1994,14(1):87-90.

[10]Nishimura S I,Nishiyama T,Murakami A. Inverse analysis of soft grounds considering nonlinearity and anisotropy[J]. Journal of the Japanese Geotechnical Society, 2005, 45(2):87-95.

[11]Hashash Y M A,Levasseur S,Osouli A,et al. Comparison of two inverse analysis techniques for learning deep excavation response[J]. Computers and Geotechnics, 2010, 37(3):323-333.

[12]黃東梅,朱樂東,丁泉順,等. 超高層建筑等效靜力風荷載的反演法[J]. 工程力學,2012,29(1):99-105.

[13]Anghileri M,Chirwa E C,Lanzi L,et al. An inverse approach to identify the constitutive model parameters for crashworthiness modelling of composite structure[J]. Composite Structures,2005,68(1):65-74.

[14]Kang Y L,Lin X H,Qin Q H. Inverse/genetic method and its application in identification of mechanical parameters of interface in composite[J]. Composite Structures,2004, 66(1):449-458.

[15]Van Beek J D,Meier B H,SchaFer H. Inverse methods in two-dimensional NMR spectral analysis[J]. Journal of Magnetic Resonance,2003,162(1):141-157.

[16]Yurko V. Inverse spectral problems for differential pencils on the half-line with turning points[J]. Journal of Mathematical Analysis and Applications,2006,320(1):439-463.

[17]尹銚,劉碧龍,白國鋒,等. 激光測振有限元反演優(yōu)化方法測量黏彈材料動態(tài)力學參數(shù)[J]. 聲學學報,2013, 38(2):172-180.

[18]常武軍,鞠玉濤,王蓬勃. HTPB推進劑脫濕與力學性能的相關(guān)性研究[J]. 兵工學報,2012,33(3):261-266.

[19]侯林法. 復合固體推進劑[M]. 北京: 宇航出版社, 1994.

[20]周光泉,劉孝敏. 粘彈性理論[M]. 合肥: 中國科學技術(shù)大學出版社,1996.

[21]趙伯華. 粘彈性泊松比與動態(tài)復數(shù)泊松比的研究[J]. 推進技術(shù),1995,16(3):1-7.

(編輯：薛永利)

Measurement of solid propellant Poisson's ratio with the method of finite element piecewise inversion optimization

ZHANG Xiao, ZHENG Jian, PENG Wei, GU Zhi-xu

(Ordnance Engineering College, Shijiazhuang050003, China)

A novel method named finite element piecewise inversion optimization was developed for solid propellant Poisson's ratio measurement, which depended on the accurate load response test and the application of the finite element software ABAQUS. Based on the inversion optimization theory, iterative model and error analysis model were built. By means of piecewise inversion, the Poisson's ratio of a certain type of solid propellant was measured in the condition of stress relaxation, the process and results were provided with error analysis as well. Using the time temperature equivalence principle, the Poisson's ratio master curve was obtained. At the last, the inversion results were compared with optical non-contact measurement results while the errors were got, which further verified the effectiveness of this method.

solid propellant；finite element；piecewise inversion optimization；Poisson's ratio measurement；stress relaxation；master curve

2015-01-20；

2015-05-23。

張曉(1990—)，男，碩士，研究方向為固體推進劑力學性能。E-mail:erebuss@outlook.com

V512

A

1006-2793(2016)04-0519-05

10.7673/j.issn.1006-2793.2016.04.013