Comparative Study of Numerical Schemes for Granular Combustion of Boron Potassium Nitrate

: Multiple applications in aerospace utilize pyrotechnic charges for their operation, and these charges are predominantly in the form of granules. One of the most used charges is boron potassium nitrate (BPN), and the present study focuses on mathematically modeling granular combustion, its experimental recreation, and carrying out a comparative study on three upwind schemes for its numerical simulation. A customized version of the seven-equation compressible multifluid formulation is presented in this paper to model granular combustion mathematically. Three upwind schemes, namely HLLC, AUSM + -up, and HLLC-AUSM, are used for the numerical comparison. Utilizing these, an axisymmetric code is developed for the comparative study. To experimentally replicate granular combustion, granular BPN is fired in a closed bomb test facility, and the experimental pressure history is used for the numerical comparisons. The developed code can adequately predict the physics of granular combustion, and all three schemes are equally capable of numerical prediction.


Introduction
In the field of aerospace, numerous applications deploy pyrotechnic charges, namely pyrotechnic initiators, gas generators, airbags, ejection recovery systems, pyrotechnic fasteners, igniters in small solid rocket motors, decoy flares, ignition cartridges, etc. Metal oxidizers are commonly used compositions for these applications.Among compositions like zirconium potassium perchlorate, aluminum potassium perchlorate, and titanium aluminum potassium perchlorate, boron potassium nitrate (BPN) is the most widely used.BPN is thermally stable, suitable for rapid initiation, and stable in a vacuum, and its burn rates are almost independent of pressure.In all of these applications, BPN is used in granular form and is initiated by passing an electric current through a bridgewire.As the electric current passes, the bridgewire becomes heated and, in turn, heats the granules to their ignition temperature and produces hot combustion products.These hot combustion products flow into the remaining granular bed, heating and igniting it.During this process, a high-pressure gradient leads the ignition front to accelerate [1].The operation of all the applications mentioned earlier is crucially decided by the granular combustion of BPN, which brings importance to the present study.In granular combustion, there are two main phases to model, viz., a solid phase and a gaseous phase.Due to the generation of very high pressures during the operation, both phases are assumed to be compressible.A closed bomb test facility is utilized to experimentally replicate granular combustion.
In 1976, Krier and Ghokale [2] published a benchmark paper on granular combustion, where they explained how to predict the pressure wave propagation and flame spread within the porous propellant charge.Kuo et al. [3] developed a theoretical model that describes flame propagation in a packed bed of granular propellant under confinement.Following these two major works, a series of numerical studies on granular combustion

Closed Bomb Test
The phenomenon of granular combustion is experimentally replicated through a closed bomb test.The cylindrical test setup comprises a canister and chamber sections, as shown in Figure 1.The canister is 60 mm long with a 20 mm diameter, and the chamber is 150 mm long with a 60 mm diameter, as shown in Figure 2. Five grams of granular BPN is used for the test, and a sample is shown in Figure 3.The granular bed is loosely packed, placed in the canister section, and occupies 29 mm out of the 60 mm.Beneath the granular bed is a bridgewire through which an electric current is passed, and the bottom of the bed becomes ignited.Nichrome wire is used as the bridgewire, and a 30 A current passes through it.The pressure build-up is measured at 20 mm from the chamber head end through an absolute pressure sensor.The sensor specifications are mentioned in Table 1.The pressure sensor data are fed through a data acquisition module (model: NI9239) and recorded using the NI Signal Express ® software (version 2015) at 25 kHz.The recorded voltage readings are then corrected for residual unbalance and converted into corresponding pressure values.Through a particle size analysis, the diameter of the granular BPN used is determined to be 897 µm.The burn rate for granular BPN is experimentally determined to be 19.759P 0.1888 mm/s with pressure in bar.The density is measured and found to be 2000 kg/m 3 .Using the NASA CEA code [18], a constant volume (uv problem) chemical equilibrium analysis is carried out for 5 g of BPN, and the peak pressure predicted is 4.147 MPa with a flame temperature of 3011.5 K, and the combustion product properties are obtained.The experimental pressure-time curve for granular BPN and the CEA predicted peak pressure are shown in Figure 4.
are obtained.The experimental pressure-time curve for granular BPN and the CEA predicted peak pressure are shown in Figure 4.          are obtained.The experimental pressure-time curve for granular BPN and the CEA pre dicted peak pressure are shown in Figure 4.

Compressible Multifluid Formulation for Granular Combustion
In granular combustion, there are mainly two phases, i.e., a gaseous phase and a solid phase.The pressure generated during the operation of actual systems is very high.Hence, it is assumed that both phases are compressible and are governed by the corresponding equation of states, namely the ideal gas equation of state for the gaseous phase and the stiffened gas equation of state for the solid phase.To model the physics of granular combustion, the two-phase compressible multifluid formulation explained by Saurel et al. [9] is customized for granular combustion.The formulation is a full non-equilibrium model and is a modified version of the seven-equation model.The model consists of a set of balance equations, viz., mass, momentum, and energy for each phase, which are the first eight equations.The solid volume fraction evolution is the ninth equation.The first nine equations used for the modeling are explained in detail by Saurel et al. [9].The following equations are incorporated into these equations to curate the model for granular combustion, viz., the conservation of the number density of granules, the condensed phase combustion products formed, and the gas species equation to account for the initial air in the system and the gaseous combustion products created.The tenth equation is the evolution of number density per unit volume of granules [19].This equation ensures the number of granules is conserved by assuming that the granules burn to a small cutoff radius and remain at that cutoff radius.Using the tenth equation, a reduction in the granular radius is found at each time step.It is observed that there is a significant amount of condensed phase products formed during the combustion of BPN when the chemical equilibrium analysis is carried out.The eleventh equation represents the conservation of condensed phase combustion product volume fraction [9].It is assumed that the condensed and gaseous phases are in equilibrium.The gas species transport equation is the twelfth equation [17], as most of the applications operate in atmospheric conditions and contain air.By including the species transport equation, the availability of air from the beginning, its heating and the later addition of gaseous combustion products could be differentiated.All of the equations are in conservation form, which makes the governing equations fully hyperbolic and a well-posed problem.
Due to the axisymmetric nature of the closed bomb test vessel, the formulation is axisymmetric, and it is written in generic form, where U represents the conserved variables, F and G are the flux terms in axial and radial directions, respectively, and S is the source term.Suffixes g, p, and I denote the gas phase, solid phase, and interface, respectively.∅, , u, v, E, P, N, and Y represent the volume fraction, density, velocity in the axial direction, velocity in the radial direction, total energy, pressure, number density per unit volume, and local mass fraction, respectively.

Compressible Multifluid Formulation for Granular Combustion
In granular combustion, there are mainly two phases, i.e., a gaseous phase and a solid phase.The pressure generated during the operation of actual systems is very high.Hence, it is assumed that both phases are compressible and are governed by the corresponding equation of states, namely the ideal gas equation of state for the gaseous phase and the stiffened gas equation of state for the solid phase.To model the physics of granular combustion, the two-phase compressible multifluid formulation explained by Saurel et al. [9] is customized for granular combustion.The formulation is a full non-equilibrium model and is a modified version of the seven-equation model.The model consists of a set of balance equations, viz., mass, momentum, and energy for each phase, which are the first eight equations.The solid volume fraction evolution is the ninth equation.The first nine equations used for the modeling are explained in detail by Saurel et al. [9].The following equations are incorporated into these equations to curate the model for granular combustion, viz., the conservation of the number density of granules, the condensed phase combustion products formed, and the gas species equation to account for the initial air in the system and the gaseous combustion products created.The tenth equation is the evolution of number density per unit volume of granules [19].This equation ensures the number of granules is conserved by assuming that the granules burn to a small cutoff radius and remain at that cutoff radius.Using the tenth equation, a reduction in the granular radius is found at each time step.It is observed that there is a significant amount of condensed phase products formed during the combustion of BPN when the chemical equilibrium analysis is carried out.The eleventh equation represents the conservation of condensed phase combustion product volume fraction [9].It is assumed that the condensed and gaseous phases are in equilibrium.The gas species transport equation is the twelfth equation [17], as most of the applications operate in atmospheric conditions and contain air.By including the species transport equation, the availability of air from the beginning, its heating and the later addition of gaseous combustion products could be differentiated.All of the equations are in conservation form, which makes the governing equations fully hyperbolic and a well-posed problem.
Due to the axisymmetric nature of the closed bomb test vessel, the formulation is axisymmetric, and it is written in generic form, where U represents the conserved variables, F and G are the flux terms in axial and radial directions, respectively, and S is the source term.Suffixes g, p, and I denote the gas phase, solid phase, and interface, respectively.H, ρ, u, v, E, P, N, and Y represent the volume fraction, density, velocity in the axial direction, velocity in the radial direction, total energy, pressure, number density per unit volume, and local mass fraction, respectively. with, Br `∅g The basic assumption is that at any instant in time, all phases are present at every cell in the domain.This requires constitutive relations to close the system of governing equations.
The total volume fraction is constrained by Equation (3).
∅ g `∅p `∅q " 1 Each phase is compressible and is governed by the corresponding equation of state.The gaseous phase is governed by the ideal gas equation of state (Equation ( 4)), and the solid phase is governed by the stiffened gas equation of state (Equation ( 5)) [9].
The total energy for each phase [13] is given by Equations ( 6) and (7).
The granules are assumed to be spherical, and the number density per unit volume [1] is given by Equation (8).
Aerospace 2024, 11, 251 6 of 16 The mass generated during the combustion of the granular bed [1] is given by Equation ( 9), and the burn rate is shown through Equation ( 10). .
. r " aP g n (10) Equations ( 11) and ( 12) compute the speed of sound in the gas and solid phases [13].
P I " ρ g a g P p `ρp a p P g ρ g a g `ρp a p U I " ρ g a g u g `ρp a p u p ρ g a g `ρp a p (14) The drag effects are accounted for by the correlation given by Saurel et al. [9].In the solid volume fraction evolution equation, the pressure of each phase relaxes each other to a common equilibrium at a rate controlled by the pressure relaxation parameter [9] in Equation (16).
where the pressure relaxation time [9] is given by Equation (17).

Numerical Schemes
To numerically predict granular combustion, an unstructured finite volume method is used for spatial discretization with second-order accurate upwind schemes for flux computation, and an explicit strong stability-preserving Runge-Kutta method in four stages with third-order accuracy (eSSPRK(4,3) by Isherwood et al. [20]) is used for temporal discretization.
Spatial discretization is used to calculate the convective fluxes and the source terms.The finite volume method, specifically with cell-centered scheme, is used.The governing equations in generic form are written as Integrating in space over the entire domain Ω (Blazek [21]): The control volume does not change with time; upon the integration of the first term, we obtain B Bt The second and third terms are the convective and source terms, respectively.By applying the Gauss divergence theorem to the second and third terms, we obtain which could be written as Among the various methods used to compute the fluxes at the cell face, upwind schemes, namely HLLC, AUSM+-up, and HLLC-AUSM, are used for the present study.The upwind schemes are explained in the following sub-sections and are used for the comparative study.All of the upwind schemes are second-order accurate, and the Venkatakrishnan limiter reported by Michalak and Ollivier-Gooch [22] is used.With the described mathematical and numerical model, an axisymmetric code is developed for the numerical study and for the prediction of granular combustion.

HLLC
A simple HLLC scheme for the compressible multifluid formulation for granular combustion is derived from the HLLC scheme mentioned by Batten et al. [23].The HLLC scheme resolves both the shock waves and contact discontinuities.It involves seven waves, i.e., the three conventional left-and right-facing waves for each phase and a contact wave.Here, each phase has three waves, i.e., left-and right-facing waves (S L , S R ), with a contact wave (S M ) separating into four states for each phase.This scheme is easy to code and is very robust.
' % E l,g ; S L,g ą 0 E l,g `SL,g pU * l,g ´Ul,g q ; S L,g ď 0 ă S M,g E r,g `SR,g pU * r,g ´Ur,g q ; S M,g ď 0 ď S R,g E r,g ; S R,g ă 0 E l,p ; S L,p ą 0 E l,p `SL,p pU * l,p ´Ul,p q ; S L,p ď 0 ă S M,p E r,p `SR,p pU * r,p ´Ur,p q ; S M,p ď 0 ď S R,p E r,p ; S R,p ă 0 (23) where S L,g " min ´ql,g ´al,g , q r,g ´ar,g ¯; S L,p " min ´ql,p ´al,p , q r,p ´ar,p ¯(24) S R,g " max ´qr,g `ar,g , q l,g `al,g ¯; S R,p " max ´qr,p `ar,p , q l,p `al,p ¯(25) S M,g " ρ r,g q r,g `SR,g ´qr,g ˘´ρ l,g q l,g ´SL,g ´ql,g ¯`P l,g ´Pr,g ρ r,g `SR,g ´qr,g ˘´ρ l,g ´SL,g ´ql,g ¯ S M,p " ρ r,p q r,p `SR,p ´qr,p ˘´ρ l,p q l,p ´SL,p ´ql,p ¯`P l,p ´Pr,p ρ r,p `SR,p ´qr,p ˘´ρ l,p ´SL,p ´ql,p

¯(27)
∅ l,g P g " ∅ l,g P l,g `∅l,g ρ l,g pS L,g ´ql,g qpS M,g ´ql,g q; ∅ l,p P p " ∅ l,p P l,p `∅l,p ρ l,p pS L,p ´ql,p qpS M,p ´ql,p q (28) Denoting k = l (left), r (right): q k,g " u k,g n x `vk,g n r ; q k,p " u k,p n x `vk,p n r,p ∅ k,g ρ k,g q k,g ∅ k,g ρ k,g u k,g q k,g `∅k,g P k,g n x ∅ k,g ρ k,g v k,g q k,g `∅k,g P k,g n r ∅ k,g ρ k,g E k,g q k,g `∅k,g P k,g q k,g ∅ k,p ρ k,p q k,p ∅ k,p ρ k,p u k,p q k,p `∅k,p P k,p n x ∅ k,p ρ k,p v k,p q k,p `∅k,p P k,p n r ∅ k,p ρ k,p E k,p q k,p `∅k,p P k,p q k,p ∅ k,p q k,p N k,p q k,p ,g P g ´∅k,g P k,g ¯nx ´∅k ,g P g ´∅k,g P k,g ¯nr ∅ k,g P g S M,g ´∅k,g P k,g q k,g 0 0

AUSM + -up
The Advection Upstream Splitting Method (AUSM) scheme, which is used to compute the fluxes, is known to give robust and accurate solutions for flows from low Mach numbers to hypersonic regimes.Liou explains the detailed evolution of AUSM schemes in [24].The AUSM + -up scheme is an extension of the AUSM family to solve flows at all speed regimes.Kitamura et al. [15] conducted a detailed survey of AUSM schemes for multifluid flows.
From Equation ( 22), we can obtain The flux, E, is split into convective flux and pressure flux: where .
m " " ρ g q g i f gas phase ρ p q p i f solid phase (37) q g " u g n x `vg n r ; q p " u p n x `vp n r (38) The flux at each face is computed using the following expression, where k stands for the gas and solid phases: where N " p0, n x , n r , 0, 0, n x , n r , 0, 0, 0, 0, 0q T The mass flux ( .m) and pressure flux ( " P) expressions for the AUSM + -up scheme are directly incorporated from the study by Kitamura et al. [15].

HLLC-AUSM
In this scheme, the gas phase flux is computed using the HLLC scheme, and the solid phase is computed using the AUSM + -up scheme.The flux computation using the HLLC-AUSM scheme is as follows.Houim and Oran have used the HLLC-AUSM scheme in most of their works [16,17].
Gas phase flux computation: Solid phase flux computation:

Computational Framework
The computational domain of the closed bomb test with the canister and chamber is shown in Figure 5. Five grams of BPN granules is filled to 29 mm in the canister with a 0.275 initial particle volume fraction, and with the assumption that the initial 5 mm is ignited, the problem is numerically simulated.The input parameters are tabulated in Table 2.The domain is filled with quadrilateral grids.It is assumed that the granules are spherical in shape.In Figure 5, the boundary marked with orange is assigned the symmetry boundary condition, and those in blue are assigned the wall boundary condition.
Solid phase flux computation:

Computational Framework
The computational domain of the closed bomb test with the canister and chamber is shown in Figure 5. Five grams of BPN granules is filled to 29 mm in the canister with a 0.275 initial particle volume fraction, and with the assumption that the initial 5 mm is ignited, the problem is numerically simulated.The input parameters are tabulated in Table 2.The domain is filled with quadrilateral grids.It is assumed that the granules are spherical in shape.In Figure 5, the boundary marked with orange is assigned the symmetry boundary condition, and those in blue are assigned the wall boundary condition.the initial ignition length assumption, the code predicts the pressure history well.Two sets of grids are taken, one with six thousand six hundred elements and another with thirteen thousand two hundred elements, for grid sensitivity studies with all three schemes, as shown in Figure 8. From this, the rest of the studies assume that 5 mm of the granular bed is ignited, and flux calculations are carried out using the HLLC scheme.
with all three schemes and compared with the experimental results, as shown in Figure 6.The numerical prediction carried out by all of the schemes complements the experimental result.The simulation is terminated when the radius of the granules reaches the cutoff value, and the pressure remains constant after that, as heat loss to the surrounding area is not considered.A sensitivity study on the assumed initial ignited length is carried out.Simulations when the initial ignited length is 2.5 mm, 5 mm, and 7.5 mm are carried out with all of the schemes.The results are shown in Figure 7.It is observed that, irrespective of the initial ignition length assumption, the code predicts the pressure history well.Two sets of grids are taken, one with six thousand six hundred elements and another with thirteen thousand two hundred elements, for grid sensitivity studies with all three schemes, as shown in Figure 8. From this, the rest of the studies assume that 5 mm of the granular bed is ignited, and flux calculations are carried out using the HLLC scheme.(c) From the beginning to the complete combustion of the granules, the centerline plot for the evolution of solid volume fraction and the reduction in the granular radius are shown in Figures 9 and 10, respectively.From Figures 9a and 10a, the spreading of the granules to the entire domain from the initial 29 mm bed can be observed.The clear From the beginning to the complete combustion of the granules, the centerline plot for the evolution of solid volume fraction and the reduction in the granular radius are shown in Figures 9 and 10, respectively.From Figures 9a and 10a, the spreading of the granules to the entire domain from the initial 29 mm bed can be observed.The clear reduction in the solid volume fraction in the initial granular bed can be seen in Figure 9a.By 650 µs, the pressure in the chamber has risen from 1 bar to 5 bar, and from Figure 10a, it is visible that the radius of the granules has not reduced much, only from 448.5 µm to 423.4 µm.A significant reduction in the radius occurs after 650 µs until complete combustion at 10 ms as the pressure in the chamber rises from 5 bar to 47 bar.This is evident from Figure 10b, and at 10 ms, the radius of the granules has reached the cutoff radius, which is 5 µm throughout the domain.
By 650 µs, the pressure in the chamber has risen from 1 bar to 5 bar, and from Figure 10a, it is visible that the radius of the granules has not reduced much, only from 448.5 µm to 423.4 µm.A significant reduction in the radius occurs after 650 µs until complete combustion at 10 ms as the pressure in the chamber rises from 5 bar to 47 bar.This is evident from Figure 10b, and at 10 ms, the radius of the granules has reached the cutoff radius, which is 5 µm throughout the domain.

Conclusions
This study focuses on experimentally recreating granular combustion through a closed bomb test, mathematically modeling the complete physics of granular combustion, numerically predicting the same, and comparing three numerical schemes.Boron potassium nitrate is used for this study, as it is widely used in various applications.A customized version of the seven-equation compressible multifluid formulation is curated and presented in this paper and is used for modeling.The finite volume method is used for spatial discretization, and three upwind schemes from the literature, namely HLLC, AUSM+-up, and a combination of the two families of schemes, referred to as HLLC- it is visible that the radius of the granules has not reduced much, only from 448.5 µm to 423.4 µm.A significant reduction in the radius occurs after 650 µs until complete combustion at 10 ms as the pressure in the chamber rises from 5 bar to 47 bar.This is evident from Figure 10b, and at 10 ms, the radius of the granules has reached the cutoff radius, which is 5 µm throughout the domain.

Conclusions
This study focuses on experimentally recreating granular combustion through a closed bomb test, mathematically modeling the complete physics of granular combustion, numerically predicting the same, and comparing three numerical schemes.Boron potassium nitrate is used for this study, as it is widely used in various applications.A customized version of the seven-equation compressible multifluid formulation is curated and presented in this paper and is used for modeling.The finite volume method is used for spatial discretization, and three upwind schemes from the literature, namely HLLC, AUSM+-up, and a combination of the two families of schemes, referred to as HLLC-

Conclusions
This study focuses on experimentally recreating granular combustion through a closed bomb test, mathematically modeling the complete physics of granular combustion, numerically predicting the same, and comparing three numerical schemes.Boron potassium nitrate is used for this study, as it is widely used in various applications.A customized version of the seven-equation compressible multifluid formulation is curated and presented in this paper and is used for modeling.The finite volume method is used for spatial discretization, and three upwind schemes from the literature, namely HLLC, AUSM+-up, and a combination of the two families of schemes, referred to as HLLC-AUSM, are used for the comparative study.With these, an axisymmetric code is developed to predict and compare granular combustion.The developed code could predict granular combustion, and it complements the experimental result.The overall computation time for the HLLC scheme is slightly longer compared to the AUSM+-up and HLLC-AUSM schemes; however, there is not much to distinguish between them.All three schemes can predict granular combustion accurately.In future work, the chemical kinetics of granular BPN can be incorporated into the model for enhanced numerical predictions.

Figure 4 .
Figure 4. Experimental pressure history and CEA predicted peak pressure.

Figure 4 .
Figure 4. Experimental pressure history and CEA predicted peak pressure.

Figure 6 .
Figure 6.Pressure history numerical and experimental comparison.
ParametersValue Make and model number Kulite; EWCTV 312(M) Pressure range 140 bar Operating temperature range 297 to 1366 K Residual unbalance 0.5 V ± 100 mV Uncertainty in pressure measurement ±0.1% of full scale (FS)
ParametersValue Make and model number Kulite; EWCTV 312(M) Pressure range 140 bar Operating temperature range 297 to 1366 K Residual unbalance 0.5 V ± 100 mV Uncertainty in pressure measurement ±0.1% of full scale (FS)
5 V ˘100 mV Uncertainty in pressure measurement ˘0.1% of full scale (FS) ´∅k ,p P p ´∅k,p P k,p ¯nx ´∅k ,p P p ´∅k,p P k,p ¯nr ∅ k,p P p S M,p ´∅k,p P k,p q k,p

Table 2 .
Input parameters for simulation.