固體推進劑藥柱局部夾雜應力場分析的邊界元法

黃拳章，陳 健，強洪夫，王學仁，岳春國，李愛華

(1. 火箭軍工程大學，西安 710025；2. 陜西師范大學，西安 710062)

0 引言

固體推進劑藥柱在澆注、固化降溫、長期貯存過程中，由于工藝、殘余應力和環境等因素的影響，不可避免地會出現夾雜、孔隙、脫粘、裂紋等缺陷形式[1-3]，這些缺陷是困擾固體火箭發動機質量、壽命和使用性能提高的重要因素。由于脫粘、夾雜和推進劑界面材料不連續，在發動機點火發射時高溫、高壓和高過載的作用下，夾雜和孔隙附近將出現局部高應力、應變區域，過高的應力、應變將會導致推進劑的脫濕、裂紋萌生和界面脫粘，而裂紋失穩擴展將會導致固體火箭發動機躥火，甚至爆炸等災難性后果[4]。因此，在固體火箭發動機藥柱結構完整性分析中，必須給予十分的關注和重點評估。為了獲得精確的應力-應變場，傳統有限元法通常采用網格局部加密的手段提高精度，由于藥柱幾何形貌復雜，且缺陷尺寸太小，勢必會帶來網格劃分的困難和巨量的計算耗費。

邊界元法具有精度高、適合復雜邊界形狀和降維求解等優點。因此，本文引入邊界元法，在固體推進劑藥柱高精度局部應力-應變分析中作為有限元分析的重要補充。

1 夾雜問題的數學模型

首先，對含有夾雜、孔隙的固體推進劑藥柱受力問題進行抽象和簡化，將藥柱抽象成一個均勻的物體域，其內部的顆粒和孔隙統一抽象為夾雜，建立如圖1所示的數學模型。圖1中，線彈性固體介質內含有n個任意形狀的夾雜或孔隙。其中，固體區域和其相應的外邊界分別為Ω和Γ，各夾雜區域和相應的邊界分別為Ω1，Ω2，…，Ωi，Ωn和Γ1，Γ2，…，Γi，Γn。另外，Γt和Γu分別為給定的面力邊界和位移邊界。此時每個夾雜不僅受外載荷的影響，還與相鄰夾雜(孔隙)之間發生相互作用，其受力狀態非常復雜。基體域與夾雜域之間相互作用，其界面上的位移和面力都是未知的，因此該問題不能直接求解。但考慮到界面位移和面力是受夾雜自身本構所控制的，因此它們并非是相互獨立的，只要找到兩者之間的函數關系，并將其作為已知邊界條件的補充條件，就能使原問題得到解決。基于此，本文建立了孔隙和固體夾雜與基體界面位移與面力的關聯矩陣，將其代入到基體子域的邊界元離散代數方程組中，以求解各夾雜邊界的位移和面力。

圖1 夾雜問題簡化模型

2 邊界元方法

2.1 基本方程

對于如圖1所示平面夾雜問題，如果不計體力，則邊界積分方程可寫成下列形式：

(1)

對于彈性力學平面應變問題，基本解和相應面力可寫為[5]

(2)

(3)

為了對積分方程進行數值求解，將所有邊界離散成若干個邊界單元，在每個邊界單元上對位移和面力分別進行插值。進一步在方程(1)中將p點依次取為離散邊界上的各個結點，整理后，可得到下列形式的線性代數方程組：

A0x0=B0y0

(4)

(5)

式(4)中系數矩陣A0和B0可用分塊形式表示：

(6)

式(6)中，分塊矩陣下標的含義為：第1和第2下標分別表示源點p和場點q所在的邊界線的類型；下標“1”表示該段邊界給定了面力邊界條件；“2”表示該邊界段給定了位移邊界條件；“3”表示該邊界為基體與夾雜的交界面；第1上標和第2上標分別表示源點和場點分別在第幾個夾雜界的面上，其中“0”代表基體的外邊界。

由于各夾雜界面的位移和面力都是未知量，使得式(4)中未知量的數目大于方程數，因此該方程不能定解。但考慮到夾雜表面位移和面力之間并不是相互獨立的，在線彈性框架內，兩者服從下列線性關系：

Tjj=DjUj

(7)

式中Tjj、Uj和Dj分別為第j個夾雜與基體界面靠近夾雜一側的面力、位移以及它們之間的關聯矩陣。

夾雜與基體界面靠近基體一側的面力T0j和Tjj是一對作用力和反作用力，則有

T0j=-Tjj=-DjUj

(8)

將式(8)代入式(4)，就可求得基體子域邊界的未知位移和面力，進而可求出基體和夾雜域內任一點的位移、應變和應力。

需要指出的是，以上公式是針對線彈性問題的，如果是粘彈性問題，則可根據彈性-粘彈性對應原理解決[6]。但對應原理需將粘彈性材料轉換成一系列的彈性材料，進而求解不同材料參數的彈性解，因此邊界元離散方程的系數矩陣將需要不斷更新，從而導致計算量成倍增加。此外，彈性-粘彈性對應原理對于非同質材料瞬態溫度問題以及邊界條件是時間的復雜函數的情況并不適用。基于此，本文采用Prony級數法描述粘彈性材料的本構關系，并將粘彈性效應以初應力的形式作為方程的右端項，進而給出時間域內粘彈性邊界元積分方程以及應力求解的迭代公式，使得邊界元方程系數矩陣在求解過程中始終保持不變，并用數值迭代方法進行求解，可同時計算出域內和邊界不同位置的應力和應變值。

不考慮體力和熱應力，則粘彈性問題的粘彈性邊界積分方程為[7]

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

式中，G為剪切模量；K=E/3(1-2ν)為體積模量；ξ和ξ′分別為受溫度影響的減縮時間，當材料不受溫度影響時，ξ=t，ξ′=t′。若溫度是連續變化的，則有

(17)

式中aT為移位因子，其大小可由實驗確定。

假設將邊界離散為n個單元，每個單元m(m>1)個節點，則令p點位于各個不同的邊界節點，則對式(9)進行數值積分后可形成關于節點未知量的方程組，如式(18)所示。

[H]2n×2n(m-1){u}2n(m-1)=[G]2n×2n(m-1){p}2n(m-1)+

(18)

(19)

將方程(18)中的未知量移到方程左邊，已知量移到右邊，可得

(20)

2.2 夾雜界面位移與面力的解耦

式(7)建立了夾雜界面面力與位移的線性關系，因此實現夾雜界面面力與位移解耦的關鍵是建立雜邊界面力與位移的關聯矩陣，這取決于夾雜材料的本構關系。下面分別按固體夾雜和孔隙兩種情況，詳細介紹夾雜界面面力與位移的解耦過程。

圖2 粘彈性邊界元迭代算法流程圖

2.2.1 固體夾雜

對于固體夾雜子域，可列出邊界積分方程如下：

(21)

將固體夾雜子域邊界進行分元離散后，由式(21)可得下列代數方程組：

(22)

由式(22)可得

(23)

2.2.2 孔隙和剛性夾雜

2.3 夾雜邊界的幾何近似

固體推進劑藥柱是通過固體推進劑攪拌澆注而成的，混進的夾雜顆粒以及推進劑配料如AP顆粒等經過大量的滾動和摩擦，棱角基本磨平，最終成球形的居多。因此，可將夾雜和孔隙近似為球形，二維情況下，可將夾雜和孔隙處理成圓形。而對于孔隙狀孔隙，則可用長細比較大的橢圓形進行近似。而為了幾何建模和網格劃分的方便，橢圓邊界可用四段兩兩相切的圓弧近似描述[8]，如圖3所示，切點分別為C1、C2、C3、C4，可證明每兩段弧的切點是唯一的。橢圓的長、短半軸分別為a和b，4個頂點依次為A1、A2、A3、A4。圓弧1和4的圓心分別為O1和O4，曲率半徑分別為R1和R4，且有圓弧1、4關于y軸對稱以及圓弧2、3關于x軸對稱。

圖3 橢圓孔的幾何近似

3 數值算例

3.1 中心含單個圓形夾雜的平面問題

考慮如圖4所示對邊受均勻壓力P作用的含單個圓形夾雜的平面應變問題，設板的邊長為a，夾雜半徑為r?；w的彈性模量E=2 GPa，泊松比ν=0.3，為了限制整體剛體位移，將平面左、右兩邊中點的豎直方向位移和上邊界中點的水平位移加以約束。

圖4 中心含單個圓形夾雜的平面

假設夾雜為線彈性材料，彈性模量E1=1.2E，泊松比ν1=ν，a=1000 mm，r=0.5 mm。由于夾雜的尺寸遠小于方板的尺寸，因此本算例可近似為無限大問題，則正確的數值解將趨近于解析解。求解時，將方板的四條邊界各均勻劃分40個線性直線型單元，夾雜邊界均勻劃分48個線性圓弧單元[9]。圖5給出了平面應變狀態下，固體夾雜界面徑向位移和徑向應力隨弧長的分布曲線，弧長沿順時針走向為正，起點弧度為π/2。由圖5可看出，邊界元解和解析解吻合很好，說明本文的算法和程序都是正確的。

(a)徑向位移

(b)徑向應力

3.2 圓筒固體火箭發動機燃燒室平面應變問題

如圖6所示，圓筒型固體火箭發動機平面應變問題。

設發動機內腔半徑為r0，殼體內外徑分別為r1和r2，殼體厚度為h=r2-r1。藥柱材料為典型的粘彈性材料，假設殼體材料為彈性材料，因為絕熱層和襯層的很薄，且材料性質與藥柱相似。為簡單起見，這里忽略絕熱層和襯層的厚度。發動機內腔受到壓力q(t)作用，如圖7所示。

圖6 圓筒發動機燃燒室平面應變模型

圖7 瞬時壓力模型

(24)

假設固體推進劑和發動機殼體的泊松比分別為ν和νk；E(t)和Ek分別為固體推進劑和發動機殼體的彈性模量；m=r1/r0為藥柱外內徑比。根據粘彈性解-彈性解對應原理進行理論解析解的計算，可得藥柱內徑r0處的粘彈性解[10-11]。

為簡單起見，先不考慮體力和熱應力，則由式(18)可知，在區域Ω1上有

(25a)

同理，在區域Ω2上有

(25b)

假設殼體和藥柱粘接完好，則殼體和藥柱在公共邊界Γ1上滿足位移和面力的連續條件：

(26)

由式(26)可知，式(25a)和式(25b)可分別改寫成下列形式：

(27a)

(27b)

綜合式(27a)和式(27b)可得到下列形式：

(28)

式(28)即是圓筒型固體火箭發動機平面應變問題的代數方程組，按照2.1節的邊界元離散方法和記憶應力迭代解法，即可求出該問題的解。此外，從該方程的形式也可看出，該方程系數矩陣具有帶狀特征，也可克服邊界元系數矩陣是滿陣的缺點，從而可節約內存并減少計算時間。

采取線性圓弧單元對圖6所示的模型進行離散，如圖8所示，1/12對稱模型，邊界Γ0、Γ1和Γ2上共劃分4個圓弧單元，一周共48個單元，為計算記憶應力和內點的應力應變，將藥柱區域和殼體劃分成四邊形等參單元(采用有限元網格)，其中藥柱內劃分28個四邊形網格，殼體內劃分8個四邊形網格，完整域內一共432個四邊形網格，模型的兩個側邊施加對稱性邊界條件，區域積分采用3×3高斯積分。

計算模型載荷參數、幾何參數和材料參數如下：

載荷參數：q0=8.9 MPa，t0=0.1 s；

模型尺寸：r0=0.3 mm，r1=1 mm，h=0.005 mm；

材料參數：Ek=206.8 GPa，νk=0.3，推進劑ν=0.495，初始模量E(0)=12.209 MPa。

松弛模量的Prony級數形式如下：

(29)

其中，τi=4×10(i-3)，剪切模量和體積模量計算式如下：

(30)

圖8 網格劃分示意圖

圖9 周向和徑向應力隨時間的變化

圖10 周向和徑向應變隨時間的變化

對粘彈性材料進行數值計算時，一般假設體積模量為常數，并取K(0)作為體積模量，這樣計算式(16)的遺傳積分時就可以簡化，只需計算記憶應力隨著剪切模量變化即可。圖9和10分別給出了周向和徑向應力、應變的解析解和數值解隨加載時間的變化曲線，從中可看出兩者吻合良好。

4 結束語

本文給出了粘彈性時域問題的邊界元法和求解含孔隙和固體夾雜問題的邊界元求解方案，通過算例驗證了方法的正確性和有效性。該方法可推廣到一般的三維問題，對于固體推進劑藥柱中含孔隙和顆粒夾雜的情況，則可精確計算夾雜界面局部的應變場和應力場，可為固體推進劑藥柱的結構完整性分析提供精確的數據。但傳統邊界元法的系數矩陣是非對稱的滿陣，求解代數方程組一般采用高斯消去法或迭代法，使得計算效率受到很大限制，進而影響了計算規模，下一步可通過引入快速多極邊界元法，以提高計算效率。