膛口流场是非定常、带有强激波间断的复杂流场.在膛内火药燃气推动作用下，加速运动的弹丸不断压缩弹前空气，形成弹前冲击波和初始流场^[1-2].在火炮射击的后效期，高温高压燃气迅速释放，进入初始流场，膛口周围出现各种复杂的物理和化学现象，弹丸在火药燃气推力作用下继续加速，膛口附近形成了复杂的瞬态膛口流场波系结构，并在远场产生冲击波和噪声^[3-4].大口径高超速平衡炮炮膛口流场的有害扰动是极其强烈的，而且对弹丸的稳定飞行影响较大，但通过试验直接测量膛口流场的特性参数非常困难，因此，通过数值模拟分析大口径高超速平衡炮的膛口流场特性，是一个经济有效的方法^[5].

目前，国内外已经对膛口流场进行了相关的试验和研究^[6-8].Jiang等^[9]用基于ALE方程及动网格技术对44 mm口径弹丸由膛内到膛外直到飞离初始流场的整个过程进行数值模拟.江坤等^[10]利用结构网格对122 mm车载火炮的膛口流场进行数值模拟，对初速为713 m/s的弹丸发射过程中的流动现象进行了分析，得到了较为可信的结果.代淑兰等^[11]采用三维非定常化学反应流控制方程组对带制退器的7.62 mm口径、735 m/s初速的弹丸膛口燃烧流场进行了数值模拟，清晰的描述了膛口流场的发展过程结构和性质及弹丸与流场的相互影响.Love^[12]通过实验和理论方法对轴对称自由射流进行了研究，分析了射流马赫数、喷管形状等因素对轴线上的膛口气流参数分布的影响.高树滋等^[13]采用MUSCL差分格式对膛口流场进行了分析, 得到了由空气冲击波和射流近场结构所组成的全流场流动图像.乐贵高等^[14]利用Osher格式及非结构网格对二、三维可压缩欧拉流动方程组对某型号火炮膛口气流现象进行分析和计算.

上述所做的研究均基于155 mm口径以下的火炮，多以常规初速或低速发射弹丸为研究对象.由于弹丸初速相对较低，多以膛口附近处开始计算，并将膛口压力或弹丸速度简化为定值.而对于300 mm大口径超高速弹丸的初始流场及火药燃气流场的射流结构的分析未见报道.本文建立该发射条件下膛口流场的二维轴对称数值仿真模型，耦合内弹道过程，研究超高速弹丸自膛底飞出膛口时的初始流场和火药燃气流场的射流结构，并对结果进行了分析，为大口径高超速火炮的设计与试验提供理论指导，为深入研究奠定良好基础.

1 数学模型 1.1 控制方程

当不考虑外加热和彻体力的影响时，笛卡尔坐标系下的二维轴对称可压缩非定常的N-S方程组为^[15]

$ \frac{{\partial \mathit{\boldsymbol{Q}}}}{{\partial t}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial x}} + \frac{{\partial \mathit{\boldsymbol{G}}}}{{\partial y}} = \frac{\sigma }{y}\mathit{\boldsymbol{S}}. $

式中：Q为守恒变矢量；F、G分别为坐标方向的通量；S为轴对称源项，具体表达式为：

$ \mathit{\boldsymbol{Q}} = \left[{\begin{array}{*{20}{l}} \rho \\ {\rho u}\\ {\rho v}\\ {\rho e} \end{array}} \right], $

$ \mathit{\boldsymbol{F}} = \left[{\begin{array}{*{20}{l}} {\rho u}\\ {\rho {u^2} + p-{\tau _{xx}}}\\ {\rho uv-{\tau _{xy}}}\\ {\left( {\rho e + p} \right)u-u{\tau _{xx}} - v{\tau _{xy}} + {q_x}} \end{array}} \right], $

$ \mathit{\boldsymbol{G}} = \left[{\begin{array}{*{20}{l}} {\rho v}\\ {\rho uv-{\tau _{xy}}}\\ {\rho {v^2} + p-{\tau _{yy}}}\\ {\left( {\rho e + p} \right)v-u{\tau _{xy}} - v{\tau _{yy}} + {q_y}} \end{array}} \right], $

$ \mathit{\boldsymbol{S}} = \left[{\begin{array}{*{20}{l}} {\rho v}\\ {\rho uv-{\tau _{xy}}}\\ {\rho {v^2} + {\tau _{yy}}-{\tau _{\theta \theta }}}\\ {\rho \left( {e + p} \right)v-u{\tau _{xy}} + v{\tau _{yy}} + {q_y}} \end{array}} \right], $

其中应力项为：

$ \begin{array}{l} {\tau _{\theta \theta }} = \frac{2}{3}\mu \left( {2\frac{v}{y} - \frac{{\partial u}}{{\partial x}} - \frac{{\partial v}}{{\partial y}}} \right), \\ {\tau _{xx}} = \frac{2}{3}\mu \left( {2\frac{{\partial u}}{{\partial x}} - \frac{{\partial v}}{{\partial y}} - \frac{v}{y}} \right), \\ {\tau _{yy}} = \frac{2}{3}\mu \left( {2\frac{{\partial v}}{{\partial y}} - \frac{{\partial u}}{{\partial x}} - \frac{v}{y}} \right), \\ {\tau _{xy}} = {\tau _{yx}} = \mu \left( {\frac{{\partial u}}{{\partial y}} + \frac{{\partial v}}{{\partial x}}} \right), \end{array} $

压力由理想气体方程给出，即

$ p = \left( {\gamma - 1} \right)\left[{\rho e-\frac{\rho }{2}\left( {{u^2} + {v^2}} \right)} \right]. $

式中：ρ为气体密度；u、v分别为x、y方向的速度分量；e为总能量; γ为气体比热比；μ为层流黏性系数；τ_θθ、τ_xx、τ_yx、τ_yy分别为不同方向的黏性力；k为热导率; q_x、q_y分别为单位质量的体积加热率; 系数σ决定流动的类型，当σ=1时，为二维轴对称流动模型，当σ=0时，为二维平面流动模型.

1.2 湍流模型

本文所利用的湍流模型为Realizable k-ε模型.该模型比起标准k-ε模型有两个主要的不同点：1)Realizable k-ε模型为湍流黏性增加了一个公式；2)为耗散率增加了新的传输方程.

引入Boussinesq的线性涡黏假设，雷诺应力表达式为

$ \tau _{ij}^R = - \overline {\rho {{u'}_i}{{u'}_j}} = {\mu _t}\left( {\frac{{\partial {U_i}}}{{\partial {x_j}}} + \frac{{\partial {U_j}}}{{\partial {x_i}}} - \frac{2}{3}\frac{{\partial {U_i}}}{{\partial {x_i}}}} \right) - \frac{2}{3}\rho k{\delta _{ij}}. $

不同的涡黏模型，其涡黏系数μ_t不尽相同，k-ε模型中$ {\mu _t} = f\left( {\frac{{\rho {k^2}}}{\varepsilon }} \right) $，其中：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;k = \frac{{\overline {{{u'}_i}{{u'}_j}} }}{2}, \\ \varepsilon = v\overline {\frac{{\partial {{u'}_i}}}{{\partial {x_j}}}\left( {\frac{{\partial {{u'}_i}}}{{\partial {x_j}}} + \frac{{\partial {{u'}_j}}}{{\partial {x_i}}}} \right)} . \end{array} $

Realizable k-ε模型中法相正应力$ \overline {{{u'}_i}{{u'}_j}} > 0 $，关于Reynolds剪切应力的Schwarz不等式：

$ \overline {{{\left( {{{u'}_i}{{u'}_j}} \right)}^2}} \le \overline {u_i^2u_j^2} . $

2 动网格与边界条件 2.1 动网格方法

为了处理计算过程中由于弹丸运动所引起的网格变化，需要采用动网格技术.动网格算法主要用来计算内部网格节点的调节，有3种算法：铺层法、弹性光顺法和局部重构法.由于本文中考虑到弹丸只沿一个方向做平移运动，所以采用铺层法即可实现对弹丸运动的模拟.铺层法中需要对网格的分割因子和合并因子进行设定.由于弹丸的运动，使得弹后网格随弹丸沿x轴向前平移，靠近膛底的一层网格被拉长，网格边长为h_s, 网格超出计算域边界部分被截断，靠近计算域边界的一层网格长度为h_c，假设给定的理想网格尺寸为h_i, 当新生成的网格尺寸满足h≥(1+c_s)h_i时，网格被分割成两个，当网格尺寸满足h≥c_ch_i时，网格将合并为一个网格.本算例中取h_i=5 mm，c_s=0.4，c_i=0.1.

2.2 网格划分与边界条件

对复杂流场进行模拟时，很多情况下很难生成单块高质量结构网格，此时可以采用分区拼接网格的方法进行处理.将整个计算域分为3个区域：弹后Ⅰ区、弹前Ⅱ区、膛口周围流场Ⅲ区，炮管总长29.6 m，管壁厚度为50 mm，如图 1(a)所示.初始时，弹丸距膛底100 mm，膛口周围流场区域长13 m，宽3 m.最小网格尺寸为2 mm，为避免非结构网格带来过多的计算量与格式耗散，本文均采用结构化网格进行划分，共28万个网格.为验证本文所采用的网格划分方法及网格尺寸合理，以总网格数35万模型作为对比，两种情况下弹底压力最大误差不超过1%，证明结论与网格划分无关.

图 1(b)给出了镜像后的初始网格划分模型.计算域中，设定身管、膛口为固壁边界条件，弹丸为移动固壁边界，弹丸在膛内运动阶段由内弹道UDF程序控制，出膛口后由六自由度UDF程序控制其运动，膛底设为压力入口边界条件，计算域边界为压力出口边界条件.壁面均假定为绝热，壁面上温度为临近壁面网格点温度.

2.3 计算模型

本文模型中，弹丸质量为160 kg，平衡体重2 000 kg，总装药量为240 kg，图 2(a)、(b)分别给出了此平衡炮内弹道过程的弹底压力随时间的变化曲线和速度随时间的变化曲线.由计算结果得出，弹丸出膛口速度为1 730.4 m/s，膛口压力为46.34 MPa, 运动时间为37.6 ms^[16-17].弹丸后效期的膛口压力随时间呈指数分布^[18-19]，其压力分布如图 2(c)所示.身管外流场区域始压力为101 325 Pa，初始温度为300 K.

3 结果与分析 3.1 数值验证

本文主要对300 mm大口径超高速平衡炮的发射过程进行了数值模拟，采用分块网格划分的整体运动处理方法及Realizable k-ε湍流模型，耦合内弹道过程及六自由度程序控制弹丸运动.为了验证本文算法的可靠性，模拟以往研究采用分层动网格技术结合Roe的迎风格式得到的某48 mm口径530.0 m/s枪炮的膛口流场特性^[20]，并得到了一致的结果，图 3、4给出了膛口流场压力等值线图与速度矢量图的对比，证明了本文算法的正确性.

3.2 初始流场分析

根据耦合内弹道过程的膛口流场数值计算，得到了大口径超高速平衡炮初始流场结构分布.弹丸在膛内推动弹前的空气柱，形成一系列压缩波，压缩波不断向膛口方向传播，由于后一个压缩波的传播速度总比前一个大，这些压缩波相互叠加形成了激波，激波在膛内逐渐加强，随着弹丸的加速运动，弹前空气被挤出膛口形成初始流场.典型的膛口流场的波系结构如图 5所示.

图 6(a)~(h)显示了超高速弹丸初始流场的形成过程.当弹丸运动至34.3 ms时，初始流场即将形成，此时弹丸距膛口5.5 m，弹丸速度为1 658.3 m/s.初始冲击波的强度主要由弹丸初速的大小决定，弹丸初速越高，弹前激波及初始冲击波就越强.当弹丸初速在1 000 m/s左右时，弹前激波到达膛口时的气体压力一般为1~2 MPa，而由图 6(a)可以看出，当弹丸以1 730 m/s超高速运动时，此压力达到了3.8 MPa左右，高于中小口径中低速发射条件下的1倍.压缩波在膛内形成轴对称涡流环，并随气流一起做轴向运动.由于波后的压力比外界压力高的多，于是形成了膨胀不足的初始流场.随着弹丸继续向前运动，初始流场逐渐形成并发展.

当弹前空气被压缩流出膛口后，气体迅速膨胀，速度增大，压力降低，并向周围传播，在身管固壁处形成涡流，涡流所在位置出现负压.由图 6(d)所示，在弹丸距膛口3.3 m处，已形成了清晰的初始流场.图 6(f)可以较清楚的观察到初始冲击波、初始瓶状激波、马赫盘以及三波点.膛口附近的涡流也不断变大并向侧前方运动.由图 6(g)可以看出，弹丸的圆柱部已运动出膛口，使膛内的弹前气体全部压出，气体不断膨胀向四周传播.被压缩的空气向前传播使得后方压力降低，周围空气流向该区域，并使弹丸周围涡流继续向前运动.

3.3 火药气体流场分析

随着弹丸运动出膛口，耦合后效期过程及弹丸六自由度运动方程，得到了火药燃气流场的分布状态.图 7(a)~(g)给出了膛口流场的波系结构及其形成过程.由图 7(a)所示，大口径超高速弹丸膛口冲击波异常强烈，高温高压火药气体首先从弹尾周围逸出，向侧前方迅速发展传播，速度高达2 500 m/s以上，而中小口径中低速发射状态火药燃气速度仅为1 500~1 800 m/s.由火药气体的速度矢量图(图 8)可以清晰的看到火药气体的运动方向.气体推动膛口周围已被初始流场扰动过的空气，形成了膛口冲击波，由图 7(b)可以看到清晰的膛口冲击波的形成.

火药气体向前运动追赶弹丸，由于弹丸初速高，2 500 m/s速度的火药气体不能包围弹丸.而对于中小口径中低速弹丸，由于弹丸弹身短、口径小、初速低，从膛口喷出的高温高压气体相对于弹丸初速高出几倍，很快将弹丸包围.图 7(c)、(d)可以看到随着弹丸的运动逐渐摆脱了火药气体对其的包围.图 7(d)可以清晰的看到弹底激波的形成，这是因为喷射出的火药气体射流大于弹丸的速度.由于初始流场主漩涡的作用，膛口冲击波波阵面在主漩涡环的中心处出现褶皱，弹底激波在膛口冲击波的作用下进一步加强并阻碍马赫盘的生长，随着弹底激波的作用越来越弱，马赫盘逐渐增大，瓶状激波不断扩大，膛口冲击波不断向远处传播并衰减.但膛口冲击波在不同方向上的衰减规律不同，膛口冲击波前方的强度要大于侧后方的强度.由图 7(e)所示，弹丸侧面形成了起始于弹丸头部并脱离弹丸侧面向外延伸的剪切层，随着火药燃气持续喷出空气接触发生复杂的变化，射流的结构更加明显.

图 7(f)可以清晰的看到膛口冲击波的三波点，膛口形成了由冠状冲击波、弹底激波、反射激波、马赫盘构成的多层次激波、间断面相互叠加的完整波系.图 9(a)、(b)给出了弹丸出膛口后弹底压力及Ma随时间的变化情况，由图 9可以看出，在1.0 ms以后火药燃气作用在弹底的压力基本为零，弹丸Ma_max可以达到4.0以上.

4 结论

1) 超高速运动的弹丸在膛内距膛口5.5 m处初始流场开始形成，距膛口3.3 m处已形成了清晰的初始流场, 且冲击波到达膛口时压力高于中小口径中低速发射状态下1倍.

2) 火药燃气流场异常强烈，火药燃气速度达到2 500 m/s以上，其波系结构更加明显，膛口射流结构更加剧烈.

3) 弹丸以1 730 m/s初速运动时弹后气体不能追赶并包围弹丸，故对弹丸运动的影响不同于中小口径中低速发射状态.