郭 攀，武文華，劉 君，吳志剛

（1.大連理工大學(xué)運(yùn)載工程與力學(xué)學(xué)部工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連116024；2.大連理工大學(xué)運(yùn)載工程與力學(xué)學(xué)部航空航天學(xué)院，遼寧 大連116024）

在固體火箭發(fā)動(dòng)機(jī)（solid rocket motor，SRM）研制、生產(chǎn)、儲(chǔ)存過(guò)程中，會(huì)產(chǎn)生如裝藥界面脫粘、澆注氣泡、拔模表面損傷、運(yùn)輸震動(dòng)裂紋和存儲(chǔ)老化裂紋等缺陷，這些小缺陷對(duì)發(fā)動(dòng)機(jī)性能的影響復(fù)雜：有時(shí)沒(méi)有明顯變化，有時(shí)引起大小不同的推力異常，導(dǎo)致裝藥結(jié)構(gòu)破壞，甚至引發(fā)爆炸事故。因此，這些缺陷對(duì)結(jié)構(gòu)的安全有重要影響。對(duì)激波與裂紋的流固耦合問(wèn)題已有大量的研究［1－7］：何國(guó)強(qiáng)［1］針對(duì)缺陷對(duì)固體發(fā)動(dòng)機(jī)性能影響的研究成果進(jìn)行較細(xì)致的綜述，C.T.Liu［2－3］通過(guò)實(shí)驗(yàn)確定了裂紋增長(zhǎng)與應(yīng)變速率之間的關(guān)系，職世君等［4］采用數(shù)值模擬研究了對(duì)流燃燒對(duì)裂紋擴(kuò)展的影響，韓波等［5］對(duì)激波掃過(guò)缺陷時(shí)裂紋壁面的變形情況進(jìn)行了數(shù)值模擬和實(shí)驗(yàn)研究，徐春光等［6］和劉君等［7］采用數(shù)值模擬研究了激波在不同尺度缺陷中的繞射規(guī)律。

強(qiáng)激波掃過(guò)缺陷時(shí)形成的流場(chǎng)極其復(fù)雜，缺陷口存在激波的反射、繞射現(xiàn)象，缺陷內(nèi)存在馬赫反射并重新生成平面激波的現(xiàn)象，在缺陷內(nèi)還存在激波誘導(dǎo)渦運(yùn)動(dòng)及激波與漩渦的相互干擾等現(xiàn)象［6－7］。研究這種激波與漩渦占主導(dǎo)地位的高度非線性非定常流場(chǎng)，探索其物理機(jī)理，進(jìn)一步了解激波在缺陷中的傳播規(guī)律、激波沖擊作用下含缺陷裝藥的動(dòng)態(tài)響應(yīng)以及應(yīng)力強(qiáng)度因子的變化規(guī)律，有重要的學(xué)術(shù)價(jià)值和應(yīng)用價(jià)值。本文中，對(duì)激波作用下固體火箭裝藥缺陷的流固耦合相互作用進(jìn)行數(shù)值模擬。

流固耦合問(wèn)題的研究方法主要有緊耦合和弱耦合。前者將流體域、固體域和耦合作用構(gòu)造在同一控制方程中，在單一時(shí)間步內(nèi)同時(shí)求解所有變量；后者在每一時(shí)間步分別采用計(jì)算流體動(dòng)力學(xué)（computational fluid dynamics，CFD）和計(jì)算固體動(dòng)力學(xué)（computational structural dynamics，CSD）對(duì)流體域和固體域依次求解，通過(guò)搭建中間數(shù)據(jù)交換平臺(tái)交換固體域和流體域的結(jié)果信息，實(shí)現(xiàn)耦合求解。具體的弱耦合包括結(jié)構(gòu)和界面共同求解的弱耦合、流體和界面共同求解的弱耦合和流體和固體都考慮界面的弱耦合。弱耦合方式立足于各自求解方法，充分發(fā)揮各自求解的優(yōu)勢(shì)進(jìn)行學(xué)科交叉，有獨(dú)特的應(yīng)用優(yōu)勢(shì)［8－9］。本文中，采用流體和固體都考慮界面的弱耦合方法分析激波作用下含缺陷固體火箭裝藥的流固耦合問(wèn)題。

1 流體計(jì)算數(shù)值方法

在任意拉格朗日－歐拉（arbitrary Lagrangian－Eulerian，ALE）有限體積描述下，二維可壓縮非定常Euler方程的積分形式可表示為：

式中：Ω（t）為運(yùn)動(dòng)的流場(chǎng)區(qū)域，?Ω（t）為其邊界；A為流場(chǎng)區(qū)域的面積，l為其邊界長(zhǎng)度；u為網(wǎng)格運(yùn)動(dòng)速度；n為控制體外邊界的外法線單位矢量；Q和F分別為守恒變量和包含相對(duì)運(yùn)動(dòng)速度的對(duì)流項(xiàng)：

式中：ρ為流體密度；u和v分別為x和y方向的流動(dòng)速度；e為流體質(zhì)量?jī)?nèi)能；p為流體壓力；U為流體在ALE坐標(biāo)系下的流體速度；un為網(wǎng)格運(yùn)動(dòng)的法向速度；nx和ny分別為n的2個(gè)分量。分別用ux和uy表示u的2個(gè)分量，則：

采用格心格式的有限體積方法對(duì)非結(jié)構(gòu)單元進(jìn)行空間離散，由格林公式求得單元內(nèi)的變量梯度，重建單元邊界兩側(cè)的變量，形成近似的一維Riemann問(wèn)題，相對(duì)通量的分裂采用Van Leer分裂。為了抑制流場(chǎng)中物理量間斷處可能出現(xiàn)的數(shù)值振蕩，采用Venkatakrishman通量限制器。網(wǎng)格單元編號(hào)為i的半離散形式為：

式中：Qi為單元i的守恒變量在單元中心處的值，Vi為單元i的體積，F(xiàn)k、nk和Sk為單元i第k個(gè)表面邊界的通量、外法向單位矢量和面積，Nf為控制體的表面?zhèn)€數(shù)。為保證時(shí)間精度，流場(chǎng)方程采用四步Runge－Kutta方程進(jìn)行求解。

2 固體計(jì)算數(shù)值方法

在裂紋的數(shù)值模擬方面，有限元法應(yīng)用最廣泛，T.Belytschko等［10－11］最早將擴(kuò)展有限元（extended finite element method，XFEM）應(yīng)用于裂紋擴(kuò)展問(wèn)題。XFEM在裂紋擴(kuò)展的數(shù)值模擬中能避免網(wǎng)格劃分，節(jié)省求解時(shí)間，因此XFEM成為目前研究的熱點(diǎn)。與常規(guī)有限元法相比，XFEM不需要重新劃分網(wǎng)格，通過(guò)富集自由度描述裂紋及其生長(zhǎng)，大大提高了計(jì)算效率。

XFEM中常規(guī)有限元的位移近似解為：

式中：Ni為第i個(gè)節(jié)點(diǎn)的形函數(shù)；si為第i個(gè)節(jié)點(diǎn)的位移。對(duì)于彈性問(wèn)題含裂紋單元的近似位移插值函數(shù)可表示為：

式中：ai為裂紋穿過(guò)單元富集節(jié)點(diǎn)位移，b（α）i（α＝1～4）為裂尖所在單元富集節(jié)點(diǎn)位移；NΓ為被裂紋穿過(guò)但不包含裂尖的單元節(jié)點(diǎn)的集合；N＾為含裂尖的單元節(jié)點(diǎn)集合；r為裂尖區(qū)域節(jié)點(diǎn)到裂尖的距離；θ為裂尖區(qū)域節(jié)點(diǎn)與裂尖的連線與裂紋的夾角；H（x）為Heaviside函數(shù)，是反映裂紋面位移不連續(xù)的附加函數(shù)，用以反映位移不連續(xù)情況，表達(dá)式為：

動(dòng)力學(xué)支配方程為：

式中：M、C 和K 分別為結(jié)構(gòu)整體質(zhì)量陣、阻尼陣和剛度陣；d（t）、d（t）和¨d（t）分別為節(jié)點(diǎn)位移列矢量、節(jié)點(diǎn)速度矢量和節(jié)點(diǎn)加速度矢量，T為時(shí)域上限，f（t）為等效節(jié)點(diǎn)力矢量。對(duì)于常規(guī)節(jié)點(diǎn)i，位移di＝［si］；對(duì)于裂尖所在單元節(jié)點(diǎn)i，di＝［si，b（1）i，b（2）i，b（3）i，b（4）i］；對(duì)于Heaviside函數(shù)加強(qiáng)的節(jié)點(diǎn)i，di＝［si，ai］。節(jié)點(diǎn)i的節(jié)點(diǎn)力矢量為：

式中：b和t分別為體力和面力；Ωc和Γc分別為單元的積分域及其邊界；dΓ為邊界的積分微元。

對(duì)于節(jié)點(diǎn)i和節(jié)點(diǎn)j，單元?jiǎng)偠染仃嘖（e）ij可統(tǒng)一表示為：

式中：下標(biāo)x和y分別表示對(duì)x和y求偏導(dǎo)。

在動(dòng)態(tài)裂紋的數(shù)值模擬過(guò)程中，時(shí)域處理上通常利用差分法諸如Newmark、Wilson－θ等時(shí)域連續(xù)的方法對(duì)問(wèn)題進(jìn)行求解［10－11］。但在高頻強(qiáng)脈動(dòng)沖擊荷載作用下，傳統(tǒng)時(shí)程方法難以捕捉波陣面的間斷，求解過(guò)程往往帶來(lái)虛假的數(shù)值振蕩。針對(duì)傳統(tǒng)時(shí)程算法求解過(guò)程的上述問(wèn)題，X.Li等［12］和 W.Wu等［13］發(fā)展了時(shí)域間斷Galerkin方法（discontinuous Galerkin finite element method，DGFEM），即在空間域有限元方法分析基礎(chǔ)上，在時(shí)域?qū)ξ灰七M(jìn)行3次Hermite插值，對(duì)位移的導(dǎo)數(shù)進(jìn)行時(shí)間的線性插值離散。已有的研究結(jié)果表明［8，12－13］，DGFEM具備自動(dòng)引入數(shù)值耗散和濾去虛假的高階模式和數(shù)值振蕩效應(yīng)的能力，特別對(duì)在脈動(dòng)和沖擊荷載作用下結(jié)構(gòu)中波傳播過(guò)程的數(shù)值模擬，與Newmark算法相比，DGFEM具有良好的捕捉波的間斷并消除虛假數(shù)值振蕩的能力。

考慮動(dòng)力學(xué)方程空間域離散（見(jiàn)（8）式），DGFEM在時(shí)域離散時(shí)允許基本未知函數(shù)φ與其導(dǎo)數(shù)ψ在離散點(diǎn)處間斷，即在時(shí)刻tn，有：

式中：φn和ψn分別為離散的基本位置函數(shù)及其導(dǎo)數(shù)，φ＋n和φ－n分別為φn的右極限和左極限，ψ＋n和ψ－n分別為ψn的右極限和左極限，Δφn和Δψn為間斷點(diǎn)處的躍變。DGFEM具體推導(dǎo)分析過(guò)程見(jiàn)文獻(xiàn)［11－12］。本文中直接給出基本求解公式：

式中：Δt為時(shí)間步長(zhǎng)，φn＝φ＋n＝φ－n，Qn為離散守恒變量，滿足：

本文中，針對(duì)沖擊波作用的間斷、作用時(shí)間短的特點(diǎn)，采用由DGFEM與XFEM結(jié)合起來(lái)的DG＿XFEM對(duì)含裂紋缺陷的固體火箭裝藥問(wèn)題進(jìn)行數(shù)值模擬。

3 模型與算例分析

模擬算例網(wǎng)格如圖1所示。初始時(shí)刻，馬赫數(shù)為2.074 61的正激波入口位于裂紋前2cm處，出口位于6m處，如圖2所示。波前氣體靜止，氣體視為量熱完全氣體，比熱比γ為1.4。參考速度為1km／s。固體網(wǎng)格劃分見(jiàn)圖3，節(jié)點(diǎn)分布為20×20，實(shí)心點(diǎn)為壓力監(jiān)測(cè)點(diǎn)所在節(jié)點(diǎn)，星號(hào)為固支邊界，空圈為裂紋面所在單元的富集節(jié)點(diǎn)，四方塊為裂紋尖端所富集節(jié)點(diǎn)。靜態(tài)強(qiáng)度因子采用平面單邊裂紋受單向拉伸問(wèn)題中的強(qiáng)度因子計(jì)算公式（KI＝Fσ（πl(wèi)a）1／2）進(jìn)行計(jì)算，其中 F＝1.12－0.23（la／lb）＋10.6（la／lb）2－21.7（la／lb）3＋30.4（la／lb）4，la為裂紋長(zhǎng)度，lb為固體板的寬度，取σ＝106。

圖1 網(wǎng)絡(luò)模型劃分Fig.1 Model and mesh

圖2 初始時(shí)刻激波流線Fig.2 Shock wave at initial time

圖3 含裂紋的固體模型及網(wǎng)格劃分Fig.3 Mesh of solid with crack

圖4 監(jiān)控點(diǎn)壓力曲線Fig.4 Pressure of monitoring points

沿固體邊界及裂紋面選擇38個(gè)測(cè)壓點(diǎn)，其中測(cè)點(diǎn)P1、P20、P38的坐標(biāo)分別為（－5cm，5cm）、（－5cm，－5cm）和（－5cm，－0.26cm）。3個(gè)測(cè)壓點(diǎn)壓力曲線如圖4所示。從圖4可以看出，P20處于波后，壓力值不變；P1在激波經(jīng)過(guò)后，壓力有明顯的波動(dòng)；P38點(diǎn)處于裂紋尖端附近，激波到達(dá)裂紋尖端附近發(fā)生多次反射，經(jīng)過(guò)多次反射，壓力最終達(dá)到穩(wěn)定。不同時(shí)刻的流場(chǎng)密度等值線如圖5所示。由圖4～5可以看出，4 000步（2 869.1μs）后，裂紋中流場(chǎng)壓力趨于穩(wěn)定。在前5 000步荷載樣本中每隔100步取樣，對(duì)裂紋的動(dòng)態(tài)響應(yīng)進(jìn)行計(jì)算，P1處的位移d隨荷載步n變化如圖6所示。分別采用位移反推法和J積分法對(duì)裂紋尖端應(yīng)力強(qiáng)度因子K進(jìn)行計(jì)算，計(jì)算結(jié)果如圖7所示，其中J積分的積分半徑1為0.4cm，積分半徑2為1cm。

從圖6～7可以看出，應(yīng)力強(qiáng)度因子與裂紋張開(kāi)位移的變化趨勢(shì)基本一致。采用位移反推法與2種積分半徑的J積分計(jì)算結(jié)果基本吻合，驗(yàn)證了該方法的可靠性。

圖5 激波掃裂紋的密度等位線Fig.5 Density contour of shock wave

圖6 P1處y方向位移曲線Fig.6 Displacement curve along yaxis at P1

圖7 應(yīng)力強(qiáng)度因子計(jì)算比較Fig.7 Contrast of the stress intensity factors

4 結(jié) 論

基于弱耦合算法實(shí)現(xiàn)了流固耦合分析程序框架，并對(duì)沖擊波作用下固體火箭裝藥缺陷動(dòng)態(tài)響應(yīng)進(jìn)行數(shù)值模擬。從結(jié)果可以看出，激波掃過(guò)裂紋的過(guò)程中，在裂紋中經(jīng)過(guò)幾次反射，使得裂紋內(nèi)壓力表現(xiàn)出振蕩的現(xiàn)象，但峰值不斷升高最后趨于穩(wěn)定，體現(xiàn)了高度非定常非線性的特點(diǎn)；從應(yīng)力強(qiáng)度因子的計(jì)算看出，裝藥缺陷裂紋的應(yīng)力強(qiáng)度因子變化與沖擊波作用下位移響應(yīng)一致，驗(yàn)證了該方法的正確性。

［1］何國(guó)強(qiáng).裝藥缺陷對(duì)固體發(fā)動(dòng)機(jī)性能影響研究［D］.西安：西北工業(yè)大學(xué)，2001.

［2］Liu C T.Investigation of the effect of predamage on the crack growth behavior in a particulate composite material［R］.La Jolla，CA，USA：American Institute of Aeronautics and Astronautics，1993.

［3］Liu C T.Crack growth behavior in a composite propellant with strain gradients：PartⅡ［J］.Journal of Spacecraft，1990，27（6）：647－652.

［4］職世君，孫冰，張建偉.固體推進(jìn)劑復(fù)合型裂紋擴(kuò)展數(shù)值計(jì)算［J］.固體火箭技術(shù)，2011，34（1）：28－31.Zhi Shi－jun，Sun Bing，Zhang Jian－wei.Numerical computation of mixed mode crack propagation in solid propellant［J］.Journal of Solid Rocket Technology，2011，34（1）：28－31.

［5］韓波，周長(zhǎng)省，陳雄.固體推進(jìn)劑裂紋內(nèi)點(diǎn)火過(guò)程流固耦合數(shù)值仿真［J］.固體火箭技術(shù)，2011，34（2）：180－183.Han Bo，Zhou Chang－sheng，Chen Xiong.Numerical simulation of fluid－solid coupling in propellant crack during ignition［J］.Journal of Solid Rocket Technology，2011，34（2）：180－183.

［6］徐春光，劉君.激波在狹縫中繞射過(guò)程研究［J］.固體火箭技術(shù)，2003，26（3）：39－41.Xu Chun－guang，Liu Jun.Study of shock wave diffraction over a crack［J］.Journal of Solid Rocket Technology，2003，26（3）：39－41.

［7］劉君，張為華，郭正.點(diǎn)火發(fā)動(dòng)機(jī)形成的受限噴流流場(chǎng)數(shù)值模擬［J］.固體火箭技術(shù)，2000，23（1）：6－8.Liu Jun，Zhang Wei－h(huán)ua，Guo Zheng.Numerical simulation of limited jets produced by ignitors［J］.Journal of Solid Rocket Technology，2000，23（1）：6－8.

［8］郭攀，劉君，武文華.爆炸沖擊載荷作用下流固耦合數(shù)值模擬［J］.力學(xué)學(xué)報(bào)，2013，45（2）：283－287.Guo Pan，Liu Jun，Wu Wen－h(huán)ua.Numrical modeling for fluid－structure interaction under blast and impact loading response［J］.Chinese Journal of Theoretical and Applied Mechanics，2013，45（2）：283－287.

［9］Felippa C A，Park K，F(xiàn)arhat C.Partioned analysis of coupled mechanical systems［J］.Computer Methods in Applied Mechanics and Engineering，2001，190：3247－3270.

［10］Belytschko T，Chen H，Xu J，et al.Dynamic crack propagation based on loss of hyperbolicity with a new discontinuous enrichment［J］.International Journal for Numerical Methods in Engineering，2003，58（12）：1873－1905.

［11］Belytschko T，Chen H.Singular enrichment finite element method for elastodynamic crack propagation［J］.International Journal of Computational Methods，2004，1（1）：1－15.

［12］Li X，Yao D，Lewis R W.A discontinuous Galerkin finite element method for dynamic and wave propagation problems in non－linear solids and saturated porous media［J］.International Journal for Numerical Methods in Engineering，2003，57（12）：1775－1800.

［13］Wu W，Li X.Application of the time discontinuous Galerkin finite element method to heat wave simulation［J］.International Journal of Heat and Mass Transfer，2006，49：1679－1684.