Scaling Turbulent Combustion Fields in Explosions

: We considered the topic of explosions from spherical high-explosive (HE) charges. We studied how the turbulent combustion ﬁelds scale. On the basis of theories of dimensional analysis by Bridgman and similarity theories of Sedov and Barenblatt, we found that all ﬁelds scaled with the explosion length scale r 0 . This included the blast wave, the mean and root mean squared ( RMS ) proﬁles of thermodynamic variables, combustion variables, velocities, vorticity, and turbulent Reynolds stresses. This was a consequence of the formulation of the problem and our numerical method, which both satisﬁed the similarity conditions of Sedov. We performed numerical simulations of 1 g charges and 1 kg charges; the solutions were identical (within roundo ﬀ error) when plotted in scaled variables. We also explored scaling laws related to three-phase pyrotechnic explosions. We show that although the scaling formally broke down, the ﬁreball still essentially scaled with the explosion length scale r 0 . However, the discrete Lagrange particles (DLP) (phase 2) and the heterogeneous continuum model (HCM) of the DLP wakes (phase 3) did not scale with r 0 , and mean and RMS proﬁles could di ﬀ er by a factor of 10 in some regions. This was because the DLP particles and wakes introduced an additional scale that broke the similarity conditions.


Introduction
Similarity theory dates back to Bridgman's work [1] summarized in Dimensional Analysis, which he published in 1922. He points out that the differential equations describing physical phenomena can be rescaled into dimensionless units, resulting in equations that are invariant over dimensionless groups. Solutions are then invariant in settings where the dimensionless groups remain constant. This is codified in the "Buckingham Π Theorem" [2]. Sedov [3] applied such concepts to mechanics in Similarity and Dimensional Methods in Mechanics (1959). His work codified and summarized similarity methods in fluid mechanics. The similarity theory was used to derive the similarity solutions of the point explosion problem first published by G. I. Taylor [4] in 1941 and later by Sedov [5] in 1946. This same similarity solution method was used by Oppenheim and Kuhl to compute all possible similarity solutions bounded by a strong shock [6], a strong detonation wave [7], and self-similar explosion waves of variable energy at the front [8]. Scaling aspects of blast waves were summarized by G. I. Barenblatt in his book Scaling, Self-Similarity, and Intermediate Asymptotics [9] published in 1996.
Here, we consider the topic of explosions from spherical high-explosive (HE) charges. We are interested in the turbulent combustion fields in their fireballs, as illustrated in Figure 1. In this work, we derive scaling laws for all fields (including the thermodynamic variables, velocities, component mass-fractions, and turbulent Reynolds stresses) on the basis of the similarity theory for explosion fields, as described by Sedov.

Relevance of This Work
One of the reasons that such scaling laws are important is because we can study the characteristics of turbulent fireballs at the laboratory scale, such as 50 g charge explosions, as reported by Glumac and Kuhl [10], then scale them up to large scale explosions, such as the 15 kg charge explosion shown in Figure 1. Many more diagnostic techniques can be applied at the laboratory scale but are impractical at the field-test scale. Such laboratory measurements include schlieren/shadow photography; holographic interferometry of particles; spectral imaging in the IR, visible, and UV wavelengths; measurements of species composition via gas chronometry; collection of carbon particles that are the source of strong Planckian radiation [10]; and optical measurements inside the fireball, as demonstrated in [10], to name just a few. In other words, using the scaling laws demonstrated in this paper, we can measure intimate details of turbulent fireballs in laboratory experiments, then confidently apply them to large-scale explosion events like that shown in Figure 1.
In contrast, the particle and mist phases of pyrotechnic explosions do not obey Sedov's scaling laws. Thus, calculations and experiments must be performed at each size of pyrotechnic charge under consideration. Even if the initial geometric charge size is scaled, the evolving flow fields of the particles and mist do not scale. Such facts have not been recognized in the past.

Gas Dynamics of Explosions
The turbulent combustion fields are governed by the inviscid conservation laws of gas dynamics. These may be expressed in strong conservation form as Mass: Momentum: Energy: Component DP:

Relevance of This Work
One of the reasons that such scaling laws are important is because we can study the characteristics of turbulent fireballs at the laboratory scale, such as 50 g charge explosions, as reported by Glumac and Kuhl [10], then scale them up to large scale explosions, such as the 15 kg charge explosion shown in Figure 1. Many more diagnostic techniques can be applied at the laboratory scale but are impractical at the field-test scale. Such laboratory measurements include schlieren/shadow photography; holographic interferometry of particles; spectral imaging in the IR, visible, and UV wavelengths; measurements of species composition via gas chronometry; collection of carbon particles that are the source of strong Planckian radiation [10]; and optical measurements inside the fireball, as demonstrated in [10], to name just a few. In other words, using the scaling laws demonstrated in this paper, we can measure intimate details of turbulent fireballs in laboratory experiments, then confidently apply them to large-scale explosion events like that shown in Figure 1.
In contrast, the particle and mist phases of pyrotechnic explosions do not obey Sedov's scaling laws. Thus, calculations and experiments must be performed at each size of pyrotechnic charge under consideration. Even if the initial geometric charge size is scaled, the evolving flow fields of the particles and mist do not scale. Such facts have not been recognized in the past.

Gas Dynamics of Explosions
The turbulent combustion fields are governed by the inviscid conservation laws of gas dynamics. These may be expressed in strong conservation form as Mass: Energy: Component DP: In the equations above, ρ is density, u represents velocity, E denotes total energy (E = e + 0.5 u·u), and Y D is the mass fraction of detonation product species. Two components are recognized: detonation products (DP), which serve as the fuel, and air (A). Their mass fractions obey the conservation relation: In pure computational cells, we have either Y D = 1 or Y A = 1; in mixed cells, Y D and Y A < 1. We assume that the components are in a state of thermodynamic equilibrium. This allows us to compute the pressure and temperature equation of states (EOS) according to Pressure: Temperature: These functions are evaluated by using the thermodynamic equilibrium code Cheetah [11,12]. Results have been tabulated for computational efficiency.
We follow the evolution of component mass fractions of fuel (detonation products): Y F , oxidizer; (air): Y O ; and products (combustion products): Y P . They were computed according to the following conservation laws: Conservation: where α denotes the stoichiometric oxidizer/fuel ratio. We use the fast chemistry, infinite Damkohler number limit (i.e., gas dynamic limit), where turbulent mixing within the cell occurs according to the MILES concept of Boris [13]: Combustion Rate: . ϕ = burn all Air : ∆ A / ∆t (Fuel − rich cell) burn all Fuel : ∆ F / ∆t (Air − rich cell) The problem is initialized by a similarity solution of Kuhl [14] for a constant velocity detonation wave propagating at the Chapman-Jouguet (CJ) speed when the detonation wave reaches the surface of the charge. See Figure 2 for the similarity profiles. (We used a perturbation on the density profile to break nonphysical numerical symmetries. We discuss HE density perturbation effects on the mean and root mean squared (RMS) profiles in the Section 6.) The detonation products expand, doing p·dV work on the surroundings, thereby creating a blast wave. The DP-air interface is unstable; perturbations grow due to vorticity creation by the inviscid baroclinic mechanism: . ω = 1 ρ 2 ∇ρ × ∇p, thereby creating a turbulent combustion layer near the edge of the fireball.
We solved the above conservation laws with an unsplit second-order Godunov method. We used the unsplit piecewise parabolic method (PPM) of Colella and Woodward [16,17] to advance the solution in time. PPM flattens slopes in cells near discontinuities to enforce monotonicity constraints (i.e., suppression of oscillations near discontinuities). This slope flattening induces a non-linear dissipation mechanism that acts on the cell level. It also reduces the scheme from second-order accurate Appl. Sci. 2020, 10, 8577 4 of 28 in smooth regions of the flow, to first order near discontinuities. This is consistent with the MILES (monotone integrated large eddy simulations) concept of Jay Boris [13], whereby mixing/dissipation only occur on the smallest grid scale. We used adaptive mesh refinement (AMR) [18,19] to follow steep gradients and mixing structures.
Appl. Sci. 2020, 10, x 4 of 30 Figure 2. Initial conditions based on the similarity solutions for a Chapman-Jouguet (CJ) detonation wave in the gas phase by Kuhl [14], and the linear velocity profile of an inertial flow expansion into a vacuum by Stanyukovich [15] for the droplet phase.
The detonation products expand, doing • work on the surroundings, thereby creating a blast wave. The DP-air interface is unstable; perturbations grow due to vorticity creation by the inviscid baroclinic mechanism: = × , thereby creating a turbulent combustion layer near the edge of the fireball. We solved the above conservation laws with an unsplit second-order Godunov method. We used the unsplit piecewise parabolic method (PPM) of Colella and Woodward [16,17] to advance the solution in time. PPM flattens slopes in cells near discontinuities to enforce monotonicity constraints (i.e., suppression of oscillations near discontinuities). This slope flattening induces a non-linear dissipation mechanism that acts on the cell level. It also reduces the scheme from second-order accurate in smooth regions of the flow, to first order near discontinuities. This is consistent with the MILES (monotone integrated large eddy simulations) concept of Jay Boris [13], whereby mixing/dissipation only occur on the smallest grid scale. We used adaptive mesh refinement (AMR) [18,19] to follow steep gradients and mixing structures.

Scaling Laws
In the inviscid gas dynamics formulation described above, variables such as density, velocity, pressure, and internal energy are related in ways controlled by their dimensional units. Sedov [3] and Barenblatt [9] show that chemical explosion fields are controlled by the non-dimensional independent variables of space, , and time, , according to where the explosion length scale is and the explosion time scale is = / (12) Figure 2. Initial conditions based on the similarity solutions for a Chapman-Jouguet (CJ) detonation wave in the gas phase by Kuhl [14], and the linear velocity profile of an inertial flow expansion into a vacuum by Stanyukovich [15] for the droplet phase.

Scaling Laws
In the inviscid gas dynamics formulation described above, variables such as density, velocity, pressure, and internal energy are related in ways controlled by their dimensional units. Sedov [3] and Barenblatt [9] show that chemical explosion fields are controlled by the non-dimensional independent variables of space, x, and time, τ, according to where the explosion length scale is and the explosion time scale is In the above, m is the mass of the charge, ∆H d represents the heat of detonation of the explosive material, and a a denotes the sound speed of the ambient atmosphere. The variable j equals 0, 1, or 2 for planarly, cylindrically, or spherically symmetric flows. Similitude theory states [3] that the non-dimensional dependent variables of U = , u, e, p, T of gas dynamics are functions of x and τ: For spherical explosions in three-dimensional Cartesian co-ordinates, Equation (13) becomes x = x/r 0 and y = y/r 0 and z = z/r 0 (22) and the spherical (j = 2) explosion length scale becomes The dependent variables then become The above is the form of the solution of the inviscid conservation laws of gas dynamics for spherical explosions. Everything scales with the explosion length scale, r 0 . Figure 3 depicts a cross-section of the temperature field of a fireball from a 1 g TNT (Trinitrotoluene explosive: C 6 H 2 (NO 2 ) 3 CH 3 ) explosion in air at 100 µs/g 1/3 . One can see a variety of scales in the mixing layer. Combustion in this inviscid formulation occurs in an exothermic sheet. Our numerical method satisfies same type of scaling laws as the fireball. Consequently, we expect the numerics to also be self-similar. We confirmed this observation by simulations at two scales.     Figure 4 compares the density fields; it shows the outer shock (blast wave) and secondary expanding shock, as well as the turbulent mixing in the fireball. Figure 5a,b shows comparisons of the vorticity fields in the fireballs, showing mean fireball radii of R FB = 30 cm/g 1/3 and R FB = 300cm/kg 1/3 , respectively. These are driven by the gas dynamic fields, which scale with the explosion length scale, r 0 .

Fireball Cross-Sections
(a) Temperature cross-section (b) Blow-up near the front         It is difficult to decern differences in the two flow fields from the above cross-section figures. Thus, we computed the maximum density difference between the two cases at each scaled time. Results are shown in Figure 6. Differences started off as roundoff errors O ∼ (10 −13 ). They grew exponentially at around 5 µs/g 1/3 , and eventually saturated at values of tens of percent.

Mean and RMS Profiles
To quantitatively compare the flow field differences in Figures 4 and 5, we azimuthally averaged the three dimensions fields in θ and φ to produce mean and RMS (root mean squared) radial profiles of the turbulent variables. The averaging method is presented in Appendix A.
It is difficult to decern differences in the two flow fields from the above cross-section figures. Thus, we computed the maximum density difference between the two cases at each scaled time. Results are shown in Figure 6. Differences started off as roundoff errors ~(10 ). They grew exponentially at around 5 μs/g / , and eventually saturated at values of tens of percent.

Mean and RMS Profiles
To quantitatively compare the flow field differences in Figures 4 and 5, we azimuthally averaged the three dimensions fields in to produce mean and RMS (root mean squared) radial profiles of the turbulent variables. The averaging method is presented in Appendix A.
Figures 7 and 8 depict azimuthally averaged profiles of the thermodynamic profiles of temperature, pressure, and density from a 1 g charge at = 0.5015 (blue curves) compared with flow fields from the 1 kg charge at the same scaled time, = 5.015 (orange curves). The mean and RMS profiles from the two scales overlayed.
(a) mean temperature: ( b) mean pressure: ̅ Figure 6. Maximum density difference between the 1 g case and the 1 kg case as a function of time.
Differences stayed at roundoff errors until 5 µs, when they grew exponentially at late times.
Figures 7 and 8 depict azimuthally averaged profiles of the thermodynamic profiles of temperature, pressure, and density from a 1 g charge at t = 0.5015 ms (blue curves) compared with flow fields from the 1 kg charge at the same scaled time, t = 5.015 ms (orange curves). The mean and RMS profiles from the two scales overlayed.
It is difficult to decern differences in the two flow fields from the above cross-section figures. Thus, we computed the maximum density difference between the two cases at each scaled time. Results are shown in Figure 6. Differences started off as roundoff errors ~(10 ). They grew exponentially at around 5 μs/g / , and eventually saturated at values of tens of percent.

Mean and RMS Profiles
To quantitatively compare the flow field differences in Figures 4 and 5, we azimuthally averaged the three dimensions fields in to produce mean and RMS (root mean squared) radial profiles of the turbulent variables. The averaging method is presented in Appendix A.         Figure 13 shows azimuthally averaged off-diagonal turbulent Reynolds stresses ( , , and ) from a 1 g charge at = 0.5015 ms, compared with flow fields from the 1 kg charge at the same scaled time, = 5.015 ms. The profiles were essentially zero, with random fluctuations due to the small ensemble size at small radii.  Figure 14 shows the maximum density histories for a 1 g charge (blue curve) and a 1 kg charge (orange curve); the curves overlayed. Maximum densities decayed because of the adiabatic expansion of the detonation products, with a spike at = 0.1 ms/g / due to the imploding shock [20].  Figure 13 shows azimuthally averaged off-diagonal turbulent Reynolds stresses (τ rθ , τ rφ , and τ θφ ) from a 1 g charge at t = 0.5015 ms, compared with flow fields from the 1 kg charge at the same scaled time, t = 5.015 ms. The profiles were essentially zero, with random fluctuations due to the small ensemble size at small radii.    Figure 14 shows the maximum density histories for a 1 g charge (blue curve) and a 1 kg charge (orange curve); the curves overlayed. Maximum densities decayed because of the adiabatic expansion of the detonation products, with a spike at = 0.1 ms/g / due to the imploding shock [20].  Figure 14 shows the maximum density histories for a 1 g charge (blue curve) and a 1 kg charge (orange curve); the curves overlayed. Maximum densities decayed because of the adiabatic expansion of the detonation products, with a spike at t = 0.1 ms/g 1/3 due to the imploding shock [20].
ppl. Sci. 2020, 10, x 12 o Figure 14. Maximum density histories for 1 g and 1 kg charges. Figure 15 presents the combustion histories for 1 g and (blue curve) and scaled 1 kg char orange curve); the curves overlayed. About 50% of the fuel mass was consumed by combustion his time.   Figure 15 presents the combustion histories for 1 g and (blue curve) and scaled 1 kg charge (orange curve); the curves overlayed. About 50% of the fuel mass was consumed by combustion at this time. Let us azimuthally average the fireball temperature field in , to find the evolution of the mean temperature field with time: Let us define the fireball radius as the radius where ( , ) = /2 (starting from outside the fireball and working inward). Then, the evolution of the fireball radius becomes Let us azimuthally average the fireball temperature field in θ, φ to find the evolution of the mean temperature field with time: Let us define the fireball radius as the radius where T FB (x FB , τ) = T max /2 (starting from outside the fireball and working inward). Then, the evolution of the fireball radius becomes Dividing by the initial nondimensional charge radius x c = r c /r 0 , one finds From experiments and numerical simulations, when the spherical fireball has expanded to one atmosphere, one obtains The largest characteristic scale in the turbulent fireball is the fireball radius, and thus turbules in the fireball scale with r FB .

Summary
We confirmed the similarity of turbulent combustion fields of fireballs from HE explosions. Figure 16 presents turbulent fireball structures from our gadynamic simulations (Figure 16a,c), which are compared with a photograph of a fireball (Figure 16b). They show similar turbulent structures that scaled with the explosion length scale, r 0 .

Summary
We confirmed the similarity of turbulent combustion fields of fireballs from HE explosions. Figure 16 presents turbulent fireball structures from our gadynamic simulations (Figure 16a,c), which are compared with a photograph of a fireball (Figure 16b). They show similar turbulent structures that scaled with the explosion length scale, r .  As expected, all gasdynamic profiles (e.g., thermodynamic profiles, combustion profiles, velocity profiles, and Reynolds stresses), both their mean and RMS values, scaled with the explosion scale r 0 = (m·∆H d /p a ) 1/3 . The fireball radius scaled according to R FB /r 0 = f (t/t 0 ). To sum up, all aspects of gaseous turbulent combustion in HE explosions scaled with the explosion radius r 0 .

Three-Phase Pyrotechnic Explosions
Next, we considered pyrotechnic explosions. Figure 17 presents a cross-section of the three-phase temperature field of a pyrotechnic explosion from our paper in the 16th International Detonation Symposium [21]. It shows the turbulent mixing and combustion of the gas phase, which reached a temperature of 2500 (yellow regions) corresponding to the adiabatic flame temperature of the HE-air system. It also shows discrete Lagrange particles (black dots) that shed micro-mist wakes (red curves). This three-phase nature of such pyrotechnic explosions is what we explore in this section. Models are described next.

Gas Phase
The conservation laws of gas dynamics are used to model the expansion of the detonation products, mixing, and turbulent combustion of the detonation products with air. The gas phase is coupled to the particles and wakes phases with drag and heat transfer terms. The conservation laws of mass, momentum, total energy, and detonation products (DP) mass fraction are where , , , , and are the gas density, velocity, pressure, temperature, and mass-fraction of DP, respectively. Here, = + • /2 is the total energy and e is the internal energy. The equation of state is specified by tables , , = ( , , ) based upon the Cheetah code [11]. Drag and heat transfer effects with the particles and wakes are included in the source terms on the right-hand side of Equations (34) and (35).

Particle Phase
The dynamics of the particle phase is modeled by discrete Lagrange particles (DLPs) and their interactions with air, namely, mass, momentum, and energy [22]:

Gas Phase
The conservation laws of gas dynamics are used to model the expansion of the detonation products, mixing, and turbulent combustion of the detonation products with air. The gas phase is coupled to the particles and wakes phases with drag and heat transfer terms. The conservation laws of mass, momentum, total energy, and detonation products (DP) mass fraction are where ρ, u, p, T, and Y D are the gas density, velocity, pressure, temperature, and mass-fraction of DP, respectively. Here, E = e + u·u/2 is the total energy and e is the internal energy. The equation of state is specified by tables p, T, a = f i (ρ, e, Y D ) based upon the Cheetah code [11]. Drag and heat transfer effects with the particles and wakes are included in the source terms on the right-hand side of Equations (34) and (35).

Particle Phase
The dynamics of the particle phase is modeled by discrete Lagrange particles (DLPs) and their interactions with air, namely, mass, momentum, and energy [22]: where x i , v i , m i , and e i denote the position, velocity, mass, and internal energy of particle i, and g denotes gravity. Mass loss, drag, and heat transfer interactions with air are specified by the following relations: . .

Wake Phase
The wake phase is modeled by the heterogeneous continuum (HC) model of Nigmatulin [27] and Collins et al. [28]: Drag and heat transfer interactions between wakes and gas are modeled by .
The same drag and Nusselt correlations of Equations (44) and (45) are used.

Scaling of Pyrotechnic Explosions
The scaling laws for the gas phase were discussed in Section 2. For pyrotechnic explosions, we must deal with all three phases. One can say that the conservation laws for phases 2 and 3 describe inertial flow, since there are no pressure gradient forces. However, there are source/sink terms on the right-hand side of the equations, expressing mass transfer, drag, and heat transfer with the gas phase.
They are functions of the Reynolds and Prandtl numbers. How these scale with change in problem scale is discussed below.
The initial conditions for a pyrotechnic charge are illustrated in Figure 2. It shows the detonation wave structure for the gas phase (as before) and the linear velocity profile for the droplet phase. This corresponds to similarity solution of inertial flow expansion into a vacuum, as published by Stanyukovich [15] in 1960. The ball radius of the droplet phase extends to 10% of the charge radius (i.e., R B = 0.1·R c ). These initial conditions scale with the explosion scale, r 0 , in other words x B = R B /r 0 = 0.1·R c /r 0 . Thus, we consider two pyrotechnic explosion scales: • A full-scale charge (the dynamics of this pyrotechnic charge system was published in the 16th Detonation Symposium [21]) where R c ≡ 100 cm with a mass of 7175 kg; • 1/3 scale charge where R c ≡ 100/3 cm with a mass of 265.7 kg.
As mentioned above [21], the pyrotechnic charge contains a liquid ball of radius R B = 0.1·R c , and thus a full-scale charge would have a ball radius of R B = 10 cm, while a 1/3-scale charge would have a ball radius of R B = 10/3 cm. The ball is made up of droplets; we estimated the droplet diameter as d i = 3 mm, which will be the same for both scales. (We assumed the liquid sphere was put under strong radial strain, creating droplets. According to Tolman (1949) [29], the droplet size is controlled by surface tension, which is a material property. Thus, the droplet size will be the same for full scale and subscale.) The number of droplets scaled as n i = (2R B /d i ) 3 , which was equal to 296,296 droplets at full scale and 10,974 at 1/3 scale. The droplet wake diameter scaled as d w = d i /100, and thus was d w = 3 µm in both cases. The wake equilibration length scaled as l w /d w = 10, 000, and thus the wake equilibration length was 30 cm in both cases. The droplet Reynolds number scaled as Re i = u·d i /υ. Assuming a characteristic velocity of 1 km/s, one finds the droplet Reynolds number of Re i = 2 × 10 5 , and the wake Reynolds number of Re W = 2 × 10 3 . These have drag coefficients of C D = 0.48 and Nusselt numbers of Nu ∼ 40 (see Appendix C for more Nusselt number correlations versus data). These results are summarized in Table 1. We note that here the particle sizes introduced an additional length scale, and thus we did not expect to maintain self-similarity.

Gas Phase Results
We investigated the scaling behavior by comparing simulations at two different scales. Pyrotechnic fireball temperature cross-sections at 1/3 scale and full scale are presented in Figure 18.
The azimuthally averaged pyrotechnic fireball temperature profiles are presented in Figure 19. Curves of the mean and RMS profiles are shown; profiles at 1/3 scale (blue curves) and full scale (orange curves) overlayed, i.e., they scaled with r 0 . It is interesting to note that the maximum mean fireball temperature in the pyrotechnic explosion was T = 1650 K, which was lower than the HE fireball case of T = 2450 K (see Figure 7a) because of heat losses to the DLPs. Pyrotechnic fireball vorticity cross-sections are presented in Figure 20. Their azimuthally averaged mean and RMS profiles are depicted in Figure 21; results for 1/3 scale (blue curves) and full scale (orange curves) overlayed, i.e., the scale with r 0 .  Pyrotechnic fireball vorticity cross-sections are presented in Figure 20. Their azimuthally averaged mean and RMS profiles are depicted in Figure 21; results for 1/3 scale (blue curves) and full scale (orange curves) overlayed, i.e., the scale with .       Pyrotechnic fireball vorticity cross-sections are presented in Figure 20. Their azimuthally averaged mean and RMS profiles are depicted in Figure 21; results for 1/3 scale (blue curves) and full scale (orange curves) overlayed, i.e., the scale with .

DLP Phase Velocity Results
Azimuthally averaged mean and RMS velocity profiles for the DLP phase of pyrotechnic explosions are presented in Figures 24 and 25. The mean and RMS velocity profiles of the DLP phase did not overlay and thus the DLP phase results did not scale with r .

DLP Phase Velocity Results
Azimuthally averaged mean and RMS velocity profiles for the DLP phase of pyrotechnic explosions are presented in Figures 24 and 25. The mean and RMS velocity profiles of the DLP phase did not overlay and thus the DLP phase results did not scale with r 0 .

DLP Phase Velocity Results
Azimuthally averaged mean and RMS velocity profiles for the DLP phase of pyrotechnic explosions are presented in Figures 24 and 25. The mean and RMS velocity profiles of the DLP phase did not overlay and thus the DLP phase results did not scale with r .

Mist Phase Results
Pyrotechnic fireball mist temperature cross-sections are presented in Figure 26. Their azimuthally averaged mean and RMS mist temperature profiles are shown in Figure 27. Results for 1/3 scale (blue curves) and full scale (orange curves) did not overlay, i.e., their profiles did not scale with r .

Mist Phase Results
Pyrotechnic fireball mist temperature cross-sections are presented in Figure 26. Their azimuthally averaged mean and RMS mist temperature profiles are shown in Figure 27. Results for 1/3 scale (blue curves) and full scale (orange curves) did not overlay, i.e., their profiles did not scale with r .

Mist Phase Results
Pyrotechnic fireball mist temperature cross-sections are presented in Figure 26. Their azimuthally averaged mean and RMS mist temperature profiles are shown in Figure 27. Results for 1/3 scale (blue curves) and full scale (orange curves) did not overlay, i.e., their profiles did not scale with r 0 .  Pyrotechnic fireball mist density cross-sections are presented in Figure 28. Their azimuthally averaged mean and RMS mist temperature profiles are shown in Figure 29. Results for 1/3 scale (blue curves) and full scale (orange curves) did not overlay, i.e., their profiles did not scale with .   Pyrotechnic fireball mist density cross-sections are presented in Figure 28. Their azimuthally averaged mean and RMS mist temperature profiles are shown in Figure 29. Results for 1/3 scale (blue curves) and full scale (orange curves) did not overlay, i.e., their profiles did not scale with .  Pyrotechnic fireball mist density cross-sections are presented in Figure 28. Their azimuthally averaged mean and RMS mist temperature profiles are shown in Figure 29. Results for 1/3 scale (blue curves) and full scale (orange curves) did not overlay, i.e., their profiles did not scale with r 0 .  Pyrotechnic fireball mist density cross-sections are presented in Figure 28. Their azimuthally averaged mean and RMS mist temperature profiles are shown in Figure 29. Results for 1/3 scale (blue curves) and full scale (orange curves) did not overlay, i.e., their profiles did not scale with .

Summary
A summary of the scaling laws for pyrotechnic explosions is given in Table 2. When implemented in the three-phase model (Equations (33)-(51)), one finds that

•
Turbulent combustion in the gas phase fireball does scale with ; • The DLP flow field does not scale with ; • The mist flow field does not scale with .

Conclusions
The gas phase of pyrotechnic explosions scaled with r , indicating that the multi-phase effects did not disrupt the scaling behavior. The DLP particle phase behaved ballistically, being only subject

Summary
A summary of the scaling laws for pyrotechnic explosions is given in Table 2. When implemented in the three-phase model (Equations (33)-(51)), one finds that

•
Turbulent combustion in the gas phase fireball does scale with r 0 ; • The DLP flow field does not scale with r 0 ; The mist flow field does not scale with r 0 .

Conclusions
The gas phase of pyrotechnic explosions scaled with r 0 , indicating that the multi-phase effects did not disrupt the scaling behavior. The DLP particle phase behaved ballistically, being only subject to drag and heat transfer without pressure (blast wave) effects, and thus did not satisfy similarity conditions. The mist phase started with zero mass and was only created by a mass sink in the DLP equations. It also was inertial (i.e., devoid of pressure effects) and only subject to drag and heat transfer effects, and thus it also did not satisfy the similitude conditions.

Initial HE Density Perturbations
When initializing the HE density profile shown in Figure 2, we introduced spatial sinusoidal variations with 1% maximum amplitude to break nonphysical numerical symmetries. This helped trigger the instabilities on the HE-air interface in order to make it less dependent on grid effects. To check whether these spatial perturbations might have any effect on scaling, we changed the sinusoidal base, and re-ran the HE fireball case shown in Figures 4-15. The result was that the profiles of the two cases overlayed, and thus scaled with r 0 .

Grid Refinement
We also explored the effect of grid refinement on scaling. We doubled the resolution and re-ran the HE fireball case shown in Figures 3-15. Table 3 summarizes the profile comparisons of 1× versus 2× zoning. We found that a factor of 2 increase in grid resolution had no effect on the mean and RMS profiles of thermodynamic variables, T, ρ, m F (t), and mean velocity-related profiles of ω, v r , v θ , v φ . The RMS profiles also scaled with r 0 , but had about 30% variation due to grid resolution. Table 3. Summary of grid resolution on mean and RMS profiles.

Conclusions
We investigated scaling behavior for turbulent combustion fields in explosions. Sedov has derived a similarity theory for explosions, stating that the flow fields U(x, y, z, τ) scale as x = r/r 0 and τ = t/t 0 , where r 0 = (m·∆H d /p a ) 1/( j+1) and t 0 = r 0 /a a . This theory was derived for the blast wave (dilatational) velocity field. However, according to the Helmholtz decomposition theorem [30][31][32], the velocity vector field u can be divided into its dilatational u ∆ and rotational u ω components (see Appendix B). The dilatational component corresponds to the blast wave solution of Sedov. The rotational component corresponds to the turbulent velocity field found in the fireball of explosions. Our algorithm numerically satisfies the same scaling behavior discretely, and thus we expect the numerical solution for inviscid gaseous explosions to show the same scaling relation. The two solutions were identical when displayed in scaled variables (i.e., the flow field profiles only differed by machine roundoff).
Numerical simulations of pyrotechnic fireballs performed at two different scales confirmed this observation; their azimuthally averaged scaled profiles overlayed, proving that they scaled with r 0 . However, the DLP flow fields and the wake flow fields of pyrotechnic explosions did not scale with r 0 because particle sizes introduced another length scale that broke similarity. For example, scaled mist temperature, density, and velocity profiles can differ by a factor of 10 in certain regions. Thus, calculations and numerical simulations of pyrotechnic explosions must be performed at each scale of interest. This has not been recognized in the past.
Author Contributions: A.K. was responsible for the gas-dynamic formulation of the three-phase models used in this research; D.G. performed the numerical simulations reported here; J.B. led the development of the high-order Godunov algorithms and Adaptive Mesh Refinement methods used in the simulations reported here. All authors have read and agreed to the published version of the manuscript.
Funding: Lawrence Livermore National Laboratory is operated by Lawrence Livermore National Security, LLC, for the U.S. Department of Energy, National Nuclear Security Administration under Contract DE-AC52-07NA27344. This work was funded by the Office of Defense Nuclear Nonproliferation Research and Development.

Conflicts of Interest:
The authors declare no conflict of interest.

Symbol
Definition ρ gas density u gas velocity vector p gas pressure E gas total energy: E = e + 0.5u·u e gas internal energy Y D mass fraction detonation products DP Y A mass fraction of air Y F mass fraction of fuel (detonation products) Y O mass fraction of oxidizer (air) Y P mass fraction of combustion products . ϕ combustion rate α stoichiometric air/fuel ratio ω vorticity vector . ω baroclinic vorticity production: . ω = ∇ρ × ∇p/ρ 2 x nondimensional radius: x = r/r 0 τ nondimensional time: τ = t/t 0 r 0 explosion length scale: r 0 = (m·∆H d /p a ) 1/( j+1) t 0 explosion time scale: t 0 = r 0 /a a m mass of the high explosives charge ∆H d heat of detonation j geometry index j = 0, 1 or 2 for planar, cylindrical, or spherically symmetric flow R nondimensional density function: R = / a U nondimensional velocity function: U = u/a a ε nondimensional internal energy function: ε = e/e a P nondimensional pressure function: P = p/p a Θ nondimensional temperature function: Θ = T/T a U ∆ nondimensional dilatational velocity: U ∆ = u ∆ /a a U ω nondimensional rotational velocity: U ω = u ω /a a x, y, z nondimensional Cartesian coordinates: x = x/r 0 , y = y/r 0 , z = z/r 0 R FB fireball radius R c high explosives charge radius R B liquid ball radius M c high explosives charge mass n i number of i droplets: n i = (2R B /d i ) 3 η i number of droplets i per unit volume: η i = 6·(2R B /d i ) 3 /π(2R B ) 3 l i droplet i equilibration length mean azimuthally averaged mean value RMS azimuthally averaged root mean squared value u ∆ dilatational velocity u ω rotational velocity D drag force . Q heat transfer n w number of particles per unit volume x i position of particle i v i velocity of particle i m i mass of particle i . s i surface recession rate of particle i e i specific internal energy of particle i d i diameter of particle i C D drag coefficient Nu Nusselt coefficient Re Reynolds number Pr Prandtl number k thermal conductivity σ density of the wake phase v velocity of the wake phase E w energy of the wake phase n w number of wake particles per unit volume: n w = σ/m w m w mass of the wake phase in cell: m w = (1/6)πd 3 w ρ w ∆ρ A mass of air in cell ∆ρ F mass of fuel in cell Subscripts 0 characteristic scale a ambient conditions FB fireball c charge i particle i D drag w wake B ball r radial velocity θ theta velocity φ phi velocity CJ Chapman-Jouguet state Given the mean Φ(R n , t), one can then define the fluctuation by Φ (R n , t) 2 = 1 δV(R n ) (R n , θ, φ, t) − Φ(R n , t) Then the root mean squared (rms) fluctuation becomes Φ (R n , t) rms = Φ (R n , t) 2 (A4) The accuracy of the azimuthally averaged variables depends on the ensemble size, N. The number of points in the ensemble is shown in Table A1.  For radii outside of 5 cm, one has 10 5 samples, which gives a good average. For radii inside 5 cm, there may not be enough points to give an accurate average.

Appendix B Helmholtz Velocity Decomposition
The complete velocity field u = u ∆ + u ω is computed by the second order Godunov scheme [16,17]. According to Equation (25), both the dilatational and rotational velocity components scale according to U ω = u ω /a a = f ω (x, y, z, τ) (A6) In particular, the above points out that the rotational velocity component u ω -i.e., the turbulent velocity field-also scales with the explosion length scale, r 0 .
Under suitable asymptotic behavior at infinity, an arbitrary velocity vector field u can be decomposed into two parts: a dilatational component u ∆ , which is irrotational, plus a rotational component, u ω , which is solenoidal: with φ and A denoting scalar and vector potentials, respectively. They satisfy the following Poisson equations: and For more details see Batchelor [30], Morino [31] and Hodge [32].