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

張　曉，鄭　堅(jiān)，彭　威，顧志旭

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

?

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

張曉，鄭堅(jiān)，彭威，顧志旭

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

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

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

0　引言

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

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

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

1　測(cè)量原理

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

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

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

2.1迭代模型

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

(1)

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

(2)

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

(3)

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

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

(4)

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

(5)

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

2.2誤差分析模型

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

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

(6)

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

(7)

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

(8)

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

(9)

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

(10)

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

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

KU=F

(11)

由于F與v無關(guān)，方程兩邊對(duì)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　測(cè)試實(shí)例

3.1試驗(yàn)部分

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

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

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

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

(17)

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

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

(18)

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

對(duì)式(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)行仿真計(jì)算，其中泊松比ν即為待反演參數(shù)。

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

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

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

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

表2　不同溫度下的松弛力

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

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

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

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

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

利用時(shí)溫等效原理，得到參考溫度為20 ℃時(shí)泊松比隨時(shí)間的變化關(guān)系，擬合成Prony級(jí)數(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時(shí)的折減時(shí)間，T0為參考溫度，aT為時(shí)溫轉(zhuǎn)換因子。

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

(21)

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

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

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

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

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

(22)

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

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

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

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

圖5　主曲線對(duì)比圖

4　結(jié)論

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

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

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

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

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

[3]蒙上陽,唐國金,雷勇軍. 材料性能對(duì)固體發(fā)動(dòng)機(jī)結(jié)構(gòu)完整性的影響[J]. 國防科技大學(xué)學(xué)報(bào),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)鐵漢,郝松林. 雙激光杠桿/電測(cè)兩用橫向變形引伸計(jì)[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é)報(bào),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]尹銚,劉碧龍,白國鋒,等. 激光測(cè)振有限元反演優(yōu)化方法測(cè)量黏彈材料動(dòng)態(tài)力學(xué)參數(shù)[J]. 聲學(xué)學(xué)報(bào),2013, 38(2):172-180.

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

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

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

[21]趙伯華. 粘彈性泊松比與動(dòng)態(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