Keywords: Partial Differential Equation PDR, Porosity Distributed Resistance RANS, Reynolds Averaged Navier-Stokes

The use of CFD (Computational Fluid Dynamics) software for the numerical prediction of difficult experiments such as the consequences of gas explosions in industrial environments remains a major challenge in process engineering. From the state of the art in this area of research where numerical simulation results have been compared to valid experimental results, it is concluded that these software can not contribute to improve safety (the differences between experimental results and numerical results are important).

However, given the potentiality of the field of Computational Fluid Dynamics, it would probably be unreasonable to marginalize it in hazard studies.

This thesis helped to define strategies for estimating the consequences of explosions by CFD.

Part of the work carried out consisted in determining the most common equations, techniques, models and methods in the most used CFD software (in the framework of study concerning the prediction of the consequences of an explosion in industrial environments).

The URANS technique (Unsteady Reynolds Averaged Navier-Stokes) was chosen for the numerical resolution of the fundamental laws of fluid mechanics. The kepsilon turbulence model and one of its variants (the low Reynolds number kepsilon model) were selected for the turbulence simulation. The modeling of reactive flows is established using the CREBCOM model (CRiteria and Experimentally Based COmbustion Model) and the EDM model (Eddy Dissipation Model). The finite volume method has been used for the discretization of continuous equations (the laws of fluids mechanics and associated turbulence and combustion models). The convective terms of these equations are solved using the numerical scheme of Roe and that of Van-Leer. The discretization of the diffusive terms was carried out using classical centered schemes. The unsteady terms are solved using the explicit Euler method. The selected meshes are of type Voronoi and of type structured and re-cutting using the AMA technique (Anisotropic Mesh Adaptation). Boundary conditions used are mainly Dirichlet type and Neumann type.

To go beyond the user aspects and have a perfect control of the software used (only pledges of a good analysis of the physical and mathematical content of CFD tools), a CFD software called MERLIN has been fully developed. It contains all the equations, techniques, models and methods selected previously and was used for all numerical simulations performed in this thesis. In order to ensure the reliability of the experiments performed with MERLIN, its verification was carried out using the MMS (The Method of Manufactured Solutions) method.

To understand the numerical representation of the physical phenomena associated with the phenomenon of the explosion, we first carried out a study on shock wave propagation in different configurations (subsonic case of the Sod'shock tube problem, reflection of an unsteady shock on a compression ramp, supersonic flow on a rising step). As a result, the accuracy of a shock wave structure predicted numerically depends on the numerical scheme and the type of mesh used. The choice of the numerical scheme and type of mesh depends on the type of shock wave to be simulated numerically. The second study performed consisted in simulating the dispersion of gas. This numerical experiment revealed that a good numerical approximation of gas dispersion is independent of the numerical scheme chosen but rather relies on the turbulence model and the type of mesh used. The choice of the turbulence model is relative to the presence or not of confinement and the type of mesh depends on the flow. The last study carried out concerns the simulation of flame propagation in different configurations (confined and unconfined environments, with and without obstacles). It demonstrates that the accuracy of a numerical prediction of the flame behavior flame is based on the choice of the combustion model and the type of mesh.

Ultimately, an optimal numerical estimation of the consequences of an explosion requires a good knowledge of the physics of the phenomena highlighted, which will help to make a reasonable choice of numerical scheme, type of mesh, model of turbulence and model of combustion.

INTRODUCTION

Virtual reality, including tools issued from the Computational Fluid Dynamics (CFD) domain, is spreading in many fields of Chemical Engineering especially in fields where "pilots" are difficult to design. This is the case for instance for some acute safety problems such as large explosions covering the entirety of a plant. The reader could find many examples of disastrous industrial explosions in the literature. Because of the large potential impact of an explosion, the user of an explosion simulation tool would like to know how reliable the results of the simulation will be.

Investigating the potential source of uncertainties of engineering CFD code in such a key aspect of chemical engineering ("industrial explosions") constitutes the red line of this PhD work.

Remember that an explosion occurs when a flammable mixture of a combustible gas in air is rapidly burned by a flame propagating in this cloud. The explosion which occurred in the oil depot of Buncefield [START_REF] Makarov | An inter-comparison exercise on CFD model capabilities to predict a hydrogen explosion in a simulated vehicle refuelling environment[END_REF] may not be the most important amongst the disastrous industrial explosions, although quite commensurate (43 injured, hundreds of homes evacuated, cost of damages amounting £1billion 1 ), but it shed light of the limitations of the state of the art about gas explosions risk assessment tools. On this industrial estate, the explosion risk was identified before the accident occurred. In many places where the congestion level is rather low, conventional risk analysis method, like the Multi-Energy method, would have significantly under-predicted the overpressures over most of the areas affected by the explosion. It was claimed that such phenomenological tools cannot take into account of the complexity of the site and that "complex" CFD (computational fluid dynamics) software would perform much better. But the poor predictability of the genuine method might also come from the complexity of the physics and it is not totally obvious that CFD codes, at least those usable for large scale explosions (e.g. in potentially complex geometries developing over tens of meters), would have the intrinsic potentiality to perform better as "benchmarking" exercises reported below tend to suggest.

CFD softwares applicable to large scale explosions are often part of commercial numerical platforms. FLACS for example, quite popular in the industrial world mainly because of its ergonomics, is a commercial numerical platform of this type. But there are others like EXSIM, REACFLOW, CAST3M, AUTOREAGAS, FLUENT ... It is worth remembering that the core (choice of equations and numerical method) of these tools was developed about thirty years ago and was more or less preserved as such until now.

The "natural" way to investigate the abilities/limitations of these tools to predict the consequences of an explosion is perhaps to make inter-comparisons on test cases or experiments. Several exercises of this type have been made over the past decades, not all of them being public (1), ( 2), (3).

1 http://news.bbc.co.uk/2/hi/uk_news/england/7777539.stm One of the first was performed on a realistic industrial situation. This was a blind inter-comparison exercise (modelers were not aware of the results of the tests but had the conditions of the experiments) concerning an explosion inside an offshore platform module 2 .

The results are presented on Figure 0-1 and reveal a very large scattering from a tool to another. Why is it like this? Based on even simple physical grounds, it can be shown that the explosion overpressure is strongly coupled to the flame dynamics (velocity, acceleration) so that the first natural conclusion is that the flame behavior may not correctly handled in the simulations.

To go deeper in this analysis, the European excellence network HySAFE organized a significant benchmarking exercise on that aspect. Garcia (3) published the results focusing on the flame propagation problem. The authors compared their simulations (CAST3M, COM3D, REACFLOW, AUTOREAGAS, FLACS, FLUENT, COBRA, CFX-4, NEWT) to the results of an experiment where an explosion was triggered inside a hemispherical balloon of 2000 m 3 filled with a homogeneous hydrogen-air mixture (Figure 012345). Mostly the designers of the codes who are in principle those managing their code at best (so no "standard users") participated in this exercise. Moreover, the experimental results were available to them before starting the simulations. Aside from the physical issues, the "user dimension" may also be tracked considering in such benchmarking exercises for instance with FLACS code, widely used in the industry. NH is such a user (well informed) and Gexcon is the designer. Using the same code, the results vary again by a factor of 2 between Gexcon and NH.

Partly based on such elements, a large national institute, (HSL), concluded that the CFD codes may not the best suited to the exercise of "explosion" prediction [START_REF] Lea | A Review of the State-of-the-Art in Gas Explosion Modelling[END_REF]. Perhaps the ground reasons are that, not only the predictability is not insured, but also that uncertainties are not mastered and that even well-informed users may not be able to make a reasonable job.

But it is certainly not good to do without (or ignore) the potential advantages of the vast domain of the computational fluid dynamics and it seems more appropriate to question the shortcomings of the present CFD solutions, in this particular context of large scale explosions simulations, and to look for potential improvements. This PhD is a contribution to this.

From the preceding analysis of past "codes benchmarking", it seems that the benchmarking exercise could be an endless task and may not be even the best way to progress towards a better "control" of the uncertainties of CFD simulations.

There is a need to look for the reasons of the limitations and to investigate their impact on the simulation work. In the present work, it is proposed to address this important question, first, from an analysis of the physics (models) embedded into the codes used to simulated large scale explosions (especially on this physics relates to the phenomena to be simulated), and second, from an investigation on the way the codes may be operated by engineers. About the embedded physics, it seems rather clear that the core of the difficulty is how the flame is propagated, especially how it interacts with the geometry and with the flow. About the "user dimension", most of the difficulty relates to the degree of awareness of the user about the underlying physics (of explosion in our case) which will drive his "modelling" choices i.e. the way the engineer conceptualizes the situation to be modelled. Some key aspects about this are raised in this work.

In chapter 1, the most important physical aspects about flame propagation and explosion development within the context of industrial explosions are presented. The fundamentals of Computational Fluid Dynamics are also recalled targeting more specifically the CFD codes dedicated to the simulation of large scale explosions. Rather than testing each code as they stand, which would have been very long and perhaps insufficient as previous benchmarking showed, a numerical toolbox, named MERLIN, was developed containing models, equations and numerical routines found in many CFD engineering codes and in particular in the "explosion" codes. This toolbox is described in chapter 2. Chapter 3 is the application of MERLIN to the simulation of the succession of events leading to a large scale explosion in the industry following rather traditional risk analysis approaches. Doing so, it is attempted to highlight the influence of the physical models contained into the software and of the modelling choices of the user.

1.1 PHENOMENA

COMBUSTION, FLAME AND PRESSURE GENERATION

Traditionally, an explosion results from a sudden energy release leading to a fast expansion of gases. This expansion is responsible of the pressure and disruption effects (Figure 1-1). An explosion can occur once an ignition source is introduced into an explosive atmosphere composed by a fuel (gas, vapour...) and an oxidizer (oxygen for example) intimately mixed (forming of a cloud) in appropriate proportions (inside the explosivity domain). In contact with "an efficient" source of ignition (strong spark for instance), the combustion starts, and hot combustion products are formed, which temperature is typically between 1000 and 2000 ° C (Lewis and von Elbe [START_REF] Lewis | Combustion, Flames and Explosions of Gases[END_REF]. These hot burnt products can then act as a «source of ignition» for the surrounding layers of mixtures so that a sort of "ignition-combustion wave" (the "flame") propagates by itself throughout the cloud. On its way, the flame turns "cold" reactants medium (20°C) in "hot" combustion products (1000-2000°C). A given portion of a cloud traversed by the flame undergoes a thermal expansion (the volume after combustion becomes 5 to 10 times larger more important) [START_REF] Daubech | Contribution à l'étude de l'effet de l'hétérogénéité d'un pré mélange gazeux sur la propagation d'une flamme dans un tube clos[END_REF].

If this transformation occurs in a gaseous volume V 0 starting at an initial pressure P 0 , it can be shown from the first principle of thermodynamics that the evolution of pressure P and volume V derives from the overall thermal power Q issued from the flame (γ is the ratio of specific heats, t the time and "e" the internal energy of the gases):

de dt = Q -P.
dV dt [START_REF] Makarov | An inter-comparison exercise on CFD model capabilities to predict a hydrogen explosion in a simulated vehicle refuelling environment[END_REF] And de dt = d(m. Cv. T) dt . [START_REF] Baraldi | An intercomparison exercise on CFD model capabilities to simulate hydrogen deflagration in a tunnel[END_REF] Where m stands for the mass of the gases and C the specific heat at constant volume. With the following additional laws:

C = R M. (γ -1)
and P. V = m. R M . T [START_REF] García | An intercomparison exercise on the capabilities of CFD models to reproduce a largescale hydrogen deflagration in open atmosphere[END_REF] (M is the molar mass of the gases), after some manipulations, we obtain:

V γ -1 . dP dt + γ. P γ -1 . dV dt = Q [4]
At each instant, the pressure increase dP depends on the value of Qf at this same time, which means that the overall shape of the pressure signal is directly related to variations of flame power Q .

Q depends on the rate at which the reactants are consumed. If at this particular time, the flame area is A 𝐟 (t), the volumetric rate per unit flame area with which the flame consumes the reactant is S u (homogeneous to a speed, this is the "burning velocity" depending on the mixture and turbulence of the flow as illustrated later) and the volumetric mass is ρ, Q reads :

Q = ρ. A (t)
. S . Q . [START_REF] Popat | Investigations to improve and assess the accuracy of computational fluid dynamic based explosion models[END_REF] Where Q is the amount of thermal energy released by the combustion per unit mass of the mixture effectively used to heat up the gases. Consider for instance the generic combustion reaction ("f'" and "F" for "fuel" and ("p'" and "P" for "products"):

ν . F + ν . O + ν . Inert → ν . P + ν . Inert. [START_REF] Lewis | Combustion, Flames and Explosions of Gases[END_REF] In most situations, the inert part constitutes the major fraction of the mixture (nitrogen for instance for combustions in air). If all the components are gaseous and the reaction fast enough, the energy is conserved (enthalpy) so that:

Q = Cp. (T -T ) = y M . ν . h _ + ν . h _ -ν . h _ = y M .
∆H . [START_REF] Daubech | Contribution à l'étude de l'effet de l'hétérogénéité d'un pré mélange gazeux sur la propagation d'une flamme dans un tube clos[END_REF] Where y is the mass fraction of the fuel in the unburnt mixture (supposed here to be the deficient species), M f the molar mass of the fuel, h f_i the standard enthalpy of formation of the species i, C the specific heat of the mixture and ΔH comb the standard heat of combustion of the fuel. T b is the final temperature of the burnt mixture and T f the temperature of the reactants (which may be different from the initial temperature if the system due to compression is any). This last equation tells that Q is an intrinsic property of the mixture, depending only on its chemical composition. Another way to write the expression is:

Qcomb = Cp. T T T -1 = Cp. T (α -1) = y M .
∆H . [START_REF] Proust | Résistance des structures de l'industrie agroalimentaires aux éffets des explosions de poussières[END_REF] Where α can be seen as the ratio of the specific masses of the reactants and burnt products (see the next section for the demonstration) or, in other words, the expansion ratio of the burnt products. Knowing that:

Cp. T = γ. R M . (γ -1) . T = γ. P ρ. (γ -1) . [START_REF] Solberg | Observations of flame instabilities in large scale vented gas explosions[END_REF] Where M is the molar mass of the unburnt mixture, expression [START_REF] Lewis | Combustion, Flames and Explosions of Gases[END_REF] finally becomes:

V γ . dP P. dt + dV dt = A (t). S . (α -1).

[10]

Thus, in any explosion phenomena, two physical domains interact strongly: the flame dynamics, represented by A (t) and S and the mechanical energy production (pressure, volume) represented by dP/dt and even dV/dt. Clearly a good estimation of the flame dynamics is required to be able to represent the consequences of explosions. This is the central piece of most explosion modelling activities.

To better illustrate this, two limiting cases might be investigated: the confined explosion as those occurring inside a strong enclosure for which the second term of the left member in [START_REF] Cooper | On the mechanisms of pressure generation in vented explosions[END_REF] is zero and unconfined explosions for which the first term is small as compared to the second one. A schematic view is presented in Figure 1-1-a.

In the "enclosure"(Figure 1-1-b), the speed of propagation of the flame is sufficiently low (<30 m / s) in order that the internal pressure remains uniform (in space 3 ). This condition is generally satisfied when the ratio between the largest and the smallest dimension of the device is less than 5 (H / D < 5) [START_REF] Proust | Résistance des structures de l'industrie agroalimentaires aux éffets des explosions de poussières[END_REF].

From expression [START_REF] Cooper | On the mechanisms of pressure generation in vented explosions[END_REF], the pressure due to the flame propagation reads:

1 P dP dt ≈ γ S A (t)(α -1) V [11] The incidence of the discharge can simply be estimated by superposing a rate of pressure decrease due to the outflow:

1 P dP dt ≈ γ Q V [12]
The volumetric outflow Q outflow is roughly:

Q = A .
2. (P -P ) ρ [START_REF] Van Den Berg | Vapor cloud explosion blast modelling ,Experimental investigation of the key parameters and blast modelling[END_REF] The maximum overpressure ΔP max = P max -P 0 is reached when expressions [11] and [START_REF] Kuhl | Al-Air Combustion in Unconfined SDF Explosions22nd ICDERS[END_REF] are equal and:

∆P ≈ ρ 2 .
S A (t)(α -1) A [START_REF] Leyer | Unconfined deflagrative explosion without turbulence[END_REF] which clearly shows that the pressure effects depend on the square of the burning velocity, expansion ratio and flame area. As shown later, the expansion ratio does not vary much and is mostly dependent on the fuel-air mixture. Similarly, the flame area is constrained by the outer walls and maximum flame area is limited by the cross section of the enclosure (≈ V 2/3 ). The more "sensitive" parameter is the burning velocity which is strongly impacted by the turbulence of the flow and, especially in a confined environment by the acoustics and acceleration, deceleration of the flow caused by the "channeling", by the walls and obstructions if any. Several of these aspects were identified in experiments ( 9), ( 7), [START_REF] Cooper | On the mechanisms of pressure generation in vented explosions[END_REF]. These phenomena can be interpreted looking the flame front as an "unstable" interface as shown later in this chapter.

Consider now the unconfined explosion situation (Figure 1-1-a). Suppose a flammable cloud in the open air, ignited at some point inside. A flame develops around, and, due to the volumetric expansion of the burnt gases, pushes away the surrounding atmosphere. This creates a pressure wave propagating outwards.

To model this, an "acoustic" model can be applied. The principles of this physical representation of the propagation of pressure waves in the environment were laid by Taylor (11) but were subjected to significant developments only in the beginning of the 1970's [START_REF] Kuhl | Al-Air Combustion in Unconfined SDF Explosions22nd ICDERS[END_REF], which led to practical tools in the 1980's ( 13), ( 14), [START_REF] Bjerketvedt | Gas explosion handbook[END_REF], [START_REF] Baker | Explosion Hazards and Evaluation[END_REF].

Leyer and Deshaies [START_REF] Deshaies | Flow Field induced by unconfined Spherical accelerating Flames[END_REF] made the most significant developments. An analytical model of the first order has been developed which predicts the propagation of the pressure wave in the environment. At a certain distance from the front, an acoustic solution is proposed:

∆p(r , t) = ρ (1 -α ) 4πr d V dt [START_REF] Bjerketvedt | Gas explosion handbook[END_REF] with [START_REF] Baker | Explosion Hazards and Evaluation[END_REF] where:

∂ V ∂t = ∂ ∂t (A × S )
 ∆p(r, t) is the pressure at time t and at the distance r ρ is the density of the atmosphere  α is the expansion rate of the combustion products  V is the volume of the "acoustic source" (flame ball)  r is the position of the observer relative to the ignition source  A is the flame area  S is the spatial flame velocity which is about S u .α. This model was validated on the basis of small-scale experiments [START_REF] Leyer | Effets de pression engendrés par l'explosion dans l'atmosphère de mélanges gazeux d'hydrocarbures et d'air[END_REF] which applicative version is the multi-energy method [START_REF] Berg | the Multi-Energy Method -a framework for vapour cloud explosion blast prediction[END_REF]. It is often admitted that the flame propagated spherically (radius R f at time t) so that with: A = 4. π. R and S = dR dt ≈ S . α [START_REF] Deshaies | Flow Field induced by unconfined Spherical accelerating Flames[END_REF] Expression [START_REF] Bjerketvedt | Gas explosion handbook[END_REF] becomes:

∆P(r , t) = ρ . (1 -α ). R r 2. S + R . dS dt . [START_REF] Leyer | Effets de pression engendrés par l'explosion dans l'atmosphère de mélanges gazeux d'hydrocarbures et d'air[END_REF] which again shows that the pressure effects are extremely sensitive to the burning velocity and even to its variations. Turbulence of the flow and flame instabilities need to be considered. There has been in particular a long lasting debate about the way obstacles interact with an outward expanding flame in the open atmosphere. Do they simply increase transiently the flame area, the flame "wrapping" around them, do they promote the turbulence in their wake due to the expanding flow or do they trigger some flame instabilities due to the deceleration of the flow upwind? This question is very important in terms of how to model them in CFD.

So not only flame propagation but more precisely, flame dynamics should be the central core of the modelling of explosions. To understand the most relevant aspects, the physics of flames propagating in premixed fuel-air medium needs to be presented.

"PREMIXED" FLAME DYNAMICS

PREMIXED FLAME AS A "COMBUSTION" WAVE

Note first that pressure effects would not exist if the burnt gases were not submitted to a volumetric expansion. As shown by expression [START_REF] Proust | Résistance des structures de l'industrie agroalimentaires aux éffets des explosions de poussières[END_REF], the expansion ratio α is linked to the combustion temperature, T b , which is also an important parameter in the flame propagation process. The second important parameter is the "burning velocity" which dictates how fast the explosion develops.

As shown before, the volumetric expansion rate of the fluid particles through the flame is a thermodynamic data that depends only on the amount of heat released from combustion, expressed, for example, through the application of the first law of thermodynamics:

α = ρ ρ ≈ T T = y M ∆H C T + 1 [19]
where  ∆H is the enthalpy of reaction at T per mole of fuel  T and T are respectively the temperatures of the combustion products and the reactants,  ρ and ρ are respectively the densities of the reactants and combustion products,  C is the average specific heat of the combustion products.

The volume expansion rate is then an intrinsic and fundamental parameter depending essentially on the composition of the mixture. For most fuels burning in air ∆H is on the order of 40 MJ/kg. At the stoichiometric composition where the maximum value of T b is expected, y f is typically 0,06 and Cp≈1100 J/kg°C. The expansion is then about 8 and T b about 2000°C. At the lower limit of flammability y f is about 0,03 and the expansion ratio drops to about 5 (T b =1200°C). Experiments and more detailed calculations confirm this (6).

The "burning velocity" is the velocity with which the flame consumes the reactants. As such it may be regarded as the volumetric consumption rate of the reactants per unit area of the flame front. Much of the present knowledge about premixed flame behavior in a number of circumstances rests upon the understanding of the burning mechanism inside the undisturbed premixed flame, sometimes called "laminar flame", although a quite idealized situation difficult to control even in laboratory condition [START_REF] Proust | Formation-Inflammation-Combustion des atmosphères explosives (ATEX) et effets associés[END_REF]. In those specific conditions, the combustion velocity is sometimes called "laminar burning velocity" or "fundamental flame speed", S . Mallard and Le Chatelier [START_REF] Mallard | Combustion des mélanges gazeux explosifs[END_REF] were the first to lay the foundations of a theory defining the laminar burning velocity (or fundamental burning velocity).

Mallard and Le Chatelier [START_REF] Mallard | Combustion des mélanges gazeux explosifs[END_REF] assumed that the flame progresses in the unburned mixture due to thermal conduction from the burnt and hot side of the flame which continuously ignites the mixture in the front face of the flame (Figure 12). They modelled the flame as consisting of two zones. In Zone I, heat is transferred by thermal conduction and the reactants are progressively heated up to the ignition point (T : ignition temperature). No chemical reaction occurs in this preheating zone. The reaction is triggered at T and proceeds until completion (T ) in zone II. The flame can propagate at a constant velocity only if the heat transmitted to the zone I is equal to the heat production rate in zone II. Because of this, Lewis and von Elbe described the flame as a "combustion wave".  The state equation:

(ρ ρ ⁄ ) = T T ⁄ [22]
where  a m is the molar concentration  Q is the reaction heat by mole of reactant  ρ is the density  ρ is the density of fresh gases and m * is the mass flux through the flame:

m * = ρ. S [23]
The Lewis number establishes a link between the energy and species balance equations so that both equations can be solved simultaneously:

Le = λ ρ. c D [START_REF] Ya | Mathematical Theory of Combustion and Explosion[END_REF] where D is the species diffusion coefficient.

Within the context of the self-heating theories, the chemical reaction term is expressed by a global Arrhenius law:

w * = Z . [a ] . exp - Ea R. T [25]
where  Ea is the activation energy  Z is the pre-exponential factor  a m is the molar concentration  N is the reaction order  R is the perfect gases constant Zeldovitch and al. [START_REF] Ya | Mathematical Theory of Combustion and Explosion[END_REF] articulate their reasoning about the behavior of the Arrhenius law of combustion reactions. They note that the activation energy is very high [Ea/R=O(10000)] so that the combustion reaction takes place with a temperature near of T on a zone with very small thickness. In a way, this situation justifies the model of Mallard and Le Chatelier with two separated zones, one chemically inert and the other where the combustion takes place. The additional point is that necessarily T is close to T . The demonstration is given below because the reasoning will be used again later.In the first zone (zone I), the energy equation is reduced to:

d T dx - m * . c λ . dT dx = 0
[26] with x = -∞, T = T ; x = 0, T = T [START_REF] Williams | Progress in knowledge of flamelet structure and extinction[END_REF] and in the zone II, with small thickness δr, the convection term is negligible so that:

d T dx + w * . Q λ = 0 [28] with x = -∞, T = T ; x = δr, T = T [29]
By integration, [START_REF] Markstein | Non-steady flame propagationOxford[END_REF] becomes:

dT dx = -2. Q λ .
w * dT [START_REF] Darrieus | Propagation d'une flamme. Essais des theories anormales de deflagration par developpement spontanee de la turbulence[END_REF] Then for [26]:

dT dx = m * . c (T -T ) λ ≈ m * . c (T -T ) λ [31]
This equation shows that the energy released by the flame contributes to preheat the reactants. Knowing that in the steady regime, λ. , = λ.

, :

m * . c (T -T ) λ = 2. Q λ .
w * dT [START_REF] Clanet | On the "tulip flame" phenomenon[END_REF] with m * = S . ρ :

S = 2. λ c J T -T with J = 1 a
. w * dT [START_REF] Eckhaus | Studies in nonlinear stability theory[END_REF] In many applications, it is sufficient to use a combustion model of order zero.

w * = Z . exp(-Ea RT ⁄ ) with N = 0 [START_REF] Batchelor | Introduction to Fluid Mechanics[END_REF] The hypothesis of high activation energy permits to express J:

J = Z a . R. T Ea .
exp -e RT [START_REF] Violeau | Les fondements de la turbulence[END_REF] so that finally

S = 2 • λ c ρ a . Z exp(-Ea RT ⁄ ) (T -T ) . R.
T Ea [START_REF] Werlé | Flow past a sphere near the critical Reynolds number[END_REF] A flame thickness characteristic scale is also defined. This quantity appears in the equation [26]: δ = λ ρc S [START_REF] Versteeg | An introduction to computational fluid dynamics: The finite volume method[END_REF] This reasoning can be extended to a number of Lewis Le ≠ 1 (but constant) and to a reaction order different to zero. For usual situations of flame propagation, N = 1 such as w * = Z . [A]. exp -.

where "A" represents the limiting reactant of the reaction (25):

S = 2. λ. Le. Z . R. T . exp - Ea RT ρ . Ea. c . (T -T ) / [38]
Maximum laminar combustion velocities of common gaseous fuels are on the order of 0.5 m/s (26).

This equation [START_REF] Tennekes | A first course in turbulence[END_REF], nowadays widely accepted, allowed to explain a number of physical behaviours of "premixed flames" including ignition, flammability limits, extinction, ... ( 27), ( 6), [START_REF] Proust | Formation-Inflammation-Combustion des atmosphères explosives (ATEX) et effets associés[END_REF],.... For the present purpose two main aspects need to be underlined:  The flame thickness is extremely small, typically a fraction of mm and represents mainly the preheating zone of the flame front. It can be shown that the chemical zone (zone II) is smaller than the preheat zone by a factor of Ea/RT b . Intuitively, the flame front would be affected by an external physical phenomenon only if it can act at the scale of the flame thickness and with a speed comparable to S lad ;

 Both parameters, flame thickness and laminar burning velocity, can be viewed as intrinsic and fundamental properties of many types of propagating flames in premixed media. Furthermore, those burning parameters are extremely sensitive to the flame temperature (T ) variations.

As explained later, if the flame front is curved towards the reactants, the heat flux from the burnt gases diverges and the heat diffused sidewise may be considered as a loss. This would reduce T and diminish S .

FLAME-ENVIRONMENT INTERACTIONS

This last point suggests that the flame front may be very sensitive to the disturbances provided by its environment.

The typical thickness of the flame is 10 -4 m (7). In practical situations (industrial safety), the scales of the geometry and of the flow disturbances are at least 2 orders of magnitude larger [START_REF] Proust | Formation-Inflammation-Combustion des atmosphères explosives (ATEX) et effets associés[END_REF] at least in industrial explosion configurations. For this reason, the flame front can be regarded as an interface (with negligible thickness) separating the unburnt from the burnt gases, propagating and separating two media of different densities. Using the conservation laws on both sides of the flame and using matching conditions at the interface, Markstein (28) formulated a "wave" model of the flame front. Then using the perturbation theory, he investigated the response of the flame front to various stimuli. The full development is rather tedious and may be found elsewhere [START_REF] Daubech | Contribution à l'étude de l'effet de l'hétérogénéité d'un pré mélange gazeux sur la propagation d'une flamme dans un tube clos[END_REF] but the major conclusions can be illustrated. For instance, it was noticed long ago that a flame front propagating in a perfectly quiescent and homogeneous mixture acquires a curved shape [START_REF] Maxworthy | Flame propagation in tubePhys[END_REF] and may even become "turbulent" (24) even though the mixture may be quiescent. Darrieus and Landau [START_REF] Darrieus | Propagation d'une flamme. Essais des theories anormales de deflagration par developpement spontanee de la turbulence[END_REF], (31) demonstrated that a propagating flame is intrinsically unstable.

A rough illustration can be obtained considering Figure 1-3 where an initially flat flame was suddenly disturbed to become wavy. Two phenomena are at work. First the flame front propagates normal to itself. Because of this, the bulges towards the reactants will grow whereas the valley will shrink leaving finally a succession of bulges. If this mechanism would play alone, the bulges would merge and the flame will flatten progressively. If fact, because of the expansion of the normal component of the flow velocity across the flame front, the flow-lines of the burnt product tend to converge and this creates a back pressure forcing those flow-lines to become parallel. This is possible only if the flame front steepens so as to increase the tangential component of the velocity at the expenses of the normal component. Finally, the bulges are continuously growing counteracting the merging mechanisms and the flame remains corrugated at least for those wavelengths where the growth rate exceeds the merging rate.

Because of this, this instability mechanism is due to the fact that the flame is a "propagating front separating two gases of different densities". In this model, the flame is unstable to all wavelength perturbations. It will amplify any disturbance.

In reality, experimental observations (confirmed later by Clanet and Searby [START_REF] Clanet | On the "tulip flame" phenomenon[END_REF])

show that the flame could be stable for at least small length scales. To account for this, Markstein [START_REF] Markstein | Non-steady flame propagationOxford[END_REF] extended the work of Darrieus-Landau [START_REF] Darrieus | Propagation d'une flamme. Essais des theories anormales de deflagration par developpement spontanee de la turbulence[END_REF], (31) and proposed to take account of the interaction between the flow and the flame structure. The underlying idea is that the burning velocity is affected by the "stretch" imparted by the flow ahead of the flame. To explain this, expression [START_REF] Tennekes | A first course in turbulence[END_REF] needs to be considered. This equation was extensively used over the last thirty years to study the behavior of disturbed laminar flames starting from the point that any perturbation inducing a change in the final flame temperature would strongly affect the burning velocity S (because of exp (-Ea RT ⁄

) in which the activation energy is assumed to be large). Within the present context, the flame stretching, imparted by the flow was invoked by Markstein [START_REF] Markstein | Non-steady flame propagationOxford[END_REF], and Eckhaus [START_REF] Eckhaus | Studies in nonlinear stability theory[END_REF] to be a potential cause for a variation of T . If the flame has a positive curvature (convex shape) with respect to the flow, the heat flux toward the reactants is not mono-directional, but has a tangential component. In the context of one-dimensional model proposed, this could be represented by a "loss" of heat for the convex parts of the flame (and a gain in the concave parts). As a first approximation, one might as well say that the flame temperature deficit associated would be proportional to this ratio:

T -T T -T = ∆T (T -T ) ≈ δ R [39]
where  T is the adiabatic temperature of the flame  R is the flow curvature radius

Using expression [START_REF] Tennekes | A first course in turbulence[END_REF], the change of S according to T can be found, using for instance a logarithmic derivative:

∆S S = - Ea 2. R. T . ∆T T = - β 2 . ∆T T [40]
where  ∆S is the variation of the burning velocities  β is the Zeldovitch number, equal to Ea R. T ⁄ .

Finally, using [START_REF] Reynolds | Effects of rotation on homogeneous turbulence[END_REF],

∆S S = - β 2 . ∆T T ≈ - β 2 .
δ R [START_REF] Candel | Recent advances in combustion modelling[END_REF] Within the limits of small disturbances, Makstein (Markstein, 1964 (28)) expressed the dependency of S as function of the curvature of the flame (relative to that of the flow) as:

S = S . 1 -L 1 R - 1 R
.

[42]

Where L , is the Markstein length proportional to the flame thickness: L = -. δ . The parameter measures the difference between the divergence of the flow and the curvature of the flame. If the former is larger than the latter, then the flame front is "stretched" by the flow and a part of the thermal energy transferred by the flame ahead in the reactant may be lost sideways. The flame velocity is reduced, and the disturbance will be dampened since the larger the disturbance the larger the reduction of the burning velocity. This would occur for wavelengths comparable to the flame thickness.

The most important implication is that the flame cannot remain passive when submitted to a perturbation coming from the flow. It will amplify or dampen it dependending on the wavelength of the perturbation relative to the flame thickness (or Markstein length). If the flame is considered as a wave, as suggested by Lewis, von Elbe, Markstein and many more, the flame cannot be simply considered as a passive interface, convected away and corrugated by the flow.

TURBULENCE AND COMBUSTION

FLOW AND TURBULENCE

Turbulence appears spontaneously in sheared region of flows, near the boundaries or in the wake of obstacles [START_REF] Batchelor | Introduction to Fluid Mechanics[END_REF]. Because the flame creates its own flow, due to the expansion of the burned gases, the influence of the turbulence on the propagation of the flame received a considerable attention in the past decades.

The phenomenology of turbulence is intrinsically linked to the nature of the Navier-Stokes equations which are on the same time dispersive (convective term) and dissipative (viscous term). Because of this it can be shown that the evolution of a small perturbation depends extremely on the initial conditions so that the resulting flow-field can be chaotic [START_REF] Violeau | Les fondements de la turbulence[END_REF].

But this can happen only when the viscous forces remain sufficiently low as compared to the inertial forces (convective term). The Reynolds number of a flow gives a measure of the relative importance of inertia forces (associated with convective effects) and viscous forces.

The initial perturbation originates in sheared zones of the flow and develops as an eddy. A typical example is the boundary near the wall in a duct or past a sphere as shown in Figure 1234. A small instability appears in the boundary flow, develops, and breaks down in a variety of flow structures having various scales. Reynolds then proposed to decompose, the instantaneous velocity u(t) into a steady mean value U and in a fluctuating component 𝑢′(𝑡) superimposed on it: 𝑢(𝑡) = 𝑈 + 𝑢′(𝑡). A turbulent flow can now be characterized in terms of the mean values of flow properties (𝑈, 𝑉, 𝑊, 𝑃 𝑒𝑡𝑐. ) and some statistical properties of their fluctuations (𝑢′, 𝑣′, 𝑤′, 𝑝′ 𝑒𝑡𝑐. ). These fluctuations are due to the passage at the measuring point of eddies of different sizes. The description of these various eddies and of their interactions was a central question that was investigated in depth during the 20 th century. It is accepted that the largest eddies may be as large as the broadest velocity gradient which is generally on the order of the flow itself. By rolling up inside those gradients, those eddies extract kinetic energy from the average flow. The characteristic velocity 𝜗 and characteristic length 𝑙 of the larger eddies ("integral scale") are of the same order as the velocity scale 𝑈 and length scale 𝐿 of the mean flow. Hence a "large eddy" Reynolds number 𝑅𝑒 = 𝜗. 𝑙/𝜈 (𝜈 is the kinematic viscosity of the fluid) formed by combining these eddy scales with the kinematic viscosity will be large in all turbulent flows, since it is not very different in magnitude from 𝑈𝐿/𝜈, which itself is large. This suggests that these large eddies are dominated by inertia effects and viscous effects are negligible. There are many mechanisms by which large eddies may give birth to smaller ones. One of them is called "vortex stretching". The presence of mean velocity gradients in sheared flows distorts the rotational turbulent eddies. Suitably aligned, eddies are stretched because one end is forced to move faster than the other. Because those large eddies are not damped by viscous effects, their angular momentum is conserved during vortex stretching. This causes the rotation rate to increase and the radius of their cross-sections to decrease. Thus, the process creates motions at smaller transverse length scales and also at smaller time scales. But in reality, other phenomena may intervene as well such as the natural instabilities (consider for instance the case of vortex rings) which break down large rotational structures. Swirling eddies convey in their periphery a strong velocity gradient in which smaller eddies may also appear. These effects cause the large eddies to break down into smaller and smaller structures ("turbulent cascade") down to an ultimate scale where the kinetic energy is dissipated into heat. The Reynolds number 𝑅𝑒 of these smallest eddies based on their characteristic velocity 𝜐 and characteristic length η is equal to 1, 𝑅𝑒 = 𝜐. 𝜂/𝜈 = 1, so the smallest scales present in a turbulent flow are those for which the inertia and viscous effects are of equal strength. These smallest scales are named the Kolmogorov microscales after the Russian scientist who carried out groundbreaking work on the structure of turbulence in the 1940's.

Finally, the kinetic energy conveyed by the fluctuations in a turbulent flow is spread over a wide range of "eddy" frequencies or wave numbers (𝜅 = 2𝜋𝑓/𝑈 where f is the frequency) as exemplified in Figure 123456, which gives the energy spectrum of the turbulence downstream of a grid.

The "spectral energy" 𝐸(𝜅) is shown as a function of the wavenumber 𝜅 = 2𝜋/𝜆, where λ is the wavelength of the eddies. The spectral energy 𝐸(𝜅) (units m 3 /s 2 ) is the kinetic energy per unit mass and per unit wavenumber of fluctuations around the wavenumber κ. The diagram shows that the energy content peaks at the low wavenumbers, so the larger eddies are the most energetic. They acquire their energy through strong interactions with the mean flow. The value of 𝐸(𝜅) rapidly decreases as the wavenumber increases, so the smallest eddies have the lowest energy content.

Figure 1-6: Energy spectrum of turbulence behind a grid [START_REF] Tennekes | A first course in turbulence[END_REF] Dimensional analysis can be used to obtain ratios of the length, time and velocity scales of the small and large eddies. This yields the following order of magnitude estimates of the ratios of small length, time and velocity scales 𝜂, 𝜏, 𝜐 and large length, time and velocity scales 𝑙, 𝑇, 𝜗 (38; 39):

𝐿𝑒𝑛𝑔𝑡ℎ -𝑠𝑐𝑎𝑙𝑒 𝑟𝑎𝑡𝑖𝑜 𝜂 𝑙 = 𝑅𝑒 ⁄ 𝑇𝑖𝑚𝑒 -𝑠𝑐𝑎𝑙𝑒 𝑟𝑎𝑡𝑖𝑜 𝜏 𝑇 = 𝑅𝑒 ⁄ 𝑉𝑒𝑙𝑜𝑐𝑖𝑡𝑦 -𝑠𝑐𝑎𝑙𝑒 𝑟𝑎𝑡𝑖𝑜 𝜐 𝜗 = 𝑅𝑒 ⁄
Typical values of 𝑅𝑒 might be 10 3 -10 6 , so the length, time and velocity scales associated with small dissipating eddies are much smaller than those of large, energetic eddies, and the difference -the so-called scale separation -increases as 𝑅𝑒 increases.

The behavior of the large eddies should be independent of the viscosity and should depend on the velocity scale 𝜗 and length scale 𝑙. Thus, on dimensional grounds we would expect that the spectral energy content of these eddies should behave as follows: 𝐸(𝜅) ∝ 𝜗 𝑙, where 𝜅 = 1/𝑙. Since the length scale 𝑙 is related to the length scale of turbulence producing processes -for example, boundary layer thickness 𝛿, obstacle width 𝐿, surface roughness height 𝑘 -the structure of the largest eddies is expected to be highly anisotropic (i.e. the fluctuations are different in different directions) and strongly affected by the flow boundary conditions.

Kolmogorov argued that the structure of the smallest eddies and, hence, their spectral energy 𝐸(𝜅 = 1/𝜂) should only depend on the rate of dissipation of turbulent energy ε (units m 2 /s 3 ) and the kinematic viscosity of the fluid 𝜈. Dimensional analysis yields the following proportionality relationship for the spectral energy: 𝐸(𝜅 = 1/𝜂) ∝ 𝜈 ⁄ 𝜀 ⁄ . Thus, for a given fluid with viscosity 𝜈, the spectral energy 𝐸(𝜅) of the smallest eddies only depends on the rate of energy dissipation and is not linked to other flow variables. The diffusive action of viscosity tends to smear out directionality at small scales. At high mean flow Reynolds numbers, the smallest eddies in a turbulent flow are, therefore, isotropic (i.e. they have the same properties in all directions).

Finally, Kolmogorov derived the universal spectral properties of eddies of intermediate size, which are sufficiently large so that their behaviour to be unaffected by viscous action (as the larger eddies), but sufficiently small that the details of their behaviour can be expressed as a function of the rate of energy dissipation ε (as the smallest eddies). The appropriate length scale for these eddies is 1/𝜅, and he found that the spectral energy of these eddies -the inertial subrange -satisfies the following relationship: 𝐸(𝜅) = 𝛼𝜅 ⁄ 𝜀 ⁄ . Measurements showed that the constant 𝛼 ≈ 1.5. Figure 1-6 includes a line with a slope of -5/3, indicating that, for the measurements shown, the scale separation is insufficient for a clear inertial subrange. Overlap between the large and small eddies is located somewhere around 𝜅 ≈ 1000.

Kolmogorov associated the dissipation rate of the turbulent energy ε (units m 2 /s 3 ) to the dimensions of the smallest eddies: 𝜀 = 𝜈 𝜂 .

Then, a fully developed turbulence cascade is totally defined when two independent parameters are given. Usually these are the turbulent kinetic energy 𝑘(≈ 𝜗 ) and the dissipation rate ε but it is sometimes 𝑘 and 𝑙 or 𝑘 and 𝜔 (𝜔 = ).

Finally, because of the existence of this eddies and "turbulent cascade" mass and momentum are diffused very efficiently. For example, a streak of dye which is introduced at a point in a turbulent flow will rapidly break up and be dispersed right across the flow. Such effective mixing gives rise to high values of diffusion coefficients for mass, momentum and heat.

TURBULENT PREMIXED FLAMES

Over the past fifty years, a major effort of understanding and modelling of the combustion of premixed gases in the presence of turbulence was accomplished (40), ( 41), ( 42)). Nevertheless, many aspects are still outstanding.

As demonstrated earlier, the combustion in a premixed flame is characterized by the laminar burning velocity and the flame thickness. The turbulence of the flowfield can be defined by the intensity of the velocity fluctuations and the integral length scale.

The turbulent combustion regimes are traditionally discussed using these parameters. In the extreme case where the size of the largest turbulent eddies is comparable or smaller than the thickness of the laminar flame 4 , some assume that the mass transfer by turbulent diffusion can occur and overlap with those of molecular diffusion. In the case of intense turbulence where u ′ is much larger than S , the combustion would rapidly have distributed throughout the explosive cloud: It is the "well-mixed reactor combustion regime". For values of u ′ comparable to S , the initial structure of the laminar flame would be disturbed by the vortices of small size that would affect the energy release rate (43) and would "thicken" the flame (Figure 1234567). Numerical simulation of these regimes was proposed introducing a heat production term in the Navier-Stokes equations in turbulent regime [START_REF] Borghi | Turbulent combustion modelling[END_REF]. But their real existence is a matter of debate.

Other situations seem more natural. When the size of the vortices is greater than the thickness of the flame, the basic structure of the turbulent flame is that of laminar flame that interacts with turbulent structures: it is the "flamelet" regime. It is admitted [START_REF] Peters | Laminar flamelet concepts in turbulent flames[END_REF], [START_REF] Williams | Combustion theory: 2nd edition[END_REF] that in most industrial situations including those relating to the explosion [START_REF] Proust | Les explosions de poussières, cours de 3 eme cycle[END_REF], the flame structure would fall into this regime.

Many studies were devoted to this regime ( 47), (45) ( 48), ( 43), [START_REF] Williams | Progress in knowledge of flamelet structure and extinction[END_REF]. One of the objectives is to link the "turbulent" burning velocity (speed with which the flame brush consumes the reactants) S t , to the fundamental burning and turbulence parameters. If it is assumed that the combustion is produced in local laminar flame fronts ("flamelets"), the main effect of the eddies is to "roll up" and corrugate the flame front, increasing its total area. In this case the ratio between S t /S lad would be equal to the flame area increase. This representation seems to be corroborated by experimental observations that show a certain similarity between the folding of a chemically inert interface between two immiscible fluids and the structure of turbulent flames ( 49), ( 50), ( 51), ( 52), [START_REF] Kwon | Flame surface properties of premixed flames in isotropic turbulence[END_REF]. But the propensity of a flame to wrinkle depends on its ability to accept an increase in its surface under the effect of the rolling up by the vortices by a mechanism called the "stretching" which involves the fundamental parameters of combustion that are S and η . Thus, rather recent studies ( 54), ( 55), [START_REF] Cambray | Length-scales of wrinkling of weakly forced unstable premixed flames[END_REF] show that S must be linked to parameters u', L, S and η . A simple dimensional analysis, supported by more phenomenological approaches [START_REF] Peters | Laminar flamelet concepts in turbulent flames[END_REF], [START_REF] Joulin | Cours dispensé à l'Ecole de Combustion[END_REF], suggests a relationship under the form:

S S = K u ′ S .
L η [START_REF] Williams | Combustion theory: 2nd edition[END_REF] where K, a and b are coefficients relatively independent of the mixing. But there may be as many correlations as authors. One reason is certainly a lack of understanding of the details of the underlying physics. Research in this area is ongoing and seems particularly active at present with the use of Direct Numerical Simulations and enhanced experimental technique (visualization).

1.2 COMPUTATIONAL FLUID DYNAMICS

BACKGROUND

CFD could be defined as the science of predicting fluid flow, heat transfer, mass transfer, chemical reactions, and related phenomena by solving mathematical equations governing these processes using a numerical method fitted to the nature of the latter equations. So there is a "physical" aspect related to the choice of the equations and a "mathematical" aspect pertaining to the way the equations are solved.

"Historically, methods were first developed to solve the linearized potential equations as those established for irrotational flows. One of the earliest types of calculations resembling modern CFD are those by the sense that these calculations used finite differences and divided the physical space in cells. Although they failed, these calculations, together with Richardson's book "Weather prediction by numerical process", and numerical meteorology. In fact, early CFD calculations during the 1940s using ENIAC (Electronic Numerical Integrator "Historically, methods were first developed to solve the linearized potential stablished for irrotational flows. One of the earliest types of calculations resembling modern CFD are those by Lewis Fry Richardson the sense that these calculations used finite differences and divided the physical space in cells. Although they failed, these calculations, together with Richardson's iction by numerical process", set the basis for modern CFD and numerical meteorology. In fact, early CFD calculations during the 1940s (Electronic Numerical Integrator Analyser and Co methods close to those in Richardson's 1922 book [START_REF] Awate | CFD -Introduction and Application[END_REF]."

Given the enormous calculation requirement of its model, Richardson proposed a technical solution that he called the "forecast-factory". The "factory" would gather 64,000 people in a stadium. Each one, using a mechanical calcula ), would perform a part of the calculation. A leader in the center, using colored signal lights and telegraph communication, would coordinate the forecast. During the 1960s the theoretical division at Los Alamos [START_REF] Harlow | Fluid dynamics in Group T-3 Los Alamos National Laboratory[END_REF] numerical methods that are still in use today, such as the method of Particle Cell (PIC), the method of Marker-and-Cell (MAC), the Vorticity Methods, the Arbitrary Lagrangian-Eulerian (ALE) and the k-epsilon s, a group working with D. Brian Spalding the Parabolic flow codes (GENMIX), the Vorticity function based codes, the SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm, the form of the k-ε equations that are used today, the , the 'Eddy break-up' and the 'presumed pdf' combustion basis of the actual Computational Fluid Dynamics was Page 42 sur 225 "Historically, methods were first developed to solve the linearized potential stablished for irrotational flows. One of the earliest types of Lewis Fry Richardson [START_REF] Hunt | Lewis Fry Richardson and his contributions to mathematics, meteorology, and models of conflict[END_REF], in the sense that these calculations used finite differences and divided the physical space in cells. Although they failed, these calculations, together with Richardson's set the basis for modern CFD and numerical meteorology. In fact, early CFD calculations during the 1940s (Electronic Numerical Integrator Analyser and Computer), used

Given the enormous calculation requirement of its model, Richardson [START_REF] Hunt | Lewis Fry Richardson and his contributions to mathematics, meteorology, and models of conflict[END_REF] factory". The "factory" would gather 64,000 people in a stadium. Each one, using a mechanical calculator ), would perform a part of the calculation. A leader in the center, using inate the forecast.

: Examples of mechanical calculators dimensional (2D) methods, using Joukowski transformations from a flow were developed in the 1930s obtains a solution for flow around a cylinder, in 1953 by using a mechanical desk calculator, working 20 hours per week for 18 months, requiring: contributed many numerical methods that are still in use today, such as the method of Particle-In Cell (MAC), the Vorticity-Stream function epsilon turbulence a group working with D. Brian Spalding [START_REF]CFD: Past, Present and Future by Brian Spalding Lecture at the Sixteenth Leontiev School-Seminar[END_REF] at Imperial the Vorticity-Stream-Implicit Method for Pressure Linked ε equations that are used today, the up' and the 'presumed pdf' combustion basis of the actual Computational Fluid Dynamics was laid.

At the beginning, CFD was performed in an academic research environment and in-house made codes were issued. Most of the commercial CFD software that are available today were issued in this way:  Fluent (UK and US)  CFX (UK and Canada)

 CFD++ (US)  Star CD (UK)  Polyflow (Belgium)  Flow 3d (US)  FLACS (Norway)  SCRYU (Japan)  Siemens PLM (Germany)  CAST3M (France)  …
In present Computational Fluid Dynamics, equations representing fluid mechanics are solved using a "numerical technique". This means that they are "discretized" in time and space and then solved step by step. Usually an elementary volume corresponding to the space step is called "cell". With this technique, an exact solution would only be possible if the time step and cell size would be smaller than the smallest physical scale. There are at least two of them: the chemical scale of the combustion in the flame, typically on the order of 10 -5 to 10 -4 m and the Kolmogorov scale of the turbulence (some mms). Obviously, at the industrial scale (tens of meters), the number of cells would be very large (10 15 at least) and incur unacceptable calculation costs. To illustrate this difficulty, 24 hours of calculation are required to solve a rather simple explosion case under a robust and fast RANS approximation over only 10 6 cells using a parallel algorithm working on five processors on a powerful PC.

Following, "direct or full numerical simulation", although feasible, is limited to very small geometries. At larger scale, small scale phenomena, typically those smaller than the cell size, are "smoothed" out and replaced by physical correlations to represent the interaction between the small scale and the large-scale mechanics. This is the "modelling" aspect of CFD. Today, usually three CFD techniques [START_REF] Eyssartier | Mise en place de l'interface de Pre-processing pour le code de calcul[END_REF] are employed:  Direct Numerical Simulation (DNS) is the most expensive in terms of calculation power. It is the simplest to design because it doesn't require any "modelling" but the discretization steps need to be very small. DNS method remains mostly restricted to research and academic calculations over simple and small geometries. Nevertheless, fruitful and meaningful "numerical" experiments can be performed with this technique and significant progresses were made in the last decades about the understanding of the turbulence and flame dynamics ;

 LES (Large Eddies Simulation) techniques appeared roughly at the turn of the century in an attempt to represent more faithfully the dynamics of the large scales of the turbulence. To do this, a low pass spatial filter is applied to remove all scales smaller than a specified length which corresponds approximately to the cell size so that only the largest scales of the flow are calculated. The incidence of filtered scales on the large-scale flow is "modeled" to introduce their influence of the transport of species and momentum. This last point is an active subject for research. Clearly, this method seems well suited to flows with large-scale unsteadiness. In practice however, this method is demanding in terms of computing resources because the calculations should be performed in three dimensions and because about 80% of the turbulent kinetic cascade should be resolved [START_REF] Pope | Turbulent Flows[END_REF] meaning typically a maximum cell size on the order of cms as far as industrial processes. Because of this, applying LES to large scale and complex geometries remains a challenge ;

 Still today, the traditional unsteady Reynolds-Averaged Navier-Stokes (uRANS) technique is largely used in industrial applications. One of the key of the success of uRANS techniques is that the physical representation of the flow corresponds to the current understanding: in particular, the turbulence of the flow is calculated as a random fluctuating velocity (average value = 0) superimposed on the mean flow following the Reynolds representation. The evolution of the turbulence in time and space is only represented by the leading parameters of the turbulent cascade (the "moments" of the random distribution) and interacts as a source or sink with the mean flow. Since only the mean flow is fully calculated, the discretization of space and time may be rather coarse and, in many situations, a two-dimension cylindrical geometry for example can mimic a three-dimensional problem which reduces even further the computational demand. Consequently, very large geometries can be considered.

CFD calculations applied to explosions were attempted rather early. Twodimensional experiments were simulated while developing advanced combustion modelling [START_REF] Barsanti | Modelling of confined turbulent explosions[END_REF], and there was a significant effort in CFD optimisation using adaptive mesh refinement [START_REF] Fairweather | A Two-Dimensional, Axisymmetric Computational Fluid Dynamic Code Capable of Predicting the Near-Field Structure of High Pressure Releases of Supercritical and Multi-Phase Carbon Dioxide[END_REF]. Recent algorithmic improvements as well as increases in computer power increased the feasibility of CFD calculations of explosion phenomena in more realistic three-dimensional configurations. Such methods were extended in the 1980's to safety applications, including explosions, with the objective of becoming a component in the range of available safety-case analysis tools [START_REF] Ledin | A review of the state-of-the-art in gas explosion modelling[END_REF].

However, two significant limitations (66) need to be considered. Even with modern computers, experience shows that each elementary volume used by a CFD code requires abound 10 3 bytes of computer memory. So that, the maximum number of cells available to represent a geometry on a powerful desktop PC is on the order of 106 . In three dimensions, this would allow approximately 100 computing cells in each co-ordinate direction, resulting in a typical cell dimension between 0.1 and 1.0 m for a typical process plant. The first consequence is that many of the objects/features within the process plant that may be important for flame acceleration (obstacles) could be smaller than this and a "trick" will be needed to represent their effect. For instance, the Porosity /Distributed Resistance (PDR) concept was developed to allow this to introduce the influence of the "subgrid" size obstacles. A second consequence is that, because the flame thickness is orders of magnitude smaller than the cell size, the physics of the combustion in the flame need also to be artificially introduced. A similarly reasoning is true also for the representation of the turbulence which usually cannot be completely resolved at the scale of the cell size. The methods used to introduce the "subgrid" phenomena are usually called "models".

In summary, CFD codes devoted to large scale explosion phenomena cannot simulate the full reality and the way the "models" of the codes (representation of the obstacles, of the turbulence, of the flame…) are chosen/implemented may have a strong impact on the simulation results. Furthermore, even admitting the computers would in the future be powerful enough to avoid using such "models", it would be very difficult to introduce all the required geometrical details. Some simplifications will be made by the user. In both case, the reality may be tweaked.

Even though some applications of LES to explosion modelling appeared recently mostly for small geometries [START_REF] Ibrahim | The effects of obstructions on overpressure resulting from premixed[END_REF], [START_REF] Molkov | LES modelling of an unconfined large-scale[END_REF] In Computational Fluid Dynamics (CFD), some forms of the conservation laws of mass, momentum, species and energy are solved. In the specific case of gaseous atmospheres, the compressibility of the medium needs also to be accounted for (Veynante, 1999 (70)). On many aspects, CFD codes applicable to explosions are amongst the most complete CFD engineering tools since they incorporate compressibility effects and chemical reactions. The conservation laws read:

 Mass balance :

∂ ∂t + ∂u ∂x = 0 [44]
u is the velocity of the flow,  is the density, t the time and x j are the space coordinate (j=1, 2,3). The Einstein notation is used where:

 =  +  + 
(the summation is applied on index "j").

 Momentum balance :

∂u ∂t + ∂u u ∂x = - ∂P ∂x + ∂𝒯 ∂x + F [45]
𝒯 is the viscous force tensor and F a "body force" (most often gravity). In practical situations, the fluids are assumed to be Newtonian, i.e. the viscous tensor is given by the Newton law:

𝒯 = μ ∂u ∂x + ∂u ∂x - 2 3
μ δ ∂u ∂x [START_REF] Proust | Les explosions de poussières, cours de 3 eme cycle[END_REF] where the laminar viscosity μ , depending on the fluid properties, is introduced. δ is the Kronecker symbol which is 1 if i = j and 0 if i ≠ j.

 Species balance (N species with k=1,…,N):

∂Y ∂t + ∂u Y ∂x = ∂𝒥 ∂x + ω̇ [47]
Y k is the mass fraction of species k, 𝒥 is the molecular diffusive flux of species k and ω̇ is the mass reaction rate of species k per unit volume. Species molecular diffusivities are generally described using the Fick law, assuming the major species with constant pressure and identical body forces.

𝒥 = μ Sc

∂Y ∂x [START_REF] Sivashinsky | Cascade renormalisation theory of turbulent flame speeds[END_REF] Sc is the Schmidt number of species k ,defined as: [START_REF] North | The fractal nature of premixed turbulent flames[END_REF] D is the molecular diffusivity of species k relatively to the major species. More complex expression may be used to describe multi-species molecular diffusion. Soret effect (species diffusion under temperature gradients) and molecular transports due to pressure gradient are usually neglected.

Sc = μ ρD
 Energy (total enthalpy) balance

∂ρh ∂t + ∂u h ∂x = ∂P ∂t + ∂ ∂x 𝒥 + u 𝒯 + u F [50]
where:

h = Y h Sc = μ ρD [51]
The specific enthalpy h of the species k is:

h = h + ∫ C , dT , [ 52
]
where the standard specific enthalpy of formation for species k and its specific heat capacity at constant pressure are respectively h and C , . Parameters u 𝒯 and u F denote respectively the power produced by the viscous and body forces. The diffusion term in the enthalpy equation is described according to the Fourier law: The Prandtl number Pr compares the diffusive transport of momentum (viscous forces) and thermal diffusivity. In equation [START_REF] Kwon | Flame surface properties of premixed flames in isotropic turbulence[END_REF], radiative heat transfert and Duffour effect (enthalpy diffusion under mass fraction gradient) are neglected because they are of a smaller order of magnitude than the effects described by the Fourier law. The Prandtl number is written as a function of the thermal diffusivity λ and the constant pressure specific heat C :

Pr = μ C λ 𝒥 = - μ Pr ∂h ∂x + Pr Sc -1 h ∂Y ∂x [54]
Then, the Lewis number Le of the species k, comparing thermal and mass diffusivity is introduced:

Le = Sc Pr = λ ρC D 𝒥 = - μ Pr ∂h ∂x + Pr Sc -1 h ∂Y ∂x [55]
Under the assumption of unity Lewis number, the enthalpy diffusive flux (equation [START_REF] Kwon | Flame surface properties of premixed flames in isotropic turbulence[END_REF]) is simplified and mass fraction and enthalpy balance equations are formally identical. This assumption is generally made to simplify turbulent flame modelling. Nevertheless, thermo-diffusive instabilities occur in premixed system when the Lewis number is lower than unity (for example for hydrogen). One direct consequence of these instabilities is an increase of the premixed flame area and of the global reaction rate. Because of compressibility, the heat release by the combustion is responsible for large expansions and movements. Consequently, an equation of state is needed to link the temperature to volume/pressure. Usually, the perfect gas law is used: [START_REF] Cambray | Length-scales of wrinkling of weakly forced unstable premixed flames[END_REF] where r = R W [START_REF] Hunt | Lewis Fry Richardson and his contributions to mathematics, meteorology, and models of conflict[END_REF] and W m is the mean molecular weight of the mixture given by:

P = ρrT
1 W = Y W [58] 1.2.

RANS AND URANS FORMALISM

The RANS technique is based on the Reynolds description of a turbulent flow-field (75), each quantity 𝑓 (especially the flow velocity) is split into a mean quantity 𝑓 ̅ and a random deviation from the mean noted f :

f = f ̅ + f with f = 0 [59]
𝑓 ̅ is a statistic average which corresponds to the quantity registered at the same position 𝑥 ⃗ and at the same time t at each experiment N realized on the same flow. It is given by:

𝑓 ̅ (𝑥 ⃗, 𝑡) = lim → 1 𝑁 𝑓 ( ) (𝑥 ⃗, 𝑡) [60]
This average satisfies to the following properties:

𝑓 + 𝑔 = 𝑓 ̅ + 𝑔̅ 𝛼𝑓 = 𝛼𝑓 ̅ 𝑤𝑖𝑡ℎ 𝛼 = 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 𝑓 ̿ = 𝑓 ̅ 𝑓 ̅ 𝑔 = 𝑓 ̅ 𝑔̅ 𝜕𝑓 𝜕𝑥 = 𝜕𝑓 ̅ 𝜕𝑥 𝜕𝑓 𝜕𝑡 = 𝜕𝑓 ̅ 𝜕𝑡
One of the difficulty of this classical Reynolds averaging technique, widely used in non-reacting fluid mechanics, is that several cross correlated terms such as 𝑢 𝑓 appear which are unknown and thus must be "modelled". A second difficulty is that in combustion problems, fluctuations of density are observed because of the thermal heat release, which generates many cross correlated terms. For instance, averaging the mass balance equation leads to:

𝜕 𝜕𝑡 + 𝜕 𝜕𝑥 ρu = 𝜕 𝜕𝑡 + 𝜕 𝜕𝑥  . 𝑢 +  ′ 𝑢 ′ = 0 [61]
where a velocity -density fluctuation correlation  𝑢 appears (similar 𝑢 𝑓 terms appear in the other conservation equations).

To limit the number of such cross correlations, Favre introduced another averaging technique (the "URANS" technique) in which all the averages are "mass weighed" f such that:

f = f + f where f = [62]
It can be shown that the Favre average of the fluctuating part is zero. With this method the cross correlated terms ρf are written ρ f = ρf and taking the example of the mass balance equation, it comes out:

∂ ∂t + ∂ u ∂x = 0 [63]
which is formally identical to the Reynolds averaged continuity equation for a constant density flow but looks simpler. This result is true for any balance equation (momentum, energy, mass fractions...). Finally, using Favre averages: With the equation of state:


P = ρ rT [67]
Those equations are systematically computed in current URANS codes devoted to compressible flows such as those considered in this work. Despite the simplifications brought by the Favre averaging, cross correlated terms nevertheless remain and need to be "modelled" as described in following paragraphs. This "modelling" work is often called "closure problems". Note that Favre averaging is mostly a mathematical simplification which physical significance may be questionable. By definition, Reynolds stresses are ρu u = ρ u u ≈ ρu u . Such parameter "measures" the forces applied on one flow line (along u i ) onto its neighbour by the transverse momentum transfer (along u j ) due to turbulence motions. Many CFD-RANS codes are based on the assumption that the Reynolds stresses can be linked to the average values of the velocities. Boussinesq observed that the turbulence (u k ') first diffuses and second is dissipated just like the molecular viscosity does in a laminar flow and proposed a similar formulation:

ρ u u = -μ ∂u ∂x + ∂u ∂x - 2 3 δ ∂u ∂x + 2 3
ρ kδ [START_REF] Wilkening | An adaptive 3-D CFD solver for explosion modeling on large scales[END_REF] where 𝜇 is a turbulent dynamic viscosity and 𝛿 is the Kronecker symbol. The last term of the equation [START_REF] Wilkening | An adaptive 3-D CFD solver for explosion modeling on large scales[END_REF] is included to ensure that the trace of the Reynolds stress tensor is equal to -2ρ k , as it should be. Parameter k is known as the "turbulent kinetic energy":

k = 1 2
u " u " + v " v " + w " w " .

[69]

Where u, v and w are respectively the velocity component following x, y and z space directions.

To our knowledge, the domain of validity of the Boussinesq model is not clearly known (i.e. demonstrated). It can be shown that it is valid between the viscous sub-layer at the wall and the transitional zone in the core flow where strictly speaking µ t =µ l . In practise however, its validity is really to be questioned since in many practical situations, very significant Reynolds stresses were measured in locations where the average velocity gradients are zero [START_REF] Schmitt | About Boussinesq's turbulent viscosity hypothesis:historical remarks and a direct evaluation of its validity[END_REF]. One of the reasons for this disagreement might be that the Boussinesq model mixes both the turbulence mechanism generation (in a shear layer as effectively represented by gradient approximation) and its transport throughout the flow which might depend significantly on the natural lifting effect of the swirls generated in the shear layers. At best, the Boussinesq model may be considered as a global average representation of the mechanical balance/interactions between the mean flow and the turbulence. [71]

A
After calculating of the averages, the following scalar equation is obtained:

∂(ρ k) ∂t + ∂ ρ u k ∂x = ρ u v ∂u ∂x -τ ∂u ∂x + ∂ ∂x τ u -ρu 1 2 u u -P u [72] -u ∂P ∂x + P ∂u ∂x .
In this equation, ρ u v represents the production of turbulent energy by the action of the Reynolds stresses ρ u v expressed by equation [START_REF] Wilkening | An adaptive 3-D CFD solver for explosion modeling on large scales[END_REF]. The second term τ is the dissipation by the molecular diffusivity. By definition, it is equal to 𝜌̅ 𝜀 where ε is the turbulence dissipation and is always positive. The third term τ u -ρu u u -P u is homogeneous to a diffusion process. It is easy to

show that τ u ≈ μ represent the molecular diffusion term whereas ρu u u can be viewed as the turbulent transport of the turbulent energy and somewhat empirically expressed as ρu u u ≈ μ . The term P u , called pressure diffusion, is not really modeled and is currently associated to the turbulent transport of k and modeled as follow: ρu u u + P u ~ where σ k is a constant of order 1.

u is the work done by the mean pressure gradient against the turbulent fluctuations and P , represents the work done by the pressure fluctuations to expand the turbulent vortices. These terms cause some serious modeling problems and are still a research issue. But since it can be shown that they are equal to zero when the flow is incompressible, in practice, they are often neglected even when the flow is compressible. So, the final form of the k equation reads:

∂(ρ k) ∂t + ∂ ρ u k ∂x = μ ∂u ∂x + ∂u ∂x - 2 3 δ ∂u ∂x + 2 3 ρ k ∂u ∂x -ρ ε + ∂ ∂x μ + μ σ
∂k ∂x [START_REF] Ibrahim | The effects of obstructions on overpressure resulting from premixed[END_REF] The problem may be "closed" if relationships are established between ε, µ t and k for instance. The turbulent viscosity cannot be a property of the fluid (unlike µ l ) but should rather depend on the turbulence characteristics. As shown in the section 1.1.3, the turbulence characteristics may be described by two independent parameters: either u' and l, k and ε or any mixture of these. A closure is partly possible on the basis of the description of the turbulence. The turbulent viscosity is proportional to the turbulence intensity times (k 1/2 ) and to the "mixing" length (l) which measures the scale (lifetime) of the turbulence so that μ = ρ k l . Based on the analysis proposed by Kolmogorov ε = C / . But one independent parameter remains.

The favored approach is to introduce another balance equation to calculate the "independent" parameter. In the famous "k-epsilon" models, the latter is often ε. As for the turbulent kinetic energy, an exact equation for Ɛ is obtained from the basic equations and some modeling is proposed for the various terms. The full equation for Ɛ can be obtained at least for the case of incompressible flows applying the following operator: Consider the right-hand side of this equation. The first term may be viewed as an accumulation of turbulence dissipation in the unit volume due to the deformation of the flow. The second term represents the "production" of the dissipation due to the vorticity of the mean flow. The third one represents the change in the intensity of the turbulence dissipation due to vortex stretching. These three terms are thought to be the source of the "turbulent dissipation" in the unit volume. The fourth term represents the final destruction of the turbulence (including the dissipation!) by the molecular viscosity. The last term of equation is the classical diffusion contribution. These terms are currently not precisely modelled but may be analysed in orders of magnitudes. The final expression reads:

∂(ρ ε) ∂t + ∂ ρ u ε ∂x = C ρ ε k μ ∂u ∂x + ∂u ∂x - 2 3 δ ∂u ∂x + 2 3 ρ k ∂u ∂x -C ρ ε k + ∂ ∂x μ + μ σ
∂ε ∂x [START_REF] Schmitt | About Boussinesq's turbulent viscosity hypothesis:historical remarks and a direct evaluation of its validity[END_REF] With constants fitted on experiments (grid turbulence):

C = 0.09 ; σ = 1.0 ; σ = 1.3 ; C = 1.44 ; C = 1.
92 [START_REF] Cinnela | Modélisation et calcul d'écoulements compressibles turbulents[END_REF] Knowing that: μ = ρ k l and that ε = C / , it follows that μ = ρ C (since l =

).

The molecular viscosity is in practice neglected because smaller by orders of magnitudes as compared to the turbulent viscosity. The merit of the k-Ɛ model is to incorporate a very significant part of the current knowledge about turbulence.

There are many variants of this ground standard k-epsilon model (78) (Realizable k-epsilon model, RNG k-epsilon model, Low Reynolds number k-epsilon model…) but the standard k-epsilon model remains the basic choice in CFD industrial explosion simulation.

The standard k-epsilon model has however many shortcomings, the first being the empirical character of the Boussinesq approximation as explained earlier. Another one is certainly its relative inability in representing highly rotating flows (wake flows). A third one is its limitation to significantly turbulent flows because the various approximations above may hold only when u' is much smaller than the mean velocity and when the turbulence cascade is sufficiently developed for the homogeneity of the turbulence to cover a significant part of the turbulence spectrum so that the turbulent diffusion could be reasonably homogeneous (μ t is isotropic). Note that, outside the boundary layers where the two above conditions are not met, the Reynolds numbers can be very large within the scope of large scale explosions. Assuming for instance a flow at a few m/s and a typical size of 1-10 m, the Reynolds number is 10 6 so that the main underlying assumptions of the k-epsilon model could be satisfied.

Nevertheless, the question of the boundary layers near the solid walls is particularly tricky and needs a special treatment. As it stands, the "k-epsilon model" would produce far too much turbulence in the boundary layers. To show this point, it is considered that in the boundary layer the viscous forces dominate the flow so that the turbulent viscosity might be smaller than the molecular viscosity. A comparison betwen μ t and μ in the boundary layer is to be performed.

To do this, the expression of the former as function of k and ε should be used so as the definition of ε. If a plane boundary layer of thickness y+ is considered with k* being the kinetic turbulent energy at y+, the turbulence velocity k' in the boundary layer should drop to zero at the wall so that at a first approximation (first order Taylor expansion of the flow velocity as function of the distance y to the wall):

k′ = k * . y y [START_REF] Cable | An Evaluation of Turbulence Models for the Numerical Study of Forced and Natural Convective Flow in Atria[END_REF] Under the same assumption:

ρ . ε ≈ μ ∂u ∂x ∂u ∂x ≈ 2. μ . k * y [79]
Finally:

μ μ = ρ C 2. μ . k * . y .
𝑦 𝑦 [START_REF] Shepherd | Turbulent transport in a confined premixed flame[END_REF] The first ratio of the right-hand side is typically on the order of 10 9 . Typical boundary thicknesses are on the order of a tenth of the size of the obstacle/confinement so, at industrial scales, about 0.1 m. Assuming a typical velocity fluctuation of 1 m/s, it turns out that middle in the boundary layer (y=y+/2), the ratio μ t /μ l would amount 10 6 which is much too large (should be close to 0). Because of this, far too much turbulence would be produced in the boundary layer, leading also to a wrong estimate of the local mean velocity gradient.

To solve this difficulty, the "classical" solution is to apply an analytical "boundary law" near the wall, between 0 and y+, and to match it with the "k-epsilon" solution applied everywhere else in the bulk flow. This solution is numerically appealing (low cost) but the performances of the simulations are sensitive to the location of the matching points, y+, which is readily understandable from the preceding discussion because of the so large dependency of μ t /μ on y+ at the limit of the boundary layer (consider for instance equation [START_REF] Shepherd | Turbulent transport in a confined premixed flame[END_REF] with y=y+).

An alternative consists in modifying the standard k-epsilon model next to the wall. In that zone, it is sufficient to set k (and the mean flow velocities) to zero while keeping its "standard" value in the bulk flow. Under these conditions equation [79], giving the evolution of the turbulent kinetic energy in the boundary layer, reduces to = (because the turbulent viscosity, being proportional to k is zero) which is physically exact. Equation [START_REF] Ibrahim | The effects of obstructions on overpressure resulting from premixed[END_REF] would result in a sort of indetermination. To solve this, we need to go back again to the definition of the turbulent dissipation rate ε = ν which may be written ε = 2ν

. To introduce these models, Jones and Launder defined a transformed variable ε * vanishing on the solid boundaries:

ε * = ε -2ν ∂k ∂y [START_REF]ANSYS Fluent User's Guide[END_REF] To introduce the modified dissipation ε * , the conservation equation for k must be changed which is done by including an additional term of dissipation which is fully active in the bulk flow,

D = ρ ε -ρ ε * .
[82]

∂(ρ k) ∂t + ∂ ρ u k ∂x = ρ R ∂u ∂x -ρ ε * + ∂ ∂x μ + μ σ ∂k ∂x -D [83]
In a similar way an additional term E is required in the equation for the modified turbulent dissipation rate:

∂(ρ ε * ) ∂t + ∂ ρ u ε * ∂x = C ρ ε * k R ∂u ∂x -C f ρ ε * k + ∂ ∂x μ + μ σ ∂ε * ∂x + E [84] with μ = C f ρ k ε [85]
and

E = 2 μμ ρ ∂ u ∂x ∂x ∂ u ∂x ∂x .
[86]

In addition to the terms D and E (which are called the low Reynolds numbers, while they intervene in the near boundary zone, where the local turbulent Reynolds number is low) the near boundary version of this model contains some damping terms, f and f : the function f that appears in the term of "dissipation of the dissipation" considers the effects of the low local Reynolds number near to the boundary, which contributes to the growth of Ɛ in this zone; the function f which intervenes in the calculation of μ is a damping function permitting to quickly decrease the turbulent viscosity in the zone of near boundary. These functions are written:

f = 1 -0.3 exp(-R ) and f = exp -2.5
1 + R 50 [START_REF] Patankar | A calculation procedure for the transient and steady-state behaviour of shell-and -tube heat exchangers[END_REF] where R = k (νε) ⁄ is the "turbulent Reynolds" number.

In many practical situations, the buoyancy effects need to be accounted for. To do this, specific source terms are added into the momentum equation along the z coordinate so that the modified balance equations read:

 Mass equation ∂ρ ∂t + ∂ρ u ∂x + ∂ρ v ∂y + ∂ρ v ∂z = 0 [88]  Momentum equation in x direction (horizontal) ∂ρ u ∂t + ∂ρ u 2 ∂x + ∂ρ u v ∂y + ∂ρ u w ∂z = - ∂P ∂x + μ l + μ t ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 + ∂ 2 u ∂z 2 [89]
 Momentum equation in y direction (horizontal)

∂ρ v ∂t + ∂ρ u v ∂x + ∂ρ v 2 ∂y + ∂ρ v w ∂z = - ∂P ∂y + 𝜇 𝑙 + 𝜇 𝑡 ∂ 2 v ∂x 2 + ∂ 2 v ∂y 2 + ∂ 2 v ∂z 2 [90]
 Momentum equation in z direction (vertical)

∂ρ w ∂t + ∂ρ u w ∂x + ∂ρ v w ∂y + ∂ρ w 2 ∂z = - ∂P ∂z + μ l + μ t ∂ 2 w ∂x 2 + ∂ 2 w ∂y 2 + ∂ 2 w ∂z 2 + ρg 𝜌 -𝜌 𝑎 ρ [91]
 Species equation

∂ρ 𝑌 ∂t + ∂ρ u 𝑌 ∂x + ∂ρ v 𝑌 ∂y + ∂ρ w 𝑌 ∂z = ∂ ∂x ρ D + μ Sc ∂𝑌 ∂x + ∂ ∂y ρ D + μ Sc ∂𝑌 ∂y + ∂ ∂z ρ D + μ Sc ∂𝑌 ∂z [92]
The turbulence model is also affected in the following way :

∂(ρ k) ∂t + ∂ρ u k ∂x + ∂ρ v k ∂y + ∂ρ w k ∂z = ∂ ∂x μ + μ σ ∂k ∂x + ∂ ∂y μ + μ σ ∂k ∂y + ∂ ∂z μ + μ σ ∂k ∂z +P + P + ρ ε [93] ∂ρ ε ∂t + ∂ρ u ε ∂x + ∂ρ v ε ∂y + ∂ρ w ε ∂z = ∂ ∂x μ + μ σ ∂ε ∂x + ∂ ∂y μ + μ σ ∂ε ∂y + ∂ ∂z μ + μ σ ∂ε ∂z +C ε ε k (P + C ε P ) -C ε ρ Ɛ k [94]
Where P incorporates the influence of the turbulence driven effects :

P = - μ ρSc g∇ρ [95]
As said previously, many versions of the k -Ɛ model exist. The main difference among these models is the way to take into account of the boundary zone. The low Reynolds terms is a solution to this difficulty (k = 0, ε * = 0). To enable these models to give correct results, this region, which is extremely sensitive for high Reynolds flows, must be sufficiently discretized which implies the use of very fine cells near the boundary. But doing so, numerical robustness problems may occur. In practice, in several "free boundary" flows like jets, turbulence is created in the mean velocity gradients which correspond to the Boussinescq approximation and implies high Reynolds flows. Because of such features, the RANS+k-epsilon modeling is expected to perform reasonably in this zone. The main remaining question is the accuracy of the model since the turbulence intensity may not be very small as compared to the mean velocity.

The interaction of the flow with obstructions is more difficult to handle for standard k-epsilon model. The Boussinecq approximation may be inadequate in stagnation zones (of the obstacles) facing the flow, where the streamlines are deflected sidewise inducing a large velocity gradient but without turbulence generation. In addition, the k-epsilon model may not perform well around the obstacles because of the curvature of the flow which contribution in the production of turbulence is neglected as explained above. This is one of the reason why the Porosity Distributed Resistance concept was developed.

SPECIES (𝐮 𝐘 𝐤 ) AND ENTHALPY (𝐮 𝐡 𝐬 ) TURBULENT FLUXES

Following the same reasoning as above, these fluxes are closed using the same Boussinescq "gradient assumption" [START_REF] Veynante | Turbulent Combustion Modeling[END_REF]:

ρ u Y = - μ Sc
∂Y ∂x [START_REF] Kolms | Computing unstructured finite element grids on free-form CAD surfaces[END_REF] where 𝜇 is the turbulent viscosity, estimated from the turbulence model, and 𝑆𝑐 a turbulent Schmidt number for species k. This assumption is at least approximate, for similar reasons than for the "Boussinesq" assumption for the momentum equation.

As for the modelling of the Reynolds stresses, theory and experiments (79), [START_REF] Shepherd | Turbulent transport in a confined premixed flame[END_REF] have shown that this gradient assumption is wrong in some turbulent premixed flames: counter-gradient turbulent transport (i.e. in an opposite direction compared to the one predicted from Eq. [START_REF] Kolms | Computing unstructured finite element grids on free-form CAD surfaces[END_REF]) can be observed in weak turbulence flames.

The molecular terms can often be neglected against turbulent transport, assuming a sufficiently large turbulence level (large Reynolds number limit). Whenever required they are modeled very simply. For example, species molecular diffusion fluxes are generally modeled as:

𝒥 = -ρD ∂Y ∂x = -ρ D ∂Y ∂x [97]
where D is a "mean" species molecular diffusion coefficient. In a similar way, the molecular heat diffusion flux in the enthalpy equation is written:

λ ∂T ∂x = λ
∂T ∂x [START_REF] Figueira Da Silva | Unstructured adaptive grid flow simulations of inert and reactive gas mixtures[END_REF] where 𝜆 ̅ denotes a mean thermal diffusivity.

THE DIFFICULTY WITH NON-TURBULENT FLOWS IN URANS

The flow might be laminar at least initially. URANS formalism implies that for any velocity, some turbulence is produced. To some extent, URANS "forces" turbulence to appear although in the real life, the flow could remain laminar (as for instance in a pipe when Re<2000).

A further difficulty is that k cannot be equal to zero in the bulk flow otherwise the calculation of μ would lead an undetermined value. Non-zero initial values need to be given (especially for k), knowing that whatever that initial value, k and Ɛ will converge to the solution in line with the bulk flow characteristics. Jones and Launder (Table 1) proposed to choose initial values not too far from the converged solutions. Starting from these initial values, the URANS calculations will adapt in the subsequent calculation steps the turbulence parameters to the calculated average flow.

Table 1 : Intensities of turbulence for different cases of turbulence ( 81), ( 82)

Flows Descriptions Intensities of turbulence High-turbulence cases

High-speed flows inside complex geometries like heatexchangers and flows inside rotating machinery (turbines and compressors).

The turbulence intensity (u'/U) is between 5% and 20%.

Medium-turbulence cases

Flows in simpler devices like large pipes, ventilation flows etc. or low speed flows (low Reynolds number).

The turbulence intensity is between 1% and 5%.

Low-turbulence cases

Nearly still flows, like external flows across cars, submarines and aircrafts. Very high-quality wind-tunnels can also reach really low turbulence levels.

The turbulence intensity is very low, well below 1%.

THE COMBUSTION TERM ( 𝝎̇𝒌 )

Modeling the combustion rate in a turbulent environment is still an open question even out of the scope of RANS techniques. In the specific situation of industrial explosion CFD simulation, it would be sufficient to estimate the mean reaction rate ω̇ . But even this is not an easy task, partly because of the relatively limited understanding of the turbulent flame phenomenology and partly because of the limitations of the CFD technology. For instance, it is difficult today to implement a flame as an interface with its own dynamics in a standard CFD environment applicable to large scale explosions [START_REF] Daubech | Contribution à l'étude de l'effet de l'hétérogénéité d'un pré mélange gazeux sur la propagation d'une flamme dans un tube clos[END_REF]. Some examples of modeling approaches used are presented hereafter.

EDDY BREAK UP MODELS

The Eddy breakup modelling ( 68), ( 83) dates back from the early seventies and originates from the work of Imperial College with Spalding. Several versions exist.

Besides the final version proposed by Spalding, a second one "Eddy dissipation model" appeared shortly after followed by the "Eddy dissipation concept".

It is extremely important to recall that it is assumed that the turbulence Reynolds number (Re = • ) is large and that the Damköhler (

Da = = • •
) is also large. The first assumption suggests that turbulent molecular mixing is much more efficient in transporting the reactants than molecular transport and the second assumption says that the reaction time is much shorter than the "mechanical time" of the turbulence, indicating that the combustion process is not affected by the turbulence and continues to proceed at the molecular level, thus depending on local properties. In support of this, Spalding showed that the properties of a turbulent flame cannot be reproduced by applying the average field characteristics (reactants mass fractions, temperature,..) on the standard Arrhenius chemical kinetic law. Spalding further assumed, using a phenomenological model, that the turbulence completely disrupts the original laminar flame (at least under the assumption of large Reynolds numbers) and transforms it into a succession of unreacted and reacted pockets of gas (and not an "averaged" mixture of them), that combustion between them is possible only at the very end of the turbulence dissipation process when the molecular phenomena become dominant again and lastly that this combustion process is instantaneous at the scale of the flow. This vision is not disconnected from the present knowledge especially insofar as the combustion occurs at the molecular level but ignores the fact that the flame is not passive against the turbulence and that, in particular, it will "select" part of the spectrum of the turbulence, smoothing out the rest.

A mathematical formulation can be extracted from Spalding's analysis. Spalding considered the "reactedness", c, of the reaction rather than the mass fractions. This "progress variable" can be based on the mass fraction or temperature. Both definitions are identical when the (turbulent) Schmidt number is unity:

c = Y -Y Y -Y ≈ T -T T -T [99]
Where b stands for burnt gas and u stands for unburned gas. The Favre averaged transport equation of c is derived from that of the fuel mass fraction and reads:

∂ρc ∂t + ∂ρ c u ∂x + ∂(ρc u ) ∂x - ∂ ∂x ρD ∂c ∂x = ω̇ [100]
Reducing to:

∂ρc ∂t + ∂ρ c u ∂x = ∂ ∂x ρD ∂c ∂x + ω̇ [101]
As a consequence of the foregoing assumptions, the combustion rate is limited by the rate with which the turbulence is being dissipated and not by the chemistry (much faster). Applied to the equation above it means that ω̇ should be equal to the volumetric rate of dissipation of the "progress variable" fluctuations (𝑐" ). Since the combustion is assumed to start only after the complete dissipation of the turbulence, the equation for 𝑐" can be derived ignoring the chemical term. The relevant equation may be obtained on the same manner than for k (turbulent kinetic energy) in which the dissipation term appearing is ν . Transposed to the progress variable equation, the dissipation rate of the fluctuations in Favre average reads D • " " = χ. In a similar way than in the balance equation for the turbulent kinetic energy, the ratio between the variance of the fluctuations and the dissipation rate gives the characteristic dissipation time : 𝜏 = and 𝜏 = " . Since the flow is assumed to be dominated by the turbulence, the latter characteristic time should be proportional to the former via a constant of order 1 : 𝜏 = 𝐶 • 𝜏 .

Rearranging the last equations, it comes : 𝜒 = 𝐶 • " . Finally, going back to the original equation for the progress variable, ω̇ should be proportional to 𝜌̅ • 𝜒 according to the initial assumptions so that:

ω̇ = -C • ε k • ρ • 𝑐" [102]
Using the properties of the Favre averages, ρ c = ρ(c -c ) = ρ c -c . Because of the assumption of infinitely fast combustion, c is either 0 or 1 and because of this, it can be demonstrated that c = c so that finally:

ω̇ = -C ρ ε k c (1 -c ) [103]
Where C EBU is a constant of order one.

This formulation is particularly appealing since no chemistry is included and seems particularly suited for RANS approaches. In presently available explosion codes, a modified version of the original Eddy Break Up model is used such as the Eddy Dissipation Model proposed by Magnussen and Hjertager (84). These authors added two further considerations: first the "reactedness" is replaced by the estimator min Y , Y -Y and second, they estimated more accurately the "turbulence free zones" in which combustion occurs, stating that they are fed by the dissipation rate of the turbulence at the Kolmogorov scale (𝜏 = where η and u' η are respectively the Kolmogorov scale and related velocity fluctuations) rather than at the integral scale. Taking also into account the real volume occupied by the Kolmogorov eddies, using the isotropic and homogeneous description of the turbulence they achieved a refined description and using the relationship linking η and u' η to k and ε, the Eddy Dissipation model applied to the fuel mass fractions reads:

ω̇ = -A • ρ • ε k • min Y , Y -Y with A = - 23.6 C • √2 Re [104]
Where Re t is the turbulent Reynolds number (based on u' and l). The ratio A is assumed to be weakly variable and a constant value is prescribed, typically 20, in many CFD software (EXSIM, early versions of FLACS…).

Nevertheless, this model would be reasonable only if the average fuel mass fractions correspond to those in the vicinity of the combustion zone i.e at the molecular scale. Local discrepancies may have a strong influence especially if non-homogeneous situations are considered. The reaction zone occupies only a tiny volumetric fraction ξ of a computational cell. This parameter was extracted from foregoing analysis but according to Magnussen, the size of the zones where the reaction can take place is in fact smaller than the Kolmogorov scale 10 .The parameter ξ is used to make a link between the average value of Y F inside the cell and the local values of the fuel mass fraction in the immediate vicinity of the reaction zone inside and outside the reaction sheet. This refinement is known as the Eddy Dissipation Concept... The final expression is:

ω̇ = - 11.2 1 -9.77Re • ρ • ε k • Y -Y * [105]
Where Y * stands for the fuel mass fraction in the reaction zone considered as a perfectly stirred adiabatic reactor. The mean mass fraction Y * can be obtained from the linear combination of properties in the fine structures and the surrounding fluid.

This last model is used with some success in "industrial combustors" to predict the yield and pollutant formation. This is certainly the most refined version of the Eddy Breakup family models.

Note however that the turbulent combustion vision contained in EDM and EDC, assuming that combustion occurs locally in non-turbulent zones does not fit with the present understanding of turbulent flame propagation especially during explosions. It is thought that the flame interacts mostly, and strongly, with the largest eddies and dampens the incidence of the smallest. So the flame can be everywhere in the domain. The original version, EBU, would on that aspect fit better because the time scale is that of the largest structures.

CREBCOM ALGORITHM

CREBCOM algorithm ( 85) is an attempt to circumvent this difficulty. It was originally developed to model non-diffusive kind of combustion, purely convection driven like "fast" flames or forest fires. In other words, small scale phenomena like diffusion are theoretically neglected. In a frame moving with the flame front, equation [START_REF] Harten | Self adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF] may be written:

ρ • u • dc dx = ώ [106]
since, from the conservation of mass, ρ • u = ρ • S = cste. This equation is valid when c>c inf (c inf is the cold boundary limit). If Δ is the distance from the point in the flame where c=c inf and c=1 (e.g. the burning zone thickness), then:

ρ • S • 1 -c ∆ = ρ • S • ρ ρ • 1 -c ∆ ≈ ω̇ [107]
The density ratio can be expressed as function of c inf and of the expansion ratio of the gases across the flame front σ so that the density ratio in equation

[107] is c • σ + (1 -c ).
Since the chemical reaction term should drop to zero when c=>1, this expression can be generalized:

ω̇ ≈ ρ • S • [c • σ + (1 -c )] • 1 -c ∆ [108]
This is a way to introduce the philosophy of the CREBCOM model. It is used with an Euler type of equation, introducing directly the CREBCOM source modelling in equation [START_REF] Harten | Self adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF] under the following form:

ω̇ = ⎩ ⎪ ⎨ ⎪ ⎧ ρC Δ 1 -c , , F , , > 1 2 0 F , , < 1 2 , [109]
where F , , , is a control parameter (c ijk the progress variable in cell i,j,k), C g (also called K 0 ) is a combustion parameter (given below) and Δ the computation cell size assumed then to be equal to the burning zone thickness.

F , , for the cell (i,j,k), is calculated as

F , , = c , , + c , , + c , , + c , , + c , , + c , , -3c , , , [110]
If this control parameter exceeds the critical quantity of 0.25, the cell is considered to be burning, otherwise it remains un-ignited. As expected, the 'combustion model parameter' C g links the flame propagation model with the estimated flame burning speed S and expansion ratio σ via,

C = (A. σ + B)S , [111]
A and B are correlation constants obtained from numerical experiments and the preferred values are respectively 0.243 and 0.375. The CREBCOM algorithm is very simple to implement in a CFD code. It has been implemented in several industrial purposes CFD codes, such as TONUS, to investigate turbulent combustion flows.

The main numerical drawback of the CREBCOM algorithm is that it involves a binary criterion function that specifies whether the control volume is burnt or not. Numerical experiments have shown that this criterion function can create numerical oscillations in the pressure which strongly affect the flow when the flame speed is low with respect to the sound speed (low Mach number regime). Also, a bad definition of the flame velocity in entrance leads to a wrong prediction to the numerical approximation.

THE  FLAME MODEL

The β flame model was proposed by Arntzen [START_REF] Arntzen | Modelling of turbulence and combustion for simulation of gas explosions in complex geometries[END_REF]. In this model, the combustion modeling is divided into two parts, flame and burning velocity modeling. It resembles the CREBCOM model inasmuch a "flame velocity" S is preset and defines a combustion parameter C g looking very similar to the CREBCOM algorithm.

Originally the β flame model was proposed to enable to spread the flame thickness over a few cells (4 cells) while preserving the preset burning velocity. But this is done by significantly modifying the transport equation of the progress variable. D T is not extracted anymore from the k-ε model at least inside the flame front (in practice D T from k-ε is replaced by 4.Δ.S t only when 0.001<c<0.999). What does it change in the behavior of the transport equation of c? The steady state version of [START_REF] Harten | Self adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF] may be considered:

ρ • S • dc dx - d dx ρ • D • dc dx = d dx ρ • S • c -ρ • D • dc dx = ω̇ [112]
At the ignition point (c inf ) and across the burning zone Δ:

1 ∆ ρ • S • (1 -c ) -ρ • 4 • ∆ • S • 1 -c ∆ = ρ • S • 1 -c ∆ ρ ρ -4 = ω̇ [113]
Knowing that the density ratio across the burning zone might be between 5 and 7, this expression is indeed very close to the CREBCOM formulation. To some extent the transformation of D T in the original RANS equations for the progress variable transforms this equation in an Euler type of formulation.

THE POROSITY DISTRIBUTED RESISTANCE TO REPRESENT OBSTACLES IN URANS

The URANS technique with the standard k-Ɛ is not appropriated in boundary layers and in rotating flows. For instance, it cannot provide correct results in the wake and in the stagnation zone of obstacles. This is one reason why the porosity distributed resistance (PDR) method was introduced [START_REF] Lea | A Review of the State-of-the-Art in Gas Explosion Modelling[END_REF]. With this method, all the geometry is represented as a porous region. In regions with no obstacles, the porosity is maximum and in others the porosity represents the "blockage" offered by the obstacles and additional turbulence source terms are added to account for the wake effect. The model can be seen as a generalization of the Navier-Stokes equations for fluid flow and of Darcy's law commonly used for flows in porous regions [START_REF] Miguel | Wind-induced airflow through permeable materials, Part I: the motion equation[END_REF]. This method was first proposed by Patankar and Spalding [START_REF] Patankar | A calculation procedure for the transient and steady-state behaviour of shell-and -tube heat exchangers[END_REF].

The PDR formulation modifies the governing equations in two ways. Only nonblocked areas are available for fluid flow and obstacles occupying a control volume give an additional flow resistance and turbulence production [START_REF] Saeter | Modelling and simulation of gas explosions in complex, geometries[END_REF]. The volume porosity, β , is defined as:

β = 1 - V V + V [114]
where V is the volume available for the fluid and V is the volume occupied by the obstacles (Figure 123456789). Similarly, the area porosity β is illustrated in the Figure 1-10, where the surface in the x-direction, Ax is partially occupied by an obstacle, A .The area porosity is defined as:

β = 1 - A A + A [115]
The definition of area porosities in the y-and z-direction is similar. Obviously β v and β x vary between 0 and 1. The PDR formulation of the transport equation for the general variable ϕ(ϕ = u,v,w,Yk,k or ε is written as follow:

∂ β ρ ϕ ∂t + ∂ β ρ u ϕ ∂x - ∂ ∂x β μ σ ∂ ∂x ϕ = S + R [116]
Where β is the area blockage ratio vector, S is the non-obstructed component of the source term (due to the mean free flow), and R is the additional component of the source term caused by the obstructions.

The turbulent viscosity (μ ) is obtained using the two equations k-Ɛ model, which has been modified to include the additional turbulence generation from the subgrid scale objects. The production rate of turbulent kinetic energy is modelled as:

S = -β ρu u : ∂u ∂x and R = C μ u A + Σ C R u [117]
where C is a constant, A is the wetted area of the obstacles per unit of volume.

The first term in equation [117] represents the production of turbulence by friction forces along the surface of the obstacles. C is a constant vector that gives the fraction of the pressure drop, in each co-ordinate direction, that contributes to the generation of turbulence kinetic energy in the wake. R is the drag force vector, and is given by :

R , = -C 1 2 ρ u u [118]
where C is the drag coefficient which depends on the shape of obstacle. In regions containing sub-grid scale obstacles the turbulence kinetic energy dissipation rate is not obtained from the standard transport equation, but is calculated from the following expression

ε = C k l [119]
where l = C D , C is a constant and D is a typical obstacle dimension.

NUMERICAL ASPECTS

Globally the system of equations to be solved can be expressed as:

∂U ∂t + ∂F ∂x + ∂G ∂y + ∂F ∂x + ∂G ∂y -S = 0 [120]
where U are the unsteady terms, F and G the convective fluxes, F and G the viscous/diffusive fluxes and S the source terms. Explicitely:

U = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ ρ ρ u ρ v ρ Y ρ h + P ρ k ρ ε ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ , F = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ ρ u ρ u + P * ρ u v ρ u Y u ρ h + P * ρ ku ρ εu ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ , G = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ ρ v ρ u v ρ v + P * ρ v Y v ρ h + P * ρ kv ρ εv ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ [121] F = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 0 τ + ρ 2 3 k -u τ -ρ u v V , -ρ u Y λ ∂T ∂x -ρu h -ρ V , Y h , μ + μ σ ∂k ∂x μ + μ σ ∂ε ∂x ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ , G = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 0 τ -ρ u v τ + ρ 2 3 k -v V , -ρ v Y λ ∂T ∂y -ρv h -ρ V , Y h , μ + μ σ ∂k ∂y μ + μ σ ∂ε ∂y ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ [122] S = ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎡ 0 0 0 ω̇ ω̇ + τ ∂u ∂x + τ ∂u ∂y + τ ∂v ∂x + τ ∂v ∂y P -ρε -C ε k P -C ε k ρε ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎤ [123]
In CFD, these equations are "discretized" in time and space and solved step by step using various "numerical schemes" built in principle to cope with the "dynamics" of the phenomena to be studied.

∂ ∂t ρ ρU E ( , ) + ∂ ∂x ρU ρU + P (E + P)U ( ) = 0 or ∂w ∂t + ∂F(w) ∂x = 0 [124]
is a combination of acoustic travelling waves. Any perturbation would lead to the emission of pressures waves. Numerically speaking, truncation for instance, may generate oscillations, possibly amplifying. Artefacts may then be generated leading to the non-conservation of the mass for instance.

To overcome this difficulty, the convection terms are solved on "finite volumes" (Malalasekera,1995 (89)) by applying on them the standard conservation laws familiar to engineers (mass, energy, species, impulse). Doing this, the conservations laws are intrinsically satisfied but at the expenses of the accuracy since the local quantities (pressure, density,..) are homogeneous inside an elementary volume. There is some "smoothing effect" rendering more difficult the follow up of discontinuities like shocks or other interfaces (like a combustion front…).

However due to its robustness, the finite volume method is used in many numerical simulation codes: Fluent StarCCM+, CFX, FineTurbo...and especially those devoted to explosions.

Mathematically, the foundations of the "finite volume method" is obtained by integrating any conservation law of a physical parameter w over a volume Ω, involving a flux F(w) through its outer surface Σ (with normal outer n ⃗) and a source term S(w). In integral form, the conservation law becomes:

∂ ∂t wdΩ + divF(w)dΩ = 0 [125]
The Ostrogradski theorem leads to:

where ∮ F. ndΣ represents the sum of the flux through Σ. Considering now that Ω is a polyhedre (cube, tetrahedre,...) called "cell" and that the flux is supposed to be constant over each face, the integral is simplified to a discrete sum over each face of the cell.

F. ndΣ = F n Σ [127]
The quantity F = F(w ) is an approximation of the flux F over a face of the cell: it is the "numerical flux" through this face. It is further admitted that w is constant (homogenous) in each cell and it is equal to an approached value of its average over the cell.

Below, the integrated conservation law for w is expressed using an Euler explicit method where ∆𝑤 is the increment of the w between t and t+Δt:

∂ ∂t wdΩ = Ω dw dt = Ω ∆w ∆t [128]
This is a "first order" Taylor approximation of the time derivative meaning that the exact derivative is known within an accuracy on the order of Δt. Finally, the integrated and discretized conservation law within the frame of the finite volume method reads:

Ω ∆w ∆t + F n Σ = 0 [129]
Applied to a 1 dimensional situation, the "cell i", represented by the coordinated x i of its center is a line bounded by the left and right hand side borders at x i-1/2 and x i+1/2 (Figure 1 -11). The size of the cell is Δx = x ⁄ -x ⁄ . Finally, over the grid 𝑥 , 𝑥 , between the times n • ∆t and (n + 1) • ∆t , the "convective" balance of any quantity w of the system may be written as:

w -w ∆t + F -F ∆x = 0 [130]
F ⁄ is an approximation of the flux F(w) at the interface x ⁄ and at the time t + nΔt. It is the "numerical flux" at x ⁄ .

This numerical flux is to be evaluated as function of the averaged values of w in the neighboring cells. The way it is done is what defines the "numerical scheme". Note that the Euler equation is a "wave equation" any local change of w will result in a transportation from cell to cell via mechanical waves. Any local change of w will result in a transportation from cell to cell via mechanical waves. To follow correctly this transportation, the time step should be smaller than the time for the wave to travel between adjacent cells borders. This constrain is expressed as the Courant-Friedrichs-Lewy (CFL) condition:

(∆t) = CFL ∆x (𝑢 + 𝑐) [131]
Where ∆x max is the maximum value of the space step and |u+c| j is the characteristic velocity of the pressure waves. The CFL condition is a necessary condition for stability. The CFL number chosen by the user must be smaller than the maximum of CFL number insuring the stability. Its value depends on the numerical scheme and typically ranges from 0.25 to 1.

There are many types of numerical schemes [START_REF] Goncalvès | RESOLUTION NUMERIQUE DES EQUATIONS D'EULER MONO-DIMENSIONNELLES[END_REF] which are briefly presented below in groups or families.

TAYLOR APPROXIMATIONS OF THE FLUXES (FCT SCHEMES)

One route is to consider that a "natural" way is to assume that the numerical flux at face j±1/2 is the arithmetic average of F(w) calculated at the cell centers located on each side of the face j±1/2 (90). For instance, for face j+1/2:

F = F w + F(w ) 2 [132]
The discretized Euler equation then reads:

w -w ∆t + F w -F w 2∆x = 0 [133]
It can be demonstrated that the preceding equation is a bad estimation of the initial equation [124]. A Taylor expansion of the later reveals that the above discretized equation differs from [124] by:

+ ∆t 2 ⎩ ⎪ ⎨ ⎪ ⎧ ∂F ∂w F w -F w - ∂F ∂w (F w -F w ) ∆x ⎭ ⎪ ⎬ ⎪ ⎫ [134]
In this expression is the Jacobian matrix A.

Lax and Wendroff corrected this deficiency by choosing a formulation of the numerical flux F defined by:

F = F(w ) -F(w ) 2 - ∆t 2 A ⁄ F(w ) -F(w ) ∆x [135]
such that after complete development of the discretized equation a term identical to equation [127] appears giving:

w = w -∆t F w -F w 2∆x + ∆t 2 ⎩ ⎪ ⎨ ⎪ ⎧ ∂F ∂w F w -F w - ∂F ∂w (F w -F w ) ∆x ⎭ ⎪ ⎬ ⎪ ⎫
.

[136]

This elegant solution is balanced by the necessity of estimating . To avoid this, Ritchmeyer proposed in 1967 a variant of the Lax-Wendroff scheme. This calculation is based on two steps:

 The first step, called the predictor step, gives an estimation of the solution at the interface 𝑥 ⁄ at an intermediary time (𝑛 + 1 2)∆𝑡 ⁄ , based on the Lax-Frederichs scheme is expressed as follows:

w ⁄ ⁄ = w + w 2 - ∆t 2
F(w ) -F(w ) ∆x [137]
 The second step, called the corrector step, is the conservative discretization of the starting equation.

w = w -∆t F(w ⁄ ⁄ ) -F(w ⁄ ⁄ ) ∆x [138]
This procedure is approximative and residual oscillations may remain especially in the vicinity of discontinuities and an artificial viscous term (91), ( 90) is added to correct this. Such schemes are called "Flux Corrected Transport" or FCT.

It is clear that the artificial viscosity "patch" is only weakly linked to any physics and seems merely a numerical trick. Doing this, this numerical scheme may also destroy the real phenomena such as the transport of information via pressure waves.

"WAVE" APPROXIMATION OF THE FLUXES (FVS AND FDS SCHEMES)

The second route is to consider that the flux on each faces of a cell results from a wave transporting w. A first level of approximation is to introduce in the numerical fluxes, a splitting method (FVS for "Flux Vector Splitting"), via which the numerical flux is directly composed of two contributions (upstream and downstream) depending on the local Mach number. A second level, more accurate, is to estimate directly the flux on each face by adding up the contribution of all the waves originating for a difference of w on a face of a cell solving locally the "Riemann problem" (FDS or "Flux Difference Splitting").

Flux Vector Splitting (FVS) methods

The differential version of Equation [124] in one dimension can be written as:

∂w ∂t + ∂F(w) ∂x = ∂w ∂t + ∂F ∂w • ∂w ∂x = ∂w ∂t + a • ∂w ∂x = 0 [139]
which can then be integrated over a finite volume Ω as shown above. But here, "a" is some velocity with which w is transported in or from Ω. This is a "wave" equation. Considering this, Steger and Warming [START_REF] Steger | Flux vector splitting of the inviscid gasdynamic equations with applications to finite difference methods[END_REF] simply realized that the flux flowing through the cell j along one particular direction (x in the present case), comes from the "upwind" flow. If for instance a is positive, "upwind" means that the flux of w comes from the cell numbered j-1 at the face numbered j-1/2 and from the cell center itself for the face j+1/2 so that:

F = (aw)
and F = (aw)

[140]

A very practical way to distinguish between fluxes coming from the positive values of x or negative was then proposed

w -w Δt + a • w -w Δx + a • w -w Δx = 0 [141]
where a + =(a+|a|)/2 and a -=(a-|a|)/2. Related expressions for F j±1/2 are given in Table 2. One difficulty is to find a correct expression for the velocity a with which w is transported. As known form the reader, it might not be simply the material velocity but can be a sound velocity or a shock velocity depending on the initial conditions. Van Leer (95) proposed to adapt the expression of "a" to the velocity of the flow, according to the Mach number, M. For large mach numbers (M>1 or <-1), "a" is the material velocity of the flow implying that waves cannot propagate faster and transport additional w. For intermediate values, the transport of w by acoustic waves is introduced ( Warming The FVS methods can be considered as a rather crude approximations of the "wave" problem. FDS formulations are much closer to the physics. But much more complicated too and computer resources demanding. The principle is to compute exactly the characteristic wave velocities "a".

For ∶ M ≥ 1 F = F F = 0 For ∶ 0 ≤ M ≤ 1 F = F -F F = ρ(u -c) 2γ 1 u -c H -uc For ∶ -1 ≤ M ≤ 0 F = ρ(u + c) 2γ 1 u + c H + uc F = F -F For ∶ M ≤ -1 F = 0 F = F Van Leer For ∶ M ≥ 1 F = F and F = 0 For ∶ -1 ≤ M ≤ 1 F = ρc 4 (M + 1) ⎣ ⎢ ⎢ ⎢ ⎡ 1 2c γ (1 + γ -1 2 M) 2c γ -1 (1 + γ -1 2 M) ⎦ ⎥ ⎥ ⎥ ⎤ and F = - ρc 4 (M + 1) ⎣ ⎢ ⎢ ⎢ ⎡ 1 - 2c γ (1 - γ -1 2 M) 2c γ -1 (1 - γ -1 2 M) ⎦ ⎥ ⎥ ⎥ ⎤ For ∶ M ≤ -1 F = 0 and F = F
In a 1-dimension situation, when a difference in pressure (for instance) appears on each side of a face of a cell, three characteristic waves will appear and propagate (this is the "Riemann" problem):  a pressure wave propagating inside the low pressure region and separating the low pressure region from an intermediate pressure region;  a rarefaction wave propagating into the high pressure region and separating the high pressure region from an intermediate pressure region;

 a contact wave in the intermediate pressure region separating the two original fluids (coming from the low pressure region and compressed to the intermediate pressure and coming from the high pressure region and expanded to the intermediate pressure) and having potentially different densities. The contributions of these waves are added up to reconstruct the numerical flux at the interface integrating for each wave (∂F(w))/∂x=a•∂w/∂x. The characteristics of the wave and of the intermediate region can be calculated exactly using the Rankine-Hugoniot relationships. The complete resolution of the Riemann problem can be found in [START_REF] Hirsch | Fundamentals of Numerical Discretization[END_REF] for instance. The Godunov numerical scheme is issued from this resolution is accurate but very demanding in computer resources since the required parameters are implicit into the equations.

Roe [START_REF] Goncalvès | RESOLUTION NUMERIQUE DES EQUATIONS D'EULER MONO-DIMENSIONNELLES[END_REF] proposed an elegant approximation of the exact Godunov scheme. Roe linearized the Riemann problem to obtain explicit relationships. The approached Riemann problem is set as:

⎩ ⎪ ⎨ ⎪ ⎧ ∂W ∂t + A ∂W ∂x = 0 W = W = (ρ , ρ u , ρ v , ρ H ) if x ϵ C W = ρ , ρ u , ρ v , ρ H if x ϵ C [142]
where A is the matrix of Roe to be looked for. With Euler equations, the flux is not varying as W 2 , so that A cannot be simply evaluated as function of an arithmetic average between the variables W in two adjacent cells. Considering the nature of the Rankine Hugoniot relationships, Roe proposed a mixture of geometric and arithmetic averages:

⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎧ ρ = ρ ρ u = u ρ + u ρ ρ + ρ v = v ρ + v ρ ρ + ρ H = H ρ + H ρ ρ + ρ a = (γ -1) H - 1 2 (u + v ) [143]
From which the characteristic wave speeds can be directly extracted as:

λ = u -a λ = λ = u λ = u + a [144]
The coordinates of the principal directions of A can also be estimated (the first line of each vector corresponds to ρ, the second to ρu,…):

K = 1 u -a v H -u a , K = ⎝ ⎜ ⎛ 1 u v 1 2 (u + v ) ⎠ ⎟ ⎞ , K = 0 0 1 v , K = 1 u + a v H + u a [145]
The intensities 𝛼 of each elementary wave are estimated as:

⎩ ⎪ ⎪ ⎨ ⎪ ⎪ ⎧ α = ρ v -ρ v -v ρ -ρ α = γ -1 a ρ -ρ H -u + u ρ v -ρ v + ρ E -ρ E + v ρ v -ρ v -v ρ -ρ α = 1 2a ρ -ρ (u + a ) + ρ u -ρ u -a α α = ρ -ρ -(α -α ) [146]
Finally the numerical flux at the interface between cell i and cell j can be computed using the Godunov expression:

F W , W = 1 2 F(W ) + F W - 1 2 λ α K [147]
The Roe scheme is particularly accurate and efficient. Nevertheless, a drawback is that rarefaction waves propagate at the same speed than shock waves with the same "jump conditions". This due to the approximation. Typically, if two adjacent states are likely to generate a rarefaction wave, the Roe solver will connect directly the two states by a discontinuous jump, decreasing the entropy… which changes the nature of the flow and may even produce an unrealistic "expansion shock". Various solutions were proposed (93) such as an "entropy" correction consisting in modifying gradually the original characteristic wave speed. A simple an efficient method was proposed by Harten. The characteristic wave speeds need to be replaced by the following expression:

= • + 𝜀 only when 𝜆 𝑢 < 2 • 𝜀
The parameter ε s , is the "correction entropy parameter". Its value of is chosen, depending on the application, between 0 and 0.5.

ACCURACY: ORDER OF THE SCHEMES

Most of the preceding numerical schemes are a "first order approximation" of the differential operator inasmuch the numerical fluxes F are a function of the variable w to the first order only. Going back to equation [124] above this means that ∂F/∂w is a constant and the Euler equation reduces to a wave equation where "a" is the wave velocity:

∂w ∂t + a ∂w ∂x = 0 F W , W = 1 2 F(W ) + F W - 1 2 λ α K [148]
when an "explicit" formulation of the derivatives is chosen (see in the next section for more explanation), the discretized version of the above equation reads:

w -w ∆t + a w -w ∆x = Ο(∆x, ∆t)F W , W = 1 2 F(W ) + F W - 1 2 λ α K [149]
Solving for w i n gives:

𝑤 = 1 - 𝑎∆𝑡 ∆𝑥 𝑤 + 𝑎∆𝑡 ∆𝑥 𝑤 F W , W = 1 2 F(W ) + F W - 1 2 λ α K [150]
Clearly, the increment in w is proportional to Δt and Δx on the "first order", thus the numerical scheme is said "first order in time and in space". Very small space and time discretisation is required to improve the accuracy of the solution.

Apart from their simplicity, a further advantage of the first order numerical schemes is that they are "TVDs" (Total Variation Diminishing) meaning that they "naturally" dampen the oscillations that may appear (depending on the numerical solver chosen) for instance around strong gradients. By definition the Total Variation (TV) in a one dimensional configuration for the discrete solution w j n at the time t=n∆t is:

TV[w ] = |w -w | F W , W = 1 2 F(W ) + F W - 1 2 λ α K [151]
The discrete solution satisfies the TVD property if between (n+1)∆t and ∆t :

TV[w ] ≤ TV[w ]F W , W = 1 2 F(W ) + F W - 1 2 λ α K [152]
A numerical scheme is TVD if the numerical solution satisfies the TVD property. This property is such that the growing of a local gradient during the integration in time is going to be compensated by the (more important) diminution of a gradient in another place of the domain. Particularly, the appearance of oscillations at the vicinity of the discontinuities is not compatible with the TVD method.

For the first order numerical schemes, it comes out from the preceding equations and for rather small values of the CFL = ∆ ∆ ≪ 1:

TV[w ] ≈ 𝑎∆𝑡 ∆𝑥 • |w -w | = 𝑎∆𝑡 ∆𝑥 • TV[w ]F W , W = 1 2 F(W ) + F W - 1 2 λ α K [153]
So first order numerical schemes are TVDs. Note that this "good" properties is closely link to one drawback of such schemes: their diffusive nature. Because of this, they tend to smooth out any rapid variations even those which may really exist like shock waves or contact surfaces, not only numerical artifacts.To improve this situation, it is therefore necessary to introduce some schemes with upper order (at least second order). The extension to an upper order can be done using two routes: either a flux development to the second order as in the FCT schemes, or incorporating a space dependency in the conservative variables (w) used to calculate the fluxes ("MUSCL reconstruction (90)"). As explained above, the development of the fluxes to a higher order does not prevent the appearance of strong oscillations in the vicinity of shocks waves/contact surfaces so that …some artificial numerical viscosity is to be introduced, significantly altering the benefit of the second order development. So the second route, the MUSCL method, is frequently used at least in commercial CFD software.

To introduce the MUSCL method, the general formulation of a numerical scheme based on finite volumes discretization is recalled:

w -w ∆t + F -F ∆x = 0F W , W = 1 2 F(W ) + F W - 1 2 λ α K [154]
where the numerical flux 𝐹 ⁄ is a function of the neighbouring states on the left and the right of the interface x ⁄ . The variables 𝑤 ⁄ 𝑎𝑛𝑑 𝑤 ⁄ are respectively the values of conservative variables on the left and the right of the interface j+1 2. For a first order numerical scheme:

w ⁄ = w and w ⁄ = w F W , W = 1 2 F(W ) + F W - 1 2 λ α K [155]
In the MUSCL method (Monotonic Upstream-centered Scheme for Conservation Law) proposed by Van Leer in 1977, it is assumed that the variables w are not constant throughout each cell but vary linearly in space within each cell. Because the values of w at the interface 𝑥 ⁄ depend both from 𝑤 and from their slope a second order development in space of w is thus implicitly introduced. Two kinds of developments are possible Figure 1-12:

 Upwind approximation : 𝑤 ⁄ depends only on w j and w j-1 located on the left of 𝑥 ⁄ , as shown in the Figure 1-12, and similarly, 𝑤 ⁄ depends only on w j , w j+1 and w j+2 on the right of 𝑥 ⁄ .

w ⁄ = w + w -w 2 w ⁄ = w + w -w 2 F W , W = 1 2 F(W ) + F W - 1 2 λ α K [156] w ⁄ = w + w -w 2 F W , W = 1 2 F(W ) + F W - 1 2 λ α K [157]
 Centered approximation : 𝑤 ⁄ and 𝑤 ⁄ depend both only on w j and w j+1 , centered on 𝑥 ⁄ , A general formulation may be preferred:

w ⁄ = w + w -w 2 F W , W = 1 2 F(W ) + F W - 1 2 λ α K [158] w ⁄ = w + w -w 2 F W , W = 1 2 F(W ) + F W - 1 2 λ α K [159]
w ⁄ = w + 1 -ϕ 4 w -w + 1 + ϕ 4 w -w
[160]

w ⁄ = w - 1 -ϕ 4 w -w - 1 -ϕ 4 w -w [161]
Where the value of parameter Φ defines the type of MUSCL reconstruction:

 upwind if ϕ = -1.  centered if ϕ = 1.
Note that under the Flux approximation of equation [124], the centered scheme reduces to a first order approximation in space of the Euler equation. Upwind schemes are definitely of second order accuracy in space and are particularly suited to the resolution of Euler type of equation since they intrinsically have the capability of catching propagating information like pressure waves. Note further that a even higher order scheme can be obtained with the MUSCL reconstruction stating for instance that Φ=1/3.

Unfortunately, second order schemes are not automatically TVD. The MUSCL reconstruction can produce unphysical oscillations around the discontinuities and slip lines. Rather than introducing an artificial and arbitrary viscosity, the slopes of the variables in the cells are bounded during the MUSCL reconstruction so that the numerical scheme can remain TVD. They limit the slope of variation of the variables by comparing them to the slope in the neighbouring cells. The role of a limiter is to restrict locally the scheme to the first order, to avoid the appearance of non-physical oscillations. The limiter are introduced in the MUSCL algorithm as an additional function 𝛹(𝑟) where r is the slope of the variable:

w ⁄ = w + 1 -ϕ 4 Ψ(r ) w -w + 1 + ϕ 4 Ψ 1 r w -w [162] w ⁄ = w - 1 -ϕ 4 Ψ(r ) w -w - 1 -ϕ 4 Ψ 1 r w -w [163]
Quantities r and r respectively represent the slopes on the left and on the right, such as:

r = w -w w -w and r = w -w w -w [164]
Many slope limiters exist some of them being listed in the Table 3.

Limiters Expressions

Minmod Ψ(r) = max (0, min(1, r))

[165] No particular limiter has been found to work well for all problems, and a particular choice is usually made on a trial and error basis. In this thesis, the Minmod and the Van Albada limiters are considered.

DISCRETIZATION IN SPACE OF OTHER TERMS

The diffusive terms are usually discretized in space using a standard "centered scheme". The spatial discretisation of sources terms does not need a specific attention.

DISCRETIZATION OF TEMPORAL DERIVATIVES

"Explicit" or "implicit" resolution should be chosen. The wave equation can be used to illustrate the difference:

∂u ∂t + a ∂u ∂x = 0 [166]
where a is the wave speed.

One possible way to discretize this equation at point i and time n is:

u -u ∆t + a u -u ∆x = Ο(∆x, ∆t)
[167] This is an explicit expression i.e. the value of u at any grid point can be calculated directly from this expression without the need for any matrix inversion. Since u at each grid point can be updated independently, these schemes are easy to implement. On the downside, it turns out that this scheme is stable only under the CFL conditions which impose severe limitations over the time step leading to expensive calculations. But because of the ease of programming and easy modifications, this is the choice that was made in this work.

In an implicit scheme, the spatial derivative term is partly evaluated at the n timeindex:

u -u ∆t + a u -u ∆x = Ο(∆x, ∆t)
[168]

In this case, u at each grid cannot be updated from the other points independently. Instead a system of algebraic equations needs to be solved in order to calculate the values at all grid points simultaneously. It can be shown that this scheme is unconditionally stable so that the numerical errors will be damped out irrespective of how large the time-step is at least for the wave equation. Nevertheless, implicit schemes are not unconditionally stable for the full Euler or Navier-Stokes equations since the non-linearities in the governing equations often limit stability. However, they allow a much larger CFL limit than explicit schemes.

But the choice between both resolutions depends partly on the temporal behaviour of the rest of the equations to be solved. For instance, considering the time dependency of source terms, especially in combustion problems where the characteristics time scales for the combustion could be very different and sometimes independent from the scales of the flow (77) derived from the wave speed reasoning used to tailor the space and time steps (CFL).

Mathematical aspects can also be considered [START_REF] Cinnela | Modélisation et calcul d'écoulements compressibles turbulents[END_REF]. When the source term is linear and negative, explicit schemes may generate oscillations. Oppositedly, for linear and positive source term, implicit schemes may generate oscillations.

In practise, a first order in time explicit discretisation is often chosen (as in the preceeding equation). The time step can be chosen as the smallest between the convective time step given by the CFL conditions and a characteristic "diffusive" time scale given by:

∆t = (∆𝑥) 𝑎 CFL • ∆x + 2 • μ u i n -u i n-1 ∆t + a u i n -u i-1 n-1 ∆x = Ο(∆x, ∆t) [169] 1.2.3.10 MESH GENERATION
It is expected that the way the mesh is being implemented would influence the calculations. Three aspects are considered in this study.

The "mesh" refers to the arrangement of the cells inside a given geometry. Each cell is defined by the position of the "nodes", featuring the angles of the cells, and the shape/length of the curve linking the nodes of the cell and shaping the "faces" through which the fluxes will be calculated.

Basically two categories of meshes can be used. In the "structured" mesh, the calculation domain is subdivided in cells having the same geometry (parallelogram in 2D and parallelepiped in 3D). The physical aspect is that of a "grid". This technique is known to be cost effective because a limited number of parameters is needed to describe the mesh. The drawback may be the difficulty to describe precisely enough the irregular geometries.

A more flexible technique is offered with "unstructured" meshes. The nodes are distributed more or less arbitrarily inside the calculation domain and linked together to produce the cells, usually triangles (2D) or pyramids (3D). It is much more adapted to any kind of geometry and local refinement is relatively easy to implement. One of the drawbacks is certainly the generation and permanent use of a matrix of connectivity, describing the location of the nodes and geometry of each cell.

For complex geometry and multiscale physics, the mesh can rapidly contain a huge number of cells rendering the calculations very long. Since the scale of the phenomena usually vary depending on the location inside the calculation domain (a finer mesh might be required in a boundary layer or in the flame front), it could be very useful to adapt the density of cells inside the geometry. It can be done using adaptative meshing (the mesh is adapted during the calculation) but in the present study a non adaptative technique was investigated (AMA method).

Generation of a structured mesh

The classical way of implementing a structured mesh is to divide the initial computational domain in several blocks which are subsequently discretized in cells. The advantage of this technique is to choose, for each block, the shape and fineness of the cells the best suited to the description of the geometry and of the phenomena.

An easy way to discretise a block is to choose the same number/distribution of nodes (corners of the cells) the nodes with segments. Each cell (tetrahedre) can be cut again by joining the corners of the cell. A typical example is shown on clearly a preferential direction of the faces of the cells which may impact the calculations as illustrated in

Generation of an unstructured mesh

In the second method, a sort of systematic procedure might be applied starting from the knowledge of the boundary of the domain. First, the global geometry of the cell is defined. In a 2D domain, cells would be triangles. The smallest angle should be larger than 5° otherwise incorrect results would be obtained given the limited accuracy of the numerical methods. Then nodes are distributed along the external edges of the domain according to the choice of the modeler and a total number of nodes is prescribed. located on an external edge of the domain and progressing along the edges counterclockwise for instance. When this first loop is over, the procedure is reapplied on the internal nodes of the first layer of ce applied on each cell generation process are the  the barycenter of all the triangles is included in the computational domain,  the intersection of two adjacent triangles contains either a nod or a common edge,  the intersection of an edge of a triangle with the boundary of the domain is empty, or includes either two of these extremities or the entire edge,  the inner angles of the triangle must be larger than typically 5°.

An easy way to discretise a block is to choose the same number/distribution of (corners of the cells) for all parallel edges. The cells are created by joining h segments. Each cell (tetrahedre) can be cut again by joining the A typical example is shown on Figure 1-13 l direction of the faces of the cells which may impact the calculations as illustrated in the last chapter of this work. unstructured mesh

In the second method, a sort of systematic procedure might be applied starting from the knowledge of the boundary of the domain. First, the global geometry of the cell is defined. In a 2D domain, cells would be triangles. The smallest angle r than 5° otherwise incorrect results would be obtained given the limited accuracy of the numerical methods. Then nodes are distributed along the external edges of the domain according to the choice of the modeler and a total number of nodes is prescribed. Then the mesh is produced starting from a node located on an external edge of the domain and progressing along the edges counterclockwise for instance. When this first loop is over, the procedure is reapplied on the internal nodes of the first layer of cells and so on. The criteria applied on each cell generation process are the following: the barycenter of all the triangles is included in the computational domain, the intersection of two adjacent triangles contains either a nod or a common rsection of an edge of a triangle with the boundary of the domain is empty, or includes either two of these extremities or the entire edge, the inner angles of the triangle must be larger than typically 5°.

Page 81 sur 225

An easy way to discretise a block is to choose the same number/distribution of The cells are created by joining h segments. Each cell (tetrahedre) can be cut again by joining the 13). Note there is l direction of the faces of the cells which may impact the Generation of the structured mesh: around of the initial geometry and projection of nodes and edges of the block on the geometry, at the centre: structured mesh generated, at the right: unstructured mesh obtained the by re-cutting of the quadrangles of the structured mesh.

In the second method, a sort of systematic procedure might be applied starting from the knowledge of the boundary of the domain. First, the global geometry of the cell is defined. In a 2D domain, cells would be triangles. The smallest angle r than 5° otherwise incorrect results would be obtained given the limited accuracy of the numerical methods. Then nodes are distributed along the external edges of the domain according to the choice of the modeler and a total Then the mesh is produced starting from a node located on an external edge of the domain and progressing along the edges counterclockwise for instance. When this first loop is over, the procedure is lls and so on. The criteria the barycenter of all the triangles is included in the computational domain, the intersection of two adjacent triangles contains either a nod or a common rsection of an edge of a triangle with the boundary of the domain is empty, or includes either two of these extremities or the entire edge, the inner angles of the triangle must be larger than typically 5°.

The second procedure used is based on an original of loops in the domain, method that could be classified between traditional methods of spatial decomposition (construction of a hierarchical tree structure for the progressive location of points and the acquisition of metric optimization methods of Delaunay contraint The second method (Voronoï) is a mixture between the preceding one and an optimization method of Delaunay method [START_REF] Kolms | Computing unstructured finite element grids on free-form CAD surfaces[END_REF]. The "loop" is an assembly of oriented segment on itself. The initial loop is built by joining all the borders of the domain. A mesh size is required for each border providing an initial distribution of the nodes. This initial loop is cut by a line on which nodes are placed acco presented just below. The resulting loop continues until all the domain is covered by lines and nodes. The construction of a cutting line from a node is performed according to the following calculating an error criterion taking the following parameters into account (by order of priority): This method is particularly suitable for hybrid meshes (triangles and quadrilaterals in 2D), for complex geometries where the convex and concave areas taken into account. A smoother transition of mesh sizes is automatically ensured from areas with weak refinement to areas with strong


The second procedure used is based on an original method of iterative re of loops in the domain, method that could be classified between traditional methods of spatial decomposition (construction of a hierarchical tree structure for the progressive location of points and the acquisition of metric optimization methods of Delaunay contraint type. method (Voronoï) is a mixture between the preceding one and an optimization method of Delaunay (97) type and is based on a "loop recutting"

. The "loop" is an assembly of oriented segments forming a line closed on itself. The initial loop is built by joining all the borders of the domain. A mesh size is required for each border providing an initial distribution of the nodes. This initial loop is cut by a line on which nodes are placed according to the algorithm presented just below. The resulting loops are cut again with a line and the process continues until all the domain is covered by lines and nodes. The construction of a cutting line from a node is performed according to the following algorithm aiming at calculating an error criterion taking the following parameters into account (by order

The distortion expected of triangular elements by measuring the gap of angles generated compared to angles multiple of π 3 ⁄ . of the cutting lines introduced by favoring the length error corresponding to the gap between the number to prescribe on the cutting line to respect certain of nodes to distribute on the domain, the latter estimated with respect to the initial distribution of nodes on the loop.

Nodes are distributed on the cutting line as function of the length of the first and last segment on this line. Those latter lengths are the average of the segments connected to first and final nodes. This method is particularly suitable for hybrid meshes (triangles and quadrilaterals for complex geometries where the convex and concave areas taken into account. A smoother transition of mesh sizes is automatically ensured from areas with weak refinement to areas with strong refinement ( Page 82 sur 225 method of iterative re-cutting of loops in the domain, method that could be classified between traditional methods of spatial decomposition (construction of a hierarchical tree structure for the progressive location of points and the acquisition of metric information) and method (Voronoï) is a mixture between the preceding one and an type and is based on a "loop recutting" s forming a line closed on itself. The initial loop is built by joining all the borders of the domain. A mesh size is required for each border providing an initial distribution of the nodes. This rding to the algorithm are cut again with a line and the process continues until all the domain is covered by lines and nodes. The construction of a algorithm aiming at calculating an error criterion taking the following parameters into account (by order This method is particularly suitable for hybrid meshes (triangles and quadrilaterals for complex geometries where the convex and concave areas can easily be taken into account. A smoother transition of mesh sizes is automatically ensured refinement (Figure 1234567891011121314).

Finally, the quality of mesh obtained looks very close to that obtained by the use of Delaunay contraint type methods for which the limitations are mainly a lack of robustness of the algorithm in preserving the border integrity (especially in the three-dimensional case) [START_REF] Hirsch | Fundamentals of Numerical Discretization[END_REF].

To further improve the quality of the mesh, a smoothing algorithm is applied simply consisting in repositioning each node in the center of gravity of the neighboring nodes.

Readjustment (or re-adaptation) of the mesh Whatever the numerical scheme employed to improve the spatial accuracy, the numerical dissipation which would result from using a too coarse mesh in areas where steep variations appear (shocks waves for instance) can totally ruin the quality of the calculated solution. Further, the discretization error is in a proportion of the local cell size via the higher order terms neglected in the expression of the finite difference approximation of spatial derivatives. In areas where the parameters of the flow vary within a length scale comparable to that of the cells, this truncation error can be commensurate. To capture this, a refined mesh is certainly required at least locally. A further aspect need to be considered, particularly important for FVS algorithm which tracks the direction of the gradient to express the numerical fluxes. In practice, these fluxes are supposed to be perpendicular to each face of a cell. Non orthogonal components are ignored. An error is likely when the orientation of the mesh elements is not locally fitted with the direction of propagation of the physical information. It is possible to limit the influence of this by gradually decreasing the cell size.

However, with speed flows containing large areas of strong discontinuity, it is not possible to perform a mesh refinement over the entire domain without quickly facing serious computer limitations. Thus, only a local refinement is performed.

Most methods from the literature were initially based on the use of a "shock indicator" (a function estimating the gradient of the different parameters in a zone).

The areas where the indicator is too large are marked. In this zone the cells are recut. An example is given by [START_REF] Figueira Da Silva | Unstructured adaptive grid flow simulations of inert and reactive gas mixtures[END_REF]. A two-dimensional flow passing around a corner is calculated using triangular cells. In each cell and for all the important variables (pressure, longitudinal velocity, radial velocity, temperature and mass fraction of species), the gradient is calculated and normalized using the maximum amplitude of variation of each variable. Where the modulus of this gradient is greater than a threshold value (set by prior numerical experiments), the triangles are marked and divided into four sub-triangles (using the midpoints of the edges of the original triangle for example), while the neighboring triangles (which are not marked) are divided into two triangles. In order to not damage the quality of the mesh, it is however necessary to apply several times a smoothing algorithm which consists of moving successively each node of the mesh at the center of the polygon formed by all triangular elements having this node in common. The overall effect of refinement improves very significantly the quality of the solution. However, if a too low threshold is used, too many elements are added, which is incompatible with obtaining acceptable calculation time. If the threshold is too high, only very strong discontinuities can be properly captured. Moreover, it is difficult to know a priori the number of smoothing iteration required: a large number of them limits the added value of the refinement, a small number of smoothing results in a "bad" mesh quality, without possible control of the error induced. More effective methods are now available which consist not only in refining the mesh in areas of strong gradients, but to de-refine in other areas according to an estimator of the interpolation error.

The AMA method ("Anisotropic Mesh Adaptation") of ( 99) was devised for this purpose. The interpolation error of the solution is estimated by calculating the second derivative of the variables in each node. In the context of the present study the density of the flow is seen as a relevant variable (captures mixing zones, slip lines, shocks and other pressure waves). According to this metrics, nodes can be added in some zones and removed in others. The new optimized mesh has not only refined areas where strong gradients are met, but also no refined-areas related to zones where solution changes very little For a prescribed number of nodes (so for a given calculation cost), this procedure minimizes the interpolation error while distributing in the most relevant way the nodes in the computational domain. Anisotropy is also explicitly accounted for locally so that the cells are progressively oriented and deformed in a controlled manner parallel to the discontinuities which limits the error in the estimation of the numerical flux. To find the optimum mesh ("convergence of the mesh"), the solution calculated on the preceding mesh is interpolated on the new mesh to initiate a new calculation and so on until the topology of mesh does not change anymore. The fact the convergence in mesh is reached does not imply the solution is correct. It is only converged and optimized for a given number of cells. To estimate the discretisation error, the total number of cells needs to be changed.

BOUNDARY CONDITIONS

Prescribing the boundary conditions of the computational domain may be challenging. The resolution of the system of differential equations inside the computational domain requires the knowledge of the solution on the boundary, at least for a number of variables of the problem.

To do this highly simplifying assumptions are often made. Obviously, the choice of the boundary conditions needs to obey the physical laws such a speed vanishing near a wall. But the way they are implemented should also be compatible with the accuracy and stability of the numerical scheme. This could be a particularly acute problem in Euler type of flow (convective dominated) where the information given at the boundary may be propagated unduly in the calculation domain because of the "wave propagating nature" of the equations. To our knowledge, there is to date no comprehensive theory to determine a priori the number and choice of physical variables required to describe correctly a free boundary for the complete system of Navier-Stokes equations.

The common choice is based on the application of characteristic theory in one dimension of space. Basically, two methods are used. A "strong" coupling of the boundary with the calculation domain according to which the values of the variables are fixed on the nodes of the boundary (referred to as the Dirichlet conditions). In this situation the boundary will be reflective to the incoming pressure waves. A "weak" coupling is also possible according to which the values of the variables at the boundary are introduced in the mathematical expressions of the fluxes at the boundary (a variant of the Neuman conditions). Doing this, the expected values of the variables at the boundary can be found in the converged solution. Nevertheless, because of the flexibility introduced in the variables at the boundary, the convergence may be more problematic. The standard Neuman condition supposes the values of the fluxes at the boundary are fixed which is more in line with the finite volume approximation and the various numerical schemes. This condition is usually the first intention.

A set of rather standard boundary conditions compatible with the present CFD context are given in Table 4.

Boundary

Options offered inlet Dirichlet type. A velocity profile can be prescribed but also the profile. An alternative is to set a mass flowrate (velocity is then deduced). Temperature and pressure are set and the specific mass is deduced (equation of state). K and epsilon values are given as

k = 3 2 (U ref T i ) 2 ; ε = C μ 3 4 ⁄ k 3 2 ⁄ l ; l = 0.07L;
Where U ref is the flow velocity at inlet, T i the intensity of the turbulence suggested in table 1.

outlet "Compatibility" : Dirichlet condition for subsonic flows. The outlet pressure is estimated solving locally the Riemann problem connecting the outside static pressure to the conditions immediately upstream from the boundary. "Continuous" : Neumann conditions for supersonic flows. The local velocity is estimated on the basis of the flowfield immediately upstream.

solid walls "Law at the wall" : the tangential velocity is calculated analytically using the "loglaw" : "Slip wall" : The velocity coordinate normal to the wall is sel to zero (Dirichlet) but not the others. The convective fluxes at the wall are set to zero (except the pressure). "Non slip" : All the velocity coordinates are set to zero on the wall. "Coupled" : to adapat the wall temperature to heating by convection and cooling by internal conduction. non-physical borders "Symmetry" : all variables at equal disatnces from the border are the same (like a mirror) "Periodic" : Dirichlet conditions at the border for all parameters. Theoretically, it would have been possible to act as a simple user of the existing CFD codes targeted in this work. Alternatively or in addition to this, a CFD platform like OPENFOAM (free CFD software) could have been used to test specific models or numerical schemes. Such questions were debated at the early beginning of this work and the choice was made to start from scratch or nearly so for the reasons presented below.

Past existing benchmarking suggested that even well informed users were never fully aware of the entire details of the programming (constants, limiters,..) of the softwares they were using, because, among other reasons, the documentation could hardly be fully exhaustive. Perhaps, even the conceptors may not even be able to predict the exact behavior of their codes because of the number superimposed layers, added capabilities over the time. There is also a risk that the comparison of the numerical solutions may turn out into a fight between commercial interests as the present team witnessed during HySAFE research program.

OPENFOAM is an alternative. However, many hands co developed this open source code and, as a result, is very bushy and it proved very tricky to trace out the details of the programming. Because a high level of understanding and control upon the tool was required to meet the objectives of the present work, it was felt that operating with OPENFOAM was not the best option.

These last criteria, deep understanding request and good control of the programming, led us to choose a full programming of the models and numerical schemes that had to be challenged. Intuitively also we thought that by "doing it by ourselves" we would have better chances to find the traps into which a CFD code developer or even more a CFD code user may fall into.

MERLIN was developed on those grounds. It evolved progressively during this PhD work from a collection of isolated routines that we programmed entirely into a CFD toolbox in which various programmation choices can be chained.

It is built inside a MATLAB platform. Obviously, this platform provides a clean and simple but also effective coding, enabling for instance a "natural" handling of fundamental linear algebra objects like matrices, and providing an extensive amount of libraries and functions available for a lot of simple and complex operations like matrix manipulation or 3D graphics and visualization. MERLIN is not intended to be an additional commercial CFD software. It is only a research object, organized as a toolbox to test the behavior of the models and numerical schemes included in a range of commercial CFD software (dedicated to industrial explosions study).

Nonetheless, proceeding this way requires a carefully "verification" of MERLIN. In the first part of this chapter, a brief description of MERLIN and capabilities is given.

The second part is devoted to the verification exercise.

DESCRIPTION OF THE TOOLBOX MERLIN

MERLIN is a distributed code made of many routines which can be chained together depending on the problem to solve. Flexibility is the main advantage and the second one is the relative ease with which errors can be detected in each routine. The drawback is a relative lack of ergonomics. But again it is not an "optimized" CFD code. The "solvers"which are matlab files (e.g. "solving" Euler equations, RANS equations, processing the data,...) are located in folders, each of them containing dedicated MATLAB routines. The linkage between the "solvers" selected by the user and input data, respectively located in folder "Userinput" and in the routine "Statupdate", is done via the head programme "Main_program" which also pilots the calculations. Postprocessing (Folder "Postprocessor") is somewhat independent and is done manually by loading the data in a Matlab routine or in Tecplot for instance. The global organisation is shown on Figure 2-1 and a summarized description of the routines and contents are presented in the Table 6.

StatsUpdate (matlab file)

Manages the computational time which depends on the number of iteration necessary to reach the convergence of the computation. Userinput (folder)

In this folder, the user pastes the solvers(matlab files) copied in the other folders ("Euler equations folder" for example) that he wants to use to perform the simulations.In general, to perform a simulation, the userinput folder must contain all the matlab files (or solvers) related to: the fluid mechanic equations (Euler, URANS), the geometry, the type of mesh(structured or unstructured mesh) or mesh topology, the cell size, the physical parameters (density,velocity,pressure,mass fraction,enthalpy,eventually turbulence parameters...), the boundary and initial conditions and the numerical schemes. In the following, more details about the characteristics of the solvers directly implicated in the simulations are provided.

USER INPUT FOLDER

In this folder, the files (also called solvers or routines) corresponding to the case to be simulated are provided. These files (Figure 2-2) contain: fluid mechanic equations, geometric parameters, cell sizes and mesh topology, physical parameters, boundary conditions, initial conditions, numerical schemes and so on. In this figure, the file "Fully3D" contains the fluid mechanic equations and all the physical parameters required to represent the physical phenomenon.The file "BoundaryUpdate" was copied from the PreProcessor folder and pasted in the UserInput folder in order to set and update the initial and boundary conditions to run the "Fully3D" case. Likwise,the MeshOperGen3D file comes from the PreProcessor folder and was pastde in the UserInput folder to create the geometry, the mesh topology and the cell sizes to run the "Fully3D" case.The file "MUSCL+Roe"(routine dedicated to the second order Roe numerical scheme) comes from the Euler equations folder and was pasted in the UserInput folder to discretize the fluid mechanic equations chosen to represent the physical phenomenon programmed in "Fully3D" file.The connexion between these files is operated via the Main_program file.

The example of the simulation of a compressible fully three-dimensional nonreacting viscous flow, contained in the file "Fully3D" is given below. First, in the MATLAB editor, "Main_program.m" file is open and the line 'uifile='Fully3D.m';' below "% User Defined Script" section is to be uncommented which will load the 3-D user defined geometry (Figure 2-3) and boundary conditions data contained once "Main_program.m" would be launched. By running "Main_program.m", the program will start the calculation and at specific iteration intervals, a calculated data will be dumped on the hard disk in excel or in ascii format (Figure 234) and the simulation will stop once a user defined criterion is met. Beforehand, the user has to implement in the "Fully3D.m" file, the following information:  Geometry: an example of geometry buit using MERLIN is presented in the next section (dimensions in MERLIN are set in meter).  Mesh factor: the value of the Mesh Factor will be used by the file "MeshOperGen3D.m" (located in the PreProcessor folder). This value defines the mesh size in regular meshes.  Mesh topology: the coordinates of the mesh are defined. The "meshing" is operated by MeshOperGen3D.m (located in the PreProcessor folder).  Boundary conditions: the boundary conditions are defined and will be operated by the program BoundaryUpdate.m (located in the PreProcessor folder).  Some other options: the purpose is to define the software options, the model-related parameters and many other parameters useful to the numerical simulation. The most frequently used options in the present work are summarized in Table 7.

Options Descriptions

Possible choices typegeoin This instruction specifies the role of "MeshOperGen3D.m" routine which interprets the user-input data by "extruding" the initial 2D Y-Z polygon defined by the user into a 3-D solid.

Extrusion

TimeSteppi-ngMethod

Indicates the method used for time discretization.

'Implicit': Implicit discretization is used (This option is not retained in this thesis). 'Explicit': The explicit method is used.

Model

Choice of the type of fluid mechanic equations used.

Euler equations in one, two or three dimensions URANS equation without combustion model and URANS equations with combustion model igraf

The number of time step intervals with which the calculated data is to be displayed with Matlab's output graphics.

To be defined by the user isave

The number of time step intervals with which the calculated data is to be stored to the hard memory.

To be defined by the user toldiv "toldiv" defines the tolerance on the overall sum of the divergence equations. It is a threshold value used to define if a simulation is converged.

To be defined by the user typcflux Choice of numerical scheme -Roe scheme(one order) -Van-Leer scheme(one ordre) -MUSCL +Roe scheme with limiters of type minmod or Van This solver contains the routines (Figure 2-5 and Table 8) which "calculate" the geometry (in one, two and three dimensions), the mesh and ascribes the boundary conditions. In order to build the geometry defined in the Figure 23456, the user defines a basic 2D polygon by setting the coordinates of points 1 to 6 in the Y-Z plane. Then the user provides the thickness of the prism which extends along the X-dimension. The original polygon is thus "extruded" along axis X. Because of this convention, the faces orthogonal to X are rectangles ("type A" faces). The other faces, "type B", are ordinary polygons. In the case of velocity boundary conditions, the user specifies the mean velocity and the profile for each face.

POSTPROCESSOR

The prostprocessor file contains the routines to visualize and to analyze the fluid dynamic results. Again either MATLAB routines or Tecplot by (using ascii data files) can be used (Figure 2-10 and Table 9). approximations in space and one order explicit in time. For the second order in space schemes, "minmod" and "van Leer-van Albada" limiters are available to preserve the TVD (Total Variation Diminishing) property. To avoid the production of expansion shocks, the Harten entropy correction is added to the Roe scheme. Simulating real geometries in a realistic way often required a very large number of cells even if the geometry seems rather simple and even though AMA meshing readjustment is used. Millions of cells are very often required. The simulation cost can rapidly be prohibitive on a standard running PC. A solution is to implement "parallel computing" which is known to be efficient in CFD.

Today it is not necessary anymore to program in view of parallelizing and to allocate explicitly specific operations/tasks to specific processors. Dedicated routines can do this. In this work "Star-P software" was used to do this. Star-P is a client-server parallel-computing platform designed to work with high level languages (HLL) such as MagnusMATLAB or Python. It is associated to a parallel HPC (High Performance Computer). Built-in tools are included to expand HLL computing capabilities via added libraries and hardware-based accelerators.

Following, the programming effort in setting up such parallel computation is low. Star-P is a global array syntax language which can be operated in parallel computing using the "App" syntax as exemplified below. To create a random matrix and take its inverse with MATLAB, the following code lines can be used:

A = rand (100,100); B = inv (A);

The same operation can be performed in parallel computing (using Star-P) using the following instructions:

App= rand (100,100*p); Bpp = inv (App);

The "*p" syntax tells the data construction routines to build the matrix on the parallel HPC (High Performance Computer) associated to Star-P instead of doing this on the computer over which the MATLAB programme is run.

SOME EXAMPLES OF MERLIN TOOLBOX USE

Some situations illustrating the capabilities of MERLIN are presented below.

COMPRESSIBLE FLOW INSIDE A MANIFOLD

This simulation represents a continuous air flow in standard conditions through a manifold (Figure 2

SHOCK WAVES

An unconfined blast wave propagation case is presented below using first a 2-D Euler solver (Roe second order in space, one order in time and explicit, minmod limiter) and then a 3-D Euler solver (Roe second order in space, one order in time and explicit, minmod limiter).

For the 2D simulation, the domain is a square with solid boundaries, 10 m length, initially at standard conditions. The blast source is modeled as a localized region initially at 10-atmosphere over-pressure and 2000°C located in the center of the box extending over 5 cells. A 2500 cells structured and regular mesh is used. Some results are shown on Figure 2

COMBUSTION

The focus is now on a methane-air laminar diffusion flame. The setup consists in a co-annular burner made up of a 1.1-cm-diameter fuel tube located inside a concentric 10.2-cm-diameter air. The air and fuel flow rates are respectively set to 1300 cm 3 /s and 9.8 cm 3 /s so that the fuel (methane) and air mean velocities are 0.103 m/s and 0.161 m/s.

VERIFICATION OF MERLIN

Some capabilities of MERLIN were illustrated above. And in the selected examples, a good agreement either with the expected physics, or experiments could be shown. In the mind of many CFD users, this would be enough to claim that the code is to some extent "checked".

In fact, the code could work properly on selected examples and fail in other circumstances. The "checking" exercise needs to be stronger and in fact covers two separate exercises "verification" first and subsequently "validation".

To use simple words, a code is "verified" if all the solvers solve correctly (within the expected accuracy) the correct equations (no bugs). A code is "validated" if the calculated solution is coherent (correct trends, within the expected error band) with the expected solution for a given problem. Usually the validation is performed against experimental results and sometimes requires some 'tuning".

Obviously, verification precedes validation. Since MERLIN is a toolbox, the validation is to be performed once a particular "simulator" is arranged to as to solve a specific practical situation: for instance, a URANS non reactive solver with the low Reynolds k-epsilon model to simulate the formation of a flammable cloud in a box. This part of the work is postponed to the next chapter. In the following, the verification exercise is presented.

Formally, verification is defined by the American Institute of Aeronautics and Astronautics as "The process of determining that a model implementation accurately represents the developer's conceptual description of the model and the solution to the model" (AIAA G-077-1998).

The strategy is to identify the errors in the model implementation (bugs, bias,..) and in the calculated solution (residues). Then two aspects are considered: the verification of the code, which results in finding and removing the errors, and the verification of a calculation which consists in checking that the observed residual error is that expected from the numerical scheme (for instance second order accuracy).

In order to verify the code, which involves removing bugs, correcting any incorrect implementations of conceptual models, errors in inputs,…the developer should normally review the coding. But this is by no means a proof especially for complex codes like CFD. Basic consistency checks (ex : mass conservation) may help but again it is insufficient. The best way is to compare the computations with series "highly accurate" test cases. The latter could be analytic or well accepted numerical solutions. Experimental results should not be used since their physical reality may be somewhat different from the computed equations so that they could not be considered as a "highly accurate" test case. A grid sensitivity study should be performed to bring out potential errors especially arising from an incorrect coding of the numerical schemes. More generally all convergence aspects need to be checked (iterative, spatial and temporal). Note that all the options of the code should be examined.

In verifying the calculations, the accuracy of a given simulation is determined through a grid convergence analysis providing an "observed order of accuracy" and grid convergence indices (GCI). The latter are compared to the expected values for the chosen numerical scheme.

Four different methods can be used to verify a code: the Method of Exact Solutions, the Method of Manufactured Solutions, the comparison to benchmark numerical solutions, and the code to code comparisons. To the opinion of the present author, the two latter methods consists more in being confident about the code to be verified than to verify it really. The Method of Exact Solutions consists in comparing the numerical solution to an exact solution obtained analytically for a specific physical problem. The shock tube problem is a well known verification case… but it is nearly the only one. It is restricted to the Euler equation which is a restricted formulation of the general fluid mechanics equation (the diffusive part of the equation doesn't exist). It is clear then the method of exact solutions cannot help much to verify all the solvers of MERLIN.

For this reason, the Method of Manufactured solutions (MMS) which is a more general approach was used in the verification exercise.

PRESENTATION OF THE MMS METHODOLOGY

This method (100) was proposed thirty years ago but was applied to CFD codes in a more extensive manner only at the turn of the century. It is a sort of inverse method in which an especially devised analytical solution is forced into the differential equations (DE) to derive analytically resulting source terms. The latter are introduced into the numerical solver of the same DE to calculate numerically the solution which is then compared to the original analytical solution.

An example is given below. The heat diffusion equation is to be solved:

∂T ∂t -α ∂ T ∂x = g(x, t) [170]
where α is the thermal diffusivity, T the temperature and g(x,t) the source term (t is time and x the space coordinate).

One of the most general solution for this kind of equation reads:

T(x, t) = T exp(t t ⁄ )sin(πx L ⁄ ) [171]
where T 0 is the amplitude, t 0 a time scale and L a spatial scale (for instance the scale of the domain). When this solution is introduced in the differential equation the adequate source terms is obtained:

g(x, t) = 1 t + α π L T sin(πx L ⁄ ) 1 t exp(t t ⁄ )
[172]

Finally, the numerical solver to be verified is asked to calculate the solution, T, for the following governing equation:

∂T ∂t -α ∂ T ∂x = 1 t + α π L T sin(πx L ⁄ ) 1 t exp(t t ⁄ )
[173]

The last step consists in comparing the numerical solution to the exact one by varying the parameters of the computation that are known to have an influence on the error like the nature of the mesh (regular or not, mixed,… which impacts the calculated fluxes) and the mesh size which impacts the discretization error. In particular, the evolution of the discretization errors when refining the mesh should decrease at a rate giving the order of accuracy. The latter should be identical to the expected order of accuracy for the numerical scheme selected in the solver.

To do this comparison, the following indicators are used:

L norm = ∑ T , -T , N ⁄ [174] L ∞ norm = max T , -T , [ 175
]
where k is an integer (1, 2, 4, 8, 16…) defining the level of mesh refinement (when k is doubled, the mesh size is doubled by two so the larger k, the coarser the mesh), n is the node number (between 1 and N including both interior and boundary nodes with the exception of the Dirichlet boundary nodes for which the discretization error is always zero by construction). T , refers to the exact solution evaluated at node n and T , the numerical solution at the same node but calculated with mesh level k. The order of accuracy can be calculated taking into account two discrete mesh levels k and k + 1 as:

p = ln L L ln(r) [176]
where "L" refers to L 2 norm and L ∞ norm. The parameter r is the grid refinement factor between k and k+1 (ratio of the mesh size which is two in our example). The drop of the discretization error should scale as of 1/r . If the refinement factor is r = 2 and the expected order of accuracy is p = 2, the error should drop by a factor of four on each successively refined mesh level. If p = p then the solver is verified.

Nonetheless, the MMS should not be chosen erratically otherwise only a partial verification would be performed or unphysical solutions emulated. There are also mathematical requirements.

The MMS has to be a continuous function with continuous derivatives to allow the formal order of accuracy to be reached even on coarse meshes. Because of this, trigonometric and exponential functions are chosen because they are continuous, and infinitely differentiable. Note also that they are natural solutions of hyperbolic problems. Their derivatives never vanish including cross derivatives. In building the manufactured solution, it is important that all the terms of the governing differential equation have similar magnitudes to avoid the fact that one term totally dominates others. This prevents the larger magnitude errors from masking errors in other terms of smaller magnitude. For instance, the Manufactured Solution of the Navier-Stokes equations has to be chosen such that the convective and diffusive terms possess the same order of magnitude.

Note that the MMS methodology is used both to verify the code and the calculations.

APPLICATION TO MERLIN

A number of code options need to be included into the verification exercise such as the equation of state, the boundary conditions, the turbulence models, the combustion models,.. The expected order of accuracy is dependent on the numerical schemes. In most simulations performed with MERLIN, the Roe scheme with the minmod limiter (for the convective terms) and with the diffusive terms discretized using the basic centered scheme was used which is designed to be of the second order in space. About the discretization in time, the Euler explicit scheme is used being first order of accuracy.

The verification procedure should be applied both for the time and space discretization because errors in both scales are intricated. A procedure was once proposed by Kamm et al. (102) and modified by the present author. The underlying idea is that the global error can be considered as a linear combination of the time and space discretization errors. The procedure used to estimate the time discretization error is presented later.

First, the starting point is to evaluate accurately the errors due to the space discretization for which a multitude of factor may play a role.

It is strongly advised to verify the code with the most complex (hybrid) meshes which include several cell topologies with various skewness, aspect ratio, curvature, and stretching…. A solver is considered as fully verified if it succeeds on a hybrid mesh. But if the solver fails on a hybrid mesh, then it can be tested on simpler meshes to find any inconsistency on a particular mesh topology or a coding mistake.

Since 2D and 3D solvers are available, 2D and 3D meshes were produced with various topologies. Although the verification process was performed for all solvers (2D and 3D, Euler, URANS,…), only the 3-D URANS case is described below. But the verification exercise was also done for the other options.

MESHES PRODUCED FOR THE VERIFICATION EXERCISE

The most general mesh type includes meshes with mild skewness, aspect ratio, curvature, and stretching. The different mesh topologies considered in this verification exercise are classified as "structured", "unstructured", and "hybrid" meshes which are a combination of structured and unstructured meshes.The different mesh level of refinement (k) and mesh types used for the verification are given in Table 10.  The 3D structured meshes used contain hexahedral cells and are fully regular (the cells are cubes) or skewed curvilinear. The 3D structured meshes are "naturally" generated using MERLIN. The skewed curvilinear mesh is specially designed to test the effects of aspect ratio, skewness, stretching and the effect of curved boundaries on the code.  The 3D hybrid meshes contain hexahedral, tetrahedral, and prismatic cells.

3D
To isolate the cell quality effects, 3D hybrid meshes which have cells close to isotropic can be used to test the code. Again, all the 3D hybrid meshes are generated using MERLIN. Because only the errors due to the space discretization are looked for, a steady state manufactured solution was chosen on the following general form:

ϕ(x, y, z) = ϕ + ϕ (x, y, z)

[177] where :

ϕ (x, y, z) = ϕ f a πx L + ϕ f a πy L + ϕ f a πz L + ϕ f a πxy L + ϕ f a πyz L + ϕ f a πzx L
where ϕ = [ρ, u, v, w, Y , p, k, ε] represent the parameters of the flow and the f (x) functions are sine or cosine functions. Again, the different constants in the above equation (Φ 0 , Φ x ,..) are chosen so that all terms in the URANS equations are on the same order of magnitude. For the same reason, the value for the molecular viscosity is set to 10 N/m so that the contributions of the inviscid and viscous terms are similar. After some trials and errors the following functions were retained: ρ = 1 + 0.15Cos(0.75πx) -0.1Sin(0.45πy) + 0.1Sin(0.8πz) + 0.08Cos(0.65πxy) + 0.05Sin(0.75πyz) + 0.12Cos(0.5πzx) u = 70 + 5Sin(0.5πx) -15Cos(0.85πy) -10Cos(0.4πz) + 7Cos(0.6πxy) + 4Sin(0.8πyz) -4Cos(0.9πzx) v = 90 -5Sin(0.8πx) + 10Cos(0.8πy) + 5Cos(0.5πz) -11Cos(0.9πxy) -5Sin(0.4πyz) + 5Cos(0.6πzx) w = 80 -10Cos(0.85πx) + 10Sin(0.9πy) + 12Cos(0.5πz) -12Sin(0.4πxy) + 11Sin(0.8πyz) + 5Cos(0.75πzx) Y = 0.75 + 0.11Cos(0.6πx) -0.075Sin(0.35πy) + 0.075Sin(0.65πz) + 0.06Cos(0.5πxy) + 0.035Sin(0.6πyz) + 0.09Cos(0.4πzx)

ϕ (x, y, z) = ϕ + ϕ (x, y, z)[F(x, y, z) -C] [178]
where m is an integer. When m=1, the Manufactured Solution is equal to a constant ϕ at the boundary condition satisfying a Dirichlet condition. For m = 2, the Manufactured Solution will satisfy both Dirichlet and Neumann (zero normal gradient) boundary conditions along the specified boundary. Note that the value of m can be different from one variable to another.

VERIFICATION PROCEDURE

The example of the 3-D URANS solver (e.g. continuity, momentum and total energy equations) is given below. The equations are tested on the structured, unstructured and hybrid meshes. As explained before, the "source terms" resulting from the introduction of the preceding functions in the URANS equations are imposed into the 3-D URANS solver to re-estimate numerically the functions. Then the later are compared to the exact analytic function and the L and L ∞ norms of the discretization error are calculated for all flow parameters. This is performed for a range of mesh refinement levels in order to extract an order of accuracy. For the numerical scheme chosen, in MERLIN a second order accuracy is expected. An observed order of 2 is shown in Figure 2-24 on the 3D skewed hybrid mesh. The solver is then verified. All the variables discretization norms approached 2 as expected, except for ρƐ the turbulent dissipation rate which drops to zero. To explore the reason for this failure of the verification on the 3D skewed hybrid unstructured mesh, tests were done on simpler meshes, as advised by the conceptors of the MMS methodology. On a highly skewed 3D structures curvilinear mesh with hexahedral cells (Figure 2-26) and on an unstructured mesh (Figure 2-27) with tetrahedral cells, the verification is successful.

of the discretization error

There is certainly an issue in the discrete formulation of some of the terms in the turbulent dissipation rate equation which does not stand into the coding, since it works as expected on certain mesh topologies, but probably on the expression of the cross-diffusion terms in the turbulent dissipation rate equation which are strongly dependent on the topology of the mesh.

EXTRACTION OF THE INFLUENCE OF THE TIME DISCRETISATION

It is more difficult to apply the verification procedure using the order of accuracy test to problems that involve both spatial and temporal discretization, especially when the spatial order is different from the temporal order. A combined spatial and temporal order verification method was developed by Kamm et al. ( 106). These authors use the Newton-type iterative procedure to solve a coupled, non-linear set of algebraic equations to calculate the coefficients and observed order of accuracies for the spatial and temporal terms in the discretization error expansion.

In this present work, a simplified approach is proposed suing only the higher order terms of the Taylor expansion, the norms (like L 2 and L ∞ ) of the discretization error can be written as:

∈ = g h + g h [179]
where p and q are the orders of accuracy in space and time that we are looking for. In this equation, g and g are coefficients, to be estimated, linking the global error to h and h , the levels of discretization in spatial and temporal scale.

The first step is to perform a spatial mesh refinement study, like above, with a fixed time step to calculate p and g using three mesh levels so that:

∈ = g h + ϕ [180]
where ϕ = g h is the fixed temporal error term. Using three mesh solutions, refined by the factor r , coarse ( r h ), medium ( r h ), and fine (h ), the observed order of accuracy p can be calculated as :

g ( r h ) -g ( r h ) g ( r h ) -g (h ) = ∈ -ϕ -∈ -ϕ ∈ -ϕ -∈ -ϕ g ( r h ) r -1 g (h ) r -1 = ∈ -∈ ∈ -∈ p = ln ∈ -∈ ∈ -∈ ln(r ) [181]
Then, the coefficient g can be calculated:

g = ∈ -∈ h [182]
As a second step, a temporal refinement study is performed on a fixed mesh to calculate q and g using three temporal discretizations, coarse ( r h ), medium ( r h ), and fine (h ) just as above.

In the third step, the spatial step size and the temporal step size have to be chosen so that the spatial discretization error term (g x .h x p

) is on the same order of magnitude than (g t .h t q ). With this precaution, the relative importance of time and space discretization are the same avoiding small errors on a scale to be masked by large errors on the other one, rendering difficult to assess the order of accuracy on the time scale 11 . It can be shown that this condition is reached when 𝑟 = 𝑟 ⁄ where 𝑟 is the temporal refinement factor and 𝑟 is the spatial refinement factor. Remember that in our case, the formal order is one in time and two in space.

Using this procedure, the unsteady time term is verified on the 3D hybrid mesh. As expected, the observed order of accuracy approached one with mesh refinement (Figure 2 Apart from the difficulty presented above about the turbulent dissipation rate, the verification was achieved on the most complex meshes. MERLIN can be now used to check (but again it is more a "validation" exercise) some models used in commercial CFD softwares used to simulated large scale explosions.

ANALYSIS OF SOME MODELLING ASPECTS IN LARGE SCALE (INDUSTRIAL) EXPLOSION SIMULATION

As shown in the introduction, the intercomparison of the results of the codes on practical cases is not satisfactory. And this is what motivated this work. The purpose is not to consider the codes in their entirety but some key modeling aspects: numerical schemes, turbulence model, obstacles, combustion modeling,…

When attempting to simulate an explosion within an industrial application context the first difficulty is to be able to simulate correctly a relatively large spectrum of interlinked physical processes. It is not enough to be capable of propagating a flame in a flammable cloud. It is necessary first to simulate the formation of the cloud since this "initial" step of the scenario of the accident has been for long time be recognized as a leading aspect. In some cases, it might also be necessary to calculate the propagation of the pressure wave away from the exploding cloud.

But certainly, one of the most challenging issue is to be able to account for very large, typically tens of meters dimensions and sometimes very complex geometries (Figure 3-1) and in particular with regards to the propagation of the flame. To overcome this difficulty, in many software devoted to industrial explosions (FLACS, EXSIM, COBRA, AUTOREAGAS…), the choice is made to model (rather than calculate) the smallest scales of the flow (typically those of the turbulence) and the smallest scales of the geometry as well so that large enough computational cells could be employed. URANS techniques appear as a good mean to do this. The various models and numerical schemes employed were presented in chapter 1 and implemented in MERLIN. Some significant phenomenological shortcomings were already identified that could easily be used to rule out one modeling strategy or another but whatever the softwares, even those still to be developed, assumptions and models would always be required to try and tackle the complexity of the real world.

The objective of this chapter is to estimate the incidence of various modeling strategies on the predictive capability within the specific context of large scale industrial explosions. This predictive capability can only be assessed within a specific context, here, the large scale industrial explosions.

Then what are the relevant "typical" situations?

The severity and extent of the explosion would unavoidably depend strongly on the size, shape and composition and turbulence level of the explosive cloud. As a first step, the latter parameters need to be estimated. As far as large scale explosions are concerned, large clouds are expected either resulting from a massive leakage from a compressed storage in the open atmosphere or from an accumulation in some semi-confined area (in a large building for instance). The fate of the subsequent explosion would also depend on the geometry of the surrounding which might incorporate some confinement and obstacles.

These various situations were studied using the CFD techniques commonly implemented in CFD softwares (Table 5), via the toolbox MERLIN, in the light of available experimental data and present phenomenological knowledge.

But before investigating the incidence of the physical modeling (flow, turbulence, combustion,…), a focus is put on the incidence of the mathematical/numerical aspects. The performances of numerical schemes are linked to the discretization of the geometry. Both aspects are discussed below.

INCIDENCE OF THE NUMERICAL METHOD

It expected that the following parameters would impact the quality of the resolution:

 The mesh topology (structured, unstructured mesh, geometry,…);  The nature of the numerical scheme. In particular, a "flux" splitting method (Van Leer scheme) and a "Riemann" solver method (Roe scheme) are considered which represent two different but relevant strategies in the present context 12 ;

The incidence of the numerical method is better investigated using purely eulerian problems since numerical schemes were developed to solve the "convective/advective" parts of the Navier-Stokes equations under the finite volume description. 12 Other solvers and were also challenged but the Roe and the van Leer schemes can trace out pressure waves, in addition of being capable of calculating steady flows, since their underlying physics is that of wave propagation. This is not the case for many other numerical schemes which were tested during the course of this thesis but then dropped.

As shown in the preceding chapter, the Euler equations the second order (in space) versions of the numerical schemes can be derived from the first order version by applying the MUSCL method (upwind). This strategy was used hereafter. As numerical dispersion may result from this transformation, front "limiters" are applied. Two of them are challenged: "minmod" and "Van Leer Van Albada" (VLVA). The time discretisation is an explicit first order scheme 13 .

For all the numerical simulations performed in this chapter, the time step is such that the CFL number is 0.8.

To assess accurately the incidence of the numerical schemes, the 1-D shock tube problem is used because steep gradients are produced which can be calculated on a theoretical basis providing a strong comparison/evaluation basis.

But the shock tube configuration is not sufficient for investigating the influence of the meshing and of the geometrical configuration which both can produce artefacts (which also depends on the numerical scheme). Following, well known 2-D configurations, with obstacles, are also considered.

INFLUENCE OF THE NUMERICAL SCHEMES: SHOCK TUBE CONFIGURATION

The shock tube (Figure 3-2) is a straight pipe equipped with a diaphragm inside (covering the full cross section). The "driving" section (on the left of the diaphragm) is filled with a high pressure ideal gas (the variables are labelled "L": P for the pressure, ρ for the density, T for the temperature and U = 0 for the initial velocityU = 0). The "working" section is filled with a low pressure ideal gas (the variables are labelled "R": pressure P , density ρ , temperature T and initial velocity U = 0). At time t = 0, the diaphragm vanishes, a process of pressure averaging appears. The high pressure gas is expanded by an expansion wave moving to the left. The expanding gas flows in the working section and produced a shock wave moving to the right and compressing the low pressure gas. A shock wave is created in the low pressure gas which is compressed. The gas coming from the driving section is separated from the working section gas by a contact surface often called "slip line" in a 1 D and 2-D situation. The shock wave and the slip line are discontinuities: across the shock wave all the variables "jump" (P(x), or ρ(x), or T(x), or U(x)). Across the slip line velocity and pressure are constant but most aften not the density. Note that across an expansion wave, the variables vary progressively so that the entropy is maintained. This is the Sod (103) problem for which exact solutions can be obtained solving the one dimensional "Riemann" problem. 13 One order in space versions of these numerical schemes were also tested but they show bad performances produced significantly smeared aout profiles. The results may even be difficult to interpret because of this.

In the present simulations, the geometry of the tube is given in Figure 3-3: 1 m long, 0.1 m wide with the diaphragm in the middle. In the driving section the specific mass of the (ideal) gas is 1 kg/m 3 , its pressure 100 kPa. In the working section, the specific mass of the gas is 0.125 kg/m 3 , its pressure 10 kPa. The initial velocity is 0 m/s.

The horizontal boundaries are solid wall with "slip wall" conditions (the normal component of the velocity is zero but the tangential velocity is maintained). Dirichlet conditions are imposed on the left to ensure the isentropic evolution is respected (this avoids a potentiel drift of the flow parameters). A compatibility condition is used on the right to provide a gradual adjustment of flow. The tube is discretized was discretised with 1305 nodes: using an unstructured Voronoï mesh and in some cases a structured/divided mesh (a structured/divided mesh is obtained by re-cutting the rectangles of the structured mesh along the diagonal).

Figure 3-4 shows the solution obtained at time 0.6 ms using an unstructured mesh and a second order numerical scheme. In this case, the accuracies of both numerical schemes are comparable and results are within a few % from the theoretical values.

Slight differences appear in the front of the waves (Figure 345) suggesting the limiter Van Leer Van Albada would be more accurate. Nevertheless, the thickness of the shock wave (Figure 3-5 right) is not zero, as it should be, but extends over 2 to 3 cells. With further refinement of the mesh size (Figure 3456), the shock wave seems better reproduced but not the expansion wave suggesting the influence of the numerical diffusion plays differently on different waves. This illustrates the difficulty to accurately simulate both the expansion and compression waves using "upwind" schemes. Although, in the present case, refinement of the mesh does not significantly improve the results, the nature of the mesh might be of some importance even in this simple configuration. The simulation presented above were performed using an "unstructured mesh" i.e. not Cartesian.

When a structured/divided mesh is used instead, the profiles downstream of the shock are dissimilar (Figure 34567) without any physical reasons behind (for instance, the velocity is larger near the top plate…).

Nothing like this appears with the unstructured mesh. This may be explained by the approximation of the fluxes, being approximated along the alignment of the faces of the nearby volumes i.e. along preferred directions in a structured mesh. A few percent difference is possible due to such details of the meshing. This point is may have much larger consequences as shown in the forthcoming.

Figure 34567: Illustration of the dissymmetry of the profiles appearing downstream of the shock with a structured mesh: Roe 2 nd order scheme with a "Minmod" limiter from the bottom (Y = 0, 005) to the top (Y = 0, 095) of the tube.

GEOMETRY INDUCED EFFECTS

REFLECTION OF AN UNSTEADY SHOCK ON A COMPRESSION RAMP: CASE OF SMR (SINGLE MACH REFLECTION) REFLECTION

The case of the propagation of a shock on a ramp is known to produce complex shock structures. On Figure 3-8, the incident shock wave propagates from left to right and impinges on a ramp. Two different reflection patterns are observed depending on the angle of the ramp () with the horizontal plane. When the latter is large enough, a "Regular" Reflection (RR) occurs according to which the incident shock (I) continues its propagation, unaffected, and a reflected bow shock (R) propagates backwards. The angle of the reflected shock depends on the intensity of the incident shock and the streamlines of the flow are bent towards the ramp. With a small enough angle, the angle of deflection of the streamlines becomes larger than that of the ramp producing a sort of stagnation zone just behind the contact point of the incident shock with the ramp. A third shock wave results, called the "Mach stem" M. This last shock propagates perpendicular to the ramp and creates downstream a layer into which the streamlines deflected by the reflected shock (R) can develop. As a result, a contact surface (G) or "slip line" appears separating two fluids having the same pressure but not the same velocities. The reflected shock, the "slip line" and the Mach stem represent the Single Mach Reflection (SMR). This problem is now well documented and is a reference case for evaluating the performance of numerical algorithms dedicated to supersonic flows. First, a shock tube problem is created in a straight section without a ramp by setting the pressures, temperatures and densities in regions 1 and 2 (Figure 3456789) to obtain the desired incident shock wave Mach number. The mesh is adjusted so that it corresponds to what is to be used as a second step on the field containing the ramp. The numerical profile is allowed to settle around of the shock discontinuity, until any oscillations related to the initialization is dampened or located sufficiently until the numerical diffusion thickness of the shock (related to the numerical scheme used) is steady. Second, this numerical solution is moved to the desired initial position in the field containing the ramp. The ramp is located sufficiently far from the inlet section so that the reflected shock cannot interact with the initial position of the shock (which could produce artefacts). angle between the direction of the ramp and the longitudinal direction, 𝑋 : relative position of the corner of the ramp relative to the inlet (distance/ length of the computational domain),𝑋 : relative position of the incident shock at the end of the simulation. Because of the excellent agreement between numerics and experiments the Takayama and Jiang (104) case (Figure 3-12 and Figure 3-13) can be considered as a reference case. Apart from the vertical incident shock, the expected "Mach stem" is visible, at right angles with the ramp on the leading edge of the shock complex, so as the bow reflected shock and the "slip line" trailing downstream and attached to the junction of the shocks ("triple point").

Figure 3-14 shows the isocontours of density fields using a van Leer numerical scheme with the minmod and VLVA limiters on a structured mesh. The reference mesh proposed by Takayama and Jiang (mesh with 255 × 256 cells) was used. But, it was verified (see later) that this level of refinement is satisfactory. In the present situation, it can be seen this anomaly is observed for both numerical schemes while the van Leer scheme is not known to be subject to this anomaly (107). It was further verified that this anomaly persists regardless of the level of refinement. Furthermore, whatever the level of entropy correction used (with the Roe solver), this anomaly remains while the entropy correction is said to make it disappear.

This was an unexpected result. To identify the reasons, the "blunt" body configuration studied by Gressier known to favour the appearance of the "carbuncle" phenomenon was numerically investigated. "Carbuncle" is an unexpected solution (numerical artefact) that appears when calculating supersonic flows over a blunt body. The strong shock solution is destabilized by the body producing a complex flow/shock structure (Figure 3-15). Using a mesh containing the same number of cells than used by Gressier (80 x 160 nodes), we simulated a particularly severe flow configuration, with an upstream flow Mach number of 10 (3480 m/s, 101325 Pa, 298 K), using the 2 nd order Roe scheme with a zero entropy correction. In Figure 3-16, the isocontour of the longitudinal velocity field is presented. It can be observed that, although the isocontours are somewhat irregular (presumably due to the mesh topology), the shock remains upstream of the blunt body without any appearance of the "carbuncle" phenomenon (Figure 3-16). Thus, the appearance of the kind of anomaly of Figure 3-14 should depend mainly on the mesh topology. Switching from a structured mesh to an unstructured mesh containing the same number of elements induces a complete disappearance of the problem (Figure 3-17).

Furthermore, the geometric characteristics of the shock and flow structure obtained are very similar to those obtained by the other numericians and by the experimentalists. Note that a relatively high level of refinement has to be used in order that the structured obtained downstream of the Mach stem is sufficiently representative. The larger the height of the Mach stem, the greater the amount of numerical instabilities which may be propagated. Reducing the angle of inclination of the ramp keeping the same incident Mach number enlarges the Mach stem. Such a typical configuration is the experimental case of Deschambault and Glass (105) are similar to those of . Note that in this experiment, the incident shock wave comes from the right and is inclined. This situation was simulated in this work (Figure 3 Some discrepancy is observed downstream from the Mach stem, in region 2. This is consistent with the acknowledged difficulty of accurately predicting the structure of both shocks and slip lines with the kind of numerical schemes implemented in the range of CFD sofwares targeted in this PhD work. Because of this, questions are raised about the ability of such numerical CFD methods to deal reliably with shocks and obstacles in which a number of shock and flow configurations can appear. This point is addressed below.

MULTIPLE REFLECTIONS

The following (standard) case was originally proposed by Woodward and Colella (108) and was used by many authors (see Dolejsi, (109) for example). In a channel (Figure 3 An exciting region appears in the corner of the step where an intense expansion occurs. The flow which is compressed upstream of the step experiences a strong expansion from a subsonic state ahead of the step to a supersonic state on the top of the step over a very small distance centered on the corner of the step. A special care is needed to mesh the domain avoiding zones with too coarse mesh, destabilizing the calculations. It appeared that the van Leer scheme was not robust enough whatever the effort to refine the mesh. Results were obtained with the Roe solver coupled with a "Minmod" limiter. A rather small entropy correction (0. As compared with the expected structure, three major differences appear. First the shock structure is shifted further downstream, second, three reflections instead of two are observed and third a Mach reflection appears on the top of step while a regular one is expected. This remains true whatever the level of the mesh refinement. Again, the corner is the center of an expansion fan through which the fluid particles are accelerated from a subsonic condition upstream from this point to a supersonic one downstream. It is known (preceding chapter) that approximation Riemann numerical schemes like Roe's version do not conserve entropy as it should in the present test case. This feature appears clearly in Figure 3-25 where the entropy drops strongly past the corner. Because of this the density of the gas drops artificially so that the flow is insufficiently accelerated (Figure 3-26). Because the intensity of the expansion fan is so strongly underestimated due to the approximation made in the Riemann type of solver implemented presently, the curvature of the first reflected shock (on the top of the channel) is too large so that the angle of reflection is affected. This affects the angles of all the subsequent reflected shocks, their number and even their nature (mach or regular reflection). This situation is likely to occur in cases where a singularity (like an edge) is located on the path of the pressure wave.

Apart for playing locally either on the entropy correction or on a more or less "manual" modification of the flow, it would be necessary to use a more accurate numerical scheme ("true" Riemann solver, higher order,…) at the expenses surely of the cost and may be also of the robustness.

IMPLICATIONS

In the preceding section, the potential incidence of some numerical aspects pertaining to large scale 'industrial" explosion CFD codes were addressed.

In the specific context of this study, a typical user will be a safety engineer, trained into the main physical aspects of explosion scenarios, but having only a general knowledge of the use of a CFD software. It means that even though he will be able to build a scenario of accident incorporating details of the conditions and to elaborate on the consequences based on engineered safety tools (using for example, the TNO guides), he will not claim to be a specialist in fluids mechanics, nor in combustion. In CFD, his culture may even be more superficial. He will know enough to understand that a mesh is required and that testing the "convergence in mesh" is necessary.

Consequently, many of the numerical and mathematical aspects of CFD will stand out of his scope to a significant extent. Unfortunately these aspects may have an impact on the simulation:

 within the finite volume approximation, the conservation laws of the fluid mechanics are applied in small volumes so that the convective fluxes are calculated along the normal direction of each face. This induces a sensitivity of the solution to the topology of the mesh used. Firstly, a preferred orientation of the faces of the cells may induce to a non-physical dissymmetrization of the velocity profiles. For sonic flows, the error may amount a few percents. The reason for this behaviour is the propagation of numerical instabilities produced in the vicinity of shock waves or discontinuities. In complex geometries, these instabilities may generate totally wrong shock structures. To avoid this, the use of a unstructured mesh is required ;  the solvers used in the considered CFD tools are chosen for their rapidity (cost) and robustness and are approximate estimators of "wave"-like behaviour of the convective part of the fluid mechanics equations. The van Leer and Roe solvers, as used in this work, were chosen as being representative of the two main families of Riemann solvers (FDS and FVS).

In the presence of intense expansion waves (expansion over hundreds of m/s), the Van Leer scheme becomes unstable. The Roe scheme combined with a "Minmod" limiter is much more robust but nevertheless ill estimates strongly the intensity of expansion waves. And there is no real simple way to solve this issue ;  because of the preceding shortcoming, they fail to properly consider the presence of singular points, like sharp angles, as soon as a fast flow is considered producing large velocity variation (typically for a mach number above 2). Many practical situations are concerned: computing a high velocity flow emerging from a pipe in the open atmosphere, computing the interaction with obstacles… in the latter situation non physical waves may even be produced so that the results may be…unreliable.

FORMATION OF A FLAMMABLE CLOUD

The focus is now on the physical aspects of large scale 'industrial" explosion CFD codes.

In this section, we challenge the ability of the URANS systems of equations associated to the k-epsilon model to predict the formation of a flammable cloud following a leakage of pressurized gas into an obstructed/confined environment. Several areas in the flow need be distinguished raising specific difficulties as explained below:  In the "transitional" zone (near-field and intermediate-field : Figure 3-27) extending up to 20 times the orifice diameter, the flow is first rapidly depressurized to the outside pressure (expansion zone near the orifice also called near field ranging from the orifice to a few orifice diameter downstream), and the core flow slows down and entrains a sufficient amount of the outside atmosphere to become further downstream a "selfsimilar" flow structure ("free jet") which is located in the far-field zone. This is an important region since all the momentum of the core flow is transferred to the atmosphere.  In the "self-similar" region (downstream from more than 20 times the orifice diameter), often referred as the "free jet". As long as the flow does not impact obstacles or walls, the momentum is conserved. Often, in commercial codes, the turbulence, which produces the plume, is modelled using the standard k-epsilon approach.  The interaction of the flow with obstructions is very challenging. First, it could be very demanding in computer resources to take into account a large number of obstacles in terms of boundary conditions, mesh size, preparation of the calculation and computing time,… Second, the k-epsilon model may not perform well around the obstacles. The PDR technique once appeared as an appealing solution to such difficulties but there is certainly a need to check the validity of this approach. Unfortunately, although the PDR technique can be programmed in MERLIN, there was not enough time and data to investigate this to a sufficient extent.

Finally, the accumulation of the flammable atmosphere in a confinment (a building in the industry) may also be a challenge because the convective part of the flow may become on the same order than the diffusive part and body forces may come into play affecting the turbulence and then the distribution of the flammable mixture. This aspect is addressed in the following. The instrumentation is displayed on a sort of vertical cross installed perpendicular to the axis of the flow. It can be placed at different distance from the orifice. The sensors (Figure 3-30) consist in 14 speed measurement probes able to capture the instantaneous axial and transverse velocities, and 9 concentration probes (oxygen depletion meters). The accuracy is about 5% of the measured value. The simulations were performed using the non-reacting flow solver of MERLIN. A cylindrical coordinates system was chosen. A real gas equation of state (Abel-Noble equation) was used because it is known that high pressure hydrogen releases cannot be accurately simulated using the ideal gas equation ( 115).

Because of its robustness (as shown in the preceding section), the 2 nd order in space, first order (explicit) in time Roe numerical scheme was used with a "minmod" limiter (entropy correction factor=0.05). The diffusive terms were discretized using the basic central differencing scheme. An upwind differentiation is used for the convective terms.

The simulation was performed in two steps. First the transitional zone is computed and second the computed solution is introduced as a starter in a larger domain to compute the plume.

THE TRANSIONAL ZONE

The physics of the expansion and intermediate zone was described earlier (111,112,113,114,115,116).

To our knowledge, computing the expansion zone remains a challenge especially when the jet is strongly under-expanded (110). This is why, with the CFd codes considered in this work, a virtual source term is calculated and replaces this zone when performing CFD simulations. This is the "notional nozzle" approach and it is based on the theory of Birch (106). With this method, instead of modelling the actual jet source, the source is modelled as a sonic flow at atmospheric pressure with the same mass-flowrate as the original high-pressure jet. Because it seems that some outside air can be entrained into the expansion zone, of the entrainment of air into the hydrogen jet inside the expansion zone, this model may lead to significant discrepancies.

The ratio between the static pressure at the nozzle (P ) and the ambient pressure (P ), called the "Nozzle Pressure Ratio" (NPR). It is the key parameter that determines the overall structure of the jet. A jet is considered highly underexpanded when NPR>10.The structure of this kind of jet is shown on Figure 3-36. o where: o M is the Mach number at the nozzle; o and P is the pressure in the tank.  D : diameter of the Mach disk, which is twice the distance between the position of the triple point and the axis. According to Antsupov (118):

D D = log P P - 3 4 .
 L : length of the subsonic zone defined as the distance between the position on the axis of the Mach disk and the point of the axis reached by the inner sonic line of the shear layer. According to Glotov(119):

L D = 1.96 × P P .
 a : initial angle between the direction of the free boundary of the jet and the axis.  R : radius of initial curvature of the free boundary of the jet,  a : angle between the direction of the reflected shock and the axis of the jet.  D : diameter of the main shock structure. According to Avduevskii(120) :

D X = 0.73 - 2.21 √Re .
Here, Re is the Reynolds number at the nozzle.

The main difficulty is the NPR ratio effects on the expansion, on mixing with the outside air and consequently on the density differences. A specific numerical investigation of this was performed hearafter with MERLIN.

The simulation is performed in a cylinder centered around the axis of the jet. The orifice is located on a face of this cylinder (Figure 3-37). Dirichlet conditions are applied at the inlet. The pressure at the orifice is 37.7bar, the temperature 298K, the hydrogen mass flow rate is 253g/s, the discharge coefficient Cd=0.94 and the volumetric concentration of hydrogen is 100%. The turbulence parameters at the orifice are initialized using the estimators of Table 4.

Around the orifice, outside air is flowing at 2.5m/s (corresponding to the average velocity of the wind during the experiments) at atmospheric pressure, 298K and with a turbulence intensity set at 1%.

A "continuous outflow" boundary condition is applied to the outlet plane and the other surfaces are specified as "non slip solid" boundaries to favour the stability of the numerical simulations.

Because shock waves are produced, an unstructured mesh is implemented. Knowing that the accuracy of the solutions depends on the mesh and on its readjustments when the structure of the flow changes, the sensitivity of the solutions to the mesh refinement was first investigated.

The evolution of the axial velocity profile is shown for example in Figure 3-41 for the different meshing strategies (unstructured) presented on Figure 3-39 and Figure 3-40.

Precautions were taken to ensure that the solution is sufficiently converged at the end of each simulation so that the readjustment of the mesh could be applied. The initial solution obtained from the initial mesh (1205741 nodes) is relatively spread (Figure 3-39 and Figure 3-41). The first application of the adaptation algorithm (setp 2 of the readjustment) ends with a slightly larger number of cells (1358974 nodes) distributed around the mixing layer, the zones of oblique shock and of the Mach disk so the calculation is more accurate and the calculated shock structure shifts upstream. At the beginning of the third step, the cells are moved and concentrated around the new geometry of the shock structure whereas a large number of elements are removed at the periphery of the mixing layer. From the end of this third step and up, the geometry of the shock structure doesn't evolve anymore and the subsequent steps mostly refine the alignment of the cells along the shock discontinuities, the sliding line and the boundary of the jet. Finally, three steps and 1000000 elements are enough. Air immediately behind the leading shock is heated up by compression. As the shock wave propagates away, it dissipates due to the flow divergence. Meanwhile, another semi-spherical shock wave is developing inside the expanding hydrogen jet. As it establishes, it intensifies due to the strong under expansion of hydrogen and its front becomes nearly flat. Finally, the well-known shock structure, Mach stem, is formed. Upstream the Mach disk, the temperature is extremely low, so as the pressure, and downstream, it is higher than ambient and the pressure too. The temperature might seem somewhat unrealistic but this may come from the limitations of the equation of state. As the jet emerges from the orifice, radial expansion waves originate around the circular edge. The expansion waves propagate downstream and reflect as compression waves at the outer boundary. The coalescence of these compression waves results in the barrel shock structure surrounding the highly supersonic region upstream the Mach disk. The flow inside the region between the flow boundary and the barrel shock is still supersonic, but its Mach number is much lower than that inside the barrel shock. Downstream the Mach disk, the pressure restores to a value slightly higher than atmosphere and then drops back to ambient at X=10D. In line with the changes in pressure, the axial velocity increases to a maximum value upstream the Mach disk but immediately downstream the flow decelerates to a rather low velocity and then accelerates again to a second peak value at X=10D. After that, it decelerates slowly due to the momentum transfer in the plume. Because of the high underexpansion, hydrogen is cooled down to a very low temperature (about 3K), then elevated to a value slightly above the exit temperature immediately behind the Mach disk. Accordingly, the temperature drops gently to a second minimum value at X=10D again due to the further mild expansion of the jet. After that point, it increases slowly back to the environmental temperature. At distances smaller than X=4D (roughly the position of the Mach disk), the hydrogen mass fraction on the centerline remains unity, then decreases gradually immediately downstream the Mach disk. The position X=10D is a critical distance corresponding by definition to (X + L ). Note that 10 is the initially (using the physical correlations) estimated value of (X + L ). It features the end of the near-source expansion of the jet. Downstream, the centerline velocity starts to decrease due to the momentum exchange with the outer flow. The flow conditions at this position are therefore considered as the most suitable for the subsequent plume simulation.

The mean axial velocity and hydrogen mass fraction at this location are plotted on Figure 3-44 also shows that the mass fraction of hydrogen is less than unity even on the centerline implying air entrainment starts well before the end of the expansion zone. This suggests that the notional nozzle approach, which neglects air entrainment, may induce significant discrepancies. In the present situation, a 50% difference in the mean axial velocity and in the H 2 mass fraction would result from the used of the notional nozzle approach.

In the present case the average velocity on the axis at X=10D is about 600 m/s (as compared to about 1000 m/s in the notional approach) and the hydrogen mass fraction is about 90% (as compared to 100 % in the notional approach).

THE PLUME

The conditions of Figure 3-44 are now used as input condition for the simulation of the plume development. The computational domain is a quarter of cylinder 120D wide (120D = 1.44m which is approximately the distance between the orifice and ground) and 1000D long (1000D = 12m which is approximately the farthest distance where experimental data are available). The mesh contains 1620000 grid cells clustered near the opening. This number of grid cells was selected as a result of a mesh sensitivity analysis comparable to that performed for the expansion zone. The boundary conditions are similarly defined also.

The centerline values of the mean axial velocity, turbulent intensity and hydrogen mass/volumetric fractions along the centerline are plotted on Figure 3-45. There is a sort of plateau up to X=25D for the mean axial velocity and hydrogen mass/volumetric fractions. Downstream, the velocity decreases quickly to 100 m/s at X=75D due to the intensified air entrainment between X=25D and X=75D induced by the growth of turbulence produced by some shear layer instability. From this point on, the velocity decreases according to the well known hyperbolic law of free jets. The evolutions of mass/volumetric fractions of hydrogen are quite similar indicating that the self similar structure of the jet is reached. When X increases, both the axial velocity and the hydrogen mass fraction decay in the axial direction and spread in the radial direction normal to the centerline due to the exchange of momentum and hydrogen with air. Note that the expected Gaussian like profiles are observed when X is larger than about 100D.

COMPARISON WITH THE EXPERIMENTAL DATA

The comparison with experiments is proposed from Figure 3-47 to Figure 3-49 concerning the axial velocity profiles, the volumetric concentration of hydrogen for positions, the reduced profiles of radial velocity and volumetric concentration respectively at x=4.5m and at x=3 m from the orifice, and the reduced profile of the intensity of turbulence at x=4.5m from the orifice. Note that data are available only after about X=100D. The numerical results are in fair agreement, within a few percents, with the experimental data indicating that the physics is well captured not only in the plume region but presumably also in the expansion zone since significant primary air induction is produced in this zone which would have a significant impact on the results if not accounted for correctly.

From the numerical results, it is observed that the expected Gaussian like profiles are observed but only when X is larger than about 100D which is commonly admitted by the experimentalist. Because of this and because also the experimental and numerical results are in agreement, the standard k-eps model seems reasonably valid even if the Taylor assumption may not seem totally valid (i.e. u' very small as compared to U).

COMPUTING THE INFLUENCE OF THE CONFINEMENT

In many practical situations, the leakage occurs in a confined environment. The formation of a potential stratified mixture is to be considered. Along the past decades, experiments were performed to investigate this aspect and provide data for modelling.

EXPERIMENTAL SETUP AND RESULTS

Among the experimental works, that proposed in GAMELAN (121), ( 122) and (123) test series are particularly well documented (Figure 3-50). The enclosure is a parallelepiped, 1.26 m high and with a 0.93m × 0.93m square section. A 90 cm wide and 18 cm high vent is located on a side wall at the top of the box. Helium is injected an injection tube pointing upwards at 210 mm from the bottom of the box. The axis of the injection is the same than that of the box.

Two experiments are considered here corresponding to two different flow-rates: 10 and 180 Nl/min respectively through injection tubes of 20 and 5 mm of diameter respectively.

Helium concentration sensors (catharometers) were placed on three vertical masts (M1, M2 and M4). All are located off the axis of the injection. The details of the arrangement are given on Figure 3-51. In this figure, the labelling of the different sensors of the mast 4 are given (the general labelling is MiNj where "j" is the number of the sensor -the sensors are numbered in ascending order going from bottom to top -and "i" indicates the number of the mast). Three dimensional simulations were performed using the "Non-reacting flows" solver of MERLIN with a buoyancy source term. The standard k-epsilon model incorporating the buoyancy driven turbulence production term was used as a first intention.

As above, the 2 nd order in space, first order (explicit) in time Roe numerical scheme was used with a "minmod" limiter. The diffusive terms were discretized using the basic central differencing scheme. An upwind differencing is chosen for the convective part.

Dirichlet conditions are given at the inlet of helium: ambient pressure and temperature, turbulence intensity = 1%, mass flow is given, volume fraction of helium = 1.

At the outlet (the vent), Dirichlet conditions are also used: ambient pressure and temperature, turbulence inflow intensity = 1%, inflow mass fraction of helium=0, flow lines normal to the vent.

On the walls, zero pressure gradient condition was applied and the "law of the wall" is used for the turbulence parameters and the tangential component of the velocity.

At the beginning of the injection, the chamber contains air only at Standard Temperature and Pressure (STP) conditions. At this moment, turbulence has not been developed in the chamber, but for calculation purposes, some turbulence parameters have to be defined to avoid divisions by zero. Later when the flow field would have generated its own (much stronger) turbulence, the influence of the initial conditions would have disappeared. When free-stream turbulence effects are negligible, the turbulent intensity is on the order of 1% of the free stream velocity U ∞ .The turbulent kinetic energy would be k ≈ 10 U ∞ and Ɛ should be set, so that the effective viscosity is of the same order than the laminar viscosity.

υ ~0(υ)~10

[183]

Thus Ɛ = C k υ ⁄ ≈ 10 U [184]
The turbulent length scale should then be set to

l = C k ⁄ υ ⁄ ≈ 10 U ⁄ [185]
When the initial velocity is zero, a satisfactorily effective viscosity can be obtained with

k = Ɛ = 10 [186]
Geometrical considerations show that a plane of symmetry is available so that only half of the domain could be considered. Furthermore, in the present configurations (Figure 3-54 and Figure 3-55), a more or less conical plume emerging from the orifice is expected impacting the upper wall to which a significant part of the momentum is transferred. Due to the buoyancy forces, the mixture impacting the wall accumulates progressively as a layer which thickness increases with time until the moment the volumetric flow through the upper vent equilibrates that of the impacting plume. If the flow-rate and initial momentum at the inlet is large enough the layer would completely fill the chamber producing an homogeneous mixture (Figure 3-54). On the other hand, for low flow-rate, a stratified regime is observed (Figure 3-55) (124).

Accordingly, the zones to be meshed more finely would be the plume and the upper part of the chamber. The opening angle of the plume was supposed to be 12° as for free jest. The other parts of the chamber can be meshed less finely. A typical mesh built on this basis is presented on Figure 3 A mesh sensitivity study was performed according to the conditions shown in Table 13. Table 13 The numerical result obtained for the lowest sensor of the mast 4 for the different meshes of Table 13 are with mesh3 which was later retained for the subsequent simulations. The numerical result obtained for the lowest sensor of the mast 4 for the different . Clearly, convergence is obtained with mesh3 which was later retained for the subsequent simulations.

Meshes

umerical results of the lowest sensor of the mast 4 obtained with flowrate

A typical sequence representing the evolution of the mixture volumetric Nl/min test case) is At t = 0.01s, the jet is formed. At t = 1 s, the gas mixture impacts the top of the box and at t = 10s, a gravity current is established and along the ceiling a layer begins to form. At t = 50s, the reversal phenomenon is achieved and produces a still front at the height of the discharge point but the mixture remains unevenly distributed in this layer. At t = 300s and whatever the next moment, the layer is now homogeneous and is steady. This ascertainment is also confirmed by the profiles of the concentration taken at locations 1, 2, 3, 4 and 5 of mast 4 (Figure 3 The simulated volumetric concentration on the upper "sensors" is close to that measured. But some stratification is calculated (-10%) whereas in the experiment, the mixture is perfectly homogeneous throughout the chamber.

The results for the second experiment (10 Nl/mn) are shown on Figure 3-60 (simulations and measurements superposed). In this case, the discrepancies are very large both in trends and in absolute values suggesting that the physics is not well accounted for. Nonetheless, the uRANS formulation in general and the k-epsilon model in particular is expected to perform correctly in this situation as suggested by the free jet situation studied before. A deeper analysis is required. Non dimensional numbers can be used and in particular the Reynolds numbers comparing the inertial and viscosity forces and the Froude numbers comparing the inertial to the buoyancy forces. Two flows need to be considered: that issued from the injection and that occurring along the boundaries of the enclosure. To calculate the Reynolds and Froude numbers of the enclosure, a convective flow velocity should be estimated. It is basically that of the plume impinging on the top of the enclosure and could be estimated assuming the momentum is preserved from the injection point to the top (to estimate the convective velocity, the buoyancy forces are ignored). The final expressions are provided in Table 14.

IMPLICATIONS

Somewhat surprisingly, even in these rather simple looking situations, obtaining a satisfactory estimation may be a challenging task. The good news is that the numerical scheme is of secondary importance contrary to what was observed as for the shock waves. Nevertheless, a refined knowledge of fluids mechanics seems required to avoid/solve the traps which may not be available to the "standard user" as defined before:  it was first shown that the "notional nozzle" approach often used to simulate underexpanded gaseous releases in CFD can be wrong to a significant extent (by 50% in the present case for significantly but realistically underexpanded gas releases). It is believed that this discrepancy is larger, the greater the underexpansion level. Otherwise, somewhat surprisingly, the assumption that u' needs to be very small as compared to U under the Taylor assumption in the context of the k-epsilon model seems to remain valid in free jets ;  when a confinement is to be accounted for, using a standard well established k-ε model, even provided with a standard log-law for the boundary walls, would ill estimate the concentration/size of the explosive atmosphere by a factor of about 30% which is considerable as far as the potential subsequent explosion effects are concerned ;

 quite a significant amount of structured knowledge (not only a culture of…) in fluid mechanics was required to find the reason for the discrepancy and the latter was called for only after consideration of the experimental results. So the confrontation with existing data is certainly desired but even finding the relevant experiments containing the relevant key physics (with regard to the practical situation to be studied) is very challenging since a deep knowledge of the physics remains compulsory. A practical worrying aspect concerns the boundary layers. Borders are everywhere in practical situations and the present illustration shows that a wrong physical representation of them may lead to severe mismatches. Certainly, many situations corresponding to the example developed in this section are concerned. It calls for a analysis of the PDR representation of obstacles, although not addressed in the present work. A further point pertains to the size of the computational cells and to the "cost" of the simulations. As shown with this example, although geometrically simple, the cell size needs to be small enough to resolve correctly the formation of the layer. The convergence in "mesh" was reached with about 10 7 cells and the CPU time is estimated at 172 hours, the pre-processing takes 1 hour and half and the postprocessing takes 2 hours (for the simulations realized with low Reynolds number k-Ɛ model). The numerical tests suggest that the accuracy is likely to drop extremely fast when "coarsening" the cells size. For instance, a decrease of only 20% in the number of cells (increase of 10% cell size) results in a drop of 5 to 10% in accuracy. So, performing reasonable simulations promises to be long and expensive.

FLAME PROPAGATION AND EXPLOSION

In this domain, the main challenge is to represent correctly the flame dynamics because, from this, depends enormously the pressure loads and final consequences (as a rule of thumb the pressure increase as function of the square of the flame velocity). In large scale explosion CFD codes, two families of combustion models are commonly implemented: EBU types of models and CREBCOM types of models. The analysis presented in the first chapter suggests that none of them can really mimic the flame behaviour. The purpose of the simulations presented hereafter is to investigate their behaviour in representative situations and guess up to which extent the computed results may be trusted in or not. The influence of the geometry particularly that of the obstacles, is investigated using the EDM and the CREBCOM combustion models. The beta flame model, which is an evolution of the CREBCOM model, could not be programmed in MERLIN.

INFLUENCE OF THE CONFINEMENT

EXPERIMENTAL SETUP AND RESULTS

The aim of the experimental work was to investigate the dynamics of an explosion starting inside a vented enclosure and propagating outwards inside the flammable cloud expelled through the vent. This is a very practical situation (125).

The experimental chamber (Figure 3-63) is a 4m 3 parallelepiped (2 m x 2 m x 1 m) with an open vent (square of 700 mm width). The flammable cloud is an homogeneous mixture of hydrogen and air (16.5% H2), ignited on the opposite side from the vent. The experimental conditions were carefully controlled using oxygen meters for the concentration, many pressure gauges inside and outside and fast cameras. The experimental results are presented on pressure data are those recorded by the pressure sensor P2 located near of the vent and by the sensor installed on a support ("lens"), outside, at 2m from the vent. The flame dynamics was extracted from the high speed films (velocity on the axis of the chamber).

To be able to see the flame, normally invisible, the flammable cloud wa with microscopic NH 4 NO (Figure 3-67). The flame front ( with a growing semi cylindrical shape until approaching the vent. At time 120/130 ms, the apex of the flame rushes inside the large vortex formed outside which expands rapidly from time 130 ms to 145 ms. The large velocity increase corresponds exactly to the ignition and burning of the vortex outside. The time period from 0 to 120 ms correspond to the development of the flame inside the chamber. There, the flame acce ignition point) to about 20 m/s at 120 ms. Immediately after, the apex of the flame enters the vent area and accelerates up to 160 m/s. It can be calculated that this velocity is closed to that estimated for th a pressure pulse is noticed corresponding to the external burning of the cloud. This pulse is produced during the rapid and nearly isotropic outward expansion of the vortex. Image processing was applied on the frame taken at 120 ms to highlight the contour of the flame front. The latter (continuous red line) is superimposed on the original frame. The flame front is very elongated but seems to be moderately corrugated. In such experimental conditions, nearly isobaric, the spatial flame speed is the expansion velocity of the flame (laminar burning velocity times the expansion ratio) multiplied by the ratio between the flame area and the vertical cross section of the enclosure. The measured flame speed is 20 m/s and the expansion velocity is 6 m/s so that the ratio between the flame area and the vertical cross section of the enclosure is about 3.3. The average area of the flame front (doted yellow lines) is mostly that of a cone with the apex on the axis of the chamber at 1.75 m from the ignition source and the base of a height of 1.65 m. The flame area /enclosure cross section ratio is then 1.9. The corrugations of the flame area represent 3.3-1.9=1.4 so 1.4 3.3 ⁄ × 100 ≈ 40 percent. In practice then, 60% of the flame acceleration seems to be due to the elongation of the flame, most probably by the flow, and 40% by some "turbulence" corrugating the front. About the scales of the corrugations, a clear large scale appears, typically 0.5 m wide and 0.15 m amplitude and apparently a much smaller scale, typically 6 to 7 cm with a smaller amplitude 3 cm. The large scale corrugations represent an increase of the flame area (as compared to that of the cone) of about 20% (multiplying by factor 1.2) and the small scale indentations about 35 % (factor 1.35). Multiplying all the contributions (elongation x large scale corrugations x small scale indentations), the flame cross section areas is 1.9 x 1.2 x 1.35 = 3.1 which is close to the value deduced from the velocities. The turbulent burning velocity results from the corrugations/indentations and represents a ratio of area of 1.2 x 1.35 = 1.6 superimposed over the "mean flame area" (here the cone). This means that the turbulent burning velocity would amount 1.6 x S lad = 1.6 m/s since the laminar burning velocity of the mixture is about 1 m/s. Note that all along the propagation of the flame in the chamber, the cloud outside in a "bubble" and not (yet) a developed jet. About the picture taken 150 ms after ignition, the pressure pulse is produced by a nearly isotropic expansion of the "bubble" of unburnt gas expelled by the vent during the propagation of the flame inside the chamber.

NUMERICAL SIMULATIONS

Three dimensional cartesian calculations were performed using the Reacting flow solver of MERLIN. The combustion models were implemented in the following way:  The CREBCOM model is implemented in the Euler equations solver as a source term in the energy, mass fraction and reaction progress balance equations. These equations are solved using a 2 nd order in space Roe numerical scheme (upwind with the minmod limiters) and first order explicit in time. The entropy correction factor is set to 0.05.

 The EDM model is used with the URANS equations as a source term in the energy, mass fraction and reaction progress balance equations. These equations, with the standard k-epsilon model, are solved using a 2 nd order in space Roe numerical scheme (upwind with the minmod limiter for the convective part and a centered scheme for the diffusive part) and first order explicit in time. The entropy correction factor is set to 0.05. The followings initial conditions are:

 Initially, an homogeneous hydrogen-air mixture (16% v/v hydrogen in air), at ambient temperature (293.5k) and pressure (1 atm) fills the enclosure. The outside atmosphere is at the same T, P conditions but air only.  Experimentally, the ignition is achieved using a pyrotechnical spark. To simulate this, a hot pocket of burnt gases (2200K) is lacated at the ignition point (opposite side from the vent). This pocket of hot gas occupies a set of cells of a total volume of 140mm 3 , representing the volume of the match (cylinder of 5mm of length and 6mm of diameter);  In the experiment, the mixture is perfectly still with no turbulence. But again the RANS calculation cannot be triggered with no turbulence. The "low turbulence" option was chosen with 1% turbulence intensity: u = 0.01 x U where U is the free stream velocity. The initial turbulent kinetic energy is then k ≈ 10 U . Ɛ should be set, so that the effective turbulence viscosity is similar to the molecular viscosity υ ~0(υ)~10 m 2 /s. Then, Ɛ = C k υ ⁄ ≈ 10 U and the turbulent length scale is l = C k ⁄ υ ⁄ ≈ 10 U ⁄ . Assuming that the free stream velocity would be on the order of the initial flame velocity, typically 1 m/s, it comes that k = Ɛ = 10 m 2 /s 2 .

The balance equation of the reaction is:

2H + O + 3.762N → 2H O + 3.762N
where water is the main product of the reaction. The basic combustion parameters are given in the Table 15  The vent is initially closed (by a fictitious membrane). An opening overpressure of 60 mbar is chosen (corresponding roughly to the sensor P2 measured overpressure after 120 ms of propagation e.g. when the flame reaches the vent ;

 At the solid wall boundaries (other than the vent), a zero pressure gradient is applied and the "law of wall" is used to estimate the turbulence parameters and the tangential component of the velocity ;  At the external boundaries of the domain, continuity conditions are applied. The computational domain is the chamber itself prolongated by six volumes (Figure 3 The domain was meshed using the multiblock decomposition approach of MERLIN. Only structured meshes were used since no significant shock waves are expected. However, the size of the cell differs from one part to another, depending on the expected flow. In Figure 3-69 for example, a finer mesh is used in the surrounding of the ignition point (up to half the volume of the chamber) to capture correctly the front of the flame during its development phase. Outside the chamber (Figure 3-70 and Figure 3-71) a fine mesh is also used to capture correctly the gas bubble formed in front of the vent and its subsequent ignition by the flame rushing out of the chamber. A coarser mesh is used otherwise around the chamber to act as "buffer zones" for the pressure/velocity waves (avoids artefacts coming from the outer boundaries. About 4892160 nodes are used including 2600000 in the chamber. This meshing strategy was obtained by trial and error.

With CREBCOM model, the cells need to be cubic otherwise the flame will be distorted. A non-homogeneous mesh decelerates the flame in the direction of the elongation of the mesh, because the reaction rate is inversely proportional to the cell size. Because of this, a regular mesh (all cells cubic and identical) was chosen for all the simulations below (4500000 cells). Significant differences appear also about the external combustion. A turbulent jet explosion is predicted by the EDM model whereas an isotropic outward expansion of a gas bubble is observed. Somewhat surprisingly the wright order of magnitude of the maximum overpressure are observed although the physics is clearly not respected. Consequently, the flame seems much more turbulent and less elongated than in the reality. The last point is about the nature of the turbulence on the flame front inside the chamber. If the small scale indentations may share some similarities with the calculated integral length scale of the turbulent flowfield (cms), the corrugations are much larger by an order of magnitude and cannot be due to the turbulence of the flow. From the analysis presented in chapter 1, both the elongation of the flame and the large scale structurations could be explained by the flame natural instabilities which are ignored in EDM models.

To run the CREBCOM algorithm, C g (also called K 0 ) has to be defined. It can be related to the turbulent burning velocity, which is proportional to the laminar burning velocity in the "flamelets" regime of the Borghi diagram. Therefore, the laminar burning velocity S L has to be defined. In this work, the laminar burning velocity evolution with the mixture ratio was given by CNRS Orleans (126). This expression takes into account the mixture composition and the initial temperature in the range of values relevant for severe accident conditions. This expression established, supposing as initial temperature 298k and initial pressure 100kPa, is:

S (m/s) = (-0.29 + 1.07ϕ + 1.44ϕ ) × (1 -χ ) [187]
where ϕ is the equivalence ratio and χ is the diluent molar fraction in the mixture.

A relationship was proposed between C g and the laminar burning velocity S for laminar flame expanding in a closed vessel by comparing the experiments to the predictions of TONUS. The expression found is:

C (m/s) = -0.0594 + 1.581(S ) + 0.4799(S ) .

[188] Another parameter to be set is the ignition threshold F which is directly linked to the progress variable c. Its choice is not straight-forward. If F is too small, the flame will spread rapidly through the volume criterion because of the numerical diffusion which is likely to "diffuse" c. Conversely, if F is too high, the flame may have difficulties to propagate. F was chosen equal to 0.5 which is in line with the practise of CREBCOM models 14 .

The numerical results are presented on Figure 3-83 and Figure 3-84 together with the experimental data. The different values of C g selected give approximately the same results and give approximately the wright order of magnitude of the maxima. Nevertheless, the experimental profiles stand out from those obtained numerically. Experimentally, the flame is laminar in the chamber while the control parameter used to approximate it numerically corresponds to turbulent propagation. As a result, the flame dynamic found experimentally in the chamber is different from the numerical prediction. Outside, a pressure peak is calculated corresponding presumably to a fast development of the flame front area rather than to an increase of the local turbulence and turbulence burning rate since no turbulence is incorporated into the CREBCOM algorithm.

Again, correct orders of magnitudes were obtained with the wrong physics. But contrary to the EDM model, some tuning (of C g ) was required. This tuning is also dependent on the cell size as shown in the first chapter.

INFLUENCE OF OBSTRUCTIONS

In most practical situations, the flame front propagates in areas where obstacles are present. It is often impossible to describe accurately all the geometry and to "mesh" it correctly. In current CFD practice devoted to large scale explosions, the overall geometry (walls, openings, partitions,…) is resolved but the obstructions smaller than the cell size are "modeled" using the PDR technique. With this technique, it is assumed that the role of the obstacles inside the cells is to produce (shear and wake) turbulence, further accelerating the flame. But obstacles do also modify the velocity field which may impact significantly the flame via the instability mechanisms. Experiments and simulations were performed to investigate this aspect.

Again, the PDR method is not used, the obstacles are introduced explicitly.

EXPERIMENTS

In the following experiments (127), the incidence of the confinement, as illustrated in the previous section, was voluntarily limited. The experimental setup is a 3D geometry with a corner ignition. The corner is formed with two vertical steel walls (length 3 m and height 1 m), covered with a flat (transparent) roof (Figure 3-87).

The two other vertical faces are sheltered using a thin plastic sheet which is simply supported so that when the flame propagates inside the chamber, they are lifted without any mechanical resistance. The inner volume is 9 m 3 . Homogeneous hydrogen-air mixtures were used and a spark is ignited in the corner formed by the two steel walls. The mixture was made visible using the same seeding technique as in the previous section and the experiments were filmed from above. Three pressure gauges were used. The obstruction consisted in vertical cylinders of different diameters and variable number.

In the specific experiment investigated here, the obstruction consisted in two rows of three cylindrical obstacles of 20 cm in diameter separated by 80 cm. The mixture contains 14.5% of hydrogen in air and is initially at rest. The results are represented on Figure 3-88 for P2 pressure sensor and for the flame velocity extracted from the video. The following interpretation of the data can be given: from time 0 to 0.1 s, the flame front develops spherically around the ignition point and the peak corresponds to the instant when the top and bottom part of the flame reach the walls. The subsequent extinction of those parts of the flame results in a decrease of the pressure. When the flame has achieved its final transformation into a quarter of a cylinder propagating outwards, the pressure stabilizes more or less (from time 0.13 s to 0.25 s) before interacting with the obstacles (first row from 0.25 s to 0.3 s and second row from 0.3 s to 0.38 s). The flame velocity was measured using the video along a straight line going at equal distances between the obstacles. To some extent, it is a mean flame velocity. The data can be correctly interpreted up to time 0.27s. At that time, the flame is between the two rows of obstacles. Before the interaction with the first row, the flame velocity is nearly constant : about 8 m/s. Somewhat surprisingly, the flame velocity increases (from 8 to 13 m/s) ahead of the obstacles (from 0.21 to 0.23 s) when the edge of the flame is still about ½ obstacle ahead. Then, the flame velocity reaches about 17 m/s at time 0.27 s. It can be shown that the oscillations in the flame speed are due to the acoustics aof the chamber not to the obstacles. Thus, the influence of the obstacles is to increase moderately the (average) flame velocity but the link with the rapid evolution of the pressure signal is not immediate. To understand better, the evolution of the flow/flame front needs to be investigated.

When the flame approaches the boundary layer of the obstacles (from time 0.21 to 0.23 s), the flow is decelerated in the boundary layers and accelerated in between the obstacles. Due to the obstruction offered by the obstacles (including the boundary layers), the flow (and thus the flame) should be accelerated by about 50% which is effectively observed. But in the meantime, the pressure decreases! The pictures indicate the flame front just wraps around the front face of the obstacle and may be extinguishes leading to a reduction of the flame front area. On the contrary, the maximum of the pressure may correspond to the reconstruction of the flame behind the obstacles leading to a rapid increase of the flame area. So, the pressure variations may rather represent the variations of the flame area and while the average flame velocity varies as function of the necking of the flow and reactivity of the mixture. The first and the second peak are due to the two successive rows of obstacles but the third one is produced when the flame rushes out (additional turbulence due to the venting of the unburnt mixture). Note also that immediately downstream of each obstacle, there is no jump in the flame velocity as would suggest a local increase of the turbulence in the wake of the obstacles.

SIMULATIONS

The same simulation strategy than in the preceding section was adopted. Differences are the reactivity of the mixture, boundary conditions, cells sizes and incorporation of the obstacles.

Only the EDM model was used because the turbulence created by the obstacle is cannot be easily incorporated in a Euler type of solver as used with CREBCOM.

The real configuration is represented in a 3D geometry, which seems reasonable.

Only the interior of the setup (size 10 m × 10 m × 10 m

At the solid walls, the "log law" is used to estimate the turbulence parameters and the tangential velocity component. The "plastic sheet" walls held at the atmospheric pressure, considered as "an open boundary"

The ignition source is match).

On Figure 3-89 are presented mesh refinments. The convergence is reached with mesh1.

Analysis of some modelling aspects in (industrial) explosion simulation

Average flame velocity (between the obstacles) and explosion overpressure as function of time (14.5 % H2 in air, 2 row of DN 200 cylinders, the first row at 1.5 m and the second 0.8 m further)

The same simulation strategy than in the preceding section was adopted. Differences are the reactivity of the mixture, boundary conditions, cells sizes and incorporation of the obstacles.

model was used because the turbulence created by the obstacle is incorporated in a Euler type of solver as used with CREBCOM.

The real configuration is represented in a 3D geometry, which seems reasonable.

Only the interior of the setup is accounted for using about 9000000 regular cells m) in the 3 m x 3 m x 1 m geometry.

At the solid walls, the "log law" is used to estimate the turbulence parameters and the tangential velocity component. The "plastic sheet" walls are assumed to be held at the atmospheric pressure, considered as "an open boundary"

The ignition source is modelled as in the preceding test case (pyrotechnical are presented the overpressure profiles obtained with different The convergence is reached with mesh1. The same simulation strategy than in the preceding section was adopted. Differences are the reactivity of the mixture, boundary conditions, cells sizes and model was used because the turbulence created by the obstacle is incorporated in a Euler type of solver as used with CREBCOM.

The real configuration is represented in a 3D geometry, which seems reasonable. 9000000 regular cells At the solid walls, the "log law" is used to estimate the turbulence parameters and are assumed to be held at the atmospheric pressure, considered as "an open boundary" (continuity). modelled as in the preceding test case (pyrotechnical the overpressure profiles obtained with different As in the preceding configuration, the pressure traces clearly different. In the early time (0 ignition point but close to all the solids walls the "log the turbulence which leads to an overestimate of both the flame velocity and pressure. Note that ta large rather constant velocity propagation phase (at 12 m/s) so that the edge of the flame front is at 2 m from the ignition point when just in the wake of the first row of obstacles.

On Figure 3-92 is shown the flow velocity ahead of the flame. The difference between this and the flame velocity ( velocity. Clearly, the turbulent burning velocity jumps from about 2 m/s, before the flame/obstacle interaction to about 10 m/s immediately after. So turbulence produced by the wake of the obstacles wh velocity increase. This does not match with the experimental observation.

Analysis of some modelling aspects in (industrial) explosion simulation

Flame dynamic during the during the explosion of a 14.5% hydrogen air mixture around obstacles using the EDM model

As in the preceding configuration, the pressure traces are resembling although In the early time (0-0.1 s), the flame develops radially around the ignition point but close to all the solids walls the "log-law" and the k the turbulence which leads to an overestimate of both the flame velocity and a large flame velocity increase occurs at time 0.17 s after a rather constant velocity propagation phase (at 12 m/s) so that the edge of the flame front is at 2 m from the ignition point when this increase occurs first row of obstacles.

is shown the flow velocity ahead of the flame. The difference between this and the flame velocity (Figure 3-91) gives the local turbulent burning velocity. Clearly, the turbulent burning velocity jumps from about 2 m/s, before the flame/obstacle interaction to about 10 m/s immediately after. So turbulence produced by the wake of the obstacles which is responsible velocity increase. This does not match with the experimental observation.

ndustrial) explosion simulation

Page 193 sur 225 of a 14.5% hydrogenmodel resembling although 0.1 s), the flame develops radially around the law" and the k-ε model boosts the turbulence which leads to an overestimate of both the flame velocity and flame velocity increase occurs at time 0.17 s after a rather constant velocity propagation phase (at 12 m/s) so that the edge of the this increase occurs, meaning is shown the flow velocity ahead of the flame. The difference gives the local turbulent burning velocity. Clearly, the turbulent burning velocity jumps from about 2 m/s, before the flame/obstacle interaction to about 10 m/s immediately after. So, it is the ich is responsible for the velocity increase. This does not match with the experimental observation. The EDM model is based on the assumption that the flame propagation is driven by the turbulence of the flow. In this particular example, the flame is not initially turbulent but the flow velocity the Reynolds number over the obstacles is about 10 been a favorable configura simulations show strong combustion accele not in line at all with the experimental observations.

The experiment do not exhibit any significant wake behind the obstacles in which the flame would rush. The flame velocity varies only slightly (due to the nec the flow because of the reduced cross section in the region of the obstacles) with no specific sign of variations commensurate to the pressure pulses. It is believed that the pressure pulse is due to the local flame area variations (elongation, reduction) while going around the obstacles. Why the turbulence does not seem active? The same values of u' and l estimated: respectively 0.05 cm and 1 m/s. This means that the turbulent diffusivity μ t /ρ is 0.05 m 2 /s. The duration of the explosion is about 0.3s so that the thickness of the turbulent boundary/shear layers (the turbulence is created in boundary/shear layers) is about (0.05 x 0.3) a very limited extent and this can expl not occur in the wake of the obstacles Analysis of some modelling aspects in (industrial) explosion simulation

Flow velocity along the axis ahead of the flame as function of time during the explosion of a 14.5% hydrogen-air mixture around obstacles using the EDM model model is based on the assumption that the flame propagation is driven by the turbulence of the flow. In this particular example, the flame is not initially velocity is large enough to become turbulent: for instance, the Reynolds number over the obstacles is about 10 5 ... A priori, it should have configuration for the considered models. As expected, the simulations show strong combustion accelerations due to the obstacles. not in line at all with the experimental observations. he experiment do not exhibit any significant wake behind the obstacles in which the flame would rush. The flame velocity varies only slightly (due to the nec the flow because of the reduced cross section in the region of the obstacles) with no specific sign of variations commensurate to the pressure pulses. It is believed that the pressure pulse is due to the local flame area variations (elongation, uction) while going around the obstacles. Why the turbulence does not seem active? The same values of u' and l T than in the previous "test case" were estimated: respectively 0.05 cm and 1 m/s. This means that the turbulent diffusivity e duration of the explosion is about 0.3s so that the thickness of the turbulent boundary/shear layers (the turbulence is created in boundary/shear layers) is about (0.05 x 0.3) 1/2 =0.1 m. So, turbulence exists but to a very limited extent and this can explain why the acceleration of the falme does not occur in the wake of the obstacles. model is based on the assumption that the flame propagation is driven by the turbulence of the flow. In this particular example, the flame is not initially is large enough to become turbulent: for instance, ... A priori, it should have the considered models. As expected, the due to the obstacles. But this is he experiment do not exhibit any significant wake behind the obstacles in which the flame would rush. The flame velocity varies only slightly (due to the necking of the flow because of the reduced cross section in the region of the obstacles) with no specific sign of variations commensurate to the pressure pulses. It is believed that the pressure pulse is due to the local flame area variations (elongation, uction) while going around the obstacles. Why the turbulence does not seem than in the previous "test case" were estimated: respectively 0.05 cm and 1 m/s. This means that the turbulent diffusivity e duration of the explosion is about 0.3s so that the thickness of the turbulent boundary/shear layers (the turbulence is created in turbulence exists but to why the acceleration of the falme does Potentially a low Reynolds k-ε model would have been able to capture this transient aspect of the turbulence flow field but a much smaller cells size is needed which is hardly feasible. The turbulence in flows impacting an obsctacles, producing some rotation (like over an obstacle/cylinder) is overestimated with kthe standard epsilon model since this model is not able to distinguish between shear and stress. A solution to this problem was once proposed by Kato and Launder (128) which suggests to replace the production term of the standard kepsilon transport equation by the following:

P = μ ∂u ∂x + ∂u ∂x ∂u ∂x - ∂u ∂x .
[189] It is more a numerical trick than a fundamental physical reasoning. Simply, this modified production term avoids excessive production of turbulent energy in stagnation regions, shock and regions with very large normal strain.

IMPLICATIONS

As suspected none of the proposed categories of combustion modelling is convenient since they ignore the physical reality of flame propagation. Nonetheless, those models may produce correct order of magnitudes of the overpressure although for wrong reasons. A reason is that the constants of the models were "tuned" on experimental results. Clearly, the EDM model seems more reliable since the prediction should not depend on the cell size and mesh topology contrary to the CREBCOM model. Because of their robustness and possibilities to capture correct orders of magnitudes, EBU categories of models have been in use for a long time and were continuously refined. A major deficiency of EBU models is the basic assumption that the turbulent flame is passive against the turbulence of the flow. As explained in the first chapter and as illustrated with the experimental configurations above, it is not true and other aspects of flame dynamics should be incorporated.

Thus, with the versions of the combustion models implemented in large scale industrial explosion CFD codes, the predictability of the simulations cannot be assessed because the underlying physical assumptions are wrong.

For all these reasons, significant developments are still required.

This PhD work originates from the scientific and technical investigations that were decided immediately after the Buncefield explosion (2005). The traditional explosion safety engineering methods failed, and alternative methods were pushed forward and especially existing Computational Fluid Dynamic Softwares dedicated to large scale explosion scenarios in industrial environments.

But can they really do better? In 2002 Health and Safety Laboratory raises some doubts. The reason for this lies perhaps in the rather poor outcomes of benchmarking exercises showing that the discrepancies between different simulations on the same cases can easily be in a factor of 4. The last benchmarking exercise was not better (129).

Nevertheless, CFD methods should do better than simple engineering tools. The objective of this thesis is mostly to find out the reason why the discrepancies are so large.

The scope is restricted to the numerical tools that can be used to simulate the propagation of the flame and the subsequent pressure effects in a large industrial environment. The size is the first challenge, typically tens of meters, and the complexity is the second one, with obstructions, confinement, transmission from a building to outside… Only a limited number of codes can do this (FLACS, EXSIM, CFX,…).

An analysis of the physics of explosion phenomena is presented in the first chapter. The discussion is centered on the dynamics of the flame. In a real environment the flame front will be submitted to many disturbances: turbulence, acoustics, acceleration/deceleration of the flow… and all of them will impact the combustion, thus not only the turbulence of the flow. The effect of the disturbances is intimately linked to the unstable nature of any flame front. In other words, the flame front is not passive against the disturbances but either filters out or amplifies them. As shown in the same chapter, all the "industrial explosion dedicated CFD softwares" use the Reynolds Average Navier Stokes technique with a k-epsilon model to calculate the flow and the turbulence. The flame is propagated using Eddy Break Up (EBU) type models or a CREBCOM type algorithm. FLACS uses the flame model, which is a CREBCOM type of model called the Beta flame model. In CREBCOM type of models, the combustion velocity is deduced from a correlation using the calculated parameters of the flow. In EBU types of models, the combustion is directly integrated in the equations. The underlying idea is of all these models is that the flame is primarily dominated by the turbulence of the flow. A discussion is proposed about the various numerical strategies to solve the RANS equations since they may have a direct impact on the results.

In the second chapter, the author explains why he chose to develop his own CFD platform rather that using as they stand the above mentioned CFD codes. The first objective was to gain in control and flexibility (have the possibility to associate as desired various numerical strategies and models). The second objective was to isolate better the "user dimension" which is known to be a main reason of "bad performances" of numerical simulations. MERLIN was thus developed during this PhD work to test the various CFD models and solvers employed in the large scale (industrial) explosion softwares. MERLIN is a CFD toolbox enabling the user to bridge openly various numerical methods (numerical schemes, mesh topology,…), systems of equations (Euler, RANS) and models (k-epsilon, low Reynolds kepsilon, EDM, CREBCOM,…). Illustrations of the capabilities of MERLIN are given showing that the toolbox may be used for other applications than safety and an original way of verifying the software is presented.

The last chapter is dedicated to the analysis of the performances of the "classical" solvers 15 used in computing large scale (industrial) explosions. When attempting to simulate an explosion within an industrial context, the first difficulty is to be able to simulate correctly a relatively large spectrum of interlinked physical processes. It is not enough to be capable of propagating a flame inside a flammable cloud. It is necessary first to simulate the formation of the cloud since this "initial" step of the scenario of the accident has been for long been recognized as a leading aspect. In some cases, it might also be necessary to calculate the propagation of the pressure wave away from the exploding cloud. And again, the geometry is usually very large and complex. The behaviour of the solvers were investigated on each step of an explosion: the expansion of a gaseous release, the formation of the plume, the accumulation in a confinement, the explosion inside and outside and the subsequent propagation of a pressure wave. The main conclusions are:  Since it is numerically difficult and resource demanding, the details of the "expansion zone" of the jet produced by a massive leakage of a pressurised combustible gas is skipped and replaced by a fictifious leakage source in the considered category of CFD codes. This is the "notional nozzle model".

It is shown that this modelling strategy may be wrong especially when the leakage is strongly underexpanded. A drift up to 50% of the flow velocity and mixing ratio is possible;

 Provided the expansion zone is correctly accounted for, the standard kepsilon model performs very well for the rest of the jet confirmed this approach is valid in shear dominated flows. Indirect observations of the flow around obstacles nevertheless highlights limits for rotating flows in which too much turbulence is produced. This deficiency is linked to the Boussinecq model. The PDR modelling of the obstacles may be a solution but could not be tested during this PhD work. The standard k-epsilon model does not perform well in boundary layer near solid wall and a special law of wall is tradionnally implemented. Implicitly, it is admitted that a steady state turbulence flowfield is reached. This assumption forces the turbulence. In highly transient situations, like during the accumulation of an explosive atmosphere in a building, the characteristics of the cloud may be ill estimated by up to 30%. This deficiency may be totally removed using a refined k-epsilon model like the "low Reynolds number" k-epsilon;  The combustion models used in the considered CFD codes simulate a physics which is not at work in the explosion configurations. Despite this, correct orders of magnitude of the maximum explosion overpressure may be obtained (by comparison with experiments) although the dynamics of the flame propagation may be very different. An explanation is that the constants of the models were tuned on experiments. There is no guarantee that in configurations totally different from the latter experiments this minimum level agreement could be kept. The conceptors of the large scale (industrial) explosion codes seem aware of this and continuously conduct calibration exercises;  Lastly the propagation of the pressure wave in the envrironment and in particular in an obstructed area was considered. First the influence of the mesh topology is important. Structured meshes should be avoided otherwise significant numerical artefacts will appear like wrong reflection patterns. Numerical schemes are often Riemann approximative schemes, like the Roe scheme. The latter is fast and robust, thus very appealing for the category of CFD software considered. Unfortunately, it is shown that the approximate nature of such schemes modifies the pressure waves so that unexisting reflections are predicted. Following the predicted pressure field in an obstructed configuration (many obstacles), is unreliable. It is believed that this trend is more pronounced the stronger the incident pressure wave.

As these tools stands presently, a very strong link should be made between the simulations and the experimental results, especially to calibrate the combustion models. The CFD software would then be more an interpolation tool with limited extrapolation (predictive) ability. This is what is done so far. Nevertheless, the key point is to be able of selecting the relevant test cases. The experimental situations need to be chosen as being representative of the physics at work.

Note that even for an expert in fluid mechanics and explosion phenomena, this is not an obvious exercise, and this puts in perspective the "user" dimension.

In the specific context of this study, a typical user will be a safety engineer, trained into the main physical aspects of explosion scenarios, but having only a general knowledge of the use of a CFD software. It means that even though he will be able to build a scenario of accident incorporating details of the conditions and to elaborate on the consequences based on engineered safety tools (using for example, the TNO guides), he will not claim to be a specialist in fluids mechanics, nor in combustion. In CFD, his culture may even be more superficial. He will know enough to understand that a mesh is required and that testing the "convergence in mesh" is necessary. Consequently, many of the numerical and mathematical aspects of CFD will stand out of his scope to a significant extent. Unfortunately, the key aspects addressed above may have a significant impact on the simulations. Will he be able to identify the situations where the standard k-epsilon model is not adequate ? Will he understand enough of the explosion phenomena to choose in the available experimental database the tests relevant to his configuration to "tune" conveniently the combustion model ?

It seems that further developments of the considered softwares would be required to release the limitations that makes a "safe use" these codes out of reach to "standard user". It seems that the URANS approach, although being approximate by certain aspects, is globally satisfactory, because containing the appropriate information (large scales), offers its robustness and makes it possible to compute in 2 dimensions. A correct representation of the boundary layers and of obstacles is certainly needed but the key point is undoubtedly the combustion model. Ideally, the flame should be modeled as it stands: a thin interface corrugated "thanks to" the perturbations. Direct numerical simulation supports this point of view (130) and modern combustion research is increasingly acknowledging this.

Other modelling strategies than URANS exists that the "Large Eddy Simulation". The applicability to routine simulation of large scale explosion in a realistic environment need to be demonstrated. This modelling strategy shares with URANS the difficulty in modelling the combustion although offering more opportunities to succeed in this area.

  Angle between the direction of the reflected shock and the axis of the jet [Specific heat capacity at constant volume [J/kg.K] C ε : Constant of the k-epsilon model [-] The standard enthalpy of formation of the fuel [J/kg] h _ : The standard enthalpy of formation of the species i [J/kg] h _ : The standard enthalpy of formation of O [J/kg] h _ : The standard enthalpy of formation of the product [J/kg] Distance y to the wall [-] y : Mass fraction of the fuel in the unburnt mixture [kg species/kg]

Figure 0 - 1 :Figure 0 - 3 :Figure 0 - 4 :

 010304 Figure 0-1: Intercomparison of maximum explosion overpressure inside an offshore module

Figure 0 - 5 :

 05 Figure 0-5: Test facility (ICT -ball of 2000 m 3 full of a stoichiometric mixture of H 2air at rest -ground ignition on the axis of symmetry).

6 :

 6 discrepancies in the levels and shape of the signals appear (Figure 0-6 a). The maximum values of the overpressure for instance vary from a code to another by a factor of 2 (Figure 0-6 b) and most of the simulated signals are very different from the measurements. So clearly, significant modelling issues of the physics remain. Signals of pressure at 2 m from the boot (CAST3M-ECA, COM3D-FZK, FLACS-GEXCON, REACFLOW-JRC, AUTOREAGAS-TNO, FLUENT-UU) (a) and report maximum simulated/ measured values in 2 m, 5 m, 8 m, 18 m, 35 m and 80 m in the boot (b).

Figure 1 - 1 :

 11 Figure 1-1: Unconfined (a) and confined (b) explosion

Figure 1 - 2 :

 12 Figure 1-2: Description of a combustion wave using temperature

Figure 1 - 3 :

 13 Figure 1-3: Deviation of the current lines through a flame front

Figure 1 - 4 :

 14 Figure 1-4: Development of turbulence in a flow past a sphere (from H. Werlé (36))

Figure 1 - 5 :

 15 Figure 1-5: Typical point velocity measurement in turbulent flow

Figure 1 -

 1 7 is an example.

Figure 1 - 7 :

 17 Figure 1-7: Representative diagram of combustion regime

  Figure

Figure 1 - 8 :

 18 Figure 1-8: Examples of mechanical calculators

Figure 1 - 9 :

 19 Figure 1-9: A control volume with a cylindrical obstacle inside. The gray area shows the volume occupied by the obstacle, 𝑉 .

Figure 1 - 10 :

 110 Figure 1-10: The surface area 𝐴 of a control volume, partially occupied by an obstacle. The bright area shows the area occupied, 𝐴 .

Figure 1 - 11 :

 111 Figure 1-11: Example of calculation domain in one dimension

Figure 1 - 12 :

 112 Figure 1-12: Upwind and centered evaluation scheme

  ) = max (0, min(1,2r) , min (2, r)) Chakravarthy Ψ(r) = max(0, min(β, r)) ; 1 ≤ β ≤ 2

Figure

  Figure At the left: Generation of a block nodes and edges of the block on the geometry, at the centre: structured mesh generated, at the right: unstructured mesh obtained the by re of the quadrangles of the structured mesh.

Figure 1 - 13 :

 113 Figure 1-13: Generation of the structured mesh: At the left: Generation of a block around of the initial geometry and projection of nodes and edges of the block on the geometry, at the centre: structured mesh generated, at the right: unstructured mesh obtained the by re of the quadrangles of the structured mesh.

  The distortion expected of triangular elements by measuring the gap of angles generated compared to angles multiple of  The length of the loops.  The adaptation error corresponding to of nodes to prescribe on the the ideal number of nodes estimated with respect to the initial distribution Nodes are distributed on the cutting line as function of the length of the first and last segment on this line. Those latter lengths are the average of the segments connected to first and final nodes.

Figure 1 - 14 :

 114 Figure 1-14: Procedure of generation of an unstructured mesh: At left: Generation of the initial loop and distribution of edge nodes, At centre: Mesh intermediate obtained during the closure process, At right

  Procedure of generation of an unstructured mesh: : Generation of the initial loop and distribution of edge nodes, : Mesh intermediate obtained during the closure process, At right: Final mesh obtained after smoothing.

  The distortion expected of triangular elements by measuring the gap of introduced by favoring the length of small between the number to respect certain criteria like e domain, the latter being on the loop.Nodes are distributed on the cutting line as function of the length of the first and last segment on this line. Those latter lengths are the average of the sizes of the Procedure of generation of an unstructured mesh: : Generation of the initial loop and distribution of edge nodes, : Mesh intermediate obtained during the closure process,

  ρ • y. τ μ ⁄ . The parameter u t is the well known shear velocity ( u τ = τ w ρ ). T w is the shear stress in the bulk flow at the limit of the boundary layer chosen by the user (30 < y < 500). The turbulence parameters in the boundary are k = τ μ ; ε = τ κ . For the other parameter, Dirichlet conditions (normal velocity=0, Temperature fixed…) can be used of Neumann type (zero pressure gradient).

Figure 2 - 1 :

 21 Figure 2-1: An overview of the folder MERLIN

Figure 2 - 2 :

 22 Figure 2-2: an overview of the content of UserInput folder dedicated to the simulation of the compressible non-reacting fully three dimensional case represented by the matlab file "Fully3D".

Figure 2 - 3 :

 23 Figure 2-3: Fully-3D simulation

Figure 2 - 4 :

 24 Figure 2-4: example of data files obtained during a simulation for a postprocessing with realized with tecplot

Figure 2 - 5 :

 25 Figure 2-5: an overview of the preprocessing solver of MERLIN

Figure 2 - 6 :

 26 Figure 2-6: computational domain example

Figure 2 - 7 :

 27 Figure 2-7: W component hat function given on outflow face

Figure 2 - 8 :Figure 2 - 9 :

 2829 Figure 2-8: staircase approximation of user defined domain boundary

Figure 2 - 10 :

 210 Figure 2-10: An overview of the Postprocessing solver

Figure 2 - 11 :

 211 Figure 2-11: An overview of the Euler equations solvers

Figure 2 - 12 :Figure 2 -

 2122 Figure 2-12: An overview of the non reacting RANS flow solvers

Figure 2 - 13 :

 213 Figure 2-13: An overview of the reacting RANS flow solvers

  -14 and Figure2-15). The height of the manifold is 2 m and the cross section is 0.5 x 0.5 m 2 . A 42600 cells structured mesh is used. There are two (Dirichlet) circular outflow orifices (0.15 m diameter) with plug flow velocity profiles and plug flow inflow (0.15 m diameter) at the basis of the manifold (flow velocity is 0.01 m/s). All other faces are solid walls with log laws (for turbulence parameters). The URANS equations with the standard k-epsilon model are used to simulate this phenomenon. A second order Roe solver (one order in time with an explicit resolution) is sued with a minmod limiter. The diffusive terms are discretized using the central differencing scheme.

Figure 2 -

 2 Figure 2-14: post representation of Dirichlet velocity boundary conditions

Figure 2 - 15 :

 215 Figure 2-15: Post-processing of a fully 3D Compressible Navier-Stokes simulation in a manifold of 2m height, 0.5 x 0.5 m 2 as cross Section and Re=32 (The diameters of the inlet/outlets are 0.15 m.). The lines in black are streamlines.

Figure 2 -

 2 Figure 2-16: zoom of stagnation point

  -18 at different times after the start of the propagation. The vertical coordinate is the overpressure in Pa at the ground level and the horizontal axis are the distances in meter. t = 2 ms t = 6 ms t = 10 ms Figure 2-18: Two-dimensional blast propagation; pressure field For the 3D case, a slightly larger domain was chosen (25 m x 25 m x 25 m) with the same initial conditions, boundaries and numerical scheme. The charge (used to model the blast source) is positioned at a corner of the block and extended over 64 cells. As previously, a 15625 cells structured and regular mesh is used. The overpressure level is the colored scale ( Figure 2-19). Clearly the results are very similar to the preceeding case suggesting that, in certain situations, a 2D simulation may be as accurate as a 3D while being much less expensive and long. t = 0 ms t = 5 ms t = 10 ms Figure 2-19: three-dimensional blast propagation; pressure field.

A

  cartesian plane 2D geometry was chosen using the compressible URANS equations, H-M combustion model (EDM), standard compressible k-epsilon model and second order Roe numerical scheme with a minmod limiter. The discretization of the equations is performed similarly to that of the compressible flow inside a manifold case. Plug flow velocity Dirichlet boundary conditions are given for the lateral air co-flow whereas Poiseuille Dirichlet boundary conditions are assigned for the central methane injector. On the outflow face at Z=0.8, compatibility boundary conditions are assigned for all the transported species. On the lateral walls, homogenous Dirichlet boundary conditions are set for the mixture fraction. A general view of the results is presented in the Figure 2-20.

Figure 2 - 20 :

 220 Figure 2-20: Contour of density of a methane air diffusion flame

The 3 -

 3 D meshes produced are presented in Figure 2-21 to Figure 2-23 and can be described as follows:

Figure 2 -

 2 Figure 2-21: 3D structured meshes: a) Cartesian and b) skewed curvilinear

Figure 2 - 22 :

 222 Figure 2-22: 3D unstructured meshes: a) unstructured mesh with tetrahedral cells, unstructured mesh with prismatic cells, c) highly skewed unstructured mesh with tetrahedral cells, and d) highly skewed unstructured mesh with prismatic cells

Figure 2 -

 2 Figure 2-23: 3D hybrid meshes: a) skewed hybrid and b) highly skewed hybrid

Figure 2 - 24 :

 224 Figure 2-24: Order of accuracy results for URANS equations on a 3D skewed hybrid mesh using a) L2 norm of the discretization error and b) L∞ norm of the discretization error

Figure 2 - 25 :

 225 Figure 2-25: Order of accuracy results on the 3D skewed hybrid mesh for k-ε turbulence model using a) L2 norm of the discretization error and b) L∞ norm of the discretization error

Figure 2 -Figure 2 - 27 :

 2227 Figure 2-26: Order of accuracy results on the 3D highly skewed curvilinear (i.e., structured) mesh with hexahedral cells for k-ε turbulence model using a) L2 norm of the discretization error and b) L∞ norm of the discretization error

  -28).

Figure 2 - 28 :

 228 Figure 2-28: Order of accuracy results for time accuracy of the unsteady flows on the 3D hybrid grid using a) L2 norm of the discretization error and b) L∞ norm of the discretization error

Figure 3 - 1 :

 31 Figure 3-1: Typical industrial geometry (Petroleum Company LP refinery at Texas City, Texas)

Figure 3 - 2 :

 32 Figure 3-2: Schematic representation of a shock tube at initial time (t=0) and waves that are propagated in the tube at time t > 0.

Figure 3 - 3 :Figure 3 - 4 :

 3334 Figure 3-3: Computation domain and boundary conditions

Figure 3 - 5 :

 35 Figure 3-5: Details of the preceding figure: head of the expansion (left) and foot of the upstream shock (right).

Figure 3 - 6 :

 36 Figure 3-6: Influence of the mesh refinement : Roe 2 nd order coupled with the "Minmod" limiter: head of the expansion (left) and foot of the upstream shock (right)..

Figure 3 - 8 :

 38 Figure 3-8: Schema of pseudo stationary reflection patterns of an oblique shock wave

Figure 3 - 9 :

 39 Figure 3-9: First initial configuration similar to shock tube problem

Figure 3 - 10 :

 310 Figure 3-10: Second initial condition in the field containing the ramp

Figure 3 - 11 :

 311 Figure 3-11: Computation domain and boundary conditions

Figure 3 -

 3 Figure 3-12: isocontours field of density simulated by Takayama and Jiang for =46°. Result taken from K.Takayama and Z.Jiang (255 x 256 cells -65280 nodes-unstructured).

Figure 3 - 13 :

 313 Figure 3-13: Numerical interferogram (left) realized by Takayama and Jiang and experimental shadowgram (right) realized by J.M Dewey in the same conditions for =46°. Results taken from K.Takayama and Z.Jiang (1997).

Figure 3 - 14 :

 314 Figure 3-14: Interaction of an unsteady shock with a ramp: isocontours field of density for the solutions obtained on a structured / divided mesh type with the different numerical schemes chosen.

Figure 3 - 15 :

 315 Figure 3-15: Topology of the carbuncle solution, temperature (left), pressure (right) (107).

Figure 3 - 16 :

 316 Figure 3-16: Supersonic flow (Mach number = 10) over a blunt body: mesh and the isocontour of the longitudinal velocity field obtained with the Roe solver without entropy correction.

Figure 3 - 17 :

 317 Figure 3-17: Interaction of an unsteady shock with a ramp: isocontours field of density for the solutions obtained on an unstructured mesh with two level of mesh refinments and the van Leer numerical scheme (continuous line), the symbols are the results of E. Ritzerfeld, H. Kleine and H. Gronig obtained with a Roe solver.

Figure 3 - 18 :Figure 3 - 19 :

 318319 Figure 3-18: Interferogram and experimental isopycnics taken from H.M.Glaz,P.Colella, I.I.Glass and R.L.Deschambault (105)

  -20), 1m high and 3 m long, an 0.2 m upward step is inserted at 0.6 m from the inlet (left). A supersonic flow of ambient air (1 bar abs, 298 K at 1048 m/s) is admitted from the left and flow towards the right-hand side.

Figure 3 - 20 :Figure 3 - 21 :

 320321 Figure 3-20: Computation domain dimensions

Figure 3 - 22 :

 322 Figure 3-22: Computation domain and boundary conditions used

  1) was applied. The AMA technique was employed and the final mesh contains 11591 nodes. The isocontour fields of density, pressure and longitudinal velocity obtained at t*=1 and t*=1.5 are shown in Figure 3-23 and Figure 3-24 respectively.

Figure 3 - 23 :

 323 Figure 3-23: Isocontours field of density for the Woodward and Colella test case at time t*=1.5 (50 levels between 0.68 and 5.5 kg /m3).

Figure 3 - 24 :

 324 Figure 3-24: Isocontours field of pressure for the Woodward and Colella test case at time t*=1.5 (50 levels between 40975 Pa and 1055077 Pa).

Figure 3 -

 3 Figure 3-25: isocontours of the reduced entropy of the supersonic flow past the corner of the step (S 0 is the specific entropy of the inlet flow) for the Woodward and Colella test case at time t*=1.5.

Figure 3 -

 3 Figure 3-26: isocontours of the flow velocity of a supersonic flow past the corner of the step for the Woodward and Colella test case at time t*=1.5.

Figure 3 - 27 :

 327 Figure 3-27: Schematic of a free jet in axial coordinate specifying the different zones encountered in a flammable atmosphere

Figure 3 - 28 :Figure 3 - 29 :

 328329 Figure 3-28: Experimental configuration : horizontal flow of hydrogen at 1.5 m above ground through a 12 mm orifice fed under about 40 bar

Figure 3 - 30 :Figure 3 -Figure 3 - 31 :

 3303331 Figure 3-30: Example of sensors positioning on the mast (in yellow velocity measurement)

Figure 3 - 36 : 7 :

 3367 Figure 3-36: Structure of highly under-expanded jet

Figure 3 - 37 :

 337 Figure 3-37: Domain of reference

Figure 3 -

 3 Figure 3-38: Computation domain

Figure 3 - 39 :Figure 3 -Figure 3 - 41 : 3 -Figure 3 - 42 :

 33933413342 Figure 3-39: Initial mesh (1205741 nodes)

Figure 3 - 43 :

 343 Figure 3-43: Computed values on the centerline of the jet (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C)

Figure 3 - 44 .

 344 Both parameters are very low outside the central core of 12D.

Figure 3 -

 3 Figure 3-44: mean axial velocity and hydrogen mass fraction profiles versus distance to the centerline at X=10D (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C)

Figure 3 -

 3 Figure 3-45: mean axial velocity, axial turbulent intensity and hydrogen mass/volumetric fraction on the centerline (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C)

Figure 3 -

 3 Figure 3-46 shows the profiles of the mean axial velocity and hydrogen mass fraction versus distance to the centerline at different X positions.

Figure 3 -

 3 Figure 3-46: mean axial velocity and hydrogen mass fraction versus distance to the centerline at different X positions (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C)

Figure 3 - 47 :Figure 3 - 48 :Figure 3 - 49 :

 347348349 Figure 3-47: Experimental and numerical profiles of main axial velocity and hydrogen volumetric fraction (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C)

Figure 3 - 50 :

 350 Figure 3-50: Experimental environment of GAMELAN (Open box with a vent in the upper part)

Figure 3 - 51 :

 351 Figure 3-51: Experimental setup scheme, top view on the left and side view on the right.

Figure 3 -Figure 3 - 53 :

 3353 Figure 3-52: Concentration (found experimentally) of helium registered by the sensor located on the mast 4 for the injection at 180Nl/min

  -56.

Figure 3 - 56 :

 356 Figure 3-56: Mesh used to simulate Gamelan

Figure 3

 3 Figure 3-54: Homogeneous layer regime

Figure 3 - 57 : 13 :

 35713 Figure 3-57: Numerical results of the lowest sensor of the mast 4 obtained with different meshes for the

Figure 3 - 58 :

 358 Figure 3-58: Concentration (v/v) of helium at different times in the plane of symmetry for the 180 Nl/min injection flowrate

  -58). These concentration profiles (Figure3-59) should be compared to the experimental results represented on Figure3-52.

Figure 3 - 59 :

 359 Figure 3-59: Volumetric concentration of helium registered by the sensors of mast 4 for the 180 Nl/min injection flowrate

Figure 3 - 60 :

 360 Figure 3-60: Volumetric concentration (experiments and simulations) of helium registered by the sensors of the mast 4 for the 10 Nl/min injection flowrate

Table 14 :

 14 Estimation of the non dimensional flow numbersClearly, the buoyancy forces play an important role especially at the lower flowrate. In this situation, the Reynolds numbers are all small suggesting the boundary layers may not be fully developed as implicitly assumed when using the log laws at the walls in the standard k-epsilon model.Because of this, the low Reynolds number k-epsilon was used and the calculations run again. The results are now presented on Figure3-61 and Figure3-62. The agreement is now excellen confirming the preceding assumption.

Figure 3 - 61 :Figure 3 - 62 :

 361362 Figure 3-61: Concentration of helium registered by the "sensors" of mast 4 (180 Nl/min helium flowrate, low Reynolds number k-Ɛ model)

Figure 3 - 63 :

 363 Figure 3-63: Experimental device to investigate the dynamics of vented explosions

Figure 3 -

 3 Figure 3-64: Overpressure profile hydrogen-air mixture

Figure 3 -

 3 Figure 3-65. The pressure data are those recorded by the pressure sensor P2 located near of the ent and by the sensor installed on a support ("lens"), outside, at 2m from the vent. The flame dynamics was extracted from the high speed films (velocity on the axis

Figure 3 - 67 :

 367 Figure 3-67: Flame propagation in the chamber. A treatment was applied to the picture to highlight the contour of the flame front. The latter is superimposed on the original frame (continuous red line).

  -68) adjoining the chamber to simulate the combustion outside. Because of the vertical plane of symmetry, only half of this is modelled. The second half of the geometry is replaced by a "symmetry" boundary condition. The outer dimensions of the domain are : (L = 6m) × (w = 2m) × (h = 3m).

Figure 3 - 68 :

 368 Figure 3-68: Computation geometry, the grey part is the chamber and the six transparent parts are the environment surrounding the chamber.

Figure 3 - 69 :

 369 Figure 3-69: Computation geometry, the white part is the finely meshed part and the blue part is meshed in a coarser manner.

Figure 3 - 70 :

 370 Figure 3-70: Meshes used for the chamber and the different adjoining this chamber.

Figure 3 - 71 :

 371 Figure 3-71: Computation domain obtained after the reunification of the different meshed parts.

Figure 3 - 72 :

 372 Figure 3-72: Comparison of the numerical simulation results obtained with different meshes (mesh1: 4512620 nodes, mesh2: 4892160 nodes, mesh3: 5523480 nodes). H-M model applied to the simulation of a 16.5% H 2 -air explosion in a 4 m 3 chamber equipped with a 0.5 m 2 vent.

Figure 3 - 73 :

 373 Figure 3-73: Pressure vector obtained at t=135ms from the simulation of 16.5% hydrogen-air mixture in a 4 m 3 chamber equipped with a 0.5 m 2 vent

Figure 3 -

 3 Figure 3-74: simulated overpressure (EDM model) signal from the explosion of a 16.5% hydrogen-air mixture in a 4 m 3 chamber equipped with a 0.5 m 2 vent

Figure 3 - 76 :

 376 Figure 3-76: Simulation of the external explosion (t=140ms) using the EDM model during the explosion of a 16.5% hydrogen-air mixture in a 4 m 3 chamber equipped with a 0.5 m 2 vent

Figure 3 - 78 :

 378 Figure 3-78: Comparison of numerical (EDM model) and experimental profiles of the overpressure signals (explosion of a 16.5% hydrogen-air mixture in a 4 m 3 chamber equipped with a 0.5 m 2 vent)

Figure 3 - 80 :Figure 3 - 81 :Figure 3 - 82 :

 380381382 Figure 3-80: Turbulent fluctuation velocity (on the axis) as function of time corresponding to the simulation performed with the EDM model (explosion of a 16.5% hydrogen-air mixture in a 4 m 3 chamber equipped with a 0.5 m 2 vent)

Figure 3 - 83 :

 383 Figure 3-83: Comparison of numerical (CREBCOM model) and experimental profiles of the overpressure (explosion of a 16.5% hydrogen-air mixture in a 4 m 3 chamber equipped with a 0.5 m 2 vent).

Figure 3 - 84 :

 384 Figure 3-84: Comparison of numerical (CREBCOM model) and experimental profiles of the axial flame velocity (explosion of a 16.5% hydrogen-air mixture in a 4 m 3 chamber equipped with a 0.5 m 2 vent)

Figure 3 - 85 :Figure 3 - 86 :

 385386 Figure 3-85: Comparison of numerical (CREBCOM model) and experimental profiles of the overpressure for different values of the combustion parameter C g (C g =K 0 ) (explosion of a 16.5% hydrogen-air mixture in a 4 m 3 chamber equipped with a 0.5 m 2 vent)

Figure 3 - 87 :

 387 Figure 3-87: Schematic representation (and photograph) of experimental setup used to investigate the propagation of a hydrogen-air flame in an obstructed media

Figure 3 - 88 :

 388 Figure 3-88: Average flame velocity (between the obstacles) and explosion overpressure as function of time (14.5 % H2 in air, 2 row of DN 200 cylinders, the first row at 1.5

  ndustrial) explosion simulation Page 191 sur 225 Average flame velocity (between the obstacles) and explosion overpressure as function of time (14.5 % H2 in air, 2 row of DN 200 m and the second 0.8 m further)

Figure 3 - 89 :Figure 3 - 90 :

 389390 Figure 3-89: Numerical profile of overpressure generated during the explosion 14.5% hydrogen different meshes

Figure 3 - 92 :

 392 Figure 3-92: Flow velocity along the axis ahead of the flame as function of time during the explosion

  ndustrial) explosion simulation Page 194 sur 225 Flow velocity along the axis ahead of the flame as function of time air mixture around obstacles

Figure 0 - 1 : 18 Figure 0 - 2 : 19 Figure 0 - 3 : 19 Figure 0 - 4 : 20 Figure 0 - 5 : 20 Figure 0 - 6 : 21 Figure 1 - 1 : 24 Figure 1 - 2 : 29 Figure 1 - 3 : 33 Figure 1 - 4 : 36 Figure 1 - 5 : 36 Figure 1 - 6 : 38 Figure 1 - 7 :Figure 1 - 8 : 42 Figure 1 - 9 : 64 Figure 1 - 10 : 65 Figure 1 - 11 : 68 Figure 1 - 12 : 77 Figure 1 - 13 : 81 Figure 1 - 14 :

 011802190319042005200621112412291333143615361638171842196411065111681127711381114 Figure 0-1: Intercomparison of maximum explosion overpressure inside an offshore module ................................................................................................................................. 18 Figure 0-2: Example of a congested geometry (MERGE project) representing a regular cuboidal pipe ....................................................................................................................... 19 Figure 0-3 : Comparison of calculated and measured maximum over-pressures for MERGE medium-scale experiments, (x) -COBRA predictions and (<>) -EXSIM predictions; a) all experiments and b) experiments with maximum over-pressures below 1.5 bar (5) ............ 19 Figure 0-4: Comparison of calculated and measured maximum over-pressures for MERGE large-scale experiments, (x) -COBRA predictions, (<>) -EXSIM predictions,(•) -FLACS predictions and (o) AutoReaGas predictions; a) all experiments and b)experiments with maximum over-pressures below 1 bar (5) ................................................................... 20 Figure 0-5: Test facility (ICT -ball of 2000 m 3 full of a stoichiometric mixture of H 2 -air at rest -ground ignition on the axis of symmetry). ............................................................. 20 Figure 0-6: Signals of pressure at 2 m from the boot (CAST3M-ECA, COM3D-FZK, FLACS-GEXCON, REACFLOW-JRC, AUTOREAGAS-TNO, FLUENT-UU) (a) and report maximum simulated/ measured values in 2 m, 5 m, 8 m, 18 m, 35 m and 80 m in the boot (b). ............................................................................................................................... 21 Figure 1-1: Unconfined (a) and confined (b) explosion ...................................................... 24 Figure 1-2: Description of a combustion wave using temperature ...................................... 29 Figure 1-3: Deviation of the current lines through a flame front ....................................... 33 Figure 1-4: Development of turbulence in a flow past a sphere (from H. Werlé (36)) ....... 36 Figure 1-5: Typical point velocity measurement in turbulent flow ..................................... 36 Figure 1-6: Energy spectrum of turbulence behind a grid (38) ........................................... 38 Figure 1-7: Representative diagram of combustion regime ................................................ 40 Figure 1-8: Examples of mechanical calculators ................................................................. 42 Figure 1-9: A control volume with a cylindrical obstacle inside. The gray area shows the volume occupied by the obstacle, 𝑉𝑆. ................................................................................. 64 Figure 1-10: The surface area 𝐴𝑥 of a control volume, partially occupied by an obstacle. The bright area shows the area occupied, 𝐴𝑆. ..................................................................... 65 Figure 1-11: Example of calculation domain in one dimension .......................................... 68 Figure 1-12: Upwind and centered evaluation scheme ....................................................... 77 Figure 1-13: Generation of the structured mesh: ................................................................. 81

Figure 2 - 1 :

 21 Figure 2-1: An overview of the folder MERLIN ................................................................

Figure 2 - 2 :

 22 Figure 2-2: an overview of the content of UserInput folder dedicated to the simulation of the compressible non-reacting fully three dimensional case represented by the matlab file "Fully3D". ...........................................................................................................................

Figure 2 - 3 :

 23 Figure 2-3: Fully-3D simulation ..........................................................................................

Figure 2 - 4 :

 24 Figure 2-4: example of data files obtained during a simulation for a post-processing with realized with tecplot.............................................................................................................

Figure 2 - 5 :

 25 Figure 2-5: an overview of the preprocessing solver of MERLIN ......................................

Figure 2 - 6 :

 26 Figure 2-6: computational domain example ........................................................................

Figure 2 - 7 :

 27 Figure 2-7: W component hat function given on outflow face............................................

Figure 2 - 8 :

 28 Figure 2-8: staircase approximation of user defined domain boundary ..............................

Figure 2 - 9 :

 29 Figure 2-9: cross section of the mesh used for the simulation of the unconfined explosion (expanding flame zone) .......................................................................................................

Figure 2 - 10 :

 210 Figure 2-10: An overview of the Postprocessing solver......................................................

Figure 2 - 11 :Figure 2 - 12 :Figure 2 - 13 :Figure 2 - 14 :Figure 2 - 15 :Figure 2 -Figure 2 -Figure 2 - 18 :Figure 2 -Figure 2 - 20 :Figure 2 - 21 :Figure 2 - 22 :Figure 2 - 23 :Figure 2 - 24 :Figure 2 - 25 :Figure 2 -Figure 2 - 27 :Figure 2 - 28 :Figure 3 - 1 :Figure 3 - 2 :Figure 3 - 3 :Figure 3 - 4 :Figure 3 - 5 :Figure 3 - 6 :Figure 3 - 7 :Figure 3 - 8 :Figure 3 - 9 :Figure 3 - 10 :Figure 3 - 11 :Figure 3 - 12 :Figure 3 - 13 :Figure 3 - 14 :Figure 3 - 15 :Figure 3 - 16 :Figure 3 - 17 :Figure 3 - 18 :Figure 3 - 19 :Figure 3 - 20 :Figure 3 - 21 :Figure 3 - 22 :Figure 3 - 23 :Figure 3 - 24 :Figure 3 - 25 :Figure 3 -Figure 3 - 27 :Figure 3 - 28 :Figure 3 - 29 :Figure 3 - 30 :Figure 3 - 31 :Figure 3 - 32 :Figure 3 - 33 :Figure 3 - 34 :Figure 3 - 35 :Figure 3 - 36 :Figure 3 - 37 :Figure 3 - 38 :Figure 3 - 39 :Figure 3 -Figure 3 - 41 :Figure 3 - 42 :Figure 3 - 43 :Figures

 211212213214215222182220221222223224225222722831323334353637383931031131231331431531631731831932032132232332432533273283293303313323333343353363373383393341342343 Figure 2-11: An overview of the Euler equations solvers ................................................. 100 Figure 2-12: An overview of the non reacting RANS flow solvers .................................. 100 Figure 2-13: An overview of the reacting RANS flow solvers ........................................... 101 Figure 2-14: post representation of Dirichlet velocity boundary conditions ..................... 102 Figure 2-15: Post-processing of a fully 3D Compressible Navier-Stokes simulation in a manifold of 2m height, 0.5 x 0.5 m 2 as cross Section and Re=32 (The diameters of the inlet/outlets are 0.15 m.). The lines in black are streamlines. ........................................... 102 Figure 2-16: zoom of stagnation point .............................................................................. 103 Figure 2-17: inflow recirculation streamlines ................................................................... 103 Figure 2-18: Two-dimensional blast propagation; pressure field ...................................... 104 Figure 2-19: three-dimensional blast propagation; pressure field. .................................. 104 Figure 2-20: Contour of density of a methane air diffusion flame .................................... 105 Figure 2-21: 3D structured meshes: a) Cartesian and b) skewed curvilinear .................... 110 Figure 2-22: 3D unstructured meshes: a) unstructured mesh with tetrahedral cells, unstructured mesh with prismatic cells, c) highly skewed unstructured mesh with tetrahedral cells, and d) highly skewed unstructured mesh with prismatic cells ............... 111 Figure 2-23: 3D hybrid meshes: a) skewed hybrid and b) highly skewed hybrid ............. 112

Figure 3 -Figure 3 - 45 :Figure 3 - 46 :Figure 3 - 47 :Figure 3 - 48 :Figure 3 - 49 :Figure 3 - 50 :Figure 3 - 51 :Figure 3 - 52 :Figure 3 - 53 :Figure 3 - 56 :Figure 3 - 54 :Figure 3 - 55 :Figure 3 - 57 :Figure 3 - 58 :Figures

 3345346347348349350351352353356354355357358 Figure 3-44: mean axial velocity and hydrogen mass fraction profiles versus distance to the centerline at X=10D (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C) . 153 Figure 3-45: mean axial velocity, axial turbulent intensity and hydrogen mass/volumetric fraction on the centerline (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C) ........................................................................................................................................... 155 Figure 3-46: mean axial velocity and hydrogen mass fraction versus distance to the centerline at different X positions (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C) ........................................................................................................................... 156 Figure 3-47: Experimental and numerical profiles of main axial velocity and hydrogen volumetric fraction (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C) ... 157 Figure 3-48: Experimental and numerical reduced profiles of velocity and hydrogen volumetric fraction taken respectively at x=4.5m and x=3m from the orifice (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C) .................................................... 157 Figure 3-49: Experimental and numerical reduced profiles of the intensity of turbulence taken at x=4.5m from the orifice (release of hydrogen through a 12 mm orifice fed at 40 bar/20°C) ........................................................................................................................... 158 Figure 3-50: Experimental environment of GAMELAN (Open box with a vent in the upper part) .................................................................................................................................... 159 Figure 3-51: Experimental setup scheme, top view on the left and side view on the right. ........................................................................................................................................... 159 Figure 3-52: Concentration (found experimentally) of helium registered by the sensor located on the mast 4 for the injection at 180Nl/min ........................................................ 160 Figure 3-53: Concentration (found experimentally) of helium registered by the sensor located on the mast4 for the injection at 10Nl/min ........................................................... 160 Figure 3-56: Mesh used to simulate Gamelan ................................................................... 162 Figure 3-54: Stratified layer regime .................................................................................. 162 Figure 3-55: Homogeneous layer regime .......................................................................... 162 Figure 3-57: Numerical results of the lowest sensor of the mast 4 obtained with different meshes for the 180 Nl/min injection flowrate ................................................................... 163 Figure 3-58: Concentration (v/v) of helium at different times in the plane of symmetry for the 180 Nl/min injection flowrate ...................................................................................... 165

Figure 3 - 59 :Figure 3 - 60 :Figure 3 - 61 :Figure 3 - 62 :Figure 3 - 63 :Figure 3 - 64 :Figure 3 - 65 :Figure 3 - 66 :Figure 3 - 67 :Figure 3 - 68 :Figure 3 - 69 :Figure 3 - 70 :Figure 3 - 71 :Figure 3 - 72 :Figure 3 - 73 :Figure 3 - 74 :Figure 3 - 75 :Figure 3 - 76 :Figure 3 - 77 :Figure 3 - 78 :Figure 3 -Figure 3 - 80 :Figure 3 - 81 :Figure 3 - 82 :Figure 3 - 83 :Figure 3 - 84 :Figures

 3593603613623633643653663673683693703713723733743753763773783380381382383384 Figure 3-59: Volumetric concentration of helium registered by the sensors of mast 4 for the 180 Nl/min injection flowrate............................................................................................ 166 Figure 3-60: Volumetric concentration (experiments and simulations) of helium registered by the sensors of the mast 4 for the 10 Nl/min injection flowrate..................................... 166 Figure 3-61: Concentration of helium registered by the "sensors" of mast 4 (180 Nl/min helium flowrate, low Reynolds number k-Ɛ model) .......................................................... 167 Figure 3-62: Comparison of the numerical results obtained with the different turbulence models for the volumetric concentration of helium registered by the "sensors" of mast 4 (10 Nl/min helium flowrate) .............................................................................................. 168 Figure 3-63: Experimental device to investigate the dynamics of vented explosions....... 170 Figure 3-64: Overpressure profile measured during the explosion of a quiet hydrogen-air mixture (16.5% of hydrogen) in a 4 m 3 chamber, provided with a 0.5 m 2 vent ................. 171 Figure 3-65: Measured flame dynamic profile during the explosion of a quiet hydrogen-air mixture (16.5% H 2 ) in a 4 m 3 chamber, provided with 0.5 m 2 vent ................................... 172 Figure 3-66 : Pictures of the flame taken at time 120 ms (a) and at time 140 ms (b) ....... 173 Figure 3-67: Flame propagation in the chamber. A treatment was applied to the picture to highlight the contour of the flame front. The latter is superimposed on the original frame (continuous red line). ......................................................................................................... 174 Figure 3-68: Computation geometry, the grey part is the chamber and the six transparent parts are the environment surrounding the chamber. ........................................................ 176 Figure 3-69: Computation geometry, the white part is the finely meshed part and the blue part is meshed in a coarser manner. ................................................................................... 177 Figure 3-70: Meshes used for the chamber and the different adjoining this chamber. ..... 177 Figure 3-71: Computation domain obtained after the reunification of the different meshed parts. .................................................................................................................................. 178 Figure 3-72: Comparison of the numerical simulation results obtained with different meshes (mesh1 : 4512620 nodes, mesh2 : 4892160 nodes, mesh3 : 5523480 nodes). H-M model applied to the simulation of a 16.5% H 2 -air explosion in a 4 m 3 chamber equipped with a 0.5 m 2 vent. ............................................................................................................. 178 Figure 3-73: Pressure vector obtained at t=135ms from the simulation of 16.5% hydrogenair mixture in a 4 m 3 chamber equipped with a 0.5 m 2 vent .............................................. 179

Figure 3 - 85 :Figure 3 - 86 :Figure 3 - 87 :Figure 3 - 88 :Figure 3 - 89 :Figure 3 - 90 :Figure 3 - 91 :Figure 3 - 92 :

 385386387388389390391392 Figure 3-85: Comparison of numerical (CREBCOM model) and experimental profiles of the overpressure for different values of the combustion parameter C g (C g =K 0 ) (explosion of a 16.5% hydrogen-air mixture in a 4 m 3 chamber equipped with a 0.5 m 2 vent) .............. 187 Figure 3-86: Comparison of numerical (obtained with the CREBCOM model) and experimental profiles of the flame velocity during the explosion of the hydrogen-air mixture (16.5% of hydrogen) in a chamber of 4m 3 , fitted with a vent on the opposite side of the ignition .................................................................................................................... 187 Figure 3-87: Schematic representation (and photograph) of experimental setup used to investigate the propagation of a hydrogen-air flame in an obstructed media ................... 189 Figure 3-88: Average flame velocity (between the obstacles) and explosion overpressure as function of time (14.5 % H2 in air, 2 row of DN 200 cylinders, the first row at 1.5 m and the second 0.8 m further) ................................................................................................... 191 Figure 3-89: Numerical profile of overpressure generated during the explosion of a 14.5% hydrogen-air mixture around obstacles using the EDM model with different meshes (mesh 1 : 9000000, mesh 2 : 9500000, mesh 3 : 8000000) .......................................................... 192 Figure 3-90: Comparison of experimental and numerical profiles of overpressure during the explosion of a 14.5% hydrogen-air mixture around obstacles using the EDM model 192 Figure 3-91: Flame dynamic during the during the explosion of a 14.5% hydrogen-air mixture around obstacles using the EDM model............................................................... 193 Figure 3-92: Flow velocity along the axis ahead of the flame as function of time during the explosion of a 14.5% hydrogen-air mixture around obstacles using the EDM model ...... 194

  Momentum equation:

			∂ Y ∂t	+	∂ u Y ∂x	= -	∂ u Y ∂x	+	∂𝒥 ∂x	+ ω̇	[65]
	 Total enthalpy equation h = h + u u 2 ⁄
	∂ρ h ∂t	+	∂ u h ∂x		= -	∂ u h ∂x	+	∂P ∂t	+	∂ ∂x	𝒥 + u τ + u F	[66]
			∂ u ∂t	+	∂ u u ∂x	= -	∂ u u ∂x	-	∂P ∂x	+	∂τ ∂x	+ F	[64]
	 Species equation (N species with k=1,…,N):

Table 2 )

 2 . FVS methods often have the great advantage of being robust and above all, simple to implement, but usually introduce too much numerical diffusion to properly capture the boundary layers or shear waves.

	Numerical	Expressions
	schemes
	Steger and

Table 2 :

 2 Examples of Flux Vector Splitting numerical schemes

	Flux Difference Splitting (FDS) methods

Table 3 :

 3 Slope limiters in 2 nd order numerical schemes

Table 4 :

 4 Typical boundary conditions 1.3 SOFTWARES COMMONLY USED TO SIMULATE LARGE SCALE (INDUSTRIAL) EXPLOSIONSAs already quoted, several CFD software are used in the domain of industrial safety particularly for large scale complex explosion phenomenon simulation. Examples are FLACS, FLUENT, CAST3M, TONUS, PHOENIX, AUTOREAGAS, REACFLOW, EXSIM… The main features of these softwares are listed in Table5.

	Software	Basic equations	Numerical technique	Turbulence model	Combustion models	Grid system	Numerical methods
	FLACS	Mass, Momentum Species Energy	Finite volume URANS	Standard k-epsilon model	Hjertager and Magnussen version of Eddy Break Up model (EDM)	-Structured and cartesian, -PDR for sub-grid obstacles	-First order Euler explicit for time discretization -Riemann solver for convection and central differencing for
					Beta flame		diffusion.
					(after 2000)
					Hjertager		-First order Euler
					and		explicit for time
	FLUENT	Mass, Momentum Species Energy	Finite volume URANS	Standard k-epsilon model	Magnussen version of Eddy Break Up model	-Structure, cartesian or curvilinear	discretization -Riemann solver for convection and central
					(EDM)		differencing for
							diffusion.
					-Hjertager		-First order Euler
					and		explicit for time
	CAST3M	Mass, Momentum Species Energy	Finite volume URANS	Standard k-epsilon model	Magnussen version of Eddy Break Up model	Structured cartesian	discretization -Riemann solver for convection and central
					(EDM)		differencing for
							diffusion.

Table 5 :

 5 Main models and methods used in the selected CFD software for large scale (industrial) explosion simulation

  Contains the routines (matlab files) which define the geometry (1D to 3D geometry), the type of mesh (structured or unstructured mesh), the cell size and the location of the boundaries.Postprocessor(folder)Contains routines to visualize and to analyze the data. The post processing is realized using MATLAB features or Tecplot by loading respectively the excel and the ascii data files.

	Parallelization(matlab file)	Contains the computer program dedicated to parallel computing.
	Euler equations (folder)	Contains Euler equations resolved in one, two or three dimensions. It was
		originally developed to investigate the behavior of the main numerical
		schemes used in finite volume method.
	Non-reacting flow (folder)	Contains URANS equations in two and three dimensions without combustion
		models.
	Reacting flow (folder)	Contains URANS equations in two and three dimensions where the
		combustion model (CREBCOM or H-M models) is added as source terms to
		the species balance equation.
	Preprocessor(folder)
	Incompressible RANS(folder)	Contains many computer programs dedicated to incompressible flow (These
		programs are not used in this thesis).
		Contains the requirements related to the use or the diffusion of MERLIN.
	Copyright(file.txt)
	README(file.txt)	It describes how MERLIN works.

Table 6 :

 6 Description of the content of MERLIN

Table 7 :

 7 Description of the basic options of MERLIN 2.1.2 PREPROCESSOR

Table 8 :

 8 Description of the PreProcessor solver of MERLIN

Table 9 :

 9 Description of the content of the PostProcessor solver of MERLIN 2.1.4 EULER EQUATIONS FOLDER This folder contains three solvers (Figure 2-11), Euler 1D, Euler 2D and Euler 3D, solving the Euler equations in one, two and three dimensions operating 4 other sub-solvers, representing different numerical schemes: Roe and Van-Leer numerical schemes with one order in space and time (explicit) approximations and MUSCL+Van-Leer, MUSCL+Roe numerical schemes with second order

Table 10 :

 10 Different mesh levels of refinement and mesh types used for the verification exercise (3D cases only)

		Mesh Topologies
	K	Structured	Unstructured	Hybrid
	16	8 × 8 × 8	320	1664
	8	16 × 16 × 16	2560	13312
	4	32 × 32 × 32	20480	106496
	2	64 × 64 × 64	163840	851968
	1	128 × 128 × 128	1310720	6815744

  Boundary conditions are always implemented which dictates to a large extent the final result and it is necessary to verify them as well. The above general MMS can be adapted to incorporate the verification of the various types of boundary conditions. The technique was developed by Bond et al. (101). The location of the boundary can always be represented by a surface corresponding to the equation: F(x, y, z) = C, where C is a constant usually set to zero. The Manufactured Solution retained for boundary conditions verification is:

	p = 100000 + 20000Cos(0.4πx) -50000Cos(0.45πy) + 20000Sin(0.85πz) -25000Cos(0.75πxy)
	-10000Sin(0.7πyz) + 10000Cos(0.8πzx)
	k = 780 + 160Cos(0.65πx) -120Cos(0.7πy) + 80Sin(0.8πz) + 80Cos(0.8πxy) + 60Cos(0.85πxy)
	-70Sin(0.6πz)
	ε = 150 -30Cos(0.75πx) + 22.5Cos(0.875πy) + 20Sin(0.65πz) + 40Cos(0.6πxy) -15Cos(0.75πxy)
	+ 25Sin(0.8πz)

Table 11 :

 11 Physical and geometric parameters used for the simulation of an unsteady shock propagating on a ramp:

𝑀 : Mach number of the shock, 𝜙:

Table 12 :

 12 -19) using an unstructured mesh with the Van Leer numerical scheme coupled with the VLVA limiter over 45360 nodes. In our work, the same representation and the same modelling technique as for the Katayama and Jiang configuration were implemented. Following, the shock wave comes from the left and is vertical. Bearing this in mind, a physically correct result is obtained showing a normal shock wave structure. A quantitative comparison between the simulations and the experiment is proposed in Table12. Comparison between some experimental values and simulated data (van Leer 2nd order in space, VLVA limiter, unstructured mesh over 45360 nodes)

	Region	ρ/ρ 0		Mach stem length/incident shock-ramp
				apex distance
		Measured	Calculated	Measured	Calculated
	1	2.71	2.7	0.275	0.275
	2	3.68	3.47	Reflected shock front-ramp apex
				distance / Mach stem -ramp apex
				distance
	3	3.33	3.21	0.38	0.38
	j	4.06	4.08

Table 15 .

 15 .

	Thermodynamic	Unburnt gases		Burnt gases
	parameters		Hydrogen(H 2 )	Air	Hydrogen(H 2 )	Air
	Molecular weight(g/mol)	0.32	24.33	0.32	24.33
	Density(kg/m	3 )	0.0838	1.205	0.0112	0.16
	Specific heat at
	Constant pressure (kj/kgK)	14.1889	1.005	16.59	1.35452

Table 15

 15

	: Thermodynamic parameters calculated at ambient temperature
	(293.15k) and pressure (1atm) for unburnt gases. The burnt gases and
	combustion products parameters are calculated at 2200k (the ignition
	temperature).
	The boundary conditions are:

Table 1 :

 1 Intensities of turbulence for different cases of turbulence (81), (82) ................... 58

Table 2 :

 2 Examples of Flux Vector Splitting numerical schemes ........................................ 72

Table 3 :

 3 Slope limiters in 2 nd order numerical schemes ...................................................... 78

Table 4 :

 4 Typical boundary conditions ................................................................................. 86

Table 5 :

 5 Main models and methods used in the selected CFD software for large scale (industrial) explosion simulation ......................................................................................... 87

Table 6 :

 6 Description of the content of MERLIN ................................................................. 92

Table 7 :

 7 Description of the basic options of MERLIN ....................................................... 95

Table 8 :

 8 Description of the PreProcessor solver of MERLIN ............................................. 96

Table 9 :

 9 Description of the content of the PostProcessor solver of MERLIN .................... 99

Table 10 :

 10 Different mesh levels of refinement and mesh types used for the verification exercise (3D cases only) .................................................................................................... 110

Table 11 :

 11 Physical and geometric parameters used for the simulation of an unsteady shock propagating on a ramp: 𝑀𝑠: Mach number of the shock, 𝜙: angle between the direction of the ramp and the longitudinal direction, 𝑋𝑖: relative position of the corner of the ramp relative to the inlet (distance/ length of the computational domain),𝑋𝑓: relative position of the incident shock at the end of the simulation. ................................................................ 127

Table 12 :

 12 Comparison between some experimental values and simulated data (van Leer 2nd order in space, VLVA limiter, unstructured mesh over 45360 nodes) ....................... 132

Table 13 :

 13 Meshes used to simulate GAMELAN ............................................................... 163

Table 14 :

 14 Estimation of the non dimensional flow numbers ............................................. 167

Table 15 :

 15 Thermodynamic parameters calculated at ambient temperature (293.15k) and pressure (1atm) for unburnt gases. The burnt gases and combustion products parameters are calculated at 2200k (the ignition temperature). ........................................................... 175

	6. LISTE OF FIGURES

http://www.fabig.com/( FABIG newsletter -issue n°22, may 1998-article R320)

because the pressure is equalised at the speed of sound, so that as long as the flame speed is much less than the speed of sound, the pressure can be assumed to be uniform.

Note that the integral scale should be comparable to the flame thickness suggesting a very small equipment.

http://www.gexcon.com/index.php?/flacs-software/article/FLACS-Overview.

http://www.ansys.com/fr_fr/Produits/Flagship+Technology/ANSYS+Fluent.

http://www-cast3m.cea.fr/.

Because the Kolmogorov scale is the inner limit of the "inertial" range of the turbulence cascade where the dissipation by viscous forces are neglected. So the inertial forces are still significant at the Kolmogorov scale.

MERLIN

If for instance the temporal discretization error term is too small when compared to the spatial discretization error term, then mistakes in the temporal discretization will not be seen on very fine meshes.

We verified that the simulations are not very sensitive to the value of F in the interval 0.2 to 0.7.

a solver is defined here as a set of system of equations, models, numerical scheme, mesh producing a result once the initial and boundary conditions of the configuration are given.

ACKNOWLEDGEMENTS

First of all, I would like to express my gratitude to my PhD thesis director Christophe PROUST for allowing me to do my thesis at UTC (Compiègne

LISTE OF TABLES