基于三維八節點內聚單元的固體火箭發動機粘接界面結構分析①

崔輝如，王佳奇，呂 軒，許玉榮

(1.陸軍工程大學 國防工程學院, 南京 210007；2.中國酒泉衛星發射中心，酒泉 732750；3.航天動力先進技術湖北省重點實驗室，湖北航天技術研究院總體設計所，武漢 430040)

0 引言

固體火箭發動機作為性能優越的動力裝置被廣泛地應用于軍事和航天等領域。高體分比、服役周期長以及豎直貯存等一系列的嚴格指標使得發動機粘接界面面臨著前所未有的考驗。

在粘接界面的監檢測方面，不少科研人員開始利用紅外熱掃描和成像法進行絕熱層/推進劑粘接界面的檢測。海軍航空大學在國內率先開展了固體火箭發動機界面實時監測系統設計與試驗相關工作。通過將傳感器埋入粘接界面處，并進行了長時間的跟蹤監測。監測結果表明，傳感器對界面的粘接性能影響較小，監測結果準確可靠。類似地，美國的空軍實驗室在絕熱層與推進劑之間內置一定數量的傳感器實時監測界面處的徑向應力以及溫度。

傳統的仿真分析過程中，將粘接界面等效處理的方法已經難以滿足現階段界面破壞過程模擬的需求。在BARENBLATT和DUGDALE提出內聚力模型的概念之后，內聚力模型儼然成為描述脫粘過程強大的分析和數值計算工具。不少研究人員針對發動機粘接脫粘問題，利用內聚力模型開展了一系列的研究。CUI等通過建立松弛內聚力和松弛參數的方法，分別創建了兩種粘彈性內聚力本構模型。仿真結果表明，這兩類模型可以很好地應用于推進劑和絕熱層之間粘接界面的分離破壞。馬曉琳等通過內聚單元構建了襯層和推進劑粘接界面的二維有限元模型，通過反演手段獲得內聚力模型參數，研究了發動機在典型工況下，推進劑和絕熱層粘接界面的力學響應。南京理工大學團隊采用內聚單元和雙線性內聚力模型，在有限元方法的基礎上分別開展了推進劑和絕熱層粘接界面的I型、II型界面分離特性研究，并取得了一系列的成果。

綜上可知，目前利用內聚力模型進行的界面力學行為模擬主要還集中在試驗件層面。為了滿足發動機粘接界面的結構完整性分析，本文從界面的力學模型入手，構建了三維八節點內聚單元，推導了單元內力矢量以及單元剛度陣，研究和建立有效的有限元仿真方法，并對某固體火箭發動機固化降溫和點火增壓兩種工況條件下的粘接界面結構進行了三維有限元分析。

1 基于松弛效應的粘彈性內聚力本構模型

1.1 PPR內聚力本構模型

PPR內聚力模型是考慮物理斷裂參數以及連續斷裂邊界的勢相關內聚力模型。利用勢函數與內聚力之間的關系，可以得到兩個方向的內聚力表達式。

(1)

(2)

式中和分別代表法向和切向的分離位移；和分別代表切向和法向內聚能；和為能量常數；和為初始斜率參數；和為模型的形狀參數，分別控制法向和切向的內聚力和分離位移曲線下降段的凹凸特性；和為模型的法向和切向的臨界位移。

和的計算與兩個方向的內聚能相關。

(1)當≠時

(3)

(2)當=時

(4)

此外

(5)

式中和為法向和切向內聚強度；和為初始斜率指針。

(6)

PPR內聚力模型的控制參數可以分成4組，即初始斜率參數、，內聚強度、，臨界位移、以及形狀參數、。

1.2 粘彈性內聚力模型

大部分的情況下，襯層內的粘接劑和固體推進劑的粘接劑采用相同的高聚物配方，是典型的粘彈性材料。研究表明，粘接界面具有典型的粘彈性效應。為此，CUI等提出了粘接界面內聚參數隨時間松弛的假設：

(1)時間相關內聚力的內在機制是內聚力模型參數的松弛，這些內聚力模型參數包括臨界位移、內聚強度以及初始斜率參數。

(2)時間相關的內聚力依然是內聚力模型中勢函數對分離位移的偏導數，這一條件在任意時刻都滿足。

(3)加載時間以及加載溫度的力學效應可以依據時溫等效原理進行相互轉化。

根據上面的基本假設，這里采用形如()=+exp(-)的標準線性固體松弛模型。為了使得系數在初始時刻為1，那么+=1。初始的內聚力模型參數乘以松弛系數就可以獲得松弛的內聚力模型參數。式(7)以I型內聚參數為例給出了松弛的臨界位移、內聚強度以及初始斜率參數的時變關系：

(7)

式中(0<<1)、(0<<1)、(0<<1)、、以及為未知的松弛參數；(0)、(0)以及(0)為=0時刻的內聚強度、臨界位移以及初始斜率參數。

考慮到加載溫度的影響，依據時溫等效原理獲得的縮減時間將會被引入到以上定義之中，即

(8)

其中

(9)

(10)

1.3 縮減時間的數值積分

本部分將具體介紹如何在增量步中計算松弛型的內聚力模型參數，主要的工作將集中在縮減時間的計算上面。接下來的部分將以內聚強度為例進行介紹。

如式(11)，縮減時間是溫度的非線性函數。通過兩邊進行求導，可得

(11)

為了計算縮減時間，這里將采用增量的方法。將計算時間區間分解成份時間子區間，即

(12)

令

(13)

假設溫度()以及()在時間區間[,+1]上線性變化，那么

()=+

(14)

其中

(15)

結合上式，可得

(16)

依據式(16)，縮減時間的增量可以表示為

(17)

計算過程中，時間增量Δ+1是已知量，與之對應的臨界位移，內聚強度以及初始斜率參數的縮減時間可以采用以上方法獲得。

2 三維八節點內聚單元

為了實現對粘彈性內聚力本構模型的應用，還需要開發相應的內聚單元。目前，為了確保有限元分析精度，三維固體火箭發動機結構完整性分析大都采用六面體單元。相應地，內聚單元主要采用三維八節點的形式。為了簡化問題，這里僅討論初始厚度為0的三維八節點內聚單元，如圖1所示。

圖1 三維八節點零厚度內聚單元應用示意圖Fig.1 Application diagram of 3D 8-node zero-thickness cohesive element

2.1 單元位移變換矩陣推導

在進行三維八節點內聚單元的推導之前，需要對單元的編號規則進行簡單的介紹。如圖2所示，單元的下表面的節點編號按照逆時針的規則分別為節點1、2、3、4。與之對應的單元的上表面的節點編號按照逆時針的規則分別為節點5、6、7、8。首先提取三維內聚單元的中面，即平面1′2′3′4′。

如圖3所示，規定為單元所處的整體坐標系，為單元參考的坐標系。假設中面各頂點的坐標在整體坐標系和參考坐標中分別表示為

(18)

依據中面的定義，頂點1′、2′、3′、4′的整體坐標可表示為式(19)～式(22)。

圖2 三維八節點內聚單元及其中面示意圖Fig.2 Diagram of 3D 8-node cohesive element and its middle surface

圖3 三維八節點內聚單元參考坐標系示意圖Fig.3 Reference frame scheme for 3D 8-node cohesive element

(19)

(20)

(21)

(22)

參考坐標系的原點可以選擇在單元內的任意一點，以1′點為例，即

=

(23)

在單元平面內選擇和軸，那么軸自然垂直于單元面，按照角點1′→2′→3′右螺旋指向的正方向。

(24)

(25)

那么，軸的方向余弦為

(26)

其中

(27)

(28)

(29)

其中

(30)

(31)

此時，方向的余弦可以由軸和軸確定，根據右螺旋法則有

(32)

坐標變換矩陣為

(33)

以上圖為例，參考坐標系內的節點位移向量為

(34)

總體坐標系內的節點位移向量為

=[]，(=1,2,…,8)

(35)

節點的位移向量滿足

(36)

為方便，令

=[]

(37)

(38)

單元的節點位移向量滿足

(39)

(40)

式中代表三階單位矩陣。

2.2 單元分離位移場推導

已知內聚單元變形后節點的參考坐標表示為

(41)

那么，參考坐標系下的分離位移為

=[]，(=1,2,3,4)

(42)

(43)

(44)

考慮到零厚度單元：

(45)

(46)

=[]

(47)

建立節點分離位移與節點參考位移間的關系

(48)

其中

=[-]

(49)

建立分離位移場與節點分離位移關系：

(50)

其中

(51)

(52)

建立分離位移場與節點整體位移間的關系，即

(53)

2.3 八節點六面體等參單元

類似四邊形的等參元，可以建立八節點六面體等參單元。對于如圖4所示的等參單元，其形函數在自然坐標系中如式(54)所示。

圖4 自然坐標系中的八節點六面體等參單元Fig.4 8-node hexahedral isoparametric element in natural coordinate system

(54)

類似地，可以通過形函數，建立自然坐標系中的點與圖5中物理坐標系中的點之間的映射關系，即

(55)

圖5 物理坐標系中的八節點六面體單元Fig.5 8-node hexahedron element in physical coordinate system

對于單元上的面積積分而言，令

(56)

(57)

利用自然坐標系有

(58)

中面上(=0)的面積積分為

(59)

(60)

(61)

故

(62)

2.4 單元內力矢量及單元剛度陣推導

依據虛位移原理，有

(63)

式中和分別為格林應變張量和第二PK應力張量；和分別為分離位移及內聚力；為單元中面；、分別為外力及節點位移矢量；為整個單元外邊界。

通過伽遼金法結合等參元理論，推導可得單元的剛度矩陣及內力矢量為：

(64)

(65)

基于以上推導并結合ABAQUS中的用戶單元模塊，將本文開發的內聚單元進行有限元應用。

2.5 單元驗證算例

驗證算例的主要對象為彈性圓管粘接結構，如圖6所示。圓管的內外徑分別為8 mm和2 mm。在圓管靠近外表面的位置，沿著圓管的徑向，設計了長度為3 mm的預制裂紋。圓管的內部受壓，壓力值為1 MPa，圓管結構的底部被固定約束。沿著預制裂紋的方向是后加工而成的粘接界面，結構受內壓過程中，裂紋將沿著粘接界面的方向擴展。材料的楊氏模量以及泊松比分別為10 MPa和0.33。PPR模型的初始參數滿足(0)=0.5 MPa，(0)=0.1672，=3.10以及(0)=2.5 mm。另外，切向的內聚參數假設和法向是一致的。

圖6 圓管結構幾何模型(單位：mm)Fig.6 Geometrical model of tubular structure (unit:mm)

為了方便分析粘彈性模型的特性，這里定義模型的參數滿足

====03

====50

(66)

網格方面，圖7針對圓管結構本文分別采用了二維四節點單元(Q4)和三維八節點單元(Hex8)。網格主要的尺寸為0.5 mm，三維網格結構的整體厚度為2 mm。這里的二維單元均采用平面應力單元。網格規模上，六面體網格的數量是四邊形網格數量的4倍。在結構的粘接界面處分別設置了經過文獻[31]檢驗的二維四節點內聚單元(COH2D4)和本文設計的三維八節點內聚單元(COH3D8)。內聚單元數量方面，六面體網格中的內力單元是四邊形網格中內聚單元數量的16倍。

(a)Quadrilateral mesh

(b)Hexahedral mesh圖7 圓管結構網格Fig.7 Mesh for tubular structure

圖8給出的是不考慮溫度變化條件下，內壓在8 s內從0 MPa勻速增加到0.8 MPa時，圓管結構頂部裂縫寬度隨時間的變化。可以發現，不同網格在平面應力下的計算結果和三維應力狀態下的計算結果是對應一致的。通過以上分析，可以認為，本文提出的內聚單元推導過程和仿真方法是準確的。

圖8 裂縫寬度隨時間分布曲線Fig.8 Crack width versus time

3 粘接界面結構完整性分析

3.1 發動機有限元模型

表1 推進劑松弛模量Prony級數系數

表2 粘彈性粘接界面內聚力模型參數

如圖9所示，在發動機推進劑和絕熱層之間對網格進行加密，盡量將內聚單元附近的網格質量對結構分析的影響降到最低。在提交計算時，根據對稱特性，對結構的環向位移進行約束。為了抑制剛體位移，對發動機殼體的一端進行約束。

(a)Overall schematic diagram

(b)Partial grid diagram between propellant and insulation(case,insulation and propellant from outside to inside)圖9 發動機有限元網格示意圖Fig.9 Finite element mesh diagram of solid rocket motor

3.2 固化降溫分析

藥柱在零應力溫度58 ℃(331.15 K)澆注成型，還需要經過長達1 d(86 400 s)的固化降溫階段，直到溫度降到室溫20 ℃(293.15 K)。在這一階段，假設發動機結構內部是沒有溫度差的。發動機內部的溫度變化曲線可以近似表示為

()=-(-)×(1-e- ×10)

(67)

其中，為初始固化溫度；為室溫；參數=8 s；時間以秒計。

固化降溫的初期，溫度變化較快，在降溫的后半段時期，溫度變化較慢。固化降溫結束時，推進劑和絕熱層粘接界面上的內聚力分布云圖如圖10所示。可以發現，在粘接界面的兩端，單元的切向和法向內聚力急劇增加。在降溫階段，藥柱熱脹系數較大，藥柱整體向中心收縮。因此，藥柱兩端的變形劇增，這就導致了粘接界面在兩端的內聚力較中間部位稍大。從發動機的縱剖面，更加直觀地表現出了這一趨勢。

此外，從圖10(a)可以看到，切向內聚力(這里的切向內聚力是發動機縱向和環向內聚力的合力)在發動機的縱向會出現明顯的極小值，這一極小值在圓管段區域。法向內聚力方面，如圖10(b)所示，由于推進劑沿著法向的裝填量有一定差異，降溫時收縮變形量不同，法向內聚力在星孔段以及翼槽段會出現一定的波動，而在圓管段，法向內聚力基本保持不變。

在固化降溫結束時，粘接界面上最大的切向和法向內聚力分別為0.03 MPa和0.065 MPa，遠小于粘接界面的內聚強度，此階段粘接界面結構破壞的可能性較低。

(a)Tangential direction (b)Normal direction圖10 固化降溫結束粘接界面內聚力空間分布云圖Fig.10 Spatial distribution of traction on the bonding interface at the end of cooling stage

3.3 點火增壓分析

發動機的點火增壓階段是發動機結構設計中另一個較為危險的階段。為了方便計算，這里直接給出發動機點火建壓過程中的-關系式：

()=(1-e-)

(68)

其中，=6 MPa為內壓峰值；時間的單位為s；參數=20 s；這里加壓時間采用0.3 s。

圖11給出了三維發動機結構在20 ℃狀態下點火增壓結束時，發動機推進劑和絕熱層粘接界面上的法向和切向內聚力分布。與固化降溫階段不同的是，在點火增壓階段，藥柱的側面也是受到燃氣壓力作用的。因此，在增壓階段，切向內聚力并沒有在界面的兩端出現極值。在點火結束時，粘接界面上最大的切向和法向內聚力分別為0.11 MPa和-0.348 MPa。此外，云圖結果顯示，粘接界面上不存在界面單元失效情況，點火增壓結束時，粘接界面結構安全。

圖12給出了發動機粘接界面上靠近模型邊界處的三個特征點上切向和法向內聚力隨時間的分布。其中，點位于圓管段的端部，點位于翼槽段的中間部位，點位于星孔段的端部。可以看到，隨著點火階段的推進，內壓逐漸增大，粘接界面上的切向和法向內聚力都將逐漸變大。此外，由于空間位置的不同，切向內聚力方面差異較大，但是法向內聚力的差異較小，這與該階段處于受壓狀態有關。

(a)Tangential direction (b)Normal direction圖11 點火增壓結束粘接界面內聚力空間分布云圖Fig.11 Spatial distribution of traction on the bonding interface at the end of ignition pressurization stage

(a)Tangential direction (b)Normal direction圖12 點火增壓階段粘接界面內聚力隨時間變化曲線Fig.12 Traction on bonding interface versus time during ignition pressurization stage

4 結論

(1)通過圓管粘接結構中的裂紋擴展算例，驗證料三維八節點內聚單元，建立了含粘接界面復合結構的有限元仿真計算方法。相比于傳統的結構分析手段，該方法可以通過內聚單元直接獲得推進劑與絕熱層粘接界面上內聚力的分布，為粘接界面的結構完整性分析提供判斷依據。

(2)針對某發動機固化降溫和點火增壓工況，開展了三維發動機粘接界面的有限元仿真分析，得到如下結論：

1)在固化階段，由于被粘接材料熱脹系數的差異，推進劑收縮的速率明顯比絕熱層高，粘接界面受到持續的分離效應。固化降溫結束時，粘接界面上最大的切向和法向內聚力分別為0.03 MPa和0.065 MPa，遠小于粘接界面的內聚強度。

2)在點火增壓階段，粘接界面的法向持續受壓，界面的法向比較安全。但是由于持續的擠壓，界面兩側的切向位移不連續加劇，導致界面的切向內聚力增加明顯，相比于固化降溫結束階段，增加了近2.7倍。點火結束時，粘接界面上最大的切向和法向內聚力分別為0.11 MPa和-0.348 MPa，不存在界面單元失效情況，粘接界面結構安全。