Application of a CFD Model in Determination of the Muzzle Blast Overpressure in Small Arms and Its Validation by Measurement

The main subject of this paper is application of a Computational Fluid Dynamics (CFD) model in determination of the muzzle blast overpressure and its physical manifestations, as well as its validation through measurements of primary parameters. Unsteady Reynolds-averaged Navier-Stokes equations (URANS) with a corresponding turbulence model were applied for numerical simulation of complex gas-dynamic process of propellant gases release from the barrel after firing. The unstructured adaptive mesh for spatial discretization was applied, as suitable model for numerical calculation and physical interpretation of these intensive dynamic processes. The provided experimental results were compared with the results of numerical simulations, which were thus validated, according to adopted minor simplifications.


INTRODUCTION
The specific phenomena that occur at the moment of gun fire are the result of high powder gas energy.The particular manifestations consist of high intensity overpressure, high temperature and high flow velocity of gas.These phenomena were able to produce manifestations, even at the distances farther than 100 calibers from the muzzle.The flow field around the muzzle of a gun barrel is complex and includes flow phenomena such as expansion waves, compression waves, shocks, shear layers, and blast waves.The gunshot is the combination of a number of acoustic waves formed as a result of four main components: the muzzle blast wave gunpowder gases that flow at supersonic velocities, the shock wave, generated due to the supersonic projectile movement, the wave formed by the air column ejected from the barrel in front of the projectile and the acoustic wave generated by impacts of weapon parts during the firing process.
In the paper of Zhangxia et al. [1], which deals with the muzzle blast formed at firing process of 35 mm antiaircraft gun, the authors stated that, while a projectile accelerates at high pressure and high temperature, the powder gases are released, which creates the noise following the muzzle blast wave.As the muzzle energy increases, the level of impulse noise also increases proportionally.The firing noises can be dissolved into muzzle blast, projectile noise, and noise of explosion when the projectile hits the target [2].
The powder gas pressure is manifested as a muzzle blast wave.The muzzle blast consists of two main properties that distinguish it from general noise: high amplitude and short duration.For large caliber guns, the positive impulse duration is a few milliseconds, while for small caliber guns it is less than 0.5 milliseconds.
The overpressure of the muzzle blast wave endangers both the shooter and the equipment, and the created cloud of dust interferes with the aiming and decreases accuracy.That is why spreading of the blast wave has to be regulated using muzzle gas devices.Most scientists use empirical and analytical methods.
The CFD calculation [1] was based on the axially symmetric, standard κ-ε turbulence model and Navier-Stokes equations.Muzzle blast overpressure values with a muzzle brake at different points of the environment were compared to the values obtained without a muzzle brake.The CFD analysis of the muzzle blast overpressure revealed that in the model with a muzzle brake, a great amount of propellant energy had been dissipated and that approximately 74% of the overpressure was reduced compared to the model without a muzzle brake.
This study contributes to better understanding of the blast wave characteristics and offers a good method to model and design muzzle brakes using CFD.The paper [3] studies the application of the axially symmetric CFD implicit model based on LES (Large Eddy Simulation) model of turbulent viscosity for designing different types of silencers, as well as the application of the FW-H (Fowcs-Williams and Hawking) acoustic model for sound waves frequency analysis.In this paper, the authors stated that the muzzle blast wave caused by the gas expansion is the main acoustic source.The internal energy formed by a chemical reaction of the propellant gases which spread rapidly from the barrel, generates muzzle blast and flash in only few milliseconds.The discharged propellant gas generates the normal wave and the blast wave.The supersonic projectile forms an acoustic shock wave that spreads away from the bullet path.The shock wave expands in a form of a cone behind the bullet, while the top of the wave, which moves at the speed of sound, goes outside the front.
As a rule, the impulse noise of small caliber weapons concentrates in the frequency range of 500-1000 Hz, while the impulse noise of large caliber weapons is in the lowfrequency range of less than 200 Hz [4].Numerous researches related to the design of suppressors are directed towards changing the frequency of gunshot noise [5][6][7].
The paper [8] studies designing of large weapon silencers and compares new data with previously obtained results.For assessment, design and simulation of the results, Gambit and Fluent software is used.The simulated results of pressure and sound pressure levels at different points inside the silencer and at different points in the environment are compared with the results at the same points obtained without the use of silencer.CFD method is used to analyze the impulse noise pressure and reduction of sound pressure at firing a large caliber weapon.Cayzac et al. used FREIN code based on 2D Euler equations and TVD (total variation diminishing) numerical scheme to simulate the precursor blast wave in the boundary layer, projectile movement and the effect of the gas brake [9].The process of separation of the APFSDS (Armour-Piercing Fin-Stabilized Discarding-Sabot) from the sealing ring of the sub-caliber projectiles is simulated [10].
The research of Kang et al. [11] studies the usage of numerical simulations for comparison of the blast wave reductions when silencers are used in weapon systems of small and large calibers.URANS equations can be used as the basic model for 3D unsteady, compressible and turbulent flow at firing a weapon.These equations can be transformed into a cylindrical coordinate system, which can be applied to supersonic waves and silencer design.Roe's upwind numerical scheme for spatial derivative and central second order schemes for determination of viscosity were applied for analysis.The multi-block mesh technique is used to describe complex muzzle geometry.Numerical grid points at the overlapped regions (part of structured mesh and transition to the unstructured mesh) overlap in order to preserve a conservative flux model in the entire mesh.The boundary conditions are applied to the wall as an adiabatic viscous model (no slip).The boundary conditions applied at the entry and the exit are based on the Riemann invariant model.The experimental and the simulation data are used to define the complex supersonic flow and the impulse noise reduction.The numerical software for the analysis of the high pressure of the supersonic blast wave is developed.Several models of silencer are analyzed through simulations and based on calculations for the three-chamber silencer, the sound reduction is found to about 42 dB.The optimal model for noise reduction at all frequencies is determined based on the actual performed tests.However, the experiments involved only small arms weapons, with no large caliber tests performed.Considering the blast wave characteristics, the model presents a good basis for predictive designing of silencers in large caliber weapons systems.
The primary objective of this paper is to make a simulation model of the blast wave concerning small arms.The simulation model is created in Ansys Fluent software package, using the RANS model for unsteady multiphase turbulent flow and achievable k-ε model of turbulent viscosity for 2D.For multiphase flow, VOF (volume of fraction) model for two phases (powder gases and ambient atmosphere) is used.The experimental parameters are measured and the mean values of the blast wave overpressure are obtained, at the chosen characteristic points in the ambient region in relation to the weapon muzzle.In the paper, the simulation and the measurements are performed for the automatic rifle М21 (Zastava Arms), caliber 5.56 × 45 mm.
The analysis of the obtained results shows that the simulation model can be successfully used to determine the intensity of the blast wave overpressure for different types of small arms.

THE MATHEMATICAL MODEL
The mathematical model is based on the Navier-Stokes differential equation system, which can be conveyed with the following equation, The Eq. ( 1) is the complete system of basic equations in the closed (conservation) form, where U, F, G, H and J are vector columns, ( ) Vector columns F, G and H in Eq. ( 1) are the flux vectors, while the Ј vector column is the source i.e. the sink.The vector column U is called the solution vector.The first elements in all the vector columns are the continuity equations.The second elements in the vector columns are motion conservation equation systems.The last elements in the vector columns are governing equations systems of energy sustainment.
If Eq. ( 1) is transformed in the time domain, a form suitable for describing unsteady changes -flows is obtained.It will be applied in this paper, The variable U in Eq. ( 7) is a solution vector, since the elements in the vector column U (ρ, ρu, ρv, etc.) are determined numerically in the time domain.The variables on the right side of the Eq. ( 7) are determined based on the known values for a previous time interval.The values ρu, ρv, ρw and ρ(e+V 2 /2) depend on ρ and they are called flux variables.The variables u, v, w and e are primitive variables.If the starting values of the primitive variables and the variable ρ are known, the rest of the values of the primitive variable can be calculated using the expressions, ( )

Multiphase Fluid Flow
For the proper mathematical interpretation of the powder gases flow in the ambient region, one should keep in mind that this is a multiphase flow, with powder gases in the first phase and atmosphere gas in the second phase.It is a common practice to have two flow phases for two ideal gases which do not interact due to short time of gas discharge [12][13][14][15][16].For multiphase flows, two basic scalar models of transferable equation are generally used: one for independent phases and the other one for the mixture.For the undetermined scalar value k in the phase l, the transferable equation (defined by the function k l φ ) for the control volumes that refer only to the l phase, is, ( ) where are α l fraction volume, ρ l density, l u l phase velocity, For the mixture, the value of the scalar function k l φ should be determined by defining all the above given expressions for all the phases.Calculation based on volume fractions of the phases or the VOF model is used for numerical calculation of timedependent solutions [14][15][16].Stability of the VOF calculations lies in the fact that the solution is independent of the formulations.In each control volume, the volume of fractions is considered unique.The values for all the valuables are the volume mean value in all the parts of the grid where the location of each phase is known.If the cells are filled with one kind of fluid, the values refer only to that particular fluid.In the boundary part, the values are determined based on the volume fraction of fluids of the two phases in it.
In other words, if the q th volume fraction of the fluid in the cell is marked as α q , then the following three conditions are possible: -α q = 0, the cell is empty (fluid q th ), -α q = 1, the cell is full (fluid q th ), -0 < α q < 1, the cell contains a fraction of the fluid q th and a fraction of one or more fluids.
Based on the local value of α q , the corresponding properties and variables will be assigned to each control volume in the domain.
The tracking of the interface S between the phases is accomplished by the solution of a continuity equation for the volume fraction of one (or more) of the [14][15][16].For the phase q th , this equation has the following form, ( ) ( ) ( ) where qp m  -mass transfer from the phase q into the phase p, and pq m  -mass transfer from the phase p to the phase q.
As a rule, the source on the right side of the Eq. ( 15), Sα q equals zero, but the values of the constants of the defined mass sources for each phase can be determined.The volume fraction equation cannot be solved for the primary phase; the primary-phase volume fraction will be computed based on the following constraint, The volume fraction equation can be solved either implicitly or explicitly by time discretization.
For the control volume model convective and diffusion fluxes should be calculated in correlation with the sources within it [14,15].For the boundary layers, the volume fractions of the phases are usually calculated using one of interpolation methods [14,15].In this case, a geometric reconstruction scheme has been applied (Fig. 1).
With geometric reconstruction schemes, the special interpolation model is applied for cells near the boundary of two phases, which is shown in Fig. 1а, the actual shape of the boundary between the phases and Fig. 1b, the layout when the geometric reconstruction method is used.The explicit and implicit schemes are used for the boundary cells which are interpolated as if they were totally filled with one of the phases (using a QUICK model scheme, Upwind of the first or second order, modified HRIC, compressive or CICSAM scheme model), before they undergo a special treatment [14][15][16].The geometric reconstruction scheme is primarily used for unstructured mesh [14][15][16].With this scheme, the boundary between two phases is interpreted by a linear inclination inside each cell.

THE SIMULATION MODEL
Ansys Fluent is a software suitable for creating simulation models of gas-dynamic processes of powder gas flow out of the barrel, during the fire of weapon.It should be pointed out that Ansys Fluent can be used to numerically calculate the change of the blast wave of the discharged powder gases in the time domain for any of the spatial discretization points.In theory, Ansys Fluent has no limitations.If the already implemented models do not yield solution accurate enough, new models can be defined.In practice, limitations are related to available space and that is why only the axial symmetric 2D model is studied in this paper.Simulation models for the powder gas discharged from the weapon barrel without a muzzle gas device are created.In order to create a mesh, a circular space of the 4500 mm diameter is used for the 2D model, as found in the literature [12,13,17].The boundary surface of the fluid entry is 4 mm diameter inside the barrel, so that the influence of the barrel edge on turbulence i.e. fluid whirling for both simulation models would not be avoided.A hybrid mesh (Fig. 2) is created.It consists of a part of structured mesh inside the barrel and a part of unstructured mesh of the ambient field.
Spatial discretization is performed using 359.044cells, with 198.562 nodes and 593.606 cell faces.
The mesh is "too rough" for this spatial discretization model, therefore, automatic pressure gradient adaptation is performed in pre-processing.An unsteady (transient) model is defined using Navier-Stokes equations based on the pressure, according to the RANS model, while the velocity is calculated as absolute.
A complete discretization of the ambient field is shown in Fig. 3.
In addition to the basic model, other equations also need to be defined.Since this is a two-phase flow, in which powder gases spread into the atmosphere, the VOF model is applied with the energy conservation equation.Powder gases are the primary phase, while air is the secondary phase.The viscosity model is also defined, and, in this case, it is the realizable κ-ε model.Due to the impulse nature of the process, ideal gases with no mutual interactions are considered.The change in the initial powder gases pressure based on the mass flow is defined with a change function, using the option UDF (user defined function).The UDF initial function was obtained based on an internal-ballistics calculation for the weapon, for the period of adiabatic propagation of the gases after the projectile left the barrel.The accuracy of the internal-ballistics calculation is checked by comparison with the experimentally measured reference values, such as the muzzle velocity and the pressure changing within the barrel in the function of time (by applying pressure sensors and test barrels).In this case, the muzzle velocity of the projectile with standard ammunition and a standard rifle was measured.The measured average muzzle velocity was 915.44 m/s, and the calculated value was 916.36 m/s.Based on this parameter it is confirmed that the accuracy of this method is very precise and the method itself is reliable.
The diagram of the pressure change of the gunpowder gases at the muzzle, based on the internal-ballistics calculation is shown in Fig. 4.
To easily calculate and avoid higher order elements of integration in the UDF in the C++ code, the constant function is set directly, as shown in Fig. 5.
The boundary conditions of the ambient field are defined based on the output pressure.The standard values of atmospheric pressure and temperature are taken.For the values inside the cells the pressure and velocity are calculated together, i.e., their values are correlated inside the numerical mesh.The pressure values are determined based on the values of volumetric forces.The density and the amount of movement are determined based on the first order upwind change scheme, while the volume of the phases is determined based on the geometric reconstruction model inside the cells.During the numerical calculation, remainder values of the polynomial functions (for which the control criterion has already been defined) are monitored.The monitoring of the pressure values at characteristic points has been defined as an additional parameter (Fig. 6).
Due to the nature of blast waves, a time step has to be short.The number of iterations increases the calculation accuracy, but if the time step is short enough, a small number of iterations can be taken in order to shorten the time needed for calculations and in order to avoid processor loading.Time discretization for simulation models is set as time step of 2×10 −7 second and five iterations within step.

The Results of Numerical Model
The values of pressure at the characteristic points are the key parameters that can be compared with the measured values.For this simulation model, the set of data is obtained for 1200 time intervals.Fig. 7 shows the function of the pressure change at characteristic points per one direction for different equidistance as a function of time.Since the moment it is formed, the blast wave spreads concentrically, [18], as it is shown in Fig. 8, [19].Fig. 9 shows visualized propagation of the powder gases blast wave using CFD simulation.The visualization shows that the simulation results correspond to the physical process to a great extent.In the simulation model, the effect of the blast wave seen in Fig. 8 is not taken into account since this influence is rather small and does not have significant effect on maximal values of the pressure at the referential points.

THE EXPERIMENTAL MEASUREMENTS
For the detection of the pressure changes in fast gasdynamic processes, it is suitable to use measuring systems based on piezoelectric effect.In order to get a complete physical pattern of this impulse phenomenon, several sensors are used as seen in the block diagram (Fig. 10), [20].Sensors PCB Piezotronics 137А23 (РР in Fig. 10) are used to measure the overpressure, while the charge amplifier 494A21 PCB Piezotronics (PS in Fig. 10) is used as an amplifying unit [20].The results are recorded in the LabView software using the NI 6008 card.
The linearity, i.e. the maximum measurement error, is defined in accordance with specifications of the measurement converter PCB 137А 24, and it was less than 1%.By adding up the overall errors given by the equipment specifications, it can be concluded that the registered value of the overpressure change is up to 1%.Therefore, these measurement system records the changes of the blast wave overpressure with high accuracy, and with this overall error, can be considered reliable.The overpressure is measured at referential points, according to the scheme identical to the simulation one (Fig. 6).For each direction, the set of data is obtained based on which mean values are adopted.Fig. 12 shows the change in the referential points per one direction.
Taking into consideration the distribution of the referential points (Fig. 6) and the phenomenon around the muzzle during firing, it can be concluded that only the overpressure blast wave is recorded.The effect of the Mach disc, which has influence at smaller angles of the projectile motion, is not taken into account.The position of the measurement sensors is given in Fig. 12.Ten measurements are taken for each referential point, at the temperature of 296 K and at the atmospheric pressure of 99600 Pa.The measurements are taken at the proving ground test facilities of Serbian Armed Forces.
The mean value of standard deviation of the maximum measured values for each referential point is 1.023%.Based on the deviations for the maximal recorded pressure values and based on the form of change in Fig. 11, it can be concluded that the given measured results accurately reflect the phenomenon of the overpressure muzzle blast wave.

THE ANALYSIS OF THE RESULTS
The quality of the described model can be assessed based on the measured parameters.As already stated, the intensities of the blast wave overpressure of the powder gases at characteristic points around the source, i.e. around the muzzle, are measured.The temperature was 296 K and the atmospheric pressure was 99600 Pa.If the mean measured overpressure values are increased for the atmospheric pressure, they can be compared to the values obtained at the points defined by CFD simulation.The comparison of the maximum pressure values at characteristic points is given in Tab. 1.The mean deviation of the simulated values of maximum pressure values at referential points in relation to the mean maximum measured values for all directions and distances is about 3.7%.The greatest deviation is for the direction 45° at the distance of 0.2 m and it is 14.8%, while all other deviations are less than 10%.
Besides, maximum pressure values, continual pressure changes at characteristic points can be also compared, as shown in Figs. 13 to 15. simplifications and neglects.The simulation model neglects the interaction between the powder gases and the environment due to the impulsive nature of the process.The influence of the projectile, which can have a minor influence at short distances and small angles, is also neglected.The transitional cones at the back cross-section of the barrel are also neglected for the purpose of the simplification.This can significantly influence the pressure intensity values for smaller equidistance and smaller angles of the propagation direction in relation to the direction in which they are discharged from the barrel.

CONCLUSION
The paper presents comparative analysis of the measured muzzle overpressure powder gases values and the values simulated around the muzzle at referential points.
The mean deviation of the simulated values of the maximum pressure values at referential points compared to the mean maximum measured values for all directions and distances is about 3.7%.
If special attention is paid to the frequency of change, wider waves with lower frequency can be noticed in the measured values, which is often due to the noise during or the effect of the environment.The noise and the error due to the wind and the vibrations are possible, because the tests are not taken indoors.The barrel geometry is significantly simplified, which can also have a considerable influence on the pressure at referential points at shorter equidistances and for smaller direction angles.Some approximations have been made in defining the pressure change function during the discharge of the powder gases from the barrel, based on the internal ballistics calculations.
Taking into consideration all these simplifications and approximations, it can be concluded that the simulation model makes a good starting point for determination of the gas dynamic parameters of the blast wave after the firing.In addition to a good agreement in intensity, frequency propagation of the pressure in relation to the distance is also in a good compliance.The visualization shows that the simulation results considerably correspond to the physical process of a weapon firing.
In addition to calculations of changes in the primary parameter (powder gases blast wave pressure), it is possible to calculate changes in all other parameters, like density, temperature, gas velocity etc.
The simulation model is open to improvements.In addition to education purposes, it can be used as a good starting point for designing of gas devices.It can be applied to both small arm and large caliber weapons.

Γ
diffusion coefficient and k l S basic expression, which must be known in advance.In this case k l φ refers only to one phase (the l phase) and it affects only the variables in the l phase.The mass flux for the l phase is defined by the equation,

Figure 1
Calculation using boundary interpolations: a) the boundary between two phases b) the boundary shown using a geometric reconstruction scheme

Figure 2 AFigure 3 A
Figure 2 A hybrid discretization mesh of a part of the ambient field around the weapon barrel

Figure 4 Figure 5
Figure 4 The function of the gunpowder gases pressure at the muzzle

Figure 6
Figure 6 Defining changes in the powder gas pressure at the muzzle

Figure 7
Figure 7 Pressure values at characteristic points for one direction The diagram shows that the function of the pressure change corresponds to the physical model of the process.The distances of the points from the source are given.The points D1, D2, D3 and D4 are at the distance of 0.2 m, 0.4 m, 0.6 m and 1.0 m from the back cross-section of the muzzle, respectively.

Figure
Figure 8The photo of the powder gases blast wave during the firing, [19]

Figure 9
Figure 9 Values of the blast wave pressure in 2D for the interval from 2.3×10 −5 s to 1.6×10 −3 s

Figure 10
Figure 10 Block diagram for measuring and collecting data.(PP -pressure converter, PS -intermediate unit, KV -coaxial lead, NP -data carrier)

Figure 11
Figure 11The recordings of the blast wave overpressure for one direction

Figure 12
Figure 12 Position of the measurement probes

Figure 13 Figure 14 Figure 15
Figure 13 Changes in the blast wave pressure in the direction of 45° in the time domain

Table 1
The Comparison of maximum blast wave pressure values at characteristic points