Numerical simulation of small pool fires incorporating liquid fuel motion

Combustion and Flame 213 (2020) 441–454

Contents lists available at ScienceDirect

Combustion and Flame journal homepage: www.elsevier.com/locate/combustflame

Numerical simulation of small pool fires incorporating liquid fuel motion Kazui Fukumoto a,b,∗, Jennifer X. Wen b, Manhou Li a, Yanming Ding d, Changjian Wang a,c,∗ a

School of Civil Engineering, Hefei University of Technology, Hefei 230009, China Warwick FIRE, School of Engineering, University of Warwick, Coventry CV4 7AL, United Kingdom c Anhui International Joint Research Center on Hydrogen Safety, Hefei 230009, China d Faculty of Engineering, China University of Geosciences, Wuhan 430074, China b

a r t i c l e

i n f o

Article history: Received 17 June 2019 Revised 17 July 2019 Accepted 29 November 2019

Keywords: Pool fire Gas-liquid approach Burning rate Thermocapillary force Large eddy simulation

a b s t r a c t For small-scale pool fires, Vali et al. [1] showed a pair of vortices in the liquid pool. The first vortex appeared just close to the sidewall of the container, and the second one emerged slightly away from the first vortex. Large-eddy simulations of small methanol pool fires coupled with liquid fuel convective flow were conducted using an in-house version of FireFOAM to investigate the above phenomenon. In this study, a three-dimensional liquid phase model is newly developed. The model incorporates the effects of thermocapillary Marangoni convection, buoyancy, shear stress, and evaporation. For the gas phase, the combustion model is the extended eddy dissipation concept model coupled with the laminar combustion model. This combustion model uses the viscous diffusion rate to consider laminar-turbulent transition. The predictions were in reasonably good agreement with the measured local mass burning rate, flame height and distributions of liquid temperature. The error of the mass burning rate was within 4%. The present predictions captured a pair of vortices in line with Vali et al.’s experiment [1]. Their sizes increased with increasing the liquid temperature. The Reynolds analogy could explain the sensible reason behind this trend. Shear stress and thermocapillary force caused convection in the liquid pool, and this convection formed a pair of vortices. Thermocapillary force was due to the different distributions of convective and radiative heat transfer. Sensitivity test for sub-models for the liquid phase demonstrated that their effects on the mass burning rate were all less than 5.1%. Conversely, the simulation assuming zero gravity only in the liquid phase resulted in almost 64% reduction in the mass burning rate. © 2019 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Liquid pool fires are one of the most fundamental fire scenarios, which attracted considerable attention from different investigators [1–7]. Weckman and Strong [2] investigated turbulence structures of a medium-scale methanol pool fire and provided several measurements. Vali et al. [3] investigated the influences of the base temperature of the liquid fuel and vessel materials. They found that the burning rate and flame height increased with an increase in the base temperature of the liquid fuel. Vali et al. [1] further studied the influences of the convection within a pool, confirming the existence of two vortices near the sidewall of the container. Sun et al. [4] studied the behaviour of a circular ring thin-layer pool fire. They found that the burning rate changed ∗ Corresponding authors at: School of Civil Engineering, Hefei University of Technology, Hefei 230 0 09, China. E-mail addresses: [email protected] (K. Fukumoto), [email protected] (C. Wang).

with an increase of the inner and outer diameters of the circular ring pool. Some investigators have also reported numerical studies. Wen et al. [5] performed medium-scale methanol pool fire simulation using Fire Dynamic Simulator (FDS) [6], capturing the pulsating behaviour of the fire. Wang et al. [7] predicted a medium scale pool fire using FireFOAM [8], a large eddy simulation (LES) based fire simulation solver within OpenFOAM [9]. The predicted temperature and velocity fields were in reasonably good agreement with McCaffrey’s measurements [10]. Maragkos et al. [11] modified FireFOAM 2.2.x to consider the effects of non-unity Lewis numbers and thermal diffusion, and they implemented the dynamic Smagorinsky model with a variable turbulent Prandtl number approach for pool fires. Moreover, several researchers performed coupled pool fire simulation with a one-dimensional (1D) liquid phase simulation. Prasad et al. [13] investigated structures and energetics of a pool fire with a gas-liquid coupled approach and achieved reasonably good agreement with the measurements for the temperature and

https://doi.org/10.1016/j.combustflame.2019.11.047 0010-2180/© 2019 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

442

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

Nomenclature arad Cp d D g h hmass hconv hvap xf L Le m˙ m m M p pdyn q q qrad qconv r rrad R Re Ro S Sc Sh t T Tb Tbot Tinter u V xf xj X Y Z

absorption coefficient (m−1 ) heat capacity at constant pressure (J/(kg·K)) pool depth (m) effective diffusion coefficient (m2 /s) gravitational acceleration (m/s2 ) heat/sensible enthalpy (J/kg) mass transfer coefficient (m/s) convective heat transfer coefficient (W/(K·m2 )) heat of vaporisation (J/kg) flame height (length) (m) representative length or length (if no subscript) (m) Lewis number mass flow rate/mass burning rate (kg/s) mass burning rate per area (kg/(s·m2 )) mass burning rate per volume (kg/(s·m3 )) molecular weight (g/mol = kg/kmol) pressure (N/m2 ) dynamic pressure (N/m2 ) heat release rate (W) heat release rate per volume (W/m3 ) radiative heat flux (W/m2 ) convective heat flux (W/m2 ) radius of pool (m) reflectivity gas constant (J/((kmol·K)) Reynolds number Criterion for flame volume Surface area (m2 ) Schmidt number Sherwood number time (s) temperature (K) boiling point (K) fuel temperature at the pool bottom (K) interface (surface) temperature (K) velocity or velocity scale (m/s) volume (m3 ) flame height (m) coordinate in j direction (m) volume fraction mass fraction mixture fraction

Greek

α dep εrad ηrad λ ρ μ σ Stefan σ Mar τ

thermal diffusivity (m2 /s) depth from interface (m) emissivity transmissivity thermal conductivity (W/(m·K)) density (kg/m3 ) viscosity (kg/(m·s)) Stefan-Boltzmann constant (W/(m2 ·K4 )) surface tension (N/m) time scale (s)

subscripts b boiling diff diffusion EDC EDC model f flame fu fuel first first cell from the interface

i, j, k n inter liq r rad conv solid SGS tan t vap

coordinate indexes normal direction interface liquid radial direction radiation/radiative convection/convective solid (pool wall) subgrid-scale tangential direction turbulence/turbulent vaporization/vapour

bar φ˜ φ¯

density-weighted average time-average

burning rate. Novozhilov and Koseki [14] developed a gas-liquid coupled approach to predict the burning rates for small to medium scale pool fires. Although the predictions were broadly consistent with several measurements, those were strongly dependent on a soot conversion factor. Sikanen and Hostikka [15] investigated the effect of in-depth radiation absorption in liquid fuels and reported on its significant effect. As described above, Vali et al. [1] reported the appearance of two vortices close to the container sidewall. The first vortex could result from two factors: (i) an upward flow owing to buoyancy force induced by the wall heated by a fire and (ii) convection caused by shear stress owing to a gas flow on the pool surface. The mechanisms of the second vortex was still not precise. Vali et al. [3] postulated that shear stress owing to gas recirculation near the pool centre dragged a liquid flow from the centre to the pool wall. The authors [16] later suggested that heat conduction from the container sidewall influenced the burning rate for small-scale pool fires. The above findings suggested that the liquid flow motion within the pool influenced the burning rate and fire dynamics. Such effects warrant further investigation. As 1D liquid phase simulation cannot capture heat conduction through the sidewall of the container, solving two or three-dimensional (3D) governing equation is necessary for the liquid phase. Hence, the first objective of the present study is developing and validating a new gas-liquid approach to better capture the burning rate, flame height, liquid temperature, and flow motion in a pool. The second is gaining insight into the formation and evolution mechanisms of a pair of vortices in the liquid pool with the 3D numerical predictions. The third is reporting the influence of sub-models (e.g. thermocapillary force and shear stress) on the mass-burning rate. To the best of our knowledge, no previous numerical studies have simulated a pair of vortices in the liquid pool [1]. The predictions by the gas-liquid approach will be compared with the experimental data of Vali et al. [1]. 2. Physical problem Figure 1 summarizes the underlying physics of small pool fires, where turbulent and laminar flames co-exist. The pool surface receives both convective and radiative heat fluxes from the flame. In the meantime, it also emits radiative heat flux. Near the interface, an entrained airflow drags a liquid flow. While the liquid fuel evaporates at the interface, a pair of vortices appears according to previous experiments [1]. Their existence enhances heat and mass transfer. Owing to the density gradient induced by the temperature, the liquid fuel with the large density moves downward, whereas that with the small density goes upward.

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

Fig. 1. Physics of small pool fires. The arrows indicate heat and momentum (or mass) transfer, respectively.

In-depth radiation occurs owing to the transparency of the liquid fuel. Convection and diffusion transfer heat from the top to bottom of the container, and a fuel delivery system supplies mass and heat from the bottom of the container. As the flame gradually heats the sidewall of the container, this heat transfer might also have some influences on fuel evaporation. The previously reported 1D evaporation models [13–15] cannot capture such complicated processes. They do not consider (g) shear stress balance, (i) liquid motion, (j) buoyancy, and (n) heat transfer from the sidewall of a vessel. However, the relative influences of these factors are unclear, and their quantification is also one of the objectives of the present study. Notably, the present study is focused on the steady-state condition, as in the experiment of Vali et al. [1]. Transient pool fires are more complex, involving (i) initial growth, (ⅱ) steady burning stage, (ⅲ) boiling burning and (ⅳ) decay stages [29]. In most laboratory pool fire tests, there would not be continuous fuel supply to maintain the pool height. The lip height and mass burning rate changed with time. In phase (ⅲ), many bubbles formed owing to evaporation would enhance heat exchange with the sidewall. Namely, to better capture the transient pool fire scenarios, a change in a lip height and bubble motion should be considered. 3. Methodology The present study uses an in-house version of FireFOAM 2.2.x, which contains the authors’ previous developments [12,17,18,21]. The momentum and continuity equations’ coupling algorithm is Pressure-Implicit with Splitting of Operators (PISO) with the outer iteration (termed PIMPLE [9]). The energy and species equations are also solved. The subgrid-scale (SGS) modelling is the dynamic Smagorinsky model [19,20] to consider laminar-turbulent transition as denoted as (b) and (c) in Fig. 1. The methodology for gasphase modelling was the same as presented in our previous stud-

443

Fig. 2. Flowchart of the solver. Top and bottom parts denote the gas and liquid phases, respectively.

ies [21]. The following sections introduce some further modifications to facilitate the present study. Figure 2 is a simple flowchart of the in-house version of FireFOAM, excluding sub-models for simplicity. The solver consists of the combination of two PIMPLE algorithms for the gas and liquid phases marked in Fig. 2. 3.1. Combustion modelling In small scale pool fires, the flame undergoes laminar to turbulent transition [22]. Considering this phenomenon is needed for modelling the combustion process as denoted as (b)–(c) in Fig. 1. In our previous study [18], the time scale of laminar combustion τ diff was estimated based on viscous diffusion. On the other hands, Chen et al.’s EDC model for LES [12] gave the time scale of turbulent combustion τ EDC. After that, their minimum was the reaction time scale. Moreover, this study adopts the infinitely fast chemistry assumption described in our recent paper [21]. 3.2. Radiation modelling The treatment of radiative heat transfer must be taken into account in this simulation as denoted as (e) in Fig. 1. The radiative heat transfer equation was solved using the finite volume discrete ordinates method. The weighted sum of the grey gas (WSGG) model [23,24] gives the gas radiative properties. The WSGG model estimates the total absorption coefficient for a non-sooty flame as

arad = −

ln(1 − εrad ) Lbeam

1/m,

(1)

where the WSGG model gives ε rad . The assumed flame shape estimates the mean beam length Lbeam . Although the banded formulations of WSGG avoided the calculation of the mean beam length, this treatment increased the computational time by 2.9 times [24].

444

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

Here, the flame shape was estimated from preliminary simulation and approximated to be of the conical shape:

3.6Vf = Sf

Lbeam =

3.6Vf

π rf (rf + Larc,f )

m,

(2)

where Vf is the flame volume, Sf is the surface area of the flame, rf is the radius of the cone, and Larc,f is the arc length of the cone. Larc , f is given as



Larc,f = Hf2 + rf2

0 . 5

m,

(3)

where rf is the radius of the cone (equal to the pool radius), Larc,f is the arc length of the cone, and the height of the cone Hf is given as

Hf =

3Vf rf2 π

m,

In the above, Vf is computed following Yang et al. [25], who used the following criterion to determine the flame border given as

Ro =

1 1 + sY˜fu /Y˜O2

,

(5)

The flame border is at 0 ≤ Ro ≤ 0.7 for methanol which is determined to fit the predicted flame height. 3.3. Liquid phase modelling 3.3.1. Governing equations and heat balance The evaporation model is developed based on that of Ali et al. [26] and Sikanen and Hostikka [15] with some modifications and extension to 3D formulations. The governing equations for the liquid phase are based upon those for the gas phase equation as follows: Continuity equation

∂ ρliq ∂ ρliq u j,liq + = 0, ∂t ∂xj

(6)

Momentum equation

  ∂ ρliq ui,liq ∂ ρliq u j,liq ui,liq ∂ ui,liq ∂ u j,liq ∂ + = + μ ∂t ∂xj ∂ x j liq ∂ x j ∂ xi  ∂ pdyn,liq ∂ ρliq 2 ∂ uk,liq − δi j − − g jx j , 3 ∂ xk ∂ xi ∂ xi pdyn,liq = pliq − ρliq g j x j

(7) (8)

∂ ρliq hliq ∂ ρliq u j,liq hliq + ∂t ∂xj

  D pliq ∂h ∂ ρliq αliq liq = + Dt ∂xj ∂xj  ∂ q liq,rad − m hvap − , ∂ xn

hliq =

293.15

C pliq dT

J/kg



Xfu,equil = exp −

hvap Mfu R



1 Tliq,inter

−



1 Tb

(11)

where Tliq,inter = Tliq,first because of the consideration of heat balance in the first cell next to the interface. Section 3.3.3 will explains this treatment. The concept of Stefan diffusion gives the mass burning rate per area m :

m = hmass m =

p¯ first Mfu ln RTfirst

m Sinter Vfirst



Xfu,first − 1 Xfu,equil − 1







kg/ s · m2 , and

kg/(s · m3 ),

(12) (13)

Eq. (13) is unit conversion because this model considers heat balance in the first cell. Sikanen and Hostikka assumed that diffusion of vaporised fuel via a stagnant film influences m . However, a thickness of the stagnant film δ film is an unknown quantity, and they modelled it by the following mass transfer correlation:

hmass =

Dvap,gas

=

ShDvap,gas

δfilm L μ Dvap,gas = m2 /s, ρ¯ Sc

Sh = 0.664Sc1/3 Re1/2 (9)

where m hvap is considered in only the first cell from the interface. Some variables exist in both the liquid and gas phases. For such variables, subscript ‘liq’ indicates the liquid phase. A variable is one for the gas phase, unless stated otherwise. hliq is defined as T

3.3.2. Evaporation The proposed approach utilises Sikanen and Hostikka’s model [15] for evaporation calculation as denoted as (h) in Fig. 1. The approach treats the volume fraction of fuel vapour above the interface under the equilibrium state Xfu ,equil as a function of Tliq,inter and Tb based on the Clausius-Clapeyron relation:

m/s, and

(14) (15)

where Sc = 0.6 [15]. Dvap , gas is the binary diffusion coefficient between a fuel vapour and gas; Meredith et al. [32] estimates Sh as

Sensible enthalpy equation



Fig. 3. Overview of heat balance in a small pool fire.

(4)

(10)

Figure 3 gives the overview of the heat balance in a small pool fire, which is considered by boundary and interface conditions. qconv is the convective heat flux, qrad is the radiative heat flux, qre is the re-radiative heat flux owing to surface emission and reflection, qsidewall is the conductive heat flux from a wall, m hvap is the enthalpy loss per time owing to vaporisation, m hliq is the enthalpy gain/loss owing to mass gain/loss per time.

Re < 5.0 × 105 ,

(16)

Sh is evaluated using field variables at the cell centre in the first cell, Re = ρ| ¯ u|L/μ, and u is the velocity vector (m/s). Sikanen and Hostikka [15] calculated Sh for Re > 5.0 × 105 . However, in the present study, Re is typically < 100, and hence Eq. (16) is used. Beji et al. [28]. estimated L as

L=

π r2 m, 2π r

(17)

Since m is the exchange of momentum at the interface, it gives the interface velocities in the normal direction un,inter and un,liq,inter . The proposed approach exchanges momentum in the vertical direction at the interface as follows:

un,liq,inter = un,inter =

m

ρliq

m/s for the liquid phase, and

m m/s for the gas phase, ρ¯

(18) (19)

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

445

Eq. (19) is the inlet condition for the gas phase, whereas Eq. (18) is the outflow condition for the liquid phase. It is noteworthy that the mass fraction of fuel in the gas phase is given comparing the mass flux and diffusion flux [18].

the equilibrium state; the gas and fuel vapour volume fractions are as follows:

3.3.3. Sensible enthalpy balance near the interface This approach assumes that evaporation occurs very near the interface following previous evaporation models [13–15,26]. The assumption is valid if the bubbling layer is so small that its thickness is not essential for evaporation calculation. This thickness is the unknown quantity, and the presumption may be invalid at the boiling burning to the decay stages [29] because of the growth of the bubbling layer. It is noteworthy that the proposed model estimates heat balance in the first cell. Its size corresponds to the bubbling layer and must be larger than the thickness of the bubbling layer. The primary reason behind this treatment is that it prevents infinite mʺ with very high heat flux at the interface. The huge mʺ could occur close to the rim, resulting in the low liquid temperature near the sidewall. Prasad et al. [13] predicted the lowest liquid temperature (= 322.5 K) at the rim owing to the heat loss caused by the vast mʺ. Conversely, Vali et al. [1] observed the highest temperature there, and it ranged from approximately 331 to 334 K. Convective heat transfer and surface emission are considered at the interface as denoted as (d) and (f) in Fig. 1. The present approach models the interface conditions for sensible enthalpy as

Xgas = 1 − Xfu,vap ,

−λliq

∂ Tliq 4 = qconv − εrad,liq σStefan Tliq ,inter ∂xj for the liquid phase, and

(20)

Xfu,vap = Xfu,equil ,

qconv = λinter

(T˜first − Tinter ) W/m2 , first /2

(22)

where first is the first cell size along the interface. The gas-phase properties are sensitive to the temperature and concentrations of the chemical species. Therefore, the present approach uses multi-component formulation for heat conductivity at the interface λinter as [27]:

λinter =

 I

1+



λI

K=I XK · 1.065IK

W/(m · K ), and

(23)

The correction factor of chemical species I and K IK is calculated as:

1

IK = √

2 2

1+

MI MK

− 12 

1+

μ  12 M  14 2 I

μK

K

MI

,

(24)

The species I and K at the interface considered only fuel vapour and gas in Eqs. (23)–(24), but the gas here does not contain the fuel vapour. The Sutherland equation [9,11] and the modified Eucken correlation [30] estimates μ and λ of the gas, and Mgas is considered the mean molecular weight. References [28,33] list several properties of the fuel vapour for methanol, ethanol, heptane and acetone. At the interface, the concentration is assumed to be

(26)

3.3.4. Shear stress balance near the interface Momentum exchange between the gas and liquid phases occurs at the interface. The interface velocities in the tangential direction for the gas and liquid phases are given based on continuity of shear stress balance [26] as denoted as (g) in Fig. 1. Solving the following equation gives the liquid interface velocity uliq,inter,tan in the tangential direction:

μinter −

˜ first,tan − uliq,inter,tan u uliq,inter,tan − uliq,first,tan = μliq,inter first /2 first,liq /2

∇tan σMar liq,first

for the liquid phase

(27)

The gas-phase interface velocity in the tangential direction uinter,tan is given as

uinter,tan = uliq,inter,tan

for the gas phase

(28)

where ∇ tan is the derivative in the tangential direction. This treatment allows to interact convection and diffusion. μinter in Eq. (27) is given by [27]:

 I

(21)

where qrad , m hvap , and m hliq do not appear in Eq. (20). Since radiative heat flux is computed based on Beer’s law, its effect does not appear in Eq. (20). The second term on the right side includes only surface emission at the interface. Section 3.3.7 explains in-depth radiation based on Beer’s law. m hvap is expressed as m hvap in Eq. (9). Since an outflow at the interface expresses the effect of m hliq , Eq. (20) does not contain it. Above balanced heat in the first cell transfers from the top to pool bottom as denoted as (l) in Fig. 1. qconv is given as

(25)

where Xfu ,equil indicates a fuel vapour at the equilibrium state.

μinter = Tinter = Tliq,inter for the gas phase

and

μI

1+



K=I

XK XI

IK

kg/(m · s ),

(29)

where the procedure of Eq. (29) is the same as Eqs. (23) and (24). 3.3.5. Sidewall condition The proposed approach takes account of conductive heat transfer from the sidewall to the liquid fuel as denoted as (n) in Fig. 1 [16]. A gas-liquid-solid coupled approach is required to calculate the wall heat flux qsidewall , and the methodology of the prediction is considered the future work. For simplicity, the present approach utilises the measured wall heat flux as the boundary condition for the sensible enthalpy of the liquid phase at the sidewall as

qsidewall =

Twall,solid − Tliq,wall

r/λwall,solid + r/λwall,liq

W/m2 ,

(30)

where Twall,solid is the measured wall temperature obtained from reference [16], Twall,liq is obtained by the predicted temperature on the wall, ࢞r = 5 × 10−4 m [16], λwall,solid = 1.4 W/(m·K) for quartz [16], and λwall,liq are the liquid heat conductivity on the wall. The no-slip condition is the sidewall treatment for the momentum equation as follows:

uliq = 0 m/s,

(31)

3.3.6. Pool bottom condition Tbot is specified at the pool bottom for the sensible enthalpy boundary condition as denoted as (m) in Fig. 1:

Tliq = Tbot

K

(32)

The mass flow rate m˙ is the same as the mass burning rate given at the pool bottom for the momentum condition in the following experiment [1], assuring mass conservation in the pool as:

m˙ = m Sinter kg/s

(33)

446

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

Fig. 4. Properties of methanol and ethanol [33].

3.3.7. In-depth radiation The proposed model takes account of radiative heat transfer in the liquid phase as denoted as (k) in Fig. 1. Following our previous study [18], qrad is computed by Beer’s law, neglecting emission inside the pool [18]:

qliq,rad = ηrad,liq qrad exp(−arad,liq dep )

W/m2 ,

(34)

where qliq,rad decreases with increasing dep because qliq,rad is absorbed in liquid fuel. 3.4. Model parameters and physical properties Unlike previous studies [13,14,28] which assumed the liquid properties to be constant, these properties are calculated as a function of T plotted in Fig. 4 [33]. It is important to note that ρ liq must be a function of T to consider the buoyancy effect in computer simulation. Eqs. (23) and (24) use λvap and μvap . LeO2 = 1.11 [31],

LeH2 O = 0.83 [31], LeCO2 = 1.39 [31], LeN2 = 1.0 [31], LeCH3 OH = 1.31 and LeC2 H5 OH = 1.72, where the Lewis numbers of methanol and ethanol were estimated by Fuller et al.’s correlation [34]. Also, Table 1 lists the input parameters used in this approach.

Table 1 Input parameters for simulations. Parameters

Values

Prt and Sct arad ,liq

0.85 1140 [1/m] [15] 0.95 [28] 0.02 (= water’s value) 0.98 (1- rrad,liq ) 337.85 K [33] 351.45 K [33] 0.7 (numerical calibration) 0.85 (numerical calibration)

ε rad,liq rrad,liq

ηrad,liq Tb Tb Ro Ro

for for for for

methanol ethanol methanol ethanol

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

447

the computational domain was created based on the experimental condition [1]. The respective dimensions in Fig. 5 were: xmax = 0.6 m, ymax = zmax = 0.27 m, r = 0.045 and δ = 0.0025 m; the coordinate origin (0, 0, 0) was set on the pool-centre surface. m˙ = 5.9 × 10−5 kg/s, the temperature at the interface was T = 337.85 K and temperature at the rim was T ≈ 385 K, where the values of m˙ and temperature at the rim were obtained from the experiment for Tbot = 268.15 K [1,16]. The present study tested two computational grids for grid sensitivity of the gas phase. The initial number of cells across the pool diameter was 32; then the first grid layer was locally subdivided in all the directions, resulting in 64 cells on the pool surface. In the y and z directions, the cell sizes were reduced to rrim = 0.751 mm close to the rim. The first grid size from the interface xf irst was 0.5 mm, with the expansion ratio ≈ 3–4% and a total of 228,288 cells. For a more precise computational grid, the number of cells in each direction was increased by 25%; as a result, there are 448,500 cells. Figure 6 shows (a) ux /q1 /5 , (b) ࢞T, (c) flame height xf , and (d) heat fluxes predicted at different grid resolutions. ࢞T is calculated as ࢞T = T – 293.15, q is given as q = ࢞hcomb × m˙ , xf is given at the highest location with the stoichiometric mixture fraction following previous studies [7,18], qnet = qtot − qre , and qre is the absolute value of re-radiative heat flux. Differences in xf , ux /q1 /5 and heat fluxes are small, whereas ࢞T is slightly different. Since the predicted xf at different resolutions are similar, the grid resolution with 228,288 cells would be sufficient for this gas-phase simulation. The following sections create a gas-phase computational grid based on this coarser one. Figure 7 shows the predicted flame volume defined as 0 ≤ Ro ≤ 0.7. The flame height agreed the tip of flame volume as confirmed in Fig. 6. An instantaneous flame shape is similar to that observed in the experiment by Hayasaka [36], referred to as the mushroomcap shape.

Fig. 5. Sketch of the computational geometry.

3.5. Numerical conditions 4.2. Grid sensitivity test for liquid phase simulation The discretization scheme for the convection term of the momentum equation was the second-order central linear scheme. Those of the mass fraction of chemical species and the energy equation of the gas phase were the central linear scheme limited by total variation diminishing (TVD). The discretization schemes for all the equations for the liquid phase were the central linear scheme. The discretization scheme for time marching was the second-order backward differential scheme. The outer iteration and number of the pressure correction were set to 7 and 2 for the gas phase simulation, respectively, while 3 and 2 for liquid phase simulation. For gas-phase combustion, the irreversible one-step chemistry model is as follows:

Cx Hy Oz + (x + y/4 − z/2)/O2 → xCO2 +yH2 O.

(35)

The small pool fires tested by Vali et al. [1] were simulated to validate the gas-liquid approach. In the experiment, Vali et al. used a circular burner made of quartz with an inner diameter = 90 mm (2r), outer diameter = 95 mm and height (pool depth d) = 12 mm and filled methanol in the burner (= container). A fuel delivery system continuously refilled the fuel to make sure that the surface height was equal to the rim height. In other words, the lip height was zero. 4. Results 4.1. Gas phase simulation This section explains the results of gas-phase simulation to obtain the essential information of the pool fire. As shown in Fig. 5,

A gas-liquid coupled simulation was conducted to test grid sensitivity. The respective dimensions in Fig. 5 were: xmax = 1.0 m, ymax = zmax = 0.405 m, d = 0.012 m, r = 0.045, and δ = 0.0025 m. In the liquid region, divisions in the x-direction were 26, 40, and 57 with the expansion ratio of 3%; the respective first grid sizes in the x-direction ࢞x were 0.15 mm, 0.301 mm and 0.0751 mm. ࢞y and ࢞z were the same as the gas phase sizes. The initial temperatures were set from the measured temperatures to save the computation time. The simulations continued until the predicted mass burning rate became the quasi-steady state. Figure 8 depicts the distribution of the predicted T and velocity vector fields with different grid resolutions in the liquid phase. The predicted mass burning rate m˙ in cases (a)–(c) in Fig. 8 are 6.0956 × 10−5 , 6.1452 × 10−5 and 6.3428 × 10−5 kg/s. Although the grid resolution slightly affects the mass burning rate, its difference between cases (a) and (b) is only about 1.3%. Hence, the resolution in case (b) is considered to be sufficiently accurate. Two vortices emerge in the predictions for cases (a) and (b), but not in case (c) which has the least grid resolution. The left vortex is the first vortex, and the right one is the second vortex following Vali et al. [1]. A difference between cases (a) and (b) is unclear, but the interface velocity is slightly different near the second vortex. The computational grid for case (b) is sufficient, but subsequent calculations use the resolution in case (a) to give better accuracy. Additional two cases show grid dependency of the first grid height from the interface in the gas phase, which is ࢞first in

448

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

Fig. 6. (a) ux /q1/5 , (b) ࢞T, (c) flame height xf and (d) heat fluxes predicted with different grid resolutions.

Fig. 7. (a) Mean and (b) instantaneous flame volumes defined as 0 ≤ Ro ≤ 0.7.

Eq. (22). Eq. (22) would be valid if the first grid size from the interface is adequately small. The respective ࢞first were= 0.421, 0.25, and 0.15 mm, and Fig. 9 plots mass burning rate m˙ . A difference in the mass burning rate between ࢞first = 0.421 and 0.25 mm seems large, but its error is 1.1%. Therefore, the following sections employ this medium grid size ࢞xfirst ,gas = 0.25 m. Determined sizes are: ࢞xfirst,liq = 0.0715, ࢞xfirst ,gas = 0.25, ࢞yliq,gas,min = 0.377 and ࢞yliq,gas,min = 0.95 mm with 533,952 cells for the gas phase and 569,088 cells for the liquid phase. 4.3. Validation of gas-liquid approach This section performs further simulations for validation of the gas-liquid approach. The dimensions in Fig. 5 were the same as

Fig. 8. Predicted distributions of temperature and velocity vector fields with different grid resolutions.

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

449

Fig. 12. Predicted and measured flame height xf vs mass burning rate m˙ .

Fig. 9. Dependency on the first grid size from the interface.

Fig. 10. Predicted and measured mass burning rate m˙ vs fuel temperature at the pool bottom Tbot .

Fig. 11. Predicted and measured flame height xf vs fuel temperature at the pool bottom Tbot .

Section 4.2. The respective Tbot were 268.15, 280.15, 297.15, and 313.15 K. Following a change in Tbot , the temperatures at the rim and Twall,solid in Eq. (30) = 385, 388, 393, and 398 K depending on Tbot obtained from the experimental data [16]. The simulations continued until reaching the quasi-steady mass burning rate. Figure 10 shows the mass burning rate vs Tbot . The mass burning rate increases with increasing Tbot because the temperature gradient in the pool decreases with increasing Tbot , resulting in the higher surface temperature. The errors between the predictions and measurements are within 4%. Figure 11 shows the predicted and measured flame height vs Tbot . The flame height is slightly overpredicted. Since a change in the mass burning rate is small, its trend is not clear. Figure 12 depicts the predicted and measured flame height vs mass burning rate and also plots McCaffrey’s correlation [10] as a reference. The predictions are slightly overestimated compared with the measurements but are close to McCaffrey’s correlation.

Figure 13 presents the normalised temperature T∗ vs normalised depth x∗ . T∗ is calculated as T∗ = (Tliq − Tbot )/(TInter − Tbot ) and x∗ is calculated as x∗ = 1 + x/d. T∗ = 0 and x∗ = 0 are the values at the pool bottom. The predicted T∗ for Tbot = 268.15–297.15 K are in reasonably good agreement with the measured data, whereas the one for Tbot = 313.15 K deviates at x∗ > 0.5 from the experimental data. Tendencies for the predicted T∗ with Tbot = 268.15–313.15 K are similar. Namely, T∗ linearly increases at x∗ = 0–0.6 or 0.7 and that is almost constant after that. Trends in the measured T∗ for Tbot = 268.15–297.15 K are similar to the predictions; however, the one for Tbot = 313.15 K slowly increases at x∗ > 0.6. Vali et al. [1] did not mention this tendency difference, and a reason behind this change was unclear. Figure 14 depicts the predicted the local (a) heat fluxes, (b) local mass burning rate per area m vs distance from the wall, and (c) local mass burning rate per area vs total heat flux. Differences in heat fluxes and m for Tbot = 268.15–313.15 K are minimal in Fig. 14(a) and (b). The predicted m close to the sidewall are slightly different for Tbot = 268.15–313.15 K. Vali et al. [1] did not measure the mass burning rate vs total heat flux. Therefore, discussion here uses different experimental data. Singh and Gollner [35] conducted tests with a porous noncombustible material soaked with fuels; they measured the local heat fluxes and m in the vertical laminar boundary layers of methanol and ethanol. It is currently challenging to simulate their configurations because of a lack of necessary information of the porous media such as heat conductivity, density and heat capacity. Assuming that the local mass burning rate is a function of the local total heat flux, their data are comparable with the present scenarios. However, a minor difference in m owing to a difference in rear side temperature (equal to a difference in Tbot ) can happen. Figure 14(c) shows the local mass burning rate vs total heat flux. The measured results are close to the predicted data, and a decrease in Tbot slightly reduces m , especially at the higher total heat flux. A difference in the predicted m between Tbot = 268.15 –313.15 K is due to m h from the rear side, as stated above. The predicted m at the wall does not follow the trend because the total heat flux and heat flux from the sidewall influence m . Figure 14(d) shows the predicted m vs total heat flux for ethanol with Tbot = 297.15 K. Although the main focus here is the methanol pool fire, Fig. 14(d) confirms the results of a minimum applicability test. The predictions are close to the experimental data. 4.4. Analysis of the liquid phase Figure 15 depicts the distributions of liquid temperature. x is the vertical direction, and y is the horizontal direction. As shown in Fig. 15, the temperature gradient is smaller for Tbot = 313.15 K. It causes a small heat loss from the pool bottom, resulting in the high

450

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

Fig. 13. Predicted and measured normalised temperature T∗ vs normalised depth x∗ along the centreline.

mass burning rate and high surface temperature. The surface temperatures are highest just near the rim owing to heat transfer from the sidewall. The surface temperatures are relatively low at −40 < y < −35 mm and are somewhat high at −30 < y < −25 mm. This non-uniform temperature distribution is due to a difference in the distributions of convective and radiative heat fluxes as observed from Fig. 14. The convective heat flux peaks near the sidewall, whereas the radiative heat flux peaks approximately around the middle coordinate between the wall and pool centre. Also at −30 < y < −25 mm inside the pool, the temperatures are slightly higher because of in-depth radiation. Figure 16 shows the predicted velocity vector fields in the liquid phase with the velocity magnitude defined as (ux ,liq 2 + uy,liq 2 )0.5 .

There is a pair of vortices near the sidewall, and their sizes increase with increasing Tbot . This trend is in line with the experiment of Vali et al. [1]. As shown in Fig. 4, μliq decreases as Tliq increases. Assuming that the definition of a Reynolds number is Re = uliq, inter ·r/μliq , the effect of inertia in the pool is stronger as the liquid temperature increases. The effect of decreasing μliq is hence similar to that of increasing uliq, inter . The interface velocities above the first vortex centre are approximately 0.02 to 0.025 m/s. Vali et al. [1] did not measure the velocity magnitude. However, they reported approximately 0.02 m/s for Tbot = 311.15 K and d = 18 mm above the first vortex centre in a subsequent work [37]. The vortex position changed depending on Tbot in the experiment, but the present simulations do not cap-

Fig. 14. Predicted local (a) heat fluxes, (b) local mass burning rate per area m vs distance from the sidewall, (c) m vs total heat flux qtot and (d) m vs qtot for ethanol.

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

451

Fig. 16. Predicted velocity vectors in the liquid phase. Fig. 15. Predicted distributions of Tliq and velocity vector fields in the liquid phase.

ture it. A possible explanation is that the predicted liquid temperature is not sufficiently precise to capture this change. As shown in Fig. 14(c), a trend in local mass burning rate vs total heat flux is very different near the sidewall. This high total heat flux causes a flow toward the pool centre direction, owing to thermocapillary force. Figure 17 shows the predicted (a) gas and liquid temperatures and (b) velocity vector fields near the interface for Tbot = 297.15 K. In Fig. 17(a), the highest liquid temperature appears near the sidewall, with the shortest distance between the flame and interface (termed a standoff distance). As shown in Fig. 14, the convective heat flux is very high near here. Near the sidewall, an airflow toward the flame arises, dragging the liquid fuel toward the pool centre. Notably, Vali et al. [1] suggested a gas circulation near the pool centre is a reason behind the appearance of the second vortex. However, the present prediction does not show this gas circulation. Figure 18 depicts the predicted temperature and velocity vector fields near the two vortices for Tbot = 313.15 K, with bounded Tliq for 333–336 K. The temperature is very high near

the sidewall, resulting in a flow toward the pool centre because of thermocapillary force. The airflow dragged a flow toward the pool centre as confirmed in Fig. 17. In other words, shear stress and thermocapillary force cause the flow from the sidewall to the centre. Furthermore, a minimal difference in Tliq appears between the first and second vortices. Therefore, a liquid flow moves from the pool centre toward the sidewall owing to thermocapillary force. The two flows collide at the low-temperature region, which causes a downward flow. This downward flow separates into two flows: the one flow toward the sidewall and the other one toward the pool centre. Besides, a relatively large amount of evaporation occurs owing to the large heat flux near the sidewall. This evaporation causes an upward flow along the sidewall, as expressed in Eq. (20). Eq. (18) gives un ,inter = 7.9237 × 10−5 m/s near the rim for Tbot = 297.15 K, where m = 0.059586 kg/(s·m2 ) and ρ liq = 752 kg/m3 obtained from Figs. 4 and 14. Thus, the effect of evaporation is much smaller than that of shear stress and thermocapillary force. Oztop and Dagtekin [39] numerically investigated a cavity flow driven by two moving walls, giving two opposite flows just like the current scenario. In their study, left and right walls moved upward

452

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

Fig. 17. Predicted (a) gas and liquid temperatures and (b) velocity vector fields near the interface for Tbot = 297.15 K. Vector lengths are not consistent between the liquid and gas phases owing to those magnitude differences.

in a closed vessel, reporting that a clockwise vortex appeared on the left side and an anticlockwise vortex emerged on the right side. Figure 19 summarizes the overview of the flow mechanisms in the small pool. Vali et al. [1] did not discuss a small difference in the temperature between the vortices. However, such a difference emerged for Tbot = 313.15 K in their experiment [1]. In the further investigation of Vali et al. [38], there was a more distinct difference. 4.5. Gas-liquid phases’ interaction Fig. 18. Predicted temperature with velocity vector fields close to two vortices for Tbot = 313.15 K, where Tliq is bounded for 333–336 K.

Fig. 19. Overview of flow mechanisms.

Figure 20(a) shows the predicted flame height and mass burning rate m˙ vs time, where (a)–(e) are also used for following Fig. 21. One cycle of oscillations of the mass burning rate and flame height are similar. This cycle is approximately 0.18 s, corresponding to (a)–(e). Their trends are slightly different, i.e. the flame height falls abruptly from (e)–(f), whereas such a change in the mass burning rate does not appear as shown in Fig. 20(a). As shown in Eqs. (11) and (12), m is a function of Tliq,inter ; m varies with time, meaning that Tliq,inter also changes with time. A significant contribution of heat transfer is convection, followed by radiation. Therefore, a flame motion such as a change in flame height would affect the mass burning rate. Further, in Fig. 20(b), a trend in the flame height between the coupled and decoupled approaches is close, and one-cycle period is similar. In other words, the oscillation of the mass burning rate does not affect well a trend in the flame height. Therefore, the flame motion preferentially influences the oscillation of the mass burning rate.

Fig. 20. Predicted (left) flame height xf and mass burning rate m˙ vs time t for Tbot = 297.15 K and (right) comparison of xf vs t using different approaches for Tbot = 268.15 K.

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

453

Fig. 21. Predicted transient flame volume and distribution of T for Tbot = 297.15 K.

bottom temperature is lower than the interface temperature, the reference case does not show the distinct buoyancy effect resulted from the density distribution. However, this density distribution inhibits thermal mixing with buoyancy. As a result, the higher interface temperature and mass burning rate appears for Case 1. Conversely, since thermal mixing is prominent under the zero-gravity condition, this effect reduces the mass burning rate for Case 6. 5. Conclusions

Fig. 22. Model sensitivity tests.

Figure 21 shows the predicted flame volume and temperature distribution between phases (a)–(e). In phase (a), the predicted flame height is the shortest, and this flame shape is similar to a ‘cone’. In phase (b), the flame near the pool surface becomes thinner than that in phase (a), and a necking phenomenon emerges. In phase (c), a mushroom cap shape appears with the highest flame height. After the flame reaches the highest point, the flame sprits up into tip and base parts in phase (d). The former disappears with convection and turbulent diffusion, and the latter becomes a similar shape to that of phase (a). These unsteady motions are close to those reported by Hayasaka [36]. 4.6. Effect of the different factors on the mass burning rate Finally, the validated model was used to investigate the different influencing factors in the mass burning rate. Case 1 was a standard condition with Tbot = 268.15 K for reference. Cases 2, 3, and 4 did not consider the effects of thermocapillary force, shear stress and pressure work term (Dpliq /Dt), respectively. Case 5 ignored heat transfer from the sidewall. For Case 6, g was zero in the liquid phase to examine the buoyancy effect. As shown in Fig. 22, the effects of thermocapillary force, shear stress, pressure work term, and heat transfer from the sidewall are all less than 5.1%. A significant difference in the predicted mass burning rate appears under the zero-gravity condition. Since the

The in-house version of FireFOAM was modified to implement a fully coupled 3D gas-liquid approach. It considered thermocapillary Marangoni convection, buoyancy, shear stress and evaporation in the liquid phase. The SGS modelling for the gas phase was the dynamic Smagorinsky model, and the combustion model was the modified eddy dissipation concept for LES with the treatment for laminar-turbulent transition. Gas-liquid phase coupled simulations were conducted for the small pool fire test of Vali et al. [1]. The predicted mass burning rate, flame height, their inter-dependence, and the liquid temperature distribution were all found to be in reasonably good agreement with the experimental measurements. The error of the mass burning rate was within 4%. The present predictions captured a pair of vortices observed by Vali et al. [1]. Their sizes increased with increasing Tbot following the experimental observations. The Reynolds analogy supported the rational reason behind this trend. Convection induced by shear stress and thermocapillary force arose a pair of vortices. The different distributions of convective and radiative heat fluxes yielded thermocapillary force. Numerical tests showed that the effects of the thermocapillary force, shear stress, pressure work term, and heat transfer from the sidewall on the mass burning rate were all less than 5.1%. Conversely, the simulation assuming zero gravity only in the liquid phase resulted in almost 64% reduction in the mass burning rate. Vali et al. also reported the mass burning rate and flame height for different boundary temperature [3], liquid pool depth [38], and wall materials [16]. These other scenarios also involved changes in the flow and thermal structures in the pool and could be a subject of future study. The validated model here holds potential applicability to the investigation into the liquid-phase transport phenomena, mass burning rates, and flame evolution of pool fires for different diameter sizes.

454

K. Fukumoto, J.X. Wen and M. Li et al. / Combustion and Flame 213 (2020) 441–454

Declaration of Competing Interest None. Acknowledgements FireFOAM 2.2.x was built and distributed by FM Global, and the authors acknowledge technical and help supports from them. The in-house version of FireFOAM used in this study firstly developed in the project funded by the National Key Research and Development Programme special for inter-governmental cooperation (No. 2016YFE0113400) and the European Commission FP7-PEOPLE2012-IIF (Grant number 328784). References [1] A. Vali, D.S. Nobes, L.W. Kostiuk, Transport phenomena within the liquid phase of a laboratory-scale circular methanol pool fire, Combust. Flame 161 (2014) 1076–1084. [2] E.J. Weckman, A.B. Strong, Experimental investigation of the turbulence structure of medium scale methanol pool fire, Combust. Flame 105 (1996) 245–266. [3] A. Vali, D.S. Nobes, L.W. Kostiuk, Effects of altering the liquid phase boundary conditions of methanol pool fires, Exp. Therm. Fluid Sci. 44 (2013) 786–791. [4] H. Sun, C.J. Wang, H. Liu, M. Li, A. Zhang, W. Zhao, C. Gao, Experimental study of combustion characteristics of circular ring thin-layer pool fire, Energ. Fuel 31 (2017) 10082–10092. [5] J.X. Wen, K. Kang, T. Donchev, J.M. Karwatzki, Validation of FDS for the prediction of medium-scale pool fires, Fire Safety J 42 (2007) 127–138. [6] Fire Dynamics Simulator (FDS) and Smokeview (SMV), https://pages.nist.gov/ fds-smv/ (Accessed 25 May 2019). [7] Y. Wang, P. Chatterjee, J.L. de Ris, Large eddy simulation of fire plumes, Proc. Combust. Inst. 33 (2011) 2473–2480. [8] FireFOAM 2.2.x, Source code, https://github.com/fireFoam-dev/fireFoam-2.2.x (Accessed 25 May 2019). [9] OpenFOAM, Source code and documentations, http://www.openfoam.com/ (Accessed 25 May 2019). [10] B.J. McCaffrey, Purely Buoyant Diffusion flames: Some Experimental Results, 20234, Centre for Fire Research, National Engineering Laboratory, National Bureau of Standards, Washington, D.C., 1979 NBSIR 79-1910. [11] G. Maragkos, T. Beji, B. Merci, Advances in modelling in CFD simulations of turbulent gaseous pool fires, Combust. Flame 181 (2017) 22–38. [12] Z. Chen, J.X. Wen, B. Xu, S. Dembele, Large eddy simulation of a medium-scale methanol pool fire using the extended eddy dissipation concept, Int. J. Heat Mass Trans. 70 (2014) 389–408. [13] K. Prasad, C. Li, K. Kailasanath, C. Ndubizu, R. Ananth, P.A. Tatem, Numerical modelling of methanol liquid pool fires, Combust. Theor. Model. 3 (1999) 743–768. [14] V. Novozhilov, H. Koseki, CFD prediction of pool fire burning rates and flame feedback, Combust. Sci. and Technol. 176 (2004) 1283–1307. [15] T. Sikanen, S. Hostikka, Modeling and simulation of liquid pool fires with indepth radiation absorption and heat transfer, Fire Safety J 80 (2016) 95–109.

[16] A. Vali, D.S. Nobes, L.W. Kostiuk, Quantifying the conduction pathways in a laboratory-scale methanol pool fire, Combust. Sci. Technol 187 (2015) 765–779. [17] C.J. Wang, J.X. Wen, Z.B. Chen, Simulation of large-scale LNG pool fires using firefoam, Combust. Sci. Technol. 86 (2014) 1632–1649. [18] K. Fukumoto, C.J. Wang, J.X. Wen, Large eddy simulation of upward flame spread on PMMA walls with a fully coupled fluid-solid approach, Combust. Flame 190 (2018) 365–387. [19] M. Germano, U. Piomelli, P. Moin, W.H. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A Fluid 3 (1991) 1760–1765. [20] D.K. Lilly, A proposed modification of the Germano subgrid: scale closure method, Phys Fluids A Fluid 4 (1992) 633–633. [21] K. Fukumoto, C.J. Wang, J.X. Wen, Large eddy simulation of a syngas jet flame: effect of preferential diffusion and detailed reaction mechanism, Energ. Fuel (2019), doi:10.1021/acs.energyfuels.9b00130. [22] A. Nakakuki, Heat transfer mechanisms in liquid pool fires, Fire Safety J. 23 (1994) 339–363. [23] T.F. Smith, Z.F. Shen, J.N. Frledman, Evaluation of coefficients for the weighted sum of gray gases model, J. Heat Transf. 104 (1982) 602–608. [24] I. Sikic, S. Dembele, J.J. Wen, Non-grey radiative heat transfer modelling in LES-CFD simulated methanol pool fires, J. Quant. Spectrosc. Ra. 234 (2019) 78–89. [25] W. Yang W, B. Wlodzimierz, Numerical study of fuel temperature influence on single gas jet combustion in highly preheated and oxygen deficient air, Energy 30 (2005) 385–398. [26] S.M. Ali, V. Raghavan, K. Velusamy, S. Tiwari, A numerical study of concurrent flame propagation over methanol pool surface, J. Heat Transf. 134 (2012) 1–9. [27] J. Warnatz, U. Maas, R.W. Dibble, Combustion, 4th ed, Springer, Berlin, 2006. [28] T. Beji, B. Merci, Development of a numerical model for liquid pool evaporation, Fire Saf. J 102 (2018) 48–58. [29] J. Chen, X. Zhang, Y. Zhao, Y. Bi, C. Li, S. Lu, Oxygen concentration effects on the burning behavior of small scale pool, Fuel 247 (2019) 378–385. [30] B.E. Poling, J.M. Prausnitz, J.P. O’Connel, The Properties of Gases and Liquids, 5th ed., McGraw-Hill, U.S., 2001. [31] M.D. Smooke, Reduced Kinetic Mechanisms and Asymptotic Approximations For Methane-Air flames, Lecture note in Physics, Springer-Verlag, Berlin, 1991. [32] K. Meredith, Y. Xin, J.D. Vries, A numerical model for simulation of thin-film water transport over solid fuel surfaces, 10th International Symposium Fire Safety Science (2011) 415–428. [33] Thermal-FluidsPedia, http://www.thermalfluidscentral.org/ (Accessed 23 May 2019). [34] E.N. Fuller, P.D. Schettler, J.C. Giddings, New method for prediction of binary gas-phase diffusion coefficients, Ind. Eng. Chem. 58 (5) (1966). [35] A.V. Singh, M.J. Gollner, A methodology for estimation of local heat fluxes in steady laminar boundary layer diffusion flames, Combust. Flame 162 (2015) 2214–2230. [36] H. Hayasaka, Radiative characteristics and flame structure of small-pool flames, Fire Technol 32 (1996) 308–322. [37] A. Vali, D.S. Nobes, L.W. Kostiuk, Characterization of flow field within the liquid phase of a small pool fire using particle image velocimetry technique, Exp Therm. Fluid Sci. 75 (2016) 228–234. [38] A. Vali, D.S. Nobes, L.W. Kostiuk, Fluid motion and energy transfer within burning liquid fuel pools of various thicknesses, Combust. Flame 162 (2015) 1477–1488. [39] H.F. Oztop, I. Dagtekin, Mixed convection in two-sided lid-driven differentially heated square cavity, Int. J Heat Mass Trans. 47 (2004) 1761–1769.