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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對于節(jié)點i和節(jié)點j，單元剛度矩陣K（e）ij可統(tǒng)一表示為：

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

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

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

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

式中：Δt為時間步長，φn＝φ＋n＝φ－n，Qn為離散守恒變量，滿足：

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

3 模型與算例分析

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

［1］何國強.裝藥缺陷對固體發(fā)動機性能影響研究［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］職世君，孫冰，張建偉.固體推進劑復(fù)合型裂紋擴展數(shù)值計算［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］韓波，周長省，陳雄.固體推進劑裂紋內(nèi)點火過程流固耦合數(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］徐春光，劉君.激波在狹縫中繞射過程研究［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］劉君，張為華，郭正.點火發(fā)動機形成的受限噴流流場數(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é)報，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.