Mathematical Modelling of Turbulent Combustion of Two-Phase Mixtures of Gas and Solid Particles with a Eulerian–Eulerian Approach: The Case of Hydrogen Combustion in the Presence of Graphite Particles

Modelling of Turbulent Combustion of Two-Phase Mixtures of Gas and Solid Particles with a Eulerian–Eulerian Approach: The Case of Hydrogen Combustion in the Presence of Graphite Particles. Abstract: The numerical modelling of turbulent combustion of H 2 –air mixtures with solid graphite particles is a challenging and key issue in many industrial problems including nuclear safety. This study presents a Eulerian–Eulerian model based on the resolution of the Navier–Stokes equations via large eddy simulation (LES) coupled with a system of ordinary differential equations (ODEs) of the detailed chemical kinetics to simulate the combustion of mixtures of gases and particles. The model was applied to predict the transient evolution of turbulent combustion sequences of mixtures of hydrogen, air and graphite particles under low concentration conditions. When applied to simulate lab-scale combustion experiments, the results showed a good agreement between experimental and numerical data using a detailed chemical kinetic model. Moreover, the model was able to predict some key experimental tendencies and revealed that the presence of a low concentration of graphite particles (~96 g/m 3 ) in the scenario inﬂuenced the hydrogen combustion dynamics for mixtures of 20% (in volume) of hydrogen in air. Under these conditions, pressure levels reached at the walls of the sphere were increased and the combustion time was shortened. The results also showed the viability of using this kind of a model for obtaining global combustion parameters such as wall pressure evolution with time.


Introduction
Combustion of gas and particles mixtures is an issue of major interest in many different fields, including industrial combustors [1][2][3], pollutant emissions [4][5][6], solid propellants [7][8][9] or accident prediction and mitigation [10,11]. In this last field, prediction of particle behaviour with and without combustion is a key topic in nuclear power plants [12,13] as well as in fusion reactors such as the International Thermonuclear Experimental Reactor (ITER) [14,15]. In this case, the presence of particles might influence the combustion dynamics during a potential accident [16,17]. Therefore, it is of outmost importance to properly predict the effects of this type of turbulent combustion sequences in presence of solid particles.
Mathematical modelling of turbulent combustion requires a proper description of two key aspects: chemistry and turbulence. Large eddy simulation (LES) is a mathematical Mathematics 2021, 9,2017 2 of 17 approach that provides a compromise between efficiency in terms of computational costs and detailed physical description of turbulence dynamics [18]. Several models are available for modelling turbulence at the subgrid scale (SGS) with LES. The Smagorinsky-Lilly [19] was the first one. It was developed for flows with homogeneous turbulence conditions. The wall-adapting local eddy model [20] is a common choice in the case of incompressible wall flows whereas the dynamic k sgs equation model has been shown to behave relatively well in compressible flow conditions [21][22][23]. Regarding chemistry, it is important to use detailed kinetic models in order to predict several mechanisms such as ignition or quenching [24]. This is critical not only in the case of hydrogen combustion [25], but also in the case of other mixtures such as hydrogen and carbon. However, in this case, the number of detailed chemistry models available is still limited. Saxena and Williams [26] proposed a chemical kinetic mechanism for the combustion of hydrogen and carbon monoxide with 13 species and 30 reactions. They showed good results in their testing. Gibeling and Buggein [27] also developed a simplified model for the oxidation mechanism of carbon monoxide in the presence of hydrogen and oxygen. It considered nine species and 12 reactions and was used in propellant applications with satisfactory results. Zhuo et al. [28] also used a chemical kinetic model with eight species and 12 reactions for modelling carbon monoxide oxidation in the presence of hydrogen for propellant applications. In the case of mixtures of gases and particles, Bournot et al. [29] simulated a reactive two-phase flow with aluminium particles for base bleed applications. In this case, the chemistry model was very simple and only considered three species for the gas phase and two species for the particle phase. The turbulent nature of the flow was modelled in a statistical way. The results permitted to identify the global flow regimes of the problem. As for carbon particles, Chelliah et al. [30,31] studied the influence of particle porosity on the combustion of graphite particles in the presence of hydrogen under quasi-steady burning conditions. They relied on the chemical kinetic models proposed by Bradley [32] for nonporous graphite particles and by Yetter et al. [33] for CO-H 2 O-O 2 . A total of five reactions were considered for the solid carbon phase and 28 for the gas phase. A total of 13 species were considered in the problem. They confirmed the suitability of this kind of kinetic models with a comparison against experimental data. As for the modelling of particle behaviour in a turbulent gas flow, results of direct numerical simulations (DNS) for canonical problems permitted to give credit to less demanding numerical models such as LES. The Eulerian treatment when solving the governing continuum equations for averaged quantities of both phases is still limited with LES and DNS. There are just a few articles that have been devoted to the implementation of the Eulerian two-fluid approach in the framework of DNS or LES [34,35]. Yeh and Lei [36] used LES to investigate the motion of particles in isotropic and homogeneous shear flows. The generated particle-statistics by neglecting the subgrid scale (SGS) effects on particles showed that LES can successfully predict second-order statistics of particle motion. Similar results were obtained by other researchers [37][38][39][40].
As previously said, in the case of a combustion sequence within a dust-laden atmosphere in an industrial environment, a mining environment [10] or an accident scenario at a fission nuclear power plant [11] or the ITER fusion reactor [17], the presence of graphite particles might influence a potential hydrogen deflagration [41]. The sequence might degenerate into a dust explosion which might increase the accident impact [17]. Therefore, it is of major importance to properly predict the effects of this type of combustion sequences [41]. Previous works have assessed the effectiveness of LES models to predict the dynamics of premixed turbulent combustion of H 2 -air mixtures [23]. In this work, we present a two-phase reacting model based on a Eulerian-Eulerian approach for both the gas and the solid phases. The model includes LES turbulence modelling and detailed kinetic schemes for the combustion chemistry of both the gas and the solid phases. These two approaches are not usually used in two-phase flow models that consider a Eulerian approach of both the gas and the solid phases. We explore the application of this model to a combustion two-phase flow problem with graphite particles in the presence of a premixed H 2 -O 2 -N 2 atmosphere. As commented above, this scenario is prototypical of several key industrial and safety problems. Specifically, we focus on the case of sequences in closed three-dimensional (3D) H 2 scenarios and investigate the effect of the presence of a low concentration of graphite dust particles in the combustion sequence. The comparison of combustion experiments from the literature with the numerical simulations permits to face the validation of the mathematical model proposed and evaluate its prediction capabilities. The article is structured as follows: firstly, the physical model and the numerical method used are presented. Later, the model is validated against ad-hoc experimental results obtained in a spherical bomb and the results obtained are discussed. Finally, some conclusions are drawn.

A Two-Phase Flow Model for Mixtures of Gases and Solid Particles
The two-phase model used in this work relies on the hypothesis of highly diluted mixtures of gas and solid particles and considers a Eulerian-Eulerian approximation of the mass, momentum and energy conservation equations for both the gas phase and the solid phase. The system of equations also includes those corresponding to the concentration of the species of each phase. The gas phase is considered to be an ideal gas initially composed of a mixture of hydrogen and air (i.e., oxygen, nitrogen and argon). The solid phase was initially considered to be a monodisperse distribution of graphite (i.e., carbon, C) particles of 35 micrometres in diameter (Sauter median diameter) and particle density of 2160 kg/m 3 . During the combustion process, these species may react and generate additional species that will be detailed in the Combustion Model section.
In this case, the system of conservation equations averaged with the FAVRE approach (filtered local volume average of the equations) [42] was as follows: where the gas mixture is formed by NGSP gases and the solid mixture is formed by NPSP solids. The subscript g indicates variables relating to the gas phase, while the subscript p refers to the solid phase; ρ g is the average gas density, → u m is the velocity of phase m (i.e., gas phase "g" or solid particle phase "p"), Y m,k represents the mass fraction of species k from phase m.
The hypothesis of a highly diluted mixture implies that void fraction (α) could be assumed to be near unity (α ≈ 1). Therefore, the concentration of particles is defined by σ = ρ g (1 − α). E m and H m indicate, respectively, the total internal energy and the total enthalpy of phase m (where "m" can be solid phase "p" or gas phase "g"). Mass fractions Mathematics 2021, 9, 2017 4 of 17 are defined so that if NGSP is the number of components of the gas mixture and NSSP is the number of solid species, the following relations are fulfilled: Note that = τ represents the stress tensor and includes the turbulent (subgrid) stress terms, → q is the heat flux vector, Y k is the mass fraction of species k, . ω k is the reaction rate of the k species and D k is the diffusion term of species k. Note that the expression of the species conservation equations includes the terms of the thickened flame model (TFM) used for modelling the turbulent combustion mechanism. These terms are explained in detail in the Combustion Model section of this work. Similarly, the energy conservation equation of the gas phase also includes the corresponding terms of the TFM in the transport terms (i.e., Q c is the heat released per unit of volume and time due to the chemical reactions, which is defined as follows: .
where ∆H f m,k is the formation heat of species k.
In the model, the pressure effect on the solid phase is negligible and, therefore, solid particles can be considered to be incompressible. The gas-particle interaction was taken into account through source terms in the mass, momentum and energy equations. → F d is the gas-particle drag force, . Q g is the interfacial heat transfer rate, Γ stands for the total mass exchange between phases, . ω m,k -for the species reaction rates. In the present approach, particle size was considered to be constant. This means that the model does not consider the change of particle diameter during the combustion process.
The equation of the state considered for the gas phase was as follows: where p denotes gas pressure, T-the gas temperature, R u -the universal constant, M kthe specific molar mass of the species k. Specific heats are temperature-dependent following the database by McBride et al. [43] and dynamic viscosity was considered to be temperaturedependent through Sutherland's formula [44]: Closure Equations for the Interphase Transport The momentum exchange between the gas phase and the solid phase was taken into account by considering the drag force acting on a particle: where C d is the drag coefficient which is a function of the particle Reynolds number Re p and A are the representative area of the particle [45]. Considering that the number of particles is σ/m p , where m p is the mass of a particle, and assuming spherical particles [46], the previous equation can be expressed as follows: The expression adopted in this work for the drag coefficient C d is the one proposed by Otterman and Levine [47] and used by Miura and Glass [48] in their work: where the Reynolds number for particles Re p = The rate of the heat transferred from the gas to a particle at its surface, . Q g , is as follows: where the Nusselt number (Nu) can be calculated as follows: where Re is the gas Reynolds number, Pr-Prandtl number. This equation is valid for Re ≤ 50,000 according to Crowe et al. [45]. It was originally proposed by Knudsen and Katz [49] and has been used by different authors [47,48].

Model of the Chemical Kinetics
In order to evaluate precisely the rates of the chemical reactions present in the problem, a system of 33 ordinary differential equations (ODE) is considered and numerically solved to calculate the concentration of the different species at each timestep. In this work, the detailed chemical kinetics model used by Chelliah [30,31] was considered to describe the combustion of H 2 and solid carbon (C) in the presence of air. It was based on the chemical kinetic models proposed by Bradley [32] for nonporous graphite particles and Yetter et al. [33] [30,31] and Yetter et al. [33] that were adopted in this work, it has been considered that atmospheric N 2 and Ar are nonreacting species that do not intervene in any of the reactions and do not undergo any oxidation process. This was done to maintain the integrity of the chemical models, without making any modifications that could alter their validity.
A total of five reactions were considered for the solid carbon phase (Table 1). Two different mechanisms were tested in this study. One considered semiglobal heterogeneous surface reactions for nonporous graphite particles (Table 1a). The other considered semiglobal heterogeneous surface reactions for porous graphite particles (Table 1b). The mechanisms assumed that the primary product was CO and that each rate w < s independent of the others. The first assumption is known to be reasonable for high-temperature oxidation of graphite (i.e., temperature at the particle surface over 800 K) [50]. Note that in Table 1a, the rate . s i is expressed in terms of k i = A i T αi exp(-E i /RT), the partial pressure P j is in Pa, T-in degrees Kelvin, whereas in Table 1b, the rate . s i is expressed in terms of . s i = W i c i B i T αi exp(-E i /RT) in kg/m 2 /s. In this formula, B i T αi is in m/s, T-in degrees Kelvin. The gas phase homogeneous reaction mechanism of CO oxidation in the presence of H 2 O considered here is the one proposed by Yetter et al. [33]. It consists of 12 species in 28 elementary reactions ( Table 2). The rate constants for this mechanism were validated for a wide range of temperatures, pressures and reactant concentrations using shock tubes and flow reactor measurements. Following Yetter et al. [33], for high-temperature oxidation of CO, the non-Arrhenius rates recommended for reaction steps 4 and 22 were implemented.  [30]). (b) Heterogenous rate constants for porous graphite oxidation (from Chelliah [30]). Here, For the species k, the total mass reaction rate . ω k was defined as the contribution of all reactions: .
where NR is the number of reactions. The Kronecker delta δ k,i permits to take into account the reaction rate if species k was involved in reaction i. Finally, the heat released by the chemical reactions was modelled by the term . Q, which included the contributions of all the reactions. It was defined as follows: where ∆H fm,i is the heat of formation of species i.

A Detailed Turbulence Model
In order to describe the turbulence phenomena during the combustion in a realistic way in the model, LES approach was considered in this work. LES models make use of the filtered local volume-averaged conservation equations (FAVRE approach) [42] to solve the flow field, being the small-scale stresses solved with a subgrid-scale model [20] due to the low dependence of these scales on the geometry. In particular, the LES dynamic subgrid-scale kinetic energy model (LES k sgs Eqn.) was used [21] in this work. This model was specially designed for compressible flows. In this model, the subgrid turbulent kinetic energy (k sgs ) is defined as followed: It is calculated making use of the transport equation: The subgrid-scale eddy viscosity is modelled as follows: where C k is a constant and V 1/3 is the local grid scale calculated from the cell volume V in each cell as V = (∆x∆y∆z) for inhomogeneous grids. The subgrid-scale turbulent stress tensor is calculated as follows: where dev S ij is the deviatoric component of the rate-of-strain tensor for the resolved scales. This way, the model relates the subgrid-scale stresses τ ij to the large-scale strain-rate tensor S ij . This LES model showed good results when dealing with compressible flows [22] and gas reacting flows [18].
Units are J, kmol, cm and K. Step

The Turbulent Combustion Model
The modelling of the combustion mechanism under the turbulent regime is a challenging physical problem that usually requires high computational costs as it must solve different time and spatial scales of a turbulent flame. In this work, the model chosen for modelling the turbulent combustion was the artificially thickened flame model (TFM). This model introduces an F factor in the energy and species equations of the gas phase that affects the thermal and molecular diffusivities (see Equation (1)). On the one hand, the F factor multiplies the pre-exponential factor of the kinetic equations; this permits to decrease the reaction rates by that factor. On the other hand, it increases the molecular diffusivity by the same factor. As the laminar flame speed is proportional to both mag- ω), the model provides a flame which propagates at the same speed. Besides, the modelled flame is F times thicker as the laminar flame thickness is a function ω . This way, the computational requirements of the mathematical model are relaxed, and less demanding grid sizes are required [51]. However, this approach modifies the physics of flame propagation since the Damköhler number is reduce [52]. This drawback is solved by taking into account the efficiency function E ∆ that takes into account the actual wrinkling of a turbulent flame by introducing subgrid wrinkling of the modelled flame. In this study, we used the efficiency function proposed by Charlette [53,54]. This function was calculated as follows: , Re ∆ is a function of the turbulent intensity u ∆ at the scale of the test filter scale ∆, the subgrid-scale turbulent Reynolds Re ∆ and the laminar flame thickness δ 0 L . The subgrid-scale turbulence intensity u ∆ was obtained from the obtained velocity resolved at the ∆ mesh scale as Reynolds number at the subgrid scale was estimated as follows: Re ∆ = u ∆ ∆ ν . The flame laminar flame thickness at each cell was estimated following Charlette [53] procedure with the relationship δ 0 Regarding the exponential factor, in this work, we used the fixed value of β = 0.5 proposed by Charlette [53].

Numerical Methods
The finite volume approach was used to numerically solve the system of equations. After testing different numerical integration strategies, the numerical procedure that provided the best results was selected for the validation of the model. The numerical schemes chosen were AUSMup-HLLC Low Mach [55][56][57], Godunov Scheme for the gas phase and a flux-difference splitting scheme (Rusanov) [58] for the particle phase, both with a Van-Leer [59] TVD scheme. Time integration was performed with the classical four-stage Runge-Kutta scheme for the fluxes, inter-phase, turbulent and chemistry source terms. Thus, the fluxes and source terms involved in each transport equation were evaluated in multiple substeps per each fluid-convection timestep. Primitive variables were then reevaluated from the intermediate conservative variables evaluated for each Runge-Kutta slope evaluation substep. Fluxes, diffusive terms and source terms were then recalculated from the corresponding intermediate primitive values. Although a higher number of operations per timestep is needed, this scheme permitted to set a higher CFL number with a more stable fluid flow behavior during the simulation, avoiding the presence of numerical instabilities that might be encountered in the simulation results [7].
A second-order linear Gauss scheme was considered for spatial discretization, including gradient, divergence and Laplacian calculations. This numerical strategy was previously tested for hydrogen premixed turbulent combustion problems with good results [23]. Regarding the solid phase, a Rusanov scheme [58] was used to evaluate the convective fluxes in the conservation equations of the solid phase. This approach was previously validated for gas-particle combustion problems [7,8,14,16,17].
The integration of the ODE equation corresponding to chemistry source terms was a stiff problem difficult to solve in a cost-efficient way. After testing different strategies, an in situ adaptive tabulation method (ISAT) with ODE(SEULEX) integration [60,61] was used to solve the system of equations of the chemical kinetics of the gas phase and to estimate the reaction rate of the k th species in the i th reaction ( . ω k,i ). ISAT is a method originally proposed for turbulent reacting flow simulations [60]. This method aims to approximate non-linear system solutions by means of a set of linear regressions of independent variables from a lookup tabulated database constructed dynamically with previous solutions (storage and retrieval method). Thus, it permits to reduce the number of ODE integrations for the chemistry set of ODE equations, being in the problem analyzed in this work one of the most computation-demanding tasks per integrated timestep. It has been reported that using this technique allows, under certain conditions, decreasing by three orders of magnitude the computer time required to compute the detailed chemistry in reacting flow computations [60]. This algorithm has been successfully applied in combustion chemistry problems involving up to 50 species [62] and was also used for premixed H 2 -air combustion with detailed kinetics and LES-TFM modelling, similar approaches to the presented in this work [23,63]. A relative tolerance of 10 −4 was set as the threshold for retraining the tabulated dataset. It is worth mentioning that this method is being extended to the applications other than the initially intended, especially for real-time predictive control [64] as an alternative to neural networks since it presents some advantages, e.g., it does not need training data before use.
The SEULEX ODE integration consists of a semi-explicit multistep method based on numerical extrapolation. An absolute tolerance 10 −9 has been set, limiting to 1000 the maximum number of iterations per chemical gas integration.
Regarding the particle phase, a first-order implicit Euler scheme was used to integrate the chemistry equations with good results. Reaction rates for each timestep were evaluated using a sub-timestep (chemical timestep) in which each reaction rate for a given species and reaction is updated by the previous sub-timestep species concentrations. The temperatures used for evaluating the reaction rates were also updated from the last sub-timestep. An initial timestep of 10 −12 s was set, thus avoiding numerical fluctuations and divergences that may lead to nonphysical results.

Results and Discussion
In the previous section, a two-phase flow model for turbulent combustion of gas and particles mixtures was presented. The use of LES, TFM and detailed kinetic schemes was explored as a way to take into account a realistic description of turbulence, flame wrinkling and reaction mechanisms in the turbulent combustion process. In order to validate the present model and evaluate its prediction capabilities, the experimental results from combustion tests presented by Sabard and Sabard et al. [65,66] were used as a reference benchmark. These experiments were performed at CNRS Orleans (France) and assessed the combustion of gas mixtures of H 2 -O 2 -N 2 and graphite (carbon) particles in a spherical bomb. The experimental facility consists of a spherical bomb with the internal radius of 125 mm equipped with different instrumentation including piezo-resistive wall pressure sensors, a Schlieren system, an electrical (spark) ignition system as well as a laser-induced ignition system (Figure 1). A detailed description of the experimental setup and methodology can be found in the works by Sabard and Sabard et al. [65][66][67]. During the experiments, different concentrations of C, H 2 and O 2 were introduced in the spherical bomb. The uncertainty in the volumetric composition of the gas mixture was below 0.3%, whereas the wall pressure measurements obtained had an experimental uncertainty smaller than 2%. The ignition of the mixture was generated with an electric spark between two electrodes. The initial pressure and temperature within the sphere were 1 bar and 298 K, respectively. The experimental conditions of the tests were as follows:  These experiments were simulated with the model presented in the previous section. The numerical domain defined for the simulations consists of an eighth part of a sphere with the radius of 125 mm. Thus, three different symmetry planes were considered. A grid sensitivity study was performed in order to check the potential influence of the spatial discretization simulated on the results. The results showed that there was a small influence of discretization in the radial direction when the mesh size was under 125 µm (i.e., 1000 radial elements). Independency of the opening angle was also assessed, reporting less than 0.3% in the variation of the maximum pressure and less than 0.115 ms differences in the rising time between the meshes with different angles. The final mesh size of 125 microns was used. Finer meshes had no influence on the chemistry mechanism and fluid fields. Regarding the initiation, it was assumed that an autoignition of the mixture which affected a small sphere with the radius of 2.5 mm initiated the sequence. This initiation was modelled as an addition of energy of 850 kJ/m 3 applied to the ignition volume in 0.1 ms. This quantity of energy was enough to initiate a laminar flame in the domain.
The results of the experimental benchmark of the two-phase model are presented in Figure 2. The figure shows a comparison of the model prediction of pressure evolution at the wall of the spherical bomb with time with the experimental data obtained in Experiment 1 and Experiment 2. For the sake of comparison, an experiment performed without particles (called in the figure H-EXP2) was also included [66,67]. Its conditions were similar to Experiment 2 except for the presence of particles. These experiments were simulated with the model presented in the previous section. The numerical domain defined for the simulations consists of an eighth part of a sphere with the radius of 125 mm. Thus, three different symmetry planes were considered. A grid sensitivity study was performed in order to check the potential influence of the spatial discretization simulated on the results. The results showed that there was a small influence of discretization in the radial direction when the mesh size was under 125 µm (i.e., 1000 radial elements). Independency of the opening angle was also assessed, reporting less than 0.3% in the variation of the maximum pressure and less than 0.115 ms differences in the rising time between the meshes with different angles. The final mesh size of 125 microns was used. Finer meshes had no influence on the chemistry mechanism and fluid fields. Regarding the initiation, it was assumed that an autoignition of the mixture which affected a small sphere with the radius of 2.5 mm initiated the sequence. This initiation was modelled as an addition of energy of 850 kJ/m 3 applied to the ignition volume in 0.1 ms. This quantity of energy was enough to initiate a laminar flame in the domain.
The results of the experimental benchmark of the two-phase model are presented in Figure 2. The figure shows a comparison of the model prediction of pressure evolution at the wall of the spherical bomb with time with the experimental data obtained in Experiment 1 and Experiment 2. For the sake of comparison, an experiment performed without particles (called in the figure H-EXP2) was also included [66,67]. Its conditions were similar to Experiment 2 except for the presence of particles.
As shown, the two-phase approach proposed with LES and a detailed chemistry kinetic model was able to predict with good results the pressure evolution with time including the fast pressure rise found at the wall between 10 and 20 ms. As indicated in Table 3, the model permitted to predict the maximum wall pressure (P max ) and the time lapsed to reach the maximum pressure (t rise ) with the relative error of 3.3% and 20.8%, respectively, for Experiment 1 and 7.8% (in P max ) and 18.2% (in t rise ) for Experiment 2. Taking into account the uncertainties linked with the ignition process (which are of the milliseconds order), the prediction of the time locus of the maximum pressure by the model could be considered satisfactory.
the model permitted to predict the maximum wall pressure (Pmax) and the time lapsed to reach the maximum pressure (trise) with the relative error of 3.3% and 20.8%, respectively, for Experiment 1 and 7.8% (in Pmax) and 18.2% (in trise) for Experiment 2. Taking into account the uncertainties linked with the ignition process (which are of the milliseconds order), the prediction of the time locus of the maximum pressure by the model could be considered satisfactory.  The model's capacity to predict particle concentration effect on the combustion sequence can also be assessed by taking into account C-EXP2 and H-EXP2. As shown in Figure 2 (down), the model was able to predict the reduction of the maximum wall pressure when the particle concentration was reduced to zero. The model also showed a slight time displacement of the curve when the particle concentration was reduced. This displacement resulted in the increase of the combustion time but with approximately the same dP/dt. This experimental tendency was also well-captured by the model qualitatively. However, the experimental increase of the combustion time was larger than the one predicted by the model.
In general terms, Figure 2 shows a good agreement of the model with the experimental data. The main deviations between the simulations and the experiments are related to the final stage of the combustion process (i.e., for t > 20 ms) when the pressure level reached is maintained. These deviations can be related to the uncertainty found in the chemical model at high temperatures and to the influence of the graphite particle size. Therefore, the range of applicability of this model can be set on the basis of the conditions used in the validation, that is, H 2 concentration in the air of 20% for mixtures of N 2 /O 2 between 2.33 and 3.76 at the initial ambient pressure and temperature, graphite particle size of the order of 35 microns and particle concentration between 0 and 97 g/m 3 .
As for the prediction of the combustion products, Table 4 shows a comparison of the concentration of CO percentage in the combustion products estimated by this model and the one experimentally measured in the test C-EXP1 and C-EXP2 [66,67]. The table also includes the numerical prediction estimated with the Cosilab software for the same tests [66]. The results also permitted to compare the effect of different modelling approaches in the oxidation mechanism considered in the solid phase (porous vs. nonporous approximation). As shown, the mathematical model presented in this study provided better predictions of the CO composition in the combustion products than the Cosilab software in the case of both porous and nonporous modelling approaches. The table shows a good estimation of the porous model proposed, with deviations of less than 1.5% in the prediction of the volumetric percentage of CO in the combustion products and less than 0.1% in the case of the nonporous approach, whereas Cosilab provided a minimum deviation of 11%. The porous model provided results slightly closer to the experiments than the nonporous model in the prediction of gas combustion products. The comparison of the simulation results with both kinetic models (porous and nonporous) and its comparison with the experimental data of C-EXP1 and C-EXP2 also highlighted that the actual particle porosity was an important factor to predict transient subproducts of the combustion sequence. In fact, the consideration of the porous or the nonporous chemical kinetics model resulted in differences in the reaction rates of two order of magnitude. For example, in the case of the reactions C + 0.5 O 2 → CO and C + CO 2 → 2 CO, there were found differences of the order of 100 between both models. This was somehow expected as the particle porosity is directly linked to the actual effective particle surface which is directly linked to the reaction rates of the solid phase (Table 1). However, the actual particle porosity in industrial and safety problems commented upon above is usually unknown and must be considered from the practical point of view as an uncertainty factor of the model. All in all, it has to be remarked that regarding the prediction of P max and t rise both, the porous and nonporous chemical kinetics models presented the similar behaviour without remarkable quantitative differences in the dynamic evolution of wall pressure with time.
As for the applications of this model, it can be used to evaluate the efficiency of different strategies for tuning the products of the H 2 -graphite combustion process to the desired conditions depending on the context.
In the nuclear safety context, a reduction of the combustion velocity and/or a quenching of the combustion process would be desirable. In this case, the strategy was to reduce the concentration of H 2 and/or O 2 in the scenario by promoting catalytic recombination that would reduce the probabilities of the deflagration-to-detonation transition. Under these circumstances, passive autocatalytic recombiners (PARs) are used to promote the reaction H 2 + O 2 → H 2 O. In this case, there are undesirable subproducts in the H 2 -graphite combustion whose concentration should be reduced to the minimum possible in order to avoid the poisoning of the catalyser. Specifically, because of the large sticking coefficient of CO compared to the other adsorbed species within a PAR and its high activation energy for desorption, the presence of CO in the mixture would poison the catalytic surface; this would prevent the desired recombination reaction in the PAR from occurring [68]. Besides, the lean limit concentration of hydrogen combustion decreases as the CO concentration increases and the flammable region widens for H 2 -CO mixtures. Thus, CO is an undesirable by-product in the overall reaction [69]. Thus, in the safety context, the present model could be a useful tool to evaluate by means of numerical simulation the efficiency of different strategies used to mitigate the potential hazard during an accident sequence. For example, it can be applied in the prediction of accident sequences in ITER or in containment buildings of nuclear power plants, to adjust the parameters of mitigation systems such as N 2 injectors or passive autocatalytic recombiners.
In the industrial context, this numerical model can be applied in syngas combustors to evaluate potential strategies to reduce the concentration of H 2 , CO and solid C at the exit of the combustor. For example, the development of IGCC technologies, involving gas-turbine combustion of syngas, derived, for instance, from air or O 2 gasification of pulverized coal, has recently promoted interest in studies of CO/H 2 combustion. In this case, this numerical model can be used to improve the efficiency of these combustors by predicting quenching, flame acceleration and/or spatial regions where the combustion process is inefficient and might lead to an increase in undesirable by-products.

Conclusions
In this work, we presented a numerical model for describing the turbulent combustion of two-phase flow mixtures of gas and particles. Specifically, we analysed the influence of the presence of solid carbon particles in the turbulent combustion of an H 2 -air atmosphere. A two-phase model was proposed to describe this reacting flow with LES and detailed chemistry. The model proposed was benchmarked against experimental combustion data obtained in a spherical bomb. The results highlighted some significant specifics:

•
In case of highly diluted mixtures of H 2 -air and graphite particles, the benchmarked results showed that LES with a detailed chemistry model were found to be an appropriate engineering approach for analysing premixed turbulent combustion of graphite-H 2 mixtures.

•
The validation against the experimental data show that the two-phase approach used in the present model based on the Eulerian-Eulerian approach seems to be accurate enough to afford this type of combustion sequences with highly diluted mixtures.

•
Under the conditions studied, the model captured well the key tendencies linked to the presence of carbon particles of the microns order. In this sense, the model was able to predict that the presence of a low concentration of carbon particles (of the order of 96 g/m 3 ) accelerated the combustion sequence, obtaining smaller combustion times than in the absence of particles and larger maximum wall pressure levels (of the order of 15%). • Classical graphite and hydrogen detailed oxidation mechanisms [30,31,33] coupled with a Eulerian-Eulerian model provided good results in the prediction of combustion products under turbulent combustion conditions. Regarding graphite combustion, the porous oxidation model provided results closer to the experiments than the nonporous model.
Future work will face the extension of this model to other metal particles such as tungsten and to engineering applications of hydrogen turbulent combustion sequences.

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

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Nomenclature
A Representative area of a particle.  . ω m,k Mass reaction rate for species k in phase m.