在过去的许多年中, 人们针对超声速气流中的氢气-空气燃烧进行了大量的研究, 其目的是了解超燃发动机及利用超燃燃烧的实验设备(如激波管)中的流动机理[1].超声速燃烧流场中流动特征时间和化学反应特征时间同属一个量级.而且, 带化学反应的超声速流场中有一种很典型的物理现象:激波诱导燃烧, 如图 1所示.激波后的高温加热预混气体.气体通过激波后并不马上燃烧起来, 而是要经过一段距离, 能量才开始释放出来.激波后到能量开始大量释放的这段距离称为诱导区.在诱导区, 温度、压力和密度近似保持不变, 但是会发生一些很重要的基元反应, 中间产物的浓度呈指数增加.诱导区之后, 储存在混合物中的能量开始释放, 温度迅速上升, 这个区域称为能量释放区.诱导区的厚度主要由诱导时间决定.激波后的温度很高、压力很大时, 由于反应诱导时间很短, 因而诱导距离也很短; 反之亦然. Lehr针对直径15 mm弹头氢气-空气爆轰开展了一系列实验研究[2].
|
| 图 1 氢气-空气激波诱导燃烧示意图 Fig.1 Sketch of shock-induced combustion by hydrogen-air |
随着计算机技术的发展, 数值方法成了模拟这种现象的主要手段[3]. Yungster等[4]模拟了Lehr实验状态中的Ma=4.18, 5.11和6.46来捕捉激波和燃烧波. Wilson等[5]也对激波诱导燃烧进行了详尽的数值模拟, 研究表明, 化学动力学模型和网格精度对模拟这种爆轰流动现象非常重要. Matsuo等[6]研究了头部半径导致的不稳定性周期性变化.
由于来流Mach数和弹头直径都是影响爆轰燃烧流场的关键参数, 因此, 本文采用数值方法对Lehr实验中的氢气-空气爆轰进行了模拟, 研究了Ma=4.18, 5.11, 6.46流动时的燃烧结构, 以及弹头半径对氢气-空气燃烧爆轰流动稳定性的影响.
1 控制方程和方法 1.1 控制方程假定化学非平衡流中的气体介质是热完全气体的混合物, 不考虑振动能的激发, 流动质量扩散采用双组元气体模型假设.在直角坐标系下, 带化学非平衡反应的三维N-S方程可以写为
| $ \frac{\partial \boldsymbol{U}}{\partial t}+\frac{\partial \boldsymbol{E}}{\partial x}+\frac{\partial \boldsymbol{F}}{\partial y}+\frac{\partial \boldsymbol{G}}{\partial z}=\frac{\partial \boldsymbol{E}_{v}}{\partial x}+\frac{\partial \boldsymbol{F}_{v}}{\partial y}+\frac{\partial \boldsymbol{G}_{v}}{\partial z}+\boldsymbol{S} $ |
其中, 守恒变量U, 对流通量E, F和G, 耗散通量Eν, Fν和Gν以及化学反应源项S的具体形式如文献[7], [8]所示.
1.2 化学及热力学模型第i种组分的化学非平衡生成源项可表示为[9]
| $ \dot{\omega}_{i}=W_{i} \sum\limits_{r=1}^{n r}\left(\nu_{i r}^{b}-\nu_{i r}^{f}\right)\left(\frac{\rho}{\bar{W}}\right)^{n}\left[ {{R}_{\text{f}\mathit{r}}}-{{R}_{\text{b}\mathit{r}}} \right] $ |
其中, Rfr, Rbr分别代表以摩尔浓度表示的第r个基元反应的正向反应和逆向反应速率(mole/(m3·s)).
| $ \left. \begin{array}{*{35}{l}} {{R}_{\text{f}\mathit{r}}}=k_{r}^{\text{f}}\prod\limits_{i=1}^{ns}{{{\left( {{\chi }_{i}}\rho \right)}^{\eta _{ir}^{\text{f}}}}} \\ {{R}_{\text{b}\mathit{r}}}=k_{r}^{\text{b}}\prod\limits_{i=1}^{ns}{{{\left( {{\chi }_{i}}\rho \right)}^{\eta _{ir}^{\text{b}}}}} \\ \end{array} \right\}\overset{当反应速率公式以{\rm{cgs}}为单位}{\mathop{\Rightarrow }}\, $ |
| $ \left\{ \begin{array}{*{35}{l}} {{R}_{\text{f}\mathit{r}}}=1000\left[ k_{r}^{\text{f}}\prod\limits_{i=1}^{ns}{{{\left( 0.001{{\chi }_{i}}\rho \right)}^{\left. \eta _{ir}^{\text{f}} \right]}}} \right. \\ {{R}_{\text{b}\mathit{r}}}=1000\left[ k_{r}^{\text{b}}\prod\limits_{i=1}^{ns}{{{\left( 0.001{{\chi }_{i}}\rho \right)}^{\eta _{ir}^{\text{b}}}}} \right] \\ \end{array} \right. $ |
氢氧燃烧机理采用文献[10]中的7组分(N2, O2, H2, OH, H, O, H2O)8步反应模型.
1.3 计算方法采用有限体积求解器进行方程求解.以三维情形为例, 在所研究的控制体积内, 将含化学反应源项的基本方程以积分形式表示为[11]
| $ \iiint\limits_{\Delta V}{\frac{\partial \widehat{\boldsymbol{U}}}{\partial t}}+\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_S (\boldsymbol{f}\cdot \boldsymbol{n}) \text{d}S=\iiint\limits_{\Delta V}{\widehat{\boldsymbol{S}}}\text{d}V $ |
六面体控制体ΔV的体积为Ω, 控制体中心取(I, J, K), 左面中心取(i, J, K), 右面中心取(i+1, J, K), 以此类推.积分方程可离散如下
| $ \begin{array}{*{35}{l}} {{\left[ \frac{\mathit{\Omega }}{\partial t}\partial \boldsymbol{\hat{U}} \right]}_{i, J, K}}+{{\left( \boldsymbol{\hat{F}}_{i+1}^{n+1}-\boldsymbol{\hat{F}}_{i}^{n+1} \right)}_{J, K}}+{{\left( \boldsymbol{\hat{G}}_{j+1}^{n+1}-\boldsymbol{\hat{G}}_{j}^{n+1} \right)}_{1, K}}+ \\ {{\left( \boldsymbol{\hat{H}}_{k+1}^{n+1}-\boldsymbol{\hat{H}}_{k}^{n+1} \right)}_{1, J}}={{\left( \boldsymbol{\hat{F}}_{\nu , i+1}^{n+1}-\boldsymbol{\hat{F}}_{\nu , i}^{n+1} \right)}_{J, K}}+\left( \boldsymbol{\hat{G}}_{\nu , j+1}^{n+1}- \right. \\ {{\left. \boldsymbol{\hat{G}}_{\nu , j}^{n+1} \right)}_{1, K}}+{{\left( \boldsymbol{\hat{H}}_{\nu , k+1}^{n+1}-\boldsymbol{\hat{H}}_{\nu , k}^{n+1} \right)}_{1, J}}+{{[\boldsymbol{\hat{S}}\mathit{\Omega }]}_{1, J, K}} \\ \end{array} $ |
此处称
| $ \begin{array}{c} \hat{u}=u n_{x}+v n_{y}+w n_{z} \\ \hat{\boldsymbol{F}}_{i+1, J, K}^{n+1}=\sigma_{i+1, J, K}\left(\boldsymbol{F}_{i+1, J, K}^{n+1} n_{x}+\boldsymbol{G}_{i+1, J, K}^{n+1} n_{y}+\boldsymbol{H}_{i+1, J, K}^{n+1} n_{z}\right) \end{array} $ |
其中, σi+1, J, K为所在面的面积, (nx, ny, nz)为该面法向量.
采用Steger-Warming格式[12-13]耦合LU-SGS[14]的全隐式数值模拟方法, 其中, 化学反应源项采用对角化隐式处理[15].
计算中采用以下边界条件:壁面为绝热壁, 且为非催化壁面.
2 结果及分析由于计算网格对超声速爆轰数值计算影响很大, 首先研究网格对弹头直径15 mm在Ma=6. 46时爆轰状态的影响, 并将计算结果与实验结果进行了对比, 验证了所用计算方法的适用性.接着研究了弹头直径对Ma=4.18, 5.11爆轰状态的影响.计算状态见表 1.
|
下载CSV 表 1 计算所需的自由流条件和弹头直径[2] Tab.1 Free stream conditions and projectile diameters for computation[2] |
选择Lehr实验中弹头直径15 mm在Ma=6.46的状态对计算方法进行验证.针对Ma=6.46计算所采用的网格如图 2所示.网格采用商业软件Gridgen生成, 图 2(a)中没有对激波捕捉进行网格加密处理, 网格为321×201×121. 图 2(b)对激波捕捉进行了网格加密处理, 网格为321×291×121, 除了法向加密外, 在激波附近对网格进行了适应性处理, 使其能够更好地对激波、燃烧波等进行捕捉.
|
| 图 2 计算网格 Fig.2 Meshes used for computation |
图 3为计算网格对Ma=6.46密度等值线的影响, 图 4为计算网格对Ma=6.46组分H2O等值线的影响.从图中可以很明显地看出, 当计算网格没有进行激波适应时, 计算结果不能很好地捕捉氢气-空气弹头爆轰流场的结构特征, 也无法捕捉流场中的诱导区域.对网格进行修正后, 可以从图中清楚地分辨出激波以及燃烧阵面, 并且激波和燃烧阵面的位置与实验测得的位置符合良好, 如图 3(b)和图 4(b)所示.在弹头头部滞止线附近, 激波和燃烧阵面重合, 沿头部直到60°左右两者开始被诱导区域分离开, 这一点与文献[12]一致.对于Ma=6.46, 由于滞止区域附近波后温度非常高, 如图 5所示, 导致其诱导距离非常小, 离开滞止线后, 诱导距离随着激波强度以及波后温度下降而逐渐增大.而且, 从图中还可以看出, 由于该状态Ma=6.46高于15 mm弹头氢气-空气混合流动的C-J Ma=5.11, 燃烧流场为稳定状态, 流场中激波和燃烧阵面都非常光滑.
|
| 图 3 网格对密度等值线的影响(z=0截面) Fig.3 Mesh effect on the density contours on z=0 section |
|
| 图 4 网格对组分H2O等值线的影响(z=0截面) Fig.4 effect on the H2O contours on z=0 section |
|
| 图 5 网格对温度等值线的影响(z=0截面) Fig.5 Mesh effect on the temperature contours on z=0 section |
图 6为Ma=6.46时沿驻点流线的压力、温度以及组分分布.从图中可以看出, 由于该状态时化学能量和动能的比值很小, 能量释放没有能够引起很强的爆轰, 热释放对流场的影响不是很大, Von Norman峰值(压力曲线上的极值点)不是很明显.沿驻点流线的温度分布曲线上可以看到一个不很明显的诱导区, 在诱导区域, 生成物浓度急剧增大.
|
| 图 6 密网格时滞止线上压力、温度及组分分布(z=0截面) Fig.6 Pressure, temperature and species along the stagnation line on z=0 section with fine meshes |
图 7为密网格时Ma=6.46状态z=0截面流线放大图.从图 7(a)可以看出, 由于氢气-空气混合物的化学反应时间尺度远大于当地的流动时间尺度, 混合物瞬间反应, 这一点不同于图 7(b).从图 7(b)中可以看出, 随着激波和燃烧阵面逐渐变得倾斜时, 燃烧阵面逐渐与流动方向趋于一致.
|
| 图 7 密网格时z=0截面流线放大图 Fig.7 Enlarged streamline on z=0 section with fine meshes |
图 8为弹头直径5 mm时, Mach数对爆轰流场结构的影响.从图中可以看出, Ma=4.18时, 驻点附近温度达到自燃温度, 氢气燃烧并向下游传播, 但由于Mach数较低, 燃烧区域局限于弹头表面附近; 随着Mach数的增大, 驻点附近的燃烧向上游传播, 并在激波后区域驻定.从图中还可以看出, 激波和燃烧阵面都非常光滑, 流场为稳定状态.
|
| 图 8 D=5 mm弹头在不同Mach数时流场特征 Fig.8 Mach effect on the flowfield characteristics with D=5 mm |
图 9为不同弹头直径时氢气-空气在Ma=5.11时的爆轰流场.从图中可以看出, 弹头直径15 mm时, Ma=5.11时的氢气-空气爆轰流场处于不稳定状态, 如图 9(c)流场组分H2O及流线图所示.随着弹头直径的减小, 流场逐渐变得稳定, 如图 9(b)和图 9(a)所示. D=5 mm爆轰流场是稳定的, 其中的激波和燃烧波阵面非常光滑, 如图 9(a)所示. 图 9(c)与图 9(b)的流动结构类似, 但前者的流场不稳定占优.因此, 弹头Mach数固定不变时, 弹头直径越小, 流动稳定性越强.
|
| 图 9 Ma=5.11时弹头直径对氢气-空气爆轰流场稳定性的影响 Fig.9 Projectile diameter effect on the flowfield stability at Ma=5.11 |
本文采用数值模拟方法对氢气-空气燃烧流场进行了研究, 重点分析了弹头飞行Mach数、弹头直径对氢气-空气燃烧爆轰流场稳定性的影响.研究结果表明:弹头直径一定时, 氢气-空气燃烧流场稳定性随着飞行Mach数的增大而增强; 弹头飞行Mach数一定时, 氢气-空气燃烧流场稳定性随着弹头直径减小而增强.
| [1] |
Wilson G J, MacCormack R W. Modeling supersonic combustion using a fully implicit numerical method[J]. AIAA Journal, 1992, 30(4): 1008-1015. |
| [2] |
Lehr H F. Experiments on shock-induced combustion[J]. Astronautica Acta, 1972, 17(4): 589-597. |
| [3] |
Moon S Y, Lee C W, Sohn C H. Adaptive finite ele-ment analysis of shock-induced combustion[R]. AIAA 2002-3857, 2002.
|
| [4] |
Yungster S, Eberhardt S, Bruckner A P. Numerical simulation of hypervelocity projectiles in detonable gases[J]. AIAA Journal, 1991, 29(2): 187-199. |
| [5] |
Wilson G J, MacCormack R W. Modeling supersonic combustion using a fully implicit numerical method[J]. AIAA Journal, 1992, 30(4): 1008-1015. |
| [6] |
Matsuo A, Fujiwara T. Numerical simulation of shock-induced combustion around an axisymmetric blunt body[R]. AIAA 1991-1414, 1991.
|
| [7] |
李俊红, 潘宏禄, 程晓丽, 等. 标模外形化学非平衡流场数值模拟[J]. 计算物理, 2015, 32(4): 395-402. Li J H, Pan H L, Cheng X L, et al. Numerical simulation on chemical nonequilibrium flowfield in standard model[J]. Chinese Journal of Computational Physics, 2015, 32(4): 395-402. DOI:10.3969/j.issn.1001-246X.2015.04.003 (in Chinese) |
| [8] |
李俊红, 潘宏禄, 沈清, 等. 超燃冲压发动机燃烧室的燃烧特性[J]. 航空动力学报, 2014, 29(1): 14-22. Li J H, Pan H L, Shen Q, et al. Combustion characteristics of scramjet combustor[J]. Journal of Aerospace Power, 2014, 29(1): 14-22. (in Chinese) |
| [9] |
Gnoffo P A, Gupta R N, Shinn J L. Conservation equations and physical models for hypersonic air flows in thermal and chemical nonequilibrium[R]. NASA Technical Paper 2867, 1989.
|
| [10] |
Anderson D A, Tannehil J C, Pletcher R H. Computational fluid mechanics and heat transfer[M]. New York: McGraw-Hill Book Company, 1984.
|
| [11] |
Clutter J K, Mikolaitis D W, Shyy W. Effect of reaction mechanism in shock-induced combustion simulations[R]. AIAA 1998-0274, 1998.
|
| [12] |
Atwood C A. Computation of viscous blast wave flow-fields[R]. NASA-CR-18808091N22512, 1991.
|
| [13] |
Campbell C H, Candler G V. A flux consistent implementation of flux vector splitting[R]. AIAA 2004-243, 2004.
|
| [14] |
Parsani M, Van den Abeele K, Lacor C, et al. Simula-tion of compressible turbulent flows with an implicit LU-SGS algorithm for high-order spectral difference method on unstructured grids[R]. AIAA 2009-4143, 2009.
|
| [15] |
Ahuja J K, Tiwari S N. A parametric study of shock-induced combustion in a hydrogen-air system[R]. AIAA 1994-674, 1994.
|