Application of Three-Dimensional Meshless Method in Muzzle Flow Field of Projectile with Large Displacement

Application of Three-Dimensional Meshless Method in Muzzle Flow Field of Projectile with Large Displacement

Liang WangRui Xue Ning Cai Wei Wu Dongliang Zhang

School of Energy and Power Engineering, Nanjing Institute of Technology, Nanjing 211167, China

Beijing Aerospace Technology Research Institute, Beijing 100074, China

Corresponding Author Email: 
wl_njit@njit.edu.cn
Page: 
41-50
|
DOI: 
https://doi.org/10.18280/ijht.390105
Received: 
15 September 2020
|
Revised: 
24 November 2020
|
Accepted: 
3 December 2020
|
Available online: 
28 February 2021
| Citation

© 2021 IIETA. This article is published by IIETA and is licensed under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).

OPEN ACCESS

Abstract: 

The mesh-less method abandons the idea of traditional mesh method. In numerical calculation, it only needs the node information, and it’s not necessary to connect nodes into mesh cells. Aiming at the problem of non-equilibrium chemical reaction flows in a three-dimensional muzzle with moving boundaries, this paper proposed a three-dimensional least squares mesh-less (3D-LSM) algorithm, and the ALE (Arbitrary Lagrangian-Eulerian) form multi-component Euler governing equations were adopted for fluid dynamics modeling; then, a finite-rate reaction model was used to deal with the chemical reactions in the flow field, a time-division method was applied to solve the stiff problem of chemical reactions, and multi-component AUFS (artificially upstream flux vector splitting) format was introduced to calculate the convective flux term of the governing equation. In addition, the weighted-point filling strategy and the dynamic point cloud method were employed to deal with the topological structure changes in the point clouds caused by the large displacement movement of the projectile in the three-dimensional muzzle flow field. Finally, the proposed mesh-less method was applied to calculate the muzzle of a 12.7mm gun, and accurately captured the bottle-shaped shock wave structure in the complex muzzle flow field and the second muzzle flame phenomenon; the calculation results were compared with the experimental results and obtained a good match, which had proved the feasibility of the proposed 3D-LSM method. The research findings of this paper provided a new solution for the simulation of high-speed non-equilibrium reaction jets with multi-component and large-displacement moving boundaries.

Keywords: 

meshless method, dynamic cloud of point, non-equilibrium reaction, AUFS, three-dimensional muzzle flow field, large-Scale movable boundary

1. Introduction

In three dimensions, the problems of non-equilibrium reaction flow with moving boundaries [1-5] widely exist in engineering practice, such as the non-equilibrium flow formed by the high-temperature environment outside hypersonic aircraft, the interference flow of the internal and external flow of the aircraft , the atmospheric layer re-entry process of the aircraft, and the gunpowder gas flow in the muzzle of the guns, etc. The problem of non-equilibrium chemical reaction flows [3, 6-11] involves multiple disciplines such as aerodynamics, air thermodynamics, aero-physics, and chemical reaction kinetics, etc. Then the numerical simulation of such problem has always been a hot spot in the field of CFD (Computational Fluid Dynamics) research.

The challenges in the numerical simulation of complex chemical reaction flows with moving boundaries are mainly reflected in [12, 13] the following aspects: the complex physical process, the stiff problem of chemical reaction calculation [14] , and the huge amount of calculation, etc. The flow field contains strong discontinuity surfaces such as moving boundaries, shock wave, and combustion wave, etc., the deformed meshes near the moving boundaries need to be reconstructed, and the judgment of intersection and penetration among points, lines and surfaces in the reconstruction process is extremely complicated. Traditional CFD methods are mostly based on meshes, but the existence of meshes has limited the application of these methods in solving flow fields with complex boundaries [15-20]. Therefore, the method that can perform numerical analysis on flow problems without the necessity of meshing has become a research hotspot and the mesh-less method [21-27] has been applied to the study of fluid mechanics ever since. After that, with the rapid development of computer hardware technology in the 20th century and the continuous emergence of high-efficiency algorithms, the research on the numerical simulation of chemical reaction flows based on mesh-less method has also made some progress, but most of the existing studies are 2D calculations or 3D non-dynamic calculations [7, 12, 28-30]. In the 1970s, Liszka and Orkisz [31], Liszka and Orkisz [32] established a prototype of mesh-less algorithm. Salam [26] calculated the movement of floating objects very well under the novel meshless local Petrov-Galerkin. Based on the finite point method, Ortega et al. [33] and Based on the Taylor–Galerkin format, Erhart et al. [34] solved the compressible flows in 2D and 3D space. Firoozjaee and Afshar [35] proposed a localized and virtualized radial basis function, and gave a local mesh-less method based on a Cartesian point-distribution scheme. The problem of supersonic flow around a plate was simulated well. Munikrishna and Balakrishnan [36] proposed the DLSM method to simulate the steady-state incompressible N-S equation and the time step could be large. Gargari et al. [37] adopted a hybrid mesh-less Cartesian point arrangement method, and calculated the turbulent flow using the RANS turbulence model of the mesh-less solver LSFD-U. In response to the incompressible Navier-Stokes equation and ALE equation, Hong et al. [38] put forward a mixed discrete least square mesh-less (MDLSM) method. Compared with the DLSM method, its accuracy is higher.

For conditions with small rigid boundary displacement and little changes in existing point cloud structure, commonly-used methods include algebraic analogy model, Delaunay graph mapping method [39] , and spring analogy method [40], etc. As for movements with larger rigid boundary displacement, the commonly adopted methods include the point cloud overlapping method [41, 42] , and the local point cloud reconstruction method [13], etc. The local point cloud reconstruction method is suitable for dealing with the unsteady flow problems with complicated boundary deformation; its idea is quite simple: delete the point cloud whose quality deteriorated due to the boundary deformation from the original point cloud topological structure and form a cavity, then use the automatic point cloud generation method to re-generate a high-quality point cloud in the cavity, and add it to the original point cloud system to continue the calculation.

To solve the stiff problem of chemical reaction calculation, Trotter [43] proposed a high-resolution detonation wave calculation program based on the decoupling algorithm combined with the structural grid adaptive method. Aiming at the problem of sufficiently small time steps, Strang [44] analyzed the classic decoupling algorithm Lie splitting method [45] and Strang splitting method [46], and the time accuracy were order 1 and order 2, respectively. Guzzo and Azevedo [47] proposed a novel mesh-based MTS method and applied the elementary reaction kinetics to improve computational efficiency. In the time-division method, the chemical reaction is a set of rigid nonlinear ordinary differential equations. At present, most of the calculations are conducted using the VODE solver [48, 49] in the Chemkin software package.

The mesh-less method has received wide attention due to its inherent superiority, flexibility, and variability. However, when dealing with complex flow fields with moving boundaries, its performance is far less mature than the unstructured mesh method; especially in terms of 3D problems, its flux calculation of points cloud is more complicated. However, few studies have been conducted on it. This paper started from the 3D dynamic point cloud processing method and established a 3D mesh-less algorithm for the numerical simulation of point cloud topological structure changes and the non-equilibrium chemical reaction flow fields, then it calculated two practical problems, and the calculation results were in good match with the experimental results, which had verified the effectiveness of the algorithm.

2. Governing Equation and Numerical Method

The ALE (Arbitrary Lagrangian-Eulerian) method [39] allows the meshes to move at any speed, and is suitable for solving flow fields containing moving boundaries with any displacement value. The mesh-less simulation algorithm for chemical reaction flow field developed in this paper is based on the ALE method, and its ALE-form Euler equation containing the chemical reaction source term can be expressed in a rectangular coordinate system as:

${{\mathbf{U}}_{t}}+{{\mathbf{F}}_{1,x}}+{{\mathbf{F}}_{2,y}}+{{\mathbf{F}}_{3,z}}=\mathbf{S}$     (1)

The specific expressions are as follows:where, the subscripts ", x", ", y", ", z" mean to calculate the partial derivatives of x, y, and z; U is a conserved variable; F1, F2, F3 are convective flux terms; S is the chemical reaction source term.

$\boldsymbol{U}=\left(\begin{array}{llllllll}\rho_{1} & \rho_{2} & \cdots & \rho_{N_{R}} & \rho u & \rho v & \rho w & \rho E\end{array}\right)^{T}$

$\boldsymbol{S}=\left(\begin{array}{llllllll}\omega_{1} & \omega_{2} & \cdots & \omega_{N_{R}} & 0 & 0 & 0 & 0\end{array}\right)^{T}$

$\boldsymbol{F}_{1}=\left(\begin{array}{c}\rho_{1}(u-\dot{u}) \\ \rho_{2}(u-\dot{u}) \\ \vdots \\ \rho_{N_{R}}(u-\dot{u}) \\ \rho u(u-\dot{u})+p \\ \rho v(u-\dot{u}) \\ \rho w(u-\dot{u}) \\ \rho E(u-\dot{u})+u p\end{array}\right)$$\boldsymbol{F}_{2}=\left(\begin{array}{c}\rho_{1}(v-\dot{v}) \\ \rho_{2}(v-\dot{v}) \\ \vdots \\ \rho_{N_{R}}(v-\dot{v}) \\ \rho u(v-\dot{v})+p \\ \rho v(v-\dot{v}) \\ \rho w(v-\dot{v}) \\ \rho E(v-\dot{v})+v p\end{array}\right)$$\boldsymbol{F}_{3}=\left(\begin{array}{c}\rho_{1}(w-\dot{w}) \\ \rho_{2}(w-\dot{w}) \\ \vdots \\ \rho_{N_{R}}(w-\dot{w}) \\ \rho u(w-\dot{w})+p \\ \rho v(w-\dot{w}) \\ \rho w(w-\dot{w}) \\ \rho E(w-\dot{w})+w p\end{array}\right)$

where, ρ is the mass density of the mixed gas, ρi is the mass density of component i; NR is the total number of components; u, v, and w are the velocity components in the x, y, and z directions of the mixed gas; $\dot{u}$, $\dot{v}$ and $\dot{w}$ are the velocity components of the node in the x, y, and z directions; p is the pressure of the mixed gas; ωi is the mass production rate of component i caused by the chemical reaction; ρE is the total energy per unit volume:

$\omega_{i}=\sum_{t=1}^{R_{N}}\left(v_{i t}^{R}-v_{i t}^{L}\right)$$\left[K_{f t} \prod_{k=1}^{N_{R}}\left(\frac{\rho_{k}}{M_{k}}\right)_{k^{L}}^{v_{k}}-K_{b t} \prod_{k=1}^{N_{R}}\left(\frac{\rho_{k}}{M_{k}}\right)^{v_{k^{R}}^{R}}\right] \cdot M_{i}$    (2)

where, $\boldsymbol{v}_{i t}^{L}$ and $\boldsymbol{v}_{i t}^{R}$ are the chemical reaction equivalent coefficients of component i on the left and right sides of reaction t; Mi is the molar mass of component i; RN is the total number of chemical reactions in the model; Kft is the forward and reverse reaction rate constant, its value is obtained by calculating the Arrhenius formula:

$K_{f t}=A_{t} T^{b_{t}} \exp \left(-\frac{E_{t}}{R_{u} T}\right)$

where, Atis the pre-exponential factor; bt is the temperature factor; Et is the activation energy; Ru is the universal gas constant; T is the temperature of the mixed gas.

The equilibrium constant [13] is used to calculate the reverse reaction rate constant, and the chemical reaction source term is handled with the finite-rate reaction model.

2.1 Spatial dispersion

In case of 3D space, we assum that the following linear relationship is introduced to calculate fluid flow variables:

$\begin{align}  & f={{a}_{1}}x+{{a}_{2}}y+{{a}_{3}}z+{{a}_{0}} \\ & =\frac{\partial f}{\partial x}x+\frac{\partial f}{\partial y}y+\frac{\partial f}{\partial z}z+{{a}_{0}} \\\end{align}$     (3)

In every point cloud structure, the flow variables of the center point i satisfy the above relationship, and the satellite points j (j=1, 2, ..., N) are also. Then:

$\left\{\begin{array}{l}f_{i}=a_{1} x_{i}+a_{2} y_{i}+a_{0} \\ f_{j_{1}}=a_{1} x_{j_{1}}+a_{2} y_{j_{1}}+a_{0} \\ \cdots \\ f_{j_{6}}=a_{1} x_{j_{6}}+a_{2} y_{j_{6}}+a_{0}\end{array}\right.$

Add the relational expression of center point i to the relational expression corresponding to each satellite point:

$\left\{\begin{array}{l}\left(f_{j_{1}}+f_{i}\right) / 2=a_{1}\left(x_{j_{1}}+x_{i}\right) / 2+a_{2}\left(y_{j_{1}}+y_{i}\right) / 2+a_{0} \\ \left(f_{j_{2}}+f_{i}\right) / 2=a_{1}\left(x_{j_{2}}+x_{i}\right) / 2+a_{2}\left(y_{j_{2}}+y_{i}\right) / 2+a_{0} \\ \cdots \\ \left(f_{j_{6}}+f_{i}\right) / 2=a_{1}\left(x_{j_{6}}+x_{i}\right) / 2+a_{2}\left(y_{j_{6}}+y_{i}\right) / 2+a_{0}\end{array}\right.$

Further, it can be written as a matrix:

$\begin{align}  & \left( \begin{matrix}   \frac{{{x}_{1}}+{{x}_{i}}}{2} & \frac{{{y}_{1}}+{{y}_{i}}}{2} & \frac{{{z}_{1}}+{{z}_{i}}}{2} & 1  \\   \vdots  & \vdots  & \vdots  & \vdots   \\   \frac{{{x}_{N}}+{{x}_{i}}}{2} & \frac{{{y}_{N}}+{{y}_{i}}}{2} & \frac{{{z}_{N}}+{{z}_{i}}}{2} & 1  \\\end{matrix} \right)\cdot \left( \begin{matrix}   {{a}_{1}}  \\   {{a}_{2}}  \\   {{a}_{3}}  \\   {{a}_{0}}  \\\end{matrix} \right) \\ & =\left( \begin{matrix}   \frac{{{f}_{1}}+{{f}_{i}}}{2}  \\   \vdots   \\   \frac{{{f}_{N}}+{{f}_{i}}}{2}  \\\end{matrix} \right) \\\end{align}$      (4)

Observing the above equation, the coefficient matrix is written as B. In the flow field dispersion, we find the number of satellite points of a node is usually greater than 4. Thus the above equations are contradictory equations. At any node i, the space derivative was calculated by the least square method:

$\left.\frac{\partial f}{\partial x}\right|_{i}=a_{1}=\sum_{j=1}^{N} b_{i j} \frac{f_{i}+f_{j}}{2}$, $\left.\frac{\partial f}{\partial y}\right|_{i}=a_{2}=\sum_{j=1}^{N} c_{i j} \frac{f_{i}+f_{j}}{2}$, $\left.\frac{\partial f}{\partial z}\right|_{i}=a_{3}=\sum_{j=1}^{N} d_{i j} \frac{f_{i}+f_{j}}{2}$

In the formula, bij is the first row of matrix $\left(\boldsymbol{B}^{T} \boldsymbol{B}\right)^{-1} \boldsymbol{B}^{T}$, cij is the second row, and dij is the third row. Then, for the node i, the convective flux term of the ALE equation can be expressed as:

${{\left. \frac{\partial {{\mathbf{F}}_{1}}}{\partial x} \right|}_{i}}+{{\left. \frac{\partial {{\mathbf{F}}_{2}}}{\partial y} \right|}_{i}}+{{\left. \frac{\partial {{\mathbf{F}}_{3}}}{\partial z} \right|}_{i}}=\sum\limits_{j=1}^{N}{\left( \sqrt{b_{ij}^{2}+c_{ij}^{2}+d_{ij}^{2}}{{\mathbf{W}}_{ij}} \right)}$     (5)

2.2 The 3D multi-component AUFS format

The artificially upstream flux vector splitting (AUFS for short) format was developed by Sun et al. [44]. This format is easy to calculate and can accurately distinguish the contact discontinuity surfaces, it doesn’t show carbuncle phenomenon under multiple dimensions, and has good robustness. For the 3D multi-component problems, the specific form of the AUFS applied to the mesh-less method is as follows:

$W_{ij}^{AUFS}=(1-M){{W}_{1}}+M{{W}_{2}}$    (6)

The fluid vectors W1 and W2, and the constant M in the formula are defined as follows:

${{\mathbf{W}}_{1}}=\left[ \frac{1}{2}\left( {{\mathbf{P}}_{i}}+{{\mathbf{P}}_{j}} \right)+\delta \mathbf{U} \right]$,${{\mathbf{W}}_{2}}=\left[ {{\mathbf{U}}^{\alpha }}\left( {{{\tilde{u}}}^{\alpha }}-{{{\tilde{\dot{u}}}}^{\alpha }}-{{S}_{2}} \right)+{{\mathbf{P}}^{\alpha }} \right.$,$M={{S}_{1}}/\left( {{S}_{1}}-{{S}_{2}} \right)$     (7)

In the formula, P and δU are defined as follows:

$\tilde{u}^{\alpha}=n_{x} u_{\alpha}+n_{y} v_{\alpha}+n_{z} w_{\alpha}$,

$\tilde{\dot{u}}^{\alpha}=n_{x} \dot{u}_{\alpha}+n_{y} \dot{v}_{\alpha}+n_{z} \dot{w}_{\alpha}$,

$\boldsymbol{P}=\left(0,0, \cdots, 0, p l_{x}, p l_{y}, p l_{z}, p \tilde{u}\right)^{T}$,

$l_{x}=b_{i j} / e_{i j} l_{y}=c_{i j} / e_{i j}$,

$l_{z}=c_{i j} / e_{i j}, e_{i j}=\sqrt{b_{i j}^{2}+c_{i j}^{2}+d_{i j}^{2}}$,

$\delta \boldsymbol{U}=\frac{1}{2 c_{i}}\left(\left(Y_{1} p\right)_{i} \quad\left(Y_{2} p\right)_{i} \quad \cdots\left(Y_{N_{R}} p\right)_{i} \quad(p u)_{i} \quad(p v)_{i} \quad(p w)_{i} \quad(H p)_{i}\right)^{T}-$$\frac{1}{2 c_{j}}\left(\left(Y_{1} p\right)_{j} \quad\left(Y_{2} p\right)_{j} \quad \cdots\left(Y_{N_{k}} p\right)_{j}(p u)_{j}(p v)_{j}(p w)_{j}(H p)_{j}\right)^{T}$

where, Yi=ρi/ρ; the values of a, coefficient M, and scaling functions S1 and S2 are obtained from the formulas below:

$\alpha=\left\{\begin{array}{ll}i & S_{1}>0 \\ j & S_{1} \leq 0\end{array}\right.$,

$M=S_{1} /\left(S_{1}-S_{2}\right), S_{1}=\left(u_{i}+u_{j}\right) / 2$

$S_{2}=\left\{\begin{array}{l}\operatorname{Min}\left(0, u_{i}-\dot{u}-c_{i}, u^{*}-c^{*}\right) \quad S_{1}>0 \\ \operatorname{Max}\left(0, u_{j}-\dot{u}+c_{j}, u^{*}+c^{*}\right) \quad S_{1} \leq 0\end{array}\right.$

where, u* and c* are the intermediate wave velocity and intermediate sound velocity derived from the isentropic relationship, and the specific expressions are:

$u^{*}=\frac{\left(u_{i}+u_{j}\right)}{2}+\frac{c_{i}-c_{j}}{\gamma-1}$, $c^{*}=\frac{\left(c_{i}+c_{j}\right)}{2}+\frac{(\gamma-1)\left(u_{i}-u_{j}\right)}{4}$

where, γ is the gas specific heat ratio.

The time term performs explicit time marching according to the fourth-order Runge-Kutta method.. The form is as follows:

$\boldsymbol{U}_{i}^{(0)}=\boldsymbol{U}_{i}^{n}$,  $\boldsymbol{U}_{i}^{(j)}=\boldsymbol{U}_{i}^{(0)}-\phi_{j} \Delta t \boldsymbol{R}\left(\boldsymbol{Q}_{i}^{j-1}\right)$$\phi_{j}=\frac{1}{4-j+1} \quad j=1,2,3,4$, $\boldsymbol{U}_{i}^{n+1}=\boldsymbol{U}_{i}^{(4)}$

$\boldsymbol{R}\left(\boldsymbol{Q}_{i}\right)=\sum_{k=1}^{N} \boldsymbol{W}_{i k}+\boldsymbol{S}\left(\boldsymbol{U}_{i}^{(0)}\right), \boldsymbol{Q}_{i}=\left(\boldsymbol{U}_{i}, \boldsymbol{U}_{i 1}, \boldsymbol{U}_{i 2}, \cdots \boldsymbol{U}_{i N}\right)$

where, Δt is the calculation time step, and its value adopts the full field time step, which can be determined by the minimum value of the local time step at each point.

For the stiff problem of chemical reaction calculation, the time-division algorithm [39, 40, 45] of the chemical reaction and flow decoupling calculation are adopted, and the time step of the chemical reaction is further subdivided. The boundaries are processed by the method of constructing the mirror point outside the flow field, and the value of the flow variable of the mirror point is determined according to the boundary type. The solid surface adopted the normal impermeable boundary condition, and the far field adopts the non-reflecting boundary condition.

3. 3D Point Cloud Generation and Dynamic Point Cloud Processing Method

3.1 Point cloud generation

This study introduces the weighed-point filling strategy [25] to the three dimensional space. The idea of the method is to set a real weight for each discrete node, and the size of the weight is represented as the radius of the sphere with the node as the center of the sphere. The filling of weighted points is to fill the calculation area with spheres of same or different radius. In actual operation, to make the method more flexible, two adjacent virtual spheres are not necessarily tangent, partial intersection or separate are also allowed. The specific processes can be divided as follows:

(1) Observing the actual problem, on the discrete boundary surface, points can be distributed based on certain principles. The distribution can be uniform or non-uniform. The weight of each boundary point i is $r_{i}=\frac{1}{2 N} \sum_{j=1}^{N} d_{i j}$, dij is the distance between point i and adjacent node j, and N is the number of adjacent points.

(2) With boundary triangle as the initial front, weighted points were filling inside. As shown in Figure 1, taking the front ABC as an example, the weight value rD of the best marching weighted point D was the average of the weight values rA, rB and rC of the three endpoints A, B, and C; according to the principle that the sphere where node D is located is tangent to the spheres where nodes A, B, and C are located, then the coordinate of weighted point D was (xD, yD, zD). βr is the weight coefficient, because two adjacent spheres were allowed to be partially intersected or separated, the purpose of non-uniform point distribution could be achieved at the same time.

$\left\{ \begin{array}{*{35}{l}}   {{\left( {{x}_{A}}-{{x}_{D}} \right)}^{2}}+{{\left( {{y}_{A}}-{{y}_{D}} \right)}^{2}}+{{\left( {{z}_{A}}-{{z}_{D}} \right)}^{2}}={{\left( {{r}_{A}}+{{r}_{D}} \right)}^{2}}  \\   {{\left( {{x}_{B}}-{{x}_{D}} \right)}^{2}}+{{\left( {{y}_{B}}-{{y}_{D}} \right)}^{2}}+{{\left( {{z}_{B}}-{{z}_{D}} \right)}^{2}}={{\left( {{r}_{B}}+{{r}_{D}} \right)}^{2}}  \\   {{\left( {{x}_{C}}-{{x}_{D}} \right)}^{2}}+{{\left( {{y}_{C}}-{{y}_{D}} \right)}^{2}}+{{\left( {{z}_{C}}-{{z}_{D}} \right)}^{2}}={{\left( {{r}_{C}}+{{r}_{D}} \right)}^{2}}  \\   {{r}_{D}}=\frac{1}{3}\cdot \left( {{r}_{A}}+{{r}_{B}}+{{r}_{C}} \right)\cdot {{\beta }_{r}}  \\\end{array} \right.$     (8)

Figure 1. Mesh-less weighted point filling

(3) Then, look for the existing weighted points around point D (the radius of the area was 1.5 rD), shown as points D1 and D2 in Figure 1, and establish the linked list of candidate points. Then, after geometric judgments such as penetration and intersection, invalid weighted points were deleted. At last, according to the existing weighted points and the overlap coefficient of the ideal marching point, the optimal marching point was determined; the overlap coefficient was defined as follows:

${{\beta }_{ij}}=\left\{ \begin{array}{*{35}{l}}   0 & \left( {{d}_{ij}}\ge {{r}_{i}}+{{r}_{j}} \right)  \\   1-\frac{{{d}_{ij}}}{{{r}_{i}}+{{r}_{j}}} & \left( \left| {{r}_{i}}-{{r}_{j}} \right|\le {{d}_{ij}}\le {{r}_{i}}+{{r}_{j}} \right)  \\   \frac{2\cdot min\left( {{r}_{i}},{{r}_{j}} \right)}{{{r}_{i}}+{{r}_{j}}} & \left( {{d}_{ij}}\le \left| {{r}_{i}}-{{r}_{j}} \right| \right)  \\\end{array} \right.$     (9)

where, ri and rj are the weight values of weighted points i and j, respectively, and dij is the distance between the two points. The overlap coefficients βij of all candidate points were calculated, and the point with the largest βij value was selected as the optimal marching point.

3.2 Dynamic point cloud processing

In a 3D space, when a rigid body in the flow field moves with a large displacement, the surrounding point cloud will be compressed or stretched, resulting in a decrease in the quality of the point cloud, and the calculation accuracy will be affected. This paper adopted the local point cloud reconstruction method to solve the problem of point cloud deformity caused by the large-displacement moving rigid body in the flow field. The local point cloud reconstruction method mainly includes the following steps:

(1) Find the internal nodes of the moving rigid body that caused point cloud deformity and form the reconstruction cavity, then establish the linked list of the initial reconstruction front. The judgment criterions for the nodes to be deleted are: 1) Detect the distance Lij between the center point i and its satellite point j, when Lij<αmin·2r0, or Lij>αmax·2 r0, point i is set as a node to be deleted, wherein r0is the average distance weight of the nodes on the moving boundary, αmin and αmax are respectively the tensile limiting coefficient and the compressive limiting coefficient; 2) Find the satellite points of center point i, and connect any three points that are mutually satellite points (S-S-S points) to form non-overlapping triangles, then check whether the spherical degree θ of the trihedral angle formed by each triangle and the center point i exceeds the angle limit, that is, whenθ<βmin·θ0, orθ>βmax·θ0, point i is set as a node to be deleted; the average angle is θ0=4π/N; N is the total number of triangles; βmin and βmax are angle limiting coefficients.

(2) Adopt the weighted-point filling strategy to re-arrange the points in the cavity and update the point cloud

With the cavity boundary front triangle as the initial front, the weighted-point filling strategy was adopted to generate points and point cloud. The quality of the cloud was improved using the following method: 1) SS cut or connect, as shown in Figure 2, point A and point B are satellite points for each other, and points C1, C2, and C3 are common satellites of point A and point B, when the intersection of line AB and triangle C1C2C3 is inside triangle C1C2C3, the satellite relationship between point A and point B can be deleted to improve the quality of the point cloud without affecting other point clouds. In a more common situation, when the number of common satellite points is greater than 3, the polygon formed by the common satellite points can be divided into multiple triangles, and then similar operations can be performed; also, the above processing can be reversed; 2) Laplace smoothing, taking any node i (xi, yi, zi) in the flow field as an example, the smoothing method can be expressed as:

$\left\{ \begin{array}{*{35}{l}}   x_{i}^{n+1}=\lambda x_{i}^{n}+\frac{(1-\lambda )}{N}\sum\limits_{k=1}^{N}{x_{k}^{n}}  \\   y_{i}^{n+1}=\lambda y_{i}^{n}+\frac{(1-\lambda )}{N}\sum\limits_{k=1}^{N}{y_{k}^{n}}  \\   z_{i}^{n+1}=\lambda z_{i}^{n}+\frac{(1-\lambda )}{N}\sum\limits_{k=1}^{N}{z_{k}^{n}}  \\\end{array} \right.$     (10)

where, n represents the number of smoothing operations, usually its value is between 3-5; λ is the smoothing factor, its value is 0<λ<1, the subscript k represents the satellite point of the center point i; N is the total number of satellite points.

Figure 2. S-S cut or connect

(3) The linear interpolation method was used to transfer the flow parameters of the new and old nodes in the flow field, delete the old nodes and point clouds from the data structure, and add the new nodes and point clouds to the data structure.

Figure 3 is an example of the local reconstruction method. The sphere moves in the flow field and its displacement is relatively large. From the point cloud after the movement, it can be seen that the reconstruction effect was relatively good. The following formula is adopted to calculate the quality of the point cloud; the minimum quality factor is greater than 0.3, and the percentage of point clouds with a minimum quality factor greater than 0.8 exceeds 65%, which has satisfied the calculation requirements.

$Q=\sqrt{\frac{R_{0}-\Delta R_{\max }}{R_{0}} \cdot \frac{\theta_{0}-\Delta \theta_{\max }}{\theta_{0}}}$

where, R0 is the maximum value of the distances between the center point and all satellite points.

Figure 3. Example of local point cloud reconstruction method

4. Numerical Examples

4.1 High-speed projectile

The projectile flies at v0=1837m/s in a premixed gas. The mixed gas consists of methane and air with an equal stoichiometric ratio. Figure 4(left) shows the shape of the projectile. Its head cone angle is 60°. The angle of attack is $\alpha$-4°. The temperature and pressure of the surrounding air are T0=298K and p0=101325Pa. The flow field is a cylinder with 809,580 nodes (Figure 4(right)). Its radius is 0.07m and Its length is 0.16m.

Figure 4. The structure of the projectile(left) and flow field node distribution (partial, right)

The calculated result is showed in Figure 5 (right). Compared with the experimental result (Figure 5 (left), the present method can correctly capture the shock waves on the upper and lower side of the projectile, which are asymmetric under the angle of attack. Figure 6 gives the flow field temperature cloud at the z=0 section and the pressure cloud at the y=0 section. The gas near the top of the projectile is compressed, causing the temperature and pressure to rise, where the temperature reaches up to 1200K, which is not enough to cause combustion. Behind the tip of the projectile, the presence of the projectile wall makes the temperature and pressure drop rapidly, then the pressure starts to rise, and the temperature slowly drops.

Figure 5. Experimental shadow photos (left) and density cloud images of XOY section (right)

Figure 6. Temperature cloud at the z=0 section and pressure cloud at the y=0 section

4.2 Calculation of muzzle flow field with high-speed moving projectile

The muzzle flow field is usually a complicated highly under-expanded combustion gas jet. This study adopted the 3D-LSM method to calculate the muzzle flow field with large-displacement moving boundaries. The calculation area of the muzzle is an axisymmetrical structure, as shown in Figure 7. The left figure shows the calculation area on the XOY section of the D=12.7mm anti-aircraft machine gun, and the right figure shows the node distribution diagram. The length of the rifled tube was 1.08m; the flow field outside the muzzle was a cylinder with a radius of 0.45m and a length of 1.1m. Since the flow field was symmetrical, 1/4 of the cylinder was taken for calculation. According to the characteristics of the multi-component muzzle flow field, the area near the flow field axis was encrypted in the process of flow field node arrangement, as shown in the right figure of Figure 7, a total of 2735806 points had been arranged.

Figure 7. Muzzle calculation area (left) and flow field node distribution (partial, right)

Figure 8 is a comparison between the calculated density shadow map and the experimental shadow map at t=3.55ms. Table 1 gives the characteristic parameter values of the shockwave bottle in the flow field at this moment. The two shapes agreed well with a small error, which had verified the feasibility of applying mesh-less algorithm in the flow field of muzzle chemical reactions.

Figure 8. Comparison of XOY section of the calculated shadow map and the experimental shadow map

Table 1. Calculated results and experimental results of characteristic parameters of muzzle flow field and the errors

 

L1/m

L2/m

L3 /m

Calculated result

0.1593

0.2878

0.1091

Experimental result

0.1683

0.2794

0.1035

Error

5.35%

3.01%

5.41%

Figure 9 gives the muzzle flow field temperature cloud at the y=0 section and the OH distribution cloud at the z=0 section at t=3.25ms, t=3.42ms, t=3.55ms, t=3.73ms, t=3.92ms,and t=4.10ms, showing the muzzle ejection process of highly under-expanded gunpowder gas in the barrel and the formed complex wave structure in the muzzle. OH is an intermediate product of the chemical reaction of H2 and O2 in the gunpowder gas. As the projectile flew away from the muzzle, the gunpowder gas expanded very rapidly, and the ejection speed changed suddenly and was greater than the speed of the projectile. The high-temperature and high-pressure gunpowder gas behind the projectile bypassed the projectile, and moved in the axial and radial directions of the muzzle, forming a spherical gunpowder gas shock wave around the projectile, and there’s bottom shockwave at the bottom of the projectile, as shown in Figure 9. At the same time, the high-temperature gunpowder gas contained incompletely combusted combustible gas H2 and CO. At the edge of the flow field, since the air had been drawn in, the gunpowder gas was ignited and undergone chemical reactions, namely the second muzzle flame, which caused the temperature of the flow field in this area to rise obviously, and the mass fraction of OH increased obviously. The gas in the core area of the flow field has the characteristics of high speed and low density. When the gunpowder gas passed through this area, it was accelerated, and its propagation speed in the axial direction was greater than that in the radial direction, forming a high-pressure air mass protruding forward, which then gradually evolved into a coronal muzzle shock wave. In Figure 9, from left to right, shows the initial shock wave, muzzle shock wave, bottom shock wave, and gunpowder gas shock wave. As the flow field continued to develop, the projectile gradually detached from the gunpowder gas jet, the bottom shock wave disappeared, the Mach disk of the gunpowder gas was no longer affected by it and began to close. At this time, the bottle-shaped shock wave was formed by the Mach plate and the intersecting shock wave, as shown in Figure 9 (c). At the same time, the jet boundary curled outwards, vortex rings appeared at the edges, and the number of vortex rings increased continuously, which accelerated the mixing of gunpowder gas and external air, and some of the external air had been drawn into the flow field, causing combustion reaction in the flow field, and the temperature rose accordingly; the OH mass fraction increased, but it’s mainly distributed in local areas where vortex rings existed. After that, the speed of the gunpowder gas further decreased and was lower than the speed of the projectile; the projectile passed through the muzzle shock wave and a shock wave was formed on its head, as shown in Figure 9 (c-f).

Figure 10 is the temperature cloud of the YOZ section at 9D on the right side of the muzzle. The temperature changes at three moments showed that, the combustion was initiated from the boundary area of the flow field. As the jet boundary curled, more air had been drawn it, the combustion spread inside, and high temperature areas increased. At the same time, the size of the shockwave bottle first increased and then decreased in the radial direction, also in Figure 9, it can also be observed that this was because the muzzle gunpowder gas density and pressure were decreasing. Figure 11 shows the pressure distribution on the projectile surface at t=3.25ms, t=3.42ms, t=3.55ms and t=3.73ms. Figure 12 shows the pressure changes (left) and temperature-OH changes (right) of the projectile surface along the x direction at the corresponding moments. In general, the pressure on the surface of the projectile decreased, when the projectile just came out of the muzzle, the pressure on its head was relatively high; in process of the projectile passing through the gunpowder gas flow field, the pressure on its tail was relatively high; when the projectile had completed passed through the flow field, the pressure on its head was relatively high. The temperature and OH distribution curves showed that when the projectile intersected with the front end of the shockwave bottle, intense combustion reaction occurred on the surface of the projectile, while at other moments, the combustion reaction was weak or there’s no combustion reaction at all.

Figure 9. Muzzle flow field temperature cloud at the y=0 section and OH distribution cloud at the z=0 section

Figure 10. Temperature cloud of YOZ section at 9D on the right side of the muzzle

Figure 11. Pressure distribution on the surface of the projectile

Figure 12. Pressure changes (left) and temperature-OH changes (right) of the projectile surface along the x direction at different time moments

5. Conclusion

The mesh-less method has received wide attention due to its inherent superiority, flexibility, and variability. However, when dealing with complex flow fields with moving boundaries, its performance is far less mature than the unstructured mesh method; especially in terms of 3D problems, its point cloud calculation and flux calculation are more complicated, therefore, few studies have been conducted on it. Targeting at the problem of topological structure changes in the point cloud caused by the large displacement movement of the projectile in the three-dimensional muzzle flow field, this paper successfully designed a 3D dynamic mesh-less algorithm. The dynamic point clouds near the large displacement boundary were processed using the local reconstruction method, and the weighted-point filling strategy was adopted to generate the new point clouds. The convective flux term of the ALE-form Euler governing equation was calculated using the AUFS format, the time-division method was introduced to solve the stiff problem of the chemical reaction calculation, and the terms in the governing equation were decoupled. The calculation of the flow field around the high-speed flying projectile with the angle of attack verifies the effectiveness of the algorithm. At last, the proposed mesh-less method was applied to calculate the high-speed non-equilibrium reaction muzzle jet flow field with high-speed moving projectile. The calculation results showed that the complex shock wave structure in the flow field had been captured, and the non-equilibrium chemical reaction was effectively simulated in the study, exhibiting the phenomenon of the second muzzle flame. The multi-scale problem of time and space in the flow field calculation had been solved, indicating that the proposed LSM method was feasible, and this paper provided a new solution for the high-speed non-equilibrium reaction jet problems.

Acknowledgement

This work is supported by the National Natural Science Foundation of China (Grant No.: 51806095), the Natural Science Foundation of Jiangsu Province of China (Grant No.: BK20181022) and the research fund of Nanjing Institute of Technology (No.: ZKJ201702).

  References

[1] Ziegler, J.L., Deiterding, R., Shepherd, J.E., Pullin, D.I. (2011). An adaptive high-order hybrid scheme for compressive, viscous flows with detailed chemistry. Journal of Computational Physics, 230(20): 7598-7630. https://doi.org/10.1016/j.jcp.2011.06.016

[2] Sun, W., Gou, X., El-Asrag, H.A., Chen, Z., Ju, Y. (2015). Multi-timescale and correlated dynamic adaptive chemistry modeling of ignition and flame propagation using a real jet fuel surrogate model. Combustion and Flame, 162(4): 1530-1539. https://doi.org/10.1016/j.combustflame.2014.11.017

[3] Santillana, M., Zhang, L., Yantosca, R. (2016). Estimating numerical errors due to operator splitting in global atmospheric chemistry models: Transport and chemistry. Journal of Computational Physics, 305: 372-386. https://doi.org/10.1016/j.jcp.2015.10.052

[4] Yee, H.C., Shinn, J.L. (1989). Semi-implicit and fully implicit shock-capturing methods for nonequilibrium flows. AIAA Journal, 27(3): 299-307. https://doi.org/10.2514/3.10112

[5] Li, Z., Ferrarotti, M., Cuoci, A., Parente, A. (2018). Finite-rate chemistry modelling of non-conventional combustion regimes using a Partially-Stirred Reactor closure: Combustion model formulation and implementation details. Applied Energy, 225: 637-655. https://doi.org/10.1016/j.apenergy.2018.04.085

[6] Saghafian, A., Terrapon, V.E., Pitsch, H. (2015). An efficient flamelet-based combustion model for compressible flows. Combustion and Flame, 162(3): 652-667. https://doi.org/10.1016/j.combustflame.2014.08.007

[7] Mei, Z., Shi, C., Fan, X., Wang, X. (2020). Numerical simulation of hypersonic reentry flow Field with gas-phase and surface chemistry models. Materials Today Communications, 22: 100773. https://doi.org/10.1016/j.mtcomm.2019.100773

[8] Sussman, M. (1994). A computational study of unsteady shock-induced combustion of hydrogen-air mixtures. In 30th Joint Propulsion Conference and Exhibit, 3101: 1-23. https://doi.org/10.2514/6.1994-3101

[9] Shuen, J.S., Yoon, S. (1989). Numerical study of chemically reacting flows using a lower-upper symmetric successive overrelaxation scheme. AIAA Journal, 27(12): 1752-1760. https://doi.org/10.2514/3.10331

[10] Suresh, A., Liou, M.S. (1992). The Osher scheme for non‐equilibrium reacting flows. International Journal for Numerical Methods in Fluids, 15(2): 219-232. https://doi.org/10.1002/fld.1650150206

[11] Harlé, C., Carey, G.F., Varghese, P.L. (2000). Analysis of high speed non-equilibrium chemically reacting gas flows. Part I. Physical modeling and sensitivity studies. International Journal for Numerical Methods in Fluids, 32(6): 669-690. https://doi.org/10.1002/(SICI)1097-0363(20000330)32:6%3C669::AID-FLD982%3E3.0.CO;2-W

[12] Wu, W., Xu, H.Q., Wang, L., Xue, R. (2014). Numerical simulation for the external combustion of base-bleed projectile using gridless method. Journal of Harbin Institute of Technology (New Series), 21(3): 95-100.

[13] Wang, L., Xu, H., Wu, W., Xue, R. (2015). Numerical simulation of muzzle flow field based on calculation of combustion productions in bore. Journal of Harbin Institute of Technology (New Series), 22(4): 72-78.

[14] Kirk, B., Stogner, R., Oliver, T.A., Bauman, P.T. (2013). Recent advancements in fully implicit numerical methods for hypersonic reacting flows. In 21st AIAA Computational Fluid Dynamics Conference, 2559: 1-30. https://doi.org/10.2514/6.2013-2559

[15] Liu, G.R., Karamanlidis, D. (2003). Mesh free methods: moving beyond the finite element method. Applied Mechanics Reviews, 56(2): B17-B18. https://doi.org/10.1115/1.1553432

[16] Jiang, T., Chen, Z.C., Lu, W.G., Yuan, J.Y., Wang, D.S. (2018). An efficient split-step and implicit pure mesh-free method for the 2D/3D nonlinear Gross-Pitaevskii equations. Computer Physics Communications, 231: 19-30. https://doi.org/10.1016/j.cpc.2018.05.007

[17] Idelsohn, S.R., Onate, E., Calvo, N., Del Pin, F. (2003). The meshless finite element method. International Journal for Numerical Methods in Engineering, 58(6): 893-912. https://doi.org/10.1002/nme.798

[18] da Silva, L.F., Azevedo, J.L., Korzenowski, H. (2000). Unstructured adaptive grid flow simulations of inert and reactive gas mixtures. Journal of Computational Physics, 160(2): 522-540. https://doi.org/10.1006/jcph.2000.6470

[19] Bussing, T.R., Murman, E.M. (1988). Finite-volume method for the calculation of compressible chemically reacting flows. AIAA Journal, 26(9): 1070-1078. https://doi.org/10.2514/3.10013

[20] Lee, J.W. (2007). Parallelized Cartesian grid methodology for non-equilibrium hypersonic flow analysis of Ballutes. Doctoral dissertation, Georgia Institute of Technology, 4548: 1-15. http://hdl.handle.net/1853/16150

[21] Ku, C.Y., Liu, C.Y., Yeih, W., Liu, C.S., Fan, C.M. (2019). A novel space-time meshless method for solving the backward heat conduction problem. International Journal of Heat and Mass Transfer, 130: 109-122. https://doi.org/10.1016/j.ijheatmasstransfer.2018.10.083

[22] Khorasanizade, S., Sousa, J.M.M. (2016). An innovative open boundary treatment for incompressible SPH. International Journal for Numerical Methods in Fluids, 80(3): 161-180. https://doi.org/10.1002/fld.4074

[23] Wu, N.J., Chang, K.A. (2011). Simulation of free-surface waves in liquid sloshing using a domain-type meshless method. International Journal for Numerical Methods in Fluids, 67(3): 269-288. https://doi.org/10.1002/fld.2346

[24] Batina, J. (1993). A gridless Euler/Navier-Stokes solution algorithm for complex-aircraft applications. In 31st Aerospace Sciences Meeting, 0333: 1-11. https://doi.org/10.2514/6.1993-333

[25] Rijas, A.S., Sriram, V., Yan, S. (2019). Variable spaced particle in mesh-free method to handle wave-floating body interactions. International Journal for Numerical Methods in Fluids, 91(6): 263-286. https://doi.org/10.1002/fld.4751

[26] You, C., Qiu, Y., Xu, X., Xu, D. (2008). Numerical simulations of viscous flows using a meshless method. International Journal for Numerical Methods in Fluids, 58(7): 727-741. https://doi.org/10.1002/fld.1760

[27] Kim, S.L., Choi, J.Y., Jeung, I.S., Park, Y.H. (2001). Application of approximate chemical Jacobians for constant volume reaction and shock-induced combustion. Applied Numerical Mathematics, 39(1): 87-104. https://doi.org/10.1016/S0168-9274(01)00054-X

[28] Choi, J.Y., Jeung, I.S., Yoon, Y. (2000). Computational fluid dynamics algorithms for unsteady shock-induced combustion, part 2: Comparison. AIAA Journal, 38(7): 1188-1195. https://doi.org/10.2514/2.1112

[29] Wang, L., Xue, R., Cai, N., Wu, W., Niu, M., Zhang, D. (2020). Application of least squares meshless method in multi-component high-speed non-equilibrium reaction jet application of least squares meshless method in multi-component high-speed non-equilibrium reaction jet. International Journal of Design & Nature and Ecodynamics, 15(4): 507-514. https://doi.org/10.18280/ijdne.150407

[30] Perrone, N., Kao, R. (1975). A general finite difference method for arbitrary meshes. Computers & Structures, 5(1): 45-57. https://doi.org/10.1016/0045-7949(75)90018-8

[31] Liszka, T., Orkisz, J. (1980). The finite difference method at arbitrary irregular grids and its application in applied mechanics. Computers & Structures, 11(1-2): 83-95. https://doi.org/10.1016/0045-7949(80)90149-2

[32] Löhner, R., Sacco, C., Onate, E., Idelsohn, S. (2002). A finite point method for compressible flow. International Journal for Numerical Methods in Engineering, 53(8): 1765-1779. https://doi.org/10.1002/nme.334

[33] Ortega, E., Oñate, E., Idelsohn, S. (2009). A finite point method for adaptive three-dimensional compressible flow calculations. International Journal for Numerical Methods in Fluids, 60(9): 937-971. https://doi.org/10.1002/fld.1892

[34] Erhart, K.J., Gerace, S.A., Divo, E.A., Kassab, A.J. (2009). An RBF interpolated generalized finite difference meshless method for compressible turbulent flows. In ASME International Mechanical Engineering Congress and Exposition, pp. 571-581. https://doi.org/10.1115/IMECE2009-11452

[35] Firoozjaee, A.R., Afshar, M.H. (2011). Steady-state solution of incompressible Navier-Stokes equations using discrete least-squares meshless method. International Journal for Numerical Methods in Fluids, 67(3): 369-382. https://doi.org/10.1002/fld.2370

[36] Munikrishna, N., Balakrishnan, N. (2011). Turbulent flow computations on a hybrid cartesian point distribution using meshless solver LSFD-U. Computers & fluids, 40(1): 118-138. https://doi.org/10.1016/j.compfluid.2010.08.017

[37] Gargari, S.F., Kolahdoozan, M., Afshar, M.H. (2018). Mixed Discrete Least Squares Meshfree method for solving the incompressible Navier–Stokes equations. Engineering Analysis with Boundary Elements, 88: 64-79. https://doi.org/10.1016/j.enganabound.2017.12.018

[38] Hong, W., Hongquan, C., Zhihua, M., Saihu, P. (2009). Gridless method for solving moving boundary problems and its dynamic cloud of points. Journal of Nanjing University of Aeronautics & Astronautics, 41(3): 296-301.

[39] Zhou, X., Xu, H. (2010). Gridless method for unsteady flows involving moving discrete points and its applications. Engineering Applications of Computational Fluid Mechanics, 4(2): 276-286. https://doi.org/10.1080/19942060.2010.11015316

[40] Anandhanarayanan, K. (2010). Development of three-dimensional grid-free solver and its applications to multi-body aerospace vehicles. Defence Science Journal, 60(6): 653-662. https://doi.org/10.14429/dsj.60.583

[41] Zhang, J., Ren, D.F., Ma, X.J., Tan, J.J., Cai, X.W. (2014). A meshless solution method for unsteady flow with moving boundary. Advances in Mechanical Engineering, 6: 209575. https://doi.org/10.1155/2014/209575

[42] Glowinski, R., Osher, S.J., Yin, W. (2017). Splitting methods in communication, imaging, science, and engineering. Springer, 627-641.

[43] Trotter, H.F. (1959). On thehe product of Semi-groups of Operators. Proceedings of the American Mathmatic Society, 10(4): 545-551.

[44] Strang, G. (1968). On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5(3): 506-517. https://doi.org/10.1137/0705041

[45] Gou, X., Sun, W., Chen, Z., Ju, Y. (2010). A dynamic multi-timescale method for combustion modeling with detailed and reduced chemical kinetic mechanisms. Combustion and Flame, 157(6): 1111-1121. https://doi.org/10.1016/j.combustflame.2010.02.020

[46] Brown, P.N., Byrne, G.D., Hindmarsh, A.C. (1989). VODE: A variable-coefficient ODE solver. SIAM Journal on Scientific and Statistical Computing, 10(5): 1038-1051. https://doi.org/10.1137/0910062

[47] Guzzo, F.R., Azevedo, J.L.F. (2010). A new particle-like method for high-speed flows with chemical non-equilibrium. Journal of Aerospace Technology and Management, 2(1): 17-32. http://dx.doi.org/10.5028/jatm.2010.02011732

[48] Sun, M., Takayama, K. (2003). An artificially upstream flux vector splitting scheme for the Euler equations. Journal of Computational Physics, 189(1): 305-329. https://doi.org/10.1016/S0021-9991(03)00212-2

[49] Fedkiw, R.P., Merriman, B., Osher, S. (1997). High accuracy numerical methods for thermally perfect gas flows with chemistry. Journal of Computational Physics, 132(2): 175-190. https://doi.org/10.1006/jcph.1996.5622