A Novel One- and Zero-Dimensional Model for Turbulent Jet Ignition

Turbulent jet ignition (TJI) is a promising combustion technology for burning highly diluted air-fuel mixtures. Computationally efficient models to assess the effect of the operating conditions and design parameters on the ignition propensity and timing are of paramount importance for the development of combustion systems employing TJI. To this end, a one-dimensional (1-D) jet model, which is based on the solution of the section integrated mass and momentum conservation equations, is derived in the present study. The model is extended with two additional transport equations for the turbulence intensity and the ignition precursor/tracer, that marks the ignition event. One-dimensional transient flamelet calculations are performed to generate tables for the ignition precursor source term that account for the turbulence and chemistry interaction. Further simplification of the model is carried out to obtain a novel penetration correlation and a computationally inexpensive Lagrangian ignition model. The extended jet model is hierarchically validated using available literature data for non-reactive and reactive jets, as well as experiments conducted in a state-of-the-art optically accessible prechamber. The derived model is able to reproduce both flow-related quantities (velocity and turbulence profiles, jet penetration) and the ignition delay time under different variations. This study also illustrates how numerical simulations in canonical configurations (one-dimensional flamelet) can be used in practical applications of TJI.


Introduction
Combustion initiation of a lean air-fuel mixture via one or multiple turbulent hot jets has been extensively investigated due to its relevance in multiple applications including explosion protection (Sadanandan et al. 2007), pulse detonation engines (Lieberman et al. 2002), wave rotor combustors (Wijeyakulasuriya et al. 2010) and internal combustion engines equipped with a prechamber (Attard et al. 2010). The mechanisms of turbulent jet ignition involve a complex interplay between chemical reactions and turbulent mixing that determines the ignition propensity, timing and location. Gussak (1975) pioneered in investigating the fundamental aspects of TJI and developed the so-called Avalanche Activated Combustion method, otherwise know as 'LAG'. This process refers to the ignition and rapid burning of ultra-lean mixtures through torch jets containing active radicals, which are generated by the incomplete combustion of a rich air-fuel mixture inside a prechamber. Significant reduction of the combustion duration and mitigation of the cyclic combustion variability (CCV) can be achieved with TJI in comparison to a conventional ignition system when applied it is applied in internal combustion engines (Gussak et al. 1979;Pitt et al. 1983). In addition to the aforementioned experimental works, early research on TJI included numerical (Goh 1971) and theoretical studies (Phillips 1963(Phillips , 1972a for the determination of the maximum experimental safe gap (MESG), which is the largest diameter between two vessels that does not allow an explosion to be transmitted from one to the other.
Optical experiments (OH* chemiluminescence and Schlieren) and numerical simulations were combined to understand the effect of temperature and mixing on the ignition location and timing of a hot jet issuing into a combustible hydrogen-air mixture (Sadanandan et al. 2007). Ghorbani et al. (2014Ghorbani et al. ( , 2015 employed a stand alone probability density function (PDF) method to numerically investigate the mechanisms of ignition of a hydrogen-air mixture by a hot turbulent jet for different thermochemical and flow conditions. Ignition was observed near the jet tip, as in the experimental work of Sadanandan et al. (2007), where mixing is sufficiently slow and does not inhibit chemical reactions. Deflagration initiation in a premixed mixture via a starting hot jet of adiabatic combustion products was numerically simulated in Iglesias et al. (2012). The reactant diffusivity is seen to have a dominant effect on the critical nozzle diameter, while the Reynolds number exerts minor influence on it. Carpio et al. (2013) observed numerically that the critical radius increases with the nozzle velocity and that the minimum critical radius corresponds to stoichiometric conditions.
The jet structure and ignition location for variations of the nozzle number and diameter was experimentally investigated by means of OH* and CH* chemiluminescence in Gentz et al. (2015a). The main outcome of this study is that the nozzle diameter which results in the fastest combustion varies depending on the mixture equivalence ratio. Biswas et al. (2016) examined the ignition mechanisms of methane-air and hydrogen-air mixtures by a hot turbulent jet for variations of pressure, temperature, equivalence ratio, orifice diameter and spark plug location. A global Damköhler number Da = u � l f s L l I was introduced, which was found to characterise the transition from extinction to ignition. The combustion mode can either resemble conventional premixed flame propagation or jet ignition, depending on the value of the Damköhler number. Xu et al. (2019) employed optical experiments (OH* chemiluminescence and Schlieren) in conjunction with 3-D Reynolds Averaged Navier Stokes (RANS) simulations to validate a new modelling approach for TJI based on the G-equation and investigate the combustion inside the main chamber. Studies on canonical configurations, such as the ones conducted in Sadanandan et al. (2007), Fischer et al. (2017a), Sidey and Mastorakos (2018), can yield further insights and enable to isolate the effect of mixing and chemical reactions on ignition initiation. A common observation is that ignition is unlikely to occur when mixing rate is very high or the hot jet temperature is low.
Numerical simulations of TJI have been limited in RANS (Xu et al. 2017(Xu et al. , 2018Chinnathambi et al. 2015;Gentz et al. 2015b;Thelen et al. 2015;Thelen and Toulson 2017;Wang et al. 2018), Large Eddy Simulations (LES) (Validi et al. 2017;Feyz et al. 2019a;Malé et al. 2019) and more recently Direct Numerical Simulations (DNS) (Validi and Jaberi 2018;Qin et al. 2018;Benekos et al. 2020). Although some preliminary lower-order modelling efforts were focused on capturing the mass entrainment rate and temporal evolution of the temperature profile in suddenly started jets (Feyz et al. 2019b(Feyz et al. , 2018, with the scope of extending the model to ignition delay prediction, a formally derived and validated low-order ignition model for TJI is missing from the existing literature. In addition, while valuable insights were obtained in canonical configurations, they did not enable quantified predictions in applications of TJI.
The objective of the present study is the development of a reduced order and computationally efficient jet model for the simulation of ignition delay by a transient hot jet. To this end, 1-D jet models from the literature are modified and extended, so that the mixing rate and ignition precursor distribution can be determined. Simulations in a 1-D reactiondiffusion mixing layer enable the computation of the ignition precursor source term, which is stored in tables and accessed during the main-simulation loop. The 1-D model is further simplified such that a penetration correlation and a simplified Lagrangian ignition model is derived.
This paper is organised as follows: The next section briefly presents the experimental set-up and the optically accessible prechamber (OPC). The derivation of the one-and zerodimensional flow and ignition sub-models is given in the subsequent sections. The Results and Discussion section is dedicated to the extensive validation of the proposed model. Concluding remarks and directions for future research are included in the final section.

Experimental Set-Up
A newly developed, constant volume divided combustion chamber, named the optically accessible prechamber (OPC) has been designed and built to allow the study of TJI under different conditions. The testrig contains two volumes, which are connected through a single interchangeable orifice with 4 mm length as shown in Fig. 1a. The smaller of the two volumes, (identified as the prechamber), has approximately orthogonal base 12.15 x 12.15 mm and a height of 25 mm. The prechamber contains an ignition source (spark plug), which allows the ignition of the mixture. The orifice connecting the prechamber to the main chamber can be changed in shape and diameter, and can fit a nozzle with maximum diameter of 4 mm. The shape of the main chamber inner surface can be seen in Fig. 1b with its dimensions given in Fig. 1a. The volume of the main chamber is 186'988 mm3, whereas the prechamber is 3'602 mm3 (1.9% of the total volume).
Both chambers have windows to allow optical access into the respective volumes, which enable optical measurements during combustion. The prechamber is fitted with two opposing windows, which allow optical access through the complete prechamber volume. The main chamber also has two opposing windows, allowing optical access from the nozzle 1 3 exit to the bottom of the chamber. A removable opposing wall is placed 40 mm away from the nozzle exit.
The operation and control sequences of the OPC allow the accurate setting of conditions in both chambers prior to ignition. The charge air pressure and temperature is feedback controlled through an intake and leakage valve system and an electrical heating system, respectively. The fuel in the main chamber is injected using a hollow-cone piezo-injector, and the mixture composition is determined using the pressure signal and a partial pressure calculation. An external mixing tank is used to prepare the mixture for the prechamber, allowing the separate determination of composition in the two chambers. The premixed mixture is injected directly into the prechamber prior to ignition. The fuel used in this study is pure methane (N45, 99,995 mol%). The operating conditions possible for the current set-up are presented in Table 1. Simultaneous high-speed Schlieren and two-dimensional OH* chemiluminescence imaging can be employed to provide the combustion and penetration characteristics of the turbulent jet emanating into the main chamber. For this purpose, two LaVision HSSX high-speed cameras were used. The frame rate was 20 kHz for the reactive experiments but further increased to 200 kHz when the main chamber was filled with air for accurate penetration measurement of the highly transient reactive jet. A Matlab routine to derive the reactive jet penetration from the Schlieren frames was developed.

Model Development
This section is further divided into three different subsections in which the flow, turbulence and ignition sub-models are derived and the fundamental modelling assumptions are presented. Whenever possible, lengthy derivations that do not belong to the core of the model are given in the corresponding Appendices.

Flow
A gaseous jet either of reactive air-fuel mixture or hot combustion products and uniform exit velocity is penetrating into a quiescent mixture (either reactive or non-reactive) as  Fig. 2 1 . The jet gradually entrains the main chamber mixture and spreads radially. The axial velocity of the jet decays from the momentum conservation requirement. A constant spreading angle j is assumed for the jet. In reality, the jet angle varies throughout the jet and in time as a function of the entrainment rate. An exact solution of the Navier-Stokes is required to fully capture this transient behaviour, which is, however, out of the scope of the present study.
In addition to the fixed jet angle the following assumptions and approximations are made: 1. No mean flow convection e.g. swirl. Therefore, the jet is symmetrical around its axis. 2. The flow is incompressible. For typical conditions encountered in prechamber gas engines, the nozzle Mach number is below one. Compressibility effects are more pronounced near the nozzle exit where the velocity is high but quickly decay further downstream. 3. The chamber pressure is assumed uniform. 4. Molecular and turbulent diffusion fluxes on the jet control volumes are neglected. This is justified on the basis of the leading order axial convention terms and the fact that the control volumes of the jet are large enough so that velocity gradients are minimal on the lateral sides. 5. Negligible viscous and turbulent forces acting on the control volumes. Despite the fact that viscous shear forces may be substantial in certain jet regions, they are contained within the control volume and, therefore, their net contribution is zero. The ones outside of the control volumes are small in magnitude, since the flow velocities and the corresponding gradients are small (Musculus and Kattke 2009). 6. The radial distribution of the ensamble averaged axial velocity 2 u follows the profile suggested by Abramovich (1963), that is where = r∕r O , with r O being the outer jet radius at a given downstream location x, and a the shape factor. The shape factor varies in the developing jet region and approaches a = 3∕2 in the far field. The methodology to derive the shape factor is detailed later. The ensamble averaged mixture fraction Z is determined as follows for the general case of non-unity Schmidt number (Sc). 7. The temperature variation across the jet is purely an effect of mixing of the two streams (main chamber and prechamber fluid) and the heat released does not result in a deviation from this temperature profile. Since the objective of the present model is to determine the ignition delay, during which there is radicals build up but no substantial heat is released, it holds sufficiently well. 8. Local density is computed as follows Sc 2 All variables that appear on this study are ensamble averaged unless otherwise noticed.

3
assuming ideal mixing between the prechamber and main chamber mixture. In Eq. (3), pc and mc is the pre and main-chamber density, respectively. The temperature ratio t r = T pc ∕T mc = mc ∕ pc is introduced (the last equality follows from the uniform pressure assumption) and Eq. (3) becomes 9. The thermochemical conditions of the two streams (pre and main-chamber) are constant during turbulent jet ignition. In reciprocating internal combustion engines, the ignition delay period is short due to the high pressure (Kyrtatos et al. 2018), resulting in approximately constant thermochemical conditions. This assumption holds well also for optical experiments in constant volume chambers. While the ignition delay might be long in optical experiments, there is no compression work done by a moving piston and, therefore, negligible pressure variation before combustion initiation.
The jet domain is discretised in cells as shown in Fig. 3. Under assumptions (1) to (5), the mass and momentum conservation equations in each cell simplify to where the term on the left hand side is the rate of change of a given quantity (either the mass of prechamber mixture (m Z ) or momentum (P)) over the whole cell, while the two terms on the right hand side represent the inflow and outflow difference from the cell surfaces.
The mass of prechamber mixture and momentum in each cell can be found from the following volume integrals is solved numerically with the known shape factor a i and the fixed parameters t r and Sc to obtain Z j cl,i . The value of the centerline velocity u j cl,i can then be found from Equations and yield the flow terms ṁ j Z,i and Ṗ j i based on the known centerline values. Equation (14) is then employed to obtain the cell integral m j+1 Z,i in the next time step j + 1.
To fully close the system of equations, the shape factor a needs to be known in each axial location. As in Musculus and Kattke (2009), it is assumed that the momentum is conserved along the jet and, therefore, the momentum flux at any axial position equals the nozzle momentum flux. Moreover, the centerline velocity and the centerline mixture fraction are set to the nozzle values; that is u cl = u j and Z cl = 1 . This results in the following non-linear integral equation through which the value of shape factor a can be derived The complete derivation of Eq. (19) and an approximate polynomial expression based on a simplified density profile are given in B. The polynomial expression derived and Eq. (19) are compared with the one suggested by Musculus and Kattke (2009) in the Results and Discussion section.

Turbulence and Scalar Dissipation Rate
The simple form for the transport equations of mass and momentum is a direct consequence of assumptions (4) and (5) that enabled only the convective terms to be retained in the final conservation equation. Nonetheless, the same assumptions do not apply for the turbulence intensity (k) transport equation. It can be easily seen from Fig. 4a that the various terms of the turbulence intensity equation, as obtained from experimental measurements of Panchapakesan and Lumley (1993) in the far field, are of the same order of magnitude. Additionally, the profile of the turbulent velocity ( u ′ ) in an axial cross-section cannot be easily parametrised. As Fig. 4b depicts, the cross sectional velocity profile has two maxima at approximately the jet half-width (symmetric around the centerline ) in the near nozzle region, which collapse to one further downstream.
Near the nozzle, the turbulence intensity is mainly determined by the production term due to the strong velocity gradient concentrated around the jet half-width. This is manifested in the shape of the turbulence intensity profile that shows a maximum located at approximately the jet half-width and very low turbulence intensity around the centerline . Downstream of the nozzle, the transition of the velocity profile results in a shift of the peak of the turbulence production term closer to the jet centerline . Turbulent diffusion in the radial direction and turbulence production by the velocity gradient in the axial direction also contribute to this trend. The turbulence intensity transport equation in a cylindrical coordinate system and for a compressible flow is as follows (Launder and Spalding 1983) where v is the enasmble averaged radial velocity, t the turbulent viscosity, k a model constant, P ij the turbulence production and the turbulence dissipation term. The effect of turbulent viscosity variation on diffusion transport has been neglected in Eq. (20). The turbulence production is given by Watkins (1977) where S ij is the strain rate of the velocity field and D denotes the velocity divergence (socalled dilatation). The production term can be written in cylindrical coordinates as follows Based on order of magnitude analysis for the turbulence terms included in Eq. (20), which is presented in C, and for reduction of the computational cost, the radial convection and axial turbulent diffusion are neglected. In addition, the production term is further simplified by neglecting terms involving the radial velocity component. This results in the following transport equation for the turbulence intensity . 4 a Budget of turbulent kinetic energy in the far field for a jet of air injected into quiescent air (experimental measurements of Panchapakesan and Lumley (1993)). The curves have been made non-dimensional by normalising with cl u 3 cl b 0.5 . b The turbulence intensity radial profile normalised with the centerline velocity at various axial locations x∕d j as measured by Iqbal and Thomas (2007) with the simplified production term P ′ ij given by Subsituting the expression for the velocity profile (Eq. 1) yields In order to compute the production term, the eddy viscosity in every jet location needs to be determined. With the Prandtl mixing length model (Pope 2001) where C t = 0.0161. For unsteady jets, it was empirically found that a better approximation of the turbulent viscosity is obtained if, in Eq. (26), the nozzle velocity is replaced with the centerline velocity u cl and the nozzle radius with the jet half-width b 0.5 . Experimental measurements of Hussein et al. (1994) for the normalised turbulent eddy viscosity ( t ∕(u cl b 0.5 ) ) suggest that it varies radially. For this reason, the turbulent viscosity profile is approximated by where The exponent of was empirically set to 2.5a in order to obtain a good agreement with the experimental measurements of Hussein et al. (1994). The dissipation term can be determined based on the integral length scale and the turbulence intensity as follows where C = 0.09 (Pope 2001). The integral length scale in each cross section is approximately proportional to the jet walf-width where C l I is a model parameter, whose value was found to be around C l I ≃ 0.6 (Wygnanski and Fiedler 1969;Bogey and Bailly 2009;Biswas et al. 2016). From the ratio of the turbulence intensity and dissipation rate, the turbulence time t can be determined where the turbulent velocity u ′ is computed assuming isotropic and homogeneous turbulence Ihme and See (2010) provide the following expression for the scalar dissipation rate in the context of LES where Z′′2 is the mixture fraction variance, Δ the LES filter width, C and C u are model parameters, which are dependent on the filter width Δ . If the filter width is set equal to the integral length scale, Eq (33) becomes where the scalar-to-mechanical time scale ratio is set to its standard numerical value C = 2 (Peters 2001;Markides and Mastorakos 2006). The mixture fraction variance Z′′2 is determined by the following transport equation (Spalding 1983;Peters 2001) where P ij,Z is the production term. In the limit of high Reynolds number 3 , Eq. (35) simplifies to The simplified production term for the mixture fraction variance is given by The solution of the turbulence intensity transport equation is based on a decomposition of the total solution procedure in a convection-production ( k cp ) and a diffusion ( k df ) step. Similar to the numerical approach for the solution of the mass and momentum conservation equations, Eq. (23) is radially integrated to yield where m k = ∫ V kdV is the total turbulence intensity inside a cell, ṁ k = ∫ A ukdA is the turbulence intesnity flow through the cell surface and Ṡ k = ∫ V (P � ij − )dV the volume integrated net turbulence production. The mean turbulence intensity inside a cell can be found from and to determine the radial profile of the turbulent kinetic energy in each cross section due to the combined effect of production and dissipation, it is assumed that The final profile of turbulent kinetic energy k df is found from the numerical solution (at each axial location) of the following 1-D partial differential equation (PDE) given the initial condition k = k cp of the previous solution step An implicit integration scheme is employed for the numerical solution of Eq. (41). Equation (36) is solved in a similar manner.

Auto-Ignition Chemistry Modelling
The analytical and detailed treatment of chemistry for turbulent jet ignition requires the solution of as many transport equations as the number of species of the chemical mechanism. Taking into account the difference in time-scales and the high non-linearity of the chemical source term, the small time-step required is expected to result in a significant computational cost even for a reduced chemical mechanism. However, for turbulent jet ignition, a two dimensional manifold can be identified to fully describe the thermochemical state in any jet location (Steinhilber and Maas 2013;Steinhilber 2015;Fischer et al. 2017b), allowing for a substantial reduction of the computation cost, since only two variables need to be determined during the simulation. From these two variables, one is used to characterise the mixing of the two streams (mixture fraction) and the other one the reaction progress (progress variable). The mixing-chemistry interaction can be accounted for with a different combustion model (cf. (Ghorbani et al. 2015;Wu and Ihme 2014)). The mixture fraction Z is defined based on the total enthalpy h, which is conserved during chemical reaction where h pc and h mc is the pre-and main-chamber mixture total enthalpy, respectively. For the present study, it is essential to provide two alternative definitions of the progress variable, whose importance will be demonstrated later. For hydrogen-air mixtures, a progress variable C is defined based on the mass fraction of the H 2 O species ( Y H 2 O ) as follows where Y mc H 2 O,u and Y mc H 2 O,b is the mass fraction of H 2 O in the non-reacted and fully burned main chamber mixture, respectively. An alternative definition of the progress variable, which depends on the mixture fraction, is the following is the mass fraction of H 2 O for a non-reacted and fully burned mixture at a given mixture fraction Z. For methane-air mixtures, Eqs. (42), (43) and (44)  The solution of Eq. (51) or (49) requires a closure of the mean reaction rate. One possibility is the use of presumed PDF for the mixture fraction and progress variable, as in (Ihme and See 2010). This approach requires an additional transport equation for the progress variable variance C′′2 . Additionally, Eq. (51) is a two dimensional PDE, since it involves a turbulent diffusion term on the right hand side and, therefore, its numerical solution can drastically increase the model computational cost. In order to determine the combined effect of reaction and turbulent diffusion, while keeping the computational cost low, it is assumed that each point of the simulation domain can be mapped to a 1-D reaction-diffusion mixing layer based on the local Z and Z . A mixing layer is characterised by a maximum scalar dissipation Z,max (the scalar dissipation rate at Z = 1) and the thermochemical conditions of the pre and main chamber boundary, at Z = 1 and Z = 0, respectively. The mixing rate varies across the mixture fraction space Z. In the jet, it is higher close to the nozzle and the jet centerline but drops towards the jet side and further downstream the nozzle. To capture this behaviour, it is assumed that the mixing rate varies linearly across the mixing layer, from Z = Z,max at Z = 1 to Z = 0 at Z = 0 (Sadanandan et al. 2007). The mixing layer that an individual point of the jet belongs to is uniquely determined from Z,max = Z ∕Z.
For a given fuel, pressure, temperature and equivalence ratio for the pre and main chamber, the initial thermochemical conditions of the two streams are determined. These mixtures define the fixed boundary conditions at Z = 0 (main chamber mixture) and Z = 1 (prechamber mixture) for the 1-D unsteady PDEs and where Le j is the Lewis number of species j. The initial condition for Eqs. (52) and (53) is obtained by pure mixing of the two streams.
The net effect of reaction and turbulent diffusion (right hand-side of Eq. (51)) can be determined from the transient solution of the 1-D mixing layer PDEs. Equation (51) can be written as where Ṡ ΔC is the source term that contains contributions from reaction and turbulent diffusion. The methodology to derive Ṡ ΔC from 1-D mixing layer calculations is detailed after the presentation of the chemistry modelling.
In order to illustrate the ignition dynamics, a zero-dimensional (0-D) constant pressure reactor is considered at the thermochemical conditions of Table 2. Chemistry is modelled with the full GRI-3.0 mechanism (Smith et al. 1999) that contains 53 species and 325 elementary reactions for both methane and hydrogen. The Cantera plugin for Matlab/Simulink (Goodwin and Moffat 2001) was employed to obtain the chemical reactions and the species thermodynamic properties. While no significant heat is released during the ignition delay time, chemical reactions do take place and a pool of radicals, which trigger the ignition event, is formed. Fig. 5a depicts the evolution of the important species mass fractions (normalised with the corresponding maximum concentration) for a reactive gas of Z = C = 0.60. It can be seen that short-lived species such as HO2 reach a peak in concentration well before substantial amount of fuel is consumed. The time evolution of the progress variable for a mixture of varying mixture fraction Z and initial progress variable C = Z is given in Fig. 5b. The mixture fraction is varied from Z = 0.00 to 0.90 in steps of 0.05, with the lines of increasing intensity corresponding to higher mixture fraction. The time instant of ignition, determined with the maximum of the mixture fraction time derivative, is shown with red dot in Fig. 5b. It can be easily seen that the progress variable varies only minor before ignition; this also justifies assumption 7 employed for the model development.
To obtain the source term Ṡ ΔC , Eqs. (52) and (53) are numerically solved from the initial condition (mixing line), and the progress variable is computed from Eq. (43). For the conditions of Table 2, this results in a different progress variable distribution for the different time instances as shown in Fig. 6a (left). The term Ṡ ΔC at a given time instant t n can then be easily obtained from and following transformed to Ṡ ΔC =Ṡ ΔC (C, Z; Z,max ) . The above calculations are repeated for different mixing rates Z,max and stored in tables.
The source term varies only minor for a small to moderate mixing rate, but drops drastically very close to the extinction limit. For this reason, the specified mixing rates are not equally distributed, but a more dense discretisation is employed close to the extinction mixing rate. The extinction mixing rate is a function of the fuel and the thermochemical conditions. Mixtures of higher reactivity (hydrogen fuel or high temperature) are more resistant to extinction and the corresponding extinction mixing rate is higher. Figure 6b provides the combined reaction-diffusion source term Ṡ ΔC in the mixture fraction-progress variable space (Z,C) at the conditions given in Table 2 for three levels of . 6 a Progress variable distribution (left) and time derivative of the progress variable (right) over the mixture fraction space for different time instances (increasing color intensity with time). b Contour plot of the combined reaction-diffusion term as a function of mixture fraction Z and progress variable C for three levels of the maximum scalar dissipation rate Z,max (a)

(b)
the maximum scalar dissipation rates Z,max . Increase of Z,max results in a narrower range (in the (Z,C) space) of sustainable chemical reaction, as well as an overall decrease of the progress variable production rate. This phenomenon is indicative of the quenching effect of a very high mixing rate, as observed in both Sadanandan et al. (2007) and Sidey and Mastorakos (2018). For the thermochemical conditions of Table 2 no ignition is observed when Z,max > 1500 1/s.

Model Reduction
In the previous sections, a computationally efficient 1-D model for TJI has been derived. The computational cost of the model for simulating the ignition by a hot turbulent jet with 1 ms duration is approximately a few minutes; ranging from 5 to 10 min depending on the exact thermochemical conditions. While its computational cost is significantly lower when compared to LES, it can be considerably reduced if the model is further simplified. The motivation for this reduction is the integration of the model in a quasi-dimensional engine simulation code for prechamber combustion with runtime requirement in the order of a few seconds.
To determine the ignition delay with the simplified model, the evolution of the ignition precursor is monitored along Lagrangian particles that are imposed at the nozzle exit at the timing of the exit of the reactive jet as Fig. 7 depicts. The axial position of each particle p is controlled via the following equation where the particle velocity u(x p,n , r p,n ) at a given location (x p,n , r p,n ) is determined from the jet velocity distribution u = u(x, r) . The expression for the jet velocity distribution is derived via the steady state solution of the momentum conservation equation in 4.
The turbulent velocity and the mixture fraction variance for the simplified jet model are computed based on the Prandtl mixing-length hypothesis (Pope 2001;Torrez et al. 2011) as follows and (56) x p,n+1 = x p,n + u(x p,n , r p,n )dt Fig. 7 Schematic of the jet with the Lagrangian particles trajectories (left) and the hypothetical evolution of different variables along a particle (right) where C u ′ and CZ ′′2 are tuning parameters. The scalar dissipation rate is computed by combining Eqs. (30), (31), (34), (57) and (58) The evolution of the variable ΔC p , through which the ignition progress variable C ign,p = ΔC p 1−Z p can be determined, is given by where the source term is a function of the local progress variable C p,n = Z p,n + ΔC p,n , mixture fraction Z p,n and scalar dissipation rate Z,p,n . The mixture fraction distribution Z(x, r) is computed assuming that Table 3 summarises the model parameters, their values, and the corresponding reference for those obtained from the literature. The parameter C u ′ of the simplified turbulence model is not directly obtained from the literature but it is adjusted such that a good agreement is obtained with the experimental measurements of Crow and Champagne (1971). The parameters CZ ′′2 and C ig,cr are calibrated to match the ignition delay time for the largest nozzle diameter of the simulations performed by Ghorbani et al. (2014).

Results and Discussion
In the problem of transient hot ignition, the flow field determines the mixture fraction and the mixing rate distribution, both of which affect the ignition precursor source term of Eq. (54) and, therefore, the total ignition delay time. For this reason, before proceeding with the validation of the ignition model, the soundness of the flow-related sub-models (60) ΔC p,n+1 = ΔC p,n +Ṡ ΔC (Z p,n + ΔC p,n , Z p,n , Z,p,n )dt is verified. A combination of data from the literature in canonical configurations and the optical prechamber is used. In this way, the model is validated with high quality data (low uncertainty and well established boundary conditions) but also under realistic turbulent jet ignition conditions. The following subsection begins with the validation of the 1-D model and proceeds with the validation of the 0-D model. Table 4 reports the conditions where the jet flow and turbulence sub-models are validated. Experimental hot-wire measurements by Crow and Champagne (1971) Table 4 (CRO1). The velocity profile gradually transitions from almost flat at x∕d j = 0.025 to approximately Gaussian downstream of the nozzle and within the self-similarity region. This behaviour is a result of the variation of the exponent a of the velocity profile expression, which is determined via Eq. (19). The exponent a decreases downstream the nozzle and reaches a = 1.5 in the self-similarity region, as shown in Fig. 8b. The model tends to slightly under-predict the velocity around the jet half-width for x∕d j = 2 and 4 . This discrepancy is linked with a lower value of a than the experimental one.

Flow and Turbulent Field
The value of a obtained via the polynomial expression provided in Musculus and Kattke (2009) is lower than the one derived via Eqs. (19) or (73). Therefore, the use of Musculus and Kattke (2009) equation for the shape factor a is expected to result in greater deviation from the experimental velocity profile. This highlights the importance of accounting for the main chamber mixture contribution in the overall jet momentum, which in comparison to Diesel jets can be significant in TJI applications. In addition, it is observed that the approximate density profile obtained from the polynomial Eq. (73) agrees well with the one derived from Eq. (19).
A steady state jet with an abruptly reducing inflow velocity is considered for the validation of the 1-D jet model. Fig. 9 compares the centerline profile of the axial velocity as predicted by the 1-D model against the DNS of Shin et al. (2017) at various time-instances after the nozzle velocity drops to zero. The steady-state jet profile predicted by the 1-D model agrees well with the one from the DNS. Additionally, the axial location where the centerline velocity of the decelerating jet approaches the velocity of the steady jet is well captured for all time instances. Nonetheless, the model tends to over-predict the centerline velocity in the near nozzle region for t∕ j ≥ 20 . This deviation can be attributed to the fact that the exponent a is determined in a steady-state manner without accounting for the actual nozzle flow. However, shortly after the nozzle velocity becomes zero, the jet core velocity drops, resulting in an overall modification of the velocity profile in the near-nozzle region and a faster deceleration of the velocity than the one predicted by the 1-D model. This reasoning is further supported by the fact that the model error is mainly concentrated up to Table 4 Operating conditions for the validation of the flow and turbulence jet sub-models a distance x∕d j < 5 , which is approximately the extent of the jet core. Neal and Rothamer (2017) proposed a dynamic equation for the determination of the shape factor, which can potentially offer better agreement with the DNS results in the near nozzle region of the jet. The source term is dependent on the scalar dissipation rate, which is linked to the turbulent mixing frequency computed from Eq. (31). To validate the proposed turbulence model we make use of available literature data. Figure 10a compares the turbulence intensity profile normalised by the centerline velocity as obtained from the 1-D model in various crosssections against the experimental measurements of Crow and Champagne (1971). While the modelled profile deviates from the experimental data, the peak turbulence location correctly shifts towards the centerline downstream of the nozzle. A strong deviation is evident just after the nozzle exit ( x∕d j = 0.025 ), where the 1-D model significantly overestimates Fig. 8 a Axial-velocity profiles obtained with the 1-D model (blue solid line) versus the measurements (red circles) from Crow and Champagne (1971) in the developing region for condition CRO1 of Table 4. For clarity, the data and curves for each successive measurement location are offset vertically by 0.5u∕u cl ; b 0.2 is the radial position from the jet centerline where the axial velocity is 20 % of its center- the turbulence intensity. It is well known that the Boussinesq approximation does not hold near the nozzle, since in reality, turbulent stresses do not develop instantaneously under the presence of flow shear. Additionally, the turbulent dissipation rate is calculated only approximately, using an imposed distribution for the integral length scale. A more accurate prediction could be obtained by solving an additional transport equation for the turbulence dissipation rate, as in the k − model. Figure 10b is a contour plot of the turbulence intensity normalised with the nozzle velocity squared. Apart from the near-nozzle area, where the model over-predicts the turbulence intensity, there is a fair qualitative agreement with the experimental measurements of Georgiadis and DeBonis (2006).
Direct numerical simulation for an unsteady jet with the nozzle velocity and mixture fraction profile of Fig. 11 (top) are employed for the validation of the penetration correlation. Initially a non reacted mixture of air and methane flows from the nozzle, which switches to adiabatic combustion products at 400 s . As can be seen in Fig. 11 (bottom) the reactive jet penetration is very well captured by the 0-D model.
Figures 12a depict the reactive jet penetration as a function of time (the time reference corresponds to the flame front exit timing) from the 0-D (solid lines) and the experiments in the optical prechamber (circles) for five operating conditions given in Table 4. The jet penetration is well predicted by the zero-dimensional model for all conditions. For the larger nozzle, the jet centerline velocity predicted by the 0-D model agrees well with the experimental measurements (Fig. 12b, left). The experimental data for the smaller nozzle (Fig. 12b, right) are noisy and do not allow a precise quantitative assessment of the model performance. However, good qualitative agreement can be seen in the far field.

Ignition Dynamics
In the final part of this section, the simplified 0-D ignition model is validated with the 3-D simulation results of Ghorbani et al. (2014) and the reactive experiments in the optically accessible prechamber. Table 5 provides the operating conditions for the validation of the ignition model, which include both steady and unsteady jets, as well as two different fuels. For the first four validation cases, the combustion products are fully burnt and their composition is the same as of the main chamber (Ghorbani et al. 2014);  Figure 13 depicts the variation of the total jet ignition delay time (left) and location (right) with the nozzle diameter for a jet of partially cooled combustion products issuing into a hydrogen-air mixture as obtained from the 0-D model (blue circles) compared to the simulations of Ghorbani et al. (2014) (red crosses). A good agreement between the simulation results of Ghorbani et al. (2014) and the 0-D model is obtained. Reduction of the nozzle diameter results in a rise of the induction time due to the increase in the mixing rate that inhibits the chemical reactions. The dependence of the ignition delay time on the nozzle diameter and the mixing rate is highly non-linear. Whereas, for large nozzle diameters, mixing has a rather limited effect on the ignition delay time, its importance increases significantly close to the extinction limit. This trend is well captured by the 0-D model, partly due to the dependency of the chemical source term on the mixing rate.
The ignition delay time predicted by the 0-D model (red stars) for the OPC conditions is depicted in Fig. 14, and compared to the one obtained by the main chamber pressure trace (blue circles) and the OH* chemiluminescence imaging (green circles). A very good agreement is obtained for all the conditions between the experimentally determined ignition delays and the ones from the 0-D model. The ignition delay as obtained by the time instant that the main chamber pressure begins to rise is higher than the one which is derived by the OH* chemiluminescence imaging. This is mainly a result of the time required for the ignition event to result in heat release and temperature rise.
An increase of the mixture temperature is associated with a reduction of the ignition delay time, due to the higher mixture reactivity. Despite the fact that higher temperature results also in more intense mixing, since the jet generated from the prechamber combustion is stronger, the effect of temperature on ignition chemistry is more dominant that the one on turbulent mixing. Similarly to the steady state simulations of Ghorbani et al. (2014), the reduction of nozzle diameter leads to a considerable increase of the total jet ignition delay time under the same thermochemical conditions. It can be also observed that the  Table 5 Operating conditions for the validation of the ignition sub-model Name Fuel ignition delay is significantly longer for methane-air in comparison to hydrogen-air mixtures, which is mainly linked with the slower methane ignition chemistry and the lower extinction mixing rate for methane-air mixtures.

Conclusions
In this work, a novel 1-D and 0-D jet model to predict the ignition delay by a hot turbulent jet issuing into a reactive mixture was developed. The basis of the model is the cross section integrated mass and momentum conversation equations. An extension of the jet model with transport equations for the turbulence intensity and the ignition precursor was employed. The chemical source term was obtained from 1-D flamelet calculations and tabulated during the main simulation loop. The model was hierarchically validated with a combination of data obtained from the literature and an optically accessible prechamber set-up. Results verified the soundness of the developed flow, turbulence and ignition sub-models with the accurate prediction of various metrics including the axial and radial velocity profiles, the jet penetration, the turbulence intensity, the ignition delay time and location. Results suggest that the ambient mixture density has an important effect on the jet penetration and, therefore, cannot be neglected when the density of the issuing and ambient fluids is comparable. In addition, the turbulence is dominated by the production terms of the transport equation, and primarily the radial velocity gradient. Preserving only the leading order terms allows the transport equation for turbulence to be significantly simplified, with a considerable reduction of the computational cost. Finally, the chemical source term as obtained by the 1-D flamelets can sufficiently account for the turbulent mixing and ignition chemistry interaction, which results in the variation of the ignition delay with the flow and thermochemical conditions.
The present model constitutes a valuable tool in the design of combustion systems employing turbulent jet ignition, which allows a fast and accurate assessment of the design and operating condition effect on the ignition propensity and timing. Future efforts will be directed towards the use of the present model as an analysis tool to interpret the different combustion modes during turbulent jet ignition. and applying the integral transformation with the variable = r∕r O yields For the flow terms it is derived in a similar manner The functions F, G, H and K are pre-tabulated before the main simulation loop for the known conditions (Sc and t r ) with Z cl ∈ [0, 1] and a ∈ [1.5, 500].

Appendix 2: Shape Factor in Developing Region
A steady jet is considered of velocity u j and density pc issuing into a quiescent mixture of density mc from a nozzle of area A NZ . The jet momentum is conserved and, therefore, the momentum flow rate in each radial cross-section is equal the one of the nozzle Additionally, within the potential core, the centerline velocity is equal to the one of the nozzle ( u cl = u j ) and the centerline mixture fraction is one ( Z cl = 1 ). Therefore, the velocity and mixture fraction profiles are given by and Substituing the mixture fraction, velocity and density profiles into Eq. (68) and applying the integral transformation with the variable = r∕r O results in the final non-linear integral equation for the determination of the shape factor a While the integral of Eq. (71) cannot by analytically determined, the computational cost of numerically solving Eq. (71) is minimal, since it is only carried out once before the beginning of the main-simulation loop. A polynomial expression for the unknown value of a can be nevertheless obtained, if instead of the actual density profile of Eq. (3) the following approximate linear profile is employed and a Schmidt number of unity is assumed ( Sc = 1 ). This results in the following polynomial expression through which the shape factor a can be determined which is similar to the one derived by Musculus and Kattke (2009).

Appendix 3: Order Analysis Of Turbulence Terms
The axial and radial velocity components in the far field of a round turbulent steady-state gas jet can be found through the solution of the momentum conservation for a self-similar velocity profile (Schlichting and Gersten 2016). Introducing the entrainment constant and the non-dimensional coordinate the axial and radial velocity component are given by The turbulent viscosity can be found from where C t = 0.0161 . The following approximate turbulent velocity profile where u ′ CL is the centerline turbulent velocity, has been determined based on a simple fit of the experimental measurements performed by Launder et al. (1973). Based on Launder et al. (1973); Lockwood and Naguib (1975) the centerline turbulent velocity in a given cross section is approximately 25 % of the centerline velocity in the far field; that is u � = 0.25u � CL .
(72) = Z pc + (1 − Z) mc = pc t r + (1 − t r )Z All terms can be easily found by combining Eqs. (74-83) and obtaining the corresponding derivatives. Normalisation of the various terms with u 3 cl ∕b 0.5 makes them independent of the axial location x, nozzle diameter d j and jet velocity u j . Figure 15a shows the terms of the k-equation (after they were made non-dimensional) divided by the maximum of the turbulence production term. It can be seen that the production term is the most dominant one, apart from very close to the jet centerline where axial convection is the leading order term. Dissipation of turbulent kinetic energy and radial diffusion are comparable in magnitude to the production term. Axial diffusion and radial convection are about one order of magnitude smaller than the production and dissipation terms and, for this reason, they can be neglected. The different terms contributing to the turbulence production are reported in Fig. 15b. It is apparent that the term P ur , associated with production of turbulent kinetic energy by the radial gradient of the axial velocity, is the dominant one around the jet half-width. Close to the jet centerline where the axial velocity component is significant, the production of turbulent kinetic energy by the axial gradient of the axial velocity component (term P ux ) becomes important.
Based on the order analysis of the terms, the turbulent intensity equation simplifies to the following equaion 1 3 which can be iteratively solved over the unknown y for each jet location x d j given j and t r .
The exponent a is determined from Eq. (19). The jet axial velocity at position (x, r) is given by Figure 16a shows the centerline velocity profile normalised with the nozzle velocity as obtained from the numerical solution of Eq. (87) (dashed lines) compared to the expression provided by Musculus and Kattke (2009) (dashed dotted lines) and Abraham and Bracco (1993) (solid lines), for a temperature ratio t r = 10 (red), t r = 5 (magenta), t r = 1 (blue) and t r = 0.1 (green); both x and y axis are in logarithmic scale. Reduction of the temperature (89) u(x, r) = u j ⋅ y(x) ⋅ 1 − r∕r O (x) a(x) 2

Fig. 16
Centerline velocity profile normalised with the nozzle velocity versus the downstream location x (normalised with the nozzle diameter) (a) and jet penetration normalised with the nozzle diameter versus time (normalised with the jet time tau j = d j /u j ) (b) for a steady state jet from Eq. (87) (dashed lines), Musculus and Kattke (2009) (dashed dotted lines) and Abraham and Bracco (1993) (solid lines), for a temperature ratio t r = 10 (red), t r = 5 (magenta), t r = 1 (blue) and t r = 0.1 (green) (a) (b) 1 3 ratio (or equivalently increase of the injected fluid density in comparison to the ambient) results in an increase of the potential core length. In the far field, the axial centerline velocity scales approximately with ∼ 1∕x for all correlations. However, differences between the various models are evident. For the high temperature ratios (corresponding to TJI relevant conditions), the correlation of the present study shows the higher penetration rate (Fig. 16b) compared to Musculus and Kattke (2009). This is mainly associated with higher value of exponent a in comparison to the correlation of Musculus and Kattke (2009). When the injected fluid has higher density than the ambient fluid, this trend is reversed and the correlation of Musculus and Kattke (2009) shows the highest penetration rate.