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

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

0 引言

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

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

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

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

圖1 夾雜問(wèn)題簡(jiǎn)化模型

2 邊界元方法

2.1 基本方程

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

(1)

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

(2)

(3)

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

A0x0=B0y0

(4)

(5)

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

(6)

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

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

Tjj=DjUj

(7)

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

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

T0j=-Tjj=-DjUj

(8)

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

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

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

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

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

(17)

式中aT為移位因子，其大小可由實(shí)驗(yàn)確定。

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

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

2.2.1 固體夾雜

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

(21)

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

(22)

由式(22)可得

(23)

2.2.2 孔隙和剛性?shī)A雜

2.3 夾雜邊界的幾何近似

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

圖3 橢圓孔的幾何近似

3 數(shù)值算例

3.1 中心含單個(gè)圓形夾雜的平面問(wèn)題

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

圖4 中心含單個(gè)圓形夾雜的平面

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

(a)徑向位移

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

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

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

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

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

圖7 瞬時(shí)壓力模型

(24)

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

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

(25a)

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

(25b)

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

(26)

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

(27a)

(27b)

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

(28)

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

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

計(jì)算模型載荷參數(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級(jí)數(shù)形式如下：

(29)

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

(30)

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

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

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

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

4 結(jié)束語(yǔ)

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