激波绕射动态裂纹的数值模拟*

罗 雷，孙振生，张世英

(第二炮兵工程大学，西安 710025)

0 引言

固体火箭发动机在生产、储存以及使用过程中，装药表面或内部不可避免的形成裂纹。在点火冲击波的作用下，这些裂纹有可能变形、扩展，进而会导致药柱结构的破坏，甚至造成发动机爆燃等恶性事故。国内外学者在点火冲击波与裂纹的相互作用以及裂纹的扩展等方面做了大量的工作，认为点火阶段压力梯度的迅速升高是导致裂纹扩展的重要因素［1-4］。点火过程中，激波在裂纹内的绕射和反射是压力梯度变化的重要原因，但对这一非定常过程的研究还非常少。因此，研究激波在裂纹中的传播过程，揭示点火冲击波与装药裂纹相互作用的规律和机理，对发动机判废准则的制定具有重要意义。

激波绕射裂纹计算的一个主要难点是对边界的处理。文献［5］在假设裂纹壁面静止的情况下对强激波绕射裂纹所形成的流场进行了详细分析，得出了许多有价值的结论。但是随着推进剂的燃烧，计算区域是不断变化的。假设裂纹壁面静止对流场结构的影响是无法估计的，而引入移动边界进行计算可有效地解决这一问题。

LevelSet作为一种有效的移动界面追踪方法，已广泛应用于多个领域［6］。文中利用高精度加权本质无振荡(weighted essentially non-oscillatory，WENO)格式［7］同Level Set方法相结合，考虑了由于燃烧造成裂纹边界的移动，对激波绕射裂纹形成的非定常流场进行了计算，并对这种激波与移动边界相互作用的非定常流场进行了详细分析，探讨了点火期间裂纹内压力突升的主要机理。并通过对静态和动态边界条件下所形成的流场结构进行比较，探讨了引入动态边界条件计算的必要性。

1 控制方程与数值计算方法

1.1 流场控制方程

二维非定常可压缩Navier-Stokes方程可写成如下守恒形式:

其中:U=(ρ，ρu，ρv，e)T;E、F 为无粘通量，Ev、Fv为粘性通量。为了使方程封闭，补充完全气体方程p=(γ－1)ρe。流场控制方程的求解采用高精度WENO格式，具体的离散方法可参考文献［7］。

1.2 界面的追踪

Level Set的控制方程如下:

其中，Level Set函数φ如下定义:

V是界面的运动速度，定义Level Set函数的单位外法向N=∇φ/|∇φ|，那么界面运动的速度V可以由式V=rbN给出。其中，rb为由于裂纹表面燃烧造成的边界点推移速度，用rb=αpn计算。其中，p为该处的静压，α与n均为与推进剂有关的常数。

1.3 Ghost Fluid 方法

Level Set最大优点的是避免了显式的追踪运动界面，大大提高了对复杂界面形状以及界面复杂运动的处理能力，但是在求解物理方程时，需要对界面附近的点做特殊的处理。为此，引入ghost网格方法，在φ(x，y，t)＞0的地方，认为是流体的真实网格，在φ(x，y，t)＜0的地方，认为是流体的 ghost网格。在ghost网格区域，采用Ghost Fluid插值。例如，标量I通过以下方程常数扩展:

而界面上满足条件I=Iw。文中将变量从φ＞0的区域扩展到φ＜0的区域，因此在上式中取减号。由于对上述方程处理方式的不同，将会得到不同的结果，在文中，按下列方式定义ghost网格上的变量:

式中:下标G表示ghost网格上的变量，下标E表示欧拉实网格上的变量，下标w表示界面上的值，τ是与n垂直的单位向量。

1.4 边界条件和初始条件

未点燃推进剂表面采用无滑移边界条件:u=0，v=0，∂p/∂n=0，T=Tw;已点燃推进剂表面则采用如下边界条件:Vτ=0，Vn=Vb，∂T/∂n=0。其中，Vτ是壁面处的切线速度，Vn是壁面处法向速度，Vb为推进剂燃烧对裂纹的喷射速度，按式Vb=rbρs/ρ计算，其中ρs为推进剂密度。因为在整个过程中，只存在燃气流同裂纹壁面的传热，所以认为在燃气流同推进剂之间的耦合边界上温度和热流密度均连续，即:TwⅠ=TwⅡ，qwⅠ=qwⅡ。其中，Ⅰ代表固体推进剂区域，Ⅱ代表燃气流动区域。在具体计算中，首先假定耦合边界上的温度分布，对推进剂内部进行导热计算，得到推进剂内部的温度分布，然后通过计算得到边界上的热流密度，再将热流密度应用于裂纹内燃气流动的计算之中，得到耦合边界上新的温度分布，以此分布作为推进剂边界的温度分布，重复上述计算过程直到边界上温度分布收敛到一个确定值。

2 算例与结果分析

文中采用图1所示的物理模型，并根据文献［8］，给定裂纹处激波前后压力比值p2/p1=10.87，波前参数取为大气参数:p1=1.01 ×105Pa，ρ1=1.225 kg/m3，T1=288.15 K，v1=0。然后根据激波关系式计算出裂纹入口处的流场条件。裂纹的宽度取为Wc=1 mm，深度h=4 mm。超压 Cp定义为 Cp=(pmaxp2)/(0.5ρ2v22)，其中，pmax为裂纹壁面处最大静压，p2为入口静压，v2为入口气流速度，ρ2为入口气流密度。

图1 发动机装药裂纹模型

图2为不同时刻裂纹内流场的等压线图谱变化。图2(a)为正激波绕射裂纹左侧拐角时的等压线分布图，激波S1绕过拐角时产生膨胀波，导致激波在裂纹进口处发生弯曲。主激波继续向右运动，在裂纹右侧壁面发生反射，在图2(b)中可以清晰地看到激波S3和壁面反射形成的反射激波S2。其中S3向裂纹顶端运动，S2向裂纹左侧壁面运动。随着时间的推移，激波S2分裂为两部分，即图2(c)中的S2和 S4。S4也向裂纹顶端运动，并且由于S4的运动速度大于激波S3的运动速度，二者逐渐发生融合形成一道更强的激波S5，同时，S4在左侧壁面反射形成S6，如图2(d)所示。S5继续向裂纹顶端运动并在裂纹壁面反射后形成激波S7，S7向裂纹入口处运动并与S2相互作用形成复杂的流场结构，造成裂纹内的压力振荡。从图中还可以看出，在激波到达裂纹顶端后，裂纹头部的燃速明显加快，说明裂纹头部压力较大，这与文献［9］的实验事实相符，也说明在裂纹不扩展的情况下，头部的钝化是必然的。

图2 激波绕射动态裂纹等压线图

图3 静止边界条件下激波绕射裂纹等压线图

图4 不同边界条件下超压计算结果对比图

图3是静止边界条件下计算得到的裂纹内流场等压线图谱变化情况，裂纹深度仍为h=4 mm，图4是两种边界条件下计算所得的裂纹顶端最大压力处超压对比图。从图中可以看出，两种边界条件下压力开始突升的时间基本相同，但在静止边界条件下计算的超压值明显的大于移动边界条件下所得的超压值。主要原因是在移动边界条件下，随着时间的推移，裂纹空腔不断增大，使得激波强度有所减弱。由于超压随着时间的变化呈现强烈的非线性，使得这种差别是无法事先估计的，而准确估计裂纹内的超压对发动机判废准则的制定具有十分重要的意义，因此必须引入动边界进行计算。图5是无量纲时间t=0.5时两种边界条件计算所得流线图，二者形成的流场结构明显不同，在静止边界条件下，裂纹内出现一系列涡结构，而在移动边界条件下，只得到两个比较明显的大的涡结构。因此，为了得到更为真实的流场结构图谱，也必须引入动边界进行计算。同时，这种引入移动边界进行计算的界面追踪方法也可以有效耦合装药燃烧和发动机内流场的计算，为全发动机三维内流场的数值模拟提供了一种很好的研究思路。

图5 t=0.5时流场中的流线图

3 结论

通过对计算结果的分析可知，用Level Set方法配合高精度WENO格式求解激波绕射动态裂纹可得到比较理想的结果。通过对激波绕射裂纹的研究，得到了激波在裂纹内部绕射和反射的流场变化，这在帮助我们认识点火冲击波与裂纹相互作用规律的同时，也可以为发动机的实际判废准则提供参考。通过对静态和动态边界条件下所形成的流场结构进行比较可知，不同的边界条件所得到的计算结果相差甚远，边界条件对流场结构的影响是无法预计的，要得到更加接近实际的裂纹内流场，必须考虑裂纹由于燃烧而产生的区域扩展。同时由于Level Set方法对复杂界面的追踪能力，也为进一步考虑裂纹的扩展提供了一种可能的发展方向。

［1］Kuo K K，Morcei J.Crack propagation and branching in burning solid propellants［C］∥Twenty-first Symposiums on Combustion/The Combustion Institute，1986:1933－1941.

［2］Lu YC，Kuo K K.Modeling and numerical simulation of combustion process inside a solid propellant crack［J］.Propellants Explosives Pyrotechnics，1994(19):217－226.

［3］崔小强，白晓征，周伟，等.点火冲击波与药柱裂纹流固耦合作用机理数值模拟研究［J］.固体火箭技术，2010，33(1):9－12.

［4］赵汝岩，于胜春，李昊，等.点火升压阶段药柱裂纹变形研究［J］.固体火箭技术，2009，32(1):43－47.

［5］徐春光，刘君.激波在狭缝中绕射过程研究［J］.固体火箭技术，2003，26(3):39－41.

［6］刘儒勋，王志峰.数值模拟方法和运动界面追踪［M］.合肥:中国科学技术大学出版社，2001.

［7］Jiang G S，Shu C W.Efficient implementation of weighted ENO schemes［J］.Journal of Computational Physics，1996，126:202－228.

［8］张文普.固体推进剂裂纹内燃烧流动规律研究［D］.西安:西北工业大学，2000.

［9］葛爱学.固体火箭发动机点火过程与装药裂纹相互作用机理研究［D］.长沙:国防科学技术大学，2004.