基于虛單元法(VEM)的固體推進(jìn)劑藥柱結(jié)構(gòu)完整性分析①

崔 輝 如

(1.陸軍工程大學(xué) 國(guó)防工程學(xué)院,南京 210007;2.陸軍工程大學(xué) 土木工程博士后科研流動(dòng)站,南京 210007)

0 引言

固體發(fā)動(dòng)機(jī)由于結(jié)構(gòu)形式簡(jiǎn)單,使用方便快捷被廣泛用作火箭彈、導(dǎo)彈和探空火箭等的動(dòng)力裝置。藥柱結(jié)構(gòu)的完整性分析是發(fā)動(dòng)機(jī)安全點(diǎn)火的重要前提。傳統(tǒng)的解析方法只能用來處理簡(jiǎn)單的藥柱結(jié)構(gòu)。試驗(yàn)方法,比如采用模擬發(fā)動(dòng)機(jī)點(diǎn)火以及冷增壓裝置等,是最能直觀判斷藥柱結(jié)構(gòu)完整性的手段,但是高昂的試驗(yàn)成本以及有限的可測(cè)力學(xué)性能參數(shù)使其不太可能大規(guī)模實(shí)現(xiàn)。數(shù)值仿真方面,國(guó)內(nèi)的大多數(shù)院校和科研機(jī)構(gòu)目前已經(jīng)可以熟練地運(yùn)用有限元手段進(jìn)行自由裝填[1]、貼壁澆注式[2]、考慮圍壓效應(yīng)[3]以及考慮界面損傷[4-5]的藥柱結(jié)構(gòu)完整性分析。此外,有研究人員還針對(duì)推進(jìn)劑的損傷[6]、老化[7]以及粘彈性泊松比特性[8],開展了一些列藥柱結(jié)構(gòu)相關(guān)的力學(xué)響應(yīng)研究。在利用有限元進(jìn)行藥柱結(jié)構(gòu)完整性分析過程中,幾何簡(jiǎn)化以及網(wǎng)格劃分處理幾乎占據(jù)了絕大部分時(shí)間[9]。另一方面,在運(yùn)輸及戰(zhàn)備服役階段,藥柱內(nèi)部會(huì)發(fā)生裂紋等缺陷,而有限元方法在處理裂紋擴(kuò)展問題時(shí)面臨的最大挑戰(zhàn)就是裂紋擴(kuò)展后單元的重構(gòu)問題[10]。因此,有必要在保證計(jì)算精度的前提下,尋求一種對(duì)網(wǎng)格形式依賴性較低的藥柱結(jié)構(gòu)完整性分析方法。

虛單元法(Virtual element method,VEM)是近10年提出來的一種適用于多邊形單元的數(shù)值仿真方法[11-13],最顯著的特點(diǎn)就是可以不需要形函數(shù)的顯式表達(dá)式,在伽遼金框架下進(jìn)行雙線性以及線性式的計(jì)算。正是在處理多邊形時(shí)不需要設(shè)計(jì)復(fù)雜的形函數(shù)表達(dá)式,VEM被廣泛地用于小變形彈性[14,15]、有限變形[16]以及彈塑性[17-19]問題的求解。國(guó)內(nèi)也有相關(guān)研究人員利用VEM開展結(jié)構(gòu)仿真工作。林姍等[20]在VEM的基礎(chǔ)上,開展了彈塑性力學(xué)問題的應(yīng)用研究,與傳統(tǒng)的有限元方法相比,VEM計(jì)算準(zhǔn)確高效,網(wǎng)格處理靈活。江巍等[21]使用VEM構(gòu)建單個(gè)塊體的位移,發(fā)展了基于VEM的非連續(xù)變形分析方法的新格式。劉傳奇等[22]對(duì)VEM的發(fā)展、優(yōu)勢(shì)、應(yīng)用等方面開展了全面的綜述。盡管VEM還處于起步階段,但是該方法在諸如非線性問題、接觸問題以及裂紋擴(kuò)展問題上展現(xiàn)出巨大的潛力。目前,國(guó)內(nèi)外尚未有利用VEM進(jìn)行固體推進(jìn)劑藥柱結(jié)構(gòu)完整性分析研究的報(bào)道。

本文首次將VEM方法擴(kuò)展到固體推進(jìn)劑藥柱的結(jié)構(gòu)完整性分析領(lǐng)域。通過開展粘彈性桿的松弛測(cè)試、薄板蠕變測(cè)試以及藥柱的點(diǎn)火增壓分析,驗(yàn)證了VEM在藥柱結(jié)構(gòu)完整性分析領(lǐng)域的可行性。本文所有的仿真計(jì)算均在課題組自主研發(fā)的發(fā)動(dòng)機(jī)結(jié)構(gòu)完整性分析軟件“CHRMULAR”上實(shí)現(xiàn)。

1 VEM方法簡(jiǎn)介

1.1 二維虛單元空間

假設(shè)空間Ω?R2被分割成一系列不重合的多邊形E的集合Ph,這里的多邊形不一定需要凸多邊形。對(duì)于任意的多邊形E,按照逆時(shí)針方向標(biāo)記其頂點(diǎn)Vi,分別為Vi(i=1,…,NV)。類似地,標(biāo)記多邊形的邊ei,分別為ei(i=1,…,NV),如圖1。

圖1 典型多邊形EFig.1 Typical polygons E

對(duì)于任意的多邊形E,定義一個(gè)局部虛單元空間Vk(E),該空間包含了所有的k次多項(xiàng)式以及其他一些函數(shù),這些函數(shù)在單元邊的約束同樣也是k次多項(xiàng)式。定義具有以下屬性的函數(shù)vh∈Vk(E):

(1)vh在單元E的任意邊e上是k次多項(xiàng)式;

(2)vh在?E上是全局連續(xù)的;

(3)Δvh在E中是k-2次多項(xiàng)式。

1.2 單元?jiǎng)偠染仃嚨挠?jì)算

對(duì)于多邊形E,其單元?jiǎng)偠染仃嚳梢岳美绽顾阕舆M(jìn)行計(jì)算:

(KE)ij=(φi,φj)0,E,i,j=1,…,Ndof

(1)

定義一個(gè)從虛單元空間Vk(E)到單元多項(xiàng)式空間Pk(E)的映射算子Π。該算子對(duì)于每一個(gè)函數(shù)vh∈Vk(E)都有如下正交關(guān)系:

(pk,(Πvh-vh))0,E=0

for allpk∈Pk(E)

(2)

式中Pk(E)為多邊形E上階次小于等于k的多項(xiàng)式空間;pk屬于多項(xiàng)式空間Pk(E)的一個(gè)子集。

利用上述算子,可以將基函數(shù)φi整理為

φi=Πφi+(1-Π)φi

(3)

那么,單元?jiǎng)偠染仃嚳梢允崂沓?/p>

(KE)ij=(Πφi,Πφj)0,E+

((1-Π)φi,(1-Π)φj)0,E+

(Πφi,(1-Π)φj)0,E+

((1-Π)φi,Πφj)0,E

(4)

利用映射算子的性質(zhì),上式中的最后兩項(xiàng)均為0,那么單元?jiǎng)偠染仃嚳梢赃M(jìn)一步表示為

(KE)ij=(Πφj)0,E+

((1-Π)φi,(1-Π)φj)0,E

(5)

進(jìn)一步地,單元?jiǎng)偠染仃嚳赏ㄟ^計(jì)算并簡(jiǎn)化成

(6)

關(guān)于輔助矩陣的數(shù)值計(jì)算,可以參考文獻(xiàn)[13]。

2 線粘彈性本構(gòu)模型

考慮泊松比ν為常數(shù),推進(jìn)劑線粘彈性本構(gòu)關(guān)系表達(dá)式為[23]

(7)

(8)

(9)

其中,σij、Sij、σkk分別為應(yīng)力張量、偏應(yīng)力張量以及球應(yīng)力張量;eij、εkk分別為偏應(yīng)變張量及球應(yīng)變張量;t為加載時(shí)間;E(t)為松弛模量,其Prony級(jí)數(shù)表達(dá)式為

(10)

為了后續(xù)計(jì)算,這里介紹推進(jìn)劑蠕變?nèi)崃康腜rony級(jí)數(shù)表達(dá)式:

(11)

VEM處理線粘彈性材料的步驟與處理線彈性方法類似,讀者可參考文獻(xiàn)[13]。唯一和處理線彈性問題不同的是,VEM需要結(jié)合牛頓-拉夫遜方法才能全面展示出材料的粘彈特性,比如蠕變效應(yīng)。此外,關(guān)于線粘彈性本構(gòu)關(guān)系的離散與數(shù)值積分方法,可以參考文獻(xiàn)[7,24]。

3 算例

3.1 單向拉伸粘彈性薄板應(yīng)力松弛分析

考慮一長(zhǎng)50 mm,寬10 mm的單位厚度長(zhǎng)方形粘彈性薄板。約束薄板左側(cè)的橫向自由度并在右側(cè)施加階躍位移,由此模擬薄板在松弛階段的粘彈性響應(yīng)。

圖2給出了薄板的網(wǎng)格以及邊界情況。薄板材料松弛模量以及蠕變?nèi)崃康腜rony級(jí)數(shù)參數(shù)見表1,泊松比取0.48。

圖2 粘彈性薄板多邊形網(wǎng)格及其邊界Fig.2 Polygonal mesh and its boundaries of viscoelastic sheet

表1 推進(jìn)劑松弛模量以及蠕變?nèi)崃縋rony級(jí)數(shù)參數(shù)Table 1 Prony series parameters of propellant relaxation modulus and creep compliance

上述問題可以簡(jiǎn)化為一維模型進(jìn)行處理。對(duì)于單軸問題,粘彈性的應(yīng)力-應(yīng)變關(guān)系可以描述為

(12)

(13)

如圖3所示,可以看到不同階躍位移水平(3.75、7.5、15 mm)下,薄板軸向應(yīng)力的數(shù)值解與解析解是一致的。上述仿真結(jié)果表明,VEM可以準(zhǔn)確地模擬粘彈性結(jié)構(gòu)的應(yīng)力松弛效應(yīng)。

圖3 不同階躍位移下薄板的軸向應(yīng)力響應(yīng)Fig.3 Axial stress response of thin plates at different step displacements

3.2 雙向拉伸粘彈性薄板蠕變分析

進(jìn)一步地,為檢驗(yàn)VEM模擬推進(jìn)劑蠕變特性的能力,設(shè)計(jì)了以下雙軸拉伸蠕變變形試驗(yàn),如圖4。該薄板長(zhǎng)100 mm,寬50 mm,厚度為1 mm。在進(jìn)行蠕變?cè)囼?yàn)時(shí),約束薄板的左側(cè)橫向位移以及下側(cè)的豎向位移,并在右側(cè)及上側(cè)表面施加大小不同的均布法向載荷。薄板的材料參數(shù)與3.1節(jié)一致。

圖4 粘彈性薄板多邊形網(wǎng)格及其邊界Fig.4 Polygonal mesh of viscoelastic sheet and its boundary

假設(shè)施加的均布載荷函數(shù)滿足

(14)

結(jié)合邊界條件和本構(gòu)關(guān)系,不難計(jì)算出平面應(yīng)力狀態(tài)下薄板兩個(gè)方向的應(yīng)變響應(yīng)為

(15)

圖5 應(yīng)變隨時(shí)間變化曲線Fig.5 Variation of strain with time

3.3 圓管形藥柱點(diǎn)火增壓分析

圖6給出了一個(gè)圓管藥柱的物理模型,藥柱的內(nèi)徑a=200 mm,外徑b=400 mm。藥柱外側(cè)受到固定約束,內(nèi)表面受內(nèi)壓p(t)=p0(1-e-kt)的作用。其中p0為1.0 MPa,壓力因子k取25,作用時(shí)間為0.3 s。藥柱的材料參數(shù)與上文一致。

圖6 圓管藥柱物理模型Fig.6 Physical model of tube grain

針對(duì)上述粘彈性平面應(yīng)變問題,藥柱結(jié)構(gòu)相應(yīng)的應(yīng)力-應(yīng)變解析解表達(dá)式為

(16)

且有

式中σr和σθ分別為徑向與環(huán)向應(yīng)力;εr和εθ分別為徑向與環(huán)向應(yīng)變;r為徑向坐標(biāo);t為加壓時(shí)間;J∞為平衡蠕變?nèi)崃?fJ(t)為蠕變函數(shù)。

分別采用商業(yè)有限元軟件ABAQUS、VEM以及解析的手段對(duì)圓管藥柱受壓?jiǎn)栴}進(jìn)行求解。在ABAQUS中,采用1975個(gè)四節(jié)點(diǎn)四邊形單元,節(jié)點(diǎn)總數(shù)為2080;VEM中多邊形的總數(shù)為1000,節(jié)點(diǎn)總數(shù)為1994。兩種數(shù)值仿真方法中總的自由度數(shù)接近。

圖7給出了藥柱內(nèi)表面徑向以及環(huán)向應(yīng)力隨時(shí)間變化的曲線,可以看到VEM方法計(jì)算出的結(jié)果與解析法以及ABAQUS方法得到的結(jié)果分布一致。另一方面,圖8給出了加載至0.3 s時(shí)刻,應(yīng)力、應(yīng)變相對(duì)誤差沿徑向的路徑分布,不難看出,兩種仿真方法的計(jì)算結(jié)果與解析解相近。

(a)Stress-time curves (b)Strain-time curves圖7 圓管內(nèi)表面的應(yīng)力、應(yīng)變隨時(shí)間變化曲線Fig.7 Variation of stress and strain on the inner surface of cylinder with time

(a)Relative error of stress curves (b)Relative error of strain curves圖8 加載至0.3 s時(shí)刻應(yīng)力、應(yīng)變相對(duì)誤差沿徑向變化曲線Fig.8 Relative error of stress and strain change along the radial direction at 0.3 s

圖9和圖10分別給出了加載結(jié)束時(shí),圓管徑向以及環(huán)向的應(yīng)變?cè)茍D,對(duì)比云圖分布情況并結(jié)合上述結(jié)果,VEM手段在計(jì)算藥柱結(jié)構(gòu)受壓時(shí)具有良好的分析精度。

(a)VEM (b)ABAQUS (a)VEM (b)ABAQUS圖9 圓管徑向應(yīng)變?cè)茍D 圖10 圓管環(huán)向應(yīng)變?cè)茍DFig.9 Radial strain contours of cylinder Fig.10 Hoop strain contours of cylinder

4 結(jié)論

本文為固體推進(jìn)劑藥柱的結(jié)構(gòu)完整性分析提出了一種全新的技術(shù)手段。以二維問題為例,介紹了VEM實(shí)現(xiàn)的基本流程。針對(duì)線粘彈性的推進(jìn)劑材料,開展了相關(guān)算例的仿真驗(yàn)證。

計(jì)算結(jié)果表明,VEM可準(zhǔn)確模擬出粘彈性材料的松弛和蠕變響應(yīng)。在點(diǎn)火增壓條件下,VEM計(jì)算出的應(yīng)力-應(yīng)變響應(yīng)與ABAQUS以及解析解吻合一致。從二維分析結(jié)果來看,VEM可以準(zhǔn)確預(yù)示粘彈性結(jié)構(gòu)在各種載荷下的應(yīng)力-應(yīng)變響應(yīng)。結(jié)合VEM對(duì)網(wǎng)格形式的低依賴性,利用VEM開展復(fù)雜藥柱結(jié)構(gòu)完整性分析工作將會(huì)大幅提高網(wǎng)格前處理的效率。此外,VEM有望進(jìn)一步應(yīng)用于藥柱的三維結(jié)構(gòu)完整性分析以及推進(jìn)劑裂紋擴(kuò)展等領(lǐng)域。