Generalised Eddy Dissipation Concept for MILD combustion regime at low local Reynolds and Damköhler numbers. Part 1: Model framework development

Moderate or Intense Low Oxygen Dilution (MILD) combustion is a fuel-flexible combustion technology featuring high efficiency and low pollutant emissions. Fundamental studies reveal that turbulence-chemistry interactions are extremely complex in MILD conditions and reactor-type approaches seem to be the adequate modelling choice. In this work we develop a generalised Eddy Dissipation Concept (EDC) adapted to MILD combustion regime accounting for finite rate chemistry. We examine two recent modifications of the standard EDC and present a generalised model. It is based on functional expressions where the model parameters are adjusted to the local conditions in terms of Reynolds and Damköhler numbers, contrary to the usually proposed ad hoc tuning of the global EDC constants. Numerical results reveal that previously presented corrections are indeed suitable for specific conditions; their appropriate combination, guided by physical premises and a scrutiny of computation results, leads to a reformulation of the EDC framework. The study consists of two parts: the model development is described here (Part 1); in a companion paper (Part 2), we present a thorough validation process performed against twelve flames issuing from the jet-in-hot-coflow burners from Delft and Adelaide, representing a wide range of operating conditions. The new, generalised model can serve as a plug and play engineering tool without complex pre- or post-processing treatment.


Introduction
New combustion concepts are emerging or getting more popular, and are intensively investigated as a result of new trends in technology development. One of them is Moderate or Intense Low-oxygen Dilution (MILD) combustion, which is a recent and promising technique for increasing thermal efficiency and decreasing pollutant emissions in combustion systems. The technique is also called Flameless Oxidation (FLOX) [1], Highly Preheated Air Combustion (HPCA) [2], High Temperature Air Combustion (HiTAC) or Colorless Distributed Combustion (CDC) [3]. Earlier it was known as Excess Enthalpy Combustion (EEC) [4]. Although some guidelines are provided there on how to define each combustion regime, this is in fact an ambiguous problem. Usually the definitions are based on predefined temperatures of inlet reactants and the temperature increase during combustion, as presented by Cavaliere and de Joannon [5], and extended to configuration of hot diluted fuel [6] by Kwiatkowski and Mastorakos [7]. The general requirements for MILD combustion are that the inlet temperature of reactants is higher than the auto-ignition temperature of the mixture and that the temperature increase due to combustion is limited [5]. Evans et al. [8] presented a new MILD combustion definition based on an equivalent activation energy and classified as autoignitive some of the flames previously regarded as MILD. Recently, Luan et al. [9] presented an interesting study identifying MILD regime in a non-adiabatic well-stirred reactor. They composed a combustion map using initial inlet temperature and oxygen concentration. Regardless of the exact definition, this mode of combustion is characterised by hardly visible flame, inherent flame stabilization, slow reaction rates, nearly-uniform temperature fields and radiation fluxes, which is desired in some industrial processes. Recirculating hot gases, due to the heat recovery, contribute to the higher efficiency. Also, reduced oxygen level and low temperature rise across combustion zones lead to the emission reduction [10]. Fundamental aspects of MILD combustion of different types of fuels with a focus on industrial applications were presented by Weber et al. [11] and some issues of its mathematical modelling were discussed by Mancini et al. [12]. Application of MILD regime to gas turbine combustors was also recently https://doi.org/10.1016/j.fuel.2020.117743 Received 8 July 2019; Received in revised form 1 February 2020; Accepted 28 March 2020 reported by Khidr et al. [13]. A comprehensive review on flameless combustion and its potential for gas turbine applications, discussing the regime definition, basic experimental investigations, computational modelling aspects, design of the combustors including future challenges has been recently published by Perpignan et al. [14].
As mentioned earlier, MILD combustion, FLOX, HiTAC or CDC are all similar concepts which may differ by the implementation, the application field or the group of researchers or engineers dealing with it. Driscoll et al. [15] mentioned that all of those regimes can be associated with distributed reactions. Swaminathan [10] pointed out that other modern combustion concepts for internal combustion engines such as Homogeneous Charge Compression Ignition (HCCI), Reactivity Controlled Compression Ignition (RCCI), and Gasoline Compression Ignition (GCI) may also share similar features. Questions of fundamental nature can be answered by analysing representative DNS cases and experimental configurations reviewed in [14,15,10]. According to Driscoll et al. [15] two types of distributed combustion burners are usually identified: premixed and non-premixed design. The first one is associated with fuel mixed with hot and diluted oxidizer injected as a premixed jet. This is realized with high inflow velocity to avoid flame anchoring and ensure auto-ignition as presented in [1]. In case of non-premixed design, usually a fuel jet is accompanied by a hot and diluted coflow and auto-ignition occurs in the shear layer between the two streams. Typical example of such configurations are the Cabra flame [16] and Jet-in-Hot-Coflow (JHC) burners from Adelaide [17] and Delft [18]. They emulate MILD conditions and will be considered further in this study as validation cases.
In the MILD conditions the chemistry and mixing time scales are comparable and, thus, the Damköhler number is usually low and often approaches unity. Galletti et al. [19] performed an analysis of turbulence and chemical timescales and confirmed that in MILD combustion the high-temperature region is characterized by Damköhler numbers lower than those observed in flame mode. Li et al. [20] pointed out that the Damköhler numbers for the MILD combustion are in the range of 0.01-5.35, concluding that this regime is controlled by both flow and chemical time scales. Recently, Zhang et al. [21] numerically investigated the effects of different injection parameters on Damköhler number. They used the EDC combustion model and pointed out that reduction of the Da value is crucial in establishing the non-premixed MILD combustion regime. Tu et al. [22] investigated MILD combustion in atmospheres of different gases; to gain quantitative insight into the reaction regimes, the Damköhler and turbulent Reynolds numbers were used to compose diagrams. Most of the previous works applied simplified methods for chemical time scale estimation. However, there are also more advanced approaches to get this quantity as shown by, e.g. Isaac et al. [23].
Recently, Swaminathan [10] reported on physical insights on MILD combustion regime from Direct Numerical Simulations (DNS). It was pointed out that fundamental aspects of this regime are still not wellunderstood. The past DNS studies have been surveyed to provide a broader perspective on MILD combustion physics. The crucial aspects which were considered concern the inception mechanism for MILD combustion, its main combustion mode, typical morphological and topological features of reaction zones and the simplest ways how to model it. Swaminathan [10] stated that the classical S-curve cannot be used to describe MILD combustion inception and a new theory based on at least two chemical timescales has to be developed. It is worth to point out the DNS studies by Minamoto et al. [24][25][26]. They investigated MILD combustion reaction zones, their structure and features. One of their observations was that the reaction zones in MILD conditions strongly interact with each other. This behaviour leads to local thickening of the reaction zone, which gives appearance of distributed combustion despite the presence of local thin reacting structures. In distributed combustion, the length scales are larger than the smallest turbulence scales, indicating non-flamelet regime. Recently, Driscoll et al. [15] discussed conditions when distributed reactions occur. Sufficient level of preheat and dilution was confirmed as a necessary requirement. They pointed out that largescale distributed combustion could not be achieved only with the increased level of turbulence without preheat and dilution as shown by Skiba et al. [27]. Additionally, reactants slowly convert into products over quite a long distance of the order of a few centimeters. High temperature of reactants and the presence of a large amount of inert gas lead to much lower increase in temperature than in traditional combustion. As a result, the temperature gradients are relatively small and thus the diffusion term becomes negligible compared to the reaction source term in the energy conservation equation [15]. Therefore, the convection-reaction balance becomes more relevant than diffusion-reaction balance typical of classical flames well described by flamelet relations.
Such observations shed some light on how the reaction zone in MILD combustion can be modelled. In particular, it was mentioned that the reactor type modelling is adequate for MILD combustion [10,26]. They also suggested that there is room for improvement of models such as the Eddy Dissipation Concept (EDC), by accounting for changes in the local reaction zone volumes, considering the effect of dilution on characteristic scales [25].
In this work, we present a comprehensive development toward an easy-to-use and simple, yet accurate, turbulent combustion model incorporating finite rate kinetics at affordable computational cost, based on the Eddy Dissipation Concept. The commonly used EDC model is revisited and its applicability is extended to low Reynolds and low Damköhler numbers corresponding to the conditions often met in MILD combustion regime. The framework of the proposed model is based on recent modifications [28,29] whose relevance was assessed on a wide range of operating conditions represented by the flames issuing from the JHC burners from Delft and Adelaide.

Numerical modelling
Numerical tools can strongly support the development of new technologies covering the full chain beginning from fundamental research on the combustion physics up to design of combustors. The use of simulations in the process of product design, manufacturing and maintenance is now being popularised, and widely implemented within the idea of the fourth industrial revolution. Recently, Raman and Hassanaly [30] reported on emerging trends in numerical simulations of combustion systems.
The performance of any combustion model strongly relies on the accuracy of the associated chemical kinetic mechanism. It is well known that using appropriate chemical mechanism is crucial for emission predictions. Perpignan et al. [14] pointed out that, in MILD combustion, the reaction pathways and the interaction with transport phenomena are different than in classical canonical premixed and diffusion flames. Detailed kinetic studies devoted to diluted combustion processes are thus needed to gain insight in this matter [5]. Thus, there is ongoing research towards optimizing chemical schemes for MILD conditions [31,32]. So far, capabilities of existing chemical kinetic models to simulate flameless combustion have been assessed in several previous works, e.g. [33][34][35].
The closure for the turbulence-chemistry interaction (TCI) is one of the key problems in turbulent combustion modelling. Veynante and Vervisch [36] classified most of the existing TCI models into three groups, based on the one-point statistics approach, geometrical analysis, and turbulent mixing description. Some of these models are dedicated only to premixed or only to non-premixed combustion, whereas some can be applied to both regimes. They are characterised by various level of complexity and accuracy. Therefore, there is a wide choice of TCI models for the standard turbulent premixed or non-premixed combustion. This is not that obvious for the flames with preheated reactants such as Jet-in-Hot-Coflow (JHC) flames. These are often lifted flames exhibiting different features with respect to standard diffusion flames, and can fall into a separate class of high temperature combustion. In MILD combustion, the mixing processes, together with finite rate chemistry, are more important than in conventional diffusion flames, which should be taken into account when selecting the combustion model. The Damköhler number is usually low and often approaches unity in the MILD conditions as the chemistry and mixing time scales are comparable. In such case the numerical modelling of TCI is a challenging task, since some fundamental questions are still open as mentioned in the Introduction.
It is essential to develop robust and accurate numerical modelling approaches for industrial applications, accounting for finite rate chemistry. Due to the straightforward inclusion of the chemical source term the transported Probability Density Function (PDF) method is worth mentioning. However, it was reported to be sensitive to the level of velocity fluctuations and very computationally intensive [29,37,38], at least for the Reynolds Averaged Navier-Stokes (RANS) simulations. Recently, Chen et al. [39] proposed a simple model based on a Perfectly Stirred Reactor (PSR) using tabulated chemistry in conjuction with prescribed PDF and compared it with the transported PDF results from Lee et al. [40]. They reported comparable or even better results for major species and temperature with lower computational cost. De and Dongre [41] assessed the performances of the EDC and PDF based models, including transported PDF method in two different approaches. The Conditional Moment Closure (CMC) method was also investigated for MILD combustion modelling, e.g. by Kim et al. [42], and for autoignitive hydrogen jet flames issuing into a hot ambient co-flow, e.g. by Tyliszczak [43]. A similar method was adopted by Labahn et al. [44,45], who employed the socalled Conditional Source-term Estimation (CSE) accounting for two mixture fractions, both for RANS [44] and LES [45], tested on the Delft-Jet-in-Hot-Coflow (DJHC) flames [18]. Both approaches, transported PDF and CMC, are regarded as advanced and accurate; however, they are characterised by a considerable computational cost.
Simplified methods are therefore needed. Among them the Flamelet Generated Manifold [46,47] approach and mixing reactor models [39,48] have been intensively developed in recent years. However, as pointed out in the Introduction, in MILD conditions the smallest turbulence scales strongly affect the reaction zones and isolated laminar flame structures are difficult to identify. Distributed reaction zone can be described by the convection-reaction balance rather than diffusion-reaction as in classical combustion [15]. In the RANS framework, some studies [33,49,50] showed the limitations of the standard flamelet approach. However, the use of the extended flamelet-like approach both in RANS and LES [46,47,51,52] has shown good performance. Ihme et al. [51] applied a three-stream flamelet/progress variable (FPV) formulation in Large Eddy Simulations (LES) and obtained satisfactory results for the Adelaide JHC. Locci et al. [53] accounted for ternary mixing and proposed a new LES model based on diluted homogeneous reactors. Colin and Michel [52] presented a two-dimensional tabulated flamelet combustion model for furnace applications using RANS, and applied it to the flameless burner of Verissimo et al. [54]. Very recently, Chitgarha and Mardani [55] used RANS to assess steady and unsteady flamelet models for MILD combustion, and compared them to EDC. They reported better accuracy with the EDC than with flamelet models, yet pointed out a higher computational cost of the former approach. It has been shown that flamelet like models can handle these regimes but they need to include several progress variables to account for the dilution etc., with a consequent increase in pre-processing and storage needs. The discussion here does not aim to assess whether flamelet or reactor models are better suited for MILD conditions. Both methods have certain advantages and specific challenges to be addressed when employed in MILD conditions.
In the present work, we consider the reactor-type modelling approach focusing on EDC, initially introduced by Magnussen and co-workers [56][57][58]. It has been extensively investigated in MILD combustion studies in recent years [29,37,38,33,59,60] and selected modifications were summarised and commented by Ertesvåg [61]. The concept of treating the reacting structures as perfectly stirred reactor is well suited for the specific MILD combustion reaction zone. The model incorporates detailed, finiterate chemistry to provide explicit expression for the turbulent combustion closure in the balance equations using the mean flow properties. It is relatively easy in implementation and very intuitive for the user. Moreover, the EDC has already been a popular choice for TCI closure in the engineering applications. It has been successfully used in a wide range of systems implementing MILD conditions, such as prototype gas turbines employing flameless oxidation [62], pulverized coal combustion [63,64], forward flow furnaces of refinery-off gas [65], hybrid solar receiver combustors [66], laboratory scale [67,68] or semi-industrial flameless furnaces [69] and other devices [70][71][72]. Recent studies on the effect of various diluents on the MILD condition characteristics both in JHC burner configuration [73] and in laboratory-scale furnace [74] also employed EDC approach. The EDC does not require defining a mixture fraction, which can be puzzling in multi-stream systems, and operates on initial and boundary conditions in physical space, which makes it straightforward to implement in the design process of novel combustion technologies.
Despite conceptual advantages of the model, EDC faces a thorny problem. Its commonly reported drawback is over-predicting the temperature values in MILD regime [28,29,37,38,59,60,72]. In [28,29] we have reviewed previous attempts to deal with this problem. A widely adopted remedy is to use a strongly modified set of constants for the fine structure mass fraction and residence time, which are key ingredients of EDC, as discussed in Section 3. Yet a case-dependent character of such procedure is not satisfactory for a general applicability of the model. Therefore, some approaches aiming at generalisation of the EDC were recently presented. Parente et al. [29] proposed functional expressions for the EDC constants dependent on dimensionless flow parameters (the Reynolds and Damköhler numbers) validated against Adelaide JHC flames. They took into account specific features of the MILD combustion mode and applied the proposed changes globally and locally. This approach was further improved by Evans et al. [75] who used more accurate approaches for the determination of the chemical time-scale. Another recent modification suggested by Aminian et al. [76] was based on the Partially Stirred Reactor (PaSR) submodel to account for finite-rate chemistry in the fine structures tested for AJHC flames. Lewandowski and Ertesvåg [28] analysed the original EDC for MILD combustion and, based on the DJHC data, proposed a variable reacting fraction to improve the model predictions. A recent EDC modification proposed by Farokhi and Birouk [77] employed a fractal modelling approach to reconsider the energy cascade model.
The model extension using locally modified EDC constants [29,75] provided considerable improvement in the predictions without ad-hoc tuning, extending the model applicability for low Damköhler number cases. However, as pointed out by Ertesvåg [61], this approach maintained the multi-step cascade model for relatively low turbulence Reynolds numbers. As will be shown in Section 4.3.2, this in fact leads to problems for flames at low turbulence conditions. We will show that a remedy for this issue is to account for a variable reacting fraction, which proved to work well for Delft JHC flames [28], although this modification alone was not sufficient to improve the predictions of the Adelaide JHC flames. The main purpose of the present work is to demonstrate, based on physical premises, that the methodology combining both modifications leads to a generally applicable EDC approach.
A thorough validation process was performed against 12 flames to prove the generality of the method for a wide range of operating conditions. For the sake of clarity, the validation results together with additional information about the simulation strategy are provided in the second part of the work [78]. This paper focuses on the methodology development presented in Section 4, supported by empirical observations made during JHC flames simulations. The experimental data come from both the Adelaide [17] and Delft [18] JHC. Calculations were run using a RANS approach in OpenFOAM. The detailed description of the studied flames, numerical setup and other aspects of the mathematical modelling are provided in the second part of this article where a comprehensive validation is performed. In brief, attention was paid to a few particular aspects. Firstly a turbulence model was selected, accounting for the round jet anomaly which is a common practice, yet for low Reynolds number jets a deviation from the known trend was reported. Multicomponent diffusion effects was included, especially in hydrogen fuelled flames. Two reduced chemical mechanisms DRM19 (19 species and 58 reactions) and KEE (17 species and 58 reactions) were employed. In addition, new expressions for inlet boundary conditions for turbulence kinetic energy and its dissipation rate for a fuel pipe inlet were implemented.

Standard EDC: detailed chemistry approach
The Eddy Dissipation Concept (EDC) of Magnussen assumes that the reactions occur where the reactants are mixed at molecular level and the turbulence energy dissipation takes place: in the so-called fine structures whose size is of the order of magnitude of the Kolmogorov scales [56,58]. These structures are not evenly distributed in space but only in some of its fraction. EDC is based on the cascade model of energy dissipation [58] from larger to smaller scales, so that relations between the scales can be described with statistically averaged (RANS) quantities. A control volume is conceptually divided into fine structures and the surroundings. The role of the cascade model is to provide information on the fine structures. These are characteristic scales that we cannot resolve but model with the use of quantities from the mean flow. It is postulated that the ratio of the mass of regions containing the fine structures and the total mass in the control volume can be expressed as [56,58] is the turbulence Reynolds number. The constants C D1 and C D2 are introduced in the eddy dissipation turbulence energy cascade model which relates the fine structures to the quantities resolved by the turbulence model. Their numerical values were set to . The mean residence time in fine structures is modelled [58] as: where m is the mass transfer rate between the fine structures and the surroundings, divided by the fine-structure mass [58]. The secondary constants C and C , which are function of C D1 and C D2 as shown in Eqs. 1,2, are introduced by some authors for convenience. The standard values are = C 2.13 and = C 0.408. The mean reaction rate is the reacted mass of species k per unit of time per unit volume of the entire fluid in the cell. Including the assumption that only a fraction ( 1) of the fine structures actually reacts, in the formulation of Magnussen [57] the mean reaction rate in the conservation equation for species k is expressed as where¯denotes the mean density, Y k is the mean mass fraction of species k and Y k is the corresponding mass fraction in the fine structures, which remains unknown. Considering that the fine structures can be treated as a homogeneous, adiabatic and isobaric reactor, the mass fractions Y k are obtained from the solution of ordinary differential equations describing such reactor. This can be either a perfectly stirred reactor (PSR) or a plug flow reactor (PFR). The influence of the reactor formulation and the inflow or initial conditions was investigated previously [28,35]. It was pointed out that the issue should be treated with care but is not likely to cause significant differences. However, recently, Bösenhofer et al. [79] presented the analysis of different fine structure treatments in the EDC for classical combustion studied in Sandia D flame configuration. Several tests of different EDC variants were performed and the influence of the use of PSR and PFR approaches was investigated. Under the conditions of high turbulence and fast chemistry they observed some differences in reaction rates of intermediate species and temperature deviations up to 10%. It was also reported that differences in the mean reaction rates obtained with PSR and PFR increased at large and low values. In the JHC flames the values of the latter quantities are relatively high, which can explain why the differences reported in [79] were not observed in our simulations.
It has been reported in several studies (e.g. [29,37,38,59,60,72]) that the standard EDC requires improvements in non-conventional regimes such as MILD combustion. This has led to the calibration of the models constants (C and C ), to deal with the problem described in Sec. 3.2. However, it quickly turned out that this is rather a case dependent remedy and thus more general approaches have been proposed (e.g. [28,29,[75][76][77]).In the following sections (Sec. 3.3 and Sec. 3.4) propositions presented by Lewandowski and Ertesvåg [28] and by Parente et al [29] are briefly recalled and discussed.

The model constants adjustments
The most common modification of the EDC for MILD combustion is the use of adjusted values of C and C . Those constants were originally defined from the constants C D1 and C D2 in the energy cascade model [58]. A fundamental reason to increase C is the fact that reaction rates are slower in MILD regime. As C scales the fine scale volume fraction, the larger (or smaller) value of the constant would lead to larger (or smaller) volume fraction of a computational cell occupied by reaction region. Considering the conservation of mass and energy, as a result of lower C , a smaller volume fraction would lead to more distributed and overall larger reaction zone in the computational domain. This argument is used in studies suggesting lower values of C in MILD combustion.
In the case of DJHC flame, De et al. [37] suggested to increase C to 3.0 or decrease C to 1.0. Later, different combinations of the modified constants were also tested on the AJHC flames [38,59,35,60]. Recently, Ertesvåg [61] presented an extensive survey of proposed modified values revealing uncertainty of the choice, which causes lack of generality of such approach. He also commented on the common EDC modifications and advocated preserving the standard set of constants or at least use caution when doing so. The major concern was related with the interconnection of the turbulence energy cascade model and the EDC constants with the turbulence models.

Variable reacting fraction
Recently, Lewandowski and Ertesvåg [28] revisited the original EDC for MILD combustion modelling revealing some important implementation details. Initially, Gran and Magnussen [56] provided functional expressions for the reacting fraction of the fine structures. It should be clarified that the expressions for were most often introduced in the EDC when fast chemistry approach was employed, e.g. Lilleberg et al. [80]. The factor was then supposed to emulate finiterate chemistry effects and one can find different versions of expressions for in the recent work of Ertesvåg [61]. Gran and Magnussen [56] did not state that = 1 when using detailed chemistry. Initially, they even suggested that different should be used for different species rather using one global value. Nevertheless, they compared EDC with detailed chemistry approach using = 1 and variable . Since the results in the tested cases were very similar they came up to the conclusion that it is convenient to set = 1. On the other hand, Shiehnejadhesar et al. [81] and Magel et al. [82] claimed that when using EDC with detailed or finite rate chemistry approach, setting = 1 indicates that the reacting fraction of the fine structures is controlled by the chemistry. Ertesvåg [61] pointed out that "the argument was that the solution of the chemical reactor with detailed chemistry would determine whether and to what extent the reaction occurred. An alternative interpretation of this could be that the quantity Y Y ( ) k o k was solved from the reactor, rather than Y Y ( ) k o k . That is, was not unity but became a part of the solution". Nevertheless, in most instances = 1 has been set when a detailed chemical mechanism is used.
Recently, Farokhi and Birouk [83] developed a hybrid EDC/Flamelet approach and investigated potential usage of non unity , employing definition of variable reacting fraction from Magnussen [84]. In [28] it has been shown for Delft JHC flames that = 1 may not be justified in the non-stoichiometric conditions with low turbulence Reynolds number and incomplete reactions. The reacting fraction of the fine structures, , was based on the global one-step irreversible reaction. Following the assumption presented by Gran and Magnussen [56], the reacting fraction may be expressed as multiplication of three factors = · · , 1 2 3 described in detail by those Authors, see also [80,28].

Locally modified EDC constants
One of the recent extensions of the EDC to MILD conditions is the correction proposed by Parente et al. [29]. They revised the energy cascade model with the assumption that the MILD combustion occurs in the so-called distributed reaction regime. In such conditions the flame structure is thickened and distributed due to the fact that the eddies with the size of Kolmogorov scale are able to enter and thicken the preheating zone and reaction region. The characteristic speed of reacting fine structures u was then assumed to be equal to the turbulent flame speed S T using the Damköhler expression The expression used for the turbulent flame speed is based on the second of Damköhlers two hypotheses. It assumes that the ratio between the turbulent and laminar flame speeds is proportional to the square root of the ratio between the turbulent and laminar diffusivities. This hypothesis is in principle applicable only to turbulence with all length scales much smaller than the flame thickness. Such a regime, however, is seldom obtained in practice because the integral length scale of practical turbulent flows is typically several times greater than the laminar flame thickness of the unburned mixture [85]. Hence, an intermediate regime is more realistic, where some turbulence length scales are smaller and some are larger than the flame thickness. Following this reasoning, the relationship we have employed for the turbulent flame speed should include an effect of increased flame surface area beside the increase of turbulent diffusivity.
It should be noted that the use of the laminar flame speed S L , which is a kinetic combustion quantity, is justified by the large degree of partial premixing occurring in MILD regime. With this assumption we can reconsider the turbulence energy cascade model [58] estimating the characteristic speed of the turbulent reacting fine structures with = u S . T Following [29], one can obtain a new expression for C D2 constant, and secondary parameter C can be expressed as a function of Da and Re T where C 0 is a proportionality coefficient equal to 0.5. The expression for C can be derived similarly and presented as . It should be however pointed out that C 0 and C 0 are not strictly set values since proportionality expressions were used during the above derivation. Thus further adjustments, as in the work of Parente et al. [29], are justified to some extent. However, the issue of how sensitive the model is to the changes in those new constants remains. Additionally, to avoid too large changes in the local parameters, the values were bounded between the standard EDC constants values and some upper and lower limit [29], so that C 0.408, 5.0 and C 0.5, 2.14 . The same proportionality coefficient C 0 has been used by Evans et al. [75] and the step by step derivation of both coefficients is shown in [86]. Further on this approach will be referred as "v2016" model version.

General remark and motivation
Now the method of coupling the two approaches presented in Sections 3.3 and 3.4 will be presented, leading to a new version of the model denoted as "generalised EDC". In [28] it was shown that the use of a variable reacting fraction approach could successfully simulate the DJHC flames, showing considerable improvement compared to the EDC approach using = 1 and the standard set of constants. Although the Delft flames operate as jets in hot co-flows to emulate MILD conditions, there has been some recent discussion whether they can be fully classified as MILD cases [8]. However, their importance on investigating a stabilization mechanism fundamentally different from that of conventional diffusion flames [14] is not questioned. A further validation study was performed on the AJHC flames [78] It was shown that the effect of employing a variable is relatively small compared to = 1, supporting the conclusion of Gran and Magnussen [56], indicating that the variable alone does not represent a general solution for MILD combustion modelling with EDC. On the other hand, the AJHC cases have been successfully simulated with the method of locally modified EDC constants by Parente et al. [29]. However, the same approach shows problems when applied to DJHC flames. In [78] it is showed that the flame is incorrectly predicted suggesting too high decrease in reaction rates due to the strongly modified constants.
Two things come out of those observations. Firstly, what are the reasons for the respective models failure or success in the considered cases? To answer this question both the Adelaide and Delft JHC flames have to be examined considering their differences and similarities. Secondly, the complementary effect of the two modifications suggests that their combination may lead to more general solution for TCI modelling of MILD combustion. In order to address these issues, a closer look to the role of the local turbulence Reynolds Re and Damköhler Da numbers is taken in the next Section.

The Re Da reaction zone map
Let us now focus on the two local distribution of the Da and Re numbers. In Fig. 1, the fields of Da and Re are presented together with the OH mass fraction map as an indicator of the reaction zone obtained with the standard EDC model using = 1, henceforth denoted as BFM2005 [57]. Then, conditioned on the OH mass fraction being greater than 2.8·10 4 , the values of Re and Da are plotted against each other for each computational cell and shown in Fig. 2 for the Delft and Adelaide burners. Such contours represent Re and Da values in the reaction zone, thus from the model point of view this can be understood as working conditions. It should be pointed out that Da in Fig. 2 has the same range in both plots, yet Re extends up to 100 for the Delft flames and up to 200 for the Adelaide flames. The first conclusion can be made that the Delft flames operate under lower Reynolds numbers than the Adelaide flames. This fact aligns with the observation of Oldenhof et al. [87] who reported that the chemical reactions in DJHC flames occur on the outer side of the interface that separates the turbulent and non-turbulent part of the flow. They explain that on average those flames are in a high-shear, high-turbulence part of flow, but do not experience this turbulence. The second clear observation is that for the lower jet Reynolds number, the maximum Damköhler numbers were higher and, especially clear for the AJHC flames, those values were shifted towards lower Re . This is reasonable, since for lower levels of turbulence we have longer Kolmogorov time scales, which leads to higher Damköhler numbers.
When it comes to the flames under the same flow conditions but different levels of dilution, as in the example of AJHC HM1-3 flames, it is observed that the Re vs Da maps are very similar but flames with higher amount of oxygen in the coflow have slightly higher Damköhler numbers in the reaction zone. This can be explained by the fact that with increasing oxygen content in the coflow, the maximum temperature increases as well, which leads to shorter (faster) chemical times and thus, higher Damköhler numbers. These quantitative observations assure that considerations based on this kind of maps are reasonable and can serve as a useful source of information. However, it should be pointed out that the maps presented in Fig. 2 are drawn up from the simulation results, thus their shapes will depend on the simulation approach.

Determination of the low Da limit
Locally modified EDC constants expressed with Eqs. (5) and (6) as functions of the two variables (Da and Re ) are presented as surface plots in Fig. 3. Those expressions are effectively implemented only in a certain range of C and C . As mentioned in Section 3.4, the value of C is bounded between 0.408 and 5.0, whereas C changes between 0.5 and 2.14. The limiting lines from Fig. 3 were projected into the Da and Re plot and are presented in Fig. 4. Above the green and blue lines, the model works using the standard EDC constants, in between those lines only one of the parameters is changed, whereas the other has the standard value. There is not much to discuss on these upper limits, since they provide a smooth and consistent transition to the standard EDC for higher Da and Re which is desired. However, setting the lower limits is not that obvious, since there is no clear physical premise. We have taken the values of = C 5.0 and = C 0.5 as suggested by Parente et al. [29], which comes from empirical observations. Moreover, a further increase in C and C would likely lead to extinction, and another good reason for setting such values as the limiting ones is that the two limits coincide on the Da and Re plot. A persisting problem is that for very low Da the EDC parameters reach extreme values. In the original work by Parente et al. [29], C was set to its standard value for < Da 0.01, to ensure ignition. The authors reported that the results were insensitive to a variation of 50% around this value. Yet, from our experience it follows that for different cases with varying Re , the change in the threshold of Da can improve the results. Our analysis ashowed that decresing the threshold of Da would lead to flame extinction at some point, which was dependent on the level of turbulence in the flame. For example, using a constant Da value of 0.04 led to the extinction of the AJHC HM1 flame at = Re 5000, whereas for the case AJHC HM1 at = Re 10000, the value of 0.03 was found to be the minimum to enable ignition. For other cases, HM2 and HM3, the limit could be decreased even more, leading to some small improvements, yet no universal constant threshold Da could be estimated. Therefore, using a Da threshold gently dependent on Re seemed to be more appropriate, leading to accurate results and stable solution in a wide range of conditions. The lower Da limit was set to be consistent with the upper limit of C (the red line in Fig. 4) as an interpolated, exponential function of Re . Below this threshold both C and C remain at their standard values.

Determination of the low Re limit
As mentioned in Section 4.2 and shown in Fig. 2, the reaction zones of DJHC flames experience much lower turbulence than in case of AJHC flames. Thus, poorer performance of locally modified EDC parameters may be somehow related to the low Re conditions and another limit is likely required. It is now important to determine where to assign such a limit and how to justify it.
In order to pursue qualitative investigations, numerically predicted radial temperature profiles were examined for each flame considering the local flow conditions. The maximum temperature values were identified and plotted in a Re vs Da map for the two investigated EDC modifications. Each maximal temperature from the simulations was compared to the experimental data and a relative error was assigned. The results from the variable approach and locally modified EDC constants have been compared for representative JHC flames. Fig. 5 shows an example of such comparison for the AJHC HM1 = Re 10000 case. The maximum temperatures for the radial distribution at four axial locations (30 mm, 60 mm, 120 mm and 200 mm) are localised and denoted with symbols on the respective Re vs Da map. Arrows point to the respective peaks in the temperature distribution plots. Fig. 5 allows to determine which values of Re and Da correspond to the maximum temperature results predicted by the two EDC variants. In Fig. 5 it is apparent that for the first two locations "variable " gave better agreement with the  experiment, whereas for the other two "v2016" performed better. It can be stated that Da is very different for both methods, yet transition in accuracy between the two corrections takes place in physical space between 60 mm and 120 mm downstream the AJHC HM1 flame. In the Re vs Da map this corresponds to Re between 22 and 38, and between 25 and 60 for variable and "v2016", respectively. Considering the three Delft and five Adelaide flames, in total 25 points with different local conditions spanning up = Re 100 were selected and compared. It was found that for < Re 23 all the "v2016" results were always outperformed by "variable " and for > Re 31 "v2016" was superior. The reason for poor performance of variable approach for > Re 31 is due to the fact that the value of was close to unity resulting in a modest effect, similar to the standard EDC predictions. However, it is not yet clear why "v2016" struggles at low Reynolds number conditions present in all Delft flames and in the upstream part of Adelaide flames. Keeping in mind that the deterioration of accuracy occurs at Re between 23 and 31, let us now consider a scenario of low turbulence, where the number of levels in the cascade is small. The cascade model is constructed as explained by Ertesvåg and Magnussen [58], but the condition that = q 4 3 is abrogated at low turbulence. If we assume that there is only one level of the cascade, we can write the energy transfer to the first level from the mean flow as and, therefore In such situation the turbulence Reynolds number is equal to Below that number there is no cascade at all. Therefore, any sophisticated corrections to the turbulence energy cascade model (and "v2016" is such correction) may lose justification and applicability. One may point out that Re was calculated using the standard C D1 and C D2 constants values. Employing a strongly modified set of constants, the value of Re when the cascade has just one level can go down much below unity. The concept of turbulence energy cascade at such low turbulence would be difficult to hold. Interestingly, the above limit falls into the empirically estimated Re range between 23 and 31 for which the model with locally modified EDC parameters starts to considerably lose accuracy.

The generalised model framework
The process of model development and differences in each variant can be visualized in Re vs Da map presented in Fig. 6. As shown in the first diagram, the standard EDC employs the assumption of = 1 and uses constant values of C and C regardless of the local conditions. The second diagram represents the "v2016" model variant behaviour as presented in [29]. The third diagram includes modifications which led to the generalized model. Taking into account the above guidelines, the constant value of = Re 28 was used as a threshold below which the EDC parameters C and C are switched back to their standard constant values. The use of variable reacting fraction ensures proper model behaviour at low Re as previously shown in [28]. For higher Re the value of is close to unity and does not interfere with the effect of locally modified EDC constants. Additionally, the low Da limit was In summary, in this model version as in "v2016" the expressions for locally modified EDC constants ensure a smooth transition to the standard EDC constants. Variable reacting fraction is used instead of = 1 assumption and stricter bounds on modified constants are set. This prevents the model from nonphysical behaviour at extremely low Re and Da . At this point the new model framework has been set. Yet, since Re and Da are the new input parameters to locally modify the model constants, the impact of their definition shall be assessed in the next section. It concerns especially Da , where the estimation of chemical time scale is needed. An extensive validation of the generalized model in the Jet-in-Hot-Coflow environment was performed and will be presented in the second part of this work [78].

Sensitivity to the chemical time scale estimation
The locally modified EDC parameters are functional expressions that depend on turbulence Reynolds number and the local Damköhler number. The model sensitivity to those local variables must be investigated. The calculation of the turbulence Reynolds number is straightforward and does not require so much attention as the local Damköhler number. It is due to the various ways of the chemical time scale estimation. Although the simulations are performed with detailed chemical mechanism, the evaluation of a characteristic local Damköhler number Da can be made with some simplifications.
Originally, Parente et al. [29] proposed the use of one-step chemistry to estimate the controlling chemical time scale c . They took a set of the reaction rate parameters for methane-air mixture given by  Westbrook and Dryer [88], which in Table 1 is denoted as set 1. We have defined chemical time scale using all five variants given in Table 1. Since the sets 1 and 2 represent first order reactions, the chemical time scale was estimated as in [29], i.e. the reverse of the rate coefficient (but not larger than 0.1) expressed with the Arrhenius equation. In case of sets 3-5, we took the definition from Turns [89], accounting for the fuel and oxidizer concentrations. Comparison of the numerical results for the flame AJHC HM1 at = Re 10000 obtained using those five chemical time scales is shown in Fig. 7. In the upstream part of the flame almost no differences can be observed, yet going downstream the flame set 5 led to higher temperatures, set 2 to lower ones and the other three sets gave approximately the same results.
It should be pointed out that even if we use global, one-step chemical mechanism to estimate chemical time scale for obtaining the local Damköhler number, the CFD computations still use a chemical mechanism consisting of 17 species and 58 reactions [90]. It is of interest to check what would be the effect of using the leading chemical time scale based on this mechanism. Evans et al. [75] recently proposed to take the maximum chemical time of only the five major species CH , H , O 4 2 2 , CO and CO 2 , showing that the same leading time-scale was obtained as when all species were considered. We have adopted this approach using the following expression: (for = … k 1, , 5) where two variants were considered to provide dY dt / k . In the first variant (KEE 5sp I), it was calculated as Y Y / k k . In the second variant (KEE 5sp II), it was evaluated from the formation rate taken from the edcSMOKE library, which handles chemistry in the edcSimpleSMOKE solver [91]. In Fig. 8 the numerical results using those two approaches are compared with chemical time scale obtained from the global chemical reaction with parameters from set 1. There is no difference between the results at the locations 30 mm and 60 mm (not shown), whereas at the location 120 mm the results obtained from global reaction deviate from the results based on the two approaches using five selected species from the KEE mechanism. At location 200 mm, significant differences are observed between all the three data series. The first variant using five species underestimates the temperature and CO 2 but provides excellent agreement with the OH peak value.  [88] given in Table 1.
The second variant highly overestimates all three scalars, whereas results with set 1 are placed in between.
Medwell et al. [92] reported that at some point downstream of the JHC flames surrounding air is entrained and is able to mix with the hot coflow stream. Then it is visibly evident that the flames are perceptibly different in structure. As a result, beyond that location the combustion regime can deviate from MILD conditions. Therefore, one should be cautious when interpreting experimental data downstream the flame. Especially, validation of numerical models should be restricted to the region not affected by the entrained air. Nevertheless, numerical analysis downstream the flame is still worth investigating to compare different models results and present their behaviour.
Recently, Evans et al. [75] claimed that in case of AJHC HM1 at = Re 10000 the local oxidant composition is controlled by the coflow stream for approximately 150 mm, whereas further downstream of the jet exit plane this role is taken over by the entrained ambient air. They also suggested that for the HM1 flame, air entrainment into the fuel jet was enhanced by reaction zone weakening, resulting in a sudden increase in temperature and the occurrence of the so-called flame brush that was also experimentally observed. Considering the above discussion, it appears that the results provided with the second variant using five species can correspond to those observations. However, the flame brush was absent in the experimental data we were provided with. It this context it is also worth to point out the recent work of Li et al. [48] who presented different methods of chemical time scale estimations applied to AJHC flames simulated with PaSR model. Nevertheless, it was is worth mentioning that estimation of the chemical time scale can be a critical point, since it can easily vary by an order of magnitude leading to correspondingly different Da values. One can regard this as a drawback of the model, yet we should remember that the role of c in this model is to estimate the value of Da . Without a clear guideline in this respect, we d we believe the most robust estimate of the chemical time scale estimated is obtained using a global, one step reaction with rate parameters from set 1, as originally suggested by Parente et al. [29]. This approach keeps implementation of chemical time scale simple and impacts only calculation of Da , while species and temperature results are still based on detailed chemical mechanism providing reasonable results.

Conclusions
This work presents the development of an accurate and efficient computational model for turbulence-chemistry interaction in MILD combustion using RANS. Having in mind the role of LES in the current state-of-art combustion modelling research, there is still room for cost-efficient RANS combustion models to be used in industrial design and development.
After in-depth analysis of the EDC approach, discussion and assessment of various model formulations, a generalised model was developed. The proposed modifications were dictated by the unique features of MILD combustion regime. The proposed new variant of the EDC is a specific combination of the revised original EDC [28] with the locally adjustable EDC constants approach [29]. To highlight the main findings of the present work, the following points can be made: • variable reacting fraction becomes influential only at very low Re . For stronger turbulence is close to unity; • the model constants are locally modified as functions of Re and Da in the ranges marked in Fig. 6; -C and C start to respectively increase and decrease from their standard values at the upper limit of Da as function of Re , -the lower limit of Da was set when the two curves for C and C overlap with each other at 5.0 and 0.5, respectively. Below that limit the standard constants values are used to ensure proper ignition. -for < Re 28 the turbulence energy cascade breaks and the expressions for locally modified EDC constants lose applicability, and standard constants values are retrieved. In this region the reacting fraction becomes lower than unity and is supposed to suppress the reaction rate overestimation; • expressions for locally modified EDC constants are dedicated for < Da 1, to ensure the fine structures size larger than the Kolmogorov length scale; • based on the above findings, a framework for the generalisation of the EDC model was presented extending its applicability for low Reynolds and Damköhler number conditions; • employing a more accurate chemical time scale estimation method is foreseen to further improve the model.
To summarise, because of the specific features of the MILD combustion, the generalised EDC model requires the use of variable reacting fraction and the ranges of the constants' modifications are now clearly specified using the functional expressions. The proposed approach adjusts the EDC parameters to the local flow conditions and does not add complexity from the user point of view, thus it can serve as a convenient "plug-and-play" engineering tool in the field of numerical combustion. The companion paper (Part 2, [78]) presents in detail the validation process against twelve flames issuing from the jet-in-hotcoflow burners from Delft and Adelaide.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.