Capturing the Swirling Vortex and the Impact of Ventilation Conditions on Small-Scale Fire Whirls

The fundamental flow structure and temperature distribution of small-scale fire whirls, including tangential and axial velocities, temperature variation, and air entrainment in the lower boundary layer, were successfully captured using a generic fire field model with large eddy simulation (LES) turbulence closure. Numerical predictions were validated thoroughly against two small-scale experimental measurements, where detailed temperature and velocity distributions were recorded. Good agreement between numerical and experimental results was achieved. Normalization was also performed to compare the numerical predictions with the empirical correlations by Lei et al. (2015) developed from medium-scale fire whirl measurements. The transient development stages of small-scale fire whirls and the impact of air entrainment on the stability of the fire whirls were also investigated based on the validated numerical results. The numerical validations showed the potential of the current LES fire field model in capturing the dynamic behaviour of the fire whirl plume and performing a quantitative analysis on its onset criteria and combustion dynamics in future.


Introduction
Fire whirls (also referred to 'fire devils') are a distinctive and potentially catastrophic form of fire, which can occur in forest fires, bushfires, or intense fires in large building voids or atrium structures. Unlike normal pool fires, the flame structure of fire whirls is sustained by a swirling vortex core along the axis of the plume [1]. The vortex core could be initiated by certain air entrainment conditions [2] and the upward movement of the heated gases [3]. Due to the circular motion within the core, the turbulent intensity, air entrainment, and the mixing between air and fuel are significantly enhanced, forming an intensified and prolonged flame structure with a burning rate that is several times higher than that of normal pool fires [4]. More importantly, fire whirls emit significant radiant heat that could further ignite surrounding combustibles, causing rapid fire spread and posing great threats to the residences and structures nearby. Due to the complex nature of fire whirls, only a few experimental measurements have been carried out over the past decades [5][6][7][8].
Early experimental studies were carried out to investigate the vorticity strength within fire whirls using a rotating cylindrical screen [9] or by channelling air entrainment [10]. More recently, several small-scale experimental studies have been carried out, attempting to correlate the burning rate and

Governing Equations
The complex whirling fluid motion and heat transfer of a turbulent reacting fire whirl are predicted by the governing equations of mass, momentum, and energy. The conservation of scalar properties, such as gas chemical species, is governed by a transport equation. For the practical simulation of a weak fire whirl, several assumptions are made, including (i) low-Mach-number flow equations are considered, (ii) thermal-physical properties are constant, (iii) the ratio between mass and thermal diffusivity (Lewis number) is unity, (iv) the variable density effect is only included in the buoyancy term in the momentum equation, where Boussinesq approximation is applied, (v) no external body force and heat source are applied towards the flow, and (vi) non-adiabatic non-premixed combustion of the energy equation is adopted. Accordingly, the following form of the governing equations were utilised in this simulation study: where ρ is the mixture density, u i is the velocity vector, and p is the pressure. As the non-premixed combustion model was adopted, the heat of formation will be included in the definition of enthalpy H; and k t is the turbulent thermal conductivity which can be calculated as k t = µ t /Pr t . µ t is the eddy viscosity and Warrantz et al. suggested that the Prandtl number Pr t is 0.7 [22]. S h accounts for source term due to radiation and heat transfer to wall boundaries. They will be further discussed at a later stage. ϕ represents the scalar quantities involved in the flow system (e.g., mixture fraction, soot particulate density, and soot volume fraction). Γ is the diffusion coefficient of ϕ.
The mixture fraction combustion model was utilized in this study. Hence, all thermochemical scalars are uniquely related to the mixture fractions. Besides, the mean scalar values are extracted from the instantaneous mixture fractions by the probability density function method. The mean fluid density, ρ, therefore, can be computed as: where p( f ) is the probability density function and f is the mixture fraction. They will be further discussed at a later stage. Using Newton's, Fick's, and Fourier's laws [23] and the assumption that the Lewis number (Le) equals to 1, the molecular stress tensor, the subgrid-scale stress, subgrid enthalpy flux term, and total enthalpy in terms of resolved quantities are given by: where σ ij is the stress tensor due to molecular viscosity, µ is the molecular viscosity, δ ij is the Kronecker delta, and τ ij is the subgrid-scale stress. H can be defined by the ideal gas law, and Y i and h 0 i T re f ,i are the mass fraction and the formation enthalpy at the reference temperature T re f ,i of species i, respectively.
The unknown subgrid-scale (SGS) terms appearing in Equations (2), (4), and (6) require closure using the SGS models. The SGS turbulent flux of the scalars is modelled by the following equation: where the molecular Schmidt number Sc ϕ for mixture fraction is the specified value of 0.7 [24], while the molecular Schmidt number for soot quantities is set to 700 [25]. Zhou et al. [26] have successfully demonstrated the feasibility of applying the Smagorinsky-Lilly model [27] to investigate unsteady low-speed buoyant jet diffusion flames. Kang and Wen [28] have also adopted the Smagorinsky-Lilly model for their studies of a small-scale buoyant fire. The subgrid-scale stress tensor can be modelled according to the Smagorinsky-Lilly formulation as follows: where S ij = 1 2 ∂u i ∂x j + ∂u j ∂x i . In this study, according to Erlebacher et al. [29], the terms involving τ kk are neglected as the flow has a low Mach number in both experiments. The Smagorinsky constant C s is taken to be 0.1. In the Smagorinsky-Lilly model, the eddy viscosity is modelled by µ t = ρL 2 s S . Combining Smagorinsky eddy viscosity and Equation (9), the subgrid-scale stress tensor can be re-written as: where L s = min(κd, C s ∆) is the mixing length for subgrid scales, κ is the von Karman constant, and d is the distance to the closest wall. ∆ is the local grid-scale and can be computed according to the volume of the computation cell V:

Subgrid-Scale Combustion Modelling
For the combustion system, as the non-premixed mixture fraction combustion model was applied in this study, the basic assumption of the mixture fraction model is that all the species mass fraction are only functions of the mixture fraction [30]. Therefore, the thermochemical state of the combustion reactions can be governed by the mixture fraction f and the mixture fraction variance f 2 .
where k is the laminar thermal conductivity of the mixture and C p is the mixture-specific heat. For large eddy simulation, the transport equation is not solved for the mixture fraction variance. Instead, it is modelled as: where C var is the dynamic stress constant that equals to 0.5 and L S is the subgrid length scale.
In LES (large eddy simulation), the flame is typically not resolved by the computational grid [24]. Therefore, the chemical reaction is assumed to occur at the flame sheet, where a state relationship exists for all reactants. Under the mixture fraction assumption, the state relationship could be represented by the chemistry equilibrium assumption and experimental state relationships established by Y. R. Sivathanu and G. M. Faeth [31]. The instantaneous temperature, species mass fraction, and density can be parameterized by the instantaneous enthalpy and mixture fraction φ i = φ i ( f , H) [32]. While, in the non-adiabatic simulation, changes in enthalpy due to heat loss must be considered. In this study, heat loss was assumed to not impact the turbulent enthalpy fluctuations. Therefore, the density-weighted scalars can be calculated as: where H is the mean enthalpy and p( f ) is the shape of the assumed probability density function (PDF). Furthermore, in the simulation, the shape of the assumed PDF could be described by the β function, and this PDF shape is given by the following function of mean mixture fraction f and mixture fraction variance f 2 . where The mathematical result of Equation (17) for the non-adiabatic system was generated and stored in a 3D look-up table for the non-adiabatic system. The mean value of each mass fraction, density, and temperature are the function of f , f 2 , and H, respectively. In the flow domain, these values can be obtained by table interpolation.

Soot Formation and Radiation Models
The soot particles contributed significantly in augmenting the overall radiative heat transfer from a fire [33]. Tracking the evolution of soot volume fraction was thereby essential in the current study. The two-equation semi-empirical soot model by Moss et al. [34], that incorporates the essential physical processes of soot nucleation, coagulation and surface growth, was adopted. This model is attractive as it can use PDF in terms of mixture fraction to predict soot formation [35].
The major parameters utilized to determine soot formation were normalized radical nuclei concentration b * nuc and soot mass fraction Y soot . The Moss-Brookes model solves two more transport equations for these two parameters.
Appl. Sci. 2020, 10, 3428 where Y soot is soot mass fraction and Pr soot is the turbulence Prandtl number for soot transport that equals to 0.4 [36]. M is soot mass concentration, b * nuc is normalized radical nuclei concentration, N is soot particle number density, and N norm equals 10 15 particles. The last term in these two equations describes the instantaneous soot mass concentration and instantaneous production rate of soot particles, respectively.
where X prec is the mole fraction of the soot precursor, which is assumed to be acetylene [37]. The mass density of soot ρ soot is assumed to be 1800 kg·m −3 , and the activation temperature T a for nucleation reaction is that proposed by Lindstedt [38]. Furthermore, in the soot mass concentration equation, X sgs is the mole fraction of the participating surface growth species. Fenimore and Jones [39] assumed the collision efficiency (η coll ) to be 0.04. C a , C β , C r , M p , C oxid , C ω , and T γ are model constants. Appropriate values for model constants are obtained from Brookes and Moss [37], C a = 54 s −1 C β = 1.0 C r = 11, 700 kg·m·kmol −1 ·s −1 C oxid = 0.015 C ω = 105.8125 kg·m·kmol −1 ·K −1/2 ·s −1 T γ = 12, 100 K (activation temperature of surface growth rate) T a = 21, 100 K (activation temperature of soot inception) M p = 144 kg·kmol −1 As in the combustion system, the flow is highly turbulent. Hence, the temperature and composition fluctuations were taken into account by the probability density function. To calculate soot concentration for the Moss-Brookes model, a two-variable PDF in terms of mixture fraction and the radiative heat loss fraction can be utilized [37]. The general expression for the mean reaction rate can be: where ω soot is the instantaneous molar rate of soot production, and p( f ) and p f 2 are the PDFs of the mixture fraction and mixture fraction variance, respectively. The luminous thermal radiation from combustion products and soot is treated by solving the filtered radiative transfer equations (FRTE) for a non-scattering grey gas using the discrete ordinates method (DOM) with the S4 quadrature scheme [40]. The grey gas assumption and the S4 approximation allow for more efficient computations. The discrete radiative transfer equation is spatially filtered, resulting in the FRTE shown below for a Cartesian coordinate system in three dimensions: In this present investigation, Equation (23) neglects the interaction of the SGS turbulence-radiation.
The blackbody radiation is defined as E b = σT 4 and σ is the Stefan-Boltzmann constant. The direction cosines ξ j , η j , and ζ j represent a set of directions for each of the radiation intensities I j that span over the total solid angle range of 4π around a point in space. The integrals over solid angles are approximated using the S4 numerical quadrature, where the maximum number of discrete ordinates is 24. The discrete version of the radiation source terms that appears in the filtered energy equation is: The filtered gas absorption coefficient K a,g for the combustion products (CO 2 and H 2 O) as well as the unburnt methane fuel can be approximated according to the weighted sum of gray gases model (WSGGM). For soot, the expression from Kent and Honnery [41], K a,s = 1864 f v T, is employed. The absorption coefficient in Equation (27) is evaluated based on the sum of the filtered gas and soot absorption coefficients, i.e., K a = K a,g + K a,s .

Experimental Arrangement of Chow and Han
Numerical simulations were carried out based on two small-scale fire whirl experiments conducted by Chow and Han [14] and Hartl and Smits [7]. A schematic view of the experiment by Chow and Han [10] is shown in Figure 1a. As depicted, a liquid propane fuel bed was placed at the centre of a 0.35 m × 0.34 m vertical shaft with a total height of 1.45 m. The top of the vertical shaft was fully opened, while the shaft was enclosed by wooden walls and a transparent plastic sheet for video footage and observations. A single ventilation gap of 3.6 cm width was opened at one side for asymmetric air entrainment to initiate the onset of fire whirl. To study the influence of air entrainment on the formation and stability of the fire whirl, the vertical gap can be fully or partially opened. For measuring the time-averaged temperature profile, a thermocouple tree was placed at the centre of the small-scale fire whirl model. The thermocouple tree measurements were taken from 0.15 m to 1.35 m above the floor level with a constant spacing of 0.15 m between each thermocouple. In this study, two numerical simulations were carried out with a liquid fuel pool area of 39 cm 2 (i.e., corresponding to an approximated fire size of 3.58 kW), while the ventilation gap was fully or partially opened as shown in Figure 1a.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 27 measuring the time-averaged temperature profile, a thermocouple tree was placed at the centre of the small-scale fire whirl model. The thermocouple tree measurements were taken from 0.15 m to 1.35 m above the floor level with a constant spacing of 0.15 m between each thermocouple. In this study, two numerical simulations were carried out with a liquid fuel pool area of 39 cm 2 (i.e., corresponding to an approximated fire size of 3.58 kW), while the ventilation gap was fully or partially opened as shown in Figure 1a.

Experimental Arrangement of Hartl and Smits
Figure 1b presents the schematic view of the experiment by Hartl and Smits [7]. As depicted, a pressure gas burner (Fisher Scientific 1201 Meker-type burner, dimethyl ether) was centred inside two, staggered, plexiglass, half-cylinders which had an inner diameter of 305 mm and a height of 890 mm. Two ventilation gaps, formed by the staggered half-cylinders, induced asymmetric air entrainment for the onset of fire whirls. As shown in Figure 1b, the exit diameter was 38.1 mm, corresponding to an approximated fire size of 5.13 kW, and the burner was sealed with epoxy to create a diffusion flame. The exhaust was situated 1.5 m above the base, which was only operated after the completion of the experiment. The flow rate of gas fuel was controlled with a needle valve, and pressure in the fuel line to the burner was maintained at 207 kPa. Furthermore, a stereo particle image velocimetry (SPIV) system, which uses two Imager sCMOS cameras, was utilized to measure the instantaneous whirl circulation, velocity distribution, and inflow velocity nearing the base. In each individual PIV image, flame positions were significantly unsteady. This wandering motion made the average velocity values different from their instantaneous values. According to previous work, the time-averaged fire whirl could be considered as the quasi-steady state and the velocity profiles could correlate to the Rankine profile [42]. Therefore, only images with u radial < (u tangential 7 ⁄ ) were collected and a subset of 127 images were averaged.

Computational Domain and Boundary Conditions
Three fire whirls were predicted in this study and the case conditions are summarized in Table 1.

Experimental Arrangement of Hartl and Smits
Figure 1b presents the schematic view of the experiment by Hartl and Smits [7]. As depicted, a pressure gas burner (Fisher Scientific 1201 Meker-type burner, dimethyl ether) was centred inside two, staggered, plexiglass, half-cylinders which had an inner diameter of 305 mm and a height of 890 mm. Two ventilation gaps, formed by the staggered half-cylinders, induced asymmetric air entrainment for the onset of fire whirls. As shown in Figure 1b, the exit diameter was 38.1 mm, corresponding to an approximated fire size of 5.13 kW, and the burner was sealed with epoxy to create a diffusion flame. The exhaust was situated 1.5 m above the base, which was only operated after the completion of the experiment. The flow rate of gas fuel was controlled with a needle valve, and pressure in the fuel line to the burner was maintained at 207 kPa. Furthermore, a stereo particle image velocimetry (SPIV) system, which uses two Imager sCMOS cameras, was utilized to measure the instantaneous whirl circulation, velocity distribution, and inflow velocity nearing the base. In each individual PIV image, flame positions were significantly unsteady. This wandering motion made the average velocity values different from their instantaneous values. According to previous work, the time-averaged fire whirl could be considered as the quasi-steady state and the velocity profiles could correlate to the Rankine profile [42]. Therefore, only images with u radial < u tangential /7 were collected and a subset of 127 images were averaged.

Computational Domain and Boundary Conditions
Three fire whirls were predicted in this study and the case conditions are summarized in Table 1. The meshing system was designed by the ICEM [43]. Referring to the experimental setup, the structure mesh system was adopted in the simulation of Chow's two cases (i.e., fully and partially opened cases) as shown in Figure 2. The total size of the computational domain was 0.5 m × 0.4 m × 1.45 m including the extension of the model, to keep numerical stability at computational boundaries, which was 0.15 m in the x-direction and 0.05 m in the y-direction. To ensure sufficient mesh resolution, the size of the meshing elements was determined through the relationship between the heat release rate, ambient density, and ambient temperature. The ANSYS Fluent package was adopted in all fire whirl simulations. For the turbulence model, the Smagorinsky-Lilly subgrid-scale (SGS) model was adopted to resolve the SGS turbulent contribution. In the combustion modelling, the combustion mixture was presented by the mixture fraction, where chemical species concentrations were evaluated based on the chemistry equilibrium mechanism. The discrete ordinates radiation model was utilized to track the radiative heat transfer between surfaces and combustion species. For the fully and partially opened vertical shaft model, the total number of mesh cells was 3,200,000. The refined boundary layer cells were positioned near the walls. For the Hartl's case, as the model included two staggered cylinders, hybrid mesh (i.e., structured hexagonal and unstructured tetrahedral) was applied to discretise the computation domain, as shown in Figure 3. The top exhaust was removed and the overall vertical height was reduced to 0.89 m to simplify experimental geometry. Therefore, the total size of the computational domain was 0.6 m × 0.6 m × 0.89 m. To enhance the numerical stability, the structured mesh was positioned at the centre of the model (i.e., with the size of 0.05 m × 0.05 m × 0.89 m), where each grid was 0.001 m × 0.001 m × 0.001 m, with a total number of 2.22 million grid cells. The outer volume in the model was meshed with the unstructured tetrahedral mesh. The total number of tetrahedra cells was 9.31 million. The pressure-based method was utilized to solve filtered governing equations by applying the SIMPLE scheme for velocity-pressure coupling. The second-order implicit was used to calculate transient formulation.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 27    Due to the heat release rate was kept relatively steady in the fully developed fire whirl stage, the influence of fuel evaporation rate by the radiation feedback can be neglected. The liquid fuel bed and the pressure gas burner, therefore, were simplified as a fixed velocity mass-flow inlet. To generate diffusion flow at the fuel inlet surface, the gas pressure of inlet was set up to 0 ps. Moreover, all wall boundaries were considered as the no-slip stationary wall. At the wall boundary, refined boundary layer cells were generated to satisfy the requirement of the alternative near-wall approach for LES based on the work of Werner and Wengle [44]. In their model, the wall shear stress was resolved by the analytical integration of the power law near-wall velocity distribution. To ensure proper entertainment into the combustion zone, a user-defined function (UDF) was implemented to consider the static pressure change along the vertical direction for the initialization and far-field pressure boundaries. The pressure change was evaluated as follows: p gauge = p 0 − ρ ambient ·g·Z.

Validation of Time-Averaged Centreline Flame Temperature
In the fire experiment by Chow and Han [14], the pool fire tests were carried out for a total burning time ranging from 150 to 180 s. Only steady-state temperature measurements were made using the thermocouple tree. For numerical validations, following our previous study [20], the time-weighted averaging calculation was applied to the predicted instantaneous temperatures with a moving average window of approximately 40 s. In addition, bare-wire K-type thermocouple trees were adopted for temperature measurements at the centre of the fire whirl. According to Chow and Han [14], the stainless steel thermocouple probes with a diameter of 1.5 mm were used throughout the experiment. Due to the lack of thermal shielding, the intense heat convection and radiation near the flame could cause an error in the actual reading at the high-temperature region (i.e., at the centre of the fire) [45]. To remedy the associated errors, Wen et al. [45] proposed a correction method to revise the temperature reading based on the thermal equilibrium assumption as follows: q convective to and f rom thermocouple = q radiative to and f rom thermocouple (25) h s,conv T g,corr − T ther = ε ther σT 4 ther − q r (26) where T g,corr is the corrected gas temperature, T ther is the thermocouple-measured value, σ is the Stefan-Boltzmann constant which equals to 1.699 × 10 −8 , q r is the local radiative flux computed by the radiation model, and h s,conv is the convective heat transfer coefficient which is calculated by an empirical correlation for a sphere in crossflow [20]. Figure 4 shows the comparison between the corrected experimental readings and predicted time-averaged temperature distribution along the centreline of the fire whirl. In general, the numerical predictions were in satisfactory agreement with the corrected experimental data. The numerical scheme successfully captured the main tendency of the centreline temperature distribution. The predicted maximum temperature of 467 • C occurred at 0.08-0.1 m above the fuel bed. As depicted, the predicted temperature profile was overpredicted around 19% at 0.2-0.4 m. The prediction errors could be attributed to the following reasons. Firstly, the equilibrium chemistry model was adopted for handling the combustion process of the fire whirl. In reality, the combustion process could become non-equilibrium due to transient turbulent mixing, leading to fuel-rich or fuel-limited conditions. The equilibrium model could, therefore, overestimate the heat generation of the combustion process, leading to overpredictions of the flame temperature. Secondly, with the limited computational resources, only thirteen species, as shown in Table 2, were considered as the intermediate species in the simulation. in the present study successfully captured the centreline temperature profile of a small-scale fire whirl. To further verify the creditability of the numerical model, the predicted time-averaged excess temperatures (∆ = − ∞ ) for Case 1 and Case 2 (i.e., fully and partially opened cases, respectively) and Case 3 (i.e., model with two half-cylinders) were compared against the empirical correlation with the medium-scale fire whirls [19]. Z was the normalized vertical height (e.g., = ⁄ , ℎ = mean flame height) as shown in Figure 5. As depicted, the numerical predictions were in satisfactory agreement with the empirical correlation in the continuous flame region. The comparison indicated the similar combustion behaviour between small-and medium-scale fire whirls within the continuous flame. Numerical predictions for Case 1 and Case 2 were in good agreement with the correlations at the intermittent flame. Predictions for Case 3 exhibited a relatively different decaying slope. This could be attributed to the different entrainment mechanisms and flame structures. In this experiment, air entrainment was initiated by two, staggered, half-cylinders, which were fundamentally different from the rectangular entrainment gap used in Chow and Han [14] and Lei et al. [19] experiments. Moreover, the vertical shaft simulation presented a cylindrical fire whirl, while this model, the half-cylinder model, presented a conical fire whirl in the simulation, and the difference in flame structures will further be discussed in a later section. According to Jiao Lei et al. [12], cylindrical fire whirls have higher imposed circulation and flame velocity, and the intensity of fuel-air mixing in the intermittent region would be sustained at the same level as in the continuous flame region. The temperature of the vertical shaft model, therefore, decreased relatively slowly in the intermittent flame and plume region.  Other intermediate species generated during the incomplete combustion process were neglected. The simplified intermediate species could pose error to the combustion heat release rate, causing overprediction of the temperature profile. Lastly, measuring error could also be a factor contributing to the discrepancies in validation. Due to asymmetric entrainment from the single ventilation gap, the fully developed fire whirl tended to skew away from the air entrainment gap. Thermocouple probes were also slightly adjusted towards the tilted flame structure (as discussed in our previous work [33]). As a result, temperature measurement was not precisely located at the centre of the experimental model. Flame temperature readings were taken slightly offset from the centre of experimental shaft posing error in the measured results. In general, the numerical scheme adopted in the present study successfully captured the centreline temperature profile of a small-scale fire whirl.
To further verify the creditability of the numerical model, the predicted time-averaged excess temperatures (∆T = T − T ∞ ) for Case 1 and Case 2 (i.e., fully and partially opened cases, respectively) and Case 3 (i.e., model with two half-cylinders) were compared against the empirical correlation with the medium-scale fire whirls [19]. Z was the normalized vertical height (e.g., Z = H/H f , where H f = mean flame height) as shown in Figure 5. As depicted, the numerical predictions were in satisfactory agreement with the empirical correlation in the continuous flame region. The comparison indicated the similar combustion behaviour between small-and medium-scale fire whirls within the continuous flame. Numerical predictions for Case 1 and Case 2 were in good agreement with the correlations at the intermittent flame. Predictions for Case 3 exhibited a relatively different decaying slope. This could be attributed to the different entrainment mechanisms and flame structures. In this experiment, air entrainment was initiated by two, staggered, half-cylinders, which were fundamentally different from the rectangular entrainment gap used in Chow and Han [14] and Lei et al. [19] experiments. Moreover, the vertical shaft simulation presented a cylindrical fire whirl, while this model, the half-cylinder model, presented a conical fire whirl in the simulation, and the difference in flame structures will further be discussed in a later section. According to Jiao Lei et al. [12], cylindrical fire whirls have higher imposed circulation and flame velocity, and the intensity of fuel-air mixing in the intermittent region would be sustained at the same level as in the continuous flame region. The temperature of the vertical shaft model, therefore, decreased relatively slowly in the intermittent flame and plume region.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 27 Figure 5. Comparison of the predicted centreline temperatures for the three fire whirls against the empirical correlation proposed by Lei et al. [19].

Validation of Time-Averaged Velocity Field
To analyse the swirling motion within the fire whirl, it was crucial to thoroughly validate the predicted velocity field with experimental data. Unfortunately, detailed velocity field measurements were not conducted by Chow and Han [14]. In this study, for the very first time, the predicted velocity field was thoroughly validated against the SPIV measurements of a small-scale fire whirl conducted by Hartl and Smits [7]. To be consistent with the experimental data, time-averaging was applied to the numerical predictions for approximately 20 s after the fire whirl had been fully developed for at least 40 s. Figure 6 shows the comparison of the measured and predicted time-averaged velocity components along a horizontal line at the centre of the fire whirl at 0.08 m above the fire bed. Although the symmetrical experimental arrangement was adopted in the SPIV measurements, it is worth noting that the measured vertical and radial velocities exhibited non-symmetric behaviour. As depicted in Figure 6, the peak vertical velocity was slightly offset from the centreline of the fire. Similarly, the two, opposite, angular velocity peaks were also slightly shifted to the left from the fire centre (i.e., negative radial locations). Moreover, a noticeable angular velocity decay was observed on the right side of the figure, while no noticeable decline was observed on the left side of the figure.
As discussed in the previous section, the motion of the fire whirl was highly unstable, causing it to wander in and out of the measuring laser sheet. This could introduce uncertainties in the timeaveraged velocities leading to the non-symmetric profiles. Furthermore, the surrounding room conditions during the experiment could also contribute to the uncertainties of the measurements.

Validation of Time-Averaged Velocity Field
To analyse the swirling motion within the fire whirl, it was crucial to thoroughly validate the predicted velocity field with experimental data. Unfortunately, detailed velocity field measurements were not conducted by Chow and Han [14]. In this study, for the very first time, the predicted velocity field was thoroughly validated against the SPIV measurements of a small-scale fire whirl conducted by Hartl and Smits [7]. To be consistent with the experimental data, time-averaging was applied to the numerical predictions for approximately 20 s after the fire whirl had been fully developed for at least 40 s. Figure 6 shows the comparison of the measured and predicted time-averaged velocity components along a horizontal line at the centre of the fire whirl at 0.08 m above the fire bed. Although the symmetrical experimental arrangement was adopted in the SPIV measurements, it is worth noting that the measured vertical and radial velocities exhibited non-symmetric behaviour. As depicted in Figure 6, the peak vertical velocity was slightly offset from the centreline of the fire. Similarly, the two, opposite, angular velocity peaks were also slightly shifted to the left from the fire centre (i.e., negative radial locations). Moreover, a noticeable angular velocity decay was observed on the right side of the figure, while no noticeable decline was observed on the left side of the figure. As discussed in the previous section, the motion of the fire whirl was highly unstable, causing it to wander in and out of the measuring laser sheet. This could introduce uncertainties in the time-averaged velocities leading to the non-symmetric profiles. Furthermore, the surrounding room conditions during the experiment could also contribute to the uncertainties of the measurements.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 27 Figure 6. Comparison of the predicted axial, tangential, and radial time-averaged velocities against the experimental measurements by Hartl and Smits [7].
For numerical predictions, as all boundary conditions were symmetric in the simulation, the predicted velocity profiles were largely symmetric with only minor fluctuating noise remaining after time-averaging. As depicted, the predicted axial velocity peak centre was more closed to the middle of the model compared to the measured result, which was slightly off from the centre. Nevertheless, the predicted vertical velocity peak value was in excellent agreement with the measurements. Similarly, an excellent comparison was also shown for the angular velocity peak values. The distance between two angular velocity peaks represents the core diameter of the fire whirl. Comparing with the experimental data, the predicted core diameter was slightly underpredicted, indicating that the intensity of the swirl motion was slightly overpredicted by the numerical scheme. This was also evidenced in the slightly overpredicted radial velocity in the comparison. The discrepancy could again be attributed to the equilibrium chemistry assumption adopted in the present numerical scheme. The assumption tends to overestimate the completeness of the combustion process, which overestimates the density contrast between the fire core and surrounding air, leading to stronger turbulence suppression. According to Lei et al. [2], the flame elongation and swirling core formation are closely related to turbulence suppression. Stronger suppression could encourage a higher swirling flame structure exhibiting a higher radial velocity in the predicted results. In general, the predicted results were in good agreement with the experimental measurement. This paper presents a complete numerical scheme to capture the swirling flame structure of a fire whirl, which is thoroughly validated against the SPIV measurements.

Normalization of Tangential Velocity and Vortex Core Radius
On the basis of many previous studies [7,14,[18][19][20], velocity profile in the plane perpendicular to the whirl axis fits the Burgers vortex assumption that the vortex consists of an inner, rigid, rotating core and an outer, free vortex. According to previous studies [7,14,[18][19][20], velocity profile in the plane perpendicular to the whirl axis aligns with the Burgers vortex assumption, where the vortex consists For numerical predictions, as all boundary conditions were symmetric in the simulation, the predicted velocity profiles were largely symmetric with only minor fluctuating noise remaining after time-averaging. As depicted, the predicted axial velocity peak centre was more closed to the middle of the model compared to the measured result, which was slightly off from the centre. Nevertheless, the predicted vertical velocity peak value was in excellent agreement with the measurements. Similarly, an excellent comparison was also shown for the angular velocity peak values. The distance between two angular velocity peaks represents the core diameter of the fire whirl. Comparing with the experimental data, the predicted core diameter was slightly underpredicted, indicating that the intensity of the swirl motion was slightly overpredicted by the numerical scheme. This was also evidenced in the slightly overpredicted radial velocity in the comparison. The discrepancy could again be attributed to the equilibrium chemistry assumption adopted in the present numerical scheme. The assumption tends to overestimate the completeness of the combustion process, which overestimates the density contrast between the fire core and surrounding air, leading to stronger turbulence suppression. According to Lei et al. [2], the flame elongation and swirling core formation are closely related to turbulence suppression. Stronger suppression could encourage a higher swirling flame structure exhibiting a higher radial velocity in the predicted results. In general, the predicted results were in good agreement with the experimental measurement. This paper presents a complete numerical scheme to capture the swirling flame structure of a fire whirl, which is thoroughly validated against the SPIV measurements.

Normalization of Tangential Velocity and Vortex Core Radius
On the basis of many previous studies [7,14,[18][19][20], velocity profile in the plane perpendicular to the whirl axis fits the Burgers vortex assumption that the vortex consists of an inner, rigid, rotating core and an outer, free vortex. According to previous studies [7,14,[18][19][20], velocity profile in the plane perpendicular to the whirl axis aligns with the Burgers vortex assumption, where the vortex consists of an inner, rigid, rotating core and an outer, free vortex. The flow will be accelerated within the inner, rigid, rotating core and maximum tangential velocity will be achieved at the edge of the rotating core. On the other hand, the flow will be decelerated in the outer, free, vortex region. In this study, the tangential velocities were normalized by maximum tangential velocity (v m ) and the radial positions were normalized by the radius of the maximum tangential velocity (r v c ). The tangential vortex core radius is defined as the location where the tangential velocity equals the maximum tangential velocity at the same level [19]. In Figures 7 and 8, the normalized tangential velocities were plotted against the normalized radial position for Case 1 and Case 3, respectively. Due to density difference, the tangential velocity will decrease along the vertical direction, the turbulence suppression will reduce, but the diameter of the vortex core will increase. As depicted in Figures 7 and 8, the position where normalized tangential velocities achieved zero again tended to be closing to the centre of vertex core in both models. Within the inner, rigid, rotating core (e.g., r/r v c ≤ 1.0), the normalized tangential velocities of the half-cylinder model and vertical shaft model fitted well by the Burgers vortex line, and both of them increased quickly. Due to the smooth chambered wall surface in the half-cylinder model, the normalized tangential velocities at different heights decreased steadily within the outer, free vortex (1 < r/r v c ≤ 9). All velocities were reasonably against the Burgers vortex assumption. While, in the fully opened vertical shaft model, as a result of the effect of asymmetric ventilation condition and shaft geometry (e.g., flat wall), tangential velocities increased significantly at the near wall region and sharply decreased to zero at the wall surface, which meant that the air movement at the near wall region was related to inlet flow induced by the vertical gap. Moreover, the normalized velocity at normalized height equaled 0.3 (z = 0.3), decayed slower than the line of the Burgers vortex function, and slightly increased at the near wall region. It indicated that the outer, free vortex at this level (z = 0.3) was dominated by the air entrainment passing through the vertical gap. of an inner, rigid, rotating core and an outer, free vortex. The flow will be accelerated within the inner, rigid, rotating core and maximum tangential velocity will be achieved at the edge of the rotating core. On the other hand, the flow will be decelerated in the outer, free, vortex region. In this study, the tangential velocities were normalized by maximum tangential velocity ( ) and the radial positions were normalized by the radius of the maximum tangential velocity ( ). The tangential vortex core radius is defined as the location where the tangential velocity equals the maximum tangential velocity at the same level [19]. In Figures 7 and 8, the normalized tangential velocities were plotted against the normalized radial position for Case 1 and Case 3, respectively. Due to density difference, the tangential velocity will decrease along the vertical direction, the turbulence suppression will reduce, but the diameter of the vortex core will increase. As depicted in Figures 7 and 8, the position where normalized tangential velocities achieved zero again tended to be closing to the centre of vertex core in both models. Within the inner, rigid, rotating core (e.g., ⁄ ≤ 1.0), the normalized tangential velocities of the half-cylinder model and vertical shaft model fitted well by the Burgers vortex line, and both of them increased quickly. Due to the smooth chambered wall surface in the half-cylinder model, the normalized tangential velocities at different heights decreased steadily within the outer, free vortex (1 < ⁄ ≤ 9). All velocities were reasonably against the Burgers vortex assumption. While, in the fully opened vertical shaft model, as a result of the effect of asymmetric ventilation condition and shaft geometry (e.g., flat wall), tangential velocities increased significantly at the near wall region and sharply decreased to zero at the wall surface, which meant that the air movement at the near wall region was related to inlet flow induced by the vertical gap. Moreover, the normalized velocity at normalized height equaled 0.3 ( = 0.3), decayed slower than the line of the Burgers vortex function, and slightly increased at the near wall region. It indicated that the outer, free vortex at this level ( = 0.3) was dominated by the air entrainment passing through the vertical gap.

Formation of Whirling Motion and Its Impact on Flame Shape and Height
With the validated temperature and velocity field, it was now possible to analyse the whirling structure and its relationship with the development of fire whirls within the two experiments. In the fire experiment by Chow and Han [13], different development stages of the fire whirl were identified based on the shape of the flame structure observed throughout the experiment. In the present numerical study, following the study by Lei et al. [19], the shape flame structure was presented by the iso-surface of the mixture fraction, where its value was equivalent to the stoichiometric ratio.
and is the air-to-fuel ratio on a mass basis. Equation (27) calculates the mixture fraction at stoichiometric conditions when = 1. Figures 9 and 10 show the predicted revolution of flame height and the instantaneous flame structure after fire ignition for Case 1 and Case 3, respectively. As depicted in the two figures, both fire whirl developments shared similar characteristics at the initial stage. The flame in both experiments first emerged as a free-standing plume fire after ignition (see the first flame height peak in both figures). Afterwards, due to the air entrainment, the fire structure became unstable and highly titled above the fuel bed. At this stage, the flame height reached its lowest point in both experiments. After the lowest point, the flame height gradually increased as the self-rotational motion started to strengthen and intensify the combustion and fire heat release rate. When the fire coupled with the self-rotational motion along the centreline above the fuel bed, the fire whirl attached to the fully developed stage, where flame structure was elongated to its maximum height (see the instantaneous flame structure at 6 and 6.4 s in both figures).

Formation of Whirling Motion and Its Impact on Flame Shape and Height
With the validated temperature and velocity field, it was now possible to analyse the whirling structure and its relationship with the development of fire whirls within the two experiments. In the fire experiment by Chow and Han [13], different development stages of the fire whirl were identified based on the shape of the flame structure observed throughout the experiment. In the present numerical study, following the study by Lei et al. [19], the shape flame structure was presented by the iso-surface of the mixture fraction, where its value was equivalent to the stoichiometric ratio.
where γ = ( f uel/air) actual ( f uel/air) stoichiometric and ν is the air-to-fuel ratio on a mass basis. Equation (27) calculates the mixture fraction at stoichiometric conditions when ν = 1. Figures 9 and 10 show the predicted revolution of flame height and the instantaneous flame structure after fire ignition for Case 1 and Case 3, respectively. As depicted in the two figures, both fire whirl developments shared similar characteristics at the initial stage. The flame in both experiments first emerged as a free-standing plume fire after ignition (see the first flame height peak in both figures). Afterwards, due to the air entrainment, the fire structure became unstable and highly titled above the fuel bed. At this stage, the flame height reached its lowest point in both experiments. After the lowest point, the flame height gradually increased as the self-rotational motion started to strengthen and intensify the combustion and fire heat release rate. When the fire coupled with the self-rotational motion along the centreline above the fuel bed, the fire whirl attached to the fully developed stage, where flame structure was elongated to its maximum height (see the instantaneous flame structure at 6 and 6.4 s in both figures).  At the fully developed stage, the flame structure is strongly influenced by turbulent suppression. According to Liu et al. [18], depending on the imposed circulation, the fully developed fire whirl could exhibit three distinctive flame shapes: conical fire whirl, cylindrical fire whirl, and curved cylindrical fire whirl. Due to the different air entrainment mechanisms, the fully developed fire whirl also exhibited two different flame shapes in both experiments. As depicted in Figure 9, the predicted  At the fully developed stage, the flame structure is strongly influenced by turbulent suppression. According to Liu et al. [18], depending on the imposed circulation, the fully developed fire whirl could exhibit three distinctive flame shapes: conical fire whirl, cylindrical fire whirl, and curved cylindrical fire whirl. Due to the different air entrainment mechanisms, the fully developed fire whirl also exhibited two different flame shapes in both experiments. As depicted in Figure 9, the predicted At the fully developed stage, the flame structure is strongly influenced by turbulent suppression. According to Liu et al. [18], depending on the imposed circulation, the fully developed fire whirl could exhibit three distinctive flame shapes: conical fire whirl, cylindrical fire whirl, and curved cylindrical fire whirl. Due to the different air entrainment mechanisms, the fully developed fire whirl also exhibited two different flame shapes in both experiments. As depicted in Figure 9, the predicted flame height decreased substantially after reaching the maximum point, indicating that the fire whirl transformed from conical to cylindrical fire whirl. Meanwhile, in the Hartl and Smits experiment, the predicted result showed that the fire whirl remained as a conical fire whirl and the flame height was retained at a similar level for an extended period.

Formation of Vortex Core Due to Asymmetric and Symmetric Air Entrainment
A closer examination of the predicted flow field could provide some insights on the whirling motion by induced air entrainment. Figures 11 and 12 show the predicted instantaneous velocity field at three levels (i.e., 0.05, 0.25, and 0.45 m) above the fire bed for the two fire experiments (e.g., Case 1 and Case 3, respectively). Figure 11a-c shows the velocity distribution at the three levels when the fire whirl was at the flame rising stage (i.e., stage 2). As depicted, at the lower levels, the entrained air entered from the vertical gap and initiated a strong rotational vortex skewing towards the bottom left corner. This rotational vortex induced a leaning flame structure towards the vertical wall. At the higher level, the velocity distribution became more uniform, indicating that the rotational vortex was mainly concentrated at the lower levels. The resultant flame height was therefore substantially reduced. In the fully developed stage (i.e., stage 3), the vortex core was in cyclostrophic balance. A clear vortex core structure can be observed at all three levels (Figure 11d-f). It is noteworthy that the vortex core was slightly offset from the centreline of the fuel bed. This could be caused by asymmetric air entrainment from the single air gap, which also agreed with the experimental observations. In stage 4, the vortex core structures remained largely similar to that in stage 3. Tangential velocities between the vortex core and surrounding walls were however increased compared with that in stage 3, indicating that the flame structure transited to the cylindrical fire whirl shape. For the Hartl and Smits experiment, a similar velocity distribution and flow structure can also be observed. Although symmetric air entrainment was adopted in the experiment, as shown in Figure 11a, the predicted vortex cores were noticeably offset from the centre of the fire bed at the lower levels. This indicated that the fire whirl was at stage 2 of its development with a leaning flame structure and relatively reduced flame height. In the fully developed stage, due to symmetric air entrainment, a clear vortex core can be observed along the centreline of the fire bed. The vortex core and strength of its tangential velocity remained in the later stage (Figure 12d-f), indicating that the flame structure remained in the conical shape throughout the experiment.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 27 flame height decreased substantially after reaching the maximum point, indicating that the fire whirl transformed from conical to cylindrical fire whirl. Meanwhile, in the Hartl and Smits experiment, the predicted result showed that the fire whirl remained as a conical fire whirl and the flame height was retained at a similar level for an extended period.

Formation of Vortex Core Due to Asymmetric and Symmetric Air Entrainment
A closer examination of the predicted flow field could provide some insights on the whirling motion by induced air entrainment. Figures 11 and 12 show the predicted instantaneous velocity field at three levels (i.e., 0.05, 0.25, and 0.45 m) above the fire bed for the two fire experiments (e.g., Case 1 and Case 3, respectively). Figure 11a-c shows the velocity distribution at the three levels when the fire whirl was at the flame rising stage (i.e., stage 2). As depicted, at the lower levels, the entrained air entered from the vertical gap and initiated a strong rotational vortex skewing towards the bottom left corner. This rotational vortex induced a leaning flame structure towards the vertical wall. At the higher level, the velocity distribution became more uniform, indicating that the rotational vortex was mainly concentrated at the lower levels. The resultant flame height was therefore substantially reduced. In the fully developed stage (i.e., stage 3), the vortex core was in cyclostrophic balance. A clear vortex core structure can be observed at all three levels (Figure 11d-f). It is noteworthy that the vortex core was slightly offset from the centreline of the fuel bed. This could be caused by asymmetric air entrainment from the single air gap, which also agreed with the experimental observations. In stage 4, the vortex core structures remained largely similar to that in stage 3. Tangential velocities between the vortex core and surrounding walls were however increased compared with that in stage 3, indicating that the flame structure transited to the cylindrical fire whirl shape. For the Hartl and Smits experiment, a similar velocity distribution and flow structure can also be observed. Although symmetric air entrainment was adopted in the experiment, as shown in Figure 11a, the predicted vortex cores were noticeably offset from the centre of the fire bed at the lower levels. This indicated that the fire whirl was at stage 2 of its development with a leaning flame structure and relatively reduced flame height. In the fully developed stage, due to symmetric air entrainment, a clear vortex core can be observed along the centreline of the fire bed. The vortex core and strength of its tangential velocity remained in the later stage (Figure 12d-f), indicating that the flame structure remained in the conical shape throughout the experiment.

Imposed Air Entrainment and Its Impact on the Vortex Core
To investigate the influence of air entrainment on the formation and stability of the fire whirl, the vertical gap could be partially opened in the Chow and Han experiment. With the same aim, numerical simulations were also carried out to capture the air entrainment mechanism and its impact on fire whirl characteristics under a partially open configuration. In the present study, the vertical gap was partially opened to 0.4 m above the floor level (Figure 1a). Figure 13 shows the comparison of the predicted vortex core region of two fire whirls under the fully and partially opened ventilation conditions. The figure provides a clear visualization of the two distinguished types of behaviour of the two vortex cores. As depicted, under the fully opened ventilation condition, the vortex core exhibited a persistent, stable, rotating flow structure. The high vorticity core appeared to be a rigidlike body established above the fire bed and extended to the full height of the vertical shaft. In contrast, for the partially opened ventilation condition, the vortex core became an unstable and weak rotating structure. The height of the vortex core also reduced significantly below the height of

Imposed Air Entrainment and Its Impact on the Vortex Core
To investigate the influence of air entrainment on the formation and stability of the fire whirl, the vertical gap could be partially opened in the Chow and Han experiment. With the same aim, numerical simulations were also carried out to capture the air entrainment mechanism and its impact on fire whirl characteristics under a partially open configuration. In the present study, the vertical gap was partially opened to 0.4 m above the floor level (Figure 1a). Figure 13 shows the comparison of the predicted vortex core region of two fire whirls under the fully and partially opened ventilation conditions. The figure provides a clear visualization of the two distinguished types of behaviour of the two vortex cores. As depicted, under the fully opened ventilation condition, the vortex core exhibited a persistent, stable, rotating flow structure. The high vorticity core appeared to be a rigidlike body established above the fire bed and extended to the full height of the vertical shaft. In contrast, for the partially opened ventilation condition, the vortex core became an unstable and weak rotating structure. The height of the vortex core also reduced significantly below the height of

Imposed Air Entrainment and Its Impact on the Vortex Core
To investigate the influence of air entrainment on the formation and stability of the fire whirl, the vertical gap could be partially opened in the Chow and Han experiment. With the same aim, numerical simulations were also carried out to capture the air entrainment mechanism and its impact on fire whirl characteristics under a partially open configuration. In the present study, the vertical gap was partially opened to 0.4 m above the floor level (Figure 1a). Figure 13 shows the comparison of the predicted vortex core region of two fire whirls under the fully and partially opened ventilation conditions. The figure provides a clear visualization of the two distinguished types of behaviour of the two vortex cores. As depicted, under the fully opened ventilation condition, the vortex core exhibited a persistent, stable, rotating flow structure. The high vorticity core appeared to be a rigid-like body established above the fire bed and extended to the full height of the vertical shaft. In contrast, for the partially opened ventilation condition, the vortex core became an unstable and weak rotating structure. The height of the vortex core also reduced significantly below the height of opening (i.e., 0.4 m).
The numerical simulation was in agreement with the experimental observations and ascertained that the impact of the imposed air entrainment was crucial on the onset of fire whirls and the stability of the vortex core.
opening configuration. For the Chow and Han experiment, within the continuous flame region, there was no substantial difference between the fully and partially opened ventilation conditions. It indicated that the air entrainment behaviour in the bottom of the inflow boundary layer was similar in both cases. In the intermittent and plume regions, the mass flow rate of air entrainment decreased significantly from 1 to 0.5 due to blockage of opening. On the other hand, the air entrainment increased gradually with the height for the fully opened case. The decreased air entrainment might not be sufficient to sustain the combustion causing instability of the flame structure. This also explained the longer burning time (indicating a lower heat release rate) as reported in the experiment.   Figure 14 shows the predicted normalized mass flow rate, where (s + 1)m f is the stoichiometric value with respect to the fuel, along with the normalized height of the fire whirl plumes. The axial mass flow rate (m f ) can be obtained by the following numerical integration: Lei et al. [19] showed that, due to the stabilizing effect of fire whirl, the shear stresses are strictly suppressed. Therefore, in a fire whirl, the mass flux ρw = ρ * w. In addition, the integral upper limit r w c is the plume radius, where the axial velocity decays to the half-maximum value at the selected level, and w is the axial velocity. In general, throughout the fire whirl plume (i.e., along the continuous flame, intermittent flame, and plume region), a higher normalized mass flow rate was imposed on the fire whirl in the Hartl and Smits experiment. This was caused by the axisymmetric two-sided opening configuration. For the Chow and Han experiment, within the continuous flame region, there was no substantial difference between the fully and partially opened ventilation conditions. It indicated that the air entrainment behaviour in the bottom of the inflow boundary layer was similar in both cases. In the intermittent and plume regions, the mass flow rate of air entrainment decreased significantly from 1 to 0.5 due to blockage of opening. On the other hand, the air entrainment increased gradually with the height for the fully opened case. The decreased air entrainment might not be sufficient to sustain the combustion causing instability of the flame structure. This also explained the longer burning time (indicating a lower heat release rate) as reported in the experiment. The impact of air entrainment on the combustion process can be evidenced in Figures 15 and 16, which show the normalized radial profile of excess temperature in the continuous and intermittent regions for three fire whirl plumes. Here, the radial profile was normalized by the vortex core radius defined by the maximum tangential velocity. As depicted, in the continuous region, one can observe a noticeable temperature dip within the vortex core (i.e., below 1) in the Hartl and Smits experiment. This indicated a fuel-rich condition within the vortex core, where the excess temperature was suppressed due to limited oxygen. It is well known that the flow within the vortex core is in cyclostrophic balance, where radial fluid motion is suppressed, causing limited mixing within the core structure. The higher the swirling of the vortex core, the stronger the suppression is imposed in the mixing. A similar temperature dip can also be observed in the fully opened Chow and Han experiment. Nevertheless, an insignificant decrease in temperature was shown for the partially opened case, indicating a relatively weak rotational motion and radial motion suppression within the vortex core. Excess temperature, caused by the radial heat transfer towards the centre of the core, exhibited a plateau-like distribution within the vortex core. As depicted, the temperature profile started to decay rapidly at the vicinity of the vortex core. The decay profiles were in satisfactory agreement with the empirical decay functions proposed by Lei et al. [19]. For the partially opened case, a noticeable fluctuation in temperature can be observed along the radial direction. This was caused by the noise in the time-averaging process, where a fire whirl structure was found unstable under the partially opened condition. The impact of air entrainment on the combustion process can be evidenced in Figures 15 and 16, which show the normalized radial profile of excess temperature in the continuous and intermittent regions for three fire whirl plumes. Here, the radial profile was normalized by the vortex core radius defined by the maximum tangential velocity. As depicted, in the continuous region, one can observe a noticeable temperature dip within the vortex core (i.e., below 1) in the Hartl and Smits experiment. This indicated a fuel-rich condition within the vortex core, where the excess temperature was suppressed due to limited oxygen. It is well known that the flow within the vortex core is in cyclostrophic balance, where radial fluid motion is suppressed, causing limited mixing within the core structure. The higher the swirling of the vortex core, the stronger the suppression is imposed in the mixing. A similar temperature dip can also be observed in the fully opened Chow and Han experiment. Nevertheless, an insignificant decrease in temperature was shown for the partially opened case, indicating a relatively weak rotational motion and radial motion suppression within the vortex core. Excess temperature, caused by the radial heat transfer towards the centre of the core, exhibited a plateau-like distribution within the vortex core. As depicted, the temperature profile started to decay rapidly at the vicinity of the vortex core. The decay profiles were in satisfactory agreement with the empirical decay functions proposed by Lei et al. [19]. For the partially opened case, a noticeable fluctuation in temperature can be observed along the radial direction. This was caused by the noise in the time-averaging process, where a fire whirl structure was found unstable under the partially opened condition.     Figure 17 shows the predicted time-averaged tangential velocity along the vertical centreline of the fire bed under the fully and partially opened ventilation conditions. As depicted, at the bottom boundary layer levels (i.e., z = 0.05-0.1 m), circulation motion was imposed by the entrained air from the opening. Due to the restricted opening (i.e., less opening area), the imposed tangential velocities from the partially opened case were slightly higher than that of the fully opened one. For the fully opened case, starting from z = 0.1 m, two noticeable velocity peaks were found near the centre. It is worth noting that the two peaks were of the similar absolute magnitude, which was roughly equal to imposed air entrainment, indicating the onset of the self-rotating vortex core. On the other hand, for the partially opened case, there was no distinct velocity peak within the bottom boundary layer. The absence of tangential peaks suggested no formation of vortex core due to the relatively weak circulation motion. Similar findings are also shown in Figure 13, where vortex core was not visible in the bottom boundary layer for the partially opened case. The onset of self-rotating vortex core was only found at a higher level (i.e., z = 0.3-0.4 m). Beyond the opening (i.e., z = 0.5-0.6 m), without the imposed air entrainment, the vortex core started to collapse, while the rotating motions dissipated to the outer vortex region. The predictions suggested that imposed circulation was crucial to maintain the stability of the vortex core. In contrast, with the consistent entrained air, the vortex core remained stable in the fully opened case.

Tangential Velocity and Stability of the Vortex Core
Appl. Sci. 2020, 10, x FOR PEER REVIEW 22 of 27 Figure 17 shows the predicted time-averaged tangential velocity along the vertical centreline of the fire bed under the fully and partially opened ventilation conditions. As depicted, at the bottom boundary layer levels (i.e., Figure 17a,b), circulation motion was imposed by the entrained air from the opening. Due to the restricted opening (i.e., less opening area), the imposed tangential velocities from the partially opened case were slightly higher than that of the fully opened one. For the fully opened case, starting from z = 0.1 m, two noticeable velocity peaks were found near the centre. It is worth noting that the two peaks were of the similar absolute magnitude, which was roughly equal to imposed air entrainment, indicating the onset of the self-rotating vortex core. On the other hand, for the partially opened case, there was no distinct velocity peak within the bottom boundary layer. The absence of tangential peaks suggested no formation of vortex core due to the relatively weak circulation motion. Similar findings are also shown in Figure 13, where vortex core was not visible in the bottom boundary layer for the partially opened case. The onset of self-rotating vortex core was only found at a higher level (i.e., z = 0.3-0.4 m). Beyond the opening (i.e., z = 0.5-0.6 m), without the imposed air entrainment, the vortex core started to collapse, while the rotating motions dissipated to the outer vortex region. The predictions suggested that imposed circulation was crucial to maintain the stability of the vortex core. In contrast, with the consistent entrained air, the vortex core remained stable in the fully opened case.

Conclusions
In this article, numerical simulations were performed to attempt to thoroughly validate a generic fire field. The generic fire field incorporated all essential physical considerations in fire dynamics, including LES turbulent closure, subgrid-scale combustion with presumed PDF approach, discrete ordinates radiation model, and two-equation soot formation model. To the best of our knowledge,

Conclusions
In this article, numerical simulations were performed to attempt to thoroughly validate a generic fire field. The generic fire field incorporated all essential physical considerations in fire dynamics, including LES turbulent closure, subgrid-scale combustion with presumed PDF approach, discrete ordinates radiation model, and two-equation soot formation model. To the best of our knowledge, for the first time, the numerical predictions were thoroughly validated against two experimental small-scale fire whirl measurements, where temperature and detailed velocity field were recorded. The predicted temperature profile and velocity field at the vicinity of the vortex core were in good agreement with the experimental results. The presented validation ascertained the capacity and potential of the adopted fire field model in capturing the swirling vortex and dynamic behaviour of fire whirls. Quantitative analyses were carried out based on the validated numerical model. The major findings are summarized as follows: 1 The centreline temperature variations within the continuous flame are consistent with the empirical correlation by Lei et al. [18] based on their medium-scale fire whirl measurements. 2 The dimensionless study shows that the tangential velocity within the vortex core follows the fundamental Burgers vortex equation. In the Chow and Han experiment, the tangential velocities start to increase and deviate from the Burgers equation at approximately three times of the vortex radius, where air entrainment from the ventilation gap is imposed. In the Hartl and Smits experiment, the tangential velocities decay rapidly at r/r v c ≥ 4 due to the boundary of the two half-cylinders. 3 The development stages of the fire whirls, including the flame rising, fully developed conical fire, and fully developed cylindrical fire stages, are successfully captured by the numerical model. 4 The swirling velocity within the vortex core suppresses the fuel-air mixing and its associated combustion process, causing a noticeable temperature dip at the fire centre. A stronger swirling motion results in a more significant temperature decrease. 5 The air entrainment is crucial in providing imposed circulation for the onset of self-rotating fire whirls. A partially opened ventilation gap significantly weakens the circulation motion, causing insufficient air entrainment within the inflow boundary layer and causing instability to the fire whirl structure. The numerical findings confirm that the strong inflow boundary layer is the key factor affecting the onset of self-rotational motion of fire whirls. Funding: The financial support provided by the Australian Research Council Industrial Training Transformation Centre (IC170100032) is gratefully acknowledged.

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

Nomenclature a, b
Exponents of the beta probability density function b * nuc Normalized radical nuclei concentration C p Specific heat (J·kg −1 ·K −1 ) C s Smagorinsky constant C α , C β , C γ , C ω , C oxid Exponential constants for soot quantities C var Dynamic stress constant D Fire bed width E b Blackbody radiation