張 澤，高 峰，郜 陽，呂 銳

（空軍工程大學，西安，710051）

0 引 言

固體火箭發(fā)動機在生產(chǎn)制造、服役及點火發(fā)射過程中受到溫度變化、振動、重力等多種載荷的作用，不可避免地會對推進劑藥柱造成一些損傷。這些損傷缺陷的存在很可能在發(fā)動機工作過程中形成“超燃燒表面”，即燃面面積大于設(shè)計值，導致發(fā)動機內(nèi)彈道性能偏離設(shè)計指標。更加嚴重的情況下，當“超燃燒表面”產(chǎn)生的瞬時壓強峰值超過發(fā)動機殼體所能承受的最大壓強時，發(fā)動機會發(fā)生爆炸。

為了評判和預估含裂紋固體火箭發(fā)動機的工作性能，國內(nèi)外學者對裂紋的擴展及裂紋內(nèi)的流場分布情況做了相應的試驗及仿真。Kumar［1］通過試驗研究了點火增壓速率、裂紋尺寸及推進劑燃速等因素之間的相互影響關(guān)系。邢耀國［2-4］對燃燒條件下裂紋及脫黏面的擴展情況及危害性做了研究。何國強［5］對燃燒增壓過程脫黏擴展條件進行了試驗分析。Kuo等［6］的研究表明，點火升壓速率只有達到或超過1 GPa/s時，點火燃氣才可能對裂紋產(chǎn)生影響，使得裂紋產(chǎn)生失穩(wěn)擴展。桂曉波［7］在冷流沖擊條件下對藥柱進行流固耦合數(shù)值分析，計算結(jié)果表明在壓力上升過程中靠近頭部位置的藥柱內(nèi)孔壁的應力和應變更大。孫博［8］對發(fā)動機點火時藥柱不同深度的裂紋進行J積分計算，總結(jié)得出裂紋的擴展受裂紋深度的影響，并存在一個臨界值，超過這個臨界值裂紋會發(fā)生擴展。

文獻［9］～［12］中采用了不同的手段對推進劑藥柱進行了流固耦合計算，本文在上述學者研究的基礎(chǔ)上，針對點火沖擊過程中推進劑藥柱未被點燃的情況，進行含裂紋缺陷發(fā)動機流固耦合仿真研究，并對裂紋內(nèi)流場壓強和速度的變化對裂紋變形的影響以及裂紋變形對流場的反作用進行了分析。

1 物理模型

圖1為BATES發(fā)動機模型，該發(fā)動機進行過上千次點火試驗，擁有較高的精度，目前已成為美國空軍火箭推進實驗室（Air Force Rocket Propulsion Laboratory，AFRPL）衡量固體推進劑性能的標準發(fā)動機［13］。根據(jù)文獻［7］和文獻［8］所得的結(jié)論建立本文的裂紋計算模型，其中裂紋寬為2 mm，深度為30 mm，距離燃燒室頭部120 mm。

圖1 含裂紋發(fā)動機計算模型Fig.1 Calculation model of cracked motor

2 數(shù)學模型

2.1 基本假設(shè)

a）將燃氣視為組分凍結(jié)的理想氣體，遵循理想氣體狀態(tài)方程；

b）點火燃氣從燃燒室頭部進入燃燒室，升壓速率為1 Gpa/s；

c）由于氣體在流動過程中溫度變化幅度較小，不考慮溫度變化對氣體各項參數(shù)的影響；

d）僅考慮點火過程生成燃氣對藥柱及裂紋的沖擊作用，不考慮推進劑的燃燒及化學反應作用；

e）將推進劑視為各向同性的線性粘彈性材料，泊松比為常數(shù)；

f）不考慮重力對計算過程的影響。

2.2 氣相控制方程

二維不可壓黏性流動N-S方程可表達為

質(zhì)量守恒方程：

動量守恒方程：

能量守恒方程：

其中，

式中ρ為氣體密度；ui為氣體速度分量；p為氣體壓強；xi為空間中點坐標分量；E為能量；qi為熱流量；τij為切應力分量。

固體火箭發(fā)動機燃燒室內(nèi)燃氣流動內(nèi)通道流，流動過程屬于高雷諾數(shù)流動，因此湍流模型采用可實現(xiàn)型k-ε兩方程模式，其中湍流能量輸運方程和能量耗散輸運方程如下：

式中μt為渦黏性；Sij為平均速度應變率張量；k為湍動能；ε為湍流耗散率；τtij為雷諾應力。

在求解流場時，采用基于密度的有限體積法離散N-S 方程，并使用SⅠMPLE 算法進行求解，對求解區(qū)域采用分塊結(jié)構(gòu)化網(wǎng)格進行劃分，并對邊界層進行網(wǎng)格加密處理?？臻g離散采用二階迎風格式，時間離散采用四階龍格-庫塔方法。

2.3 固相控制方程

依據(jù)有限元原理，結(jié)構(gòu)動力學方程為

式中 [M]為總體結(jié)構(gòu)質(zhì)量矩陣；[C]為總體結(jié)構(gòu)阻尼矩陣；[K]為總體結(jié)構(gòu)剛度矩陣；{d?}為節(jié)點加速度列矢量；{d?}為節(jié)點速度列矢量；z5pd5vd為節(jié)點位移列矢量；[F(t)]為節(jié)點等效載荷列矢量。

對于固體區(qū)域，采用任意拉格朗日-歐拉（Arbitrary Lagrangian-Eulerian，ALE）有限元法對結(jié)構(gòu)區(qū)域進行離散。對計算區(qū)域采用四邊形結(jié)構(gòu)網(wǎng)格劃分，為了提高裂紋附近計算精度，對裂紋附近網(wǎng)格進行加密處理。在計算中，將前端面和外壁面約束固定，后端面和內(nèi)壁面為流固耦合面，采用Newmark時域連續(xù)方法對動態(tài)問題進行求解。圖2為裂紋區(qū)域網(wǎng)格示意。

選取HTPB推進劑為研究對象，將推進劑作為黏彈性材料處理，材料密度為1 600 kg/m3，泊松比為0.495，表1為松弛彈性常數(shù)。

表1 Prony級數(shù)參數(shù)Tab.1 Prony parameters

粘彈性材料剪切松弛核函數(shù)G(t)和體積松弛核函數(shù)K(t)采用Prony級數(shù)表示為

式中G∞，Gi為剪切模量；K∞，Ki為體積模量；為各Prony級數(shù)分量的松弛時間。

G(t)和K(t)可用松弛模量E(t)表示為

將松弛模量E(t)表示為Prony級數(shù)形式［14］：

2.4 流固耦合計算方法

流固耦合計算方法通常分為強耦合和弱耦合。強耦合方法在同一控制方程內(nèi)同時對流體域與固體域進行求解。弱耦合方法分時間步對流體域和固體域分別進行求解，通過兩個區(qū)域之間的網(wǎng)格插值實現(xiàn)耦合計算。相比于強耦合方法，弱耦合更容易實現(xiàn)計算收斂，并且節(jié)省計算資源，因此本文選取弱耦合的方法進行計算。

由于采用弱耦合方法進行計算，流體域和固體域之間網(wǎng)格節(jié)點信息不匹配，在進行插值運算時需要尋找與兩區(qū)域耦合面上距離最近的節(jié)點進行匹配。節(jié)點之間建立匹配之后，將流體區(qū)域的壓強等信息傳遞到固體區(qū)域，然后求解固體區(qū)域的節(jié)點位移、速度并通過插值傳遞至流體區(qū)。接著采用動網(wǎng)格技術(shù)分別對流體域與固體域的網(wǎng)格進行更新，實現(xiàn)耦合面的移動，完成一個時間步的運算，按上述步驟進行迭代運算，直至收斂。

為了保證計算精度，需要對收斂進行判斷，應力收斂判斷表達式為

式中σk+1，σk分別為k+1和k迭代步耦合界面處的結(jié)構(gòu)應力；εu為設(shè)定的應力收斂值；ε0為一個非常小的常數(shù)。速度和應變的收斂判定與應力形式相同，達到收斂標準后即認為計算結(jié)束。

3 計算結(jié)果與分析

點火藥被點燃后，生成的燃氣從燃燒室頭部向噴管出口方向欠膨脹噴出，在燃燒室內(nèi)原有氣體的阻滯作用下，形成激波以及膨脹波，并在燃氣流過區(qū)域迅速形成高壓。隨著燃氣不斷向前運動，高壓區(qū)域不斷擴大，燃氣前峰到達噴管后，與噴管壁面形成反射波，沿與主通道燃氣相反方向運動，并與主通道燃氣相互作用，圖3為不同時刻燃燒室內(nèi)壓強分布云圖。

圖3 不同時刻燃燒室內(nèi)壓強分布云圖Fig.3 Pressure contours in the combustion chamber at different moments

燃氣在流動過程中，在激波的沖擊以及燃氣壓力作用下，推進劑出現(xiàn)變形并產(chǎn)生應力。在t=0.15 ms時刻，燃氣前峰到達裂紋處，與裂紋壁面碰撞形成壓縮波，如圖4所示。隨著燃氣的繼續(xù)運動，裂紋腔內(nèi)開始注入燃氣，裂紋入口處壓強開始快速升高，因此在該時刻，裂紋內(nèi)最大壓強出現(xiàn)在入口處，并且隨著裂紋深度的增加壓強逐漸降低。不同時刻裂紋內(nèi)部沿裂紋方向壓強分布如圖5所示。在t=0.25 ms時到達燃氣末端與壁面碰撞后被反射形成反射波，在裂紋末端形成高壓區(qū)域，該時刻裂紋腔內(nèi)最大壓強出現(xiàn)在裂紋末端處，沿裂紋向出口方向逐漸降低。在高壓的作用下，裂紋末端產(chǎn)生的形變大于裂紋內(nèi)其他區(qū)域，形成一個鼓形區(qū)域，裂紋變形呈現(xiàn)內(nèi)寬外窄的形狀。隨著壓強反射波繼續(xù)向裂紋入口處運動，在t=0.4 ms時刻裂紋中部壓強開始升高，此時壓強最大區(qū)域出現(xiàn)在裂紋中部，向兩邊依次遞減，并且由高壓產(chǎn)生的鼓形區(qū)域也隨著反射波向裂紋開口處運動，裂紋基本呈現(xiàn)內(nèi)外同寬的形狀。在t=0.6 ms時壓強反射波到達裂紋入口處時，形變最大區(qū)域也出現(xiàn)在裂紋入口處，此時裂紋呈現(xiàn)內(nèi)窄外寬的形狀，由于裂紋出口的擴張，裂紋腔內(nèi)的反射波釋放至燃燒室主流區(qū)域，因此出口處壓強略低于裂紋內(nèi)部壓強，裂紋腔內(nèi)壓強分布逐漸趨于均勻。不同時刻裂紋位移、應力及應變云圖分別如圖6～7所示。

圖4 燃氣前峰與裂紋作用壓強云圖Fig.4 Pressure contour of gas front and crack action

圖5 不同時刻裂紋內(nèi)壓強分布曲線Fig.5 Curves of pressure distribution in cracks at different moments

圖6 不同時刻裂紋位移云圖Fig.6 Crack displacement contours at different moments

圖7 不同時刻裂紋應力云圖Fig.7 Crack stress contours at different moments

圖8 為裂紋腔內(nèi)不同時刻燃氣速度分布。在t=0.15 ms時刻，燃氣開始注入裂紋，裂紋入口處速度較大，隨著深度增加燃氣速度逐漸趨于0。此后裂紋內(nèi)燃氣速度不斷提高，直至t=0.25 ms 時壓強波到達裂紋末端后形成反射波，由反射波產(chǎn)生負的壓強梯度，使裂紋內(nèi)燃氣速度開始減小。當t=0.6 ms時，整個裂紋內(nèi)部燃氣速度幾乎為0，這是因為裂紋內(nèi)壓強升高，主通道內(nèi)燃氣進入裂紋受阻。

圖8 不同時刻裂紋內(nèi)燃氣速度分布曲線Fig.8 Distribution curve of gas velocity in crack at different moments

上述仿真計算所得流場壓強及速度變化規(guī)律，與Kumar［15］在試驗及仿真中所得結(jié)論基本相符，驗證了本文所采用的計算方法的合理性。在本文的研究基礎(chǔ)上可對工程中所遇到的含裂紋缺陷的固體火箭發(fā)動機進行仿真計算，從而對其工作性能進行初步判斷，以便為其延壽或判廢工作提供參考。

4 結(jié) 論

本文研究得出以下結(jié)論：

a）點火生成燃氣進入燃燒室后，推進劑藥柱受到激波以及高壓的沖擊作用發(fā)生變形，變形導致流場邊界發(fā)生變化，從而影響燃氣氣流的流動，因此在發(fā)動機流場計算中不能忽視固體區(qū)域的變形。

b）燃氣進入裂紋后膨脹加速，在裂紋前端達到最大速度，在向末端運動的過程中不斷減速，到達裂紋末端后發(fā)生反射形成反射波，造成裂紋內(nèi)部壓強進一步升高，導致裂紋變形加劇，反射波在運動過程中使裂紋最大變形區(qū)域不斷向入口處擴展，最終形成內(nèi)窄外寬的擴張型裂紋。

c）通過流固耦合方法計算了裂紋在點火燃氣作用下產(chǎn)生的形變，并且分析了裂紋附近產(chǎn)生的應力?？稍诖嘶A(chǔ)上進一步采用斷裂力學的方法對裂紋的擴展問題進行判斷，從而為含裂紋缺陷的固體火箭發(fā)動機的工作性能的評估工作提供數(shù)據(jù)參考。