沖擊作用下推進劑變形的流固耦合分析方法

金銘君，李 強

(西北工業(yè)大學 燃燒、熱結(jié)構(gòu)與內(nèi)流場重點實驗室，西安 710072)

沖擊作用下推進劑變形的流固耦合分析方法

金銘君，李 強

(西北工業(yè)大學 燃燒、熱結(jié)構(gòu)與內(nèi)流場重點實驗室，西安 710072)

固體火箭發(fā)動機的點火過程是一個復(fù)雜多變的理化過程，具有時間短、升溫、升壓梯度大等特點。針對固體火箭發(fā)動機點火過程中的裝藥結(jié)構(gòu)完整性問題，文中建立了一套用于分析沖擊作用下固體推進劑變形現(xiàn)象的仿真模型。采用RANS和ALE方法，分別對流體域和固體域進行求解，以兩場獨立交叉耦合迭代的模式實現(xiàn)了仿真過程。以一個推進劑冷流沖擊實驗作為算例，對仿真模型進行了驗證，計算值與測量值間誤差不超過10%，仿真模型計算可靠，具有向固體火箭發(fā)動機實際點火過程拓展的價值。

固體火箭發(fā)動機；點火過程；推進劑變形；流固耦合

0 引言

固體火箭發(fā)動機的點火過程是一個短暫的快速升壓過程。其中，高溫高壓燃氣對固體推進劑裝藥的沖擊作用，會導致裝藥結(jié)構(gòu)的變形、破損和移位，甚至導致發(fā)動機失效爆炸。因此，為確保整個發(fā)動機的工作性能及可靠性、安全性，對發(fā)動機點火瞬態(tài)過程開展裝藥結(jié)構(gòu)完整性及瞬態(tài)流暢變化研究是非常有必要的。一直以來，固體火箭發(fā)動機的研制和定型主要以試驗為主，周期長、耗資大，加之燃燒室內(nèi)工作環(huán)境惡劣，點火過程又極為短暫，且復(fù)雜多變，對其進行實驗研究是極為困難的。近年來，隨著各種高精度、高效率計算方法的提出和完善，以及高速計算機的快速發(fā)展，計算機數(shù)值仿真技術(shù)憑借其高效、經(jīng)濟、安全等優(yōu)勢，在火箭發(fā)動機的研制中，占據(jù)了愈發(fā)重要的地位。

國內(nèi)外大量研究表明，升壓梯度是破壞裝藥結(jié)構(gòu)完整性的重要因素[1]。在點火階段，過高的壓強峰或升壓速率會對推進劑藥柱產(chǎn)生沖擊，造成藥柱變形或移位，而藥柱變形或移位又使得流動區(qū)域發(fā)生了變化。對于這樣一個典型的流固交互問題，采用流固耦合法對其分析，無疑是最佳選擇。1982年，Tante G L將流體和結(jié)構(gòu)之間的耦合關(guān)系作為動力學模型，驗證了燃燒室和藥柱的相對運動，建立了用于分析SRM點火爆炸原因的二維非定常模型。1991年，Titan IV SRMU地面試車失敗，掀起了針對點火沖擊影響的流固耦合方法研究熱潮。1991年，Johnston采用二維流場耦合結(jié)構(gòu)變形模型，計算了內(nèi)壓作用下藥柱變形及其對流場的影響，給出了較為安全的藥柱臨界模數(shù)[2]。此后，許多研究小組開展了此方面的研究工作，并得到了有價值的結(jié)論。國內(nèi)報道的主要有海軍航空工程學院飛行器工程系通過MPCCI連接FLUENT和ABAQUS軟件，對SRM快速升壓過程的研究[1]；以及南京理工大學機械工程學院利用CFX和ANSYS軟件，對固體推進劑裂紋內(nèi)點火過程的流固耦合研究等[3]。

本文著眼于具有良好適應(yīng)性的高精度SRM點火過程流固耦合方法研究，提出了一種適用的流固耦合計算分析模型。由于實際發(fā)動機點火過程的推進劑應(yīng)變和變形數(shù)據(jù)無法通過實驗獲得，因此本文以一個推進劑冷流沖擊實驗作為算例，對計算模型進行了校驗，通過對比計算結(jié)果與實驗結(jié)果中測試段內(nèi)的壓強曲線和推進劑變形最大位移的大小，來對仿真模型的準確性做出評價。

1 計算模型

1.1 流動計算模型

針對氣流對推進劑沖擊過程的流動特性，對流動模型做如下假設(shè)：

(1)計算中所涉及流體為有粘可壓縮牛頓流體，且可視為各向同性的連續(xù)介質(zhì)；

(2)冷流沖擊過程中流體溫度脈動較小，忽略由溫度脈動帶來的粘性系數(shù)的脈動；

(3)流場運算中忽略重力影響。

將密度ρ、速度ui、壓力p及溫度T經(jīng)上式處理后，代入有粘可壓縮N-S方程組，即可得到可壓縮雷諾控制方程組，進行必要簡化后，有[8]

(1)

(2)

(3)

1.2 湍流模型

本文采用渦粘性方法中的一方程模型：Spalart-Allmaras模型(S-A模型)來求解雷諾平均方程。

對于一方程模型的一般形式，有

(4)

式中SP為生成項；SD為耗散項；D為擴散項。

(5)

其中

這里，d為到物面的最近距離。對于式中其他參變量的關(guān)系，有

g=r+Cw2(r6-r)

上述各式中，常系數(shù)取值如表1所示。

表1 S-A模型中常系數(shù)取值

1.3 推進劑結(jié)構(gòu)力學計算模型

對推進劑結(jié)構(gòu)模型做如下假設(shè)：

(1)認為沖擊時間短暫，固體推進劑的粘彈特性不能完全體現(xiàn)，將其近似為大變形彈性模型；

(2)將推進劑材料視為均質(zhì)的各向同性Kirchhoff材料，忽略重力影響。

對于上述假設(shè)，材料本構(gòu)關(guān)系[9-10]滿足：

Sij=Cijk1Ek1=λEkkδij+2μEij

(6)

其中，Sij為基爾霍夫應(yīng)力；Eij為格林應(yīng)變；Cijk1為彈性模量的四階張量；λ和μ為拉梅(Lamé)常數(shù)，可表示為

λ=υE/(1+υ)(1-2υ)

μ=E/2(1+υ)

1.4 非線性有限元方法

采用ALE(Arbitrary Lagrangian-Eulerian)有限元方法，對結(jié)構(gòu)域進行離散。在ALE方法中，由于質(zhì)量守恒被強加作為偏微分方程，所以需要建立一種弱形式[10]。將ALE域劃分為單元，定義所有非獨立變量為單元坐標的函數(shù)，對于單元e，其ALE坐標為

χ(ξe)=Φe(ξe)=χ1N1(ξe)

其中，ξe為單元e的坐標。對于網(wǎng)格運動，有

其中，x1(t)為節(jié)點的運動。對于網(wǎng)格速度，有

對于有限元矩陣方程，有連續(xù)方程：

其中,Mρ、Lρ、Kρ分別為容量、轉(zhuǎn)換和散度矩陣：

對于動量方程，有：

其中M和L分別是廣義質(zhì)量和傳遞矩陣，對應(yīng)于在參考構(gòu)型描述下的速度fint和fext分別為內(nèi)力和外力向量：

1.5 耦合計算

目前，流固耦合問題的求解方式有兩類：兩場交叉迭代和直接全場同時求解。后者計算量巨大，且具體實現(xiàn)十分復(fù)雜，故本文采用兩場分開進行交叉迭代計算，各自單獨使用獨立解算器和網(wǎng)格進行解算的方式。在兩場數(shù)據(jù)交換的處理上，一般是依靠在兩套非匹配網(wǎng)格間進行雙向插值來完成的，其精度受兩插值網(wǎng)格間重合度等諸多因素影響。本文采用在交界面處通過插值重構(gòu)一層虛網(wǎng)格進行數(shù)據(jù)傳遞和相關(guān)變量解算的耦合方法，數(shù)據(jù)通過最小二乘法在兩場和虛擬網(wǎng)格間進行傳遞，減少了來自雙向插值帶來的誤差；交界面節(jié)點力可在虛網(wǎng)格上直接進行計算，節(jié)省了計算時間和計算量。設(shè)上角標nf為流體解算器中的時間步，ns則為結(jié)構(gòu)解算器中的時間步；下標f代表流體解算器中的變量，s代表結(jié)構(gòu)解算器中的變量，f-s則代表交界面處的變量。

通常固體膨脹波速cd比流體中波速c大，出于對算法穩(wěn)定性的考慮，固體解算器中的時間步長ts通常小于流體解算器中的時間步長tf。因此，兩場求解器中分別取其最佳時間步長進行交叉迭代，其具體流程和數(shù)據(jù)傳遞如圖1和圖2所示。

2 實驗

2.1 實驗方法

該實驗采用一套冷流來流模擬裝置來模擬發(fā)動機點火過程中壓強迅速升高的燃氣流，對縮比的推進劑進行沖擊實驗。實驗開始前，先在高速攝影儀驅(qū)動程序界面上對推進劑試件取一系列距離已知的標定點；實驗時，通過高速攝影儀來獲得實驗過程的影像資料，并在測試段首尾兩端配置壓力傳感器來對內(nèi)腔壓強進行監(jiān)測；實驗結(jié)束后，對影像資料中標定點間的距離變化情況進行分析，獲得推進劑形變、形變速度隨時間變化的曲線[11-13]。

2.2 實驗裝置

實驗裝置主要由3部分組成：氮氣恒壓供氣系統(tǒng)，推進劑變形測量段及數(shù)據(jù)采集系統(tǒng)。氮氣恒壓供氣系統(tǒng)由4個氮氣氣瓶和穩(wěn)壓裝置組成，可提供2～4.5 MPa的沖擊壓強。實驗測量段內(nèi)腔尺寸為240 mm×200 mm×140 mm，內(nèi)有固定槽，可保證在測量過程中試件沒有移動；前端為30 mm厚的有機玻璃，在保證測量段密閉性的同時，又兼顧了圖像采集系統(tǒng)對視野清晰度的要求；實驗采用兩路測壓，測壓孔配置在測量段后蓋板上，測量點P1臨近氣流入口處，測量點P2臨近氣流出口位置。數(shù)據(jù)采集系統(tǒng)分為圖像采集系統(tǒng)和壓強測量系統(tǒng)兩部分。圖像采集部分采用美國Phanton高速攝影儀(V4.3)，所采用的傳感器分辨率為800×600，采樣率為300幀；壓強測量部分采用Dewetron并行數(shù)據(jù)采集系統(tǒng)。

2.3 實驗方案

實驗對不同沖擊壓強下推進劑的變形情況進行了測量，為了減小實驗誤差，對每個工況進行了3次實驗，最終結(jié)果為3次實驗結(jié)果的平均值。本文采用前文所述計算模型，對其中沖擊壓強為3.3、3.6、4.1 MPa的工況開展了數(shù)值計算工作。

3 結(jié)果對比

3.1 實驗結(jié)果

圖4給出了在沖擊壓強為3.3、3.6、4.1 MPa的工況下，氣流入口和出口處壓強隨時間變化曲線。其中，曲線p1代表氣流入口處壓強，p2代表氣流出口處壓強。

在各個工況下，沖擊壓強達到最大值需要約0.12 s，與實際發(fā)動機點火過程的壓強建立時間大致相當，可認為該冷流推進劑沖擊實驗具備一定對真實發(fā)動機點火過程的參照性。

在各個工況下，氣流入口與出口處存在明顯的壓強差，分析認為這是造成推進劑變形的主要原因之一。推進劑形變的形狀與氮氣沖擊的位置有關(guān)，在該實驗中，沖擊位置位于推進劑的右側(cè)，推進劑受沖擊后向左凹陷，位移最大處為氣流沖擊位置的中央處，如圖5所示。各工況下入-出口壓強差、推進劑最大位移及形變速度如表2所示。

沖擊壓強/MPa前后壓差/MPa推進劑形變大小/mm推進劑形變速率/(mm/s)3．30．1610．70846．9363．60．2360．81753．1304．10．3291．19259．792

3.2 計算結(jié)果

參照上述實驗，依據(jù)其設(shè)計參數(shù)配置物理模型和邊界條件，開展了對推進劑冷流沖擊的數(shù)值計算工作。流動求解器采用了有限體積方法，空間對流項的離散采用了HLLC格式，粘性項差分采用WENO格式，時間項上采用四階Runge-Kutta積分方法。在計算中，參照實驗測量到的推進劑位移最大位置，設(shè)置了位移監(jiān)測點。由于實驗結(jié)果分析中只給出了氮氣沖擊造成的推進劑變形的最大位移結(jié)果，而未能給出具體位移情況隨時間變化的過程，考慮到氮氣的沖擊作用只是一個脈沖過程，因此在數(shù)值計算中，在氣流入口位置給定一個單脈沖過程的沖擊壓力作為初始邊界條件。計算域分為16個分區(qū)，流體計算網(wǎng)格總數(shù)約為28萬，其計算結(jié)果如下：

首先，由表3可知，各工況下最大位移的計算值和測量值間誤差不超過10%，可認為計算結(jié)果較為準確。其次，在圖6中給出的推進劑最大位移隨時間變化曲線可看到，在沖擊過后最初的0.02 s內(nèi)，3個工況下推進劑均沒有發(fā)生明顯的形變。此后，推進劑開始變形，且形變量逐漸變大。對比圖4可知，在沖擊最初的0.02 s內(nèi)，雖然氮氣氣源和實驗器內(nèi)的絕對壓差很大，但實驗器內(nèi)氣流入口和出口間的相對壓差卻較小，因此推進劑的變形很不明顯；而在0.02 s之后，隨氣流入口和出口間相對壓差的逐漸增大，推進劑變形逐漸明顯，且在壓差達到最大值的時刻，推進劑的位移變化幾乎同時達到最大值。最后，在沖擊過后的降壓階段，由于計算邊界無法準確模擬氮氣壓強衰減的過程，可認為此后的計算結(jié)果沒有意義。

沖擊壓強/MPa推進劑最大位移/mm計算值實驗值3．30．7210．7083．60．8920．8174．11．0861．192

由于實驗條件所限，實驗結(jié)果并沒有給出推進劑變形大小隨沖擊壓強漸變的具體過程，因此在計算結(jié)果中，只有推進劑形變的最大值和實驗結(jié)果具有比照作用，但其得到的過程值可認為是具有一定預(yù)驗價值的。通過對比計算結(jié)果與實驗結(jié)果中3個工況下的推進劑變形最大位移以及測試段內(nèi)升壓過程的壓強曲線，驗證了仿真模型的可行性和準確性。

4 結(jié)束語

(1)為了對氣流沖擊作用下推進劑的變形現(xiàn)象開展數(shù)值仿真研究，本文提出了一套基于RANS和ALE方法的流固耦合仿真模型。

(2)通過在交界面構(gòu)造虛擬網(wǎng)格的方法，實現(xiàn)了兩場間的數(shù)據(jù)交互，既減少了來自雙向插值帶來的誤差，又節(jié)省了計算時間和計算量。

(3)以推進劑冷流沖擊實驗作為算例，對仿真模型進行了校驗，計算結(jié)果與實驗結(jié)果誤差不超過10%，驗證了仿真模型的可靠性和準確性。

[1] 于勝春, 趙汝巖, 許濤, 等.固體火箭發(fā)動機快速升壓過程的流固耦合分析[J].固體火箭技術(shù), 2008, 31(3):232-235.

[2] Johnston W A. Flow-structural analysis of the Ariane-5 solid rocket motor during ignition transient[R]. AIAA 2011-4056.

[3] 韓波, 周長省, 陳雄.固體推進劑裂紋內(nèi)點火過程流固耦合數(shù)值仿真[J].固體火箭技術(shù), 2011, 34(2):180-183.

[4] 閻超.計算流體力學方法及應(yīng)用[M].北京: 北京航空航天大學出版社, 2006.

[5] Mark Salita. Modern SRM ignition transient modeling(part 1): introduction and physical models[R]. AIAA 2001-3443.

[6] 王新月, 胡春波, 張堃元, 等.氣體動力學基礎(chǔ)[M].西安: 西北工業(yè)大學出版社, 2006.

[7] David C Wilcox. Turbulence modeling for CFD[M]. California: DCW Industries, Inc., 1993.

[8] 陳矛章.粘性流體動力學基礎(chǔ)[M].北京: 高等教育出版社, 1993.

[9] Toro E F. Riemann solvers and numerical methods for fluid dynamics: a practical introduction[M]. Berlin:Springer-Verlag Berlin and Heidelberg GmbH & Co. K, 1999.

[10] Ted Belytschko, Wing Kam Liu, Brian Moran. Nonlinear finite elements for continua and structures[M]. London: John Wiley & Sons Ltd., 1998.

[11] Soares D Jr, Rodrigues G G, Gonc alves K A. An efficient multi-time-step implicit-explicit method to analyze solid-fluid coupled systems discretized by unconditionally stable time-domain finite element procedures[J]. Computers and Structures, 2010, 88(5-6):387-398.

[12] Luciano Garelli, Rodrigo R Paz, Mario A Storti. Fluid-structure interaction study of the start-up of a rocket engine nozzle[J]. Computers & Fluids, 2010, 39(7):1208-1218.

[13] Guerraiche Bilal. Invstigation on mechanical behavior of composite propellant during ignition process[D].西安: 西北工業(yè)大學, 2013.

(編輯：崔賢彬)

A solid-fluid interaction model for propellant deformation under impact condition

JIN Ming-jun，LI Qiang

(Science and Technology on Combustion, Internal Flow and Thermal-Structure Laboratory,Northwestern Polytechnical University, Xi'an 710072, China)

The ignition of SRM is a complicated and changeful physicochemical process, which has the features of extremely short duration and particularly high temperature & pressure gradient. An applicable solid-fluid interaction simulation model was established to analyze the phenomenon of propellant deformation under impact effect for the purpose of dealing with grain structure integrity problem during SRM ignition transient. RANS and ALE method were used independently in fluid and solid sub-domain. The interaction between two sub-domains was accomplished by the transmission of interface forces. The model was validated with a cold flow test. The results of simulation agree well with test data. It shows that the model has a relatively high accuracy and could be used to investigate the ignition process in a practical SRM.

solid rocket motor；ignition transient；propellant deformation；solid-fluid interaction

2016-02-25；

2016-10-27。

金銘君(1988—)，男，碩士生，從事固體火箭發(fā)動機點火過程流固耦合研究。E-mail:jackking413@163.com

V512

A

1006-2793(2017)02-0158-06

10.7673/j.issn.1006-2793.2017.02.005