Numerical Study of Hydrogen Auto-Ignition Process in an Isotropic and Anisotropic Turbulent Field

: The physical mechanisms underlying the dynamics of the ﬂame kernel in stationary isotropic and anisotropic turbulent ﬁeld are studied using large eddy simulations (LES) combined with a pdf approach method for the combustion model closure. Special attention is given to the ignition scenario, ignition delay, size and shape of the ﬂame kernel among different turbulent regimes. Different stages of ignition are analysed for various levels of the initial velocity ﬂuctuations and turbulence length scales. Impact of these parameters is found small for the ignition delay time but turns out to be signiﬁcant during the ﬂame kernel propagation phase and persists up to the stabilisation stage. In general, it is found that in the isotropic conditions, the ﬂame growth and the rise of the maximum temperature in the domain are more dependent on the initial ﬂuctuations level and the length scales. In the anisotropic regimes, these parameters have a substantial inﬂuence on the ﬂame only during the initial phase of its development.


Introduction
Turbulent combustion processes exist widely in various technical applications, e.g., gas turbines combustion chambers, burners, car engines. Numerical simulations of turbulent reactive flows are one of the most challenging tasks in contemporary CFD (Computational Fluid Dynamics) both in the academic and industrial science centres. The turbulent flame being the result of the two-way interaction between chemical reactions and turbulent structures is characterized by significant oscillations of all flow variables. The combustion strongly changes the density and viscosity associated with temperature variations and these variables strongly affect the velocity field. The opposite impact is equally important, i.e., the turbulence affects the combustion process by enhancing the convective and diffusive mixing mechanisms. In this work we concentrate on the auto-ignition phenomenon for which the mixing process is particularly important, especially in non-premixed conditions where the regimes favorable for ignition as well as an initial flame development are conditioned by the mixing intensity [1,2]. Modeling the ignition process is one of the most difficult tasks in CFD research of reactive flows. It requires the use of Direct Numerical Simulations (DNS) or Large Eddy Simulations (LES) with advanced combustion models, e.g., in non-premixed conditions the Conditional Moment Closure (CMC) [3][4][5], Eulerian stochastic fields (ESF) [6][7][8] and Partially Stirred Reactor (PaSR) type models [9,10] have been proven to predict the ignition process and flame stabilization relatively well. Recently, LES method is used in investigations of very complex combustion problems including the flame initiation by a spark ignition and auto-ignition in gasoline engines [11][12][13] and combustion instabilities (e.g., thermoacoustic oscillations, blow-off) in liquid-fuelled burners [14,15] or gas turbines [16,17].
In the present work we apply LES to study a very fundamental issue of the autoignition process, i.e., we analyze to what extent this phenomenon is influenced by the isotropy of turbulent velocity field. Homogeneous isotropic turbulence (HIT) is an idealized model of turbulent flow that has attracted the attention of researchers for many years. A relative simplicity of this configuration caused that for a long time it has been also the focus of research related to fundamental problems of auto and forced ignition phenomena. A prominent example of such investigations is a study of Mastorakos et al. [18] who performed DNS of the auto-ignition process in a 2D mixing layer between methane and air. The most important findings of this research were that the ignition is a localized phenomenon and that turbulent flows ignite earlier than laminar ones. A DNS study of a 2D hydrogen/air mixing layer performed later by Im et al. [19] showed that, on average, when a characteristic time-scale of turbulence (τ t ) was shorter than a laminar ignition time (τ ign,l ), i.e., τ t /τ ign,l < 1, the ignition was delayed, whereas for τ t /τ ign,l > 1 it occurred earlier. More recently, a 3D hydrogen/air configuration was analyzed by Doom and Mahesh [20]. Their DNS results showed that when the flame kernels collide they do not always extinguish as reported in [18] for 2D cases. It was observed that the flame kernels could also combine and spread with/without increasing flame intensity. Compared to the auto-ignition phenomenon a situation in spark-ignition problems seems even more complex. Here, the born flames can have very complex structures, e.g., a tribrachial one as shown by Chakraborty and Mastorakos [21], or their occurrence can be significantly conditioned by shear stresses, which arise between flow streams [22]. The presence of strong shear stresses causes that an initially assumed isotropic field evolves and becomes highly anisotropic leading to significantly different regimes than in the HIT conditions.
Although in the real world anisotropic turbulent flows are common the fundamental research devoted to them is rare and definitively less advanced compared to the HIT studies. The anisotropic conditions occur for instance in shear layers of jet flows, however, in these cases the auto-ignition process is studied mainly from the point of view of fuel/oxidizer parameters including their compositions, temperatures, speeds or turbulence intensities [23][24][25]. Knowing the fact that the ignition process originates in small turbulent scales one should assume that their anisotropy can affect both the flame kernel occurrence as well as its further development. The objective of this paper is to investigate a relation between an anisotropic flow field and the nature of the auto-ignition process. We consider the homogeneous anisotropic turbulence (HAT) with various turbulent intensities and length scales. The flow configuration constitutes an academic type problem for which there are no experimental data and the investigations have an exploratory character. They are performed using LES combined with the Probability Density Function (PDF) approach based on the Eulerian stochastic fields (ESF) method. The LES-ESF approach has been proved to work well for ignition problems, e.g., in the study on the auto-ignition of a hydrogen jet performed by Jones and Navarro-Martinez [8].

Mathematical Formulation
In the LES framework, the reacting low speed flows are governed by the conservation equations of mass, momentum and energy given as: where the bar and tilde symbols denote spatially filtered and Favre-filtered quantities [26] respectively. In these equations, u represents the velocity vector, ρ is the density, p is the hydrodynamic pressure, Y α represents the mass fractions of species α = 1, . . . , N and h denotes the total enthalpy. The symbols τ and D α = µ/ρSc, D = µ/ρPr represent the viscous stress tensor and the mass and heat diffusivities where Sc = 0.7, Pr = 0.7 are the Schmidt and Prandtl number. The sub-grid scale stress tensor is computed via an eddy viscosity concept as τ SGS = 2µ t S, where S is the strain rate of the resolved velocity field and µ t is the sub-grid viscosity computed using the model of Vreman [27]. The terms D SGS α = µ t /ρSc t , D SGS = µ t /ρPr t , are the sub-grid diffusivities where Sc t , Pr t are the turbulent Schmidt and Prandtl numbers, both set to 0.9 [28]. The turbulent Schmidt number quantify the ratio of rates of the momentum and mass transfer due to the eddy motion whereas the turbulent Prandtl number defines the relation between the momentum eddy diffusivity and heat flux produced by turbulent eddies. Equations (1)-(4) are complemented with the equation of state p 0 = ρR T, where p 0 , R and T are the thermodynamic pressure, gas constant and temperature, respectively.
The chemical source termsẇ α represent the net rate of formation and consumption of the species. A highly non-linear nature of this term means that sub-grid fluctuations cannot be ignored. In the present work the scalar Equations (3) and (4) are replaced by an equivalent evolution equation for the density-weighted filtered PDF function, which is solved using the stochastic field method proposed by Valiño [6]. Each scalarφ α is represented by 1 ≤ n ≤ N s stochastic fields ξ n α such that The stochastic fields evolve according to: where the total diffusion coefficient is defined as Γ = D α + D SGS α , the micro-mixing time scale equals to τ =ρ∆ 2 /(µ + µ t ) [8] with ∆ = V 1/3 cell being the LES filter width and dW represents a vector of Wiener process increments different for each field.

Numerical Method
The computations have been performed using an in-house LES code SAILOR based on the low Mach number approach [29]. The filtered transport equations are solved on a Cartesian half-staggered mesh [30][31][32]. The time integration is performed by means of an operator splitting strategy where the transport in physical space and chemical terms are solved separately. The convection and diffusion terms of the governing equations are calculated using a predictor-corrector approach with help of the 2nd order methods, Adams-Bashforth and Adams-Moulton, respectively. For the time integration of term representing the chemical source the VODPK (Variable-coefficient Ordinary Differential equation solver with the Preconditioned Krylov method) solver [33] is applied. The spatial discretization is performed by the 6th order compact difference method for the Navier-Stokes and continuity equations. The transport equations of scalars, i.e., species mass fraction and total enthalpy are discretized using 5th order WENO (Weighted Essentially Non-Oscillatory) [34] scheme. The chemical reactions are computed using CHEMKIN interpreter [35]. There are a number of studies that evaluate hydrogen combustion using different chemical kinetics. Recently, Bazooyar and Darabkhani [36] compared the statistical errors of the combustion mechanisms in a numerical study of premixed hydrogen turbulent flame in a thermo-photovoltaic combustor. In the present investigation, we apply a detailed hydrogen oxidation mechanism proposed by Mueller et al. [37] with 9 species and 21 reactions.
The SAILOR code has been thoroughly verified in previous studies devoted to nonreactive flows [38,39] and number of combustion problems in jet flows [23,40,41] and mixing layers [22,42,43]. In all these cases it turned out to be very accurate. Moreover, thanks to the applied high-order discretization schemes grid-independent results could be obtained on relatively coarse meshes. The same is expected to carry over to the simulations reported in the present work, as the complexity of the cases analyzed here is similar as in the problems analyzed previously. However, to verify the sensitivity of the results to the settings of the combustion model, in the selected cases the simulations have been performed applying 4 and 8 stochastic fields.

Problem Description
We consider an auto-ignition problem in a stationary mixing layer formed between a hydrogen/nitrogen mixture (0.1 H 2 / 0.9 N 2 ) at 300 K and hot air (1000 K). Figure 1 shows the computational domain (L x × L y × L z = 0.025 × 0.04 × 0.025 m 3 ) in which in the vicinity of its centre (y = 0) isotropic (HIT) or anisotropic (HAT) initial velocity fluctuations (u ) are generated. The iso-surfaces presented in Figure 1 correspond to sample initial velocity fluctuations in the 'x' direction equal to −0.5 m/s (light grey) and 0.5 m/s (dark grey). The flow is periodic in the 'x' and 'z' directions and the no-slip walls are assumed at y = ±L y /2. The fuel concentration is given by C(y) = (tanh(2y/δ) + 1)/2, where δ denotes the initial vorticity thickness assumed to be equal 0.5 mm. ; initial temperature (contours at the side boundary) and anisotropic turbulent field (iso-surfaces of u 1 ).

Initial Conditions
The turbulent velocity fields are generated by a method proposed by Orszag S. [44] that was originally used for the generation of the HIT fields. In the first step it relies on generation of an initial random velocity fluctuation (u i,ran ) where the subscript 'i' stands for the 'x', 'y', 'z' direction. These random fields have the zero mean value and assumed RMS level. Then, they are transformed into the Fourier space where the particular velocity components are mutually correlated (u i,ran → u i ) according to an assumed energy spectrum and requiring that ∇ · u = 0. Due to these operations the RMS values of u i,ran and u i are slightly different. In the present research we extend this procedure such that both HIT and HAT initial velocity fields could be generated. In the region of the mixing layer ±10δ the energy of the initial velocity is defined through the Passot-Pouquet kinetic energy spectrum defined as [45]: where k represents the wave number, k i,0 is an adjustable wave number allowing generation of the velocity fluctuations with a required Taylor length scale λ i for a particular component of the velocity fluctuations. Their global RMS value is defined as u RMS = ∑ 3 i=1 u i,RMS /3. In the HIT conditions the Taylor length scale is defined as λ = u 2 RMS /(∂u RMS /∂x) 2 where u RMS is arbitrary velocity component and 'x'-arbitrary direction. In the case of HAT the global Taylor length scale is defined, as λ = ∑ 3 i=1 λ i /3 with λ i referring to particular directions. The initial HIT and HAT velocity fields are obtained by specifying either the same or different u i,RMS and λ i in particular direction.
The turbulent Reynolds numbers can be computed either based on u i,RMS and λ i as Re λ i = u i,RMS λ i /ν, with ν-the air kinematic viscosity, or based on the turbulent quantities that can be treated as representative for a given test case (u RMS , λ), i.e., Re λ = u RMS λ/ν. In the present work we consider four HIT and four HAT cases with Re λ ≈ 50, 100, 150, 200 of which all parameters are given in Table 1. Note, that due to above-mentioned procedure (u i,ran → u i ) the final values of Re λ slightly diverge from the assumed one. The last column shows the turbulent time scale defined as τ t = λ/u RMS to which we will refer discussing the impact of τ t /τ ign,l on the ignition process. For HAT cases τ t remains on an almost constant level regardless of Re λ , whereas, in the HIT conditions its values pronouncedly depend on Re λ .

HIT
Case  Moreover, it can be seen that the impact of Re λ on the initial solution is large. In the case HAT-050 where λ i are small the size of the turbulent structure is also small compared to the case HAT-200. As will be shown later these differences have a significant influence on flame development. At this point, one should mention that in all HIT and HAT cases the length scales λ i are approximately ten times smaller than the sizes of the computational domain. Hence, the impact of the boundary conditions (walls, periodicity) on the evolving flow structures is very small and does not affect the ignition phenomenon.

Computational Details
All the computations have been performed on the mesh consisting of 192 × 128 × 192. In Figure 1 every four mesh line is shown. The mesh nodes are distributed uniformly in the periodic directions (∆x = ∆z = 1.31 × 10 −4 m) and are compacted in the 'y' direction. The minimum cell size is equal to ∆y min = 5.59 × 10 −5 m and its average value in the region −10δ ≤ y ≤ +10δ is ∆y −10δ,+10δ = 1.11 × 10 −4 m. The simulations were conducted with a constant time step ∆t = 1 × 10 −6 s corresponding to the CFL number approximately equal to 0.5, and lasted for the time period allowing the flame development across the whole domain.
When the simulations start the flow evolves and the initially imposed turbulent structures combine and develop. In LES, they are accurately represented on the numerical mesh provided that the Pope's criterion is fulfilled [46], i.e., the ratio of the sub-grid (k SGS ) to the total kinetic energy (k TOT ) should be less than 0.2. The former is computed as k SGS = µ 2 t /(ρC∆) 2 (C = 0.1 [47]). Figure 4 shows the contours of k SGS /k TOT at the time instances before and during the ignition phase (t ≈ 1.0 ms). It can be seen that the values of k SGS /k TOT are very small and there are only a few places where they reach the level of 0.1. Additionally, in this figure we show the isolines of the ratio of the sub-grid to the molecular viscosity (µ t /µ). The values of µ t are approximately two orders of magnitude smaller than µ. All this means that the impact of the modeled sub-grid scale on the solution is very small and the computational mesh is appropriate for the present problem. Nevertheless, to show the level of dependence of the results on the mesh density, we performed the test computations on the mesh consisting of 70% more nodes, i.e., 224 × 160 × 224 nodes.

Auto-Ignition Scenario and Initial Flame Development Phase
After the simulation starts the fuel mixes with the hot air. The strength of the mixing process depends on the initial flow field but in all analyzed cases it leads to the conditions favorable for the auto-ignition, which is discussed below for the case HAT-200. In agreement with the findings of Mastorakos et al. [18] the ignition spots have the highest probability to appear in the places with the mixture composition corresponding to the most reactive mixture fraction (ξ MR = 0.0178) and where the scalar dissipation rate is low. In this work the mixture fraction ξ is computed following the Bilger's definition: where Y's are the elemental mass fractions, W's are the atomic weights and the subscripts, f and o, correspond to the fuel and air, H and O refer to the hydrogen and oxygen atoms, respectively. The scalar dissipation rate is defined as Figure 5 shows the iso-surface of the most reactive mixture fraction coloured by the scalar dissipation rate at the time moments t = 0.9 ms and t = 1.0 ms from the beginning of the simulation. It can be seen that there are preferred locations (small χ) in which the ignition spots can appear. Comparing the presented solutions, one can notice that during 0.1 ms the flow evolves only slightly, however, in this time the flame kernels significantly grow and propagate towards the stoichiometric mixture (ξ ST = 0.225). At the later time the flame develops along ξ ST and eventually expands across the mixing layer. Figure 6 shows semi-transparent iso-surfaces of three different temperature values and the vorticity module (Ω) contours in the time moments corresponding to the results presented in Figure 5. The high temperature regions seem to be distributed randomly, however, their localisations correlate well with low values of χ (Cf. Figure 5). It can be seen that between t = 0.9 ms and t = 1.0 ms the initially local flames have connected and spread almost along the entire 'x-z' cross-section of the computational domain. The view along the 'x' direction ( Figure 6c,d) shows that the regions of high temperature induce large local velocity gradients that manifest by a significant growth of the vorticity module (white contours). Its contours in the 'x-z' cross-section plane indicated by the green vertical line are shown in Figure 7. The presented results correspond to the time moments t = 0.6 ms, t = 0.9 ms and t = 1.0 ms. It can be seen that before the auto-ignition occurs (Figure 7a) and during its initial phase ( Figure 7b) the vorticity values are comparable and low, in general. However, at the later time ( Figure 7c) the thermal expansion caused by an increasing flame volume accelerates the flow in the layers attached to it. In effect, in a direct flame vicinity the vorticity grows significantly. The contours of its high values (Ω ≈ 1.0) correspond to the places in which the analyzed cross-section plane crosses the high temperature iso-surface.

Analysis of the Flame Volume and Maximum Temperature Growth
The above-depicted auto-ignition process is qualitatively the same in all cases considered. In this section, we analyze how the differences in the initial solutions affect the moment of auto-ignition and the subsequent evolution of the flame kernel. As the reference solution we assume the auto-ignition process in the laminar conditions in which the velocity fluctuations are equal to zero and the mixing between the fuel and air is driven by the diffusion mechanism. The auto-ignition time corresponds to the moment when the maximum temperature in the domain reaches the level of 1% above the oxidizer temperature that is commonly used criterion [25,48]. In the reference laminar case it is equal to τ ign,l = 0.86 ms. Figure 8 shows the temporal evolutions of the mass fraction of HO 2 , the maximum temperature in the computational domain and the flame volume V f for the cases HAT-200 and HIT-200 as the function of a normalised time τ = t/τ ign,l . The radical HO 2 is the so-called preignition species whose appearance precedes the auto-ignition process and the temperature rise. The flame volume is a measure of the flame growth at its later phase. This parameter is defined as the sum of the volume of the computational cells in which the temperature is larger than 1200 K, i.e., V f = ∑ V T>1200K cell . The results presented in Figure 8 were obtained applying 4 and 8 stochastic fields (N s in Equation (5)). It can be seen that their number has only a minor influence on the solutions. The existing differences are very low and only quantitative, which confirms that the applied combustion model provides the results that can be regarded as independent of N s . Regarding the temporal evolution of HO 2 shown in Figure 8a one can see that the turbulence speeds up the preignition phase. A sudden rise of HO 2 around τ = 1.0 is related to the auto-ignition phenomenon accompanied by a steep growth of the maximum temperature (T max ) presented in Figure 8b. A close inspection on the ignition times shows that HIT-200 is characterized by a shorter ignition delay compared to HAT-200, though the differences are very small (≈ 1%τ ign,l ). The exact values are τ = 0.932 for HIT-200 and τ = 0.945 for HAT-200 when 4 stochastic fields are applied. Apparently, the conditions favorable for the ignition occurs in both these cases with nearly the same probability. Much more pronounced are the differences in the maximum temperature at the later stage when the flame starts to grow and propagate. For instance, at τ = 1.1 it can be seen that in the case HIT-200 T max is 230 K higher. However, when the flames are well established the difference in T max vanishes and at τ = 1.25 T max in HIT-200 is only 40 K larger than in HAT-200. The temporal evolution of V f is presented in Figure 8c. It can be seen that V f starts growing immediately after the temperature rise. Initially, its increase in HIT-200 and HAT-200 as well as for the laminar initial conditions is almost the same. Substantial differences are seen only when T max values in all the cases become similar (τ ≈ 1.25). From this moment the flame expands significantly faster in the case HIT-200 and at the end of the simulation time its volume is approximately 20% larger than in HAT-200. This means that the significant differences in the initial flow conditions have a larger impact on the flame propagation phase than on its initialization and very early development. Figure 9 shows the temporal evolution of the maximum temperature for all the cases from Table 1. The results corresponding to HIT and HAT are normalized such that the values T max are divided by T max,lam obtained for the laminar initial conditions at the end of the simulation time (τ = 1.45). It enables the comparison of the shapes of profiles and allows us to better quantify the impact of Re λ and the initial regimes. It can be seen that the auto-ignition always occurs faster for the turbulent initial conditions. To some extent, this seems to contradict with the findings of [19] who observed that for τ t /τ ign,l < 1 (the cases HIT-050, HIT-100, see τ t values in Table 1) the auto-ignition in turbulent flow was slightly retarded. There, however, the initial conditions characterized a strong turbulence that could have induced large values of the scalar dissipation rate delaying the ignition. Among the HIT cases T max grows significantly faster during the entire simulation time when the Reynolds number is large (Re λ = 150 and Re λ = 200). In the HAT conditions T max is much less dependent on Re λ and the turbulence level has a negligible impact on the ignition delay. This is because in the HAT cases, the turbulent time scale is nearly constant (τ t /τ ign,l = 5.8 − 6.1).
(a) (b)  Figure 9 additionally shows the solutions obtained for the lowest Reynolds number on the 70% denser mesh. Analysis of these results aims to complement the discussion on the very low level of the sub-grid kinetic energy presented in Section 2.5. It confirms the accuracy of the performed simulations. It can be seen that the use of significantly more mesh points virtually does not affect the obtained profiles. The difference in predicting the auto-ignition time is less than 0.5%τ ign,l , whereas the difference in respect to T max during the flame development phase does not exceed 1% T max,lam . Taking into account that these results were obtained for the cases at Re λ = 50 with very tiny initial flow structures it can be deduced that the basic mesh is appropriate also for the cases with larger vortices (Re λ = 100, 150, 200).
The normalised flame volume presented in Figure 10a (V f /V f,lam ) is found to be analogously affected by the initial conditions. The differences between the particular solutions are substantial and grow with time. In the HIT conditions the mixing process turns out to be much more intense than in the HAT regimes. This enhances the flame propagation across the mixing layer, and in effect, at the end of simulations the flame volumes in the isotropic regimes are significantly larger. Even for the smallest Re λ (HIT-050) V f is larger than in any HAT case. In the case HIT-50 V f is 25% larger than V f,lam , whereas in the case HAT-050 V f and V f,lam are comparable. Figure 10b presents the instantaneous variation of the flame volume dV f /dτ. It shows that the anisotropy causes a faster flame growth at an early stage of the flame kernel propagation (τ = 1.04-1.1) and the steepest change of the profile is found for the case HAT-050. In the later time (τ = 1.12-1.22) dV f /dτ is almost constant for all HAT cases that means that in this period of time the flame volumes grow in the time linearly. It can be seen that the values of dV f /dτ for HAT-100, HAT-150 and HAT-200 are practically the same and visibly larger than for HAT-050. This suggests that for the HAT conditions there is some threshold value of Re λ (or Re λ i ) above which the flame growth is the fastest and only weakly dependent on Re λ . The situation is significantly different in the HIT regimes.
The differences in dV f /dτ in particular cases are seen over the entire simulation time. Moreover, the profiles have a significantly different character with a well-marked maximum at τ ≈ 1.25. This corresponds to the time moment when the maximum temperatures in the domain for all HIT cases reach virtually the same level (see Figure 9a) and when the values of the flame volumes start to visibly diverge (see Figure 10a). The largest values of dV f /dτ are in the cases HIT-100 and HIT-150. Unlike as in the HAT conditions, this could suggest that when the initial turbulent field is isotropic there is an optimal moderate value of Re λ at which the flame expansion is the most dynamic.  Figure 11 shows the contours of the mixture fraction ξ and isolines of the selected values of the heat release in the central 'y-x' cross-section plane. The subfigures on the left half of Figure 11 correspond to the HIT conditions and on the right half to the HAT regimes. The presented results correspond to two time moments: (i) τ = 1.1 (upper row), when for the HAT cases dV f /dτ reaches the constant level; (ii) τ = 1.23 (bottom row), when for the HIT cases dV f /dτ are maximal. It can be seen that in the HIT conditions the distribution of the mixture fraction is much more wrinkled compared to the HAT regimes, especially in the cases HIT-050 and HIT-100. The heat release isolines are locally torn and the flame sources are dispersed. In effect, in these conditions the flame volume initially (τ = 1.1-1.23) grows slower than in the HAT regimes (Cf. Figure 10). The situation changes at the later time when the wrinkling translates to a larger surface area where the flame can appear (ξ MR ) and along which it can propagate. The increase of Re λ in the HIT conditions smoothes the mixture fraction contours, i.e., the flow structures become larger and local small wrinkles vanish. In the HAT conditions the change of Re λ causes the opposite behavior. In the case HAT-050 the values of Re λ 1,2 ('x' and 'y' direction) are small (see Table 1) and the mixing zone between the fuel and oxidiser is almost flat. The heat release isolines are continuous and extend over almost the entire 'z-x' plane. Unlike as in the HIT conditions, the increase of Re λ makes the mixture fraction distribution more complex. This is the effect of mutual interactions between the flow scales along with the particular directions. The large initial value of u 3,RMS corresponding to Re λ 3 induces velocity fluctuations both in the 'x' and 'y' direction leading to pronounced wrinkling of the flame surface.

Conclusions
The large eddy simulation method combined with Eulerian Stochastic Fields was applied to study the auto-ignition mechanism in the stationary mixing layer with the initial isotropic (HIT) and anisotropic (HAT) velocity field. The dynamics of the turbulent flames was compared to the laminar case including the entire combustion process starting from the auto-ignition, through the propagation of the flame kernels, up to the time moment when the flame developed along with the mixing layer. The ignition delay was found to be shorter for all the turbulent cases and the minor impact of the initial Reynolds number (Re λ ) was observed only for the HIT regimes. Contrary to the ignition times, the evolution of the flame temperature, the flame size and its temporal growth were found to be significantly dependent on Re λ , both in the HIT and HAT conditions. For its higher values (Re λ > 50) the flame volumes turned out to be at least 25% larger than in the laminar conditions. At the same time, the differences between the flames in HIT and HAT conditions were observed at the level of 22-34%, with the larger impact of Re λ in the HIT cases. It was found that the anisotropy of the initial velocity field seeds up the initial flame kernel growth but when the flame is developed it spreads over the domain slower than in the isotropic turbulent field. Analysis of the time evolution of flame volume variation has shown that in the HAT regimes its profiles characterize a plateau over some period of time, whereas in the HIT conditions distinct maxima occur. The flame surface waviness in the HIT regimes was found to be definitively more complex in the direction in which the velocity fluctuations and turbulent length scale were larger. It was observed that the increase of Re λ in the HIT conditions smoothers the flame surface, whereas in the HAT cases it makes the flame surface more wavy.
The performed research showed that depending on the type of initial conditions the flame dynamics can be significantly different and can be more or less dependent on Re λ . In the HAT conditions the Re λ values were set by changing the velocity fluctuations and turbulent length scales in each direction. In the future investigations, it would be advisable to consider the cases in which one of these parameters is kept constant while the other changes such to have the resulting Re λ the same. Such an analysis would give separate information to what extent the flame evolution depends on the turbulence intensity/length scale in the particular directions (perpendicular or parallel to the flame surface).

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.