Cléante Langrée

Alberto Victor

Leandro

Emmanuel Hector

Fernando

Conclusion and perspectives

During regulatory hazard studies, governmental or private organizations are called upon for their expertise in predicting the consequences of accidental scenarios. During these hazard studies, effect distances must be determined. Several methods and tools are used to determine these effect distances. First based on phenomenological methods, more complex numerical modeling has become more affordable and applicable to real industrial situations. Computational fluid dynamics (CFD) is one of these methods.

One of the most devastating accident scenarios is the unconfined vapour cloud explosion (UVCE). Following a loss of containment of a flammable fluid on an industrial site (gasoline, natural gas, hydrogen, etc.), a gas plume can form and mix with the surrounding air to form a flammable mixture. This flammable mixture can potentially spread over a large area. Upon contact with an energy source within the area where the cloud has spread, a combustion reaction can be initiated and propagate throughout the cloud until all fuel has been consumed. The propagation of this turbulent flame front will generate overpressure in the near and far field that can cause heavy human and material damage.

This thesis investigates the simulation of UVCE employing CDF. The main objectives are to improve the skills of organizations like INERIS in this field, and to increase the confidence in predictions produced by CFD methods. This manuscript is organized thematically and according to the path followed during these three years of study.

First, a review of turbulent combustion models is proposed. Several modeling approaches and techniques are discussed. The aim was to identify the advantages and limitations of each approach for a better evaluation in coherence with the final objective of an UVCE simulation. The UVCE is a very particular case of turbulent combustion, since it is on a very large scale and in a completely open environment. After a review of turbulent combustion models, a so-called geometrical approach is selected for the CFD simulations. The global system of equations is closed by modeling the subgrid flame wrinkling via a turbulent flame speed correlation. It should be noted that this geometrical approach allows to distinguish the chemical and the topological contribution to the propagation of the turbulent flame front. This differentiation allows, in particular, to integrate, relatively easily, a partial premixed effects.

A description of OpenFOAM, the numerical tool used for the application of this theoretical model is then presented. OpenFOAM is an open source C++ platform for modeling and simulation including CFD. It is used for a wide variety of problems and includes many solvers for different configurations: two-phase flow simulation, reactive flows, compressible/incompressible flows, etc... This code is based on the finite volume method described in the chapter 3. This chapter concludes with a description of the solver used in the simulations of this thesis.

Cléante Langrée

UVCE Simulation

Résumé -Français

Lors d'études de danger réglementaires, des organismes gouvernementaux ou privés sont sollicités pour leur expertise dans la prédiction de conséquence de scénarios accidentels. Lors de ces études de danger, des distances d'effets doivent être déterminées. Plusieurs méthodes et outils sont utilisés pour déterminer ces distances d'effets. D'abord basée sur des méthodes dites phénoménologiques, une modélisation numérique plus complexe est devenue plus abordable et applicable aux situations industrielles réelles. La mécanique des fluides numérique fait partie de ces méthodes.

Un des scénarios accidentels les plus dévastateurs est l'explosion de nuage de gaz non confiné (UVCE). À la suite d'une perte de confinement d'un fluide inflammable sur site industriel (essences automobiles, gaz naturel, hydrogène, etc.), un nuage de gaz peut se former et en se mélangeant à l'air ambiant former un mélange inflammable. Ce mélange inflammable peut potentiellement se répandre sur une large zone. Au contact d'une source d'énergie au sein de la zone où le nuage s'est dispersé, une réaction de combustion peut s'enclencher et se propager à travers tout le nuage jusqu'à ce que tout le carburant ait été consommé. La propagation de ce front de flamme turbulent va générer une surpression en champ proche et lointain pouvant causer de lourds dégâts humains et matériels.

Cette thèse se propose d'étudier la simulation d'UVCE par des méthodes de mécanique des fluides numérique. Les objectifs principaux sont l'amélioration des compétences d'organisme tel que l'INERIS en la matière, ainsi que l'augmentation de la confiance dans des prédictions produites par de telles méthodes. Ce manuscrit est organisé de manière thématique et selon le cheminement suivi lors de ces trois années d'étude.

Premièrement, une revue des modèles de combustion turbulente est proposée. Plusieurs approches et techniques de modélisation sont évoquées. L'intérêt de cette revue était d'identifier les caractéristiques avantages et limites de chaque approche afin de les évaluer en cohérence avec l'objectif final d'une simulation d'un UVCE. L'UVCE est un cas très particulier de combustion turbulente, puisqu'à très large échelle et dans un environnement complètement ouvert. Après cette revue des modèles de combustion turbulente, une approche dite géométrique est sélectionnée pour la simulation d'UVCE. Le système global d'équations est fermé par une approche de modélisation basé sur le plissement de sous maille, lié à une corrélation de vitesse de flamme turbulente. Il est à noter que cette approche géométrique permet de distinguer la contribution chimique et la contribution topologique à la propagation du front de flamme turbulente. Cette différenciation permet notamment d'intégrer, relativement facilement, un prémélange partiel.

Cléante Langrée

UVCE Simulation February 8, 2022 iv Une description de l'outil numérique utilisé pour l'application de ce modèle théorique est ensuite présenté. OpenFOAM est une librairie C++ open source pour la simulation numérique de mécanique des fluides. Il est utilisé pour une grande variété de problématiques et inclut de nombreux solveurs pour de nombreuses configurations différentes : simulation d'écoulement biphasique, réactifs, compressible/incompressible, etc... Ce code est basé sur la méthode des volumes finis décrite au chapitre 3. Ce chapitre se conclut par une description du solveur utilisé dans le cadre des simulations de cette thèse.

Un test de ce modèle a été réalisé en deux étapes. Premièrement une configuration 1D laminaire est étudiée. Au cours de cette étude, les caractéristiques de flamme laminaire telles que la vitesse de flamme laminaire et la température de flamme laminaire sont retrouvées. Ensuite, une configuration 3D sphérique est étudiée afin d'évaluer l'impact d'obstacles sur la propagation d'un front de flamme prémélangée turbulent. Un fort impact de la turbulence initiale et de la présence d'obstacle est retrouvé. Différentes corrélations de vitesse de flamme turbulente pour la fermeture du terme de plissement de sous-maille ont été testées.

Enfin, un cas d'application est proposé comme chapitre final du manuscrit. Ce cas d'application est la simulation numérique d'une configuration expérimentale réalisée en 2014 sur le site expérimental de l'INERIS. L'objectif est de comparer les prédictions du modèle sélectionné au chapitre 2, numériquement appliqué selon les outils présentés au chapitre 3 et testé lors du chapitre 4 aux résultats obtenus expérimentalement. Les données considérées sont principalement l'historique de position du front de flamme, ainsi que les signaux de pression. La vitesse de flamme est globalement reproduite numériquement. Néanmoins, un phénomène de brusque changement de vitesse est constaté expérimentalement, sans être reproduit numériquement. Une discussion sur cette divergence est tenue au cours du chapitre 5. Les signaux de pressions obtenus numériquement sont globalement cohérents avec les niveaux de surpression observés expérimentalement et cohérents avec l'historique de flamme numérique associé. En conclusion, ces résultats paraissent prometteurs et encouragent la poursuite d'étude sur l'apport potentiel de la CFD aux problématiques de sécurité industrielle.

vi

A test of this model was performed in two steps. First, a 1D laminar configuration is studied. During this study, laminar flame characteristics such as laminar flame velocity and laminar flame temperature are investigated. Then, a 3D spherical configuration is studied to evaluate the impact of obstacles on the propagation of a turbulent premixed flame front. A strong impact of the initial turbulence and the presence of obstacles is found. Different turbulent flame velocity correlations for the closure of the subgrid wrinkling factor have been analysed.

Finally, new CFD case is proposed in the final chapter of the manuscript. This application case represents a numerical simulation of an experimental configuration performed in 2014 on the INERIS experimental site. The objective is to compare the predictions of the model selected in chapter 2, numerically applied according to the tools presented in chapter 3 and tested during chapter 4 to the experimental results.The data considered are mainly the historical position of the flame front, as well as the pressure signals. The flame speed is globally recovered numerically. Nevertheless, a phenomenon of abrupt change of velocity is observed experimentally, without being reproduced numerically. A discussion on this divergence is held during the chapter 5. The pressure signals obtained numerically are globally consistent with the experimentally observed overpressure levels and consistent with the associated numerical flame history. In conclusion, the CFD results are very promising and encourage further studies on the potential contribution of CFD to investigate problems related to safety in industrial environments.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

Chap. 1 | Introduction 1.1 Context of this PhD

Industrial sites have the potential to be the set of both spectacular and destructive events. Industrial catastrophes like the ones of Bhopal [START_REF] Broughton | The Bhopal disaster and its aftermath: a review[END_REF], Seveso [START_REF] Homberger | The seveso accident: its nature, extent and consequences[END_REF], or AZF [START_REF]First lessons of the Toulouse ammonium nitrate disaster, 21st September 2001, AZF plant, France[END_REF], are well known because of their consequences and the media impact they had on public opinion. The events previously cited are all different in the physical description of the incident, but they all have in common a human cost, both immediately and in the long run and an enormous material cost.

A company owning an industrial plant where one or several processes imply hazardous (i.e. toxic and/or flammable) materials, is responsible for the safety of its operations. A simple management method is the risk analysis. It can take several forms (HAZOPs, Preliminary Risk Analyses, ...) but always aim at identifying: initiating events (leak, presence of an ignition source, ...), central critical events (explosion, loss of confinement, ...) and barriers, for prevention and mitigation effects. A given central critical event can lead to several dangerous phenomena (explosion, fire, toxic material release, ...).

In order to get a hierarchy of the dangerous phenomena at the end of the risk analysis, a scoring matrix is established by the team in charge of this task. This team in general gathers at least people responsible of the process operations and of the safety of the plant. The scoring matrix is based on a scale with a few levels for the intensity (I i ) and the probability of occurrence (No). Each box of the Ni x No scoring matrix corresponds to a couple intensity/probability and all the boxes should be sorted by the team as: 'acceptable', 'acceptable under condition' or 'not acceptable'. An example of the matrix is shown in Fig. 1.1. In this one, the color red stands for 'not acceptable', the orange one for 'acceptable under condition' and green means 'acceptable'. A main output of the risk analysis is a scoring matrix filled with the identified dangerous phenomena. Depending on the criteria set by the owning company or by the regulations (for example: 'no unacceptable dangerous phenomenon in the matrix', maybe completed with 'at most two acceptable dangerous phenomena under condition in the matrix'), the scoring matrix itself is declared 'acceptable' or 'unacceptable'. In this latter case, some simple principles can be used to displace the problematic dangerous phenomena towards the green boxes:

Introduction

• reducing the stored inventory of a product involved in a dangerous phenomenon with high effects distances, • displacing the capacity susceptible to undergo a loss of confinement as far as possible from the site limits, • increasing the level of reliability of technical barriers designed to prevent a given dangerous phenomenon to reduce its probability of occurrence, • ... Nevertheless, such solutions cannot be systematically retained due to technical and economical constraints which are specific to the plant.

As seen previously, the risk analysis, which is the core of the risk management strategy of a given plant requires the evaluation of characteristic effects distances. In a context of safety, an overestimation of the distances is always wished. Nevertheless, this overestimation should be as slight as possible. Important effects distances can easily make a scoring matrix unacceptable, preventing operations on the industrial plant or lead to install costly and overdimensioned mitigation devices on site or increase the protection of surrounding habitations. These aspects make a precise quantification of the effects distances necessary.

The current thesis will be focused on a typical scenario for an industrial plant storing or using gaesous or liquid flammable products: the UVCE. This acronym commonly used in risk analyses stands for 'Unconfined Vapor Cloud Explosion'.

UVCE

3

• Mixing with oxygen from the air to form a flammable volume • Inflammation of the flammable cloud • Propagation of a flame front. Associated with the expansion of the burnt gases, the flame acts like a piston on the surrounding fresh gases and can be at the origin of the formation of an air pressure wave.

An UVCE is one of the most devastating accidents that can occur in an industrial facility [START_REF] Proust | Formalisation du savoir et des outils dans le domaine des risques majeurs (eat-dra-76) les explosions non confinées de gaz et de vapeursω uvce[END_REF]. As described before, such a phenomenon first requires the formation of a gaseous flammable cloud. Mixing between the flammable product and air can be direct if the product is gaseous in the pressure and temperature of the storage or follow a phase of thermodynamic flash quite close to the release point or a vaporization if a liquid pool is formed on the ground. The flammable product is characterized by flammability limits: the Lower Flammability Limit (LFL) gives the minimum quantity of flammable product that should be mixed with air to get a flame if an ignition source is present. The Upper Flammability Limit (UFL) is the maximum amount of the product that can be burned in air. LFL and UFL depend on the initial temperature and pressure of the mixture. Some LFL and UFL values are given for products used in industrial processes in Table 1 1.2 -Flammability limits of some flammable products encountered in the industrial processes. Source [START_REF] Petit | Ed911 -les melanges explosifs -1. gaz et vapeurs[END_REF].

For a given cloud, the flammable cloud is the portion in which the product mass fraction value is between LFL and UFL. If an ignition source, powerful enough, is located in this part a premixed flame can be created. This structure is a thin interface that separates burnt and fresh gases. The interface propagates towards the latter. This propagation in the flammable cloud is susceptible to increase the pressure of the burnt gases and to create a zone of compressed gases in front of the flame. The propagation flame speed is piloted by several local parameters: the length on which the flame can propagate, the composition of the flammable mixture, the distribution of the flammable product in the flammable cloud, the presence of obstacles and of confinement ... The explosion can also produce a large amount of burned gases able to impact structures and people around through thermal convective and radiative effects. Gugan [START_REF] Gugan | Unconfined vapour cloud explosions[END_REF] studied several accidental ignitions of flammable clouds. He separated fireballs from explosions, the former giving only thermal effects, the latter leading to pressure and thermal effects. Generally the pressure effects are much more closely studied for explosions as their distance effects are more important than the thermal effects ones.

It should be emphasized that the pressure effects generated by the flame propagation directly depend on the history of the flame front. Two propagation regimes are possible: the deflagration and the detonation. When the ignition energy of the flammable mixture is low enough (lower than the energy of direct detonation initiation for the considered mixture), the initiated flame front propagates around a kernel of burned gases at a speed lower than the speed of sound in the fresh gases. In practice, this speed is about a few dozens or Introduction hundreds of meters/s. An example of pressure wave generated by a deflagration is shown below. It is characterized by a positive and negative overpressures, both being of the same order of magnitude and with a specific duration. In the framework of industrial risks, the main parameter of interest is the positive overpressure peak. The positive pressure impulse can also be of interest for computing mechanical structures responses to the pressure wave impact. If the deposited energy is high enough, a direct initiation of a detonation is possible. If the speed of a deflagration flame front increases sufficiently, a phenomenon known as DDT (Deflagration-to-Detonation Transition) can occur. With a detonation, the front separating burnt and fresh gases is coupled with the pressure wave, now a shock wave. The propagation of the interface is explained by successive auto-ignition events of the fresh gases compressed by the shock wave. The characteristic speed of propagation is about a few thousands of meters/s. It should be kept in mind that the magnitude of the overpressure peaks reached in a given flammable cloud is higher that the ones related to deflagrations. Furthermore, the shape of the generated pressure wave is different as it integrates a shock.

Other consequences of a release of a flammable product can be pool fires [START_REF] Jolly | Formalisation du savoir et des outils dans le domaine des risques majeurs (eat-dra-76) modélisations de feux industriels ω 2[END_REF], or jet fires [START_REF] Jolly | Formalisation du savoir et des outils dans le domaine des risques majeurs (eat-dra-76) feu torcheω 8[END_REF]. The first phenomenon follows a jet explosion if the flammable product release goes on after the explosion or if the flammable cloud is ignited as it is created. The second phenomenon is related to the ignition of vapors above a flammable liquid pool. These phenomena are feared because of their potentially important thermal, convective and radiative, effects. The pressure effects linked to these phenomena are negligible.

In the case of an UVCE, safety distances are generally characterized by computing the distances at which specific overpressure magnitudes are reached.

In France, for industrial plants that have to provide a Danger Study to the public administration, the distances at which the overpressure effects thresholds given in Table 1.3 are reached for the explosion scenarios have to be supplied.

Cléante Langrée

University of Rouen -CORIA February 8, 2022 In industrial history, several typical cases of UVCE can be cited to assess for their destructive potential.

One of the investigation methods often used to characterize the pressure source is to relate the observed damages on material targets (windows, cars, walls, ...) to those that can be generated by a given mass of TNT. This method is not perfect, it assumes the pressure decay related to a condensed phase detonation may differ from the one generated by a gaseous deflagration, but it gives orders of magnitude for the strength of the pressure source.

Below are detailed some major industrial accidents involving an UVCE.

Port Hudson, 1970 [START_REF] Burgess | Detonation of a flammable cloud following a propane pipeline break: the december 9[END_REF] The 9th of December 1970, at Port Hudson (USA), a 200mm diameter propane pipeline under 66 bars broke. The leakage with a flow rate around 0.08 m 3 /s during 24 minutes led to the formation of an approximately 500m diameter, flat propane-air cloud. It covered a hilly land with woods and crops. It is suspected ignition occurred in a concrete room. The observed damages allowed to estimate an energy release equivalent to 50 tons of TNT with an overpressure level of 1 bar in the flammable cloud. A detonation is considered as possible. Flixborough, 1974 [START_REF] Venart | Flixborough: the explosion and its aftermath[END_REF] This accident occurred the 1st of June 1974 at Flixborough (UK). It resulted from a break of a 700mm cyclohexane pipeline under 10 bars. 100 tons of cyclohexane were released. It is assumed the release took the form of a highly turbulent jet and the flame developed in a place containing numerous process pipes. The evaluation of the far-field damages allowed to estimate an energy release equivalent to 15-20 tons of TNT. An overpressure level of 1 bar was also estimated at the limit of the cloud.

Cléante Langrée

University of Rouen -CORIA February 8, 2022 Buncefield, 2005 [START_REF] Johnson | The potential for vapour cloud explosions -lessons from the buncefield accident[END_REF] The 11th of December 2005, at Buncefield (UK), an overflow of an gasoline tank by its roof and to the release of 180 tons of product. The released product fell into the retention bund, a large part of product going out of this latter. The dispersion conditions formed a flammable cloud on approximately 120 000 m 2 , with a height of a few meters. According to the post-event investigation, overpressures about 2 bar could be reach within the flammable cloud. At distances lower than 500 m but out of the cloud, damages corresponded to an explosion of 7 tons of TNT. Further, the equivalence was assumed to be about dozens of tons.

Discussions are still open about the nature of the flame propagation. Some works explain the effects of the explosion with a peculiar kind of deflagration [START_REF] Atkinson | Review of vapor cloud explosion incidents[END_REF] whereas others prefer the assumption of a DDT occurrence in a part of the cloud covering dense vegetation [30] [31]. All these examples of industrial explosions show that:

• flammable clouds with a characteristic dimension about a few hundreds of meters can be generated in case of a massive leak, • overpressure about 1 bar or even higher can be reached in the flammable cloud, • damages can be observed hundreds of meters away from the explosion center.

Introduction

Furthermore, in some cases, the interaction between obstacles (pipes, trees) is suspected to have boosted the flame and increased pressure effects in the flammable cloud.

The ignition mode is often described as it can play a role in the flame acceleration. Indeed, if ignition occurs in a zone in which turbulence is pre-existing, the flame front can be quickly accelerated. If a set of obstacles is located close to the ignition site, a quick flame would interact with it, promoting a higher flame speed after flame/obstacles interaction than the one that would be obtained if ignition occurred in a cloud at rest. Pre-existing turbulence could be due to: the jet of flammable products at the breech or a secondary explosion. If ignition happens in a confined area, a jet of fresh products can be generated if there is a permanent opening or a surface which resists less to explosion pressure than the others. This jet being followed by a flame, a turbulent ignition of the unconfined flammable cloud is possible.

According to other feedbacks from accidental explosions (St Herblain in 1991 [START_REF] Lechaudel | Assessment of an accidental vapour cloud explosiona case study : Saint herblain[END_REF] and La Mede in 1992 [START_REF] Michaelis | Methods applied to investigate the major uvce that occured in the total refinery's fluidcatalytic cracking unit at la mede[END_REF]) with important pressure effects, important directional effects can be observed:

• structures damages in certain directions from the assumed ignition center are explained by overpressures about a few hundreds of mbar, • in other directions, at the same radius, the observed damages prove the overpressure was not higher than dozens of mbar.

It should be kept in mind that the example mentioned corresponds to historic accidents. UVCE characterized with a lower scale (for example a flammable cloud with a width of dozens of meters) could also lead to severe structural damages and threaten people or workers life.

Usual calculation methods for the explosion

Many entities such as industrial companies of the Oil and Gas branch (Total Energies, Shell ...), consultants (DNV-GL, GexCon, ...) and public institutions like INERIS led studies or took part in research programs to gain knowledge about UVCE and propose specific modeling tools.

The main source of knowledge concerning industrial explosions is simply of the feedback of such events [START_REF] Lechaudel | Assessment of an accidental vapour cloud explosiona case study : Saint herblain[END_REF][START_REF] Michaelis | Methods applied to investigate the major uvce that occured in the total refinery's fluidcatalytic cracking unit at la mede[END_REF]. Experimental campaigns are often led in the wake of investigations to test one or several hypotheses proposed by the explosion experts to explain how an explosion took place [START_REF] Davis | Investigation techniques used to determine the massive vapor cloud explosion at the buncefield fuel depot[END_REF][START_REF] Daubech | Propagation of a confined explosion to an external cloud[END_REF]. It should be noted that the scale of the experiment is a parameter of importance as physical phenomena relative to a scale can explain flame acceleration. Most of the time, if the experimental campaign is performed at small-scale [START_REF] Daubech | Un-ignited and ignited high pressure hydrogen releases: concentration -turbulence mapping and overpressure effects[END_REF], it will be highly instrumented whereas the sensors used for large-scale experiments will be less numerous [START_REF] Puttock | Explosion assessment in confined vented geometries. solvex large-scale explosion tests and scope model development[END_REF].

While gaining knowledge, predictive models were proposed. They can be phenomenological or rely on CFD. Without detailing each one of the available model, here's a brief review of a few of them.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

1.2 UVCE 9 Equivalent TNT method
This method is based on the hypothesis that the overpressure field of any explosion could be reproduced by the detonation of a certain amount of TNT. The equivalent TNT mass M TNT is found by :

M TNT = a E gas E TNT (1.1)
with E gas being the potential energy that the flammable mixture can release, E TNT is the energy released by 1kg of TNT and a is the explosion yield.

The maximum overpressure is then determined via a pre-determined abacus. This method cannot be used to quantify the peak overpressure in the space crossed by the flame propagation. The physical meaning of the explosion yield is different, depending on the representation of E gas that can be estimated as :

• the total mass of the flammable product accidentally released, • or the mass which is contained in the flammable cloud.

In both cases, the value of the yield parameter is most often estimated with a statistical analysis from past industrial accidents.

Introduction

Congestion Assessment Method [START_REF]Developments in the Congestion Assessment Method for the prediction of vapour-cloud explosions[END_REF] This method recommends to identify the portions of space that are congested due to the presence of obstacles. The important flame speeds and overpressures are assumed to be reached in these zones due to a loop 'flame propagationflow of fresh gases ahead of the flame front -turbulence generation in this zoneflame acceleration'.

The maximum overpressure has to be quantified thanks to a decision tree [START_REF] Puttock | Fuel gas explosions guidelines -the congestion assessment method[END_REF]. Such a tree could lead to ask the opinion of an explosion specialist numerous cases. Once the pressure is known, it has to be adjusted with a factor depending on the flammable product.

The long-range pressure decay is computed using curves which were obtained from measurements carried out in the MERGE project [START_REF] Mercx | Overall final report of the merge project[END_REF]. The measured curves were generated with quasi-hemispherical deflagrations propagating at varying speed.

Baker-Strelhow method [START_REF] Baker | Recent developments in the Baker-Strehlow VCE analysis methodology[END_REF] The Baker-Strehlow method shares characteristics withe the Multi-Energy method, detailed in the next part. Thus, consequences of an accidental explosion are assumed to depend on the energy potentially at stake in the flammable cloud but also on the presence of solid obstacles in the flame path.

First, flame speeds have to be estimated with recommendations. They can take the form of a table that provides orders of magnitude depending on three parameters: the type of confinement (1-D, 2-D or 3-D) and its impact on flame propagation, the density of obstacles (weak, average or strong) and the reactivity of the mixture (weak, average or strong).

Once the violence of the explosion is known, a pressure decay curve has to be selected in the set proposed by Strehlow [START_REF] Strehlow | Blast wave generated by spherical flames[END_REF]. These curves were obtained with a simulation tool. A set of flames propagating spherically with a fixed speed was computed for speeds ranging from a few m/s to about 2000 m/s. This latter situation mimicking a stable detonation.

According to Strehlow, its approach supplies overpressures levels which are stronger than those that would be obtained with a varying flame speed, if the chosen flame speed corresponds to the observed maximum flame speed. Van Wingerden [START_REF] Van Wingerden | Investigation into blast produced by non-steady flames[END_REF] limits this statement to rapid deflagrations ie with a propagation speed higher than 340 m/s. For slower flames, the results seem to depend more on flame acceleration and on the place where it occurs. [START_REF] Van Den | The multi-energy method -a framework for vapour cloud explosion blast prediction[END_REF] It should be kept in mind that an explosion generates high overpressures if the propagation speed is high enough. As seen previously, a flame propagating in a gaseous flammable mixture accelerates if this latter contains repeated obstacles and partially confined spaces. Without obstacles or confined spaces in an ignited flammable cloud, the observed overpressures are quite weak (a few mbar).

Multi-Energy Method

Cléante Langrée

University of Rouen -CORIA February 8, 2022

1.2 UVCE 11
Thus, assuming a flammable cloud covers free or weakly congested zones and highly congested zones, only the latter are responsible of high ovepressures generation when crossed by the flame. The explosion can be seen as a set of elementary explosions occurring in the zones composing the flammable cloud.

The modeler has to identify each elementary explosion and quantify their effects separately. In this end, a violence index has to be chosen for each explosion. It ranges from 1 to 10. 10 corresponds to a detonation and the other indices to deflagrations, the corresponding speed increasing with the index value.

Once the explosion index and energy are known, the pressure decay can be deduced from an abacus. It was built solving numerically Euler equations, considering constant-speed hemispherical charges explosions. Figure 1.8 -Abacus related to the Multi-Energy method giving the overpressures generated by deflagrations at constant flame speed for hemispherical charges located on the ground [START_REF] Van Den | The multi-energy method -a framework for vapour cloud explosion blast prediction[END_REF].

CFD

In Fluid Mechanics, partial differential equations are used to describe the behavior of any fluids in motion, called the Navier-Stokes equations. Although they are the result of years of research and efforts by many famous physicists like Leonhard Euler (1707-1783) or Joseph-Louis Lagrange (1736 -1813), these equations are named after Claude-Louis Navier (1785 -1836) [START_REF] Navier | Mémoire sur les lois du mouvement des fluides[END_REF] and George Gabriel Stokes (1819Stokes ( -1903) ) [START_REF] Stokes | printed in 1847. On the theories of the internal friction of fluids in motion and of the equilibrium and motion of elastic solids[END_REF]. Mathematically, those equations remain unsolved analytically. But with the increase in computational power, approximated solutions can be found with iterative methods on points of a grid that discretize the computational domain. Thus complex behavior of fluids can be simulated. Those equations previously cited will be described in the next chapter of this thesis.

Introduction

To apply CFD to combustion modeling, enriched physics to the equations must be added to take into account the chemical reaction of combustion. In addition to the equation of continuity (conservation of mass) and momentum (conservation of linear momentum), equations for the conservation of both energy and species mass fraction must be solved [START_REF] Veynante | Turbulent combustion modeling[END_REF].

In the industrial area, CFD simulation of premixed flames has been used largely for engine's combustion chamber studies [START_REF] Richard | Towards large eddy simulation of combustion in spark ignition engines[END_REF][START_REF] Gicquel | Large eddy simulations of gaseous flames in gas turbine combustion chambers[END_REF] and explosions [START_REF] Hansen | Using computational fluid dynamics (cfd) for blast wave predictions[END_REF]. One of the first CFD modeling explosion was performed in the early 80s [START_REF] Hjertager | Computed simulation of turbulent reactive gas dynamics[END_REF]. The geometry was a tube of 50 m 3 with inner rings and filled with a methane-air or propane-air stoichiometric mixture.

In order to increase the scale of the computational domain and deal with complex obstacles sets, the Porosity Distributed Resistance (PDR) approach was proposed [START_REF] Hjertager | Computer modelling of gas explosion propagation in offshore modules[END_REF]. It is a supplemental modeling layer to classic CFD modeling. The obstacles that are not resolved by the mesh are accounted for in the transport equations trough the notion of area and volume porosity. This kind of approach was used by many CFD codes dedicated to explosion modeling in a safety framework.

Additionally to this kind of modeling, CFD modeling usage with no PDR appeared more and more in the past few years for medium scale explosion [START_REF] Bauwens | Vented explosion overpressures from combustion of hydrogen and hydrocarbon mixtures[END_REF][START_REF] Quilattre | Simulation aux grandes echelles d'explosions en domaine semi-confine[END_REF]. This is partly explained by the increase in computational power.

Phenomenological tools can give a result very quickly, simplifying the explosion characterization to the evaluation of an energy and a pressure violence or an explosion yield. Nevertheless, these estimations may not be always straightforward to perform for the modeler.

CFD can theoretically provide more precision in the computation of large-scale UVCE as it intrinsically integrates a description of the geometry of the flammable cloud and of the elements it covers (obstacles, confined zones, ...). The effects of the preferential propagation flame direction and of the flame speed history on the pressure field could be accounted for. Although the computational power keeps on increasing with time it remains limited. Then, it is not currently possible to resolved all the physics involved in large-scale UVCE and modeling choices have to be made, notably for modeling turbulence, combustion/turbulence interactions. Questions still arise concerning the physics that should be necessarily modeled in explosions for predicting accurately the pressure effects.

Objective and approach of the study

The risk management engineer would like to get a CFD tool that could predict industrial explosions, whatever their scale and their intensity, based on current knowledge. The theoretical limits (scale of the explosion, maximum flame speed accounted for, for example) of such a tool should be identified. The tool should be able to evolve with the expected increase in computational power in the years to come. At last, the engineer wishes to benefit from the best ratio "accuracy of results"/"time needed to get the results".

Cléante Langrée

University of Rouen -CORIA February 8, 2022

1.3 Objective and approach of the study 13

The aim of this thesis is to evaluate the performance and applicability of CFD methods in a context of risk assessment missions, based on two major conditions, fidelity of the results and CPU cost of the simulation. Modeling choices and numerical approaches will be discussed keeping in mind those two objectives.

The open source CFD code OpenFoam will be used, with the final goal of developing methods and solver for risks management engineers. This study will be organized in three parts. First a review of the turbulent combustion modeling. Turbulent combustion modeling is a very old and vast research topic and in order to push forward the UVCE modeling a general but precise review of the models, their limits and capability will be conducted. The theoretical models will be split in two parts, the "general" turbulence modeling, the turbulent combustion models which are specific to combustion.

Two CFD studies are carried out for the phenomenon of flame/obstacle interaction. Such a phenomenon is often present in industrial explosions.

In a first study, a configuration whose spatial scales are directly those that interest us (several meters) is discussed. However, the geometry is simplified in order to be able to make a thorough statistical analysis of the behavior and structure of the spherical flames propagating within an obstacle network. These obstacles are of toroidal shape in order to respect the sphericity of the problem. In this part our objective is a first application of the different models presented in the first part. It will also give a first look at the flame structures propagating in an obstacle network, whether they are aligned one behind the other or offset.

An experiment measured at INERIS and designed to study flame/obstacle interaction is finally modeled. A methane/air cloud initially at rest covers a three-dimensional module made of pipes. After ignition on one side of the cloud, a flame crosses the module. The characteristic length scale of the module is the meter. The comparison between modeling and the measurements will focus on the flame front position with time in its main direction of propagation, based on a movie of the flame and on the pressure field. This latter is experimentally described by a few sensors located inside and outside of the flammable cloud that recorded pressure signals.

Introduction Cléante Langrée

University of Rouen -CORIA February 8, 2022

Chap. 2 | Theoretical basis for the numerical simulation of turbulent combustion

As previously explained, UVCE are accidental scenario where a dispersed gas cloud is lit and burnt, causing thermal and pressure effects to the surrounding environment.

In order to study and then model this phenomenon, several theoretical basics need to be detailed. In the first part of this theoretical part, the focus will be put on theoretical description of the physical phenomenon at play. In a second part, the existing models for numerical simulations of turbulent combustion will be presented and discussed. A special focus will be put on the ability of the different models to take into account our application case specificity.

Physical description

On the first part of this theoretical chapter, the main physical phenomenon at play during an UVCE will be described. First combustion, turbulence and their interaction will be described and then the global phenomenon of UVCE will be detailed step by step.

Combustion

Chemical description

Combustion is a physical phenomenon characterized by a series of elementary chemical reactions that can be summarized by a global irreversible chemical reaction of the form :

Fuel + Oxidizer -→ Product(s) + (Heat) (2.1)
As recalled in the well known fire triangle (see Fig. 2.1), combustion needs three elements to happen : An oxidizer, a fuel and an energy source. When these three elements are present at the same instant, the combustion reaction can happen. From this description two types of combustion can distinguished : Theoretical basis for the numerical simulation of turbulent combustion

• When the fuel and oxidizer are mixed before the injection of the energy in the system that will trigger the combustion reaction, the combustion is called premixed. In this case, the premixed flame, will be the interface between the "fresh" mixture and the burnt product of the reaction. • When the fuel and oxidizer are not mixed before the chemical reaction occurs, it's diffusion combustion. In this case, the diffusion flame will be the interface between the oxidizer and the fuel. In both cases, the dynamic of the flame front will be dictated by the flux of the two phases : burnt and unburnt for a premixed flame, oxidizer and fuel for a diffusion flame.

As previously mentioned, an UVCE is by definition, the propagation of a premixed flame front. The focus will then only be put on premixed combustion.

For premixed combustion, under the influence of an initial energy input (spark, autoignition), a mixture of fuel and oxidizer (also known as fresh gases), will react together to form different combustion products (burnt gases) and heat. As a result, the combustion reaction is classified as exothermic. In the case of premixed combustion, several conditions are required for the start and maintenance of the combustion reaction. The proportion of fuel and oxidizer in the mixture must be within the gas mixture's flammability limits. The nature of the fuel, the pressure and the temperature all influence the upper and lower limits of a quantity we'll call equivalence ratio. The latter being the product of the mass fractions of the fuel and the oxidizer with the stoechiometric mass ratio. When the reaction is complete, the latter corresponds to the ratio between the mass of oxidant and the mass of fuel. This is the amount of oxidizer needed to burn a given amount of fuel.

When a hydrocarbon is burnt in the presence of air, the equation becomes in the stoechiometric case :

C X H Y + � X + Y 4 �� O 2 + 3.76N 2 � → XCO 2 + Y 2 H 2 O + 3.76 � Y 4 + X � N 2 (2.2)
Cléante Langrée University of Rouen -CORIA February 8, 2022 The physical phenomenon of combustion is the succession of elementary chemical reactions. A mixture of fuel and oxidant, also called fresh gas, will react together under the effect of an initial and then self-sustained energy input to form different combustion products (burnt gases) and heat.

A laminar flame is a reactive front in which an oxidation-reduction reaction takes place between a fuel and an oxidizer. This reaction can be written in the form of a simple equation but in reality can be decomposed into several elementary reactions. Several reactive schemes exist to account for this combustion reaction. Going from a simple equation to the most complete reactive scheme, for example, the GRI 3.0 mechanism of methane combustion comprising 325 elementary reactions involving 53 different species. It can be noted that for a simple one-dimensional simulation of laminar flame, the choice of the reactive scheme modeling the reaction has a very strong impact on the evolution of the flame velocity as a function of the local equivalence ratio.

Structure of a premixed laminar flame front

The surrounding aerothermochemistry drives the structure of a premixed flame. It drives it through mass and heat transport, convection and chemical reactions and kinetics. Mallard and Le Chatellier's theory [START_REF] Mallard | Recherches expérimentales et théoriques sur la combustion des mélanges gazeux explosives[END_REF] described it as follows : The heat and reactive species produced locally by the reaction diffuse towards the fresh gases. Thus, the flame propagates step by step, gradually bringing the adjacent fresh gases to the conditions of ignition. The intermediate zone between the fresh and burnt gases can be separated into two zones : the preheating zone and the reaction zone (see Figure 2.2). In the first one, the heat flow from the burnt gases raises the temperature of the fresh gases. In this preheating zone, the chemical reactions can be considered negligible compared to the thermal diffusion effects. On other hand, in the reaction zone of thickness δ r , the heat release peak is due to chemical reactions (Eq. 2.1). These two zone leads to the post-chemical equilibrium zone. In practice, two distinct physical phenomena can drive the diffusion of species in the flame front : Theoretical basis for the numerical simulation of turbulent combustion

• The Fick diffusion is due to concentration gradients,

• Soret diffusion is due to temperature gradients.

The Fick diffusion mode is dominant in most combustion phenomena and the Soret diffusion mode is only significant in front of the Fick diffusion mode for very light molecules (e.g. H 2 ) or very heavy molecules subjected to a strong temperature gradient.

Characteristic speeds of a premixed combustion

The concept of velocity is crucial to understanding premixed combustion phenomena. There are a variety of definitions and methods for measuring it experimentally [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]. Local velocities, which are linked to a flame front isotherm, and global velocities, which are linked to the spatial profile of a species produced or consumed in the flame front, can both be defined. These velocities are related to the flame kinematic properties or the integral of the reaction rate of a reactive or produced species, respectively.

In the context of our study, laminar flame speed is defined as the velocity at which the flame consumes the unburnt gases, in a 1D premixed flame. In other word, in a 1D premixed, S u is the speed at which the flame front propagates, in the reference frame of the unburnt gases. Outside the 1D unstretched laminar flame framework, defining flame speed can be more difficult though. For example, a 2D planar flame in pre-existing flow will have more characteristic identifiable "speeds", more or less related to each other.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

2.1 Physical description 19

• The absolute flame speed S a is the scalar product between the flame surface normal vector and the "displacement" vector (called w in fig. 2.3) between the position of the same point on the flame surface at two consecutive instant. It is the speed at which the flame front propagates in a fixed reference frame. • The displacement flame speed S d defined as the speed at which the flame front propagates with respect to the flow. It is therefore related to the absolute flame speed :

S d = S a -u • n (2.3)
The velocity S u is affected by flame front deformations such as stretching (more detailed therein after in a dedicated paragraph). The notations S 0 u will be used to designate the laminar propagation and burning speed of a plane, adiabatic and non stretched flame. This characteristic speed S 0 u however is not trivial to obtain. It can be theoretically derived from the conservation equation applied to a 1D premixed flame front. Such mathematical derivation [START_REF] Williams | Combustion Theory[END_REF] using asymptotic analysis can lead to the conclusion that laminar flame speed is proportional to the square root of the thermal diffusivity and Arrhenius pre-exponential factor. However, experimental configuration have been designed in order to study and determine laminar flame speeds of various fuel, such as the study of spherical expanding flames [START_REF] Varea | Experimental analysis of laminar spherically expanding flames[END_REF].

Premixed flame front thickness

The premixed flame thickness δ can be defined in various ways. The two most commonly used definitions are depicted in Figure 2

Theoretical basis for the numerical simulation of turbulent combustion

A first definition was given by Zeldovich [START_REF] Zeldovich | The theory of combustion and detonation[END_REF] :

δ p = λ ρ u C p S 0 L = D t h S 0 L (2.4)
where D th denotes the thermal diffusivity of fresh gas, such as :

D th = λ ρ u C p (2.5)
Nonetheless, this definition of flame thickness falls short of experimental measurements and it necessitates knowledge of the thermal properties of the fresh gas D th and the laminar burning rate S 0 L .

The second most used is the flame front's thermal thickness, which is defined by Spalding [START_REF] Spalding | Some fundamentals of combustion[END_REF] as :

δ th = T b -T u |∇ (T ) | max (2.6)
Where T b and T u are the temperatures of the burnt and unburnt gases, respectively. Compared to δ p , this definition is closer to experimentally available flame thickness measurements. It characterizes the heat release zone, yet it requires knowledge of the temperature profile.

Finally, note however, that it is possible to define a so-called "total" flame thickness, δ tot , as the interval between which the reduced temperature (

T -T u T b -T u
) is between 0.01 and 0.99. However this quantity has no relation to the species or temperature profiles across the flame front.

Flame stretch

The stretch represents the rate at which the flame area grows locally as a function of time. It allows to characterize the combined effect of the curvature and the flow's tangential stresses. Its definition reads [START_REF] Williams | A review of some theoretical considerations of turbulent flame structure[END_REF] :

κ = 1 δ A ∂δ A ∂t (2.7)
with δ A an infinitesimal surface element.

The deformation imposed by the total stretch is schematically shown in Figure 2.5 on an infinitesimal element δ A of a surface A propagating with the velocity S f towards fresh gas at velocity u.

The total stretch can be decomposed into two components based on the geometry of the flame and the flow using a kinematic approach [START_REF] Matalon | On flame stretch[END_REF] :

κ = κ curvature + κ shear = S f • nC + κ shear = S n f • C + κ shear (2.8)

Cléante Langrée

University of Rouen -CORIA February 8, 2022 where n is the exponent, C is the curvature-related factor and K shear is the shear-related factor and are defined as :

C = ∇ • n = - � 1 R 1 + 1 R 2 �
(2.9)

K shear = -n • Ē • n = (u • n)C + ∇ τ • u τ (2.10)
Where 1/R 1 and 1/R 2 are two characteristic curvature radii, R 1 being the largest curvature radius and R 2 being the curvature radius obtained in the normal direction n. Lastly, Ē stands for the stress tensor, which can be written as 1 2

� U + U T � .
The first term on the right-hand side of Eq. (2.10) corresponds to normal stresses, while the second term corresponds to tangential stresses.

Finally, the total stretch can be written as : [START_REF] Candel | Flame Stretch and the Balance Equation for the Flame Area[END_REF] used the transport equation of the normal unit vector n on an element of surface A to demonstrate these two contributions to the total stretch. As a result of analyzing the terms in Eq. 2.11 and Fig. 2.6, the flame responds to stretching in the following way : Theoretical basis for the numerical simulation of turbulent combustion

κ = S n f ∇ • n + (u • n)∇ • n + ∇ τ • u τ = (S n f -n n )∇ • n + ∇ τ • u τ = (S d )∇ • n + ∇ τ • u τ (2.11) Candel & Poinsot
• The normal propagation can either increase or decrease the flame area depending on the positive or negative curvature. • The velocity projected in the normal direction n of the fresh gas can be positive or negative depending on the curvature. • The tangential component of the positive or negative velocity diverges, resulting in a tangential force that either extends or decreases the flame surface.

Figure 2.6 -Curvature effects on flame front. From [START_REF] Lefebvre | Analyses théorique, numérique et expérimentale de la détermination de la vitesse de combustion laminaire à partir de flammes en expansion sphériques[END_REF].

Instabilities of laminar flame front

Laminar flame fronts, considered as interface between burnt and unburnt mixture, are subject to instabilities. Different types of instabilities can be identified and related to a physical phenomenon affecting the laminar flame front and its propagation : gravity, thermal diffusivity etc...

Gravitational instabilities

When two fluids of different densities are superimposed, gravity instabilities occur. As shown in Figure 2.7, they cause fluid motion that favors or opposes the movement of the flame front. When the cooler, denser gases, are above, flame propagation is favored because they are drawn to the flame front moving in the opposite direction. The gravitational effect stabilizes the situation. If the fresh gases are below, on the other hand, gravity pulls the flame front in the direction of propagation of the flame front. The gravity effect is destabilizing in this case. In fact, these two limit cases can occur in the case of a spherically expanding flame. Fresh (denser) gases are located above hot (lighter) burnt gases in the upper half-sphere. In the opposite direction of the gravity field, the flame front propagates upwards. As a result, the fresh gases are drawn downwards, directly meeting the flame front. Therefore, this configuration is stable. The fresh gases are placed below the burnt gases in the lower half-sphere. Fresh gases will "leak" from the flame front as it propagates in the same direction as gravity.

Nonetheless, the characteristic times required for the establishment of these instabilities are very long in comparison to those required for chemistry or flame front propagation. Only flames with slow propagation speeds, i.e. those on the verge of extinction, are susceptible of such a behavior. Gravity has no effect on burning velocities greater than 15 cm/s, according to an experimental study conducted by Ronney and Wachman [START_REF] Ronney | Effect of gravity on laminar premixed gas combustion I: Flammability limits and burning velocities[END_REF].

Hydrodynamic instabilities

Darrieus [START_REF] Darrieus | Propagation d'un front de flamme[END_REF] and Landau [START_REF] Landau | On the theory of slow combustion[END_REF] were the first to formalize knowledge of hydrodynamic instabilities in combustion. The Darrieus-Landau (or hydrodynamic) instabilities, named after them, result from the rapid acceleration of gases across the flame front. In fact, fresh gases have a much higher density than burnt gases (a ratio of about 7 for most hydrocarbons). The speed at which the gases cross the flame front is then multiplied by the ρ b ρ u ratio, known as the expansion factor. An infinitely thin flame front, which is initially flat and undergoes a small amplitude perturbation at a given time, exhibits wrinkles distinguished by their curvature (by convention, the curvature is negative at a concave flame front with respect to fresh gas and positive in convex areas). Combustion instabilities can occur as a result of streamline deflection at the flame front. Indeed, a reactive flow is distinguished by a sharp difference in density between Theoretical basis for the numerical simulation of turbulent combustion fresh and burnt gases. As the streamlines approach a convex zone, they diverge to ensure mass conservation and the velocity of the fresh gases decreases locally. When approaching a concave zone, on the other hand, the streamlines converge and the velocity of the flow increases, while the unstretched laminar combustion velocity remains constant (see Figure 2.8). As a result, the wrinkles caused by a disturbance continue to grow, making the flame front intrinsically unstable. These instabilities become more prominent as the laminar flame thickness decreases (i.e., when the density gradient increases).

Thermo-diffusive instabilities

In addition to the Darrieus-Landau instabilities caused by the expansion of the burnt gases, instabilities caused by the imbalance between thermal and mass diffusivity across the flame front must be considered for slightly wrinkled flames (Figure 2.9). To characterized this imbalance, the Lewis number Le, ratio of thermal over mass diffusivity, is introduced. When the species diffusion is weaker than the heat diffusion (Le > 1), the temperature drops behind the convex part of the flame front due to a decrease in mass diffusion of the minority species and an increase in heat loss, resulting in a local decrease in combustion rate. The opposite phenomenon is observed in the concave parts of the flame: the local combustion temperature rises, resulting in a local increase in combustion speed. As shown in Figure 2.9, these two effects tend to smooth the flame. When the Lewis number is less than one, the opposite phenomenon occurs, resulting in the formation of surface asperities on the flame front.

In the case of an expanding flame, the time it takes for cells to appear is proportional to the Péclet number, which is the ratio of the flame's radius to thickness (P e = r δ ). When the radius is large enough, cell structures appear when the Péclet number exceeds a critical value (P e crit = r crit δ ). Bauwens et al. [START_REF] Bauwens | Experimental study of sphericalflame acceleration mechanisms in large-scale propane-air flames[END_REF] and Kim et al. [START_REF] Kim | Effect of propagation behaviour of expanding spherical flames on the blast wave generated during unconfined gas explosions[END_REF] studied the growth of a flame in a large volume in order to see the full development of hydrodynamic instabilities. Kim et al. demonstrated that, regardless of the mixture, the critical Péclet (P e crit ) increases with the Markstein number (Ma). The same trends were observed in flame speed measurements for propane/air mixtures by Bauwens et al [START_REF] Bauwens | Experimental study of sphericalflame acceleration mechanisms in large-scale propane-air flames[END_REF]. They emphasized the oscillatory nature of the velocity acceleration after the onset of cells in their work. They interpreted these oscillations as the result of a cycle of cell growth and subsequent stabilization. The first cells will form on the flame's surface and then grow. When the cells have grown largely enough, new smaller cells will form on the surface. They then used the fractal theory developed by Gostintsev et al. [START_REF] Gostintsev | Self-similar propagation of a free turbulent flame in mixed gas mixtures[END_REF] to increase the flame surface, i.e. a power law in which the radius is proportional to the flame propagation time raised to the power : r ∝ t n with n = 1.2.

Obstacle induced instabilities

The propagation of a flame in an enclosure generates acoustic waves, which can interact with the flame front via multiple instability mechanisms after reflection from walls or obstacles. In the presence of obstacles or a complex geometry, strong instabilities such as Kelvin-Helmholtz instabilities (shear instabilities), Rayleigh-Taylor instabilities (density difference instabilities), or Richtmyer-Meshkov instabilities develop (Rayleigh-Taylor instability perturbed by a shock wave). Besides, fast flames can produce shock waves, which can reflect of surfaces and interact again with the flame.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

2.

.9 -Thermo diffusive instabilities. From [START_REF] Goulier | Comportements aux limites de flammes de prémélange hydrogène/air/diluants. Étude de la transition flamme laminaire-flamme turbulente[END_REF].

When compared to the instabilities caused by geometry, the Darrieus-Landau and thermo-diffusive instabilities have a minor impact (shock wave, fluid dynamics). When the flame speed is low, these instabilities prevail. The Kelvin-Helmholtz, Rayleigh-Taylor and Richtmyer-Meshkov instabilities have a significant impact on flame development and are the primary causes of flame acceleration when the flame propagates through an obstruction.

Turbulence

Turbulence characterizes a flow in which the mean motion is disturbed by dissipative structures. These latter make the behavior of the fluid unpredictable and are characterized by an energy spectrum over a wide range of scales . As it is well known, they are linked to the Reynolds number, which compares inertial forces to viscous forces and characterizes the transition from a laminar to a turbulent regime :

Re = | u | l ν (2.12)
Where l and u are the characteristic length and speed of the considered flow, respectively, and ν being the kinematic viscosity.

Theoretical basis for the numerical simulation of turbulent combustion

Turbulent flows are distinguished by large variations in the velocity field in both time and space. The presence of a continuous spectrum of vortex structures convected by the mean flow is an essential component of turbulent flows. These vortices interact strongly with one another through an energy cascade process that increases mass, momentum and heat transfer when compared to a laminar flow. Following the Kolmogorov theory [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers[END_REF][START_REF] Kolmogorov | Dissipation of energy in the locally isotropic turbulence[END_REF], there are three distinct areas that can be identified (see Fig. 2.10) :

• The integral zone has the lowest wave numbers and the largest and most energetic structures associated with the integral length scale l t . The structures in the integral zone have length scales l t and velocity scales u � (l t ) that are comparable to those used to define the Reynolds number of the flow and are not affected by viscous effects. This zone, which exhibits inhomogeneous and anisotropic behavior, is primarily controlled by the geometric characteristics of the flow. To characterize this zone, the turbulent Reynolds number Re t is introduced :

Re t = u � (l t )l t ν (2.13)
• In the second zone, large vortices become unstable and split into smaller vortices due to the process of "energy cascade". Because of viscous forces, there is no energy dissipation effect here. The energy is transferred directly from the larger structures to the smaller ones, according to the famous k -5/3 law for isotropic homogeneous turbulence. At different scales, the flow of energy from one scale to another is constant.

It is given by the turbulent kinetic energy dissipation rate. This dissipation rate � is estimated for a scale length r as the ratio of the turbulent kinetic energy u � 2 (r ) and the time scale r /u � (r ) associated with the vortex :

� = u � (r ) 2 d r u � (r ) = u � (r ) 3 r (2.14)
• The dissipative zone: it has the highest wave numbers and its scale range is limited by the Kolmogorov scale l K , which is the smallest vortex scale in the flow. It is the scale at which inertial and viscous forces compensate for each other, so that epsilon dissipation converts turbulent kinetic energy k into heat due to the fluid's kinematic viscosity ν. The Kolmogorov length scale l K and characteristic velocity u K are calculated as follows :

l K = � ν 3 � � 1 4
(2.15)

u K = (ν�) 1 4
(2.16)

Turbulence -Flame interaction

Turbulent combustion can be represented as the interaction of a laminar flame and turbulence. Indeed during its propagation the flame front will interact with its environment and especially the turbulence of this environment. This interaction can be characterized using the scales of laminar flames : its laminar flame velocity and its laminar flame thickness and to the scales of turbulent structures : the turbulent kinetic energy and the turbulent dissipation (via which a spatial scale can be characterized for the turbulence).

Cléante Langrée

University of Rouen -CORIA February 8, 2022 Turbulent structures mainly have two effects on the flame front: wrinkling and thickening. The nature of this effect will depend on the scale ratio between the local turbulence in contact with the flame front and that of the flame front. Larger turbulent structures will wrinkle the flame front and smaller ones will thicken it.

The propagation of the flame front, via the thermal expansion of the fuel gas, will generate a flow of fresh gas upstream of the flame front. This flow, associated with possible obstacles, will naturally generate turbulent intensity by its acceleration upstream of the flame front propagation.

By taking into account these different effects, flame-turbulence interactions can be categorized in relation to the space and velocity scales of the flame front and turbulence. Two dimensionless numbers can be introduced to characterize this interaction. The Damkholer number defined as the ratio of the turbulent and chemical time scales and the Karlowitz number defined as the ratio of the chemical and Kolmogorov time scales. The classification of flame-turbulence interaction regimes are represented schematically in a Borghi diagram [START_REF] Borghi | On the structure and morphology of turbulent premixed flames[END_REF], as depicted in Fig. 2.11, which identifies 4 main flame-turbulence interaction regimes.

In the case of a UVCE the turbulent flame interaction can be characterized. The turbulent scales with which the flame interacts are multiple: atmospheric turbulence, jet turbulence and drag turbulence (due to obstacles). These different scales will evolve during the propagation of the flame front due to the acceleration of its propagation. Nevertheless, a characterization of the flame turbulence interaction, mainly based on the scales of atmospheric turbulence, had been proposed by Proust [START_REF] Proust | Formation, inflammation, combustion des atmosphères explosives (atex) et effets associés[END_REF] for the situation of industrial UVCE. It is displayed in the Borghi diagram presented below.

Unconfined Vapour Cloud Explosion (UVCE)

As previously said, an Unconfined Vapour Cloud Explosion occurs when a flammable mixture is spread from a source, forms an explosive cloud and explodes. A flammable mixture is defined as a mixture of a fuel and an oxidant, the two reactants necessary for the combustion reaction to occurs. Eq. 2.2 is the combustion equation of hydrocarbon as fuel and air as oxidant. As stated earlier, for the mixture to be flammable, fuel and oxidant have to be under specific conditions. Under the Lower Flammable Limit (LFL), there's isn't enough fuel, above the Upper Flammable Limit (UFL) there is too many fuel. Some examples of flammability limits for various fuel are represented in figure 2.12. The first step of an UVCE is the formation of a flammable cloud across the plant. It begins with the loss of confinement of a flammable product (see Fig. 2.13). Multiple phenomenon can lead to such a loss :

0
• Breach in a pressurized pipe/tank; • Overflowing of a storage.

Note that depending on the timing of the unwanted events occurring during a loss of confinement, it can lead to several accidental scenarios, including pollution, pool fire, jet fire etc. This loss of flammable product can be two-phase, or fully gaseous or liquid. The formation of the flammable cloud that will "explode" happens when the gaseous part and/or the evaporating liquid part mixture with surrounding air to form a cloud within the inflammability limits. It is important to note that the type of accidental scenario leading to the formation of the flammable cloud has a great impact on the UVCE itself. The impact on the cloud geometry, internal turbulence and concentrations heterogeneity can lead to very different UVCE depending on the context of the fuel release.

Flammable mixture ignition

Once the flammable mixture is spread, the chain reaction which is the UVCE begins with the ignition of the mixture at one point in the cloud. This ignition will give the amount of energy needed to engage the combustion reaction at the given location, if the fuel and air at the given location are inside the flammability limits. And then, the reaction will give, via exothermic heat, enough energy to the surrounding part of the cloud for the combustion reaction to happen. In an industrial context, several energy sources can ignite the flammable mixture : Theoretical basis for the numerical simulation of turbulent combustion

• A heated surface (from any heating machine of an industrial process)

• An electric arc • A sparkle (from the collision between to metallic objects for instance)

The ignition of a flammable mixture is in itself a complicated question which has been a full subject of research for years now. In the case of our study, as the main point of interest is predicting the consequences of an UVCE, the ignition will be considered as the starting point of the flame front propagation, without considering any failed, or too complex cases of inflammation.

Flame front propagation

Once the mixture is lit, the combustion reaction spreads across the flammable cloud in the form of the propagation of a flame front. The flame front being a diffusive interface between the cold fresh gases (flammable mixture) and hot burnt gases (combustion products).

Flame front starts propagating laminarly, which means that the flame front will be nearly unperturbed. At this point the flame front is propagating at laminar flame speed, which is only affected by the chemical kinetic of the reaction and the curvature/stretching of the said front. Very quickly, the unperturbed flame front starts being affected because of intrinsic flame instabilities. Those instabilities are only dependent on the nature of the fuel/air mixture and do not depend on the environment (like its turbulence). Such flame fronts are called "cellular" flame fronts and depending on the nature of the fuel/air mixture, they can vary in stability and in propagation speed.

The flame front will also interact with turbulence, on time scales depending on the turbulence itself (very quickly for high turbulent flows, quite slowly for low turbulent flow). By interacting with the turbulence, the flame front will be called turbulent and its propagation speed will be higher than its laminar flame speed propagation.

During all those transformations, one important rule will be at play. The flame front propagation speed is proportional to its surface :

S T S U = A T A (2.17)
With S T and S U being the turbulent and laminar flame speeds respectively and A T and A being the surface of the perturbed and unperturbed flame front respectively. In fact, most of macro-scales phenomena leading to an increase in flame front surface, will be considered as accelerating phenomena.

Cléante Langrée

University of Rouen -CORIA February 8, 2022 Theoretical basis for the numerical simulation of turbulent combustion

Flame front propagation

When a flame spreads through a flammable mixture, two mechanisms amplify the pressure increase. First, the flame propagation speed is directly linked to the pressure induced by the propagation. Deshaies and Leyer [START_REF] Deshaies | Flow field induced by unconfined spherical accelerating flames[END_REF] and Leyer [START_REF] Leyer | Effets de pression engendres par l'explosion dans l'atmosphere de melanges gazeux d'hydrocarbures et d'air[END_REF] gave a general description of the pressure evolution in which ∆P is directly correlated to the flame propagation speed (

dr f (τ ) dt )
and its acceleration (

d 2 r f (τ ) dt 2 ) : ∆P (r, t) = ρ 0 � 1 - 1 α � � 2 r f (τ ) r � dr f (τ ) dt � + r 2 f (τ ) r . d 2 r f (τ ) dt 2 � (2.18)
With r f being the flame radius, ρ 0 the atmosphere density, α the burnt gases rate of expansion (i.e. ratio of unburnt over burnt density) and τ = tr c 0 , with c 0 the sound velocity in the atmosphere.

Second, the containment of the flammable mixture plays a role as well. If the flammable mixture (and therefore the deflagration) is confined, even partially, the pressure can build-up inside the confinement zone, and then releases in an instant if the containment is broken. In the most serious accidents, a combination of these two mechanisms is usually responsible for the overpressure generated.

Flame front acceleration

In the cases of UVCE occurring on an industrial site, several phenomenological accelerating factors can be identified. First the interaction between flame and turbulence is a well known accelerating factor (fig. 2.15) [START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Lee | Validation of analytical expressions for turbulent burning velocity in stagnating and freely propagating turbulent premixed flames[END_REF]. In UVCE scenario, mainly two types of turbulence will be at play : the atmospheric turbulence of the industrial environment, and the turbulence of the flow upstream of the flame front induced by the flame propagation.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

2.1 Physical description 33 Some significant phenomena in UVCE applications are the effects of the obstacles on the flame speed. Obstacles, and more generally, obstructions in the path of a propagation are another well known accelerating factor [START_REF] Zeeuwen | Experimental investigation into the blast effect produced by unconfined vapour cloud explosions[END_REF][START_REF] Van Wingerden | Flame propagation in the presence of repeated obstacles: influence of gas reactivity and degree of confinement[END_REF]. Obstacles can accelerate the propagation of a flame front by deforming it and thus increasing its surface. Another effect, is the turbulence induced from interaction with the unburnt gases flow generated by the flame front propagation.

Mixture inhomogeneities between the fuel and oxidizer has been shown experimentally to be another factor that increases the flame propagation speed [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][START_REF] Zhou | Flame propagation characteristics in a heterogeneous concentration distribution of a fuel-air mixture[END_REF]. Most common inhomogeneities encountered in UVCE scenarii are stratified mixture. Such stratified mixture are likely to form in the case of a cloud formation by evaporation, or even in any leak scenario (cf fig. 2.13) if the cloud is left unlit long enough so that a stratification appears naturally by gravity. Inhomogeneities are also likely to appears within a flammable cloud formed by a highly turbulent jet (in breaches scenario, cf fig. 2

.13).

Flame propagation regime

Depending on the explosion conditions (confinement, obstacles, fuel type. . . ), two propagation regimes are possible :

• Deflagration : "Autonomous subsonic propagation mode of the reaction in a combustible environment (ideally pre-mixed) thanks to its coupling with the heat and material transport mechanisms. The deflagration mode is characterized by a decrease of the pressure and the density at the same time as an acceleration of the gases crossing the reaction zone." • Detonation : "More or less autonomous propagation of a combustion zone coupled to complex shock waves which precede it, taking place with a speed higher than the speed of sound in relation to the reactive medium. The detonation mode is characterized by an increase in pressure and density at the same time as a deceleration of the gases with respect to the reaction zone through which they pass. In this mode, the propagation speed of the reaction zone must be supersonic with respect to the fresh gas."

Due to the level of complexity that detonation brings in numerical simulation, mainly with the numerical management of the induced shock, this study will focus purely on deflagration.

Deflagrations, main subject of our investigation, are complex multi-physical and multi-scale phenomena involving chemical reactions, species and heat transport, turbulence, flame/obstacle interaction, acoustics and so on, all of which are strongly coupled. Each of these phenomena must be correctly reproduced when studying explosions numerically if the parameters of critical interest are to be reproduced : either directly or by modeling. The overpressure field generated by the propagation of the flame front is very different from one combustion regime to another. Therefore, when the flame front propagates at a speed inferior to the speed of sound (i.e. deflagration), the overpressure generated via the flame front propagation begins to attenuate itself via its propagation towards infinity (i.e. atmospheric pressure). On the other hand, when the flame front propagates at a speed superior to the speed of sound, the overpressure is "pushed" without having time to attenuate, as it propagates at an inferior speed towards infinity. Typical pressure profile generated by a flame front propagation is shown in figure 2.16. From [START_REF] Proust | Formalisation du savoir et des outils dans le domaine des risques majeurs (eat-dra-76) les explosions non confinées de gaz et de vapeursω uvce[END_REF].

Theoretical basis for the numerical simulation of turbulent combustion

Transition from deflagration regime to detonation regime are often referred as DDT. Various parameters can affect occurrences of DDT, such as the obstruction levels or the type of fuel. As a major turning point in an UVCE scenario, DDT are extensively studied, both experimentally and numerically. In our study the focus will be put on deflagrations only, because studying detonation would bring too much complexity (shock capturing, shock-flame interactions etc...). In deflagration cases the reaction rate is slower than in a detonation case. The head wave is not a shock wave and is followed by a depression which can also cause damages. Considering a spherical system, the pressure field generated by the flame front propagation decreases with the radius.

Modeling

Now that the theory behind UVCE has been described, the next step is to model the different physical phenomena involved.

Reactive flow : Navier-Stokes equations

The Navier-Stokes equations are the classic equations to represent the conservation of mass, momentum and energy of a fluid. In the case of a reactive flow, an equation for the conservation of mass species fraction is also present, to take into account the species transformation induced by the combustion reaction.

∂ρU ∂t + ∇ • (ρUU) = -∇p + ∇ • τ (2.

20)

Species mass fraction equation

∂ρY k ∂t + ∇ • (ρUY k ) = -∇ • J k + ωk (2.

21)

Enthalpy equation

∂ρH ∂t + ∇ • (ρUH) = ∂p ∂t + ∇ • � J h + τ • U � + Q (2.22)
The above system (2.19-2.22) of partial differential equations (PDE) must be associated with the definitions of the unknown variables. First, the shear stress tensor τ for a Newtonian fluid reads :

τ = µ∇U + µ∇U T -µ 2 3 (∇ • U)δ ij (2.23)
To get the laminar viscosity, the Sutherland transport model which expresses the viscosity as a function of the temperature and two model constants A s and T s :

µ = A s √ T 1 + T s T (2.24)
The laminar enthalpy and species flux are expressed with a Fick's law, assuming that the laminar flows are only caused by species and enthalpy gradients.

J k = µ Sc k ∇Y k (2.25) J h = µ P r ∇h + N � k=1 � P r Sc k -1)h k ∇Y k � (2.26)
Several non-dimensional numbers are formulated in those equations, the Schmidt number (ratio of momentum to mass diffusivity), the Prandlt number (ratio of momentum to thermal diffusivity) and the Lewis number (ratio of thermal to mass diffusivity) :

Sc k = µ ρD k (2.27)
Theoretical basis for the numerical simulation of turbulent combustion

P r = µC p λ (2.28) Le k = Sc k P r = λ ρC p D k (2.29)
Finally, Q represents the chemical source term, i.e. the heat release by the combustion reaction.

Q = - � i h i ωi (2.30)
To invert enthalpy in order to retrieve the temperature the janaf model is used [START_REF] Burcat | Thermochemical Data for Combustion Calculations[END_REF]. With the Janaf thermophysical model, the temperature is computed by solving a c P = f (T ) equation, knowing the heat capacity at pressure constant and the model constants a n :

c P = R � � � (a 4 T + a 3 ) T + a 2 � T + a 1 � T + a 0 � (2.31)
Finally the state equation is used, linking pressure and density :

p = ρ R � M T (2.32)

Turbulence modeling

Turbulent energy spectrum : RANS/LES/DNS Modeling turbulence in Navier-Stokes equations can be done at several levels (fig. 2.17). Resolution methods are called DNS (Direct Numerical Simulation) when the Navier-Stokes equations are exactly solved. In DNS, all the scales of turbulence that can be captured by the mesh are captured . In our context, DNS is too costly in terms of mesh refinement and computational resources, and thus will be avoided in our investigations.

Another approach consists in modeling the turbulent movement of the flow on a coarser mesh compared to DNS. This is called RANS (Reynolds Average Navier Stokes). Rather than trying to solve all the turbulent scales via the direct resolution of the equations, the average fields of interest are solved and the fluctuating part is modeled in order to consider the turbulence of the flow. Note that the fields of interest are averaged over the physical simulation time. However, there is a derived URANS (Unsteady Reynolds Average Navier-Stokes) approach, which allows an unsteady RANS simulation.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

Figure 2.17 -Turbulent spectrum representation of each turbulence modeling method.

Figure 2.18 -Classification of turbulence modeling methods. Adapted from [START_REF] Sagaut | Multiscale and multiresolution approaches in turbulence-LES, DES and Hybrid RANS/LES Methods: Applications and Guidelines[END_REF].

Theoretical basis for the numerical simulation of turbulent combustion

An intermediate approach consists in the resolution of large scales of turbulence, easier to capture and in the modeling of smaller scales, more homogeneous between them and complicated to solve (because involving a finer mesh by definition). This is called LES, or Large Eddy Simulation. The LES solved equations are said to be filtered, because a filter is applied to solve/model turbulent structures according to their spatial scales.

There are variations and different approaches of these 3 families of resolution, nevertheless it forms a continuum of resolution approach varying in precision and cost of calculation (fig. 2.18). Note that the averaged and or filtered equations are mathematically similar, with the only difference being the formalism differentiating the terms filtered/solved for LES and average/fluctuating for RANS.

Favre averaged Navier-Stokes equations

Reynolds Averaging

In order to average the Navier-Stokes equations, a Reynolds operator is applied to the fields of interest of the equations: velocity, density, energy variable and mass fractions of the species. This leads to the following :

φ = 1 T � T φ(x i , t) dt (2.33)
with

Φ = Φ + Φ � (2.34) Φ � = 0 (2.

35)

Favre Averaging

In our case, the density cannot be considered constant, and a different kind of averaging operator, called Favre averaging need to applied. It is defined as follows :

� φ = ρΦ ρ = 1 ρT � T ρ(x i , t)φ(x i , t
∂ρ ∂t + ∇ • � ρ � U � = 0 (2.40) ∂ρ � U ∂t + ∇ • � ρ � U � U � = -∇ • � ρ � U �� U �� � -∇p + ∇ • τ (2.41) ∂ρ � Y k ∂t + ∇ • � ρ � U � Y k � = ∇ • � ρ � U �� Y �� k � + ∇ • J k + ωk (2.42) ∂ρ � H ∂t + ∇ • � ρ � U � H � = ∂p ∂t + ∇ • � ρ � U �� H �� � + ∇ • � J h + τ U � + Q (2.43)
This decomposition into mean and fluctuating parts brings new terms into the equation. It is these new terms, not directly solvable, that will be modeled and that will account for the turbulent character of the flow.

Several unknown terms remain:

• the Reynolds stress tensor : � U �� U ��

• the species and enthalpy turbulent flux

: � U �� Y �� k � U �� H ��
• the species chemical reaction rate : ωk

As far as turbulence is concerned, two main terms are to be modeled. These are the turbulent transport terms and the turbulent stress tensor or Reynolds tensor. In order to model turbulent transport terms, it is common to use a gradient diffusion hypothesis or Fick's law. This approach consists in modeling turbulent transport as a function of the gradient of the transported field. Note that in some specific cases, counter-gradient diffusion can be observed. However, in our study this gradient diffusion law will be used to solve the turbulent transport terms.

� U �� Y �� k = - µ t S c kt ∇ • � � Y k � (2.44) � U �� H �� = - µ t P r t ∇ • � � H � (2.45) Cléante Langrée UVCE Simulation February 8, 2022 40
Theoretical basis for the numerical simulation of turbulent combustion

The Reynolds stress tensor is based on the linearity or Boussinesq hypothesis which expresses :

ρ � U �� U �� = -µ t � ∇ � U + (∇ � U) T - 2 3 ∇ � ∇ • � U � Ī� + 2 3 ρk Ī (2.46)
The turbulent viscosity ν t , corresponding to the contribution of the flow turbulence to the fluid viscosity, is involved in all these functions. The models of turbulence are presented below provide a model for this turbulent viscosity, which will allow us to solve for the turbulent transport terms and the Reynolds stress tensor.

In the case of the (first-moment closure) RANS models, the turbulent viscosity is expressed as a function of two characteristic magnitudes of turbulence. In order to characterize the turbulence of the flow, its energy and spatial scales need to be characterized. This approach is called a 2-equation approach because it is required to define two variables that will characterize the turbulence and provide the model with an equation for each of these terms. The resolution of these two equations will allow us to go back to the turbulent viscosity and thus the closure of the Navier-Stokes equations will be ensured.

Two-equations URANS models k-�

A commonly used family of models are the k � models [START_REF] Launder | Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc[END_REF]. In these models the turbulent kinetic energy k represents the energy scale of the turbulent structures and the turbulent dissipation �, represents the dissipation rate of this turbulent kinetic energy and thus indirectly, the spatial scale of the turbulent structures. The two transport equations of k and � are defined in order to solve these terms. This model is widely used at high Reynolds for its simplicity of implementation. Nevertheless, it can have difficulties accounting for turbulence in the near-wall regions.

µ t = ρC µ k 2 � (2.47) ∂ρk ∂t + ∇ • � ρ � Uk � = ∇ • �� µ + µ t σ k � ∇k � + P k -ρ� (2.48) ∂ρ� ∂t + ∇ • � ρ � U� � = ∇ • �� µ + µ t σ � � ∇� � + C �1 � k P k -C �2 ρ � 2 k P k (2.49)
This model includes several constants with usual values attributed to them : C �1 = 1.44, C �2 = 1.92, C µ = 0.09, σ k = 1.0, σ � = 1.3. A production term P k is also defined as :

P k = -ρ � U �� U �� ∇ • � U (2.50)
Cléante Langrée University of Rouen -CORIA February 8, 2022

2.2 Modeling 41 k-ω
Another family of models are the k-ω models [START_REF]Turbulence modeling for CFD[END_REF][START_REF] Wilcox | Formulation of the kw turbulence model revisited[END_REF]. These models carry turbulent k energy, like the k-� models, but use the specific dissipation rate ω rather than the � dissipation. Similar to the k-� model, solving 2 additional equations closes the Navier-Stokes PDE system. It can be noted that the k-ω model allows a better resolution in close wall, as well as for low Reynolds numbers and can also be enhanced as for instance with the Shear Stress Tensor (SST) version [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF] :

µ t = ρ a 1 k max(a 1 ω, SF 2 )
(2.51)

∂ρk ∂t + ∇ • � ρ � Uk � = ∇ • ((µ + σ k µ t ) ∇k) + P k -β * ρkω (2.52) ∂ρω ∂t + ∇ • � ρ � Uω � = ∇ • ((µ + σ ω µ t ) ∇ω) -βρω 2 + 2ρ(1 -F 1 )σ ω 2 1 ω ∇k • ∇ω (2.53)
With F 1 and F 2 being sensor function, with F 1 = 1 in near wall regions, F 1 = 0 away from the surface, F 2 = 1 for boundary layer flow and F 2 = 0 in free shear layers. This model also includes several constants with usual values attributed to them : a 1 = 0.31, σ k = 0.85, β * = 0.09, β = 0.075, σ ω = 0.5, σ ω 2 = 0.856. A production term P k is also defined as :

P k = mi n(τ ∇ (U) , 10β * kω) (2.54)

Turbulent combustion modeling

At the end of the favre averaging of Navier-Stokes equations, several terms remain unknown :

• the Reynolds stress tensor : � U �� U ��

• the species and enthalpy turbulent flux

: � U �� Y �� k � U �� H ��
• the species chemical reaction rate : ωk A model has been provided for the Reynolds stress tensor and the species and enthalpy fluxes. The last section will be about the modeling of the last term : The chemical reaction source term.

Theoretical basis for the numerical simulation of turbulent combustion

Turbulent combustion models

No model approach

The flame-turbulence interaction should represent this chemical source term in the species mass fraction conservation equation. Numerically this chemical source term represents the average over a mesh, whereas in reality the flame front represents a discontinuity that can cross the meshes. It would therefore be possible in theory to develop the reaction rate obtained via mean Arrhenius law. This modeling of the average source term, by decomposing the terms using a Reynolds mean, brings out new terms too complex to model or obtain numerically. An approximation neglecting these new terms in the equation is in most cases too imprecise, which in theory could be acceptable for a sufficiently fine flame front resolution, but which in our cases seems completely unrealistic.

Other methods to model this average source term are therefore needed. In turbulent conditions 3 main approaches can be identified to model this chemical source term : turbulent mixture type analyses, statistical analyses and geometric analyses. Note that there are also other approaches based on the mean of the Arrhenius source term seen above. Among these methods are thickened flame models [START_REF] Ducros Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF], which consist in introducing into the equation system a thickening term F , aimed at thickening the flame front in order to obtain a finer relative resolution. This manipulation allows the hypothesis seen above to be acceptable. This kind of model has several advantages, as it allows to solve each mass fraction of each species. However, this model requires a basic resolution fine enough not to thicken the flame too much and to remain on a faithful model. In our case, this model does not seem to be very suitable because at large scale, the thickening factor would become too important. Nevertheless some similar configuration studies have used this model and obtained conclusive results [START_REF] Quilattre | Simulation aux grandes echelles d'explosions en domaine semi-confine[END_REF][START_REF] Vermorel | LES of explosions in venting chamber: A test case for premixed turbulent combustion models[END_REF].

Progress variable

For the approaches covered next, the notion of progress variable needs to be introduced. In turbulent combustion, the evolution of combustion can be modeled not by following the evolution of mass fractions, but by following a progress variable. This progress variable quantifies the progress of the combustion reaction in the flammable mixture, and goes typically from 0 (unburnt) to 1 (fully burnt). This variable of progress can be defined according to the temperature (eq. 2.55) or according to the mass fraction of fuel (eq. 2.56) or oxidant. A transport equation for this progress variable can be defined and will be in fact the equation that will be solved, putting aside the equation of the mass fractions specific to each species.

c = T -T u T b -T u (2.55) or c = Y F -Y u F Y b F -Y F (2.56)
Cléante Langrée University of Rouen -CORIA February 8, 2022

Modeling 43

In premixed combustion, the progress variable propagates through the mixture. Its propagation can be modeled as follows :

∂ρc ∂t + ∇ • (ρUc) = ∇ • (α∇c) + ωc (2.57)
If eq.2.57 is averaged :

∂ρ� c ∂t + ∇ • � ρ � U � c � = ∇ • � α∇� c -ρ � U �� c �� � + ωc (2.58)

Statistical approaches

These approaches consist in defining at each point of the domain a probability density function (PDF) in order to determine the state of the mixture. In fact, by considering that the mixture can be either burnt or unburnt, the PDF can be expressed as the weighted sum of two Dirac functions centered at 0 and 1 respectively. With the help of these PDF, the fields of interest can be reconstructed for the simulation.

As depicted in Fig 2 .19, in order to determine these PDF, one can assume basic profiles (e.g. Gaussian, beta. . . ) or solve a transport equation. The mean chemical source term at position x and at time t can then be expressed as the integral of the source term as a function of driving variables times the probability density of the driving variable at position x and at time t (eq. 2.59). Models based on turbulent mixing assume that the combustion reaction is only driven by the mixture induced by turbulence. Initially developed for diffusion flames (or not premixed), the two phases considered are the fuel and the oxidizer. In this framework, the combustion reaction is considered to be infinitely fast and as soon as the reagents are mixed in proportions that allow the combustion reaction to take place, and burnt gas will be produced. In the framework of premixed combustion, a similar approach can be used. However the two phases considered here will be unburnt mixture and burnt gases. The propagation of the flame front is then only driven by the mixture between thoses two phases, with the mixing being driven mainly by turbulence. Among these models is the Eddy Break Up (EBU) model [START_REF] Magnussen | On mathematical models of turbulent combustion with special emphasis on soot formation and combustion[END_REF][START_REF] Magnussen | On the structure of turbulence and a generalized eddy dissipation concept for chemical reaction in turbulent flow[END_REF], which proposes to model the chemical source term of the transport equation of the progress variable as the product of a constant of the density of the turbulent characteristic time and the variance of the progress variable (eq. 2.60).

ωc = C EBU ρ � k � c(1 -� c) (2.60)
One can also cite the Bray-Moss-Libby analysis [START_REF] Bray | Unified modeling approach for premixed turbulent combustion-Part I: General formulation[END_REF], which in the same way proposes to model the chemical source term (eq. 2.61); this time only explaining the constant and linking it to a numerically determined mixing progress variable.

ωc = 2 2c m -1 ρ� c(1 -� c) � k (2.61)
The main disadvantage of these methods in our case is that combustion is only driven by turbulence. Nonetheless, a premixed flame front propagation starts with a laminar phase. In the absence of turbulence during this phase, the formulation of this model will not allow the flame speed to be reproduced correctly. Moreover, these models do not allow direct consideration of inhomogeneities.

Geometrical approaches

They are based on the geometric analysis of the flame front to estimate the average reaction rate. Based on the hypothesis of an infinitely fine flame front, i.e. for the interface between fresh gas and fuel gas, there are two ways to deal with this problem: the solution of a scalar flame front tracking field equation and the modeling of the source term via a flame surface density (see Fig 2.20).

The first approach, whose most well-known representative model is the G-equation [START_REF] Williams | Turbulent combustion[END_REF], which no longer solves an equation for the progress variable, but an equation for the G field instead. This one can be defined in several ways. The most common is to define it as the distance to the flame front. This modeling introduces a new chemical source term that is defined as the product of the fresh gas density, the turbulent flame velocity and the G gradient.

Cléante Langrée

University of Rouen -CORIA February 8, 2022 In order to solve all the unknowns in the model, a turbulent flame velocity correlation must be provided. This makes the flame front propagation results dependent on the turbulent flame velocity correlation used. Nevertheless, the reactive source term depends only on turbulence, since it is dependent on a turbulent flame velocity model.

Another type of geometrical analyses, is to define the flame area density Σ as the amount of flame area in a given volume (in m 2 m -3 ) (illustrated in fig. 2.20). Flame surface density models then assume that the average chemical source term is equal to the reaction rate per flame surface, multiplied by the flame surface density :

ωc = ΩΣ (2.62)
A common assumption, is to assimilate the reaction rate per flame surface to the unburnt gases density times the laminar flame velocity towards the unburnt gases. Therefore :

ωc = ρ u S u Σ (2.63)
This definition has the advantage of clearly separating the chemical and geometric contributions in the reaction source term. To close the system an algebraic expression or transport equation for the flame surface density must be provided.

For our UVCE application case, a flame surface density model seems preferable for the reasons mentioned above. It allows a wide range of applications because it does not require direct resolution of the flame front. Moreover, it permits to take into account the chemistry and the inhomogeneities of richness directly via the laminar flame velocity. These models correspond to our field of application in turbulent flame interaction regimes: flamelet models.

Theoretical basis for the numerical simulation of turbulent combustion Ξ model

One can go further with this kind of analysis by defining subgrid flame wrinkling Ξ as the ratio of flame real surface density to the amount of flame surface density captured by the mesh, assimilated to the gradient of the progress variables. One can therefore write the progress variables source term as follows : In order to close the system of equations, an algebraic expression or transport equation for flame wrinkling must then be provided. An equation has been proposed by Weller [START_REF] Weller | The development of a new flame area combustion model using conditional averaging[END_REF] :

ωc = ρ u S u ∇ (c) Ξ (2.
∂Ξ ∂t + U s • ∇Ξ = GΞ -R(Ξ -1) + (σ s -σ t )Ξ (2.65a) G = R Ξ eq -1 Ξ eq (2.65b) R = 0.28 τ η Ξ * eq Ξ * eq -1 (2.65c) Ξ eq = 1 + 2(1 -b)(Ξ * eq -1)
(2.65d)

Ξ * eq = 1 + Ξ coeff � u � S u R η (2.65e)
Note that in Weller's model, a change of variable has been done. Instead of characterizing the advancement of the combustion reaction with a progress variable c (eq. 2.55), a "regress" variable b = 1 -c is used. In that framework, the chemical source term becomes :

ωb = -ρ u S u ∇ (b) Ξ (2.66)
Cléante Langrée University of Rouen -CORIA February 8, 2022

Modeling 47

Eq. 2.65a can be simplified by giving Ξ an algebraic expression, where Ξ = Ξ eq for example :

Ξ = 1 + 2(1 -b)Ξ coeff � u � S u R η (2.67)
with R η the Kolmogorov Reynolds number, u � the sub-grid turbulence intensity and Ξ coeff a model constant equal, in our case, to 0.62. Now looking at a conservation equation across the flame front :

ρ u S u = � +∞ -∞ ωY F dx = � 1 0 ωb �∇b� db (2.68)
In addition, when applied to a turbulent flame, it can be stated that :

ρ u S t = � 1 0 ρ u S u Ξ�∇b� �∇b� db (2.69) S t S u = � 1 0 Ξdb (2.70)
Now, if Ξ is replaced in 2.70 with the expression in 2.67 :

S t S u = � 1 0 1 + 2(1 -b)Ξ coeff � u � S u R η db (2.71)
that is to say :

S t S u = 1 + Ξ coeff � u � S u R η (2.72)
Consequently, with a value of Ξ coeff = 0.6, the turbulent flame speed correlation from Gulder [START_REF] Gülder | Turbulent premixed flame propagation models for different combustion regimes[END_REF] can be obtained. Logically, other correlations can also be used, such as the Zimont correlation [START_REF] Zimont | Gas premixed combustion at high turbulence. turbulent flame closure combustion model[END_REF][START_REF] Zimont | Modelling turbulent premixed combustion in the intermediate steady propagation regime[END_REF]. In any case, a generic formula for turbulent flame speed velocity can be written as :

S T S u = 1 + C( u � S u ) Cu ( l t δ f ) C l ( p p r ef ) Cp R Cη η (2.73) With C, C u , C l , C p , C η different correlation constants.
Different correlations exist for this general formula with coefficient adaptations and will be detailed later.

Theoretical basis for the numerical simulation of turbulent combustion

These models, based on subgrid flame wrinkling, allow for the natural recovery of the laminar flame velocity in the absence of turbulence since the flame surface density of an undisturbed flame is equal to one. Moreover, these models allow a simplified consideration of mixture inhomogeneities. Indeed, the separation of the chemical and geometric contributions in the chemical source term allows to vary the local laminar flame velocity and thus to take into account the inhomogeneities of equivalence ratio in the flame. In order to fully take into account the mixture inhomogeneities, a mixing term and a corresponding transport equation are also introduced. In order to link local equivalence ratio and laminar flame velocity, a laminar flame velocity correlation is used. More details about this mixture variable will be given in chapter 3.

Turbulent flame speed model Similar to the interest in determining the consumption rate of reactants in laminar flow, describing and determining the consumption rate of a flame in turbulent flow is of both fundamental and practical importance. Nevertheless, this problem here is significantly more complex and less well defined than the problem of determining the velocity of a laminar flame. While the consumption velocity of a laminar flame is solely determined by the diffusive and chemical properties of the flame, the velocity of a turbulent flame is also determined by the features of the turbulence of the flow as well as its possible interaction with the internal structure of the flame, as previously seen.

The flame front in the flamelet regime can be locally associated with a premixed laminar flame that is simply stretched and distorted by turbulence. Turbulence main effect on combustion is the wrinkling of the flame front by large turbulent structures, resulting in an increase in the effective flame area A T . As a result, the rate of reactant consumption increases, accelerating the propagation speed of the flame front. The flame front is assumed to propagate locally at the laminar flame speed S L in the flamelet regime. The average turbulent flame front then moves at a turbulent velocity S T equal to the laminar flame velocity multiplied by the ratio of the total pleated area A T to its equivalent average non-pleated area A.

Several expressions have been proposed in the past to estimate turbulent flame velocity, resulting in various models which are typically dependent on several parameters of the turbulent flow (u", l t , etc.) and the flame (S 0 L , δ 0 L , Le, etc.). The problem is still open today and is of great interest in the context of numerical simulations of turbulent reactive flows because many models can be used to simulate the effect of turbulence on the flame [START_REF] Flohr | A turbulent flame speed closure model for LES of industrial burner flows[END_REF] :

S t S u = 1 + f (Re, Da, P r ) (2.74)
with, the following definition for the Reynolds, Damköhler and Prandtl numbers

Re = u � l t ν , Da = τ t τ f , P r = ν D (2.75)
It exists dozen of turbulent flame speed correlations [91, Tab. 1 p. 744]. They are based on the same flamelet hypothesis, but have some intrinsic specificities based on the primary configuration they were aiming at modeling. Available CFD codes, both in industry and academical field, use different correlations for large scale deflagrations. A comparative study of these correlations and their response to our points of interest is presented in chapter 4. However, the main used relations, that will be tested in the present work, are given :

Gulders [91]
The Gülders correlation used in OpenFoam solver XiFoam :

S t S u = 1 + C G � u � S u � 1/2 Re 1/4 t with C G = 0.62 (2.76) Zimont [92]
Zimont correlation is the turbulent flame speed correlation used in ANSYS-Fluent.

S t S u = 1 + C Z u � Da 1/4 with C Z = 0.52 (2.77)
FLACS [START_REF] Gexcon | Flacs-cfd v21.1 user's manual[END_REF] The CFD code FLACS use a specific turbulent flame speed correlation [START_REF] Gexcon | Flacs-cfd v21.1 user's manual[END_REF].

S t S u = 1 + C F � u � S u � 0.

Conclusion

To sum it all up, let us recap our choices of modeling for our application case. As said previously, a geometrical approach is chosen for the turbulent combustion model. In particular the progress variable equation will be modeled via a wrinkling factor equation. This model have been chose for several reasons : First this is the modeling choice more able to handle low resolution of meshes. As a matter of fact, resolving totally the flame front requires to have several mesh sizes in the flame front thickness. It means for example for laminar flame fronts, a mesh size of the order of 10 -4 m, which is unrealistic in our case. Some modeling choices attempt to overcome this issue without giving up the total flame front resolution, the Thickened Flame Model for example. But this model seemed also unrealistic in the sense that it would require a thickening factor so high that the results could not be considered with much confidence.

Theoretical basis for the numerical simulation of turbulent combustion

Second, for our cases, a discrimination between 3 phases is needed : the fresh flammable mixture, the burnt gases and the inert ambient air. Because meshing the domain as far as possible is needed (first for assessing the overpressure effect in far field and second, for numerical boundary condition issues), the model has to take into account the finite volume of the initial flammable cloud and not propagate the flame front further, where there is no flammable mixture to burn.

Finally, the model has to be consistent throughout all steps of propagation of an UVCE : from the initial spherical laminar propagation to the late turbulent fully deformed flame front. For all those reasons, the geometrical approach was a suitable option. First because it allows a discrimination between the chemical and the "geometrical deformation" contribution to the flame speed. This allow to easily introduce the inhomogeneities (which will only impact the chemical contribution) and to treat the initial stages of the propagation (i.e. when the geometrical contribution is negligible). Though the latter assumption is highly based on the infinitely thin flame front analysis, it is reasonable in the application range of UVCE with initial flammable mixture at rest. Therefore, this is the equation system we're are solving. In the next chapter the numerical methods and how the theoretical set of equations are defined in the numerical solver used will be presented.

Chap. 3 | Numerical tools for UVCE simulation 3.1 Introduction

The partial differential equations discussed previously must be transformed into a set of algebraic equations to be resolved. Numerical solution techniques are divided into three categories : finite difference, finite element and finite volume. Each has the goal of converting differential equations to algebraic equations. Distinctions between these three methods rely on the approach use to transform differential equations into algebraic equations.

(1) Finite Difference.

In the finite difference method the representation of any function is based on the knowledge of point values. Then, partial derivatives are substituted by a series expansion representation [START_REF] Ferziger | Computational methods for fluid dynamics[END_REF], commonly Taylor series, which is truncated typically at first or second order even if some method used higher order approximation. It is quite easy to apply this strategy to simple geometry with regular shape. However, it increases rapidly the number of term to be considered in the Taylor series for more complex geometries with irregular mesh and with the dimension of the problem. Thus, the finite difference approach is fairly efficient and effective on structured grids. The downside is that, unless specific precautions are taken, conservation is not enforced. Furthermore, with complex flows, the constraint to simple geometries is a considerable disadvantage. (2) Finite Element.

When using the finite element method, the governing partial differential equations are integrated over an element that correspond to a piece of volume completed with an approximation of the integration at a given order that depends of the number of point used to characterize the element. An important advantage of finite element methods is the ability to deal with arbitrary geometries at a specified order. However, depending of the integration approximation the algebraic equations can become quite complex leading to a difficult physical interpretation on each term.In addition, it is quite difficult to ensure a mathematical conservation for key quantities [START_REF] Peter | Aspects of conservation in finite element flow computations[END_REF]. (3) Finite Volume.

The finite volume method (FVM) describes any function by the value of its integration over predetermine volume. It is appropriate for complex geometries after a meshing process that aims to split the total volume in small pieces of volume (cells) that are interconnected trough small faces. By concept, the approach is conservative as long as

Cléante Langrée UVCE Simulation February 8, 2022 52
Numerical tools for UVCE simulation the surface fluxes (which represent convective and diffusive fluxes) are subtracted and added to the neighbors cells [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. Since PDEs used for CFD are conservation laws, the FVM is an appropriate candidate for tackling CFD challenges. Still it needs to evaluate fluxes at cell face boundaries which can be quite complex for high order method.

Let us summarize the main goals of the numerical approach :

• Solving the Navier-Stokes equations completed by the energy equation and the mass fraction equation of considered species • Unsteady, three-dimensional resolution on complex geometries • Need to fulfill the conservation of mass, momentum and energy • Complex mesh with cell volumes that can have complicated shapes • Computation of large size with respect to the smallest length scales : flame thickness, Kolmogorov length scale ...

In this context, the set of equations presented above will be solved by applying the finite volume method on complex meshes. The OpenFoam solver has been selected for several reasons :

(1) OpenFoam is competitive technologically with all other commercial solutions currently available, the quality of the results is similar. (2) Open source philosophy with full access to source codes. Moreover, It is easy to set up collaborations because of the availability of the solver. (3) Easily scriptable, which means that many processes can be completely automated.

Changes in geometry or flow conditions can be quickly checked and compared to previous results once a case has been set up for the first time.

This chapter summarizes the different techniques of discretization of the geometry (mesh) and the equations (FVM) which are implemented to solve the propagation of a premixed and partially premixed flame in a complex environment.

After a brief presentation of the general numerical aspects, the xiFoam solver is detailed so that the reader can fully understand the whole work flow from the equations established in the previous chapter to the numerical results ontained and presented in the next chapter.

OpenFoam : A brief description

In last decades, numerical computations in fluid mechanics, using computational tools such as OpenFoam, have developed considerably, thanks to increasingly powerful computing resources that allow the numerical solution of partial differential equations with great accuracy. However, despite the great progress, a major challenge remains : the simulation of complex real problems. For most of these problems, there is no analytical solution and the validation and verification of the results becomes a difficult but necessary task. OpenFoam is a simulation tool mainly focused on the numerical solution of fluid mechanics equations based on finite volume library. Distributed under the GNU/GPL open source license, OpenFoam is widely used in both academia and industry.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

Finite Volume Method 53

Written on the C++ language, OpenFoam is an object-oriented software with all the properties of this programming paradigm such as encapsulation, polymorphism and inheritance [START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF]. In addition, OpenFoam contains a very large number of libraries for solving numerical fluid mechanics problems as well as a variety of specialized data structures for handling different fields (e.g., scalar, vector, or tensor fields) for different types of physical properties.

Different solvers developed on OpenFoam cover a wide range of applications with several solvers dedicated to compressible and incompressible flows.

More specifically, many solvers are also available for other types of physical phenomena such as combustion, heat transfer, multiphase flows, particle tracking flows, electromagnetism, fluid-structure interactions, etc.

OpenFoam allows combining its wide range of solvers with another wide range of boundary conditions, which include standard types of boundary conditions such as Dirichlet, Neumann, symmetric and periodic types.

Eventually, after discretizing the set of equations, the solution of the problem consists in optimizing and solving a linear system based on matrix representation. Such systems in OpenFoam uses iterative methods which are more efficient in term of computing resources that direct resolution for large system. The solvers use preconditioning, so that the convergence of the preconditioned system is much faster than for the original, and use smoothing, which are transformations to reduce the iteration-mesh dependence. Some options are : ICC, conjugate gradient based solver with incomplete Cholesky preconditioning; BICCG, solver with diagonal preconditioning and incomplete LU decomposition; PCG, bi-conjugate gradient preconditioning solver for symmetric LDU matrices that uses runtime selectable preconditioning, etc.

One of the objectives of this work is to determine and, if necessary, adapt within OpenFoam, the best combination of all these parameters (solver, boundary conditions, resolution methods, etc.) to capture the physics of the propagation of flames throughout obstacles.

Finite Volume Method

The principle of discretization is to transform the partial differential equations considered in the previous chapter into a system of algebraic equations. It is this new system of equations which, once derived within the geometry, gives the spatio-temporal variations of the quantities of the system, under certain conditions which will be developed below.

The discretization process can be divided into two steps, the first is the decomposition of the domain into a set of elementary cell volumes, the so-called control volume (CV). The second is the integration of the problem equations over all control volumes. This step implies that the equations of the problem have been previously transformed into a system of linear equations which require the CV to be small enough to allow for linearization process.

Numerical tools for UVCE simulation

Domain discretization

Each control volume is constructed around a point P at its center of gravity as shown in figure 3.1. Let N be the center of the neighboring control volume. d is defined as the vector linking P to N . Let S f be the unit vector orthogonal to the common face S f of the two control volumes. The variables such as velocity and pressure are defined as the mean value over the CV and associated to the center P of the CV control volume.

This distribution simplifies the implementation in the code and minimizes the amount of information needed related to the CV geometry. Some variables need to be defined at the surface, for example for flux calculations; these will be denoted with the index f .

Q P = Q = 1 V P � V P Q(x)dV (3.1)
All the control volumes in the domain form the mesh on which the equations are discretized.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

Finite Volume Method 55

From a practical point of view, the mesh containing all the CVs is generated from a specific software. This subject alone is the subject of extensive research because the quality of the mesh determines the quality of the results. Indeed, without a high-quality mesh to capture the solution's complexities, results may be unreliable or prohibitively expensive. Furthermore, suitable meshes for huge problems are frequently too large for traditional workstations, necessitating mesh development and solution on large-scale distributed memory systems. Several point have been considered :

(1) Having a well-defined, simple, clean, and watertight geometry is really the difference between a high-quality mesh and a bad mesh. If possible, geometries should be with no hole and free of irregularities like intersections and sharp outcroppings. With this consideration in mind, as will be shown later at the end of this chapter, the final complex simulation of this work was therefore performed a first time with a simplified geometry and a second time with a more complex geometry presenting angles difficult to be caught by the solver.

(2) Maintaining low value the cell's skewness ratio is crucial for accuracy and quality. Preserving the skewness ratio of each cell in complex geometries might be challenging, if not impossible. Different scenarios necessitate and dictate different skewness ratios, but in this work, strong cell distortions have been avoided as much as possible by using adapted meshing tools according to the configuration.

In the next chapter, the SALOME-mesh tool was used to generate an unstructured mesh (tetra) with very few distortion on a spherical geometry. In the following chapter a combination of the blockMesh and snappyHexMesh tools was used to generate a hexa dominant mesh in the domain with strong rectilinear singularities but with tetra cells for the complex areas near obstacles.

A brief presentation of the meshing tools that have been used in this work is now proposed. blockMesh For simple geometries, the mesh generation utility blockMesh can be used. The blockMesh utility generates graded and curved edges to generate parametric meshes based on blocks. A dictionary file entitled blockMeshDict is contained in the constant/polyMesh directory of a case and is used to construct the mesh. blockMesh reads this dictionary, constructs the mesh, and writes the mesh data to the same directory's points and faces, cells, and boundary files. The blockMesh algorithm works by breaking down the domain geometry into a set of one or more three-dimensional hexahedral blocks. Straight lines, arcs, and splines can be used as block edges. The mesh is ostensibly specified as a number of cells in each block direction, which provide enough information for blockMesh to build mesh data. This meshing tool produces high-quality meshes and is ideal for very simple geometries. The effort and time necessary to build up the dictionary increases dramatically as the geometry becomes more complex.

In general using blockMesh is only the first step for addressing a complex case before of a more sophisticated mesh construction that combines the tools blockMesh + snappyHexMesh.

Cléante Langrée

UVCE Simulation February 8, 2022

Numerical tools for UVCE simulation

Because the background mesh for snappyHexMesh is usually a single rectangular block, blockMesh is straightforward to use in this framework. An example of mesh obtained with blockMesh is shown in figure 3.2. Figure 3.2 -Mesh generated by blockMesh. snappyHexMesh For complex geometries, the mesh generation utility snappyHexMesh can be used. It generates 3D meshes containing hexahedra and split-hexahedra from a triangulated surface geometry in Stereolithography (STL) format. The mesh is generated from a dictionary file named snappyHexMesh located in the system directory and a triangulated surface geometry file located in the directory constant/triSurface. snappyHexMesh has two main fonction : refine certain region of an existing mesh, and insert a structure inside an existing mesh (i.e. adapt the mesh around it). The following step are necessary to create a mesh with snappyHexMesh, they may be seen in figure 3.3 :

(1) Creation of a background mesh with blockMesh .

(2) Definition of geometry.

(3) Creation of a castellated mesh around the body / obstacle. (4) Creation of a snapped mesh or a mesh that fits the body.

(5) Meshing of boundary layers or the addition of layers close to the surfaces.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

SALOME

The SALOME software is the result of a collaboration between the CEA, EDF and Open Cascade. This software was designed to be a 3D CAD tool for the creation of closed objects, i.e. objects that facilitate mesh creation for the study of internal or external flows, unlike software developed for mechanical purposes that limited from the meshing point of view, such as Solidworks or CATIA. Moreover, SALOME is opensource, thus our whole work flow : pre-processing, simulation and post-processing is based on open source software. SALOME offers many features. However, in the context of this work, two major aspects have been used :

(1) create, modify, repair and clean CAD models (geometry module) (2) generate mesh from CAD models; edit meshes; check mesh quality; import and export mesh data to an appropriate format to be used by OpenFoam (mesh module).

Figure 3.4 shows an example of CAD (obstacles) and mesh generated with SALOME. For the sake of clarity, the complete 3D mesh is not shown. Only the cells triggered by the passage of the flame are shown. The objective is to have the best possible mesh with very little distortion and elongation of the meshes. In the case shown in the figure, it is clear that the software has generated a homogeneous mesh of good quality (from a geometric point of view). This configuration will be detailed in the next chapter.

Equations discretization

A transport equation of any quantity Q (i.e. velocity, mass, energy ...) is typically of the form :

∂ρQ ∂t + ∇ • (ρUQ) -∇ • ρΓ Q ∇ (Q) = S Q (Q) (3.2)
It can be integrated over any control volume V P :

� V P ∂ρQ ∂t dV + � V P ∇ • (ρUQ) dV - � V P ∇ • ρΓ Q ∇ (Q) dV = � V P S Q (Q)dV (3.3)
Then, the Gauss theorem, or divergence theorem, is used to rewrite the equation. The Gauss theorem can be expressed, for any vector a as :

� V P ∇ • (a) dV = � δV a.ndS (3.4)
with δV the external closed surface of the volume V P , and n the outward pointing normal vector of the considered surface element. Therefore, eq. 3.3 can be rewritten as :

� V P ∂ρQ ∂t dV + � δV ρQU.ndS - � δV ρΓ Q ∇ (Q) .ndS = � V P S Q (Q)dV (3.5)
Cléante Langrée University of Rouen -CORIA February 8, 2022

3.3 Finite Volume Method 59

Surface and Volume integrals approximations

In equation 3.5, there are volume and surface integrals that need to be formulated based control volume in order to be able to solve the equation. In OpenFoam CV are general polyhedra thus surface integral are split each surface of the polyhedra :

� δV ρQU.ndS = � f nS f .(ρuQ) f (3.6) � δV ρΓ Q ∇ (Q) .ndS = � f nS f .(ρΓ Q ∇ (Q)) f (3.7)
Notice that mean values at surface, such as (Q f ), require an approximation based on the known mean values of this quantity over the control volume. A linear combination based on both neighbor CVs is generally used :

(Q) f = f x (Q) P + (1 -f x )(Q) N (3.8)
In addition source term integral is split in two parts to make apparent the linear dependency on the transported variable Q :

� V P S Q (Q)dV = S c V P + S p V P Q P (3.9)
With S c , S p corresponding respectively to the further explicit and implicit part of the source term. With those approximation our discretized equation can be written :

� V P ∂ρQ ∂t dV + � f nS f .(ρuQ) f - � f nS f .(ρΓ Q ∇ (Q)) f = S c V P + S p V P Q P (3.10)
Once our equation has been discretized, surface interpolation and temporal derivative terms still need to be addressed. Derivations of the equation and numerical schems are fully described in Jasak 1996 [START_REF] Jasak | Error analysis and estimation for the finite volume method with applications to fluid flows[END_REF].

Temporal scheme

For the temporal derivative term ∂Q ∂t in eq. 3.10, several expressions, that are already available in OpenFoam can be used.

Numerical tools for UVCE simulation Euler

The Euler method is the simplest basic time integration method. It simply approximate the time derivative based on the previous time step value Q 0 , and the current time step value Q itself :

∂Q ∂t = Q -Q 0 ∆t (3.11)
It's a first order, bounded, implicit method.

backward Following the idea of the Euler method, the backward scheme extend it to two previous time step (Q 0 and Q 00 , making it a second order "version" of the Euler method.

∂Q ∂t = 1 ∆t ( 3 2 Q -2Q 0 + 1 2 Q 00 ) (3.12)

Crank-Nicolson

The Crank-Nicholson method is a blended explicit-implicit approach second order scheme. It is based on the Euler approach, but the rest of the equation is evaluated explicitly for one half and implicitly for the other half. Note that this approach is coded in OpenFoam with a blending coefficient ψ, to mix this approach with the implicit Euler method. ψ = 1 corresponds to pure Crank Nicholson scheme whereas ψ = 0 corresponds to pure Euler. A value of ψ = 0.9 is advised for simulation where pure Crank Nicholson can become too unstable.

Spatial schemes

As for the convective term in eq. 3.10, a formulation of the value Q f , i.e. the value of the quantity Q, at the face f, is needed.

Linear scheme

The linear scheme interpolates linearly the face value from the 2 neighboring cells, as represented in fig 3.5.

Q f = f x Q P + (1 -f x )Q N (3.13)
Cléante Langrée University of Rouen -CORIA February 8, 2022 with

f x = f N P N (3.14)
It is a second order, unbounded scheme.

Upwind difference scheme

The face value can also be expressed as the value of the neighboring cells from which the general flux come from. For example, in fig 3.5, if the general flux was going from cell P to cell N, then φ f would be equal to φ P . In practice it can be written as :

Q f = � Q P if F ≥ 0 Q N if F < 0 (3.15)
It is a first order, bounded scheme. It can be extended to second order scheme with a weighted upwind interpolation for the linearUpwind variant.

Boundary conditions

Finally, to solve the set of equations, boundary conditions needs to be set. As seen previously, each solution for each cell is highly dependent on the solution of the adjacent cells, notably because of the flux between those two cells. But cells that have a face in contact with the border of the numerical domain require a particular treatment. Therefore, those boundary fluxes must either be known or be calculated via values inside the calculation domain completed by other information that aim to represent the boundary conditions. Two main types of boundary conditions exists and are overly used in industrial applications :

Numerical tools for UVCE simulation

• Dirichlet conditions : Fixes the value of the considered field at the boundary • Neumann conditions : Fixes the value of the gradient of the considered field Other conditions like symmetrical conditions are also commonly used. In our case, the boundaries of our simulation domain will be of three types : the ground, the industrial obstacles (mainly metal) and the open atmosphere. Note that the free, open atmosphere is a challenging boundary condition to set up, as it is expected to not reflect any wave propagating outward of the computational domain (as it does not represent any change of condition of propagation).

Solid conditions : ground and obstacles

Near wall boundary conditions are a vast topic of innovation. Several models have been developed to capture the dynamic of a flow near a wall, for each quantity characteristic of this flow (e.g. Temperature, velocity ...). An important tool for studying the interaction between the main turbulent flow and walls, is the determination of y+. It is a quantity that can be calculated based on the velocity and viscosity of the flow, to represent how refined the mesh is close to the wall. Therefore it represents how relevant/needed is the use of wall functions.

OpenFoam has a lot of wall functions coded, some used in our simulation, even if most times, the study of the y + make it not clearly relevant.

Atmospheric conditions

The atmospheric boundary conditions is a challenging part in open simulation. The ideal for this kind of boundary condition would be to be fully non reflective, especially for the pressure field. Such fully non reflective boundary condition exists but is often complex to implement, and for our simulations more classical approaches will be used.

A usual trick consists in "pushing" the boundary as far as possible to avoid the reflected waves to impact the phenomenon understudy. But as pressure waves propagates at the speed of sound (≈ 340m/s). Thus the distance to avoid these reflected waves increases as the total simulation time increases. Therefore increases unreasonably the computational cost of the simulation.

OpenFoam has some "pseudo non reflective" boundary condition coded, or non specific conditions that can be used. For both pressure and velocity, waveTransmissive will be used, as well as totalPressure for pressure in combination to fluxCorrectedVelocity for velocity.

Cléante Langrée

University of Rouen -CORIA February 8, 2022 The waveTransmissive is coded in OpenFoam. It is based on the resolution of a simple conservation equation of the considered variable. This conservation law is applied to any field Q, with a wave velocity, which is calculated to be the flux velocity plus the speed of sound :

w P = φ P |S f | + � γ ψ P (3.16)
With φ P being the patch face flux, S f the patch face area vector, γ the ratio of specific heats and ψ P the patch compressibility. It gives us the following equation to solve for each field Q, in our case, pressure :

DDt(W , Q) = ∂Q ∂t + W • ∇ (Q) = 0 (3.17) totalPressure + fluxCorrectedVelocity
Another way of setting up pseudo non reflective boundary condition is a set of boundary conditions were both pressure and velocity are fixed to values calculated by : p P = p 0 -0.5ρ|U| 2 (3.18)

U P = U c -n(n.U c ) + nφ P |S f | (3.19)
With U p being the velocity at the patch, U c the velocity in cells adjacent to the patch and n the patch normal vector.

Pressure-velocity coupling

Solving equation 2.20 and 2.19 simultaneously is not an easy task, both in compressible case and incompressible cases. The computational cost is high, and for an application like ours a resolution algorithm that can find an approximate solution by solving the transport equations one by one (the so-called segregated approach) will be used.

In order to do so an equation for pressure has to be derived.

Momentum equation (2.20) can be expressed in a semi discretized form as :

M.U = -∇p (3.20)
Numerical tools for UVCE simulation with M the [n × n] matrix of coefficient for the n cells. For each direction i related to

x i :     M 11 M 12 ... M 1n M 21 M 22 ... M 2n ... ... ... ... M n1 M n2 ... M nn     .     (U x i ) 1 (U x i ) 2 ... (U x i ) n     =      ( ∂p ∂x i ) 1 ( ∂p ∂x i ) 2 ... ( ∂p ∂x i ) n      (3.21)
In-cells values can be separated from neighboring values :

M.U = A.U -H(U) (3.22)
With A the matrix of coefficient from in-cell value and H(U) the matrix of coefficient from neighboring cells and source terms. A being diagonal matrix by definition, it can easily be inverted. Therefore :

U = A -1 H(U) -A -1 ∇p (3.23)
With H the operator defined by definition :

H(U) = A.U -M.U (3.24)
Then ψ the compressibility is defined :

∂ρ ∂t = ψ ∂p ∂t (3.25)
And this expression can be injected in the continuity equation (2.19) :

ψ ∂p ∂t + ∇.(ρ(A -1 H(U) -A -1 ∇p)) = 0 (3.26)
Which give a Poisson equation for the pressure :

∇.(ρA -1 ∇p) = ψ ∂p ∂t + ∇.(ρ(A -1 H(U))) (3.27)
Two equations are therefore obtained, one for velocity and one for pressure. It exists several algorithm to solve those two equations sequentially until a convergence is reached.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

3.4 Pressure-velocity coupling 65

Semi Implicit Method for Pressure Linked Equations -SIMPLE

The SIMPLE algorithm for resolving the pressure velocity coupling goes as follows :

• Solve the momentum equation (2.20). The pressure is taken from the previous time step. Note that the velocity field solution at this step does not satisfy the continuity equation (2.19). This step is often called "momentum predictor step" • Solve the Poisson equation (3.27) using the predicted velocity field. This loop is often called SIMPLE loop, or outer loop.

Pressure Implicit with Splitting Operator -PISO

The PISO algorithm for resolving the pressure velocity coupling goes as follows :

• Solve the momentum equation (2.20). The pressure is taken from the previous time step. Note that the velocity field solution at this step does not satisfy the continuity equation (2.19). This step is often called "momentum predictor step" This loop is often called PISO loop, or inner loop.

PIMPLE : PIso + sIMPLE

Those two way of solving the pressure velocity coupling can be combined, by having both an outer and an inner loop. The inner loop aims to ensure pressure convergence thus to determine the pressure that fulfill the pressure-velocity coupling. Whereas the outer loop is more dedicated to velocity convergence solving the non linear term of the momentum equation.

Note that, by definition, running a PIMPLE solver with one inner loop is the same as running a SIMPLE algorithm, and running a PIMPLE solver with one outer loop is the same as running a PISO algorithm.

Numerical tools for UVCE simulation

XiFoam solver

The great advantage of OpenFoam is the possibility to retrieve the code of transport equation that are solved. The code is based on OpenFoam library that represent tensor calculus. In the following the algorithm and the code are explained together with the corresponding transport equations.

General scheme

XiFoam is a combustion solver based on the b-Xi model [START_REF] Weller | Application of a flame-wrinkling les combustion model to a turbulent mixing layer[END_REF]. In XiFoam, Navier Stokes equations are solved (eq. 2.40), (eq. 2.41), (eq. 2.43) in addition to a "b" equations (eq.

Momentum equation

The Favre averaged momentum equation (eq. 2.41) is solved.

∂ρ � U ∂t + ∇ • � ρ � U � U � = -∇ • � ρ � U �� U �� � -∇p + ∇ • τ + F (3.29)
The body force source term is neglected.

∂ρ � U ∂t + ∇ • � ρ � U � U � + ∇ • � ρ � U �� U �� -τ � = -∇ (p) (3.30)
The viscous stress tensor τ is defined as :

τ = µ∇ � U + µ∇ � U T -µ 2 3 (∇ • � U)δ ij (3.31)
And representing the turbulent Reynolds stress with the hypothesis of turbulent viscosity µ t and µ ef f = µ + µ t :

-τ = ∇ • � µ t ∇ � U � -∇ • � µ ef f ∇ � U � -∇ • � µ((∇ � U) T - 2 3 tr (∇ � U)I) � (3.32)
Moreover the Reynolds stress tensor is defined as R = � U �� U �� . In XiFoam, ρRτ , is implemented via divDevRhoReff in ReynoldsStress.C (lst. The MRF.DDt(rho, U) line is from the OpenFoam module Moving Reference Frame, for rotating engine. As our reference frame is steady, this term is not activated. Also the kinetic energy K is calculated. It will be used in the transport of the energy variable.

Progress variable equation

The turbulent flux of c is modeled as :

ρ � U �� c �� = α t ∇ (� c) (3.33)
From (eq. 2.58) with (eq. 3.33) and α ef f = α + α t is obtained :

∂ρ� c ∂t + ∇ • � ρ � U � c � = ∇ • � α∇ (� c) -ρ � U �� c �� � + ωc (3.34) ∂ρ� c ∂t + ∇ • � ρ � U � c � = ∇ • (α ef f ∇ (� c)) + ωc (3.35)
A change of variable is then applied :

b = 1 -c (3.36)
Cléante Langrée University of Rouen -CORIA February 8, 2022 Therefore :

∂ρ � b ∂t + ∇ • � ρ � U � b � = ∇ • � α ef f ∇ � � b �� + ωb (3.38)
ωb is modeled as :

ωb = -ρ u S u Ξ | ∇ � b | (3.39)
The flame flux term is defined as : The stCorr is a correlation term used during the igniting time of the simulation when using a special module coded in XiFoam. As this method of ignition is not used in our simulations, it will not be discussed here. Outside of the ignition process, stCorr = 1.

Φ st = ρ u S u Ξn f (3.
n f = ∇ � b | ∇ � b | .S f (
| ∇ � b | 2 = ∇ � b • ∇ � b = ∇ � b • ∇ � b + � b∇ • � ∇ � b � -� b∇ • � ∇ � b � (3.42)
With, for any scalar f , and any vector A :

∇ • (f A) = ∇f • A + f ∇ • (A) (3.43)
Therefore :

| ∇ � b | 2 = ∇ • � � b∇ � b � -� b∇ • � ∇ � b � (3.44) | ∇ � b |= ∇ • � � b ∇ � b | ∇ � b | � -� b∇ • � ∇ � b | ∇ � b | � (3.45) ρ u S u Ξ | ∇ � b |= ∇ • � � bρ u S u Ξ ∇ � b | ∇ � b | � -� b∇ • � ρ u S u Ξ ∇ � b | ∇ � b | � (3.46) ωb = ∇ • � Φ st � b � -� b∇ • (Φ st ) (3.47)
And finally : In OpenFoam, the laminar flame speed is retrieved via correlation. The Gulders correlations are used :

∂ρ � b ∂t + ∇ • � ρ � U � b � + ∇ • � Φ st � b � -� b∇ • (Φ st ) -∇ • � α ef f ∇ � � b �� = 0 (3.
Su 0 (T, p, φ) = W × φ η × e -ξ(φ-σ) × ( T T 0 ) α × ( p p 0 ) β (3.49)
Where W , η, σ, α, β are the Weller model constants and φ, the equivalence ratio, being recovered from f t and ( Y air Y f uel ) sto , the stoichiometric air-fuel ratio.

Listing 3.9 -Flame speed correlation -equivalence ratio -Gulder.C In addition to the progress variable, a "total fuel mass fraction" f t can be transported. Note that f t represent both unburnt and burnt fuel mass fraction. "Burnt fuel" may seems odd at first, but what it really represent is literally the part of the fuel that has already burned, that is to say, the part of the combustion product mixture, that came from the fuel (and not from the oxydant).

∂ρf t ∂t + ∇.(ρ � Uf t ) = ∇.(µ∇f t ) (3.50)
With µ the dynamic viscosity. The effects of equivalence ratio inhomogeneity on the flame front propagation speed is taken into account by the S u term in the b equation with the Gulder flame speed correlation (3.49) ( [START_REF] Weller | Application of a flame-wrinkling les combustion model to a turbulent mixing layer[END_REF]). From this total fuel mass fraction, others volumes scalars are derived to caracterize the mixture. In XiFoam they're named ox (oxydant mass fraction), pr (products mass fraction), f u (unburnt fuel mass fraction) and f r es the residual fuel at a certain stoichiometric ratio. They are linked so that :

φ = ( Y air Y f uel ) sto f t 1 -f t (3.
f u = bf t + (1 -b)f r es (3.52) o x = 1 -f t -(f t -f u )( Y air Y f uel ) sto (3.53) p r = 1 -f u -o x (3.54)
Cléante Langrée University of Rouen -CORIA February 8, 2022 With f r es being the residual fuel mass fraction left in the burnt mixture.

f r es = max(f t - 1 -f t ( Y air Y f uel ) sto ; 0) (3.55)
In summary, mixture composition can be recovered from the two parameters b and f t as follow (in the case of a stoichiometric mixture) :

Energy equation

In XiFoam, once the mixture is ignited, an energy equation (eq. 2.43) is resolved for both unburnt gases (lst. 3.15) and the total mixture (lst. 3.14). Note that here, the equation solved by our solver is slightly different from the one presented in 2. Instead of considering the non chemical enthalpy, the absolute enthalpy is transported. The non chemical enthalpy presented in 2 is the sum of the sensible and kinetic energy. Here, the absolute energy we're considering is sensible and chemical energy. Therefore the variation of the kinetic energy K is introduced in the equation, while the variation of chemical energy is incorporated into the total enthalpy h t . The exact solved equation is :

∂ρ � h t ∂t + ∇ • � ρ � U � h t � + ∂ρ � K ∂t + ∇ • � ρ � U � K � - ∂p ∂t -∇ • � α ef f ∇ • � � h t �� = 0 (3.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

3.5 XiFoam solver 77

Conclusion

Turbulent combustion can be handled in many different ways using OpenFoam. This chapter focuses on the XiFoam solver, which is the application of the model selected in the previous section. But other model that won't be presented here are also usable, such as reactingFoam. reactingFoam is based on a more detailed description of the chemistry, by tracking directly the species mass fraction and detailing the reaction scheme used. For our purpose and final application using reactingFoam directly would have required to resolve all the scales of the problem that is impossible due to large scale of the computational domain with respect to the smallest length scale induce by the flame structure and the turbulent scales.

Notably, XiFoam includes a separate module for ignition. This module has not been used, but instead a "burnt volume deposit" method, which will detailed on the next section. The XiFoam ignition module proved to be rather complicate to tune and operate, and as the focus was put on the propagation of the flame front, and not that much on the dynamic of the ignition, a more straight forward approach has then been chosen.

In the next section, this solver will be evaluated in test configurations to validate that it can reproduce correctly the behavior of a turbulent large scale flames. One test case will be a simple 1D laminar test case, to ensure that laminar combustion is accurately represented. The other will be a more complex 3D large scale configuration, with the introduction of complex geometry.

Chap. 4 | Flames expanding through obstacles 4.1 Introduction

In this chapter, calculations have been performed in academic settings in order to test both the capabilities of the physical and numerical models selected in this work. First, a configuration of a 1D laminar flame has been considered in order to validate the capabilities of the numerical model to capture the properties of the methane flame, which is the object of this study. In a second step, the focus has been put on the case of a spherically expanding turbulent flame in a medium with and without obstacles. The structure of the flame within the obstacles and the capabilities of different correlation models are highlighted.

Laminar premixed flames

First, 1D laminar flames simulations were performed with the XiFoam solver. The objective is to select appropriate numerical models and to test their limitations before moving on to larger scale simulations. The main interest here is to be able to capture the main properties that the prevention and risk industry is interested in : (1) the flame speed and (2) the pressure variations. Indeed, capturing the laminar flame velocity will demonstrate a good consideration of the reactive aspects in the description of the flame based on the progress variable selected. It should also be kept in mind that the different turbulent correlation models for combustion that will used are based on a correct estimation of the laminar flame speed from which the turbulent flame speed is estimated.

Then calculations are necessary to refine the numerical parameters, such as the boundary conditions and especially the pressure conditions which can be particularly difficult to set if pressure waves are not correctly evacuated from the calculation domain as they do in the physical domain. Some physical properties inherent to the description of the flames cannot be respected. For example, the thickness of the flame. However, this is not the objective. Indeed, the flame will be grid-thickened in order to be able to work with a mesh large enough to be used in geometries of the order of several meters. "Thickening" the flame is therefore part of the Flames expanding through obstacles modeling process. But if one can estimate the average position of the flame front over time (and thus its speed) as well as the pressure effects that result from its passage, very valuable information can be derived from these calculations. This could allow the implementation of safety measures and the design of industrial sites that limit the harmful effects of a blast to a minimum. One dimensional simulations have been performed in order to evaluate the core properties of the planar propagating flame issued from the models described in chapter 2. From a practical point of view, the domain is a volume of 30 cm long, discretized by steps of length ∆x according to the x direction (see fig. 4.1). The y and z directions have the same discretization step but one cell only is considered since the calculation is uni-dimensional. The discretization step ∆x is variable according to the different convergence tests performed. Its basic value is set to ∆x = 0.15 cm but evolved from ∆x = 0.05 cm to ∆x = 1 cm.

1D Configurations

In order to evaluate XiFoam results from the thermo-chemical point of view, an accurate resolution of the corresponding flames is performed with the Cantera library. Cantera is an open-source suite of tools for problems involving chemical kinetics, thermodynamics and transport processes. It allows to analyze free flame propagation in any fuel-air mixture 1D configuration. The flow is solved using mass conservation, energy, state equations and species diffusion equations. The solver uses a Newton-Raphson integration to solve coupled differential equations. The solver first solves the domain based on a basic isotropic grid and then it refines it according to slope and curve parameters of the solution variables. The refinement is done till the tolerance criteria is met. The flame velocity is assumed zero and other velocities are calculated relatively to this speed. This transformation makes the transient phenomenon, a steady state one. Then, the results are transformed back to original frame of reference with respect to the ground. Note that in our simulation performed with XiFoam, a different reference frame is used since in this case, it is the fresh gases velocity that is set to 0 and not the flame velocity. A change in the reference frame of the results is therefore carried out before they are compared.

Cléante Langrée

University of Rouen -CORIA February 8, 2022 Flames expanding through obstacles

Ignition setup

Contrary to some specific premixed configuration such as spark ignition engines, ignition and initial flame kernel development is not a vital process in this work. Indeed, the main objective is to recover the propagation of the flame over long distances and therefore, the initial energy deposit is less important than in the case of engine studies. However, the initialization phase of the ignition must be done correctly because it can generate spurious pressure waves that will deeply modify the first moments of the flame. Moreover they will be particularly difficult to evacuate if their amplitude is too high.

OpenFoam proposes an ignition library : ignite, based on the deposit of a fixed amount of energy over a given time and in a prescribed volume. In this particular configuration, this means varying the regress variable b from 1 to 0 in a given volume and time. However, the tests that have been carried out with this methodology did not satisfy us. Indeed, non physical results appear at first (see fig. 4.2-(a,b)) : for example an adiabatic temperature superior to 3500 K. Even if these non-physical data are very quickly eliminated after the first moments of the simulation, this is not very satisfactory. Moreover, this generates pressure waves of such a level that their evacuation is difficult. This has an impact on the flame instabilities that appear. It is clear that this ignition module needs improvement and it has been decided not to use it.

The methodology used, very classical in the field of premixed flames, consists in initializing directly an established flame profile at time 0. The ignition, which does not interest us specifically, is therefore not taken into account and the simulation starts from a core of burnt gas has been deposited in a small area of 0.5 cm at the left (x = 0) of our computational domain. A seen in fig. 4.2-(c,d), temperature and velocity profiles are more coherent that the ones obtained with ignite, even if in this case a large step length is used.

It will be shown in the next section that the flame has a behavior and a propagation speed that corresponds well to what is expected from the kinetics of methane.

Flame velocity

The basic configuration with a length of 30 cm is used here to study the propagation speed of the laminar methane/air premixed flame.

From the numerical results, two measurement methods of the velocity have been used. The first one is based on the instantaneous position of the flame front. In this case, the front was defined as the position of the maximum of the temperature profile gradient. Thus, by defining :

x f = P os � Max � | ∂T ∂x | �� . (4.1)
Then, it is possible to estimate

S L = dx f dt . (4.2)
Cléante Langrée University of Rouen -CORIA February 8, 2022

Flames expanding through obstacles

This first estimation of the velocity is visible in fig. 4.3, square symbols.

The difficulty of such a measurement comes from the discretization of the domain, which induces a discontinuous signal for the position and thus a particularly jerky velocity information. However, it is always interesting to evaluate a visual method of analysis as experimenters can do. In the case of numerical calculations, a reliable method is to rely on the equivalent distance traveled by the burnt gases. It is possible to estimate the volume of gas burned in the domain :

V f = � V bdv = bAdx , (4.3)
with A the constant surface normal to the direction of propagation. The equivalent distance of the flame front is then introduced to the origin R f such that the volume R f • A is equal to the volume of burned gas V f . Finally :

R f = � X bdx , (4.4)
and consequently, In order to compare these results with validated thermophysical data, flame speeds have calculated with the GRIMech 3.0 mechanism and the CANTERA library. As mentionned in chapter 2, this detailed mechanism for methane-air combustion is comprised of 325 reactions and 53 species. It is originally an optimized mechanism designed to model natural gas combustion.

S L = dR f dt , ( 4
It turns out that the comparison of the three curves shows that the OpenFoam solver gives good results for the estimation of the propagation speed of a laminar flame. Let's not forget two important elements in the OpenFoam case : (1) on the one hand a regular mesh has been kept, which does not refine at the flame level as Cantera can do. On the other hand, the propagation mechanism is based only on the evolution equation of the progress variable and not on those of the 53 species of the GRIMech 3.0 mechanism.

The right order of magnitude of speed is observed whatever the equivalence ratio of the mixture. The method by tracking the position of extremum of the temperature gradient slightly overestimates the velocity and gives rather noisy results. On the other hand, the method based on the volume of gas burned in the domain gives a clear curve with a slight underestimation of the velocity but which remains very acceptable. In order to see how well the thermochemical data are mimicked, the adiabatic temperature of the burnt gases obtained

Cléante Langrée

University of Rouen -CORIA February 8, 2022 Then, an overestimation of the temperature can be observed around the stoichiometry and in rich mixtures. Let us remember that in this work equivalence ratios that varies from 0 to 1 will be considered. As will be seen in the last chapter of this thesis, the final configuration studied is the premixed combustion of a stoichiometric methane/air mixture.

It is clear that the simplified model does not perfectly match the results of a 325 reaction kinetic but for our purposes this difference is reasonable. As a matter of comparison, fig. 4.5 is particularly instructive. It was generated entirely using the Cantera solver and complex methane kinetics. The aim was to see the impact of the definition of the air mixture on the air/methane combustion, i.e. the molar proportion of nitrogen and oxygen that defines the 'air' mixture. Then, this air was mixed with methane and the flame speed analysis was performed. The velocity envelope for a molar fraction of 18 % to 25 % for oxygen has been plotted. Focusing on the maximum velocity, results range from 0.24 m/s to 0.58 m/s.

Finally a mesh convergence analysis has been performed on this 1D flame calculation. fig. 4.6 shows the reference velocity obtained with Cantera and compares it with 5 simulations performed with XiFoam for a mesh step ∆ x varying from 0.5 mm (600 cells) to 1 cm (30 cells) in the 30 cm long domain. It is particularly interesting to note the robustness of the numerical model to correctly capture the flame velocity while the number of cells is getting smaller. Let us remind that this ability to keep being consistent, even when mesh resolution is coarse, was a key criterion in our selection of a model. Indeed, the fact that the distant objective of the work is to be able to calculate premixed or partially premixed flames on very large scales with a sometimes large mesh, must not be forgotten.

Flames expanding through obstacles

Pressure and boundary conditions

In order for the reader to be able to reproduce the results presented here, it is important to give some information about major numerical features. The most important of them concerns the evacuation of the pressure waves. Indeed, the ignition and propagation of the flame will generate pressure waves of both physical and numerical origin. In both cases, it is necessary to succeed in evacuating these pressure waves at the exit limit of the domain.

For the pressure, a boundary condition called waveTransmissive was used. Boundary conditions close (but not identical) to the Navier-Stokes characteristic boundary conditions (NSCBC) applied to the pressure and velocity are approximated by waveTransmissive [START_REF] Piscaglia | Development of a non-reflecting boundary condition for multidimensional nonlinear duct acoustic computation[END_REF]. This condition was chosen, because of its capacity to modify the reflectivity of the walls. It necessitates the input of four different parameters :

• gamma, the ratio of specific heats (γ = 1.4),

• fieldInf, denotes the far-field value to be applied to the pressure (always atmospheric pressure in our case), • lInf, the distance to which the far-field condition should be applied, • value denotes the initial field pressure (always atmospheric pressure in our case).

The effectiveness of this boundary condition can be seen in fig. 4.7, where the pressure profile of the 1D flame is shown in a case where the atmospheric pressure is imposed at the exit of the domain and in another case where this same pressure is imposed at 10 domain lengths thanks to the variable lInf. It is clear that waveTransmissive achieve the evacuation of

Cléante Langrée

University of Rouen -CORIA February 8, 2022 pressure waves. Note that some reflective waves may remain, but their scale is smaller than the one of the fixed value boundary condition configuration. However, it should be noted that the setting of the distance parameter is not trivial, especially when the configuration is based on a non-structured mesh.

That concludes the 1D study using the XiFoam solver. The results obtained so far were needed to step up in the simulation of premixed turbulent combustion. Now that the solver has been proven to be at least capable of handling premixed combustion, recovering key premixed flame characteristics and even showing some consistency with mesh size variation, the next "test" configuration will be studied.

Expanding spherical flame

In this configuration, the flame propagates in a homogeneous reactive mixture among a spherical domain. The aim of this study is to assess for initial turbulence and obstacle presence impact on the flame propagation speed. Several studies in the literature worked on similar aims [START_REF] Di Sarli | Using Large Eddy Simulation for understanding vented gas explosions in the presence of obstacles[END_REF][START_REF] Di Sarli | Sub-grid scale combustion models for large eddy simulation of unsteady premixed flame propagation around obstacles[END_REF][START_REF] Coates | Numerical investigation of the effect of obstacle shape on deflagration to detonation transition in a hydrogen-air mixture[END_REF], but rarely on a fully opened configuration. Here the core interest of this configuration is to have a reference configuration, as close as possible to our goal configuration (industrial scale UVCE) to look at, without being limited by the comparison with experimental data.

Initially, a zero average velocity is set in the domain but several levels of turbulence are considered. The reaction is initiated in the center of the chamber by the deposit of a core of burnt gas. The flame propagates uniformly in all the space and its average properties such as its speed, the pressure jumps, the deformation of the flame front can then be determined at each instant.

Flames expanding through obstacles

The expanding flame is an unsteady flame. Its evolution is the result of different flame/turbulence interactions over time [START_REF] Abbasi | A scheme for the classification of explosions in the chemical process industry[END_REF]. Three main stages can be distinguished in the propagation of an expanding turbulent flame :

1. The ignition of the reactive mixture and the formation of a flame core. At this stage, the evolution of the flame is mainly governed by stretching and by heat loss by conduction. The flame front is "smoothed" by positive stretching effects. This stage is not considered in this work since a core of burned gases mimicking the early stage of the flame is considered at time 0. 2. As the flame expands, the flame is first affected by the largest frequencies of the turbulence spectrum, i.e., by the small spatial scales [START_REF] Bradley | The problems of the turbulent burning velocity[END_REF] (smaller than the diameter of the flame core). 3. As the flame spreads, it is influenced by the smallest frequencies, which correspond to the turbulence's largest spatial scales. [START_REF] Bradley | The problems of the turbulent burning velocity[END_REF] defines an effective turbulent velocity, which is lower than the turbulent velocity u � measured in the combustion chamber when there is no flame and increases as the flame propagates.

Due to the presence of the obstacles, there is first appearance of shear in the flow. This induces an elongation of the flame front, locally slowing and meanwhile accelerating in jets between the obstructions. Consequently, the flame surface is increased and distorted, which leads to a faster combustion. However, this stretching of the flame can also alter the behavior inside the reaction zone and provoke local quenching and possible partial extinctions. Finally, new structures are likely to appear behind the obstacles, which lead to classic interactions between vortices and a flame front, i.e. to an increase of the combustion. The reference geometry is a sphere of 6 meters diameter meshed in an unstructured way. Two topologies of torus-shaped obstacles have been selected : A staggered 'stag' configuration, presented in Fig. 4.8-(a), and an aligned 'lin' configuration (Fig. 4.8-(b)). In addition, two values were considered for the diameter of these tori : 4 cm and 12 cm. Up to a given major radius, the tori are repeated every 18 degrees of angle and they are all aligned with each other (lin) or an offset by a 9 degree rotation is applied every other row of torus (stag). The first row of torus has its centers placed at 0.8 m from the center of the domain. Then, the following series are placed every 30 cm, until 6 rows of obstacles have been placed in the domain, the last row being at 2.3 m from the center.

Configurations

The unstructured mesh is refined around the obstacles in order to correctly capture the interactions between the flow and the wall, in particular the wake induced by the fresh gases pushed by the flame. In these configurations the mesh size is of the order of 15 million cells.

The domain is initially filled with a methane mixture at stoichiometry. The ignition is done thanks to a core of burnt gas of 10 cm radius placed in the center of the domain. With the first tori placed at a radius of 0.8 m, the flame begins to propagate in an obstacle-free zone, allowing it to have an established spread before it begins to interact with the obstacles.

Calculations are done with the version v2012 of the OpenFoam library, using more specifically the XiFoam solver. The various tested configurations are listed in Table 4 First, a spherically expanding flame in an unobstructed environment is considered. A visualization of the flame front can be seen in Fig. 4.9 for several times. It is possible to notice the perturbation of the flame surface due to the Darrieus-Landau, i.e. hydrodynamic instabilities of combustion. In premixed combustion, these instabilities favors the formation of corrugated flames with relatively sharp edges directed towards the consumed gas. Fig. 4.9 presents the view from the fresh gases. Indeed, the shape of the instabilities shows a rounded aspect on the fresh gas side and sharp troughs towards the burnt gas. It is possible to observe identical structures on the experimental flame surfaces observed in [START_REF] Ahmed | Simulation of Spherically Expanding Turbulent Premixed Flames[END_REF].

Because of the spherical symmetry of this configuration the physical characteristics of the flame have been averaged to obtain the statistical data as a function of the radius. In this framework, Fig. 4.10-top shows the average temperature profiles at four given times. The envelopes of shaded color around each curve represent the deviation of the considered variable at the same instant. Note that the position of the flame is represented at the four instants thanks to a red vertical bar. It is possible to notice the thickening of the mean temperature because of the development of the instabilities.The more the instabilities are important the wider the envelope is.

The Fig. 4.10-bottom represents the average of the gas velocity module through the flame. The profile of the fresh gas velocity as well as the burnt gas velocity can thus be observed.

Flame radius

The turbulent combustion velocity is defined as the mass reduction velocity of the fresh gas. It is then necessary to define a specific mean radius to define the turbulent burning velocity of an expanding turbulent flame with a wrinkled flame front. The choice of the radius is critical to capture the correct flame velocity.

Bradley et al. [START_REF] Bradley | Turbulent burning velocity, burned gas distribution, and associated flame surface definition[END_REF] proposed to take the radius that corresponds to the equivalent turbulent flame radius r f where the volume of burned gas v b outside the sphere by this radius is equal to the volume of fresh gas inside this same sphere. By defining r o , the largest radius encompassing the maximum position of the burned gases and r i the minimum radius that does not contain any fresh gases, the average flame radius r f is such that :

� r f r i (1 -c) dv = � ro r f c dv (4.6)
Cléante Langrée University of Rouen -CORIA February 8, 2022

Cléante Langrée

University of Rouen -CORIA February 8, 2022

Flame expansion 93

This is generally the recommended method. However, if obstacles are present in the domain, it is no longer possible to implement it in a generic way. Another possible way is simply to calculate the radius of the equivalent sphere. In this case, the data from the mesh easily gives the inverse function r e = f (V b ) including the presence of obstacles. The volume of burnt gases V b is defined by the integral of the progress variable c over the whole domain,

V b = � Ω cdv .
Finally, it can be interesting to get closer to experimental visual methods to examine the propagation speed of the flame front. These methods are based on the visibility limit of the flame. The positions of the flame front that are farthest from the center of the domain is selected. It corresponds to r o , the 'outer' radius defined above. Figure 4.11 -Flame radius evolution depending on the selected method to compute flame front position (r f : burnt/fresh gas equilibrium, r e : volume equivalent sphere, r o : outer radius) -case free-01_00.

It is interesting to consider the three possible definitions of the flame radius simultaneously and to compare their evolution in an unobstructed case. The result can be seen in Fig. 4.11. Naturally, the outer radius r o that marks the farthest distance from the flame to the center of the sphere is the largest. The radius r f recommended by Bradley [START_REF] Bradley | Turbulent burning velocity, burned gas distribution, and associated flame surface definition[END_REF] because of its balance between fresh and burning gases volumes comes next. Finally, the smallest radius is the equivalent radius r e which corresponds to the radius of the sphere of the same volume than the volume of all the gases burned in the field. The important point is that the evolution of these three radii follows a similar slope, which means that all three correctly capture the velocity of the flame. Since the Bradley method is particularly difficult to implement in the presence of obstacles and the outer radius method may undergo abrupt changes in the position of the flame tips, the 'equivalent radius' method was chosen in this work to mark the position of the flame. Of course, the volume of the obstacles was taken into account in the corresponding cases.

Flames expanding through obstacles Figure 4.12 -Flame radius evolution depending on the mesh refinement, from 700k cells up to 7M cells -case free-01_00.

Mesh convergence

The complete domain considered in this work is a 6 meter diameter sphere in order to evaluate models that can be deployed on an industrial scale. As previously stated, the flame is described by a progress variable whose source term is associated with correlation models for turbulent combustion. The calculation of the evolution of the progress variable c allows capturing the flame surface at the scale of the considered mesh. However, the correlation model must be able to take into account the surface production (and thus the variation of the flame speed) due to the fluctuations associated to scales smaller than the mesh size. The difficulty in this endeavor is to successfully capture a flame velocity that is not mesh-dependent. Therefore, the evolution of the unobstructed flame has been calculated for 4 different meshes. The considered meshes have a number of cells varying from 700k to 7M meshes. The result can be seen in Fig. 4.12. Even if a slight variation of the flame position can be observed as a function of the global number of cells in the domain, it is clear that the flame position is little affected by the variation of the cell sizes. This shows that the Gulder correlation model is appropriate for the large scale calculations that are the long term objectives of this work. This point is confirmed in the following analysis, which shows that the turbulence level in the domain has a much larger impact.

Turbulence

In the case of a flame propagating in a turbulent environment, both the burning velocity and the turbulent flame thickness increase. Indeed, the increase in velocity fluctuations initially leads to an increase in the folding of the flame front and thus an increase in the reaction surface. Damköhler [START_REF] Damköhler | Zeilschrift für Elektrochemie und angewandte physikalische Chemie[END_REF] proposed the principle of flame front wrinkling as the main

Cléante Langrée

University of Rouen -CORIA February 8, 2022 mechanism controlling turbulent flames. Developed from his work and many researchers who followed the basics he stated, models of correlations between laminar and turbulent flame velocities are the basis of many current combustion models, including the three presented in chapter 2 (Gulder, Zimont and FLACS) and used in this work.

However, increasing turbulence does not lead necessarily to an unlimited increase in turbulent combustion velocity. Depending on the circumstances, it is possible to observe a slowing down of the increase in turbulent burning velocity and even a decrease in velocity. This effect is referred to in the literature as "bending". According to Peters [START_REF] Peters | The turbulent burning velocity for large-scale and small-scale turbulence[END_REF], the "bending" evolution observed corresponds to the transition from the corrugated flamelets to the thin reaction zones regimes. The flame surface does not increase indefinitely due to the folding of the flame front. Indeed, two zones close to the flame front can merge by contact and thus decrease the flame surface and the combustion speed [START_REF] Filatyev | Measured properties of turbulent premixed flames for model assessment, including burning velocities, stretch rates, and surface densities[END_REF]. This bending effect may be observed in Fig. 4.13 where the global turbulent burning velocity, i.e. averaged over the whole flame has been plotted. It is defined by :

S t = ρ b ρ u dr e d t . (4.7)
According to accurate DNS analysis carried out in [START_REF] Nivarti | Direct Numerical Simulation of the bending effect in turbulent premixed flames[END_REF], the observed bending effect derives from a shift in balance towards mechanisms that favor flame surface area destruction. In conditions that generate the bending effect, these mechanisms tend to preserve the reaction layer, ensuring the validity of Damköhler's hypothesis and flamelet models as shown here since the Gulder model is able to capture it. The influence of turbulence intensity on the propagation of expanding turbulent flames is investigated in Fig. 4.14 for two different initial turbulence intensity levels in combination with a mesh sensitivity analysis. The objective here was to verify the maintenance of mesh convergence for different levels of tubulence. Clearly, as shown in Fig. 4.14, turbulence does have a strong impact on the flame speed. Nevertheless, the combustion correlation model seems to be able to handle the variation of turbulence in the domain regardless of the mesh size. The Gulder model was used in this example but the conclusion can be extended to the Zimont or FLACS model. Now that the basic behavior of this large flame (final diameter : 6 meters) have been analyzed in an unobstructed environment, it is time to move on to the heart of this study, which is to analyze the impact of the presence of obstacles on the propagation of the turbulent flame.

Flames expanding through obstacles

Impact of obstacles

Accidents associated with gas or even dust explosions usually take place in an industrial environment surrounded by obstacles of widely varying size and shape. The presence of these obstacles has a great influence on the speed of the flame during its propagation. The first effect is an acceleration of the flame, which generates even greater damages associated to the flame propagation but also to the increase in the resulting pressure wave or an even faster transition from deflagration to detonation. If the acceleration of the flame within the obstacles is a known phenomenon observed in many studies, the ability to take into account all the phenomena that come into play remains limited.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

Flame expansion 97

The propagation of the flame generates a high velocity of fresh gas upstream of the flame, which is pushed around the obstacles. A turbulent wake is naturally formed and it is possible to ask which process dominates between the acceleration of the flame due to the increase in turbulence behind the obstacles or due to the deformation and stretching of the flame front by the wakes.

Obstacle layout

As detailed previously, two different arrangements of obstacles were used. In the first case, the obstacles are aligned one behind the other (in-line arrangement) in six rows, while in the second case, the obstacles are staggered every other row. Moreover, two different radii of obstacles were considered : 2 cm and 6 cm. The combination of arrangement and radius lead the following configuration names : lin2, lin6 and sta2, sta6. fig. 4. [START_REF] Pitsch | Cefrc combustion summer school -modeling turbulent combustion[END_REF] shows two examples of the propagating flame front in the case of the staggered obstacle arrangement for the two different diameters : sta2, sta6. This arrangement implies that the flame front, initially developed in an unobstructed space, systematically impacts an obstacle. On the other hand, the in-line arrangement means that the flame front developing in the free space between two obstacles is free to evolve throughout the propagation of the flame and will be much less affected by any wake effect.

Nonetheless, it must be kept in mind, that for a given radius of obstacle, the volume obstruction rate remains the same regardless of the obstacle arrangement. The flame propagates through the obstacles, which naturally generates a strong deformation of its front along the azimuthal plane (plane perpendicular to the obstaclesfig. 4.16 -left). However, it is possible to observe flame instabilities in the elevation plane (plane parallel to the obstacles, fig. 4.16 -right). As for the free flame of fig. 4.9, this is a Darrieus-Landau type instability.

In fact, the flame front is not at all the same depending on the arrangement of the obstacles. When the obstacles are staggered, a "flower" type propagation with successive curves formed by the crossing of the obstacles can be observed. On the other hand, in the case of obstacles in line, the flame presents a very regular structure with a regular, almost flat front between each obstacle. However, it is particularly interesting to note that the evolution of the flame surface in both cases is almost identical. Indeed, fig. 4.18 represents the evolution of the flame surface in the five main configurations. It turns out that for a given obstacle diameter (i.e. a given rate of obstruction), the curves of evolution of the flame surface coincide. First, the free flame (free-01_00) evolves until it reaches a surface of 400 m 2 , which makes it the least developed of all. Then in the case of the 2 cm radius obstacles, the flame surfaces evolve in a similar way from end to end to reach a maximum of 800 m 2 . Finally, in the case of the 6 cm radius obstacles, the flame areas are initially very close and end up with a (relatively) small gap : 1300 m 2 and 1550 m 2 respectively for the staggered and aligned cases.

The flame area directly reflects the amount of fuel consumed and therefore the speed of the flame. This similarity in the evolution of the flame area translates directly into a similarity in the equivalent flame radius over time (see fig. 4.19). If one observes the slope of the equivalent radius, the free flame initially shows a natural acceleration and then the speed reaches a maximum before slightly decreasing. As already introduced, it is the bending phenomenon described by Peters et al. [START_REF] Peters | The turbulent burning velocity for large-scale and small-scale turbulence[END_REF]. On the other hand, the four cases with obstacles show a continuous acceleration of the flame, rather weak with obstacles of 2 cm radius but quite visible with those of 6 cm radius. The decrease of the flame radius at the end of the simulation is directly linked to the absence of obstacle in the far field. In summary, the presence of obstacles causes a continuous increase in the speed of the flame. This acceleration is directly linked to the size of the obstacles (or obstruction rate), while the arrangement of the obstacles has only a weak impact on the flame surface and, therefore, propagation. It may be interesting to look at the structure of the flame to understand this last point.

Flame structure

One way to study the flame structure from a statistical point of view is to start by analyzing the evolution of the position of the flame surface, which is defined as the value 0.5 for the progress variable in this work. Though this value may seem arbitrary, it has been observed that changing it has only a marginal impact on the presented and discussed results.

Flames expanding through obstacles

In order to obtain more information about the flame surface topology, the orientation of the flame surface is also observed. For this purpose the angle β have been defined as the angle between the normal to the flame surface and the vector issued from the origin. Thus, a perfect sphere surface normal present a null β angle, which will increase with any variation of the surface orientation. The PDF of the angle beta β for the case of the free flame has been plotted in fig. 4.21. Here again, it is possible to observe a unimodal distribution, but with a right skewed form. Indeed, the mass of the distribution is concentrated on the smallest angles (i.e., the smallest deviation of the surface since the maximum position varies from 15 o to 25 o . This mostly reflects a moderate deformation of the flame compared to an equivalent sphere. Nevertheless, the tail of the distribution on the right gives us valuable information about the fact that the flame forms deep troughs towards the burnt gas and in these cases the angle of the normal can rise up to 60 o . On the fresh gas side, the flame has to push a heavy gas which will tend to dampen its possible deformations while on the burnt gas side, the high speed of the latter will promote the formation of these steep 'hollows'. Similar conclusions have been drawn by Ahmed et al. [START_REF] Ahmed | Simulation of Spherically Expanding Turbulent Premixed Flames[END_REF].

Staggered vs in-line obstacles

Now that the statistical data of the evolution of the free flame can be used as a reference, let us observe the evolution of the PDF of the flame position in the cases with obstacles. Fig. 4.22 represents the evolution of the PDF of the flame position in the sta6-04_00 case (staggered 6 cm radius obstacles). It is clear that as soon as the flame reaches the obstacles, it loses its unimodal shape for a multimodal shape with a growth of the number of peaks with time.

Cléante Langrée

University of Rouen -CORIA February 8, 2022 This three-peak structure is also found in the case where the obstacles are aligned, as depicted in fig. 4.23, orange curve. It is also worth mentioning that, in both configurations, a fourth peak in the probability of the flame front is present in the burned gases. This corresponds to the combustion of small pockets of fresh gas that have not yet burned completely around the obstacles of the previous row; their PDF is however three to four times lower, the appearance of pockets not being systematic around every obstacles.

Even if the general shape remains the same, the probability of presence at the level of the central peak (behind the obstacles) is more important in the staggered case (Fig. 4.23, blue curve, 2 < r < 2.2) because of the 'flower' shape of the flame (Fig. 4.16). This shape is directly due to the systematic encounter of the flame with an obstacle in this configuration. On the other hand, in the aligned case, part of the front never encounters any obstruction, which is why probability of larger radius is more important on the right side of the curve, at the limit of the fresh gas (Fig. 4.23, orange curve, 2.3 < r < 2.5). of having a 'lateral' front is much lower (Fig. 4.24, orange curve). It even shows a plateau between 10 o and 40 o . On the other hand, the flower shape of the flame in the staggered case induces a more regular progression of the curve and a higher probability of having a front directed towards the fresh "outside" gas. Going back to the in-line case, after the 'attack' front, which is almost oriented according to the radius, it is possible to observe (Fig. 4.17) a particularly developed 'lateral' front, which will propagate between the bars and eliminate the pockets of fresh gas not consumed by the attack front that did not encountered any obstacles. This front is at the origin of the main probability peak of the orange curve of fig. 4.24. This peak exceeds in probability, and in angle, the one of the flame of the staggered case (70 o against 60 o ).

Finally, the third zone of the PDF P (β) concerns angles greater than 90 o , i.e. flame fronts oriented towards the interior of the domain. In both the staggered and in-line cases, the shape and probability level are similar. This is in fact what could be called the "return" front. Its origin is the flame front which, after having bypassed the obstacle, at a more or less important distance from it, will propagate towards the "back" to consume the remaining pockets of fresh gas.

Obstruction rate

Apart from the arrangement of the obstacles, another major parameter that will influence the propagation of the flame is the size of the obstacles. Fig. 4.25 represents the PDF P (r ) of the position of the front of three flames, all with an equivalent radius of 2 m : the free flame (without obstacles), the flame with staggered 2 cm radius obstacles and finally the flame with staggered 6 cm radius obstacles. It is this last flame that was studied in the previous section for comparisons of flame structure as a function of obstacle arrangement. First, it is possible to recognize the aforementioned quasi-Gaussian shape of the free flame distribution (blue curve). In this configuration, only the Darrieus-Landau instabilities are responsible for the flame shape. It is also important to keep in mind that the equivalent radius that serves as a localization reference (Fig. 4.11) is different from the average localization of the b = 0.5 isosurface that was used to study the flame structure, i.e., the PDF average is not equal to the equivalent radius.

The green curve in fig. 4.25 represents the probability of presence of the flame front in the case of staggered 6 cm radius obstacles. This case has already been discussed previously with the three main peaks corresponding to the three stages of flame progression in the obstacles : (1) upstream obstacle front (2) downstream obstacle front and (3) advanced front that has not been influenced by the obstacles. In the case where the obstacles are arranged in the same way but with a radius 3 times smaller, i.e. a volume 9 times smaller, the shape of the PDF P (r ) presents a new flame structure (orange curve of the fig. 4.25). Indeed, the PDF is formed by a main peak with a high probability of presence at the level of the obstacles, followed by a plateau which corresponds to a rather uniform distribution of the flame over about 30 cm. Because of the lesser presence of obstacles, this shape (orange) can be seen as an intermediate step between the Gaussian like distribution (blue) of the free flame and the three-peak distribution (green) of the case with large obstacles. Therefore, it is clear that the presence of obstacles 'spreads' the flame front and implies a dissimilarity of distribution. It should be kept in mind that several phenomena are in competition during the propagation of the flame. Three of them can be highlighted : [START_REF] Trelat | Impact de fortes explosions sur les batiments representatifs d'une installation industrielle[END_REF] turbulence which will locally increase the flame area as well as (2) hydrodynamic instabilities and finally (3) the presence of obstacles. However, it may be interesting to note that if the

Cléante Langrée

University of Rouen -CORIA February 8, 2022 As for the position distribution, the flame surface orientation distribution P (β) of the case sta2-02_00, visible in fig. 4.27, can be seen as the intermediate distribution of the free flame free-01_00, which presents a uni-modal PDF with a right skewness (blue curve), and the case with large obstacles sta6-04_00 (green curve), which has been described previously. For the case of interest, the shape of the distribution shows a more important 'attack' area (angles between 0 o and 45 o ). This is due to three reasons : firstly, the impact of the 'small' obstacles on the flame structure is much less important and the flame is less deformed and tends to remain 'oriented' towards the outside. Secondly, the free space between the obstacles is larger. Therefore, even if the flame deforms at the obstacle, it will tend to damp this deformation, thus decreasing the angular deviation. Finally, as mentioned above, the thinness and the distance between the obstacles will also allow the natural instabilities of the flame to develop. These instabilities, in the case of the free flame, maintain a flame orientation with a relatively small angle. Finally, the last part of the P (β) distribution, which concerns angles greater than or equal to 90 o (the 'return' front), is much weaker and has a similar shape to that of the free flame.

Eventually, the last part of the P (β) distribution which concerns angles equal or greater than 90 o is much weaker and has a similar shape to that of the free flame. Indeed, the bypassing of small obstacles remains statistically a minor phenomenon compared to the volume considered.

Flames expanding through obstacles

Combustion correlation model

A final set of analyses was considered. It should be seen as an introduction to further work. It concerns the correlation models for turbulent combustion that have been introduced earlier. Many researchers ( [START_REF] Dahoe | On the Application of the Levenberg-Marquardt Method in Conjunction with an Explicit Runge-Kutta and an Implicit Rosenbrock Method to Assess Burning Velocities from Confined Deflagrations[END_REF] and references therein) have attempted to establish correlations, mostly empirical, expressing the turbulent combustion rate as a function of the characteristics of the reactive gas flow. In the simplest models, it is considered that the turbulent combustion velocity is directly related to the velocity fluctuations, however more complex models express this correlation as a function of the Lewis number, the Damköhler number or a Reynolds number, with the aim of establishing a unified law to describe a large number of turbulent premixed flames.

However, all these models are based on the assumption of an equilibrium between the reaction front and the turbulent flow. More clearly, they do not take into account the possible presence of obstacles which will modify this equilibrium or even prevent its establishment. fig. 4.28 shows the evolution of the equivalent flame radius in the sta6 configuration (staggered obstacles, 6 cm radius). Three correlations are used : Gulder's correlation (by default implemented in OpenFoam), Zimont's correlation (which can be found in softwares such as Ansys FLUENT) and finally FLACS correlation, present in the eponym software, which is very much used in the context of industrial safety. The three curves represent the position of the flame for the three correlation models. It is interesting to note that in the free propagation zone before the first obstacle, the average flame position is the

Cléante Langrée

University of Rouen -CORIA February 8, 2022 On the other hand, as visible in fig. 4.28, after passing the first obstacle, symbolized by the gray horizontal line located at r = 0.8 m from the center, the three models have a different behavior and the equivalent radius of the flame varies. The Gulder model presents the fastest flame, followed by the 'Zimont' flame and finally the 'Flacs' flame. However, it should be noted that if the flame speeds differ, the overall behavior of the flame remains similar.

As it have been done for the other analyses, it is possible to plot the statistical structures of the flames when they reach the equivalent radius of 2 m. In fig. 4.29, is plotted the radius distribution of the flame position obtained with the Gulder model that have already been analyzed (blue curve) with its three characteristic peaks corresponding to the flame front in front of the obstacle; to the one downstream of the obstacle and finally to the flame front that, having been little or not at all disturbed by the obstacle, has propagated further into the free zone. It is interesting to note that the structure of the flames obtained with the two other correlation models present significant differences. The flame from the Zimont model (in orange) is quite close to the Gulder flame in the first two peaks. However, the statistics for the presence of the flame front in the free zone are much lower. On the other hand, they are more important than those of Gulder upstream and downstream of the obstacle. In the case of the Flacs model (green curve), an opposite phenomenon can be observed with a flame that is prevalent in the free zone and very present in front of the obstacle. The first upstream peak on the other hand is almost non-existent in this example. It has been verified that this is not a slight shift in the progression of the flames presented but a different structure depending on the models. Indeed, as shown in fig. 4.30, the flames progression at the obstacles is different, which implies these differences in the statistical structure of the flame. In fact, with a 'lateral' front prevailing for the Flacs correlation, especially in the range of angles between 45 o and 90 o , this angle implies that the flame has a strong tendency to propagate diagonally in the wake of the obstacle. The efficiency of the 'lateral' front implies a high consumption of fuel downstream of the obstacle, which limits the formation of a 'return' front (angles > 90 o ) as shown in the statistics of the green curve for this angle range. On the contrary, a more important 'cleaning' work carried out by a 'return' front is observed in the case of Zimont and Gulder correlations.

Flames expanding through obstacles

In conclusion, if in the case of free propagation flames, these correlation models have the same behavior. It is unsurprising since they already demonstrated their efficiency in many other works. Nonetheless, in the presence of obstacles, some discrepancies appear. This may be related to local modifications of the flow that affect the model parameters. At this stage, it is impossible to say if one model is more suitable than another and it would be necessary to rely on direct numerical simulations (DNS) or experimental results. Concerning DNS, the scales of the order of ten meters and more concerned by the safety studies cannot be reached. And to obtain a validation of a model of the order of a millimeter or a centimeter would not mean it is valid for domains a hundred or a thousand times larger -even if precious information could be gathered. The experimental point of view is the most promising however with many blocking points and notably the difficulty of measurements in the presence of obstacles. Indeed, if it is straightforward to obtain results on the global position of the flame

Cléante Langrée

University of Rouen -CORIA February 8, 2022 front, the adaptation or the development of adapted models requires to know the evolution of the physical parameters in the heart of the configuration, which is difficult to reach by the measurement tools. The aspect of the behavior of turbulent correlation models within an obstacle should be the subject of future studies.

Conclusion

In a first part, a simple 1D configuration have been studied with the OpenFoam solver to calculate the propagation of premixed flames. This allowed us to evaluate the performance of the solver in a purely thermochemical framework. It turns out that the flame velocity is correctly estimated for different equivalence ratio (bell-shaped curve) and in particular for the lean or stoichiometric mixtures that interest us.

Then, a configuration developed to analyze the propagation of a spherical flame was specifically developed. Analyses have been carried out based on the temporal evolution of the flame, on the geometric position of its front and finally on the statistical data of its structure. The following findings could be put forward :

• Different shapes are observed for the flame front, depending on the arrangement of the obstacles, however the time evolution of the flame surface is almost identical in all cases.

Flames expanding through obstacles

• The impacts of the presence of the obstacles are mainly twofold : i) they reduce the development of the instabilities, ii) they modify the flame structure by increasing its surface. All-in-one, the corresponding overall impact tends to be a strong continuous increase of the flame speed. • Though the evolution of the surface remains the same, important modifications occur due to the interactions with the obstacles. First, the flame tends to loose its unimodal shape for a multimodal shape with a growth of the number of peaks with time (three peaks are observed after the first line of obstructions). Secondly, both configurations show the appearance of a fourth peak in the probability of the flame front (with a magnitude three to four times lower), which corresponds to combustion of small pockets of fresh gas at the rear of the obstacles. Thirdly, the PDF is more important at the central peak in the staggered case. Fourthly, in spite of keeping an unimodal shape as in the obstructions-free case, the PDF of the flame surface orientation presents modified maxima and three mains zones. Thus, there appears a dissimilarity of distribution. • Concerning the turbulent combustion models, it has been clearly highlighted that the assumption of equilibrium between the reaction front and the turbulent flow can raise some issues. More important, significant and important differences are obtained in the flame structure depending on the correlation used for the turbulent combustion rate (as a function of the features of the reactive gas flow). This would clearly need further investigations to identify the origin of such divergences, and to bound the correctness of each model in case of flame propagation in presence of obstacles.

Chap. 5 | Industrial Application

Introduction

The majority of methods to evaluate the consequences of an explosion are based on experimental results [START_REF] Paris | Évaluation des effets d'une explosion de gaz à l'air libre. base documentaire[END_REF]. However, these empirical methods quickly find their limits in the case of gas explosions. Indeed, they do not take into account many factors influencing the violence of the explosion such as the nature of the gas, the turbulence, the confinement by the walls or the obstruction generated by the obstacles. This is why a lot of research has been undertaken to better understand the mechanisms involved and thus improve the methods for predicting the associated effects. The availability of ever more powerful computers has allowed the development of specific modeling tools, in particular through CFD calculations which now make it possible to take into consideration a large part of the physical phenomena related to a gas explosion. The characterization of the pressure levels undergone by individuals or by structures can now be characterized thanks to the quantified risk analysis approach.

Industrial Application

The work presented in this manuscript is a step towards the objective of CFD modeling of deflagrations and detonations for UVCE at the core of an industrial site. The current power of computers is still limited compared to the scales to be considered, but this objective is not very far from being achieved in a few years. As part of this work, these future configurations are are being prepared by studying the physical and numerical models that will be able to meet this challenge. This work follows an incremental progression, as shown in figure 5.1. After evaluating the laminar flame kinetics using a 1D configuration of the order of the centimeter, the impact of obstacles on a deflagration have been studied in the case of a turbulent flame propagating in a medium at the scale of the meter. The "academical" configuration used allowed us to extract converged statistical data. This gave us a certain amount of information with respect to various physical phenomena and in particular, the turbulence, the instabilities of flame and the confinement of the latter within an obstacle network. This chapter is the last step of the work: it involves making qualitative and quantitative comparisons, at the scale of several meters, of these models with experimental measurements such as those used by industrial safety specialists to develop empirical estimates. The experimental UVCE configuration that interests us is a simplified version but representative of an industrial accident. This kind of experimental UVCE have been triggered in a controlled environment so that they will not harm anyone or valuable material. The environment is usually equipped with many sensors (for pressure, temperature, gas concentration etc...) and the operators can apply some parametric variations to fully characterize various physical phenomena.

EXJET experiments

The experimental work considered here has been carried out by INERIS and GDF-SUEZ ( [START_REF] Sail | Review of knowledge and recent works on the influence of initial turbulence in methane explosion[END_REF]). The aim of the experiments was to study methane jet explosion in a congested environment. The congestion module has a medium size (3 m × 1 m × 0.45 m) obstacle module, consisting of 296 steel tubes of 2 cm diameter. Each obstacle is separated by a distance of 12 cm from the others. The sketch of the experimental setup is visible figure 5.2 and a portion of the obstacle module is visible in figure 5.3. . In the framework of the EXJET experimental setup, several experiments have been carried out :

1. Ignition of a steady methane-air mixture : this setup is the first study step. It consists of applying a thin film of polyethylene over the metal cage formed by the network of obstacles. In this protected space, a stoichiometric methane / air mixture is confined before ignition is caused by the operator. The flame propagates within the domain and the polyethylene sheet rises and then moves away as the pressure increases. This configuration allows to focus on the effects of obstacles on combustion without interfering with other parameters such as partial mixing with air or varying levels of turbulence. These two points concern the following two configurations.The mixture is ignited by an electrical ignition device (location 1 in Fig. 5.2) located inside an aluminum box (19 cm × 15 cm × 13 cm) with a lateral opening. The ignition box is visible in the left part of the figure 5.3. The box is slightly shifted from the module y axis (Fig. 5.2-top green square) so that the primary flame, exiting from the box, will ignite the flammable mixture along the main axis of the module. 2. Methane jet dispersion within the module (no ignition) : it allows to mimic the methane/air mixing in case of an accidental release situated at 2.35 m from the obstacles, whose length in the direction of the jet is 3 meters.This configuration allows to produce initial experimental conditions close to the reality of an industrial accident. First, it generates a partially premixed mixture due to the diffusion and mixing of the methane jet with the free atmosphere. Moreover, the level of turbulence is correctly represented. 3. Methane jet ignition within the module : it correspond to the ignition of the previous mixture within the obstacles in an open environment.

Industrial Application

In order to proceed logically after the statistical studies of the previous chapter, a numerical simulation of the EXJET configuration has been proposed. The configurations 2 and 3 are quite possible from a domain dimension point of view. However, correctly capturing the formation of a mixing zone downstream of a turbulent jet is a particularly difficult task that involves significant numerical work on injection, mixing layers and turbulent transport within the jet. Such a simulation would have led to results that could not be validated because neither the velocity profiles in the jet and obstacle grid nor the local equivalence ratio of the mixture have been identified experimentally.

Thus, the first setup with a controlled environment has been reproduced. The experimental measurements of this configuration are based on two types of measurement :

• the measurement of pressure variations : several pressure probes are positioned in and around the obstacle module as in Fig. 5.2. Six KISTLER 0 -2 bar piezoresistive differential pressure sensors mounted on lenticular mounts were used to measure incident airborne overpressure. These sensors have an uncertainty range of 0.1 percent of full scale or 2 mbar . The overpressure sensors are all positioned in the same way.

Image analysis

A first step was to process the raw images obtained by the fast camera. Initially in black and white with a low contrast, each image was processed with the Photoshop software and an algorithm based on artificial intelligence. A colorized image has thus been obtained, which allows to differentiate the five elements of the scene : the tarpaulin stretched in the background, the cage formed by the obstacles, the polyethylene film, the flame and the grassy foreground. All the colors respect reality except for the flame which has been represented in white for a better contrast. The figure 5.4 shows both initial and final experimental images.

After the image processing, a first visual analysis of the ignition and propagation of the flame in the field has been performed.

The first image corresponds to a time of 0.008 s after the ignition. It is the first sign of flame appearance coming out of the ignition module and reaching the obstacle module (Figure 5.5 -top). A red circle has been added around the white dot that is the flame at this time.

At time 0.02, corresponding to the central image of the figue 5.5, the point has evolved into a white plume oriented to the right of the domain (direction y of the configuration). The flame has started its propagation process and it can be seen crossing the first row of obstacles at this moment.

It is important to note that the light intensity of the white color representing the flame has not been manipulated from one image to another. Since this scale is constant on the whole 2000 image constituting the database, the flame front in the third image of the figure becomes almost invisible. It is however possible to see a white fog that has also spread according to the height and depth of the cage. In particular between the first two rows of obstacles (counted according to the horizontal direction y ). This moment was chosen because a careful analysis of the image shows the first sign of the flame passing after the second row of bars. Its perceived position is surrounded by a red dashed ellipse. Depicted in fig. 5.5, 5.6 and 5.7 are instant images taken during the experiment at 6 moment in time. Those 6 instantaneous pictures illustrate issues and phenomenon that have to be taken into account in order to interpret the experimental data.

After the ignition phase, the second set of selected images in figure 5.6 shows a first phase of propagation in the domain. The most notable feature of this phase is the low luminosity of the propagating flame. The front is visible but it requires an effect to differentiate it from the other elements of the image. This propagation takes place in all directions. According to the main axis y of the obstacles, according to the depth of the image, which cannot be perceived and finally according to the vertical. The ultra-thin film of polyethylene rises up in front of the flow. Slightly at first at time 0.09 s then in a more pronounced way. The film finally leaves the ground as shown in the image corresponding to the time 0.11 s. The edge of the film is indicated by a green dashed line. Finally, the last series of photographs in figure 5.7 shows the second and last phase of flame propagation. In this phase, the flame becomes clearly visible, showing a marked intensity. On the figure, it is highlighted by a red dashed line. The flame propagates beyond the whole obstacle, in an area outside our study area. The polyethylene film continues to rise without any apparent blockage.

Flame velocity

Two methods were used to measure the flame front velocity. In cases where the flame front is clearly visible, INERIS determined the position of the flame front using automated processing of the raw camera images (Fig. 5.4-left). This consists of determining the position of the flame front by subtracting two images obtained at consecutive times. Using the differential obtained, it is possible to determine the maximum position of the flame front along the y propagation direction. However, this method could not be implemented due to the high visibility of the flame. The image processing software was no longer able to distinguish it from the digital noise. The 'human' operators of the experiment therefore directly determined the position of the flame front by eye. They then obtained the yellow triangles curve of figure 5.8. Please note that the time is not the time of ignition in the ignition box but rather the time of the appearance of the flame in the obstacle network.

It is possible to recognize two main phases in the progression of the flame front : a first phase -quite irregular -from 0 to 0.11 s with a rather low propagation speed (around 7 m/s) and then a second phase from 0.12 s to 0.2 s with a much cleaner signal and a constant and stronger propagation speed within the obstacles (around 25 m/s).

Two elements first caught our attention. First, the unevenness of the signal at the beginning of the propagation and second, the clear change in the speed of propagation of the flame around 0.11 s. Regarding the first point, the diagnosis is clear: the lack of contrast in the raw image makes it particularly difficult to determine the position of the flame front, which is almost invisible to the naked eye. Therefore, a second extraction of the flame front position from the AI colorized images has been performed. It is the blue square curve of the same figure 5.8.

With this new position curve determined with a much more visible front, the signal obtained is much cleaner. This inflection is still observed in the curve around 0.1 s even if it is less sharp than before. The major interrogation remains : what is the origin of this inflection ? It is possible to record the position from which the inflection takes place, this corresponds to a distance around 50cm from the origin (located at the beginning of the obstacle network). This corresponds to more than four rows of obstacles.

From the results obtained in the previous chapter, the passage of a freely expanding flame arriving in a zone with an obstacle network have been observed to provoke an immediate increase in the speed and thus, an inflection of the curve of the flame front position. However, once the first row of obstacles was crossed and the flame had adapted to the new environment, the change in velocity remained very small and no inflection was visible in the calculations.

Cléante Langrée

University of Rouen -CORIA February 8, 2022 This is not what is observed here. The inflection takes place in the heart of the obstacle network. A possible explanation would be that the numerical model used in the previous chapter was incomplete and that a physical phenomenon not taken into account, such as the production of turbulence energy at the surface of the flame or particular wake/flame interactions, cause a marked increase in the speed of the flame. But in this case, how the quasi-constant speed after the inflection, when dozens of obstacles are crossed, can be explained ?

In our opinion, one of the possible reasons for the speed variation at this moment is to be correlated with the experimental device. Indeed, as previously described, the combustible mixture is initially prepared under a polyethylene film which is maintained around the obstacle network by a specific frame a little wider than the latter. Only gravity holds the film that covers the top and edges of the cage, but the weight and total surface area of this film is not known. Even if this film is as light as possible, it may be important in the initial phase of flame propagation. There are two main phases of the film: the first one where it starts to lift but does not break the separation of the mixture from the outside air, and the second one where the edge of the film that was lying on the ground lifts off and then lifts completely during the flame propagation.

The instant when the film lifts off the ground corresponds to the instant of the inflection in the flame speed and it is difficult not to correlate the two even if the reason for this acceleration cannot be clearly determined. A priori, two explanations can be proposed : (1) a variation in the propagation of pressure waves under the film when it rises or (2) a mixing effect with the fresh outside air that is able to circulate in the obstacles. However the mixing being initially to the stoichiometry, the impact could be done via an increase of the turbulent mixing between the fresh and burnt gases. This unconfined mixing, now that the film is lifted, would have stimulated the general propagation of the flame. The temporal order of magnitude of the flame propagation is 1/10 of a second. If the speed that the air would have to reach to cross an inter-obstacle distance of 12 cm is considered, an order of 1.2 m/s is obtained for the air velocity, which is not at all impossible depending on (1) the external atmospheric conditions during the experiment and (2) the impact of the polyethylene film lifting on the flow.

Pressure signal

Measurements of the pressure signal within the obstacle network were made during the experiments. The three positions of the sensors according to the main study direction y are: 0.28 cm, 0.98 cm and 1.96 cm (Fig. 5.2-L3, L4 and L5).

It is important to note that these sensors are located at ground level at a zero height. This information is important because the passage of the flame is measured according to its maximum y position which corresponds in fact to a height of half the height of the obstacle module (around 20 cm). It is possible to notice that the signals do not show a marked difference in their shape after 0.1s, instant when the film lifts. However, the maximum overpressure increases by 50 % from 4 to 6 mbar . With the little experimental data in our possession it is difficult to affirm that this is a direct consequence of the film lift off but it is important to underline it nevertheless, just like the change in flame speed.

Numerical Simulation

The objective of the study is to take the models and methods presented in chapter 2 & 3 and tested in chapter 4 and to apply them for the simulation of the EXJET configuration. As previously said, the main points of interest here are the ability of the model to reproduce the flame propagation speed, and to predict a correct overpressure signal, coherent with the deflagration associated.

Cléante Langrée

University of Rouen -CORIA February 8, 2022 The first step is to consider the different elements of the geometry in order to generate a suitable mesh. Two different geometries have been worked on : the one that represents exactly the experimental configuration as represented by the planes and a simplified version. Indeed, the initial configuration, which seems particularly simple at first sight, presents some difficulties which make the calculations more difficult. The obstacle network is composed of steel tubes of 2 cm diameter which are placed one on top of the other for an obvious ease of construction. Such a geometry have been used in the simulations, but with the drawback of having to refine the mesh at the level of the points of contact between the bars, which only touch each other at one point when there are two tubes and at three points when there are three tubes (see figure 5.10). Note also that the micro-cavities formed at the intersection of the bars can be difficult to apprehend by the meshing tool and generate errors. Besides, a high precision of flow in the crossover zones is not of interest with respect to the objectives of the simulations which wish to capture the right flame speed according to the global congestion rate. Therefore, the obstacle network has been modified in order to avoid the formation of these difficulty points. By doing so, a lot of computational resources is saved by avoiding too much meshing of secondary areas. Both meshing configuration are depicted in fig. 5.11.

Ignition box

The second element of the simulation that have been simplified is the flame ignition box. From an experimental point of view, this one is designed in such a way that the flame front that crosses the first row of obstacles along the main direction of the study (y ) is as little disturbed as possible and that the flame propagation is correctly established according to the mixture richness parameters. This is why the ignition is carried out in a lateral box as shown in figure 5.12. As for the crossing of the bars, an ignition by implementing the complete lateral box has been tested. No major difference on the flame propagation in the obstacle network could be observed. The only requirement is that the origin of the times must be reconsidered. Since using the ignition time is no longer an option. That is why the instant when the flame crosses the first bar of the obstacle network is now the zero time.

Boundary and initial conditions

The boundary conditions are the same as those used in the study of the spherically expanding flame in the previous chapter. They are summarized in the table 5.1. Concerning the initial settings, the velocity field is zero as well as the relative pressure everywhere in the domain. Then three phases can be distinguished. First the flammable mixture, inside the polyethylene film is represented by a mixing variable f t = 0.055 (the methane mass fraction of the mixture, at stoichiometry) and a regress variable b = 1, as the mixture is still not burnt. Then, as described previously, a sphere of burnt gases is introduced in the first phase to ignite the mixture. In this sphere, b will be initialized at 0, as the sphere is already burnt and f t , being a mixture variable, is unaffected by the combustion and is initialized to 0.055. Finally, the ambient air phase is initialized by setting to the domain outside the polyethylene film a value of regress variable b = 1, and mixture variable f t = 0.

Global flame propagation

Therefore all the physical and numerical models presented and tested in the previous sectio have been used to simulate the EXJET configuration. In a first step it is possible to compare the numerical and experimental flame fronts from a qualitative point of view. Figure 5.14 shows the flame front at time 0.03 seconds, just as the flame has crossed the first row of bars and the inter-obstacle space that follows. At this moment, the flame still has its spherical shape both on the experimental and numerical images. A second encouraging point is the general shape of the flame above the obstacle. Indeed, the pear shape that can be observed numerically (bottom) is also found at the level of the polyethylene film (top). This film does not delimit the flame front, but the moving fluid that pushes the film away is directly the result of the advancing front. It is therefore quite possible that the shape of the film is a replica of the shape of the front at a given distance and the flame front shape obtained with the simulation is quite similar.

Industrial Application

In figure 5.13 two examples of temperature and pressure profiles through the flame have been plotted. On this curve, the distance from 0 to 3 m is where the obstacle is located. It is possible to see along an arbitrary line that the temperature variation remains very close to a classical profile with a strong gradient at the flame. Similarly, the pressure profile shows a clean curve with an amplitude of 8mbar of variation. The pressure peak, from a spatial point of view, is located just upstream of the flame front. It is clear from these pictures that the overpressure detected by the probes is "generated" by the propagation of the flame front, as the pressure signal takes the "N" shape described theoretically in chapter 2.

Another phenomenon previously mentioned can be identified in these two figures. At 0.025 s after the ignition, the pressure on the burnt gases side is significantly above atmospheric pressure (3 mbar ). As said previously, at that stage of the propagation, the flame front is still quasi spherical and a central symmetry can be considered. On the other hand at 0.1 s, the flame front is not spherical anymore and has reached the end of fresh gases in the opposite direction to the obstacle module. This means that burnt gases has a way of evacuating the pressure generated by thermal expansion on this side of the burnt gases volume. This results in a return to atmospherically pressure in the burnt gases, which can be seen in figure 5.13-right.

Flame velocity

Let us first look at the position of the flame front over time. The curves of the flame front position according to the images used by the fast camera have already been seen in figure 5.8. For this comparison, the curve from the clearest data (blue squares) is used. It has been plotted again as a continuous blue line in figure 5.16.

The curve obtained by the CFD simulation is plotted in orange. Two fuchsia colored lines were drawn at the end of each curve. The two lines have exactly the same slope corresponding to a similar flame speed equal to 25 m/s. This numerical result is a particularly encouraging one with respect to the different choices and developments that have been made so far.

On the other hand, the numerical simulation have not captured any inflection of the flame front position in the first third of the obstacle block, as what have already been discussed in the description of the experimental results. Indeed, the numerical calculation shows a flame that crosses the obstacles at an almost constant speed all along. The same behavior was observed when the spherical flame was spread in the torus network in the previous chapter. This shows that at this stage the numerical model does not take into account the phenomenon causing this speed variation. Indeed, the polyethylene film is totally ignored in the numerical calculation. If, as it is supposed, the presence of this film is at the origin of the inflection, it is normal not to capture it. Nevertheless, it is particularly positive to find the right flame speed between the calculation and the part of the experiment that totally rejected the film upwards. As a result, its presence no longer has an effect on the longitudinal propagation of the flame. It is possible that it still has an effect on the vertical propagation of the flame, but this is not considered because of the difficulty of delimiting the position of the flame front behind the polyethylene film.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

Industrial Application

Then, the flame front dynamics is evaluated in the obstacle module, based on the correlation model for turbulent flame speed. The result can be seen in figure 5. [START_REF] Ahmed | On the stress-strain alignment in premixed turbulent flames[END_REF]. The blue curve represents the Gulder correlation that used by default.

The green curve represents the FLACS correlation, used with the coefficient C F = 1.536, determined in the previous chapter thanks to the simulation of the spherical flame expanding in the torus network. It is clear that the flame, after having covered half of the obstacle modulus, has a similar behavior to the one obtained with the Gulder model. However, this simulation was not completed due to a local numerical problem. It has been chosen not to modify the numerical implementation which would have allowed us to finalize the calculation because it was already obvious that the flame would not capture a possible inflection point.

Eventually, the red curve represents the Zimont model with a coefficient C Z = 0.832, also determined in the previous chapter. As previously observed in figure 4.28 the turbulent flame speed obtained with this correlation is lower than that of Gulder. A modification of the coefficient could correct this, however the Zimont model shows a slight correction of the flame position from time 0.16 s which makes it coincide with the experimental curve. This last point is not specifically determinant but it is interesting to note that as for the Gulder model, the flame speed obtained numerically is very close to the experimental speed, which is encouraging.

In fact, the three correlation models provide very similar behaviours. This is quite natural since they come from the same family of models and use the same basic parameters. Indeed, let's recall that the turbulent flame speed is directly related to the velocity fluctuations by a relationship similar to the following equation [START_REF] Flohr | A turbulent flame speed closure model for LES of industrial burner flows[END_REF]:

S t S l
= 1 + f (Re, Da, P r ) , (5.1)

It seems that changing the proportion of one of these parameters has little impact on the behavior of the model flame. It would be necessary to consider models based on another paradigm. However, two points must be emphasized : (1) the flame speed obtained numerically with this formulation is very close to the final experimental speed and (2) no certainty about the origin of the inflection point can be stated but the presence of the polyethylene film remains one of the most likely causes. A flame speed in agreement with the experiment is obtained from the moment the film is lifted. Therefore, before anything else, a new measurement campaign or experimental comparison in similar situations should be considered, but with different obstruction rates to confirm our conclusions.

Pressure levels

After the position of the flame front over time, the amplitude of pressure variation at the flame passage is the second major parameter to capture when studying a large-scale deflagration. The temporal pressure signal have already observed in figure 5.9 at the three sensors located in the obstacle module and the spatial pressure signal in figure 5.13 when plotting different profiles from the numerical calculation. Finally, let us recall that in figure 5.8 the passage of the flame have been confirmed thanks to the experimental pressure signal.

Cléante Langrée

University of Rouen -CORIA February 8, 2022 A first very positive and encouraging observation is that the level of the pressure signals obtained from the numerical calculation is identical to that of the experiment. This is very good news because obtaining the amplitude of this variation is one of the key objectives of the simulation. On the other hand it is clear from the three figures that the two experimental and numerical pressure signals are out of phase in time. The offset is around 0.07 s for the first sensor, 0.05 s for the second and 0.03 s for the third sensor.

However, if one looks closely at the curves in figure 5.18, it appears that the experimental signal raises questions. Let's start with the first probe located at 0.28 cm inside the obstacle module (i.e., two inter-bar distances). The corresponding figure shows the two experimental and numerical signals as well as the position of the flame as seen on the experimental images. Since the experimental and numerical flames are located at the same place for this instant, both flame zones are superimposed on the hatched area. If the numerical pressure signal corresponds well to the passage of the corresponding flame with a peak in the flame zone, this is not the case for the experimental signal. The peak of the experimental signal is in fact located at time 0.1 s, i.e. well after the passage of the flame. It is difficult to interpret the origin of this shift, but one hypothesis can be put forward: this sensor is located in the

Cléante Langrée

University of Rouen -CORIA February 8, 2022

Conclusion 135

zone where the flame is still propagating in all directions in a quasi-spherical manner, including in the direction of the negative y's, in the mixture of fuel included in the thickness between the bars and the film. As a result, it is possible that the pressure of the burnt gases continues to increase after the flame has passed.

For the sensor placed at 0.98 m from the origin, the pressure signal can be seen in the figure 5.19 as the blue curve. In this case the blue band representing the position of the flame has a logical position. It is positioned a few fractions of a second before the pressure peak. However, let us not forget that the position of the flame front is not taken at ground level but at the center of the obstacle module. In addition, the position of the flame front is visually integrated along the depth of the module, whereas the pressure sensor measures a signal at a given position and thus, local fluctuations may occur.

Similarly, figure 5.20 compares the position of the flame front (blue band) with the pressure signal (blue curve) and a consistent result can be observed, again with a slight offset, this time with the flame downstream.

In both cases, the pressure signals corresponding to the numerical simulation were also plotted (orange curves) with the passage of the corresponding flames (orange bands). As already underlined, signals equivalent in amplitude are found. They are also very close from a temporal point of view. The difference is quite important and has two causes: (1) the shift of the average flame position already discussed above in figure 5.16 and (2) the fact that the experimental and numerical signal sensors capture a local instantaneous signal. Even without the flame shift, it is not possible for the two instantaneous signals to coincide perfectly. Statistical comparisons should be considered, but only one set of experimental measurements was available to us.

Conclusion

In this chapter, the comparison between a numerical calculation and an experimental configuration have been carried out directly, dedicated to the study of methane deflagration within an obstacle network.

First, several tests were performed to estimate the importance of some geometrical data on the main characteristics of the flame. It turns out that it was possible to simplify the way the bar crossings are taken into account as well as the experimental ignition box. This allowed us to focus on the main focus of the study: the interactions of the flame with the obstacle module.

In a first analysis of the experimental results, an inflection in the flame position curve have been observed, corresponding to two characteristic speeds of the flame: around 7 m/s on a first distance of 80 cm and around 25 m/s on the remaining 2,2 m of obstacle.

Industrial Application

It seemed particularly difficult to interpret the origin of this inflection so an AI algorithm has been used to reconstruct the experimental images with a clearer view of the flame. A similar curve, although less noisy, was found. From a numerical point of view, this inflection has not been observed in the case of the spherical flame propagating in the torus network and it was not observed in the EXJET case of this chapter either. Its origin could therefore be due to a runaway phenomenon that is not taken into account in the models or to a difference in initial data between the experiment and the calculation. Both are possible.

Another element raises questions : if the presence of obstacles generates a phenomenon that accelerates the flame (for example, the increase in the level of turbulence in the wakes), why is the acceleration not progressive ? On the observed curve, there is a clear swing at a given moment and then a constant speed. The hypothesis that this inflection of speed could be due to the initial presence of the polyethylene film has to be formulated. This film is placed on the metallic frames, completely covering the obstacle module initially and is lifted by the flame's blast. However, the moment when the film leaves the ground to be lifted completely is the one which corresponds to the inflection point. More experimental data are needed to support this hypothesis but it is clear that the simulation could not have taken into account the presence of this film.

Because of the inflection, it seems appropriate to discuss the time origin that has been selected for all the results presented. Indeed, the experimental zero time corresponds to the ignition time in the dedicated box as seen in figure 5.3. However, in order to concentrate on the aspects of the obstacle/flame interactions, the numerical ignition procedure has been simplified, as seen in figure 5.12. It was therefore necessary to choose a numerical time scale consistent with the experimental observations and the flame entrance in the obstacle module has been a base for this.

The figure 5.16 shows well, for small times, two concomitant numerical and experimental flames. However, due to the inflection of the experimental flame, a shift appears. Another solution would have been to select a time origin allowing to match the position of the flame at the output, but this would not have given a better insight.

In the absolute, the most important thing is to be able to find the flame speed and the pressure signal level. This was achieved in both cases. Indeed, the set of physical and numerical models which were selected allowed us to find two major information for the study of the explosions: the velocity of the flame and the amplitude of the pressure waves. This result is particularly encouraging for the development of numerical tools for the study of deflagrations and explosions.

Cléante Langrée

University of Rouen -CORIA February 8, 2022

Chap. 6 | Conclusion and perspectives

This study was set in the context in the general context of industrial risk assessment and prevention. In that context public actors like INERIS, HSE and also private actors like Total, DNV-GL and Air liquid invest every year to improve their skills and capability. A key skill in that regard, is the ability to precisely predict both the causes that can trigger hazardous events and the outcome of those events, should they ever occur. In order to be able to act more efficiently on both the prevention side and safety side, a better understanding of the prediction methods and tools is required. In this thesis, the focus has been put on one kind of hazardous event, particularly damaging, the Unconfined Vapour Cloud Explosion. Example of occurrences of such devastating events have been provided in the introduction section 1.

The more damaging effects of an UVCE is the overpressure effects in close and far range. Overpressure of few mbar can give a wide range of damaging effects to both structures and human, such as windows and wall destruction and irremediable wound to a human body. Thermal effects, whistle present, are not considered to be the main source of damage in this case. To evaluate those overpressure effects, experimental studies have been conducted for years. The INERIS participated in several experimental campaign in this perspective, DRA72, EXJET03 and AIRRE to name a few. These studies and those which came before, permitted to established the main parameters that will impact the overpressure effects consecutive to an UVCE, namely turbulence, presence of obstacle in the flammable mixture and mixing inhomogeneity. The present study proposed to study the first two, and focus on their impacts on our simulation of UVCE.

In the safety industry, both phenomenological tools, and CFD methods are presently used for this purpose. Historically, phenomenological tools have been the first technical solution for the need to predict UVCE consequences. However the adaptation of CFD methods from other fields (aeronautics, engine science) and the increase in CPU availability has made the use of CFD more and more relevant and reliable to handle this issue. FLACS (GexCon) have been pioneer in that regards, with solid experimental validation as early as the 1990s. It stays the most used tools for CFD simulation of UVCE, along another commercial code like ANSYS-Fluent (ANSYS). In the last decades, open sources code like OpenFoam started being used in UVCE numerical simulations. It present the main advantages of being easier to modulate and developed into, and benefiting from a large community of both users and developers. In that sense, using OpenFoam in this study served the purpose of improving the understanding of numerical methods and theoretical used for CFD simulation of UVCE. This improvement is expected to lead to a gain in confidence in the prediction capacity of CFD tools in the context of a risk analysis.

Conclusion and perspectives

The first step of this study was to review the vast turbulent combustion modeling topic. As presented in chapter 2, many approaches and models exist for modeling turbulent combustion, and many other for modeling turbulence itself. The focus of this review was to identify the features, strength and weaknesses of every approach, and relate them to our particular case of application. The UVCE being a really extreme case of turbulence combustion, at the opposite of the typical case of application of turbulent combustion models : the engine. As a matter of fact, an engine is a, relatively small scale combustion, happening in a closed vessel, where an UVCE is a large scale combustion happening in an open environment. The geometry, while potentially complex on both case, is also very different. In that regards a geometrical approach, with a sub grid wrinkling factor model has been used, for its theoretical relevance in every step of the UVCE process (from the quasi laminar spherical flame front to the fully perturbed flame front) and its capability to easily take into account mixing inhomogeneities. Turbulence has been modeled using RANS 2 equations models.

Those choices of modeling had been tested first in two configurations. First a simple 1D configuration to check that the typical characteristic of laminar combustion, like laminar flame speed and burnt gases temperature were recovered. Second a more complex 3D spherical configuration, on which interaction with turbulence and obstacle have been tested. As expected a strong impact of initial turbulence and obstruction on the flame propagation speed has been observed. Moreover, several turbulent flame speed correlation have been tested. Unsurprisingly, those already solidly validated correlation showed expected responses to both turbulence and obstacle, although different from one another. At the scale considered, both DNS and experiment presents difficulties to conclude on which model is more coherent. Finally, continuing towards a "real case" simulation, an experimental configuration of UVCE has been simulated. This experimental configuration has been set up and realized at this experimental facility of the INERIS. It consist in a deflagration across an obstacle module, with pressure probes across the field, both inside and outside the obstacle module. The main goal of this simulation was to compare experimental and numerical flame front propagation and pressure fields. Where flame front propagation and induced pressure fields were globally well reproduced, some aspects of the experimental UVCE were not captured. First, a sudden increase of flame propagation speed is experimentally observed, and not recovered in numerical results. A few hypothesis to explain this, not self evident, increase of flame propagation speed have been formulated. Other numerical work on experimental configuration were some of the potential experimental bias of the EXJET configuration can be diminished could provide definitive answer on that matter. Second, pressure signal, while conserving the order of magnitude, could improve its fitness to experiment in shape.

Nevertheless, it seems to us that those results are encouraging for the pursuit of work in numerical simulation of UVCE with OpenFoam. While already showing good results, the open source nature of the code make it easier to investigate and implement "on demand" improvements. Such leads for improvement can already be proposed. First, while the turbulent combustion model was the real focus of the study, further work on turbulence modeling could provide more precise capture of the phenomenon at play. For instance, it's expected that the fresh gases flow induced by the flame propagation should generate wake turbulence when interacting with obstacle upstream. Looking at figure 6.1 an increase in TKE before the passage of the flame front is noted. But it's : 1) really late, meaning that the flame front has to be relatively close to induce the effect and 2) negligeable in comparison to the increase of TKE done by the passage of the flame itself. Given eq. 2.67, TKE value impacts the

Cléante Langrée

University of Rouen -CORIA February 8, 2022 most the flame propagation within the flame front (i.e. when ∇ (b) � = 0). One hypothesis would be that the deflagration is too slow, and the obstacle not important enough to induce significant wake turbulence. Simulation of a more powerful UVCE could still trigger this kind of effects with the model as it is used in this thesis. In the latest month of this PhD thesis, some turbulence models with integration of a turbulence production terms has been explored, but with no conclusive results to show for it. Nevertheless, this work could be continued to conclude on the relevance of such a modeling choice. Moreover, while only methane/air mixtures have been considered in this study, the enlargement of the study to other fuel could prove to be insightful for the general simulation of UVCE with OpenFoam. While also not conclusive enough to be shown, some simulation of hydrogen/air have been realized during this thesis. In addition to being a trendy topic, hydrogen/air flame, due to the extreme behavior of hydrogen in comparison to other hydrocarbon fuels, could provide extreme situation in which to test the modeling choices made in this thesis

3. 3 5 . 8

 358 Process of the snappyHexMesh software. Except near the obstacle, the mesh obtained will be mostly hexahedral. A last optional stage is not represented : the addition of boundary layer cells around the obstacle. . . . . . . . . . . . . .3.4 Propagating flame (white) around obstacles (purple). The mesh structure is apparent on the flame front. Mesh has been generated with the module netgen (3D non-structured) of the SALOME software. . . . . . . . . . . . . . 3.5 Schematic representation of face interpolation. . . . . . . . . . . . . . . . . . . 4.1 Schematic representation of the uni dimensional configuration. . . . . . . . . . . 4.2 Comparison of the flame ignition modes. The ignite module (a and b) generates non-physical velocities and temperatures. The burnt gas deposit (c and d) shows correct velocity and temperature profiles. . . . . . . . . . . . . . . 4.3 Flame velocity vs equivalence ratio. Comparison of the XiFoam solver using a progress variable transport equation with Cantera using the GRIMech 3.0 mechanism. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Adiabatic temperature vs equivalence ratio. Comparison of the XiFoam solver using a progress variable transport equation with Cantera using the GRIMech 3.0 mechanism. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Impact of the air mixture definition (molar proportion of oxygen vs nitrogen) on the velocity bell shape curve. Obtained with Cantera. . . . . . . . . . . . . . 4.6 Grid Convergence : impact of the grid step on the premixed flame velocity. . . . 4.7 Impact of the waveTransmissive boundary condition on the pressure field. . . 4.8 Different arrangements of obstacles. Staggered : the flame always meets an obstacle or aligned : the flame spreads freely between a series of aligned obstacles.88 4.9 Tubulent flame surface evolution in the case without obstacle (case free-01_00), time evolution from top to bottom and left to right, from 0.04 s to 0.14 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Average temperature (top) and velocity norm (bottom) profiles over timecase free-01_00. The colored areas represent the standard deviation of the corresponding color curve. The vertical red bars indicate the flame position. . . . 4.11 Flame radius evolution depending on the selected method to compute flame front position (r f : burnt/fresh gas equilibrium, r e : volume equivalent sphere, r o : outer radius) -case free-01_00. . . . . . . . . . . . . . . . . . . . . . . . 4.12 Flame radius evolution depending on the mesh refinement, from 700k cells up to 7M cells -case free-01_00. . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Turbulent burning velocity S t evolution and bending apparition in the obstacle free environment -case free-01_00. . . . . . . . . . . . . . . . . . . . . . . . 4.14 Flame radius evolution depending on the mesh refinement and initial turbulence level.96 4.15 Flame propagation, staggered obstacles (2 and 6 cm radius), equivalent flame radius : 2 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Azimuthal and elevation cuts of the sta6-0400 flame front, equivalent flame radius : 2 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Azimuthal and elevation cuts of the lin6-0500 flame front. The elevation cut plan is located between two ranks of obstacles, equivalent flame radius : 2 m. 4.18 Flame surface evolution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.19 Impact of the various obstacles on the flame expansion. Grey lines represent the radius of the first (bottom line) and last (top line) obstacles encountered by the flame during its propagation. . . . . . . . . . . . . . . . . . . . . . . . . 100 4.20 Time evolution of the PDF of the flame radius. . . . . . . . . . . . . . . . . . . 101 Time evolution of the PDF angle between the flame normal vector and the sphere normal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Time evolution of the PDF of the flame radius, sta6-04_00 case. . . . . . . . 4.23 Comparison of the PDF of the flame radius : staggered (sta6-04_00) vs in line obstacles (lin6-05_00) case. . . . . . . . . . . . . . . . . . . . . . . . . . 4.24 Comparison of the PDF of the flame surface orientation : staggered (sta6-04_00) vs in line obstacles lin6-05_00 case. . . . . . . . . . . . . . . . . . . . 4.25 Comparison of the PDF P (r ) flame position : free (free-01_00) vs staggered 2 cm radius obstacles (sta2-02_00) vs staggered 6 cm radius obstacles (lin6-05_00). All the flames have the same equivalent radius : 2 m. . . . . . . 4.26 Azimuthal and elevation cuts of the sta2-02_00 flame front, t = 0.095 s. . . . 4.27 Comparison of the PDF P (β) of the flame surface orientation : free (free-01_00) vs staggered 2 cm radius obstacles (sta2-02_00) vs staggered 6 cm radius obstacles (lin6-05_00). All the flames have the same equivalent radius : 2 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.28 Evaluation of various turbulent combustion correlation models. . . . . . . . . . . 4.29 Comparison of the PDF P (r ) of the flame position using three different combustion correlation models. All the flames have the same equivalent radius : 2 m.110 4.30 Comparison of the PDF P (β) of the flame surface orientation using three different combustion correlation models. All the flames have the same equivalent radius : 2 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Different scales of turbulent combustion simulation. From left to right : Numerical temperature profile for 1D laminar flames with Cantera and Open-Foam ; Global view of a numerical flame front propagating in the 3D spherical configuration ; Numerical flame front propagating in the numerical simulation of the EXJET experimental configuration ; Aerial view of the Buncefield industrial site after the Buncefield disaster. . . . . . . . . . . . . . . . . . . . . 5.2 Locations of the release orifice, pipes, igniter and pressure sensors used by INERIS and GDF SUEZ for methane jets within a 3D array of 20 mm diameter tubes. ([18]). Up : top view of the obstacle module ; Bottom : side view of the obstacle module. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Photograph of the "left" side of the experimental configuration : the steel frame support of the polyethylene sheet, the ignition box and the tubes. . . . . . 5.4 Example of the processing of the raw image provided by the high speed camera by the Photoshop software and an algorithm based on artificial intelligence allowing the recognition of the different elements of the image. . . . . . . . . . 5.5 Instant capture by high frequency camera of the early stages of the propagating flame front. Instant depicted are 0.008 s (top) 0.02 s (middle) and 0.06 s (bottom). Flame kernel is identified with plain red circle, area where the flame front is identified but barely visible is identified with dashed red circle. . . . . . . 5.6 Instant capture by high frequency camera of the propagating flame front. Instant depicted are 0.07 s (top) 0.09 s (middle) and 0.11 s (bottom). Polyethylene film border is identified with dashed green line. . . . . . . . . . . . 5.7 Instant capture by high frequency camera of the early stages of the propagating flame front. Instant depicted are 0.12 s (top) 0.16 s (middle) and 0.2 s (bottom). Flame front is identified with dashed red line and polyethylene sheet border is identified with dashed green line. . . . . . . . . . . . . . . . . . . . . . Experimental evolution of the position of the flame front after ignition. Yellow triangles: first extraction from the raw images, blue square: extraction from the AI colorized images, orange lines : pressure signal from the flame (ground sensors). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Experimental temporal pressure signal at 3 different probe locations inside the obstacle module: y = 0.28 cm (blue), y = 0.98 cm (orange), y = 1.96 cm (yellow). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Arrangement of the crossing between the obstacles. The simplified layout allows an economy in the number of meshes. . . . . . . . . . . . . . . . . . . . 5.11 Paraview screenshot of the two mesh configurations. . . . . . . . . . . . . . . . 5.12 Simplification of the ignition process. . . . . . . . . . . . . . . . . . . . . . . . . 5.13 Temperature and pressure spatial profiles over the main y axis of the obstacle module. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.14 Comparison of numerical (bottom) and experimental (top) flame front, experimental time 0.03 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.15 Comparison of numerical (bottom) and experimental (top) flame front, experimental time 0.16 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Evolution of numerical and experimental flame front position. The fuchsia lines have similar slopes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.17 Evolution of numerical flame front position for 3 different turbulent flame speed correlation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.18 Numerical (orange) and experimental (blue) pressure signal at 0.28 meters from the ignition side of the obstacle module. . . . . . . . . . . . . . . . . . . . 5.19 Numerical (orange) and experimental (blue) pressure signal at 0.98 meters from the ignition side of the obstacle module. . . . . . . . . . . . . . . . . . . . 5.20 Numerical (orange) and experimental (blue) pressure signal at 1.96 meters from the ignition side of the obstacle module. . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 -

 11 Figure 1.1 -Example of a scoring matrix.

Figure 1 . 2 -

 12 Figure 1.2 -Example of a pressure wave produced by a deflagration. From [1].

Figure 1 . 4 -

 14 Figure 1.4 -Aerial view of the devastation after the Port Hudson accident.

1.2 UVCE 7 Figure 1 . 5 -

 715 Figure 1.5 -Aerial view of the devastation after the Flixborough accident.

Figure 1 . 6 -

 16 Figure 1.6 -Aerial view of the devastation after the Buncefield accident.

Figure 1 . 7 -

 17 Figure 1.7 -Overpressure against radius abacus from TNT equivalent method [3].

Figure 2 . 1 -

 21 Figure 2.1 -Fire triangle, or combustion triangle.

Figure 2 . 2 -

 22 Figure 2.2 -Structure of a laminar premixed flame, from [5].

Fig. 2 .

 2 3 illustrate new definitions of "flame speeds".

Figure 2 . 3 -

 23 Figure 2.3 -Flame speed definitions.

  .

  4.

Figure 2 . 4 -

 24 Figure 2.4 -Flame thickness definitions. From [6].

2. 1 description 21 Figure 2 . 5 -

 12125 Figure 2.5 -Flame curvature definitions. From [7].

Figure 2 . 7 -

 27 Figure 2.7 -Gravity induced instabilities : stable configuration (left) and unstable configuration (right). From [8].

Figure 2 . 8 -

 28 Figure 2.8 -Darrieus Landau instabilities. From [9].

Figure 2 .

 2 Figure 2.10 -Kolmogorov's turbulent energy cascade.

Figure 2 . 11 -

 211 Figure 2.11 -Borghi diagram, modified by Proust [11].

Figure 2 . 13 -

 213 Figure 2.13 -Schematic representation of leak scenarios.

  (a) Dispersion of the flammable cloud. (b) Ignition. (c) Hemispherical flame propagation. (d) Transition between hemispherical and azimutal propagation. (e) Azimutal propagation.

Figure 2 . 14 -

 214 Figure 2.14 -Stages of an UVCE from [12].

Figure 2 . 15 -

 215 Figure 2.15 -Turbulence intensity effects on the flame velocity. From [13].

Figure 2 . 16 -

 216 Figure 2.16 -Overpressure wave in a detonation scenario (a) and a deflagration scenario (b).From[START_REF] Proust | Formalisation du savoir et des outils dans le domaine des risques majeurs (eat-dra-76) les explosions non confinées de gaz et de vapeursω uvce[END_REF].

Figure 2 . 19 -

 219 Figure 2.19 -Illustration of PDF modeling for turbulent combustion, from [15].

2. 2 45 Figure 2 . 20 -

 245220 Figure 2.20 -Illustration of a flame surface crossing a control volume[START_REF] Ahmed | On the stress-strain alignment in premixed turbulent flames[END_REF]

64 )Figure 2 .

 642 Figure 2.21 -wrinkling concept representation. Dashed line : flame front, as captured by the mesh.

Figure 3 . 1 -

 31 Figure 3.1 -Schematic representation of a control volume. From [17].

3. 3

 3 Finite Volume Method57 (a) Background mesh. (b) Geometry. (c) Castellated mesh. (d) Snapped mesh.

Figure 3 . 3 -

 33 Figure 3.3 -Process of the snappyHexMesh software. Except near the obstacle, the mesh obtained will be mostly hexahedral. A last optional stage is not represented : the addition of boundary layer cells around the obstacle.

Figure 3 . 4 -

 34 Figure 3.4 -Propagating flame (white) around obstacles (purple). The mesh structure is apparent on the flame front. Mesh has been generated with the module netgen (3D non-structured) of the SALOME software.

3. 3 Finite Volume Method 61 Figure 3 . 5 -

 36135 Figure 3.5 -Schematic representation of face interpolation.

3. 4

 4 Pressure-velocity coupling 63 waveTransmissive

2 . 58 )

 258 , with b = 1 -c and c, the progress variable. All those equations are solved following the algorithm shown in (fig. 3.5.1). The equation for b and Xi together with the energy equation are solved at the same level as the momentum equation in the Pimple, Piso or Simple loop. This aims to enforce the coupling between the reaction rate and the flow dynamics.

  3.3) and TensorI.H (lst. 3.4), both OpenFoam source files. Listing 3.3 -divDevRhoReff in OpenFoam -ReynoldsStress.C fvc::laplacian ( this ->alpha_*rho*this ->nut(), U, "laplacian(nuEff ,U)" ) + fvc::div(this ->alpha_*rho*R_) -fvc::div(this ->alpha_*rho*this ->nu()*dev2(T(fvc::grad(U)))) -fvm::laplacian(this ->alpha_*rho*this ->nuEff(), U) Listing 3.4 -dev2 in OpenFoam -TensorI.H inline Tensor <Cmpt> dev2(const Tensor <Cmpt >& t) { return t -SphericalTensor <Cmpt >::twoThirdsI *tr(t); } And then it's implemented in XiFoam in UEqn.H (lst. 3.5) : tUEqn.ref(); UEqn.relax(); fvOptions.constrain(UEqn); if (pimple.momentumPredictor()) { solve(UEqn == -fvc::grad(p)); fvOptions.correct(U); K = 0.5*magSqr(U); }

3. 5

 5 XiFoam

40 )Listing 3 . 6 -

 4036 Φ st calculation -bEqn.H surfaceScalarField phiSt("phiSt", fvc::interpolate(rhou*StCorr*Su*Xi)*nf);

3 . 41 )

 341 Listing 3.7n f calculation -bEqn.H volVectorField n("n", fvc::grad(b)); volScalarField mgb(mag(n)); dimensionedScalar dMgb = 1.0e-3* (b*c*mgb)().weightedAverage(mesh.V()) /((b*c)().weightedAverage(mesh.V()) + small) + dimensionedScalar("ddMgb", mgb.dimensions(), small); mgb += dMgb; surfaceVectorField SfHat(mesh.Sf()/mesh.magSf()); surfaceVectorField nfVec(fvc::interpolate(n)); nfVec += SfHat*(fvc::snGrad(b) -(SfHat & nfVec)); nfVec /= (mag(nfVec) + dMgb); surfaceScalarField nf((mesh.Sf() & nfVec)); n /= mgb;

48 )

 48 It is implemented in XiFoam in bEqn.H (lst. 3

inlineListing 3 . 10 -

 310 Foam::scalar Foam::laminarFlameSpeedModels::Gulders::SuRef ( scalar phi ) const { if (phi > small) { return W_*pow(phi, eta_)*exp(-xi_*sqr(phi -1.075)); Flame speed correlation -Pressure/Temperature -Gulder.C inline Foam::scalar Foam::laminarFlameSpeedModels::Gulders:static const scalar Tref = 300.0; static const scalar pRef = 1.013e5; return SuRef(phi)*pow((Tu/Tref), alpha_)*pow((p/pRef), beta_)*(1 -f_*Yres); } Numerical tools for UVCE simulation 3.5.5 Model for mixture inhomogeneity

51 )

 51 Listing 3.11 -total fuel mass fraction equation -ftEqn.H if (composition.contains("ft")) { volScalarField& ft = composition.Y("ft"

3. 5 solver 75 Listing 3 . 12 -

 575312 XiFoam XiFoam composition calculation -inhomogeneousMixture.C scalar fu = b*ft + (1.0 -b)*fres(ft, stoicRatio().value()); scalar ox = 1 -ft -(ft -fu)*stoicRatio().value(); scalar pr = 1 -fu -ox;

13 -

 13 Residual burnt fuel mass fraction left calculation -basicCombustionMix-tureI.H inline Foam::tmp<Foam::volScalarField > Foam::basicCombustionMixture::fres ( const volScalarField& ft, const dimensionedScalar& stoicRatio ) const { return max(ft -(scalar(1) -ft)/stoicRatio.value(), scalar(0)); }

Figure 4 . 1 -

 41 Figure 4.1 -Schematic representation of the uni dimensional configuration.

4. 2

 2 Laminar premixed flames 81 (a)-Initial velocity, ignite module. (b) -Initial temperature, ignite module. (c) -Initial velocity, burnt gases core deposit. (d) -Initial temperature, burnt gases core deposit.

Figure 4 . 2 -

 42 Figure 4.2 -Comparison of the flame ignition modes. The ignite module (a and b) generates non-physical velocities and temperatures. The burnt gas deposit (c and d) shows correct velocity and temperature profiles.

4. 2 premixed flames 83 Figure 4 . 3 -

 28343 Figure 4.3 -Flame velocity vs equivalence ratio. Comparison of the XiFoam solver using a progress variable transport equation with Cantera using the GRIMech 3.0 mechanism.

Figure 4 . 4 -

 44 Figure 4.4 -Adiabatic temperature vs equivalence ratio. Comparison of the XiFoam solver using a progress variable transport equation with Cantera using the GRIMech 3.0 mechanism.

. 5 )

 5 visible in fig. 4.3, triangle symbols.

4. 2 premixed flames 85 Figure 4 . 5 -

 28545 Figure 4.5 -Impact of the air mixture definition (molar proportion of oxygen vs nitrogen) on the velocity bell shape curve. Obtained with Cantera.

Figure 4 . 6 -

 46 Figure 4.6 -Grid Convergence : impact of the grid step on the premixed flame velocity.

Figure 4 . 7 -

 47 Figure 4.7 -Impact of the waveTransmissive boundary condition on the pressure field.

Figure 4 . 8 -

 48 Figure 4.8 -Different arrangements of obstacles. Staggered : the flame always meets an obstacle or aligned : the flame spreads freely between a series of aligned obstacles.

  .

Table 4 . 1 -Flames

 41 Numerical configurations. Combination of the tested relation for the turbulent velocity with the obstacles modes : no obstacle (free), staggered obstacles (Fig. 4.8-(a)) with 2 or 6 cm radius tori (sta2 and sta6), aligned obstacles (Fig. 4.8-(b)) with 2 and 6 cm radius tori (lin2 and lin6).

4. 3 expansion 91 Figure 4 . 9 -

 39149 Figure 4.9 -Tubulent flame surface evolution in the case without obstacle (case free-01_00), time evolution from top to bottom and left to right, from 0.04 s to 0.14 s.

Figure 4 . 10 -

 410 Figure 4.10 -Average temperature (top) and velocity norm (bottom) profiles over timecase free-01_00. The colored areas represent the standard deviation of the corresponding color curve. The vertical red bars indicate the flame position.

Figure 4 . 13 -

 413 Figure 4.13 -Turbulent burning velocity S t evolution and bending apparition in the obstacle free environment -case free-01_00.

Figure 4 . 14 -

 414 Figure 4.14 -Flame radius evolution depending on the mesh refinement and initial turbulence level.

Figure 4 . 18 -

 418 Figure 4.18 -Flame surface evolution.

Figure 4 . 19 -

 419 Figure 4.19 -Impact of the various obstacles on the flame expansion. Grey lines represent the radius of the first (bottom line) and last (top line) obstacles encountered by the flame during its propagation.

  Free flame structure First, The probability density function (PDF) of the distance of the flame surface (as seen in fig. 4.9) to the center of the spherical domain can be plotted. The average value of this radius is very close to the value of the equivalent radius used so far. Fig. 4.20 represents the temporal evolution of the position of the free flame over time. It appears that the PDF has the form of a unimodal, Gaussian like, distribution. This PDF is very steep due to the low spread of the flame right after ignition. Then it widens with time due to the appearance and growth of combustion instabilities (see fig.4.9).

Figure 4 . 20 -

 420 Figure 4.20 -Time evolution of the PDF of the flame radius.

Figure 4 . 21 -

 421 Figure 4.21 -Time evolution of the PDF angle between the flame normal vector and the sphere normal.

4. 3 expansion 103 Figure 4 . 22 -

 3103422 Figure 4.22 -Time evolution of the PDF of the flame radius, sta6-04_00 case.

Figure 4 . 24 -

 424 Figure 4.24 -Comparison of the PDF of the flame surface orientation : staggered (sta6-04_00) vs in line obstacles lin6-05_00 case.

Figure 4 . 25 -

 425 Figure 4.25 -Comparison of the PDF P (r ) flame position : free (free-01_00) vs staggered 2 cm radius obstacles (sta2-02_00) vs staggered 6 cm radius obstacles (lin6-05_00). All the flames have the same equivalent radius : 2 m.

4. 3

 3 Flame expansion 107 obstacles increase the flame surface, they tend to decrease the development of instabilities by their interactions with the flame. By observing the elevation section of fig.4.26 -right, multiple 'burgeoning' structures can be noticed. They are much smaller in the case with obstacles 9 times larger (Fig.4.16 -right). However, even if the obstacles slow down the development of instabilities, the modification of the flame structure that they induce is largely prevalent, hence the strong increase in the propagation speed.

Figure 4 . 26 -

 426 Figure 4.26 -Azimuthal and elevation cuts of the sta2-02_00 flame front, t = 0.095 s.

Figure 4 . 27 -

 427 Figure 4.27 -Comparison of the PDF P (β) of the flame surface orientation : free (free-01_00) vs staggered 2 cm radius obstacles (sta2-02_00) vs staggered 6 cm radius obstacles (lin6-05_00). All the flames have the same equivalent radius : 2 m.

4. 3

 3 Flame expansion 109 same. This shows that the three models, based on the same assumptions, yield similar results as soon as the model parameters are correctly set. In these simulations, the default Gulder parameter was used (C g = 0.62) and the Zimont and FLACS parameters are C z = 0.96 and C f = 1.536, respectively.

Figure 4 . 28 -

 428 Figure 4.28 -Evaluation of various turbulent combustion correlation models.

Figure 4 . 29 -

 429 Figure 4.29 -Comparison of the PDF P (r ) of the flame position using three different combustion correlation models. All the flames have the same equivalent radius : 2 m.

4. 4 111 Figure 4 . 30 -

 4111430 Figure 4.30 -Comparison of the PDF P (β) of the flame surface orientation using three different combustion correlation models. All the flames have the same equivalent radius : 2 m.

Figure 5 . 1 -

 51 Figure 5.1 -Different scales of turbulent combustion simulation. From left to right : Numerical temperature profile for 1D laminar flames with Cantera and OpenFoam ; Global view of a numerical flame front propagating in the 3D spherical configuration ; Numerical flame front propagating in the numerical simulation of the EXJET experimental configuration ; Aerial view of the Buncefield industrial site after the Buncefield disaster.

Figure 5 . 2 - 115 Figure 5 . 3 -

 5211553 Figure 5.2 -Locations of the release orifice, pipes, igniter and pressure sensors used by INERIS and GDF SUEZ for methane jets within a 3D array of 20 mm diameter tubes. ([18]). Up : top view of the obstacle module ; Bottom : side view of the obstacle module.

Figure 5 . 2 (Figure 5 . 4 -

 5254 Figure 5.4 -Example of the processing of the raw image provided by the high speed camera by the Photoshop software and an algorithm based on artificial intelligence allowing the recognition of the different elements of the image.

Figure 5 . 5 -Figure 5 . 6 -

 5556 Figure 5.5 -Instant capture by high frequency camera of the early stages of the propagating flame front. Instant depicted are 0.008 s (top) 0.02 s (middle) and 0.06 s (bottom). Flame kernel is identified with plain red circle, area where the flame front is identified but barely visible is identified with dashed red circle.

5. 3 Figure 5 . 7 -Figure 5 . 8 -

 35758 Figure 5.7 -Instant capture by high frequency camera of the early stages of the propagating flame front. Instant depicted are 0.12 s (top) 0.16 s (middle) and 0.2 s (bottom). Flame front is identified with dashed red line and polyethylene sheet border is identified with dashed green line.

Figure 5 . 9 -

 59 Figure 5.9 -Experimental temporal pressure signal at 3 different probe locations inside the obstacle module: y = 0.28 cm (blue), y = 0.98 cm (orange), y = 1.96 cm (yellow).

Figure 5 .figure 5 . 9 .

 559 Figure 5.8 shows the moments when the passage of the flame is measured at the three sensors. The orange crosses are a rough estimation of the flame front position (peak of the pressure signal) while the horizontal orange lines mark the moments when the pressure increases slightly in front of the flame to this peak. The pressure signals can be seen in

Figure 5 . 10 -

 510 Figure 5.10 -Arrangement of the crossing between the obstacles. The simplified layout allows an economy in the number of meshes.

1 -Figure 5 . 11 -

 1511 Figure 5.11 -Paraview screenshot of the two mesh configurations.

Figure 5 .

 5 Figure 5.15 below shows the same comparison at time 0.16 s. From the simulation point of view (bottom), it is possible to see the flame front disturbed by the bars in the obstacle network. Its position is very close to that of the experimental image (top), this is an important first positive observation.A second encouraging point is the general shape of the flame above the obstacle. Indeed, the pear shape that can be observed numerically (bottom) is also found at the level of the polyethylene film (top). This film does not delimit the flame front, but the moving fluid that pushes the film away is directly the result of the advancing front. It is therefore quite possible that the shape of the film is a replica of the shape of the front at a given distance and the flame front shape obtained with the simulation is quite similar.

Figure 5 . 13 -

 513 Figure 5.13 -Temperature and pressure spatial profiles over the main y axis of the obstacle module.

5. 4 Figure 5 . 14 -

 4514 Figure 5.14 -Comparison of numerical (bottom) and experimental (top) flame front, experimental time 0.03 s.

Figure 5 . 15 -Figure 5 . 16 -

 515516 Figure 5.15 -Comparison of numerical (bottom) and experimental (top) flame front, experimental time 0.16 s.

Figure 5 . 17 -

 517 Figure 5.17 -Evolution of numerical flame front position for 3 different turbulent flame speed correlation.

5. 4

 4 Numerical Simulation 133 A quantitative comparison of the pressure signals from the simulation and the experiment is performed for the three sensors included in the obstacle module. Figures 5.18

  , 5.19 and 5.20 show these signals for probes place on the y axis at 0.28 cm, 0.98 cm and 1.96 cm from the first tube of the module (Fig.5.2 -L3, L4 and L5).

Figure 5 .

 5 Figure 5.18 -Numerical (orange) and experimental (blue) pressure signal at 0.28 meters from the ignition side of the obstacle module.

Figure 5 .

 5 Figure 5.19 -Numerical (orange) and experimental (blue) pressure signal at 0.98 meters from the ignition side of the obstacle module.

Figure 5 .

 5 Figure 5.20 -Numerical (orange) and experimental (blue) pressure signal at 1.96 meters from the ignition side of the obstacle module.

Figure 6 . 1 -

 61 Figure 6.1 -Turbulent kinetic energy profile along the obstacle module axis.

Table of Contents Page 1

 ofPage Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Domain discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3.2 Equations discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.3.3 Surface and Volume integrals approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.3.4 Temporal scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.3.5 Spatial schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

	x
	1.1	Context of this PhD		1
	1.2	UVCE		2
		Introduction		51
	3.2	OpenFoam : A brief description	52
	3.3	Finite Volume Method		53
	3.3.1 3.4	Pressure-velocity coupling		63
	Cléante Langrée	UVCE Simulation	February 8, 2022

3.3.6 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.4.1 Semi Implicit Method for Pressure Linked Equations -SIMPLE . . . . . . . . . . . . . . . . . 65 3.4.2 Pressure Implicit with Splitting Operator -PISO . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

TABLE OF CONTENTS

 OF 3.4.3 PIMPLE : PIso + sIMPLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Example of a scoring matrix. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Example of a pressure wave produced by a deflagration. From [1]. . . . . . . . . 4 1.3 Characteristics of a shock wave produced by a detonation. From [2]. . . . . . . . 5 1.4 Aerial view of the devastation after the Port Hudson accident. . . . . . . . . . . 6 1.5 Aerial view of the devastation after the Flixborough accident. . . . . . . . . . . . 7 1.6 Aerial view of the devastation after the Buncefield accident. . . Darrieus Landau instabilities. From [9]. . . . . . . . . . . . . . . . . . . . . . . . 2.9 Thermo diffusive instabilities. From [10]. . . . . . . . . . . . . . . . . . . . . . . 2.10 Kolmogorov's turbulent energy cascade. . . . . . . . . . . . . . . . . . . . . . . 2.11 Borghi diagram, modified by Proust [11]. . . . . . . . . . . . . . . . . . . . . . 2.12 Flammability limits. Midbar : stoichiometry. . . . . . . . . . . . . . . . . . . . . 2.13 Schematic representation of leak scenarios. . . . . . . . . . . . . . . . . . . . . 2.14 Stages of an UVCE from [12]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.15 Turbulence intensity effects on the flame velocity. From [13]. . . . . . . . . . . . 2.16 Overpressure wave in a detonation scenario (a) and a deflagration scenario (b). From [12]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.17 Turbulent spectrum representation of each turbulence modeling method. . . . . 2.18 Classification of turbulence modeling methods. Adapted from [14]. . . . . . . . . 2.19 Illustration of PDF modeling for turbulent combustion, from [15]. . . . . . . . . 2.20 Illustration of a flame surface crossing a control volume [16] . . . . . . . . . . . 2.21 wrinkling concept representation. Dashed line : flame front, as captured by the mesh.

	xii			LIST OF FIGURES
	3.5	XiFoam solver		66
	1.1 Cléante Langrée	UVCE Simulation	February 8, 2022

. . . . . . . . . 7 1.7 Overpressure against radius abacus from TNT equivalent method

[START_REF]U.s. department of the army, s"structure to resist the effects of accidental explosions[END_REF]

. . . . . . . 9 1.8 Abacus related to the Multi-Energy method giving the overpressures generated by deflagrations at constant flame speed for hemispherical charges located on the ground

[START_REF] Van Den | The multi-energy method -a framework for vapour cloud explosion blast prediction[END_REF]

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Fire triangle, or combustion triangle. . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Structure of a laminar premixed flame, from [5]. . . . . . . . . . . . . . . . . . . 2.3 Flame speed definitions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Flame thickness definitions. From [6]. . . . . . . . . . . . . . . . . . . . . . . . 2.5 Flame curvature definitions. From [7]. . . . . . . . . . . . . . . . . . . . . . . . 2.6 Curvature effects on flame front. From [6]. . . . . . . . . . . . . . . . . . . . . 2.7 Gravity induced instabilities : stable configuration (left) and unstable configuration (right). From [8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Schematic representation of a control volume. From [17]. . . . . . . . . . . . . 3.2 Mesh generated by blockMesh. . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 1 .

 1 [START_REF] Trelat | Impact de fortes explosions sur les batiments representatifs d'une installation industrielle[END_REF] gives and example of a qualitative scale for the intensity.

	Intensity	Qualitative scale
	I1	Intensity of the dangerous phenomenon limited to a few meters
	I2	Intensity of the dangerous phenomenon limited to 20 m
	I3	Intensity of the dangerous phenomenon limited to 50 m
	I4	Intensity of the dangerous phenomenon out of the plant perimeter

Table 1 .

 1

1 -Example of a qualitative scale for the intensity.

Table

  .2 at atmospheric conditions.

	Flammable product	LFL (in vol. with air) UFL (in vol. with air)
	Methane	5 %	15 %
	Propane	2.2 %	10 %
	Hydrogen	4 %	75 %
	Gasoline (octane index about 55)	1.4 %	7.4 %
	Ethylene oxide	3 %	100 %

  As shown in Fig 2.14 an UVCE scenario can be divided into three elementary steps :

	2.1 Physical description					29
	Step by step description
	Flammable cloud formation
	H 2
	CH 4
	C 3 H 8 C 2 H 6 "Fuel"
	C 4 H 10
	C 6 H 12
		20	40	60	80	100
		"Fuel volume fraction"
	Non-flammable mixture	Flammable mixture
	Cléante Langrée	University of Rouen -CORIA	February 8, 2022

Figure 2.12 -Flammability limits. Midbar : stoichiometry.

  • Use the pressure from the previous step to compute a new velocity field with equation (3.23) • Repeat from the momentum predictor step until the pressure and velocity field both satisfy the continuity (2.19) and momentum equation (2.20).

•

  Solve the Poisson equation (3.27) using the predicted velocity field . • Use the pressure from the previous step to compute a new velocity field with equation (3.23) • Repeat from the Poisson equation for pressure, with the new velocity for last step,

until the pressure and velocity field both satisfy the continuity (2.19) and momentum equation

(2.20)

.

  As the PIMPLE algorithm is used, a Poisson pressure equation (eq. 3.27) is derived from the momentum (eq. 2.41) and continuity equation (eq. 2.40).

	76	Numerical tools for UVCE simulation
	fvScalarMatrix heauEqn
	(
	fvm::ddt(rho, heau) + mvConvection ->fvmDiv(phi, heau)
	+ (fvc::ddt(rho, K) + fvc::div(phi, K))*rho/thermo.rhou()
	+ (
	heau.name() == "eau"
	? fvc::div
	(
	fvc::absolute(phi/fvc::interpolate(rho), U),
	p,
	"div(phiv,p)"
	)*rho/thermo.rhou()
	: -dpdt*rho/thermo.rhou()
	)
	-fvm::laplacian(turbulence ->alphaEff(), heau)
	==
	fvOptions(rho, heau)
	);
	Listing 3.15 -Total energy equation -EaEqn.H
	fvScalarMatrix EaEqn
	(
	p,
	"div(phiv,p)"
	)
	: -dpdt
	)
	-fvm::laplacian(turbulence ->alphaEff(), hea)
	+ fvOptions(rho, hea)
	);
	3.5.7 Pressure equation
	Listing 3.16 -Pressure equation -pEqn.H
	fvScalarMatrix pEqn
	(
	fvm::ddt(psi, p)
	+ fvc::div(phiHbyA)
	-fvm::laplacian(rhorAUf , p)
	==
	fvOptions(psi, p, rho.name())
	);
	Cléante Langrée	UVCE Simulation	February 8, 2022

[START_REF] Williams | Combustion Theory[END_REF]

Listing 3.14 -Unburnt gases energy equation -EauEqn.H fvm::ddt(rho, hea) + mvConvection ->fvmDiv(phi, hea) + fvc::ddt(rho, K) + fvc::div(phi, K) + ( hea.name() == "ea" ? fvc::div ( fvc::absolute(phi/fvc::interpolate(rho), U),

Cléante LangréeUniversity of Rouen -CORIA February 8, 20222.1 Physical description

Remerciements

List of Figures

As it has been observed in fig. 4.16 and fig. 4.17, the flame structure in the staggered and aligned cases is different. As detailed above, the general shape of the probability of presence of the flame front P (r ) is similar in both cases due to the presence of the obstacles, which all have the same size and obstruction rate. However, variations in the visible probability of presence can be observed because of the staggered arrangement of every other row. Consequently, it is now interesting to study the orientation statistics of the different flame fronts during propagation. As observed previously in fig. 4.21, in the free case, the PDF of the β angle of the normal to the flame front with the vector issued from the origin has a unimodal shape (maximum around 20 o /30 o ) with a right skewed shape. Over time, the shape remains the same with a decrease in the level of its maximum which, moreover, shifts up to 30 o . Finally, a widening of the curve linked to the development of flame instabilities is also perceptible. In the configurations with obstacles, the shape of the β angle distribution is presented in fig. 4.24 for the staggered and aligned cases.

Several salient points can be observed : first of all, a unimodal shape remains as for the free flame. However, its maximum probability, in the example chosen here, is between 60 o and 70 o . That is to say, much more than for the free flame. This is the direct result of the presence of obstacles that cause major perturbations in the orientation of the flame front. More important, the general shape of the curve also gives us some indication. It can be broken down into three zones : the first zone between 0 o and 45 o , the zone of maximum probability between 45 o and 90 o degrees and finally the tail of the PDF beyond 90 o . The 0 -45 o zone represents the "attack" front, which is mainly oriented towards the 'outside' fresh gas. It is the flame front that progresses between the obstacles and mainly promotes the expansion of the flame. It is interesting to note that due to the very geometric and less disturbed shape (Fig. 4.17) of the flame within the obstacles in line, the probability

Cléante Langrée

University of Rouen -CORIA February 8, 2022