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

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

?

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

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

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

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

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

0　引言

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

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

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

1　測量原理

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

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

2　數(shù)學(xué)模型

2.1迭代模型

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

(1)

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

(2)

分析整個加載過程，將測量時間點(diǎn)對應(yīng)的泊松比 寫成向量形式為

(3)

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

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

(4)

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

(5)

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

2.2誤差分析模型

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

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

(6)

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

(7)

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

(8)

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

(9)

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

(10)

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

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

KU=F

(11)

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

(12)

整理式(12)得

(13)

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

(14)

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

K=?VBTDBdV

(15)

(16)

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

3　測試實(shí)例

3.1試驗(yàn)部分

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

采用標(biāo)準(zhǔn)啞鈴型試件，標(biāo)距70mm，截面10mm×10mm，用微機(jī)控制五頭電子式萬能試驗(yàn)機(jī)進(jìn)行單向應(yīng)力松弛試驗(yàn)。

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

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

(17)

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

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

(18)

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

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

(19)

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

3.3有限元模型及仿真

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

圖2　有限元模型

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

反演變量設(shè)置見表1。

表1　反演變量設(shè)置

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

試驗(yàn)測得的不同溫度下的松弛力見表2。

表2　不同溫度下的松弛力

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

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

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

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

以上結(jié)果表明，隨著溫度的升高和載荷作用時間的增加，泊松比的值增大，且逐漸趨向于0.5。這是由于一方面隨著溫度的升高，推進(jìn)劑中存在的主價力、次價力、界面之間的化學(xué)鍵力和物理吸附作用力等降低，推進(jìn)劑的應(yīng)變能力提高[19]；另一方面，隨著加載時間的增長，推進(jìn)劑的應(yīng)變表現(xiàn)出延遲性。運(yù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)換因子。

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

(21)

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

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

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

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

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

(22)

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

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

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

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

圖5　主曲線對比圖

4　結(jié)論

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

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

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

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

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

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

[4]王君祺,職世君,張建偉. 非均布瞬態(tài)內(nèi)壓作用下固體推進(jìn)劑藥柱泊松比隨機(jī)結(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]郝松林. 高填充復(fù)合推進(jìn)劑的初始粘彈性泊松比[J]. 固體火箭技術(shù),1996,19(4):46-50.

[9]趙伯華. 固體推進(jìn)劑粘彈泊松比的研究[J]. 北京理工大學(xué)學(xué)報,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]黃東梅,朱樂東,丁泉順,等. 超高層建筑等效靜力風(fēng)荷載的反演法[J]. 工程力學(xué),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)力學(xué)參數(shù)[J]. 聲學(xué)學(xué)報,2013, 38(2):172-180.

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

[19]侯林法. 復(fù)合固體推進(jìn)劑[M]. 北京: 宇航出版社, 1994.

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

[21]趙伯華. 粘彈性泊松比與動態(tài)復(fù)數(shù)泊松比的研究[J]. 推進(jìn)技術(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—)，男，碩士，研究方向?yàn)楣腆w推進(jìn)劑力學(xué)性能。E-mail:erebuss@outlook.com

V512

A

1006-2793(2016)04-0519-05

10.7673/j.issn.1006-2793.2016.04.013