Three-Dimensional Numerical Modeling of Internal Ballistics for Solid Propellant Combinations

The processes that take place within the Internal Ballistics cycle of an artillery round are highly influenced by geometric effects. They are also highly affected by the presence of a combination of energetic materials, such as the propellant, igniter, primer, and the combustible cartridge cases. For a more realistic simulation of these phenomena, a multidimensional and multicomponent numerical model is presented, based on adaptations and improvements of previous models of conservation equations, maintaining a two-phase, Eulerian–Eulerian approximation. A numerical method based on Finite Volumes and conservative flux schemes (Rusanov and AUSM+), with the ability to predict detonation effects, is proposed. As a result, a versatile 3D numerical code was obtained that was tested in the simulation of artillery firing with conventional and modular charges (MACS). Results show the code is able to characterize the heat and mass transfer of the different energetic materials during the combustion of the propellant and the cartridge cases, the gas expansion, and the projectile acceleration.


Introduction
The international panorama regarding propelling charges for artillery howitzers has changed in recent years with the invention, development, and qualification of Modular Artillery Charge Systems (MACS).
Year after year, the armies of more and more countries adopt this projection system for their ground artillery units, totally or partially, due to the great versatility it offers. Traditionally, in field artillery, since the large-scale introduction of smokeless powders at the beginning of the 20th century and, especially since the World Wars, the generalized use of projection charges in cloth bags became a global phenomenon. The method to achieve different range zones using these projection charges consists of the use of different models of conventional charges (M4A2, M119A1, M203, etc.), which contain different propellants in different quantities, in accordance with the weapon system. Another option, when conventional charges are composed of several individual bags, is the addition/elimination of some of these units before introducing the charge into the chamber. On many occasions, these remaining bags must be eliminated on the battlefield or proving ground itself, with the several issues that this entails. On the contrary, MACS are composed of modules formed by a solid casing, made of a combustible energetic material (combustible cartridge case), which always contains the same amount of propellant, and a central hole that allows the rapid transmission of the initiating flame throughout the chamber [1,2].
From a logistics point of view, MACS are very advantageous. Since all modules are identical (for uni-modular charge systems [3]), MACS makes it easier to plan operations UXGun3D can simulate Internal Ballistics problems with a number of different energetic materials, as in the case of MACS, from a fully three-dimensional perspective. The starting hypothesis is that all the components of the charges can be modeled from the Eulerian approximation. The mathematical and numerical models used by the software are presented below, followed by application of the models to the resolution of a MACS case. Finally, the main conclusions of the study are presented.

Physical Model and Numerical Methods
The Internal Ballistics problem of a field-artillery shot with modular charges is particularly complex since, in addition to the complications inherent in multiphase, transient deflagration/detonation modeling, the presence of multiple energetic materials is added, including the igniter, the combustible cartridge cases, and the propellants. That is, for a bi-modular system, modules can contain two main types of gunpowder. Even for a uni-modular system, there might be a separated kind of gunpowder located around the central hole, with greater vivacity than the main propellant. Therefore, for generalization, a system of conservation equations that is valid for a number n of energetic materials is proposed below.
The added value of the mathematical model presented in this work relies on its capability to account for the thermo-physical properties and combustion laws of different propellants. In other words, the numerical model presented in this study is able to account for different granular solid materials in specific locations within the gun chamber and analyze their transient combustion and combined interaction.
On the other hand, the propellant is confined in modules, which in turn are surrounded by air within a cylindrical spatial domain (i.e., the gun chamber). The volume of the chamber is limited by the breech and the base of the projectile, which eventually shift during the burning of the charge, due to the increase in pressure. Therefore, the geometry of the computational domain needs to be variable. In addition, the change in geometry is coupled with the variables of the problem.
Therefore, the mathematical model adopted includes the conservation equations of mass, momentum, and energy for the gas phase (variables denoted with the subscript "1") and the solid phase (variables with the subscript "2"), composed of n types of materials. This results in the following equations: ∂ ∂t (α 2 ρ 2 U 2 ) + ∇ α 2 ρ 2 U 2 U 2 + p 2 I − p 2 ∇(α 2 ) = −Γ c U 2 + F D ∂ ∂t (H s1 ) + U 2 ·∇(H s1 ) = κ 1 . q t,1 (. . .) ∂ ∂t (H sn ) + U 2 ·∇(H sn ) = κ n . q t,n (7a) where Equations (1)-(4a) represent the gas phase sub-system, and Equations (5)-(9a) represent the solid sub-system. Symbols are defined in the Nomenclature section. The gas phase has unique mass, momentum, and energy conservation equations for the m gaseous species, Equations (1)-(3), while n = m − 1 equations are considered for the conservation of the mass fractions Y 1j , Equations (4) and (4a). Combustion modifies the relative concentration of each species through the chemical source terms ω n . For the solid phase, one mass and one momentum conservation equation is included, Equations (5) and (6), but n = m equations model the energy conservation of the m solid species, to account for the different thermal properties of the energetic materials involved. Equations (8) and (8a) are the balanced equations for the cell-averaged burning distance, d qj (also known in the literature as surface regression length) of solid j, whose definition and limitations are not general, but depend on the form functions and geometry of the solids. Therefore, through the coupling to d qj , the equations include constitutive laws for the averaged reduction of the size of all propellant grains contained in each fundamental volume cell due to the combustion process. Finally, Equations (9) and (9a) are the n = m − 1 conservation equations for the mass fractions of solid species Y 2j . The ignition of the propelling charge takes place by means of the primer, composed of a pyrotechnic mixture that at the initial moment is degraded by an exothermic reaction into the ignition gases. This chemical reaction is not modeled. Instead, in the cells of the computational domain next to the chamber breech, a mass rate of ignition gases (Γ ign ), a momentum associated with these gases (Γ ign U ign ), and the thermal energy of ignition gases (Γ ign Q ign ) are added to Equations (1)-(3). Mass and heat transfer processes inherent to the system of equations transport this energy into adjacent cells, and this initiates the combustion of the propellant.
In this model, an assumption has been made that each of the gaseous species corresponds to the combustion gases of each of the solid species. In this way, when a solid component degrades into gas, at a rate equal to −ω j , it disappears from the solid phase, and Γ c = ∑ n j=1 ω j is added to the mass conservation equation for the gas phase, while ω j is added to each one of the gas species. The reaction rate of each solid j is given by the burning velocity . r j and the number of solid particles per unit volume, so the total mass rate source term is defined by The Nobel-Abel equation of state is used to model the thermodynamic behavior of the gas phase as an approximation to real gas: where γ 1 is the ratio of specific heat of the mixture of n gases, and η 1 is the co-volume of the mixture calculated by η 1 = ∑ n j=1 (η j Y 1j ). The solid phase is assumed to be incompressible, so no equation of state is needed for the solid. This hypothesis is one of the identifying characteristics of Gough's model [15], and therefore it is assumed also in the conservative version of the model developed by the authors in a previous work [10] and used in this paper. Of course, considering the solid phase incompressible is a simplifying hypothesis, but it has proven to be valid for Internal Ballistics problems [17]. It is also very useful from the numerical point of view. Even so, the compaction of the granulated solid propellant j can be modelled by the intergranular stress Λ j (α 2 , ρ 2 ), defined as The average intergranular stress for all the solids in a control volume can be calculated where σ j is S pj multiplied by the number of particles of propellant j per unit volume: This way, although different pressures are defined for each phase, the model can be seen as a one-pressure model, since p 2 = p 1 + Λ(α 2 , ρ 2 ). The intergranular stress models, more specifically, the response of a bed of granulated solid, against variations in the volumetric fraction of the solid phase α 2 = 1 − α 1 , in comparison with a characteristic porosity value, called critical porosity, α c , which is constant for each propellant. When the void fraction α 1 is greater than the critical porosity, the grains are not compacted, and p 2 = p 1 . On the contrary, when there is compaction, the intergranular stress term acts by increasing the pressure of the solid phase, as a source of momentum, aiming to release the stresses in the solid. The Λ formulation of Equation (12) is taken from Gough [15], but other definitions were discussed and tested in a previous work [10]. A further discussion about the effects of intergranular stress in the modeling of propellant compaction is given by Nussbaum [17]. An experimental correlation by Ergun [18] for packed solid beds is used to model the interphase drag F Dj for solid species j: where we assume The total drag force for all the solid species, averaged in a control volume, is simply the contribution of the n solid species: Similarly, the interfacial heat exchanged between phases is modeled by Newton's law of cooling, extended to the average grain surface in a control volume: where T 2j,s is the temperature at the surface of solid j, the heat transfer coefficient is h j = Nu · λ j · S pj / 6V pj , and λ j is the thermal conductivity of the gas j. An experi- mental correlation is used for the Nusselt number: Nu = 2 + 0.4Re 2/3 p,j Pr 1/3 j . Note that T 2j,s = T 2 in Equation (17). In fact, T 2j,s is calculated by a polynomial law for temperature distribution within the solid grain [17], presented and discussed in a previous work [10].
In this case, the approach developed by Gough is used for the averaging of the continua within the spatial domain of the chamber [15]. A more detailed description of the physical model used can be found in [10].

Numerical Methods
The numerical method used to solve the system of Equations (1)-(9a), with the closure laws of Equations (10)- (17), is based on Finite Volumes. First, the system of governing equations can be rearranged and expressed in vector form, as a function of the vector of conserved variables U: where H(U) = [F G H] is the tensor of convective flux, and S(U) is the vector of source terms, which includes some non-conservative terms: . . .
In order to obtain a numerical solution to this problem, a splitting method can be used in such a way that the system is solved in different steps. In the first place, the homogeneous part of the system, the advection problem, dU/dt + ∇ · H(U) = 0 is integrated over a control volume Ω: where A Ω is the boundary of Ω, andn is the vector normal to surface A Ω . Considering the first integral as a time-rate of change of the averaging of the conserved variables U and the boundary A Ω formed by N surfaces so that A Ω = ∑ N s=1 A Ω(s) , Equation (21) can be written as follows: ∂U ∂t Adopting U n+1 j − U n j /∆t for the time derivative of the conserved variables, and considering the rotational properties, the surface integral of the flux is approached by This gives the explicit finite-volume scheme used in this work for the convective part of the system of equations: where A Ω(s) is the area of the surface s = (1, . . . , N) bounding the volume Ω, and T −1 s is the inverse of the rotation matrix. Later, in a second step, the ODE problem corresponding to the vector of source terms is solved: dU/dt = S(U). That is, the update of the vector of conserved variables at time n, in each of the cells of the discretized domain, is carried out through a succession of operators that provide the solution of conserved variables at time n + 1 through a split integration: where L ∆t h represents the integration operator of the hyperbolic advection problem, and L ∆t s is the solution operator of the ODE system. The numerical schemes used to compute the convective fluxes are AUSM+ [19] and AUSM+ up [20] for the gas phase. An extension of the Rusanov scheme [21] was adapted to fit with the particularities of the solid phase [10], and it is also implemented here. To reconstruct the values of U in the cell faces from the values at the centers, a second-order, centered, Gauss linear method is used (i.e., no TVD options are employed, in order to avoid numerical diffusivity).
For the source term integration, the 4th order Runge Kutta method with standard coefficients is used. In this way, although more operations are performed by the code, it allows setting a high CFL value, around 0.70. The code also includes an advanced source-term treatment (ASST), which could serve as a relaxation method in cases when the source-term problem becomes stiff, developed by García-Cascales et al., which is described in detail in [22]. The gradients of α 1 and α 2 of the non-conservative terms included in the S(U) vector are solved by a centered Gauss linear method.
In this new software, UXGun3D, OpenFOAM version 5 was used as the basis for programming. Although UXGun was programmed using object-oriented Fortran, UXGun3D was developed based on OpenFOAM, due to the facilities it provides in handling vector fields because of the libraries available to assist in numerical resolution by finite volumes.

Results and Discussion
This section contains the application of the 3D numerical model to a typical Internal Ballistics problem, based on a 155 mm howitzer firing a conventional charge. A comparison is made between the experimental results measured in the proving ground and the 2D and 3D numerical results, for validation and verification of the code.
At the end of this section, the numerical model is applied to the resolution of a case with MACS.

Conventional Charge Firing
Conventional cloth-bag propelling charges can be modeled as a one-dimensional problem with reasonable accuracy. In fact, for those properly working charges, when combustion takes place without detonation or pressure waves, even zero-dimensional codes are valid to give an approximation for the pressure-time curve, projectile acceleration, and muzzle velocity. This was already addressed with the validation of the 1D model implemented in the UXGun code [10].
The leap forward with the new 3D numerical model (which is implemented in the UXGun3D code) allows more complex scenarios to be faced. The numerical extension of the model to 3D and the use of a new solver developed ad hoc, make a prior verification necessary to ensure this new software can model correctly the shots with charges loaded in bags, with a fully 3D approximation. This version of the code also allows 2D and 2D axisymmetric simulations, just by the use of a mesh with these features and proper boundary conditions in the symmetry axis and in the front/back planes. The 2D and 2D axisymmetric results can be compared with those of the fully 3D simulations, as is done later in this section.
The case selected for this verification is the M4A2 propelling charge, commonly called "white bag". This charge, the main specifications of which can be found in the open literature, is designed to be fired from a 155 mm howitzer and contains approximately 6 kg of M1-type propellant [23]. It is composed of a base load and four bags of different weights, functioning as increments to achieve a wider range. In addition, between the base load and the increments, there is a flash reducer pad that can contain potassium nitrate or potassium sulfate [23]. The charge does not have any perforations or central tube to distribute the flame of the initiator, which is made up of a single sachet of a nitrocellulosebased gunpowder called Clean Burning Igniter (CBI), located at one end of the chamber next to the breach [23]. The ignition occurs in that region, initiated by the primer, which transmits the activation energy to the CBI bag. Later, the flame progressively spreads in a deflagration to the rest of the charge. Despite being a relatively small load, if there are deficiencies in the initiation or in the fire transmission, pressure wave phenomena can occur.
For the case tested, a charge of 6.2 kg of M1 propellant was considered, with a density of 1589 kg/m 3 , manufactured in the form of cylindrical pellets (length 11.46 mm, diameter 5.51 mm) with seven perforations. The igniter contains 14.5 g of CBI. We assume the charge can be three-dimensionally modeled as two sachets (the base load and the rest), separated by the flash reducer pad, which we consider inert.
The charge is fired from an M114 weapon with a chamber volume of 15.04 dm 3 and a barrel length of 2.87 m, with a 43.09 kg M107 projectile. The intrusion of the projectile into the chamber is 115 mm. The burning law of the M1 gunpowder was found through experimental closed vessel tests, following the procedure described in STANAG 4115 [24], obtaining coefficients a = 0.00281 m/(s·MPa n ), n = 0.684, and b = 0. With those coefficients, the burning law of the propellant is assumed to follow Vieille's Law in the following form: .
In Figure 1, a diagram is shown of the charge assembly and the modeled geometry. As can be seen, the igniter bag is represented as a cylinder of constant thickness and volume equal to the lenticular bag. Similarly, the inert sachet is modeled as a propellant-free space of the same volume. We consider the charge to have a uniform diameter of 14 cm and a length of 17 cm for the base charge sack and 28.5 cm for the set of increment bags.
Mathematics 2021, 9, x FOR PEER REVIEW 10 of 23 potassium sulfate [23]. The charge does not have any perforations or central tube to distribute the flame of the initiator, which is made up of a single sachet of a nitrocellulosebased gunpowder called Clean Burning Igniter (CBI), located at one end of the chamber next to the breach [23]. The ignition occurs in that region, initiated by the primer, which transmits the activation energy to the CBI bag. Later, the flame progressively spreads in a deflagration to the rest of the charge. Despite being a relatively small load, if there are deficiencies in the initiation or in the fire transmission, pressure wave phenomena can occur.
For the case tested, a charge of 6.2 kg of M1 propellant was considered, with a density of 1589 kg/m 3 , manufactured in the form of cylindrical pellets (length 11.46 mm, diameter 5.51 mm) with seven perforations. The igniter contains 14.5 g of CBI. We assume the charge can be three-dimensionally modeled as two sachets (the base load and the rest), separated by the flash reducer pad, which we consider inert.
The charge is fired from an M114 weapon with a chamber volume of 15.04 dm 3 and a barrel length of 2.87 m, with a 43.09 kg M107 projectile. The intrusion of the projectile into the chamber is 115 mm. The burning law of the M1 gunpowder was found through experimental closed vessel tests, following the procedure described in STANAG 4115 [24], obtaining coefficients = 0.00281 m/(s·MPa n ), = 0.684, and = 0. With those coefficients, the burning law of the propellant is assumed to follow Vieille's Law in the following form: In Figure 1, a diagram is shown of the charge assembly and the modeled geometry. As can be seen, the igniter bag is represented as a cylinder of constant thickness and volume equal to the lenticular bag. Similarly, the inert sachet is modeled as a propellant-free space of the same volume. We consider the charge to have a uniform diameter of 14 cm and a length of 17 cm for the base charge sack and 28.5 cm for the set of increment bags. When the load is modeled from a 1D view, the bag is considered to occupy the entire width of the gun chamber, since it is not possible to model the directions perpendicular to the axis of symmetry. In 2D and 3D modeling, this changes. Figure 2 shows the numerical results in 2D Axisymmetric and 3D, compared with the experimental pressure-time curve measured by Monreal-Gonzalez et al. [10]. In plots for '2D Axisymmetric' and '3D Axisymmetric' simulations, the load is considered to be concentric with the chamber's axis of symmetry, while in 'Full 3D', the effect of the load resting on the lower part of the chamber is considered. As shown in Figure 2a, the difference is practically negligible. On other hand, Figure 2b shows the differential pressure curves (instantaneous pressure difference between the breech and the opposite end of the chamber) over time. This When the load is modeled from a 1D view, the bag is considered to occupy the entire width of the gun chamber, since it is not possible to model the directions perpendicular to the axis of symmetry. In 2D and 3D modeling, this changes. Figure 2 shows the numerical results in 2D Axisymmetric and 3D, compared with the experimental pressuretime curve measured by Monreal-Gonzalez et al. [10]. In plots for '2D Axisymmetric' and '3D Axisymmetric' simulations, the load is considered to be concentric with the chamber's axis of symmetry, while in 'Full 3D', the effect of the load resting on the lower part of the chamber is considered. As shown in Figure 2a, the difference is practically negligible. On other hand, Figure 2b shows the differential pressure curves (instantaneous pressure difference between the breech and the opposite end of the chamber) over time. This differential pressure is one of the key design variables. Strong negative pressure differences can cause damage to the weapon; therefore, it is important to ensure its correct modeling. However, it was not possible to measure the pressure difference experimentally, since a single pressure gauge was used in the experiments [10].
Mathematics 2021, 9, x FOR PEER REVIEW 11 of 23 differential pressure is one of the key design variables. Strong negative pressure differences can cause damage to the weapon; therefore, it is important to ensure its correct modeling. However, it was not possible to measure the pressure difference experimentally, since a single pressure gauge was used in the experiments [10].
(a) (b) Figure 2. Experimental and numerical pressure evolution with time: (a) pressure evolution at the breech; (b) differential pressure between breech and port. Figure 3 shows the gas volumetric fraction contours for the three geometric approximations previously specified. From the qualitative point of view, the analysis of the numerical results shows how, as the combustion process progresses, the solid phase disappears to become gas, which diffuses throughout the domain, transporting both mass and heat. Although the combustion takes place simultaneously over the entire surface of the charge, since the initiation takes place in the left zone (next to the breech), the combustion process proceeds fundamentally from left to right, until the entire domain is homogenized.

2D Axisymmetric
3D Axisymmetric Full 3D   Figure 3 shows the gas volumetric fraction contours for the three geometric approximations previously specified. From the qualitative point of view, the analysis of the numerical results shows how, as the combustion process progresses, the solid phase disappears to become gas, which diffuses throughout the domain, transporting both mass and heat. Although the combustion takes place simultaneously over the entire surface of the charge, since the initiation takes place in the left zone (next to the breech), the combustion process proceeds fundamentally from left to right, until the entire domain is homogenized.
Mathematics 2021, 9, x FOR PEER REVIEW 11 of 23 differential pressure is one of the key design variables. Strong negative pressure differences can cause damage to the weapon; therefore, it is important to ensure its correct modeling. However, it was not possible to measure the pressure difference experimentally, since a single pressure gauge was used in the experiments [10].
(a) (b) Figure 2. Experimental and numerical pressure evolution with time: (a) pressure evolution at the breech; (b) differential pressure between breech and port. Figure 3 shows the gas volumetric fraction contours for the three geometric approximations previously specified. From the qualitative point of view, the analysis of the numerical results shows how, as the combustion process progresses, the solid phase disappears to become gas, which diffuses throughout the domain, transporting both mass and heat. Although the combustion takes place simultaneously over the entire surface of the charge, since the initiation takes place in the left zone (next to the breech), the combustion process proceeds fundamentally from left to right, until the entire domain is homogenized.

2D Axisymmetric
3D Axisymmetric Full 3D  A quantitative analysis of the results, obtained with the 3D model (i.e., UXGun3D) called 'Full 3D' simulation (see Tables 1 and 2), shows that the relative error of the predicted maximum pressure within the gun chamber (relative to experiments) is less than 1%, while the relative error is less than 1.5% for the maximum velocity of the projectile (i.e., muzzle velocity). Table 1 shows the numerical results for absolute pressure measured by a computational probe located at the end of the chamber next to the breech, as well as the velocity of the projectile. Experimental results by Monreal-Gonzalez et al. [10] were measured by an internal pressure gauge that could freely move along the chamber, so it was impossible to know the exact location of the sensor at any time. In addition, the only experimental velocity available is the muzzle velocity at the bore, measured by Doppler radar. That is the maximum velocity, as the projectile accelerates continuously along the tube. Absolute and relative errors between experiments and simulations are gathered in Table 2. All of them are below 1.3%. Based on these results, we can consider that the performance of UXGun3D for this classic Internal Ballistics problem is satisfactory. Calculations were performed on a shared-memory parallel system with two Intel ® Xeon ® E5-2665 CPUs with eight cores each, and up to 128 GB of RAM. For proper parallel performance, it is worth noting that 3D domains must be decomposed appropriately in order to have balanced sub-mesh sizes, even when the domain grows, due to the projectile's motion. This balance is achieved by decomposing the processor domains longitudinally, to include in all the sub-domains the whole length, including the projectile's moving boundary, where the addition of cell layers is carried out. Note that if these sub-domains are divided transversally, only one domain processor is affected by the new computational cells, and thus this processor will become a bottleneck (similar or even lower performance as in a serial execution). The simulation time required for M4A2 tests is approximately 5 min in the case of 2D domains and 10 min for 3D computations.

Mesh Sensitivity Analysis
A study was carried out to analyze how the spatial discretization size affects the solution, with the same conditions as in the reference case in the previous section. The M4A2 propelling charge firing was simulated using three meshes with different resolutions: a coarse mesh with a cell size of d x = 9 mm, a reference case using the same mesh as in the previous simulations (d x = 6 mm), and a finer mesh, with a cell size of d x = 3 mm. The results (Figure 4 and Table 3) show that, regarding maximum chamber pressure, there is no clear trend. Rather, it is concluded that the maximum pressure is not sensitive to the size of discretization. The size of the element does influence the speed of the process such that, with the finest mesh, the pressure rise occurs earlier, and this is delayed when the mesh becomes coarser. The main conclusion is that little spatial resolution is required, except to correctly define the initial conditions (charge geometry). It is crucial to ensure that the same propellant mass is introduced, regardless of the discretization size.
Mathematics 2021, 9, x FOR PEER REVIEW 13 of 23 except to correctly define the initial conditions (charge geometry). It is crucial to ensure that the same propellant mass is introduced, regardless of the discretization size.

Capabilities of the Model for the Simulation of Modular Charges Firing
We now address the application of the mathematical model to the resolution of a case of MACS firing. In this case, the internal ballistics problem becomes more complex due to the reasons raised in Section 2. The approach is based on a six-module uni-modular MACS (all modules are identical), with a combustible case composed almost entirely of nitrocellulose, which means that it burns completely without leaving residues, and a triple base propellant (nitrocellulose, nitroglycerin, and nitroguanidine) extruded in the form of rosette-shaped pellets with 19 perforations, as shown in Figure 5. The propellant fills the modules in a disorderly way; that is, the grains are not aligned but loaded in bulk inside the cases. The modules contain a small flash reducing bag, which is considered inert for simplicity, and a central hole for flame transmission. Around this, there is an igniter impregnated by black powder.
Solid propellant grains (both rocket motor grains, and pellets for artillery propellants) can be assumed to burn in parallel layers. Therefore, it is not a bad approximation to consider that the normal burn distance is similar on all surfaces around a grain. The burning rate of the propellant (̇= d /d ) normal to the grain surface can be determined for any particular propellant composition by experiments (mainly closed-vessel tests). We assume the burning rate can be modeled by Vielle's Law, as in Equation (26).
The burning rate ̇ is independent of grain geometry, but it has a decisive influence on the burning mode of grains. For the type of propellant considered, the surface of the grain increases with the burning distance, and the volume decreases. Therefore, the burning behavior of the propellant is said to be progressive. These shape functions, and , were obtained from [25].

Capabilities of the Model for the Simulation of Modular Charges Firing
We now address the application of the mathematical model to the resolution of a case of MACS firing. In this case, the Internal Ballistics problem becomes more complex due to the reasons raised in Section 2. The approach is based on a six-module uni-modular MACS (all modules are identical), with a combustible case composed almost entirely of nitrocellulose, which means that it burns completely without leaving residues, and a triple base propellant (nitrocellulose, nitroglycerin, and nitroguanidine) extruded in the form of rosette-shaped pellets with 19 perforations, as shown in Figure 5. The propellant fills the modules in a disorderly way; that is, the grains are not aligned but loaded in bulk inside the cases. The modules contain a small flash reducing bag, which is considered inert for simplicity, and a central hole for flame transmission. Around this, there is an igniter impregnated by black powder.
Solid propellant grains (both rocket motor grains, and pellets for artillery propellants) can be assumed to burn in parallel layers. Therefore, it is not a bad approximation to consider that the normal burn distance d qj is similar on all surfaces around a grain. The burning rate of the propellant ( . r j = dd qj /dt) normal to the grain surface can be determined for any particular propellant composition by experiments (mainly closed-vessel tests). We assume the burning rate can be modeled by Vielle's Law, as in Equation (26).
In order to facilitate the loading of the six modules into the gun's chamber, the combustible case of modules has a geometry with recesses and flaps, so that one module is easily coupled to the next. For 3D modeling, the geometry was simplified as shown in Figure 6a. This modification was carried out while maintaining the internal volume of the module and, therefore, the amount of original propellant and igniter.

Dynamic Mesh Generation
One of the fundamental issues in the three-dimensional modeling of Interior Ballistics is the need for a computational domain whose geometry is variable with time. Indeed, the gun chamber, which is the domain where the problem takes place, is delimited by the cylindrical surface of the chamber, the breech on one side, and the base of the projectile on the other, as represented in Figure 6b. The forcing bands on the projectile ensure the sealing between the projectile and the rifled bore of the barrel. When the pressure increases above the resistive pressure caused by the rubbing of the projectile in the rifled bore, the projectile begins to move and the gas-solid mixture in the chamber occupies that free space. Therefore, the pressure produced by combustion inside the chamber precisely determines the volume of the domain.
This means that, in simulation, the geometry of the computational domain, where the governing equations of the problem are solved, depends in turn on the solution of the problem. This coupling has been solved by programming in OpenFOAM a particular dynamic mesh generator for this specific problem. dq dq dq The burning rate . r j is independent of grain geometry, but it has a decisive influence on the burning mode of grains. For the type of propellant considered, the surface of the grain S p increases with the burning distance, and the volume V p decreases. Therefore, the burning behavior of the propellant is said to be progressive. These shape functions, S p and V p , were obtained from [25].
In order to facilitate the loading of the six modules into the gun's chamber, the combustible case of modules has a geometry with recesses and flaps, so that one module is easily coupled to the next. For 3D modeling, the geometry was simplified as shown in Figure 6a. This modification was carried out while maintaining the internal volume of the module and, therefore, the amount of original propellant and igniter. In order to facilitate the loading of the six modules into the gun's chamber, the combustible case of modules has a geometry with recesses and flaps, so that one module is easily coupled to the next. For 3D modeling, the geometry was simplified as shown in Figure 6a. This modification was carried out while maintaining the internal volume of the module and, therefore, the amount of original propellant and igniter.

Dynamic Mesh Generation
One of the fundamental issues in the three-dimensional modeling of Interior Ballistics is the need for a computational domain whose geometry is variable with time. Indeed, the gun chamber, which is the domain where the problem takes place, is delimited by the cylindrical surface of the chamber, the breech on one side, and the base of the projectile on the other, as represented in Figure 6b. The forcing bands on the projectile ensure the sealing between the projectile and the rifled bore of the barrel. When the pressure increases above the resistive pressure caused by the rubbing of the projectile in the rifled bore, the projectile begins to move and the gas-solid mixture in the chamber occupies that free space. Therefore, the pressure produced by combustion inside the chamber precisely determines the volume of the domain.
This means that, in simulation, the geometry of the computational domain, where the governing equations of the problem are solved, depends in turn on the solution of the problem. This coupling has been solved by programming in OpenFOAM a particular dynamic mesh generator for this specific problem. dq dq dq

Dynamic Mesh Generation
One of the fundamental issues in the three-dimensional modeling of Interior Ballistics is the need for a computational domain whose geometry is variable with time. Indeed, the gun chamber, which is the domain where the problem takes place, is delimited by the cylindrical surface of the chamber, the breech on one side, and the base of the projectile on the other, as represented in Figure 6b. The forcing bands on the projectile ensure the sealing between the projectile and the rifled bore of the barrel. When the pressure increases above the resistive pressure caused by the rubbing of the projectile in the rifled bore, the projectile begins to move and the gas-solid mixture in the chamber occupies that free space. Therefore, the pressure produced by combustion inside the chamber precisely determines the volume of the domain. This means that, in simulation, the geometry of the computational domain, where the governing equations of the problem are solved, depends in turn on the solution of the problem. This coupling has been solved by programming in OpenFOAM a particular dynamic mesh generator for this specific problem.
Since the geometry of the part of the chamber next to the breech is invariable, the best way to tackle the problem of dynamic mesh is considered to be generating mesh cells in the new regions that the projectile leaves free. Thus, the code is prepared for the generation of mesh from the coupling of frustoconical figures (since, as can be seen in Figure 7, the diameter of the domain is variable and axisymmetric). Each conical frustum is defined by an initial radius, a final radius, and two positions of the horizontal axis calculated with respect to the origin. This defines the geometry of the outer boundary of the domain as a function of the position of the projectile. Another variable defines the number of mesh elements in each section, so that the meshing resolution is similar to that of the rest of the domain, and finally, a variable defines the type of the projectile and collects its geometry from a library, where different types of projectiles are defined.
Mathematics 2021, 9, x FOR PEER REVIEW 15 of 23 Since the geometry of the part of the chamber next to the breech is invariable, the best way to tackle the problem of dynamic mesh is considered to be generating mesh cells in the new regions that the projectile leaves free. Thus, the code is prepared for the generation of mesh from the coupling of frustoconical figures (since, as can be seen in Figure 7, the diameter of the domain is variable and axisymmetric). Each conical frustum is defined by an initial radius, a final radius, and two positions of the horizontal axis calculated with respect to the origin. This defines the geometry of the outer boundary of the domain as a function of the position of the projectile. Another variable defines the number of mesh elements in each section, so that the meshing resolution is similar to that of the rest of the domain, and finally, a variable defines the type of the projectile and collects its geometry from a library, where different types of projectiles are defined.  When, in a time step, the pressure of the gas phase is greater than the resistive pressure of the projectile (p 1 > p res ), a new row of cells is generated, the size of which corresponds to the movement of the projectile in a given ∆t. The new cells are initialized with the values of the conserved variables (for each phase: pressure, temperature, velocity, etc.) in the adjacent cell. These values were calculated in the immediately previous time step. Of course, the speed at which the mesh is generated depends on the instantaneous speed of the projectile, in accordance with Newton's Second law as applied to its center of gravity: According to Equation (27), the projectile experiences an acceleration (i.e., du proj /dt > 0) only if the pressure at the projectile's base p m is greater than the resistive pressure P res , which models the wall friction between the projectile and the riffled bore of the howitzer. This value is not a constant, but P res is a function of the x proj coordinate of the projectile position. The resistance is greater in the first stages when the projectile is swaged into the rifling and decreases progressively as it advances along the tube. When the projectile reaches the coordinate corresponding to the muzzle, the domain has reached its maximum possible extension, and the simulation is concluded. The projectile's velocity at that point is the predicted muzzle velocity.

Case Definition and Numerical Results
The case presented here, as a sample of the application of the 3D mathematical and numerical model to the Interior Ballistics problem of modular charges, can be considered an illustrative example with non-real data, but consistent and coherent with what is expected for a MACS of 155 mm. The case is based on the design of the modular charge developed by the company EXPAL Systems S.A., but the geometry of the charge, the weapon, and the propellant parameters have been modified for confidentiality reasons and to safeguard industrial secrecy.
A first simulation, which addresses the firing of five modules in a 155/39 howitzer (Figure 8), allows us to observe the combustion sequence predicted by the code. In particular, contours of the void fraction α 1 are shown, in a representation of a longitudinal section of the fluid domain (since there is axial symmetry, only one of the two halves is represented, for the sake of clarity). It can be seen in the first frame (Figure 8a), which represents the initial conditions, that α 1 = 1 in the central zone (next to the symmetry axis), in the region near the walls of the chamber, and around the bottom of the projectile. Inside the modules, however, α 1 = 0, since the regions occupied by the propellant also contain air, in the gaps between the grains and inside the perforations. In these regions, the value of the void fraction is the minimum possible, the critical porosity, which corresponds to α 1 = α crit = 0.305. Of course, this approximation is the result of the Eulerian conception of the system of equations. Later, once the sequence begins and the combustion progresses (Figure 8b,c), it can be seen how the boundary between the solid and the gas is blurred. In observing the numerical results, it was determined that this phenomenon does not correspond only to the effect of combustion, but also to the numerical diffusion itself that takes place in the solid phase, and this does not respond to the physics of the problem.
Following these results, the numerical method was revised; in particular, revisions were made to the approximate Riemann solvers that were being used to solve the numerical fluxes for the solid phase. The problem came from the coupling of the eigenvalues of both systems of equations. Specifically, in the previous numerical approach, the eigenvalues of both subsystems of equations (gas phase vs. solid phase) were coupled. That is, both in the flux schemes of the gas and of the solid phase, the eigenvalues used corresponded to the gas phase. This aided in the numerical resolution of the problems studied thus far [9,10], with smaller propellant grains and without solid surfaces (combustible cases); in addition, it is a useful procedure for small-caliber shots, where the propellant is formed of tiny grains that move quickly along with the gas phase. In the modular charge problem, however, it was concluded that the coupling of the eigenvalues causes this unwanted diffusive effect. Once the eigenvalues of the flux schemes of each phase are decoupled (i.e., for each phase, the respective eigenvalues are used), the result is clearly better, less diffusive, more realistic concerning the physics of the problem. The characteristic time corresponds now to the combustion process. However, a diffusion effect is still observed in the solid phase that, analyzing the system of equations, we can conclude that it comes from the definition of the pressure in the solid phase, as the sum of the gas pressure plus intergranular stress of Equation (12). In fact, intergranular stress acts as a momentum-generation term when the initial solid phase velocity, 2 , is zero. This is precisely an effect sought in Interior Ballistics modeling, as it takes into account the compaction effects of the propellant: when the propellant grains are compressed by any cause in a region, their elasticity tends to release the tension, displacing contiguous propellant grains in the normal direction, similar to pressure. However, for the modeling of combustible cases, this effect is undesired, since the cases do not behave as a granular solid, but rather as a flexible solid wall, until it breaks or burns. Therefore, the intergranular stress term should be eliminated (Λ 2 = 0) to make pressures equal 2 = 1 , for the material of the case. This is not possible with the current system of governing Equations (1)-(9n), since there is a single momentum equation for all species in the solid phase (propellant, case, and initiator). Therefore, as an intermediate solution, it was considered here that 2 = 1 for the entire solid phase. In practice, this is equivalent to an incompressible solid phase, which is an approach used in other Interior Ballistics models. The implementation and testing of a dedicated subsystem of equations for each solid species, a task that would considerably increase the number of equations of the system, is left for future work.
All in all, we can conclude from the results (Figures 9-12) that by eliminating the intergranular stress term, for this specific case of MACS simulation, there is an improvement in the modeling of the initial moments of the combustion sequence, and the numerical result is less diffusive. To illustrate this, the firing of a six-module MACS in a 155/52 howitzer (chamber volume of 23 L) is simulated. The simplified geometry of Figure 6(a) is considered for the modules, whose interior is filled with 2.5 kg of extruded propellant in the shape of a 19-perforation rosette. For the case presented, a standard triple-base propellant composition M30 [26] was used, with a heat of explosion of 974 cal/g, specific heat ratio of 1.2385, covolume of 1.057·10 −3 m 3 /kg, and density of 1660 kg/m 3 . For the burning rate, modeled after Equation (26), we consider the coefficients = 4.597·10 −7 m/(s·Pa n ), = 0.652, and = 0 [26]. The results of this simulation are represented in Figure 9, where, employing contour plots of the gas volume fraction 1 , the evolution of solid concentration as the combustion process progresses is observed (considering that 2 = 1 − 1 ). In this figure (and also in Video S2 of Supplementary Materials, where the results are represented with a higher temporal resolution for 1 , pressure, temperature, and velocity), it is observed that, when the solid begins to burn through the area corresponding to the combustible case, 1 increases as the combustion source terms convert the solid phase Once the eigenvalues of the flux schemes of each phase are decoupled (i.e., for each phase, the respective eigenvalues are used), the result is clearly better, less diffusive, more realistic concerning the physics of the problem. The characteristic time corresponds now to the combustion process. However, a diffusion effect is still observed in the solid phase that, analyzing the system of equations, we can conclude that it comes from the definition of the pressure in the solid phase, as the sum of the gas pressure plus intergranular stress of Equation (12). In fact, intergranular stress acts as a momentum-generation term when the initial solid phase velocity, U 2 , is zero. This is precisely an effect sought in Interior Ballistics modeling, as it takes into account the compaction effects of the propellant: when the propellant grains are compressed by any cause in a region, their elasticity tends to release the tension, displacing contiguous propellant grains in the normal direction, similar to pressure. However, for the modeling of combustible cases, this effect is undesired, since the cases do not behave as a granular solid, but rather as a flexible solid wall, until it breaks or burns. Therefore, the intergranular stress term should be eliminated (Λ 2 = 0) to make pressures equal p 2 = p 1 , for the material of the case. This is not possible with the current system of governing Equations (1)-(9a), since there is a single momentum equation for all species in the solid phase (propellant, case, and initiator). Therefore, as an intermediate solution, it was considered here that p 2 = p 1 for the entire solid phase. In practice, this is equivalent to an incompressible solid phase, which is an approach used in other Interior Ballistics models. The implementation and testing of a dedicated subsystem of equations for each solid species, a task that would considerably increase the number of equations of the system, is left for future work.
All in all, we can conclude from the results (Figures 9-12) that by eliminating the intergranular stress term, for this specific case of MACS simulation, there is an improvement in the modeling of the initial moments of the combustion sequence, and the numerical result is less diffusive. To illustrate this, the firing of a six-module MACS in a 155/52 howitzer (chamber volume of 23 L) is simulated. The simplified geometry of Figure 6a is considered for the modules, whose interior is filled with 2.5 kg of extruded propellant in the shape of a 19-perforation rosette. For the case presented, a standard triple-base propellant composition M30 [26] was used, with a heat of explosion of 974 cal/g, specific heat ratio of 1.2385, covolume of 1.057 × 10 −3 m 3 /kg, and density of 1660 kg/m 3 . For the burning rate, modeled after Equation (26), we consider the coefficients a = 4.597 × 10 −7 m/(s·Pa n ), n = 0.652, and b = 0 [26]. The results of this simulation are represented in Figure 9, where, employing contour plots of the gas volume fraction α 1 , the evolution of solid concentration as the combustion process progresses is observed (considering that α 2 = 1 − α 1 ). In this figure (and also in Video S2 of Supplementary Materials, where the results are represented with a higher temporal resolution for α 1 , pressure, temperature, and velocity), it is observed that, when the solid begins to burn through the area corresponding to the combustible case, α 1 increases as the combustion source terms convert the solid phase into gas. In this process, temperature and pressure increase accordingly (see Video S2). The combustion process takes place on a time scale appropriate to the Internal Ballistics of a shot of the same characteristics. In fact, projectile motion starts at t ≈ 6 ms ( Figure 11) and is apparent to the naked eye (Figure 9 and Supplementary Video S2) at t ≈ 7 ms. The entire ballistic cycle, from the ignition to the exit of the projectile through the muzzle, lasts 16.25 s (in howitzers of this type, the shot can be between 15 and 20 ms [25]). The outcomes of 2D Axisymmetric and 3D simulations ( Figure 10) are practically identical, because the modules fill the whole diameter of the chamber and, therefore, in the 3D case, there is an almost concentric symmetry, like in the 2D Axisymmetric case.
Mathematics 2021, 9, x FOR PEER REVIEW 18 of 23 into gas. In this process, temperature and pressure increase accordingly (see Video S2).
The combustion process takes place on a time scale appropriate to the internal ballistics of a shot of the same characteristics. In fact, projectile motion starts at ≈ 6 ms ( Figure 11) and is apparent to the naked eye ( Figure 9 and Supplementary Video S2) at ≈ 7 ms. The entire ballistic cycle, from the ignition to the exit of the projectile through the muzzle, lasts 16.25 s (in howitzers of this type, the shot can be between 15 and 20 ms [25]). The outcomes of 2D Axisymmetric and 3D simulations ( Figure 10) are practically identical, because the modules fill the whole diameter of the chamber and, therefore, in the 3D case, there is an almost concentric symmetry, like in the 2D Axisymmetric case. The results for muzzle velocity and chamber pressure obtained with the last simulation are those represented in Figures 11 and 12. The pressure-time curve in Figure 12, with a maximum pressure of 349.25 MPa, corresponds to the expected curve in a modular charge shot, if the propellant used were an M30 composition (more energetic and faster than usual compositions of MACS) and the muzzle velocity, 940 m/s (i.e., the velocity of the projectile at the bore), corresponds to the nominal velocity of the system. into gas. In this process, temperature and pressure increase accordingly (see Video S2). The combustion process takes place on a time scale appropriate to the internal ballistics of a shot of the same characteristics. In fact, projectile motion starts at ≈ 6 ms ( Figure 11) and is apparent to the naked eye ( Figure 9 and Supplementary Video S2) at ≈ 7 ms. The entire ballistic cycle, from the ignition to the exit of the projectile through the muzzle, lasts 16.25 s (in howitzers of this type, the shot can be between 15 and 20 ms [25]). The outcomes of 2D Axisymmetric and 3D simulations ( Figure 10) are practically identical, because the modules fill the whole diameter of the chamber and, therefore, in the 3D case, there is an almost concentric symmetry, like in the 2D Axisymmetric case. The results for muzzle velocity and chamber pressure obtained with the last simulation are those represented in Figures 11 and 12. The pressure-time curve in Figure 12, with a maximum pressure of 349.25 MPa, corresponds to the expected curve in a modular charge shot, if the propellant used were an M30 composition (more energetic and faster than usual compositions of MACS) and the muzzle velocity, 940 m/s (i.e., the velocity of the projectile at the bore), corresponds to the nominal velocity of the system.

Conclusions
Relying on an Interior Ballistics model for numerical simulation is always a great help in the process of designing and characterizing propelling charges. If this mathematical and numerical model can predict three-dimensional effects in the combustion process, this would be an advantage, and even more so for the design of modular charges where

Conclusions
Relying on an Interior Ballistics model for numerical simulation is always a great help in the process of designing and characterizing propelling charges. If this mathematical and numerical model can predict three-dimensional effects in the combustion process, this would be an advantage, and even more so for the design of modular charges where  The results for muzzle velocity and chamber pressure obtained with the last simulation are those represented in Figures 11 and 12. The pressure-time curve in Figure 12, with a maximum pressure of 349.25 MPa, corresponds to the expected curve in a modular charge shot, if the propellant used were an M30 composition (more energetic and faster than usual compositions of MACS) and the muzzle velocity, 940 m/s (i.e., the velocity of the projectile at the bore), corresponds to the nominal velocity of the system.

Conclusions
Relying on an Interior Ballistics model for numerical simulation is always a great help in the process of designing and characterizing propelling charges. If this mathematical and numerical model can predict three-dimensional effects in the combustion process, this would be an advantage, and even more so for the design of modular charges where geometric effects are of utmost importance. To do this, the mathematical models underlying the simulation tool must be reliable, validated, and provide accurate results.
With this objective, a model of conservation equations that are capable of describing the transient 3D combustion behavior of propellants inside the chamber of an artillery gun was presented in this work. Both the governing equations and the closure laws were described. In addition, the numerical methods used in their resolution were presented. This mathematical model was implemented in the new software UXGun3D, which was first tested by simulating the firing of a conventional projection charge, obtaining very satisfactory results. The study shows it is feasible to simulate with a Eulerian-Eulerian model of conservation equations the problem of modular charge firing. The results obtained are in line with the expected nominal values. The main challenges are the presence of various energetic materials (propellant, initiator, combustible case) and, mainly, the different mechanical and chemical behavior of the combustible cartridge cases. The propellant responds well to this type of model; however, for the combustible cases, diffusion is observed that is not physical for this type of material and that comes from the generation of momentum from the intergranular stress term. To overcome this drawback, simulations were carried out, avoiding the intergranular stress term and obtaining satisfactory results. As for the numerical methods, the approach proposed, based on a formulation of the system of equations in terms of conserved variables, and the numerical flux schemes proved to be a good choice for the solution of this complex problem. However, a non-physical numerical diffusion was observed where the eigenvalues were coupled in the numerical flux schemes. It was found that decoupled eigenvalues are preferable, so that the eigenvalues of each phase subsystem of equations are used for their own numerical flux schemes. This way, the final results are fully satisfactory and provide the expected pressure and muzzle velocity values.
Future work will focus on implementing an even more detailed model, in which a dedicated momentum equation is available for each species of the solid phase, so that the intergranular stress term can be enabled for gunpowder, while being deactivated for the solid case. Likewise, the work may be continued by carrying out more experimental tests.