Numerical analysis of heat transfer characteristics of spray flames impinging on a wall under CI engine-like conditions

Design of Compression Ignition (CI) engines with improved thermal eﬃciencies needs better understanding of the heat transfer mechanism from spray ﬂame to the combustion chamber wall. In this regard, heat transfer occurring during the interaction between impinging spray ﬂame and wall, under CI engine-like conditions, is investigated in this study using 3-Dimensional numerical simulations based on an Eulerian– Lagrangian framework. Simulations are performed for different fuel spray injection velocities (which are representative of different fuel injection pressures in CI engines), to examine their inﬂuence on the heat transfer between impinging spray ﬂame and wall. To couple the convective and radiative heat transfer at the wall surface with the conduction heat transfer occurring within the ﬁnite thickness wall, Conjugate Heat Transfer (CHT) is incorporated in the simulations. A Non-Adiabatic Flamelet/Progress Variable (NA-FPV) approach is employed as the combustion model of n -dodecane, which is considered to be the fuel for liquid spray. Dynamics of the liquid ﬁlm formed on the wall surface by impinging fuel droplets are captured using a particle-based approach. Contribution of radiative heat ﬂux is taken into consideration using the Discrete Ordinates (DO) method. Results indicate that the total heat ﬂux (sum of convective and radiative heat ﬂuxes) at the wall surface increases with the fuel injection velocity. It is observed that the total wall heat ﬂux is largest in the stagnation zone where the spray ﬂame impinges directly on the wall surface, while the radiative heat ﬂux at the wall surface becomes larger as the distance from this stagnation zone increases. Additionally, it is found that the inﬂuence of fuel injection velocity on the radiative heat ﬂow rate at the wall surface is rather small. This radiative heat ﬂow rate when expressed as a percentage of the total wall heat ﬂow rate, ranges from ≈ 18% to 30% (depending on the 3 cases investigated), indicating that its contribution cannot be neglected for the CI engine-like conditions under which the present simulations are performed. Furthermore, to characterize the heat transfer occurring during spray ﬂame-wall interaction process, correlations between the Nusselt number Nu (corresponding to the wall heat loss) and Reynolds number Re (of the ﬂow ﬁeld) of the form Nu ∝ Re n , are analysed and compared with that of a previous experimental study to assess their applicability. It is found that, depending on how the Nusselt number Nu is deﬁned (either using the total wall heat ﬂux or the convective heat ﬂux), the value of the correlation index n changes. When Nu is calculated based on the total wall heat ﬂux (which includes the contribution from the radiative heat ﬂux), the value of n is found to be 0.49 which is close to the correlation index value of n = 0 . 4 reported in the recent experiments performed at Toyota Central R & D Labs., Inc.


Introduction
According to the United Nations Environment Programme (UNEP) annual emissions gap report [1] , global greenhouse gas emissions should be cut drastically in the next decade to limit the adverse consequences of climate change.Therefore, to mitigate the challenges of climate change, it is necessary to suppress the anthropogenic CO 2 emissions, most of which arise from the burning of fossil fuels.In this regard, many governments across the globe are expected to set ambitious new commitments to slash their CO 2 emissions, as part of the 2015 Paris Agreement [2] at the 21 st Conference of the Parties (COP21) to the United Nations Framework Convention on Climate Change (UNFCC).However, all energy outlooks [3,4] acknowledge that combustion will play an indispensable role in the energy sector for the foreseeable future.The road transportation sector which is vital for any economy, currently accounts for approximately 16% of the global CO 2 emissions [5] .Under these circumstances, further improvements in the thermal efficiencies of Compression Ignition (CI) engines (which are a type of widely-used Internal Combustion engine) are urgently needed.One way to achieve this, is to reduce the heat loss through the engine cylinder/combustion chamber wall, which accounts for nearly 30% of the fuel energy released in the form of combustion [6] .Consequently, good understanding of the wall heat loss mechanism under different combustion conditions such as engine size, fuel injection rate, etc., is paramount for developing fuel-efficient CI engines.If the heat loss through the combustion chamber wall could be accurately predicted in the early design stages, then the optimal combustion conditions for reducing heat loss can be devised, thereby achieving significant savings in the engine design cycle time.Normally, heat loss through the wall in Internal Combustion (IC) engines is estimated based on the principle of convective heat transfer, where the dimensionless average Nusselt number Nu (corresponding to the average heat flux at the wall surface assuming quasi-steady state conditions) is a function of two other dimensionless parameters, viz. the Reynolds number Re and the Prandtl number P r.In this regard, since the 20th century, many researchers have performed IC engine experiments, using both Spark Ignition (SI) and CI engines [7][8][9][10][11][12] , and proposed empirical equations/correlations (derived from their experimental data) of the form: Nu ∝ Re n (1) under the assumption that the Prandtl number P r is constant.In Eq. (1) , n is the correlation index.Once such a correlation for Nu is known, it is possible to directly calculate the value of Nu , and then using this knowledge of Nu , the average convective heat transfer coefficient h heat can be estimated.The average wall heat flux q w can then be computed easily from Newton's law of cooling, which is expressed as: where T is the representative gas-phase temperature and T w is the wall temperature.To this extent, Woschni [10] applied the analogy of non-reacting fully-developed turbulent flow in pipes (i.e., forced convection assumption) to the flow in engines, and reported a correlation between Nu and Re of the form Nu ∝ Re 0 . 8, similar to Eq. (1) .Therefore, Woschni determined the value of the correlation index n in Eq. ( 1) to be 0.8.However, the applicability of Woschni's correlation during the combustion period inside an IC engine has not been verified, even though it is still being used for zero-dimensional cycle simulations [e.g., 13 , 14] .In the same manner, other researchers have also proposed values for the correlation index n in Eq. (1) based on their experiments (using forced convection assumption), such as Oguri ( n = 0 .5 ) [8] , Annand ( n = 0 .7 ) [9] , Sitkei and Ramanaiah ( n = 0 .7 ) [11] , and Hohenberg [12] ( n = 0 .8 ).
Recently, Inagaki et al. [15] investigated the effects of bore size and fuel injection conditions on the heat loss through combustion chamber wall, experimentally and theoretically from the viewpoint of fuel efficiency.Follow-up experiment performed by Inagaki et al. [16] aimed to propose a correlation for the heat transfer coefficient, by applying the similarity of spray characteristics to diesel spray combustion, and tried to elucidate the dependence of heat loss through the wall (during spray flame-wall interaction process) on the representative velocity.From their experiments [15,16] , a power law dependence of Nu on Re of the form Nu ∝ Re 0 . 4was deduced, indicating the value of the correlation index to be n = 0 .4 .
However, in a more recent experimental study conducted by Mahmud et al. [17] , a correlation index value of n = 0 .8 was found (i.e., Nu ∝ Re 0 . 8).Hence, it is evident from all the above-mentioned experimental investigations that, there is an appreciable variation in the value of the correlation index n , indicating the lack of a clear consensus on the relation between Nu and Re corresponding to wall heat loss during spray flame-wall interaction process.Although measurement errors might occur while performing experiments, the variations in the value of correlation index n proposed in the aforementioned studies [8][9][10][11][12][15][16][17] , could also arise due to the way in which Nu and Re are defined in the respective studies.This is owing to the physical constraints and restrictions of a given experimental setup, which could make it very challenging to take into account the detailed spatial distributions and temporal evolution of flame temperature, spray flame-wall contact area, gasphase velocity and wall heat flux, while performing measurements.Another reason could be the fact that, when wall heat flux is measured in an experiment, it might implicitly include the contribution of radiative heat transfer.Under CI engine-like conditions, the contribution of radiative heat transfer from spray flame to the wall heat flux could be non-negligible, and can account for 5% to 15% of the total wall heat flux as reported by Yan and Borman [18] .So, when such heat flux measurements are used to calculate the Nusselt number Nu (a parameter fundamentally defined in the context of convective heat transfer alone), significant variations in the value of n could occur.For these reasons, numerical simulations can prove to be effective tools for elucidating such complex two-phase reacting flow problems.Considerable progress has been made on numerical simulation and modelling of two-phase flows over the decades, from the earlier studies by Bray and co-workers [19][20][21] to the more recent ones [e.g., 22 , [23][24][25][26][27][28][29][30][31][32][33][34][35][36][37] , 38] .Modelling and simulation of spray combustion under Diesel engine-like conditions have also been performed [e.g., 39 , [40][41][42][43][44][45][46] , 47] .For premixed flames impinging on walls and interacting with walls, there exist some early theoretical [e.g., 4 8 , 4 9 , 50] and experimental [e.g., 51 , 52 , 53] studies by Bray and co-workers, and several numerical investigations by other researchers [e.g., 54 , [55][56][57][58] , 59] .However, when it comes to turbulent spray combustion, it is an intricate two-phase phenomenon in itself, involving turbulent mixing, evaporation of the liquid fuel droplets, their dispersion and convection by the carrier gas-phase and combustion of the evaporated fuel, all of which occur simultaneously while influencing each other.Furthermore, the simulation of a spray flame impinging on a wall adds another level of complexity (due to fuel droplets impacting on the wall surface and forming a liquid film on it, and flame-wall interaction), is difficult to model and is demanding in terms of computational resources.Hence, studies involving numerical simulations of spray flames impinging on walls are scarce to the best of the authors' knowledge, with one example being that of Hori et al. [45] .In their study, Hori et al. [45] numerically analysed Diesel spray flames impinging on a flat wall using Reynolds Averaged Navier-Stokes (RANS) simulations, to investigate the relationship between fuel injection pressure and wall heat loss.But, they did not take into account the temporal and spatial variations of temperature inside the wall, nor did they consider the influence of radiative heat transfer on wall heat loss.
In this study, 3-Dimensional numerical simulations of spray flames impinging on a flat wall under CI engine like conditions, are used to analyse the influence of fuel injection velocity (and hence the fuel injection rate) on the wall heat flux (convective, radiative and the sum of both these heat fluxes) and the wall heat flow rate.The differences between simulations in this study and those of previous studies are that, radiative heat transfer and the dynamics of liquid fuel film formed on the wall surface (due to impinging fuel droplets) are taken into account.Additionally, Conju-gate Heat Transfer (CHT) is also considered, wherein the heat conduction occurring within the finite thickness wall is coupled with the convective and radiative heat transfer at the wall surface, to track the temperature variations inside the wall.To the best of the authors' knowledge, such a comprehensive analysis of spray flames impinging on a wall in CI engine-like environment, using numerical simulations has not yet been performed.

Numerical method and model formulation
In the present numerical simulations of spray flames impinging on a wall, neither subgrid-scale turbulence model nor turbulencechemistry interaction model is used, since fine grids with coarse DNS-like resolution are employed in each simulation (discussed later).Simulations of the two-phase reacting flows are based on an Eulerian-Lagrangian framework [22,24,26,29,32,60] in which, the continuous gas-phase is modelled in the classical Eulerian coordinate framework, while the dispersed-phase fuel droplets are tracked individually in the Lagrangian framework as point-masses.Additionally, liquid film formation on the wall surface due to fuel droplet impingement, is captured using a particle-based approach (Lagrangian framework).As mentioned previously, CHT is incorporated in these simulations to couple the spray combustion calculations with the conduction heat transfer occurring inside the solid wall.Details of the various numerical approaches mentioned above, along with information about the models for liquid fuel evaporation, combustion reaction, soot formation and radiation are presented in the following.

Governing equations of gas-phase
Conservation equations of mass, momentum, enthalpy, mixture fraction and reaction progress variable (RPV) are solved for the gasphase, which is treated as an Eulerian continuum, in the following manner: In Eqs.
(3) - (7) , ρ is the density, u is the gas-phase velocity vector, τ is the viscous stress tensor, P is the pressure, h is the specific enthalpy, Z is the mixture fraction, C is the reaction progress variable (RPV) defined as the sum of mass fractions of combustion products, i.e.
where c p and λ are the specific heat capacity at constant pressure and thermal conductivity, respectively, of the gas-phase.Q rad is the radiative heat source term whose calculation method is described later.Source terms S ρ , S ρu , S ρh and S ρZ appearing in Eqs.(3) -( 6) , take into account the interactions between the gas-phase and the liquid-phase, allowing for two-way coupling between both the phases.These source terms are expressed as: S ρZ = S ρZ,d + S ρZ,l (12) where terms with the subscript d (i.e., S ρ,d , S ρu ,d , S ρh,d and S ρZ,d ) represent the interactions between the dispersed-phase (contribution from each individual fuel droplet) and gas-phase, while the terms with subscript l (i.e., S ρ,l , S ρu ,l , S ρh,l and S ρZ,l ) represent the interactions between the liquid fuel film formed on the wall surface and the gas-phase.

Combustion model
Combustion of vaporized n -dodecane (fuel used in this study) is modelled using the Non-Adiabatic Flamelet/Progress Variable (NA-FPV) approach [61,62] .NA-FPV is a type of flamelet approach in which the influence of heat loss due to latent heat of liquid fuel evaporation, cooled wall and radiative heat transfer, on the variation of each physical quantity is accounted for.With the NA-FPV approach, the conservation equations for mass, momentum, enthalpy, mixture fraction and RPV, i.e., Eqs.(3) - (7) , are solved in the physical space, while the mass fractions of various chemical species Y k and the reaction rate of C ( ˙ ω C ) at each grid point in the gas-phase's computational domain, are extracted from a precomputed database called the flamelet library.A reduced set of control variables, viz.mixture fraction Z, reaction progress variable C and enthalpy loss h are used for extracting Y k and ˙ ω C from the flamelet library.Here, h is the difference between the enthalpy h 0 when there is no heat loss (adiabatic condition) and the enthalpy h in Eq. ( 5) when heat loss is considered.Contrary to the conventional FPV approach, the NA-FPV approach stores flamelet data with heat loss in the flamelet library, which enables accurate prediction of chemical reactions in non-adiabatic combustion fields involving latent heat of evaporation of fuel droplets and liquid fuel film, heat loss through the wall, and heat loss due to radiative heat transfer from soot and gas-phase species.The diffusion flamelet database (flamelet library) is generated using FlameMaster [63] by solving the flamelet equations for 1-Dimensional counter diffusion flame, using a detailed chemical kinetic mechanism of n -dodecane consisting of 255 chemical species and 2289 reactions [64] .Details of the procedure for generating the flamelet library are described in Moriai et al. [61] , Kishimoto et al. [62] .

Governing equations of dispersed-phase fuel droplets
Dynamics of the dispersed-phase fuel droplets (treated as spherical point masses) are described by a set of Lagrangian equations governing the evolution of the fuel spray.Evaporating fuel droplets of the dispersed-phase are individually tracked using a Lagrangian framework [22,29,60,65] by solving the equations for droplet trajectory x d , velocity u d , temperature T d , and mass m d in the following manner: ) Here, c p,d is the specific heat capacity of fuel droplet, τ d is the droplet response time in Stokes regime, g is the acceleration due to gravity, u , T , Pr g and Sc g are the local gas-phase properties, viz.gas-phase velocity, temperature, Prandtl number and Schmidt number, respectively, at the droplet position.Similarly, Nu g and Sh g are the local gas-phase Nusselt and Sherwood numbers, respectively, at the droplet position, which are calculated using the Ranz-Marshall correlations [66,67] .L V is the latent heat of vaporization at the droplet temperature T d , and B M is the Spalding mass transfer number.f 1 and f 2 are the correction factors for Stokes drag and interphase heat transfer of evaporating fuel droplets, respectively [29,60,68,69] .A non-equilibrium Langmuir-Knudsen evaporation model [60,65,70,71] is used to capture the evaporation of dispersed-phase fuel droplets in Eq. ( 16) .Details of this evaporation model and the calculation procedures for the above-mentioned quantities appearing in Eqs. ( 14) -( 16) are described elsewhere [22,26,29,68,69] , and interested readers are directed to these studies for further information.By solving the system of Lagrangian Eqs. ( 13) -( 16) , the source terms accounting for the interactions between dispersed-phase fuel droplets and the gas-phase (i.e., S ρ,d , S ρu ,d , S ρh,d and S ρZ,d ) in Eqs. ( 9) - (12) , can be computed using the Particle-Source-In-Cell (PSI-Cell) approach [72] as follows: where V is the volume of the computational grid-cell within which a certain fuel droplet resides, h d is the specific enthalpy of a fuel droplet, and N is the number of fuel droplets within a given computational grid-cell.In the classical PSI-Cell approach, the contribution of a fuel droplet to the above source terms is interpolated to each of the 8 neighbouring grid points of the gas-phase's computational cell (i.e., when the geometry of the Cartesian grid-cell is cuboidal) within which that fuel droplet resides.The interpolation of each droplet's source term contribution is achieved using the conventional trilinear interpolation technique, considering the relative positions of the droplets within a computational cell.However, this simple trilinear interpolation technique can be applied only when the point mass assumption for a droplet is true, meaning the grid-cell size must be about 10 times larger than the droplet size [73,74] .For the computational grid used in the present simulations, the grid-cell size is of the order of droplet sizes or even smaller than the droplet sizes in certain regions of the computational domain (especially in the regions closer to the wall).Applying the conventional trilinear interpolation technique to the droplet source terms in such regions would lead to erroneous estimates of droplet dynamics, because with this interpolation technique, the range of influence of the evaporating fuel droplets on the gas-phase depends on the grid-cell size.To overcome this drawback, a Gaussian distance function technique proposed by Borghesi et al. [75] is used for interpolating the droplet source term contributions, in the regions where the grid-cell sizes are equal to or smaller than twice the initial Sauter Mean Diameter (SMD) of the injected fuel spray.This technique allows for the range of influence of the fuel droplets on the gas-phase to be kept constant regardless of the computational cell size.Although the Gaussian distance function interpolation technique is computationally more demanding than the trilinear interpolation technique, it has been used in the above-mentioned regions of the computational domain to capture droplet dynamics accurately.Similar techniques for interpolating the dispersed-phase source terms to the gas-phase have already been applied successfully by other researchers [76,77] .Details of the implementation of the Gaussian distance function technique are presented in Haruki et al. [27] , Borghesi et al. [75] .

Spray-wall interaction modelling
When liquid fuel is sprayed directly into the cylinder of a CI engine, there is a possibility of some of the fuel spray impinging on the piston and cylinder walls (before vaporization and mixing are finished), resulting in the formation of liquid films.In this study, the spray-wall interaction phenomena is modelled using a combination of the spray-film interaction model proposed by Stanton and Rutland [78] , and the wall film dynamics model proposed by O'Rourke and Amsden [79,80] , which are briefly introduced below.

Spray-film interaction model
The behaviour of fuel droplets impinging on a wetted wall surface can be classified into 4 regimes, namely: Stick, Rebound, Spread and Splash, depending on the Weber number W e of the incoming dispersed-phase droplets [78] , which is defined as: (21) where, ρ d is the droplet density, u d,n is the incoming droplet's velocity component perpendicular to the wall before impingement, and γ d is the liquid droplet's surface tension.It is assumed that the solid wall surface is covered by a thin film of liquid fuel when the droplet impact occurs.The transition criteria for each regime is decided by the Weber number W e , which is calculated using the incident dispersed-phase droplet's wall-normal velocity component u d,n , as shown in Eq. (21) .As suggested by the name of each regime, the outcome of the collision of a droplet with the film surface can be easily discerned.Droplets undergoing Sticking, Spreading and Splashing phenomena usually contribute to the formation and growth of liquid film on the wall surface, and further details of each impingement regime are discussed in Stanton and Rutland [78] .

Wall film dynamics model
The wall film dynamics model proposed by O'Rourke and Amsden [79,80] is a particle-based numerical method used to model the dynamics and evaporation of liquid fuel films formed on the walls of IC engines.Dominant physical phenomena to be considered when tracking the motion of liquid film on the wall surface are summarized in Fig. 1 .Using the above assumptions and setting the wall velocity u w = 0 (since the wall is kept stationary in the present simulations), the governing equations for the conservation of mass, momentum and energy of the liquid film formulated by O'Rourke and Amsden [79,80] , can be re-written as: where ρ l is the liquid film's density (temperature-dependent), h l is the film thickness, u l is the mean film velocity, T l is the mean film temperature, T s is the gas-film interface temperature, T w is the wall surface temperature, τ w is the shear stress acting on the gas-side of the wall film (at the gas-film interface), and t is the unit tangent to the wall in the direction of u l .μ l (T l ) , λ l (T l ) and I l (T l ) are the liquid film viscosity, thermal conductivity, and internal energy, respectively, at the mean film temperature T l , and g is the acceleration due to gravity.˙ M imp is the liquid mass source per unit area due to impingement of dispersed-phase fuel droplets, ˙ M v ap is the film evaporation rate/area, P l is the film pressure arising from the impingement of fuel droplets, ˙ P imp is the momentum source/area due to collisions of the fuel droplets with the film, n is the unit vector normal to the wall surface, and δ p l is the pressure gradient across the film (in the direction perpendicular to the wall) to satisfy the condition of u l • n = 0 .Finally, c v l is the liquid specific heat capacity at constant volume (temperature-dependent), and ˙ Q imp is the energy source term due to impingement of fuel droplets.
In this model, the particle numerical method [79,80] is used to track the motion of the wall film, wherein the wall film is assumed to be constituted of several discrete particles as illustrated by the schematic diagram in Fig. 2 .In this particle method, the dispersed-phase fuel droplets impinging on the solid wall/liquid film and undergoing either Stick, Spread or Splash phenomena, are converted into wall film particles.α denotes the computational cell face on the wall surface (see Fig. 2 ) over which the discrete film particles reside.Calculation method of the wall film source terms corresponding to mass, momentum, and energy, viz.˙ M imp , ˙ P imp and ˙ Q imp , respectively, which arise from the impingement of fuel droplets on the wall film surface, over any given computational cell face α, are detailed in O'Rourke and Amsden [79] , 80 ].
However, there is an important distinction between the numerical solution procedure applied by O'Rourke and Amsden [79] , and that applied in our simulations pertaining to: (1) the gas-side heat flux to the liquid film ˙ Q gas , (2) the gas shear stress acting on the gas-film interface τ w , and (3) the film particle evaporation rate ( ˙ M v ap ) p (subscript p denotes a wall film particle).In the RANS simulations performed by O'Rourke and Amsden [79] , τ w , ˙ Q gas and ( ˙ M v ap ) p are computed using wall functions developed in the context of the K/ turbulence model.However, in this study, τ w and ˙ Q gas are calculated directly from the real-time solutions of the numerical simulations, i.e., using the velocity and temperature gradients, respectively, on the gas-side of the gas-film interface, which can be readily computed in the simulations, since fine grid spacing in the wall-normal direction with DNS-like resolution is used.On the other hand, the film particle evaporation rate ( ˙ M v ap ) p is calculated in this study as [27,81] : where ρ is the gas-phase density above the film surface, D F is the diffusion coefficient of the fuel vapour calculated using Eq. ( 8) , n α is the unit vector normal to the wall surface on cell face α, and (Y F,g ) p is the equilibrium fuel vapour mass fraction at the gas-film interface calculated as: where (X F,g ) p is the fuel vapour mole fraction at the gas-film interface, P sat,p is the saturated vapour pressure corresponding to the gas-film interface temperature T s,p , and W F and W are the molecular mass of fuel vapour and the gas-phase's average molecular mass, respectively.In this study, the liquid film evaporation rate given in Eq. ( 25) is modelled under the assumption of vapourliquid equilibrium at the gas-film interface.The physical meaning of Eq. ( 25) is that, film vaporization occurs due to the difference between the mass fraction of fuel vapour at the gas-film interface and in the bulk gas above the film surface.
For further details pertaining to the wall film dynamics model, its sub-models and numerical implementation, interested readers are directed to the original works of O'Rourke and Amsden [79,80] .Finally, using the solution of the wall film dynamics model, closures for the source terms required to couple the wall film with the gas-phase, i.e., S ρ,l , S ρu ,l , S ρh,l and S ρZ,l , appearing in Eqs. ( 9) -( 12) can be computed as: where N p,α is the number of film particles residing on the cell face α (on the wall surface), and m p , u p and h p are the mass, velocity and specific enthalpy of wall film particle p, respectively.The contribution of each wall film particle to the above source terms used for phase coupling, is then interpolated to the grid points of the gas-phase above the film surface, using the Gaussian distance function technique mentioned in Section 2.2 .

Radiation and soot formation models
Models used in the present numerical simulations to account for radiation heat transfer and soot formation, are briefly described here, since extensive literature on these models already exist.For the computation of radiation heat transfer, the Radiative Transfer Equation (RTE) shown in Eq. ( 32) below, is solved by using the Discrete Ordinates (DO) method with S8 quadrature set [82] .
Here, I λ ( r , s ) is the spectral radiation intensity, which is a function of the position vector r and the direction vector s , I bλ ( r ) is the blackbody spectral intensity at the temperature of the medium, and κ λ is the spectral Planck mean absorption coefficient.By integrating I λ ( r , s ) obtained from Eq. ( 32) over wavelength λ and solid angle , the net radiant power per unit volume ( ρQ rad ) can be calculated as shown below: where Q rad is the radiative heat source term appearing in the gasphase enthalpy conservation Eq. ( 5) .Thus, the RTE is coupled with the gas-phase enthalpy conservation equation through Q rad .The radiative wall heat flux q w,rad is calculated as shown in Eq. ( 34) below: where, n is the wall-normal vector, and I absorbed and I emitted are the intensities of radiation absorbed and emitted by the wall surface.Dominant radiating entities considered in this study are soot and gas-phase species viz.CO, CO 2 , and H 2 O, all of which are assumed to be gray (Planck mean absorption coefficient of each gas-phase species is calculated based on the data from HITEMP 2010 highresolution database [83] ) to keep computational costs realizable.The calculation methodology for the Planck mean absorption coefficients of gas-phase species (CO, CO 2 , and H 2 O) is described in Chmielewski and Gieras [84] , and that for the Planck mean absorption coefficient of soot is described in Byun and Baek [85] .
Soot formation is predicted using the Moss-Brookes model [86,87] in the present simulations, in order to account for the radiative heat transfer from soot.In this model, the transport equations for the soot number density and mass density are solved for predicting soot formation due to inception, coagulation, surface growth, and oxidation as described in Haruki et al. [27] , Watanabe et al. [87] , Muto et al. [88] .Precursor models based on Polycyclic Aromatic Hydrocarbons (PAH) [89] are employed for the computation of inception process.Mass fractions of C 2 H 2 , C 6 H 5 , C 6 H 6 , H 2 and OH, required for the computations of different processes in the soot model are obtained from the flamelet library [27,87] .Although this soot formation model cannot provide the soot particle size distribution, it is computationally less expensive than the other more sophisticated models.

Implementation of Conjugate Heat Transfer (CHT)
To couple the heat transfer occurring at the wall surface from the spray flame, liquid film formed on the wall surface, and radiation, with the unsteady heat conduction occurring inside the stationary solid wall, it is necessary to solve the 3-Dimensional transient heat conduction equation shown in Eq. ( 35) below, along with the application of appropriate boundary conditions at the gaswall and film-wall interfaces.Here, ρ s , c s and λ s are the temperature-dependent density, specific heat capacity and thermal conductivity, respectively, of the wall material (iron), and T s is the temperature inside the wall.By solving the above heat conduction equation, the temporal and spatial variations of temperature inside the wall can be computed.Total heat flux at the wall surface q w,total is calculated as the sum of the convective heat flux q w,con v and the radiative heat flux q w,rad as: q w,total = q w,con v + q w,rad (36) where the convective heat flux at the wall surface q w,con v is calculated using Eq. ( 37) based on Fourier's law, and the radiative heat flux at the wall surface q w,rad is calculated from Eq. ( 34) using the Discrete Ordinates method [82] .

Interface boundary conditions
The schematic diagram in Fig. 3 illustrates the three phases (gas, liquid film and solid wall) among which heat transfer occurs.For coupling the heat exchange between the liquid film and the wall, the heat flux balance equation at the liquid film-wall interface must be solved implicitly, under the constraint that the liquid film temperature T l and the wall temperature T s at the film-wall interface are equal, as shown below: In the above equation, T l,c represents the mean film temperature at the center of the film, and T l,w represents the liquid film temperature at the wall surface.A similar heat flux balance equation must be solved implicitly for the gas-wall interface, under the constraint of equal gas-phase temperature T g and wall temperature T s at the gas-wall interface, in order to couple the heat exchange between the gas-phase and wall: + q w,rad = −λ s ∂T ∂x s ( at the wall surface ) T g = T s (39) In Eqs. ( 38) and ( 39) above, the subscripts l, g and s denote the physical quantities of the liquid fuel film, the gas-phase and the solid wall, respectively.The liquid film is thin enough to assume that its transmissivity is 1.0, and therefore, all the radiation penetrates the liquid film and reaches the wall.

Computational configuration and solution algorithm
Schematic of the computational domain used in the present simulations is depicted in Fig. 4 .The domain's dimensions are the same as those of the combustion chamber apparatus used in a previous experimental study by Inagaki et al. [15] .For the spray flame calculations, the computational domain measures 25 mm × 65 mm × 75 mm in the x -, y -and z-directions, respectively.The spray nozzle is located at the top as shown in the figure, and the point on the nozzle axis at the nozzle exit plane is the origin of the computational domain.The spray flame simulations are performed on a non-uniform staggered Cartesian grid consisting of 720 × 560 × 560 grid points in the x -, y -and z-directions, respectively, so the total number of grid points used for the gas-phase in each simulation case is approximately 226 Million.For the spray flame computations, minimum grid spacing of x = 2.9 μm, y = z = 40 μm are used, and the grid resolution is adequate in the core regions of the jet spray flames, near the wall and other regions where combustion occurs (discussed in Section 3.1.2).In the experiment [15] , firstly, premixed Hydrogen-Air mixture was charged into the combustion chamber and ignited using a spark plug.Subsequently, the flame propagated throughout the chamber creating a high pressure and temperature gas.The resulting ambient pressure and temperature inside the chamber were 3.0 MPa and 1600 K, respectively.Hence, ambient gas temperature and pressure of 1600 K and 3.0 MPa, respectively, are used in the simulations.The ambient gas initially consists of N 2 , O 2 , and H 2 O, whose volume fractions are 0.60, 0.21, and 0.19, respectively.
Dimensions of the solid wall's computational domain are 4 mm × 65 mm × 75 mm in the x -, y -and z-directions, respectively.This means that the wall thickness is 4 mm and it is uniformly discretized using 400 grid points in the wall-normal x -direction, yielding a grid spacing of x wall = 10 μm.The grid spacing and the number of grid points in the y -and z-directions are the same as those of the spray flame's computational domain (i.e., 560 grid points each in the y -and z-directions).No-slip velocity boundary condition is imposed on the upper wall surface which is exposed to the spray flame and the liquid fuel film.The initial temperature distribution inside the wall varies linearly in the x -direction, from 481 K at the upper surface (flame side) to 463 K at the bottom surface as depicted in Fig. 5 .The bottom surface of the wall is kept isothermal at 463 K throughout the duration of each simulation (similar to the experiment), and the wall surface emissivity ε is assumed to be unity ( ε = 1 .0 ).In Fig. 4 , the stagnation point "O" on the upper surface of the wall, refers to the point where the nozzle/spray axis intersects the wall surface (i.e., radial location r = 0 ).
At the start of each simulation, n -dodecane fuel droplets with an initial temperature of 330 K are injected into the computational domain from the nozzle exit plane at the top of the domain.The nozzle diameter is 0.083 mm, which is same as that of the experiment.Table 1 summarizes the cases simulated in this  study in terms of their fuel injection parameters.In order to investigate the influence of fuel injection velocity U in j on the wall heat loss, three cases with U in j values of 100 m/s, 150 m/s and 216 m/s are considered.The injection Mach number M in j based on the fuel injection velocity U in j and the speed of sound in ambient gas C amb ≈ 795 m/s, for each case is listed in Table 1 .The maximum Mach number observed in each case's simulation is also provided in Table 1 .The time variations of injected fuel mass in all three cases are depicted in Fig. 6 , where the slope of each line corresponds to the fuel injection rate listed in Table 1 .In all three cases, the simulations are performed until 0.64 mg of fuel is injected, as shown in Fig. 6 .The spray cone angle used for injecting the fuel droplets into the computational domain in each case, is calculated using the following empirical formula obtained from experiments [90] : where c = 0 .00121 , D n is the nozzle diameter in meters ( D n = 0 .083 × 10 −3 m), L n is the nozzle length in meters ( L n = 0 .68 × 10 −3 m), P a is the ambient pressure in Pascals, P in j is the fuel injection pressure in Pascals, μ a is the ambient gas-phase dynamic viscosity [ Pa s ], and ρ a is the ambient gas-phase density [kg/m 3 ].The fuel injection pressure P in j for a given fuel injection rate (or fuel injection velocity U in j ) can be calculated from Hiroyasu's equation [91] : U in j = 0 . 39 2(P in j − P a ) /ρ d (41) Distributions of the fuel droplet diameter d are prescribed using the Rosin-Rammler [92] volume-based Probability Density Function (PDF) f v (d) , shown in Eq. ( 42) below.These size distribution PDFs for the fuel sprays that are injected from the nozzle exit, are obtained from the best-fit curves to the experimental data In Eq. ( 42) , N and X r are the function parameters obtained from experimental measurements [15] .The minimum and maximum injected droplet diameters are d min = 1 μm and d max = 80 μm, respectively.Fig. 7 shows both the volume-based and number-based droplet size distributions of the injected fuel spray, corresponding to the different fuel injection velocities U in j of the three cases.It is evident from Fig. 7 (a) that, the droplet size distribution becomes biased towards the smaller droplet diameters as the fuel injection velocity U in j increases, and hence the reduction in the SMD of the injected fuel droplets with increasing U in j in Table 1 .
The fuel injection parameters listed in Table 1 are the same as those in the experiments of Inagaki et al. [15] .However, it can be seen that when the fuel injection velocity is changed, other parameters such as fuel injection rate, spray cone angle and SMD also change.This is a limitation of the fuel injector used in the experiment, wherein for a fuel injector of a given nozzle diameter, varying the fuel injection velocity will inevitably change the other parameters.So, it is very difficult to isolate the effect of fuel injection velocity alone, in the experiment.Main factors governing the heat transfer occurring during the interaction between impinging spray flame and wall, are the flame temperature and the velocity with which the flame impinges on the wall.The fuel injection velocity directly influences the flames impingement velocity, and therefore it is obviously an important parameter.On the other hand, spray cone angle and SMD do not have strong effects on the impingement velocity of the flame.However, we do acknowledge that the fuel injection rate and SMD affect the heat release rate, and hence the flame temperature.But, as will be shown later, the fuel injection velocity does not influence the flame temperature appreciably in the near wall region, in the present simulations.In the experiments, only the fuel injection velocity was varied and its effect on heat transfer to the wall was investigated.The simulations in this study are mimicking the experiments and therefore, we have only varied the fuel injection velocity and the other parameters which inevitably change as a consequence of varying the fuel injection velocity, have been matched with those of the experiments [15] and used as it is.
The simulations are performed using an in-house thermal flow analysis code called FK 3 [93] .The FK 3 code employs a pressurebased semi-implicit solver for compressible flows, whose algorithm consists of a fractional-step method [94] .The FK 3 code has been successfully applied for simulating various types of gaseous and spray combustion problems, and has also been extensively validated [93] .In the present simulations, the spatial derivatives of the convective terms in the momentum equation are evaluated using a 4 th order accurate finite difference scheme [95] that conserves kinetic energy locally, and the spatial derivatives of the convective terms in the governing equations of the scalar quantities are evaluated using the 5 th order WENO (Weighted Essentially Non-Oscillatory) scheme [96] .Time integration of the convective terms is performed using a 3 rd order explicit TVD Runge-Kutta scheme.Gas-phase thermophysical properties ( c p , c v , λ and μ) taking temperature dependence into account are computed according to CHEMKIN [97,98] .Thermopysical properties of liquid fuel (i.e., dispersed-phase droplets and wall film particles) are computed using temperature polynomial functions generated from the NIST database [99] .The solid wall's material is iron and its thermophysical properties, viz.density, specific heat capacity and thermal conductivity, are computed using temperature polynomial functions generated from the database of Valencia et al. [100] .Each simulation was performed using 1600 cores on the HPE SGI8600 supercomputer at the Central Research Institute of Electric Power Industry (CRIEPI), Japan.The wall clock time and CPU hours for each case are listed in Table 1 .

Combustion field
Analyses of the results pertaining to the spray combustion fields, obtained from the numerical simulations are now presented.Simulation results of Case 3 are also compared with available experimental measurements.

Features of wall impinging spray flames
Instantaneous distributions of gas-phase temperature in the central x − y plane and at progressing time instances, are shown for all the three cases in Figs.8-10 .Each time instance corresponds to a specific value of injected fuel mass, since (injected fuel mass) = (fuel injection rate) × (elapsed time).The white entities in these figures represent the dispersed-phase fuel droplets and the particles belonging to the liquid fuel film formed on the wall surface.It is evident from these figures that the gas-phase temperature in the spray jets' core regions, where majority of the fuel droplets exist, is relatively low owing to the latent heat of fuel droplet evaporation.Injected fuel droplets evaporate as they flow downstream and the evaporated fuel mixes with the surrounding hot ambient gas, and gets ignited by it.Subsequent combustion of the fuel vapours occurs as is evident from the high temperature regions in the figures.With increasing downstream distance from the nozzle exit, fuel droplet size decreases and a considerable reduction in the number density of fuel droplets is also observed, as a consequence of evaporation and combustion.In all the three cases, the vaporized fuel auto-ignited at some distance away from the spray nozzle exit, and no special treatment was applied for igniting the vaporized fuel.As time advances, it can be confirmed that the spray flames impinge on the flat wall surface at time instances t ≈ 0 .38 ms, t ≈ 0 .27 ms and t ≈ 0 .19 ms, respectively, in Case 1, Case 2 and Case 3.These time instances correspond to an injected fuel mass value of m in j = 0 .15 mg ∼ 0.16 mg, among the three cases.After impinging on the wall, the flames start spreading in the direction parallel to the wall surface.
Since no significant qualitative differences can be discerned visually, among the combustion fields of the three cases in Figs.8-10 , the combustion field of Case 3 ( U in j = 216 m/s) alone will be analysed in the following for the sake of brevity.Instantaneous distributions of the iso-surfaces of gas-phase temperature at 1700 K (coloured in yellow) and 2200 K (coloured in red) are shown in Fig. 11 , along with the dispersed-phase fuel droplets (coloured entities) which are coloured by their respective evaporation rate ˙ m d , for Case 3. The distributions are shown at the time instance when injected fuel mass m in j = 0 .64 mg.The red iso-surface for 2200 K envelopes the yellow iso-surface for 1700 K in most regions, which is why the overall colour of the temperature iso-surfaces appears to be orange.The figure reveals that the evaporation rate of the fuel droplets becomes higher in the streamwise direction as they move into the hotter regions of the combustion field.For the time instance at which these temperature iso-surfaces are shown, the spray flame has already impinged on the wall surface and started m d , and gas-phase temperature T .Radial distributions of all these quantities are plotted for the streamwise location of x = 15 mm from the nozzle exit, and at the time instance corresponding to m in j = 0 .64 mg [see Fig. 10 (d)].Here, "spatial averaging" implies averaging performed along the circumference of a circle of a given radius r, whose center is situated at r = 0 (i.e., at the axis of the nozzle/spray jet).Fig. 12 (a) shows that the mass fraction of gaseous fuel Y C 12 H 26 is higher in the central region of the spray jet, due to the high evaporation rate of fuel droplets ˙ m d in that region, and decreases towards the ambient.On the other hand, the mass fraction of oxidizer O 2 (i.e., Y O 2 ) is larger outside the jet, and decreases to a small value towards the axis of the spray jet (i.e., at r = 0 ).Therefore, the radial distributions of the mass fractions of both fuel ( Y C 12 H 26 ) and oxidizer ( Y O 2 ), decrease towards the region of the mixing layer where the reaction zone represented by the high value of reaction rate ˙ ω C exists, thereby reaction does not occur in the center of the jet where many fuel droplets are present, because the high evaporation rate ˙ m d in this region leads to fuel-rich conditions that are unconducive for combustion to occur.On the contrary, combustion occurs at the outer edge of the free shear layer where the most favourable mixing of fuel and oxidizer can be achieved, such that the mass fractions of fuel and oxidizer are closer to their stoichiometric proportions.Furthermore, from Fig. 12 (b) it can be confirmed that the mass fraction of OH (an intermediate radical species of combustion reaction) attains maxima very close to the reaction zone, and so does the gas-phase temperature T , due to the heat released by combustion.
Next, to investigate the combustion characteristics and the chemical flame structure in the near wall region of Case 3, the streamwise distributions of spatially averaged values of m d and T , at the radial location of r = 5 mm from the spray axis and the time instance corresponding to m in j = 0 .64 mg [see Fig. 10 (d)], are shown in Fig. 13 .In other words, these streamwise distributions represent the variations of the aforementioned quantities with height H above the wall surface (where H = 0 represents the wall surface).In Fig. 13 (a), it is observed that ˙ m d is zero at all heights above the wall surface in the near wall region, because there are no fuel droplets or liquid film particles at this radial location ( r = 5 mm).However, gaseous fuel exists at this location's near wall region which is a result of the convection of vaporized fuel from the center region of the spray jet.Additionally, it can be seen that Y C 12 H 26 has maximum value very close to the wall surface (i.e., H = 0 ), and a reaction zone representative of a diffusion flame similar to that in Fig. 12 (a), is formed between the maximum values of Y C 12 H 26 and Y O 2 in their respective streamwise distributions.Once again, the blue bands in Fig. 13  ω C occurs.Even though the peak reaction rate in the near wall region is an order of magnitude smaller compared to that at the streamwise location of x = 15 mm [ Fig. 12 (a)], its presence alone indicates that chemical reactions are active in the vicinity of the wall surface.Such occurrences of heat release close to the wall surface will influence the velocity and thermal boundary layers, and the flow-field in the near wall region in general.Furthermore, in Fig. 13 (b), a phenomenon is observed in the vicinity of the wall surface (indicated by the black rectangle very close to H = 0 ) in which, the mass fraction of CO 2 ( Y CO 2 ) increases towards the wall, while that of CO ( Y CO ) decreases towards the wall.This is caused by the oxidation of CO into CO 2 according to the following exothermic Occurrence of the above oxidation reaction can be attributed to the fact that, the reduction in gas-phase temperature T due to the wall heat loss, works in favor of this exothermic reaction to proceed (in order to mitigate the effect of wall heat loss based on Le Chatelier's principle).However, the amount of heat loss to the wall is considerably larger than the amount of heat generated by this exothermic reaction, due to which T near the wall surface decreases despite the progress of the above exothermic reaction.
To examine the behaviour of turbulence in the flow field, the instantaneous distributions of the gas-phase velocity magnitude in the central x − y plane, and the iso-surfaces of the second invariant of the velocity gradient tensor Q are shown in Fig. 14 for Case 3. Here, Q is defined as: where, − is the 1 st invariant, S i j is the strain rate tensor and ω i j is the rotation rate tensor.Q represents the turbulent flow structure with Q > 0 indicating the vorticity-dominated regions of the flow field, and Q < 0 indicating the strain rate-dominated regions.In Fig. 14 , the iso-surfaces of Q are shown for a positive value (i.e., Q > 0 ), which implies that the turbulent structures in the vorticity-dominated regions are depicted.Initially, when the fuel spray is injected from the nozzle into the quiescent gas-phase, the momentum exchange occurring between the dispersed-phase fuel droplets and the gas-phase is responsible for imparting velocity to the gas-phase.From the gas-phase velocity magnitude distribution in Fig. 14 , formation of the mixing layer between the spray jet axis and the ambient gas, which is caused by the velocity gradients between the spray jet and the quiescent ambient gas, is evident.After the spray flame collides with the wall surface, its mostly streamwise momentum (before collision) gets converted into momentum in the wall-parallel direction in the near-wall region.
The distribution of Q iso-surfaces shows that turbulence is generated in the free shear layer (or mixing layer) of the spray jet flame, and develops in the downstream direction.However, turbulence starts decaying in the near wall region (with increasing radial distance from the spray jet axis in the wall-parallel direction) after the spray jet flame impinges on the wall surface, and spreads over the wall.The streamwise distributions of time-averaged gas-phase temperature and viscosity in the near wall region are plotted in Fig. 15 , at various radial locations for Case 3. It is evident that with increasing radial distance from the spray jet axis ( r = 0 mm), the gas-phase temperature increases and so does the viscosity, because viscosity is a temperature-dependent property.The gas-phase temperature increases with radial distance because chemical reactions are active in the vicinity of the wall surface (as already shown in Fig. 13 ), and progress as the radial distance from the spray jet axis increases.Therefore, as the viscosity becomes larger with increasing radial distance, it will play a major role in decaying the turbulence in the boundary layer formed on the wall surface.Furthermore, there is also evidence of turbulence decay brought about by chemical reactions (heat release) occurring in the vicinity of wall, from a previous DNS study on boundary layer flashback in a turbulent channel flow [101] , by our group.

Grid resolution requirements
As mentioned in Section 2 , neither subgrid-scale turbulence model nor turbulence chemistry interaction model is employed in the present simulations of wall impinging spray flames.Computational grid used for the simulations of all the three cases is representative of that used for coarse DNS.At this stage, the grid resolution requirements are analysed for the simulation of Case 3, since the fuel injection velocity is highest in Case 3 among the three cases investigated in this study.Additionally, results of a grid sensitivity test conducted on Case 3 are also presented.It should be emphasized that the present simulations have been considered in the context of turbulent spray combustion.When performing coarse DNS of turbulent spray combustion, the grid size should be capable of resolving the Kolmogorov length scales.But, at the same time, due to the two-way coupling strategy retained between the gas and dispersed phases using the PSI-Cell approach [72] , the grid size must ideally be 10 times larger than the droplet size in order to accurately capture the droplet evaporation dynamics [73,74] .For the present spray flame simulations to match the experiments, the fuel droplets are injected into the computational domain in a polydisperse manner using droplet size distribution PDFs (see Section 2.6 ), and the SMD of the injected fuel droplets ranges from 18.5 μm in Case 1 to 10.3 μm in Case 3.This puts a constraint on how small of a grid size can be used and hence, we have avoided using smaller grid sizes that are of the same order of the droplet diameters, in the core regions of the spray jet flame (for x < 20 mm) where majority of the fuel droplets are present.Otherwise, the predictions for droplet dynamics will be significantly erroneous.The wall-normal grid size x in the near wall regions becomes very fine, and can be of the order of droplet sizes or even smaller than the droplet sizes.The fuel droplets existing in those locations are treated using the distance function technique explained in Section 2.2 .This technique is computationally expensive and therefore, has been applied only in the near wall regions where x becomes equal to or smaller than twice the SMD of the injected fuel droplets.
Additionally, the radial distributions of η and t η at various streamwise locations in the near wall region for Case 3, are presented in Fig. 17 .In the central region of the spray jet flame, closer to the spray jet axis ( r = 0 mm) where the temperatures are relatively low, and many fuel droplets exist, the Kolmogorov length scales η are smaller.However, η and t η increase with radial distance from the spray jet axis at all streamwise locations, due to increasing temperature (caused by heat release from combustion).At the radial locations r > 1 mm, the grid sizes ( x , y and z) are either less than or equal to twice the Kolmogorov length scale, which is often used in DNS of turbulent combustion [29,102] , and within the range recommended by Pope [103] .Therefore, resolution of the grid is adequate for performing the simulations of Cases   2 .As evident from Table 2 , the time step value t used in the simulation is sufficiently small to ensure accurate temporal resolution of all the relevant time scales.
Next, to check whether the grid resolution is adequate for resolving the velocity and thermal boundary layers formed on the wall surface, results of a grid sensitivity test conducted on Case 3 are discussed.For this grid sensitivity test, another simulation was performed using the same conditions, computational domain con- figuration and fuel spray injection parameters as those in Case 3. But, the grid size x , in the wall normal x -direction is made much finer in the near wall region, i.e., for x > 20 mm (grid sizes in the y -and z-directions are kept the same).This finer grid simulation is designated as Case 3'.The simulation of Case 3' is performed using a Cartesian grid consisting of 840 × 560 × 560 grid points in the x -, y -and z-directions, respectively, compared to 720 × 560 × 560 grid points in the original Case 3. The minimum grid size in the x -direction x min = 0 .93 μm in Case 3' (compared to x min = 2 .9 μm in Case 3).
Fig. 18 shows the comparisons between results computed in Case 3 and Case 3', for the streamwise profiles of mean wallparallel radial velocity component V rad , temperature and mixture fraction Z, close to the wall surface and at different radial locations ( r = 1 mm, 2 mm and 3 mm) from the spray jet axis.The results do not show any drastic differences and are mostly similar between the original Case 3 and the finer grid Case 3'.There are some minor discrepancies though (especially in the temperature profiles), which are mainly due to the fact that the number of snapshots in time available for averaging in Case 3' were fewer than those in Case 3. The computational cost of simulating Case  3' was exorbitant, and we were able to perform the simulation of Case 3' for a physical time of 0.4 ms only, compared to 0.75 ms for Case 3.However, the results of Case 3 are nearly the same as those of Case 3', and hence, the resolution of the original grid of Case 3 is sufficient for the analysis presented in this study.

Comparisons with experiment
Results obtained from the simulation of Case 3 ( U in j = 216 m/s) are now compared with measurements obtained from an experiment conducted by Toyota Central R&D Labs., Inc., for the same conditions as that in Case 3's simulation, except for the fuel.Experimental data for the other two cases were unavailable due to proprietary reasons, so we are only able to show the comparisons for Case 3. Additionally, we would like to emphasize that the fuel used in the experiment was Diesel, while that used in the simulations is n -dodecane.Since data of the different hydrocarbon constituents of Diesel used in the experiment, and their exact proportions were unknown, n -dodecane has been chosen as the fuel for the simulations.Fig. 19 illustrates the instantaneous distributions of computed soot radiation emission E soot (line of sight integrated emissive power) for Case 3 at progressing time instances, along with the luminous flame photographs at the same time instances obtained from the experiment conducted by Toyota Central R&D Labs., Inc.Here E soot is calculated as: where κ soot is the Planck mean absorption coefficient of soot and σ = 5 .669 × 10 −8 W/(m 2 K 4 ) is the Stefan-Boltzmann constant.In the above Eq.( 46) , the term κ soot σ T 4 represents the radiant power emitted by soot per unit volume, and E soot is obtained by integrating this term in the z-direction, which implies that E soot is the line of sight integrated soot emissive power.The brightness of the luminous flame in the experimental flame photographs, can be considered to be proportional to the soot emissive power E soot , however, it is not an accurate representation of the flame's radiant power.Therefore, the comparison shown in Fig. 19 is for qualitative purpose only.From the figure, it can be confirmed that after the evaporated fuel auto-ignites, combustion occurs and then the spray flame impinges on the wall.Additionally, the distributions of computed soot emissive power E soot show reasonable qualitative agreement with the experimental luminous flame photographs, in  terms of the temporal evolution of the flame and its brightness caused by the soot radiation emission.
Quantitative comparisons are also performed between the simulation and experiment as shown in Fig. 20 .In Fig. 20 (a), the time variation of temperature at the stagnation point O on the wall surface (point on the wall surface where the spray jet flame impinges directly, also refer Fig. 4 ), obtained from the simulation of Case 3 is compared with experimental data.The simulation result shows an overall good agreement with measurements, however, the stagnation point temperature is slightly under-predicted by the simulation for time after start of injection (ASOI) greater than 0.6 ms.Similar comparison is performed between the computed and measured time variations of wall heat flux at the stagnation point as shown in Fig. 20 (b).The simulation result shows features similar to that of the experiment, like the initial steep rise of the wall heat flux when the spray flame comes in contact with the stagnation point, followed by the wall heat flux getting stabilized and becoming nearly constant in time.In this sense, the simulation is able to capture the time variation trends of the stagnation point's wall heat flux reasonably well.However, the simulation under-predicts the wall heat flux for time ASOI greater than 0.5 ms.This is a direct consequence of the under-predicted temperature in Fig. 20 (a), because the wall heat flux result of the simulation in Fig. 20 (b), has been computed using the wall-side temperature gradient at the wall surface (i.e., by applying Fourier's law of heat conduction).The most likely reasons for the discrepancies in the results of the stagnation point's temperature and wall heat flux in Fig. 20 are as follows: 1.As mentioned earlier, the fuel used in the simulations is ndodecane while that used in the experiment is Diesel (a blend of various hydrocarbons).Hence, relevant fuel vapour properties such as density, thermal conductivity, specific heat capacities, viscosity, calorific value, etc. are different, and are bound to influence the results.Thermophysical properties and transport coefficients of the fuel droplets and the liquid fuel film formed on the wall surface, will also be different from those of Diesel.2. In Fig. 20 (b), it can be seen that the stagnation point's wall heat flux in the experiment is zero at time ASOI = 0, but that in the simulation is non-zero at time ASOI = 0.This is because in the simulation, the initial gas-phase temperature is set to be 1600 K uniformly throughout the computational domain, and the initial temperature distribution inside the wall is assumed to vary linearly from 481 K at the top surface to 463 K at the bottom surface (refer Section 2.6 ).However, at time ASOI = 0 in the experiment, the gas-phase temperature in the near wall regions might not have been 1600 K, and the temperature profile inside the wall might not necessarily be linear.Since detailed descriptions of the initial spatial distribution of gas-phase temperature in the near wall region, and the initial temperature distribution in the wall are unknown from the experiment, the aforementioned assumptions had to be made for the initial conditions in the simulations.As these initial conditions will also dictate how the solution evolves over time, they constitute another source of the discrepancies in Fig. 20 .
Next, the ignition delay time obtained from the simulation is compared with the measured value.In the simulation, ignition delay time τ id is defined as the time after start of fuel injection to the time when the maximum rate of maximum temperature rise occurs in the domain, i.e., the following mathematical expression: where τ SOI is the time at start of injection (SOI) and for the present simulation τ SOI = 0 , and T max is the maximum gas-phase temperature in the combustion field, i.e.T max (t) = max (T (x, y, z, t )) .The ignition delay time τ id computed for Case 3 is 22.14 μs τ id = 34 .9 μs for Case 1, and τ = 26 .43 μs for Case 2), while τ id measured in the experiment is 60 μs (data provided by Toyota Central R&D Labs., Inc.).The simulation of Case 3 yields faster ignition of the evaporated fuel compared to the experiment, since the NA-FPV approach is used as the combustion model.Although the NA-FPV approach is unable to accurately reproduce τ id of spray combustion, the predicted ignition delay time (22.14 μs) is of the same order of magnitude as that of the experiment (60 μs).Moreover, the main objective of this study is to analyze the wall heat transfer characteristics when spray flames impinge on a wall, so in this context the NA-FPV approach is adequate.Another plausible reason as to why the simulation yields shorter ignition delay time, is the assumption that a spray is already established at the nozzle exit (whose size distribution is given by fitting the measured data with a Rosin-Rammler PDF, see Section 2.6 ).In other words, the atomization process at the nozzle exit is assumed to be instantaneous and its time scale is not considered in the simulation.
Fig. 21 depicts the computed time variation of the flame lift-off length in Case 3. In the simulation, flame lift-off length (LOL) is defined as the streamwise distance from the spray nozzle exit to the location of 14% of the maximum OH radical mass fraction corresponding to the quasi-steady state [44] .From the result in Fig. 21 , the quasi-steady LOL is calculated, which comes out to be approximately 3.441 mm.The predicted flame lift-off length is in close agreement with the experimental value of 4.4 mm (data provided by Toyota Central R&D Labs., Inc.), although it's slightly underpredicted.Such an under-prediction is attributed to the limitations of the applied combustion model, i.e., the NA-FPV approach.This flame LOL analysis is presented nonetheless, to check how close the prediction is to the experimental value.

Wall heat flux 3.2.1. Dependence of total wall heat flux on fuel spray injection velocity
Temporal variations of the total wall heat flux at different radial locations on the wall surface from the stagnation point O (i.e., r = 0 , 2, 4, 6, 8 and 10 mm), are presented in Fig. 22 for all the three cases.Here, the radial location of r = 0 mm on the wall surface corresponds to the stagnation point O (refer Fig. 4 ).Total wall heat flux is the sum of the convective and radiative heat fluxes at the wall surface, and the injected fuel mass is representative of time.From Fig. 22 (a) to (d), it can be seen that at the radial locations of r = 0 , 2, 4 and 6 mm on the wall surface, the total wall heat flux increases rapidly after the flame reaches the respective radial locations in each case.Thereafter, the total wall heat flux gradually approaches a quasi-steady state and becomes roughly constant in each case.The general trend observed from this figure is that, at any given radial location, the total heat flux is larger for the case with the higher fuel injection velocity.At the radial location of r = 8 mm in Fig. 22 (e), the total wall heat flux has not yet reached a steady value in any of the cases.On the other hand, as evident from Fig. 22 (f), the flames are yet to arrive at the radial location of r = 10 mm in all the three cases, which is why the total heat flux has not increased in any of them.Therefore, in the following, wall heat fluxes (total, convective and radiative heat fluxes) in the region corresponding to 0 mm ≤ r ≤ 6 mm are discussed.
Radial distributions of time-averaged total heat flux (quasisteady condition) at the wall surface are presented in Fig. 23 for all the three cases.Here, time-averaging of the total heat flux at a given radial location, is performed after the total heat flux reaches a roughly steady value at that particular radial location, as mentioned previously for Fig. 22 .Additionally, to analyze the individual contributions of the convective q w,con v and radiative q w,rad heat fluxes to the total wall heat flux, breakdown of the total heat flux at various radial locations are shown in Fig. 23 for Case 3, in the form of a stacked column bar graph.The percentage above each stacked column represents the ratio of the radiative heat flux to the total heat flux, at the respective radial location on the wall surface in Case 3. The other two cases exhibited similar qualitative trends, in terms of the radial distributions of the contributions from q w,con v and q w,rad to the total wall heat flux, and hence breakdowns of the total heat flux in Cases 1 and 2 are not shown here.From Fig. 23 , it is found that the total wall heat flux is maximum at the stagnation point O ( r = 0 mm), and decreases with increasing radial distance from the stagnation point in all the cases.Furthermore, with increasing fuel injection velocity U in j , a marked increase in the total wall heat flux is observed.This is why at any given radial location on the wall surface, the total heat flux in Case 3 ( U in j = 216 m/s) is the largest, followed by the total heat flux in Case 2 ( U in j = 150 m/s), and then the total heat flux in Case 1 ( U in j = 100 m/s) which has the smallest value.The instantaneous distributions of wall surface temperature at the time instance cor- responding to injected fuel mass m in j = 0 .64 mg in each case, are depicted in Fig. 24 .In each case, it can be seen that the temperature at and around the vicinity of the stagnation point O is the highest, and it decreases with increasing radial distance.Moreover, the higher the fuel injection velocity, the higher the wall surface temperature.Thus, the wall surface temperature trends are consistent with those of the total wall heat flux in Fig. 23 .
Since the total heat flux is expressed as the sum of the convective and radiative heat fluxes, the features of radiative heat flux in Case 3 shown in Fig. 23 are analyzed next.It is observed that the radiative heat flux accounts for a small fraction (2.9%) of the total wall heat flux at the stagnation point, but the magnitude of the radiative heat flux as well as its contribution to the total wall heat flux, increase with radial distance from the stagnation point ( r = 0 mm).The same features were observed for the radiative heat flux in Cases 1 and 2, and the differences among the three cases are only quantitative.To understand why the radial distribution of radiative heat flux exhibits a trend opposite to that of the total heat flux, the streamwise distributions (i.e., variations with height above the wall surface) of spatially averaged gas-phase temperature at different radial locations of r = 0 mm and r = 6 mm, are analyzed in Fig. 25 .Spatial averaging of temperature in this figure bears the same meaning as the spatial averaging of the quantities presented in Figs. 12 and 13 .From Fig. 25 , it is evident that in each case, the gas-phase temperature in the near wall region is higher at the radial location of r = 6 mm, compared to the radial location of r = 0 mm (i.e., above the stagnation point).This happens because, combustion reaction progresses as the radial distance from the spray jet axis increases, due to the availability of vaporized fuel that gets convected from the central region of the spray jet flame, to the radial locations farther away from the spray axis ( r = 0 mm) in the wall-parallel direction, as seen in Fig. 13 .Consequently, the radiation intensity and hence, the radiative heat flux are smaller at the stagnation point, and the radiative heat flux at the wall surface becomes larger with increasing radial distance from the stagnation point.However, the major contribution to the total wall heat flux comes from the convective heat flux rather than the radiative heat flux, so the features of convective heat flux are examined next.

Dependence of convective heat flux on fuel spray injection velocity
Fig. 26 depicts the time variations of convective heat flux q w,con v at the radial locations of r = 0 , 2, 4, 6, 8 and 10 mm on the wall surface, for all the three cases.The temporal variation trends of q w,con v are the same as those of the total wall heat flux observed in Fig. 22 , such at the radial locations of r = 0 , 2, 4 and 6 mm, the convective heat flux rises rapidly after the flame arrives at the respective radial locations, and subsequently attains a roughly steady value in each case.Furthermore, at any given radial location, the convective heat flux is largest in Case 3 ( U in j = 216 m/s) followed by Case 2 ( U in j = 150 m/s) and then Case ( U in j = 100 m/s), which is similar to observation made in 22 for the total wall flux.As with the total flux, the convective heat flux in all the cases is reach steady at the radial location of r = 8 mm, as depicted in Fig. 26  Analyzing the radial distributions of time-averaged convective heat flux (quasi-steady condition) at the wall surface for all the three cases in Fig. 27 , in the same manner as those of the timeaveraged total heat flux in Fig. 23 , the exact same trends can be seen to emerge.The convective heat flux is largest at the stagnation point O and reduces with increasing radial distance from the stagnation point, in each case.Comparing the convective heat flux at any given radial location on the wall surface among the three cases, it is clear that the higher the fuel injection velocity U in j , the larger the convective heat flux.Therefore, the reason as to why the total wall heat flux exhibits characteristics similar to those of the convective heat flux is that, the contribution of convective heat flux to the total wall heat flux is dominant compared to that of the radiative heat flux.
Next, the reason for the radial distribution trend of convective heat flux at the wall surface, and the reason for the increase in convective heat flux with increasing fuel spray injection velocity are investigated.The wall-normal turbulent heat flux arises from the activation of turbulent heat transport in the wall-normal direction by coherent vortical structures near the wall surface [104] .Therefore, an intuitive way to understand the role of turbulent heat transport in promoting the convective heat flux, is to analyze the sweep and ejection events close to the wall surface.Firstly, to understand the dependence of convective heat flux on the fuel injection velocity, the joint probability density function (PDF) distributions of temperature fluctuation T and the wall-normal ( xdirection) component of velocity fluctuation u , close to the wall surface above the stagnation point (i.e., at r = 0 mm) are illustrated in Fig. 28 for all the three cases.The first quadrant of this joint PDF distribution marked by the red square in Fig. 28 (c) corresponds to sweep, and the third quadrant marked by the green square corresponds to ejection.Sweep is the wallward interaction motion that supplies hot gas to the vicinity of the wall surface (i.e., when u > 0 and T > 0 ), and ejection represents the roll-up motion that causes the relatively cold gas in the vicinity of the wall surface to get entrained away from the wall surface (i.e., when u < 0 and T < 0 ).From Fig. 28 , it is clear that the higher the fuel injection velocity U in j , the stronger the correlation between u and T .A stronger correlation between u and T indicates larger contributions of the sweeps and ejections to the turbulent heat transport in the wallnormal direction, which results in the enhancement of convective heat flux.
Next, to understand why the convective heat flux is maximum at the stagnation point and decreases with increasing radial distance from the stagnation point, Fig. 29 shows the joint PDF distributions of u and T close to the wall surface at different ra-  dial locations of r = 0 mm, 3 mm and 6 mm, for Case 3. From the figure, it is evident that the correlation between u and T is the strongest at r = 0 mm, i.e., above the stagnation point O, and this correlation becomes weaker with increasing radial distance from the nozzle axis ( r = 0 mm).Therefore, wall-normal turbulent heat flux is greatly enhanced in the vicinity of the stagnation point, due to the largest contributions of sweeps and ejections to the wallnormal turbulent heat transport.This is the reason as to why the convective heat flux is maximum at the stagnation point, and decreases with increasing radial distance from the stagnation point.

Wall heat flow rate
Description of the wall heat flow rate Q w for all the three case is now presented.Wall heat flow rate Q w is defined as the amount of thermal energy or heat being transferred to the wall per unit time, and is calculated by integrating the wall heat flux over the wall surface area A through which heat transfer occurs, as follows: (48) 3.3.1.Dependence of wall heat flow rate on fuel injection velocity Fig. 30 depicts the temporal variations of total heat flow rate, convective heat flow rate, radiative heat flow rate, and the ratio of radiative heat flow rate to total heat flow rate at the wall surface.Here, total heat flow rate is simply the sum of convective and radiative heat flow rates.In Fig. 30 (d), the ratio of radiative heat flow rate to total heat flow rate is shown only for the time after the spray flame impinges on the wall surface, and the convective and total heat flow rates begin to increase in each case.In terms of the injected fuel mass m in j , the time around which these heat flow rates begin to increase roughly corresponds to m in j = 0 . 2 mg.From Fig. 30 , it can be seen that the total, convective and radiative heat flow rates increase steadily with time in all the three cases, because once the flame has impinged on the wall surface, the flame and the hot gas mixture begin to propagate in the wall-parallel direction with time, thereby increasing the wall surface area through which heat loss occurs in each case.Comparing the total and convective heat flow rates in Fig. 30 (a) and (b), respectively, among the three cases at any given time instance, Case 3 ( U in j = 216 m/s) exhibits the highest values of these heat flow rates followed by Case 2 ( U in j = 150 m/s), and then Case 1 ( U in j = 100 m/s) which has the lowest values of these heat flow rates.On the other hand, the radiative heat flow rates of all the three cases shown in Fig. 30 (c) are not significantly different from one another, and approach nearly the same value with time.Consequently, the ratio of radiative heat flow rate to total heat flow rate in each case, shown in Fig. 30 (d) becomes roughly constant with time, with Case 1 having the highest ratio followed by Case 2 and then Case 3. Fig. 31 shows the breakdown of the total heat flow rate Q w,total in terms of its contributions from the convective Q w,con v and radiative Q w,rad heat flow rates for all the three cases, at the time instance when the injected fuel mass m in j = 0 .64 mg in each case.It is evident that, the higher the fuel injection velocity, the larger the total and convective heat flow rates through the wall surface.In Section 3.2.2, the influence of fuel injection velocity on the convective heat flux at the wall surface was discussed, and it was found that the convective heat flux increases with increasing fuel injection velocity.As shown in Fig. 31 , the convective heat flow rate accounts for about 70% ∼ 82% of the total heat flow rate (depending on the case), while the contribution of radiative heat flow rate to the total heat flow rate remains nearly the same in all the three cases (magnitude-wise).Therefore, the total heat flow rate also increases with increasing fuel injection velocity.
On the contrary, the influence of fuel injection velocity on the radiative heat flow rate seems rather insignificant from Fig. 31 .This is due to the fact that, the fuel injection velocity does not have an appreciable influence on the flame/gas-phase temperature, as shown previously in Fig. 25 , and thus, the intensity of radiation emitted by soot and the gas-phase species does not change significantly among the three cases.Consequently, the ratio of radiative heat flow rate to total heat flow rate is larger when the fuel injection velocity is smaller.In the present simulations, this ratio lies in the range of 18.1% (for Case 3) to 29.7% (for Case 1), which is in good agreement with the values measured in previous experiments [105,106] .

Correlation between Nu and Re
The relation between Nusselt number Nu and Reynolds number Re corresponding to heat loss at the wall surface, during spray flame-wall interaction is now examined.In an attempt to estimate the heat loss through the walls of IC engines, previous experimental studies [7][8][9][10][11][12][15][16][17] have investigated the relation between Nu and Re based on the principle of convective heat transfer, under the assumptions of quasi-steady conditions and a constant Prandtl number.These investigations led to the development of empirical correlations for calculating the heat transfer coefficient, and showed that Nu has a power law dependence on Re as expressed in Eq. ( 1) .The present simulations mimic the phenomenon of spray flames impinging on a wall under CI engine-like conditions, occurring during the combustion period of a CI engine.It should be noted that the correlation shown in Eq. ( 1) , is applicable for pure convection heat transfer problems and not for situations in which radiative heat transfer is also involved, because strictly speaking, the Nusselt number Nu is a parameter that provides a measure of the convection heat transfer occurring at a surface, and radiation heat transfer is not a function of Re .However, in the experiments, the measured wall heat flux might implicitly include the contribution of radiation heat transfer.Furthermore, as described in Section 3.3.1 , the contribution of radiative heat flow rate to the total heat flow rate at the wall surface varies from 18.1% to 29.7% (depending on the case), and as also pointed out by some previous experimental studies [18,105,106] , such contributions from the radiation heat transfer certainly cannot be neglected.
Therefore, in the following, two types of Nusselt number are defined and their dependence on Re is analyzed.The first Nusselt number Nu total is calculated based on the average total heat flux at the wall surface (i.e, convective heat flux + radiative heat flux) in the quasi-steady state, while the second Nusselt number Nu con v is calculated based on the average convective heat flux in the quasisteady state (i.e., the conventional definition).These two Nusselt numbers and the Reynolds number Re are expressed as: where h heat is the average heat transfer coefficient, H is the characteristic length, T is the average gas-phase temperature (in quasisteady state) in the near wall region and is used as the characteristic temperature, T w is the wall temperature, q w,total and q w,con v are the average total wall heat flux and convective heat flux (in quasisteady state), respectively, U is the characteristic velocity, and ρ, μ and λ are the gas-phase density, dynamic viscosity and thermal conductivity, respectively.Following previous experiments [16,17] , the characteristic length H used in Eqs. ( 49) -( 51) is set equal to the spray impingement distance value of 0.025 m (i.e., the distance from the spray nozzle exit to the wall surface) in this study.In the previous experiments [16,17] , the local gas-phase velocity and thermophysical properties such as ρ, μ and λ could not be measured directly.Hence, values of the thermophysical properties corresponding to either air [17] or the ambient conditions [16] were used to calculate Nu and Re .To get an estimate for the characteristic flow velocity U, Mahmud et al. [17] used the waveform signals of wall heat flux measured by the thermocouple heat flux sensors embedded in the flat wall.On the other hand, Inagaki et al.
[16] used the flame snapshots in the near wall region at different time instances (obtained from high-speed imaging) to estimate U.However, in this study, the maximum value of average radial velocity (wall-parallel component) in the velocity boundary layer is chosen for the characteristic velocity U Eq. (51) .λ in Eqs. ( 49) and ( 50) is the average gas-phase thermal conductivity at the wall surface.Furthermore, for the present analysis, the values of gas-phase density ρ, dynamic viscosity μ and characteristic temperature T are obtained by temporal and spatial averaging.Here, spatial averaging implies averaging over a cylindrical volume of radius r = 6 mm, with the axis of this imaginary cylinder being coaxial with the spray nozzle axis.Different values for the height of this imaginary cylinder in the wall-normal x -direction (starting from the wall surface) are used for the spatial averaging of these quantities, such as 0.5 mm, 1 mm, 2 mm and 3 mm.This is done in order to check the dependence of Nu -Re correlations, on the height of the imaginary cylinder used for the spatial averaging of ρ, μ and T , since some of these quantities ( ρ, μ and λ) are defined differ- ently in the previous experimental studies [16,17] , as mentioned above.However, in this study, complete descriptions of the spatial and temporal variations of all the quantities are readily available from the simulation results.Hence, temporal and spatial averaging of the relevant quantities can be performed in the near wall regions of interest, and used for calculating N u total , N u con v and Re in Eqs. ( 49) - (51) .It should be noted that in the present simulations, the reaction zones in the near wall region exist up to a height of 1 mm ∼ 2 mm above the wall surface.Fig. 32 illustrates the dependences of Nu total and Nu con v on Re by plotting the pairs of ( Nu total , Re ) and ( Nu con v , Re ) values calculated for each simulation case.The sub-figures (a) to (d) in Fig. 32 show the relation between Nu total and Re , as well as the relation between Nu con v and Re , when spatial averaging of the quantities ρ, μ and T are performed using different heights of the imaginary cylindrical volume explained above.The solid lines represent the best fit lines for the correlation: Nu total = C total Re n total , while the dashed lines represent the best fit lines for the correlation: Nu con v = C con v Re n con v .Here, C total and C con v are the regression coefficients, and n total and n con v are the correlation indices of the respective correlations.It can be seen that the value of the correlation index n total ranges from 0.48 to 0.50, indicating that n total does not vary drastically depending on the height of the aforementioned cylindrical volume used for the spatial averaging of density ρ, dynamic viscosity μ and characteristic temperature T , in the near wall region.Considering the dependence of Nu con v on Re in Fig. 32 , as in the case of correlation index n total (corresponding to the relation between Nu total and Re ), the correlation index n con v corresponding to the relation between Nu con v and Re , also does not show any significant variation depending on how the spatial averaging is performed in the near wall region, i.e., sub-figures (a) to (d) in Fig. 32 .The value of n con v ranges from 0.55 to 0.57.The regression coefficients ( C total and C con v ) and correlation indices ( n total and n con v ) of the best fit lines for the respective Nu -Re correlations, are summarised in Table 3 for the corresponding sub-figures in Fig. 32 .From Fig. 32 , the mean value of n total is 0.49 and the mean value of n con v is 0.56, indicating that the value of correlation index n total is smaller than that of n con v .The reason for this is that, Nu total is calculated based on the total heat flux at the wall surface, which includes the contribution from the radiative heat flux as well.As described previously in Section 3.3.1 , the radiative heat flow rates were nearly the same among all the three cases investigated in this study, which means that the average radiative heat flux at the wall surface is approximately constant regardless of Re .So, when a correlation of the form shown in Eq. ( 1) is applied to Nu total and Re as in the previous experimental studies, the value of its correlation index n total will deviate from the correlation index value defined in the conventional sense (i.e., for pure convection heat transfer).The extent to which the value of n total deviates, will depend on the ratio of radiative heat flux to the total heat flux, which in turn depends on the combustion conditions such as flame temperature, ambient pressure, fuel type, etc.Furthermore, the present value of correlation index n total = 0 .49 is not that far off from the value of n = 0 .4 reported in the experiments by In-

Table 3
Correlation indices and regression coefficients of the best fit lines for correlations: Nu total = C total Re n total and Nu con v = C con v Re n con v in Fig. 32  ) and its contributions from radiative ( Q w,rad ) and convective ( Q w,con v ) heat flow rates, at the time instance when injected fuel mass m in j = 0 .64 mg in each case.The percentages in each stacked column indicate the ratios of radiative and convective heat flow rates to total heat flow rate.
agaki et al. [15,16] .Therefore, it becomes necessary to distinguish between the contributions of convective and radiative heat fluxes, if reliable correlations between Nu and Re are to be formulated.

Conclusions
In this study, a 3-D numerical analysis of n -dodecane wall impinging spray flames under CI engine-like conditions, was carried out to investigate the combustion characteristics, and the effect of fuel spray injection velocity on the wall heat flux (convective, radiative, and total heat fluxes).The dependence of wall heat flow rate (convective, radiative and total heat flow rates) on fuel injection velocity was also examined.The present simulations incorporated models for considering radiative heat transfer from soot and gas-phase species, and liquid fuel film formation on the wall surface.Additionally, Conjugate Heat Transfer (CHT) was also incorporated in the simulations for coupling the heat transfer occurring at the wall surface from the gas-phase (convective and radiative) and the liquid fuel film, with the heat conduction occurring inside the solid wall.Simulations were performed for three cases with fuel injection velocities of 100 m/s (Case 1), 150 m/s (Case 2) and 216 m/s (Case 3).The main findings of this research are as follows: 1. Total wall heat flux is largest in the vicinity of the stagnation point and reduces with increasing radial distance from the stagnation point in each case.This feature is consistent with that of the convective heat flux, which has the major contribution to the total wall heat flux, and is caused by the fact that the sweep and ejection events which characterize the turbulent heat transport near the wall surface, are largest near the stagnation point and decrease with increasing radial distance from the stagnation point in each case.Furthermore, the convective heat flux (and hence the total wall heat flux) increases with the fuel injection velocity, because with increasing fuel injection velocity, these sweep and ejection events become more frequent and enhance the convective heat flux.2. Radiative heat flux at the wall surface becomes larger as the radial distance from the stagnation point increases in all the cases, because the flame temperature increases with increasing radial distance from the stagnation point in each case.The influence of fuel injection velocity on the radiative heat flow rate is found to be rather insignificant.This is due to the fact that the radiative heat flux is determined by the flame temperature, and the flame temperature is not significantly affected by the fuel injection velocity among the three cases.However, because the total heat flow rate increases with the fuel injection velocity, the ratio of radiative heat flow rate to the total heat flow rate decreases from 29.7% in Case 1 to 18.1% in Case 3, meaning this ratio decreases with increasing fuel injection velocity.3. The contribution of radiative heat flow rate to the total heat flow rate is non-negligible for the conditions under which each case is simulated.So, when the average total wall heat flux is used for constructing a correlation between Nu and Re of the form Nu ∝ Re n (based on the principle of convective heat transfer), corresponding to heat loss through the wall, the value of its correlation index n will be different from that of the con-ventional Nu -Re correlation applicable for pure convection heat transfer.Therefore, it is necessary to clearly distinguish between the contributions of convective and radiative heat flow rates, in order to construct a robust correlation between Nu and Re that can be used for calculating the heat transfer coefficient accurately.
This work has investigated the dependence of heat loss through the wall on fuel injection velocity (and hence the fuel injection pressure), and attempted to elucidate the relation between Nusselt number Nu and Reynolds number Re relating to wall heat transfer during spray flame-wall interaction process.However, there are several other parameters which influence the wall heat loss in CI engines, such as fuel temperature at injection, spray nozzle diameter, impingement distance and fuel spray impingement angle, radiative properties of the wall surface, bore size, operating pressure and temperature, etc.In order to develop robust and improved correlations between Nu and Re that take into account the influences of the aforementioned parameters on wall heat transfer, much more work is needed.Such correlations can be used for accurate predictions of the heat transfer coefficient in CI engines, which will have strong implications on designing lowemission CI engines with better thermal efficiencies.The above issues are open questions in the fields of thermal, combustion, fluid and automotive engineering, and researches pertaining to these issues will form part of our future work.With the ever-increasing capabilities of computing resources around the world, the application of numerical simulation tools such as Large-Eddy Simulation (LES) and Direct Numerical Simulation (DNS), open up new possibilities for investigating various mechanisms associated with complex spray combustion problems involving flame-wall interactions in detail.One aspect that is integral to spray combustion problems and of particular interest to the combustion engineering community, is the liquid atomization process.The computational cost of accurately reproducing fuel atomization process (both primary and secondary atomization) using fully-resolved simulations has been quite prohibitive until now.Moreover, for precise evaluation of the effects of spray nozzle parameters, fuel injection pressure, etc. on liquid fuel atomization and spray combustion itself, the use of presumed droplet size distribution PDFs or even existing liquid atomization/breakup models may not be adequate for certain applications.Therefore, numerical simulations of spray combustion incorporating well-resolved methods that can accurately capture the physics of liquid atomization process are strongly desired.Such simulations will soon be possible by using recently developed state-of-the-art supercomputers.Findings from these numerical investigations will no doubt complement those of experimental studies which are ever so indispensable, and will be essential for addressing challenges related to energy-efficiency, pollution and CO 2 emissions control.
and ˙ ω C is the reaction rate of C (explained later).D Z and D C are the diffusion coefficients of Z and C, respectively, and are assumed to be equal to the gaseous thermal diffusivity D h under the unity Lewis number assumption ( Le = 1 .0 ) as:

Fig. 2 .
Fig. 2. Schematic of liquid fuel film formed on the wall.A particle-based approach is used for calculating the wall film dynamics.

Fig. 3 .
Fig. 3. Schematic of conjugate heat transfer among gas, film and wall.

Fig. 4 .
Fig. 4. Schematic of the computational domain used in the simulations.

Fig. 7 .
Fig. 7. Droplet size distributions of injected fuel droplets in all three cases.

Fig. 8 .
Fig. 8. Instantaneous distributions of gas-phase temperature along with spray and film particles (particle sizes are exaggerated) in the central x − y plane for Case 1, at progressing time instances corresponding to increasing values of injected fuel mass m in j .Only some fuel droplets are shown in these temperature distributions for the sake of clarity.

Fig. 9 .
Fig. 9. Instantaneous distributions of gas-phase temperature along with spray and film particles (particle sizes are exaggerated) in the central x − y plane for Case 2, at progressing time instances corresponding to increasing values of injected fuel mass m in j .Only some fuel droplets are shown in these temperature distributions for the sake of clarity.

Fig. 10 .
Fig. 10.Instantaneous distributions of gas-phase temperature along with spray and film particles (particle sizes are exaggerated) in the central x − y plane for Case 3, at progressing time instances corresponding to increasing values of injected fuel mass m in j .Only some fuel droplets are shown in these temperature distributions for the sake of clarity.
(a) and (b) indicate the height above the wall surface where the peak reaction rate ˙

Fig. 11 .Fig. 12 .Fig. 13 . 2 →
Fig. 11.Instantaneous distributions of the iso-surfaces of gas-phase temperature at 1700 K (coloured in yellow) and 2200 K (coloured in red), along with the dispersedphase fuel droplets coloured by their respective evaporation rate ˙ m d , at the time instance corresponding to injected fuel mass m in j = 0 .64 mg for Case 3. (Dimensions shown in all Cartesian directions are in meters).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 14 .
Fig. 14.Instantaneous distributions of gas-phase velocity magnitude in the central x − y plane, and iso-surfaces of the second invariant of the velocity gradient tensor Q for Case 3.
Firstly, the Kolmogorov length scale η and the Kolmogorov time scale t η are estimated as η = (ν 3 / ) 1 / 4 , and t η = (ν/ ) 1 / 2 , where ν and are the mean kinematic viscosity and the mean dissipation rate of turbulent kinetic energy, respectively, at a given location in the flowfield.Evolutions of the computed Kolmogorov length scale η and Kolmogorov time scale t η along the free shear layer of the spray jet flame in Case 3 are depicted in Fig. 16 .It can be seen that, the size of the Kolmogorov length scale η gradually decreases with increasing streamwise distance from the nozzle exit, due to increasing turbulence and lower temperatures in the central regions of the spray jet flame (caused by the latent heat of fuel droplet vaporization).However, η becomes nearly constant for x > 18 mm, and the Kolmogorov time scale t η first decreases with streamwise distance from the nozzle exit and then increases slightly.The Kolmogorov length scales are smaller than the local grid size of x = y = z = 40 μm, but these grid sizes are a necessary trade-off when performing spray combustion simulations.

Fig. 15 .
Fig. 15.Streamwise distributions of time-averaged temperature and dynamic viscosity of gas-phase in the near wall region, at various radial locations for Case 3.

Fig.
Fig. Evolutions of Kolmogorov length η and time t η scales along the free shear layer of the spray jet flame in Case 3.

Fig. 17 .
Fig. 17.Radial distributions of Kolmogorov length η and time t η scales, at various streamwise locations in the near wall region for Case 3.

Fig. 18 .
Fig. 18.Comparisons of the streamwise distributions of mean radial velocity V rad (wall-parallel component), mean temperature and mean mixture fraction close to the wall surface, at different radial locations between Case 3 and Case 3'.

Fig. 19 .
Fig. 19.Comparison between (a) computed soot radiation emission distributions (line of sight integrated) of Case 3 in the x − y plane at different time instances, and (b) Luminous flame photographs of experiment [15] conducted by Toyota Central R&D Labs., Inc. at the same time instances.

Fig. 20 .
Fig. 20.Comparisons of the variations of (a) temperature of the stagnation point on wall and (b) total wall heat flux at the stagnation point in Case 3, with experimental measurements[16] made by Toyota Central R&D Labs., Inc.

Fig. 21 .
Fig. 21.Time variation of the flame lift-off length obtained from the simulation of Case 3.

Fig. 22 .
Fig. 22. variations of total heat flux at the wall surface at different radial locations for each case.Here, injected mass of fuel is representative of time.

Fig. 23 .
Fig. 23.Radial distributions of time-averaged total heat flux at the wall surface for all three cases, with the stacked column bar graph showing the contributions of convective heat flux q w,con v and radiative heat flux q w,rad to the total heat flux in Case 3. The percentages above the stacked columns for the different radial locations on the wall surface, indicate the ratio of radiative heat flux to total heat flux at the respective radial locations in Case 3.
(e).Similarly, as shown in Fig. 26 (f), the convective heat flux has not increased in any of the cases at r = 10 mm, because the flames have not yet reached this radial location.Hence, the following discussion pertains to convec-tive heat flux at the wall surface, in the region corresponding to r = 0 mm ∼ 6 mm.

Fig. 24 .
Fig. 24.Instantaneous distributions of wall surface temperature in the y − z plane for each case, at the time instance when injected fuel mass m in j = 0 .64 mg .

Fig. 25 .
Fig. 25.Streamwise distributions of spatially averaged temperature at different radial locations of (a) r = 0 mm and (b) r = 6 mm , for each case at the time instance when injected fuel mass m in j = 0 .64 mg .

Fig. 26 .Fig. 27 .
Fig. 26.variations of convective heat flux at the wall surface at different radial locations for each case.Here, injected mass of fuel is representative of time.

Fig. 28 .
Fig. 28.Joint probability density function distributions of wall-normal component of velocity fluctuation u and temperature fluctuation T for (a) Case 1 ( U in j = 100 m/s), (b) Case 2 ( U in j = 150 m/s) and (c) Case 3 ( U in j = 216 m/s), at the radial location of r = 0 mm (above the stagnation point O).

Fig. 29 .
Fig. 29.Joint probability density function distributions of wall-normal component of velocity fluctuation u and temperature fluctuation T , at different radial locations of (a) r = 0 mm , (b) r = 3 mm and (c) r = 6 mm for Case 3 ( U in j = 216 m/s).

Fig. 30 .
Fig. 30.Time variations of (a) total heat flow rate, (b) convective heat flow rate, (c) radiative heat flow rate, and (d) ratio of radiative heat flow rate to total heat flow rate.Here, injected mass of fuel is representative of time.

Fig. 32 .
Fig. 32.Dependence of the total Nusselt number Nu total and convective Nusselt number Nu con v on the Reynolds number Re .Spatial averaging of the relevant quantities is performed over a cylindrical volume of radius r = 6 mm, with height of the cylinder above the wall surface equal to (a) 0.5 mm, (b) 1 mm, (c) 2 mm, and (d) 3 mm.Data points for the relation between Nu total and Re are represented by circles, and those for the relation between Nu con v and Re are represented by triangles.Meaning of symbol colours: Black = Case 1, Red = Case 2, Blue = Case 3. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Cases simulated in this study and their respective conditions and fuel injection parameters.

Table 2
Various time scales encountered in the simulation of Case 3.
Injection time scale, t in j = D n /U in j 0.384 μs Minimum evaporation time scale, min( t e v ap = ρ/S ρ ) ≈ 7 μs Ignition delay time, τ id 22.14 μs Minimum Kolmogorov time scale, min( t η ) 4.67 μs Simulation time step value, t ≈ 0 .01 μs . Total heat flow rate ( Q w,total