固體推進(jìn)劑藥柱局部夾雜應(yīng)力場分析的邊界元法

黃拳章，陳 健，強(qiáng)洪夫，王學(xué)仁，岳春國，李愛華

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

0 引言

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

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

1 夾雜問題的數(shù)學(xué)模型

首先，對含有夾雜、孔隙的固體推進(jìn)劑藥柱受力問題進(jìn)行抽象和簡化，將藥柱抽象成一個均勻的物體域，其內(nèi)部的顆粒和孔隙統(tǒng)一抽象為夾雜，建立如圖1所示的數(shù)學(xué)模型。圖1中，線彈性固體介質(zhì)內(nèi)含有n個任意形狀的夾雜或孔隙。其中，固體區(qū)域和其相應(yīng)的外邊界分別為Ω和Γ，各夾雜區(qū)域和相應(yīng)的邊界分別為Ω1，Ω2，…，Ωi，Ωn和Γ1，Γ2，…，Γi，Γn。另外，Γt和Γu分別為給定的面力邊界和位移邊界。此時每個夾雜不僅受外載荷的影響，還與相鄰?qiáng)A雜(孔隙)之間發(fā)生相互作用，其受力狀態(tài)非常復(fù)雜?；w域與夾雜域之間相互作用，其界面上的位移和面力都是未知的，因此該問題不能直接求解。但考慮到界面位移和面力是受夾雜自身本構(gòu)所控制的，因此它們并非是相互獨立的，只要找到兩者之間的函數(shù)關(guān)系，并將其作為已知邊界條件的補(bǔ)充條件，就能使原問題得到解決?；诖耍疚慕⒘丝紫逗凸腆w夾雜與基體界面位移與面力的關(guān)聯(lián)矩陣，將其代入到基體子域的邊界元離散代數(shù)方程組中，以求解各夾雜邊界的位移和面力。

圖1 夾雜問題簡化模型

2 邊界元方法

2.1 基本方程

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

(1)

對于彈性力學(xué)平面應(yīng)變問題，基本解和相應(yīng)面力可寫為[5]

(2)

(3)

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

A0x0=B0y0

(4)

(5)

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

(6)

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

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

Tjj=DjUj

(7)

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

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

T0j=-Tjj=-DjUj

(8)

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

需要指出的是，以上公式是針對線彈性問題的，如果是粘彈性問題，則可根據(jù)彈性-粘彈性對應(yīng)原理解決[6]。但對應(yīng)原理需將粘彈性材料轉(zhuǎn)換成一系列的彈性材料，進(jìn)而求解不同材料參數(shù)的彈性解，因此邊界元離散方程的系數(shù)矩陣將需要不斷更新，從而導(dǎo)致計算量成倍增加。此外，彈性-粘彈性對應(yīng)原理對于非同質(zhì)材料瞬態(tài)溫度問題以及邊界條件是時間的復(fù)雜函數(shù)的情況并不適用?；诖耍疚牟捎肞rony級數(shù)法描述粘彈性材料的本構(gòu)關(guān)系，并將粘彈性效應(yīng)以初應(yīng)力的形式作為方程的右端項，進(jìn)而給出時間域內(nèi)粘彈性邊界元積分方程以及應(yīng)力求解的迭代公式，使得邊界元方程系數(shù)矩陣在求解過程中始終保持不變，并用數(shù)值迭代方法進(jìn)行求解，可同時計算出域內(nèi)和邊界不同位置的應(yīng)力和應(yīng)變值。

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

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

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

(17)

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

假設(shè)將邊界離散為n個單元，每個單元m(m>1)個節(jié)點，則令p點位于各個不同的邊界節(jié)點，則對式(9)進(jìn)行數(shù)值積分后可形成關(guān)于節(jié)點未知量的方程組，如式(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)建立了夾雜界面面力與位移的線性關(guān)系，因此實現(xiàn)夾雜界面面力與位移解耦的關(guān)鍵是建立雜邊界面力與位移的關(guān)聯(lián)矩陣，這取決于夾雜材料的本構(gòu)關(guān)系。下面分別按固體夾雜和孔隙兩種情況，詳細(xì)介紹夾雜界面面力與位移的解耦過程。

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

2.2.1 固體夾雜

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

(21)

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

(22)

由式(22)可得

(23)

2.2.2 孔隙和剛性夾雜

2.3 夾雜邊界的幾何近似

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

圖3 橢圓孔的幾何近似

3 數(shù)值算例

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

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

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

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

(a)徑向位移

(b)徑向應(yīng)力

3.2 圓筒固體火箭發(fā)動機(jī)燃燒室平面應(yīng)變問題

如圖6所示，圓筒型固體火箭發(fā)動機(jī)平面應(yīng)變問題。

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

圖6 圓筒發(fā)動機(jī)燃燒室平面應(yīng)變模型

圖7 瞬時壓力模型

(24)

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

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

(25a)

同理，在區(qū)域Ω2上有

(25b)

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

(26)

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

(27a)

(27b)

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

(28)

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

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

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

載荷參數(shù)：q0=8.9 MPa，t0=0.1 s；

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

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

松弛模量的Prony級數(shù)形式如下：

(29)

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

(30)

圖8 網(wǎng)格劃分示意圖

圖9 周向和徑向應(yīng)力隨時間的變化

圖10 周向和徑向應(yīng)變隨時間的變化

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

4 結(jié)束語

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