1 引言

对膛口流场的研究涉及气体动力学、气动声学等多个学科, 因此膛口流场的机理研究是一个具有挑战性的课题。在膛内火药燃气推动作用下, 加速运动的弹丸不断压缩弹前空气, 形成初始流场^[1-2]。高温高压火药燃气在弹丸出炮口后迅速喷出, 受到弹丸的阻碍并与初始流场融合, 在膛口处形成具有复杂波系的膛口流场^[3-4]。大口径高超速火炮膛口流场的有害扰动极其强烈, 但由于炮口焰、炮口烟、电磁干扰等的存在, 使得瞬态变化的膛口流场的激波结构很难测量和计算^[5]。

人们对膛口流场形成的发展机理、膛口冲击波场的分布规律、膛口气流对弹丸的后效作用等方面都有了较为清晰的认识^[6]。Love^[7]等通过实验和理论方法对轴对称自由射流进行了研究, 分析了射流马赫数、喷管形状等因素对轴线上的膛口气流参数分布的影响。李鸿志^[8]等分析了膛口冲击波与射流的特征及形成机理, 提出了膛口冲击波的变能量、动球心的球形物理模型, 并总结出膛口冲击波远场传播规律; 姜孝海^[9]等用基于ALE方程及动网格技术对44 mm口径弹丸由膛内到膛外直到飞离初始流场的整个过程进行数值模拟。江坤^[10]利用结构网格对122 mm车载火炮的膛口流场进行数值模拟, 对初速为713 m·s^-1的弹丸发射过程中的流动现象进行了分析, 得到了较为可信的结果。代淑兰^[11]采用三维非定常化学反应流控制方程组对带制退器的7.62 mm口径、735 m·s^-1初速的弹丸膛口燃烧流场进行了数值模拟, 清晰的描述了膛口流场的发展过程结构和性质及弹丸与流场的相互影响。郭泽庆^[12]等以7.62 mm弹道枪为研究对象, 基于ALE方程及分区结构贴体网格对膛口初始流场及后效期过程进行了试验分析。Jiang^[13]基于二维结构网格分别用TVD格式和耗散控制格式离散Euler方程, 对轴对称膛口流场进行了数值模拟, 考虑了圆柱形弹丸对流场的影响, 并分析了弹丸与膛口气流间的相互作用。张焕好^[14]基于三维Euler方程, 结合Roe格式与结构化动网格技术对含有制退器的20 mm枪炮膛口流场进行了数值模拟。

前人所研究的重点多为燃气流场的特性分析, 且多以155 mm口径以下, 1000 m·s^-1以下初速的枪炮为研究对象。而对于300 mm大口径1730 m·s^-1超高速弹丸的初始流场对火药燃气射流结构的影响分析未见研究。本研究采用有限体积法并结合Realizable k-ε湍流模型, 使用结构动网格技术, 建立两种二维轴对称数值仿真模型, 研究超高速弹丸含初始流场和简化的不含初始流场情况下的火药燃气流场结构特性, 得到了该发射条件下初始流场对火药燃气流场分布的影响, 为大口径高超速火炮的设计与试验提供理论指导, 为深入研究奠定良好基础。

2 数值方法 2.1 控制方程

当不考虑外加热和彻体力的影响时, 二维轴对称可压缩非定常N-S方程组为^[15-16]:

$ \frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial {\rm{t}}}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial x}} + \frac{{\partial \mathit{\boldsymbol{G}}}}{{\partial y}} + \mathit{\boldsymbol{Q}} = \frac{1}{{Re}}(\frac{{\partial {\mathit{\boldsymbol{F}}_v}}}{{\partial x}} + \frac{1}{{{y^\sigma }}}\frac{{\partial {y^\sigma }{\mathit{\boldsymbol{G}}_v}}}{{\partial y}} + {\mathit{\boldsymbol{Q}}_v}) $

(1)

$ \frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial {\rm{t}}}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial x}} + \frac{{\partial \mathit{\boldsymbol{G}}}}{{\partial y}} + \mathit{\boldsymbol{Q}} = \frac{1}{{Re}}(\frac{{\partial {\mathit{\boldsymbol{F}}_v}}}{{\partial x}} + \frac{1}{{{y^\sigma }}}\frac{{\partial {y^\sigma }{\mathit{\boldsymbol{G}}_v}}}{{\partial y}} + {\mathit{\boldsymbol{Q}}_v}) $

(2)

$ \mathit{\boldsymbol{U}} = {\left[{\rho \;\;\rho u\;\;\rho v\;\;\rho e} \right]^T} $

(3)

$ \mathit{\boldsymbol{F}} = {\left[{\rho u\;\;\rho {u^2} + p\;\;\rho uv\;\;\rho e + pu} \right]^T} $

(4)

$ {\mathit{\boldsymbol{F}}_v} = {\left[{0\;\;{\tau _{xx}}\;\;{\tau _{xy}}\;\;u{\tau _{xx}} + v{\tau _{xy}} + {q_x}} \right]^T} $

(5)

$ \mathit{\boldsymbol{G}} = {\left[{\rho v\;\;\rho uv\;\;\rho {v^2} + p\;\;\rho e + pv} \right]^T} $

(6)

$ {\mathit{\boldsymbol{G}}_v} = {\left[{0\;\;{\tau _{yx}}\;\;{\tau _{xy}}\;\;u{\tau _{yx}} + v{\tau _{yy}} + {q_y}} \right]^T} $

(7)

$ \mathit{\boldsymbol{Q}} = \frac{\sigma }{y}{\left[{\rho v\;\;\rho uv + p\;\;\rho {v^2}\;ev + pv} \right]^T} $

(8)

$ {\mathit{\boldsymbol{Q}}_v} = \frac{\sigma }{y}\left[{0\;\;0\;\;-{\tau _{\theta \theta }}\;\;0} \right]{^T} $

(9)

其中应力项为:

$ \left. \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} \right\} $

(10)

系数σ决定流动的类型, 当σ=1时, 为二维轴对称流动模型, 当σ=0时, 为二维平面流动模型^[11]。

其中, ρ为气体密度, kg·m^-3; u、v分别为x、y方向的速度分量, m·s^-1; e为总能量, kJ, 其表达式为: $ e = \frac{p}{{\gamma-1}} + \frac{1}{2}\rho \left( {{u^2} + {v^2}} \right) $, γ为气体比热比。Re为雷诺数, τ_θθ、τ_xx、τ_yx、τ_yy分别为不同方向的粘性力, N; μ为粘性系数, 有层流粘性系数和湍流粘性系数组成。q_x、q_y为单位质量的体积加热率。

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

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

(11)

2.2 湍流模型

本研究所采用的湍流模型为Realizable k-ε模型^[17]。该模型比起标准k-ε模型有两个主要的不同点: ① Realizable k-ε模型为湍流粘性增加了一个公式; ② 为耗散率增加了新的传输方程。

引入Boussinesq的线性涡粘假设, 雷诺应力表达式为:

$ \begin{array}{l} \tau _{_{ij}}^{^R} =-\rho {u^\prime }_i{u^\prime }_j\\ = {\mu _t}(\frac{{\partial {\mathit{\boldsymbol{U}}_i}}}{{\partial {x_j}}} + \frac{{\partial {\mathit{\boldsymbol{U}}_j}}}{{\partial {x_i}}}-\frac{2}{3}\frac{{\partial {\mathit{\boldsymbol{U}}_i}}}{{\partial {x_i}}})-\frac{2}{3}\rho k{\delta _{ij}} \end{array} $

(12)

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

其中:

$ k = \overline {\frac{{{u^\prime }_i{u^\prime }_j}}{2}} $

(13)

$ \varepsilon = v\overline {\frac{{\partial {u^\prime }_i}}{{\partial {x_j}}}(\frac{{\partial {u^\prime }_i}}{{\partial {x_j}}} + \frac{{\partial {u^\prime }_j}}{{\partial {x_i}}})} $

(14)

其中Reynolds剪切应力的Schwarz不等式:

$ \left( {\overline {{u^\prime }_i{u^\prime }_j} } \right){^2} \le \overline {u_{_i}^{^2}u_{_j}^{^2}} $

(15)

3 边界条件与计算模型 3.1 网格划分

本研究涉及到两种计算模型, 第一种计算域不考虑初始流场的影响, 为弹丸从炮口处开始运动, 如图 1a所示。将整个计算域分为三个区域:弹后Ⅰ区、弹前Ⅱ区、及膛口周围流场Ⅲ区。膛口周围流场区域长13 m, 宽3 m。最小网格尺寸为2 mm, 位于膛口附近。计算域边界处网格相对稀疏。均采用结构化网格进行划分, 共10万个网格, 如图 1b所示。

第二种计算域为弹丸自膛底运动出炮口整个过程, 炮管总长29.6 m, 管壁厚度为50 mm, 如图 1c所示。膛口周围流场区域与第一种模型相同, 长13 m, 宽3 m。最小网格尺寸为2 mm, 位于膛口附近。计算域边界处网格相对稀疏。共28万个网格, 如图 1d所示。

3.2 边界条件

不考虑初始流场的计算模型中, 设定膛口处为压力入口边界条件, 使用自编的UDF程序对入口压力进行设定, 膛口周围流场边界为压力出口边界条件。设定身管、炮口为固壁边界条件, 弹丸运动由六自由度程序控制, 设定对称轴为轴边界。壁面均假定为绝热, 壁面上温度为临近壁面网格点温度。身管外流场区域初始压力为1.01×10⁵ MPa, 初始温度为300 K。

考虑初始流场的计算模型中, 弹丸在膛内运动阶段由内弹道程序控制, 出膛口后由六自由度程序控制, 设定膛底为压力入口边界条件, 设定身管、炮口为固壁边界条件, 弹丸为移动固壁边界条件。

3.3 内弹道与后效期

以300 mm大口径超高速平衡炮为研究对象, 弹丸质量为160 kg, 总装药量为240 kg, 平衡体重2000 kg, 图 2a、图 2b分别给出了内弹道过程的弹底压力和速度随时间的变化规律。弹丸炮口初度为1730.4 m·s^-1, 炮口压力为46.34 MPa, 运动时间为37.6 ms。

针对本算例, 后效期炮口压力随时间变化公式如下^[18-20]: p=46.3×10⁶e^-24.886t

弹丸后效期的炮口压力随时间呈指数分布规律, 其压力分布如图 2c所示。

3.4 网格无关性验证

对于网格收敛性测试, 本研究采用网格数为六十万的无初始流场模型作为对比, 以弹底压力为参考值, 如图 3所示, 当两种不同网格数量下弹底压力相差很小时, 可以认为计算结果与网格无关, 由图 3可以看出, 两者误差小于1%, 由此可以认为本计算结果与网格无关。

4 结果分析 4.1 初始流场对燃气射流的影响分析

模拟结果基于对有无初始流场的对比上, 图 4a~图 4d显示了高超音速弹丸初始流场的形成过程。弹丸在膛内推动弹前空气, 形成了一系列压缩波, 压缩波不断向炮口方向传播。随着弹丸的加速运动, 弹前空气被挤出炮口形成初始流场。当弹前空气被压缩流出膛口后, 气体迅速膨胀, 压力降低, 速度增大, 并向周围传播。弹丸的圆柱部运动出炮口时, 弹前气体全部被压出, 被压缩的空气向前传播使得后方压力降低, 周围空气流向该区域。

随着弹丸运动出炮口, 弹后高温高压火药气体迅速喷出, 形成了复杂而强烈的火药气体射流结构, 图 5a~图 5d和图 6a~图 6d分别给出了无初始流场和有初始流场情况下的火药燃气射流波系结构。

高温高压火药气体首先从弹尾周围逸出, 向侧前方迅速发展传播, 由图 5a~图 5d所示, 当无初始流场扰动时, 火药燃气流场波阵面近似为球形, 火药气体不能追赶并超越弹丸, 火药燃气刚喷出时速度约为2169 m·s^-1。弹丸侧面形成了明显的起始于弹丸头部脱离弹丸侧面向外延伸的剪切层。因为炮口喷射出的气体速度大于弹丸的运动速度, 弹丸底部形成了清晰的弹底激波。弹底激波在膛口冲击波的作用下进一步加强并阻碍马赫盘的生成, 随着弹底激波的作用越来越弱, 马赫盘逐渐增大, 瓶状激波不断扩大, 膛口冲击波不断向远处传播并衰减。炮口形成了由冠状冲击波、弹底激波、马赫盘等构成的完整波系。

由图 6a~图 6d可以看出有初始流场的火药燃气流场结构与无初始流场影响下的燃气流场结构明显不同。火药燃气推动炮口周围已被初始流场扰动过的空气, 形成了炮口冲击波。在炮口初始流场的影响下, 火药燃气追赶并超越弹丸, 在弹体周围形成了复杂的火药燃气团, 火药燃气刚喷出膛口时速度约为2359 m·s^-1, 比无初始流场影响下的速度提高约200 m·s^-1。随着弹丸的运动, 弹丸逐渐摆脱了火药气体对其的包围。随着弹丸的运动, 弹底激波的作用越来越弱, 马赫盘逐渐增大, 瓶状激波不断扩大, 膛口冲击波不断向远处传播并衰减, 同样形成了由冠状冲击波、弹底激波、马赫盘等构成的完整波系。

在初始流场的干扰下, 火药燃气流场的各参数也明显不同。图 7给出了计算模型中监测点的分布状况, 其坐标分布见表 1。为了研究无初始流场干扰下和有初始流场干扰下的火药燃气流场参数变化, 在弹前膛口周围处监测了四个点的压力和温度变化, 表 2给出了监测点在有初始流场和无初始流场情况下的最大压力和最大温度。

图 8a~图 8d和图 9a~图 9d给出了有无初始流场情况下的监测点压力和温度的对比曲线。由图 8可以看出, 四个监测点的最大压力随膛口距离的增加而减小。对于监测点1, 有初始流场扰动情况下的初始压力比无初始流场扰动情况下的初始压力低, 而监测点2、3正好相反。这是因为监测点1位于膨胀波影响的区域内, 膨胀波的影响使得此处的压力低于环境压力。而监测点2、3位于压缩波影响的区域内, 压缩波使得此处压力高于环境压力。而此时监测点4位于初始流场没有影响到的区域内, 所以有无初始流场扰动情况下的初始压力相同。值得注意的是, 有初始流场扰动情况下的最大压力低于无初始流场扰动情况下的最大压力, 这是因为初始流场推动了膛口前方空气, 形成了负压, 因此当火药燃气射流传播至此时, 有初始流场扰动的最大压力低于无初始流场扰动下的最大压力。

图 10为初始流场的温度等值线图。由图 10可以看出, 经初始流场扰动过的弹前空气的温度分布有了明显的变化, 无初始流场下的膛外温度为室温300 K, 而有初始流场的膛外温度一般在1000~1500 K左右。因此, 经初始流场扰动过的膛口前方温度高于没有初始流场扰动过的温度。

4.2 初始流场对弹丸运动的影响

由于初始流场的存在, 使得弹丸出炮口前的流场区域已经受到复杂的干扰, 火药燃气对弹丸底部的作用也受到了影响。图 11a和图 11b分别给出了无初始流场和有初始流场情况下的弹底压力和马赫数分布情况。

由图 11可以看出, 无初始流场时后效期弹底最大压力为40.3 MPa, 有初始流场情况下为38.9 MPa。无初始流场情况下, 弹底马赫数最高达到4.13, 在初始流场的影响下, 弹底马赫数出现波动, 最大马赫数为4.09。

初始流场波系强度较弱, 与火药燃气流场相差一个数量级, 但初始流场的存在影响了火药燃气对弹丸的作用, 使得弹底压力减小, 马赫数减小。

5 结论

针对300 mm大口径1730 m·s^-1超高速平衡炮, 建立两种数学模型, 探讨了无初始流场和有初始流场下火药燃气流场的异同, 以及初始流场的存在对弹丸运动的影响。本研究编写内弹道程序控制弹丸在膛内的运动, 编写后效期程序对炮口压力进行设定、编写六自由度程序控制弹丸在膛外的运动。

(1) 当没有初始流场影响时, 火药燃气射流呈球状, 燃气不能追赶上弹丸。而在初始流场的影响下, 火药燃气可以追赶上弹丸, 但由于初速高, 并不能超越并包围弹丸。

(2) 初始流场的存在使得火药燃气刚喷出膛口时的速度增高。

(3) 初始流场的干扰使得燃气射流的最大压力降低近一倍, 温度增加1000 K以上, 弹底压力和弹丸马赫数稍有减小。