基于虛單元法(VEM)的固體推進劑藥柱結構完整性分析①

崔 輝 如

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

0 引言

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

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

本文首次將VEM方法擴展到固體推進劑藥柱的結構完整性分析領域。通過開展粘彈性桿的松弛測試、薄板蠕變測試以及藥柱的點火增壓分析,驗證了VEM在藥柱結構完整性分析領域的可行性。本文所有的仿真計算均在課題組自主研發的發動機結構完整性分析軟件“CHRMULAR”上實現。

1 VEM方法簡介

1.1 二維虛單元空間

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

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

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

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

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

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

1.2 單元剛度矩陣的計算

對于多邊形E,其單元剛度矩陣可以利用拉普拉斯算子進行計算:

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

(1)

定義一個從虛單元空間Vk(E)到單元多項式空間Pk(E)的映射算子Π。該算子對于每一個函數vh∈Vk(E)都有如下正交關系:

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

for allpk∈Pk(E)

(2)

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

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

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

(3)

那么,單元剛度矩陣可以梳理成

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

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

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

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

(4)

利用映射算子的性質,上式中的最后兩項均為0,那么單元剛度矩陣可以進一步表示為

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

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

(5)

進一步地,單元剛度矩陣可通過計算并簡化成

(6)

關于輔助矩陣的數值計算,可以參考文獻[13]。

2 線粘彈性本構模型

考慮泊松比ν為常數,推進劑線粘彈性本構關系表達式為[23]

(7)

(8)

(9)

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

(10)

為了后續計算,這里介紹推進劑蠕變柔量的Prony級數表達式:

(11)

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

3 算例

3.1 單向拉伸粘彈性薄板應力松弛分析

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

圖2給出了薄板的網格以及邊界情況。薄板材料松弛模量以及蠕變柔量的Prony級數參數見表1,泊松比取0.48。

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

表1 推進劑松弛模量以及蠕變柔量Prony級數參數Table 1 Prony series parameters of propellant relaxation modulus and creep compliance

上述問題可以簡化為一維模型進行處理。對于單軸問題,粘彈性的應力-應變關系可以描述為

(12)

(13)

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

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

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

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

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

假設施加的均布載荷函數滿足

(14)

結合邊界條件和本構關系,不難計算出平面應力狀態下薄板兩個方向的應變響應為

(15)

圖5 應變隨時間變化曲線Fig.5 Variation of strain with time

3.3 圓管形藥柱點火增壓分析

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

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

針對上述粘彈性平面應變問題,藥柱結構相應的應力-應變解析解表達式為

(16)

且有

式中σr和σθ分別為徑向與環向應力;εr和εθ分別為徑向與環向應變;r為徑向坐標;t為加壓時間;J∞為平衡蠕變柔量;fJ(t)為蠕變函數。

分別采用商業有限元軟件ABAQUS、VEM以及解析的手段對圓管藥柱受壓問題進行求解。在ABAQUS中,采用1975個四節點四邊形單元,節點總數為2080;VEM中多邊形的總數為1000,節點總數為1994。兩種數值仿真方法中總的自由度數接近。

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

(a)Stress-time curves (b)Strain-time curves圖7 圓管內表面的應力、應變隨時間變化曲線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時刻應力、應變相對誤差沿徑向變化曲線Fig.8 Relative error of stress and strain change along the radial direction at 0.3 s

圖9和圖10分別給出了加載結束時,圓管徑向以及環向的應變云圖,對比云圖分布情況并結合上述結果,VEM手段在計算藥柱結構受壓時具有良好的分析精度。

(a)VEM (b)ABAQUS (a)VEM (b)ABAQUS圖9 圓管徑向應變云圖 圖10 圓管環向應變云圖Fig.9 Radial strain contours of cylinder Fig.10 Hoop strain contours of cylinder

4 結論

本文為固體推進劑藥柱的結構完整性分析提出了一種全新的技術手段。以二維問題為例,介紹了VEM實現的基本流程。針對線粘彈性的推進劑材料,開展了相關算例的仿真驗證。

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