CFD simulation of deflagrations in open congested

environment for industrial accident application

Cléante Langrée

To cite this version:

Cléante Langrée. CFD simulation of deflagrations in open congested environment for industrial ac-
cident application. Fluid mechanics [[Link]-ph]. Normandie Université, 2021. English. �NNT :
2021NORMR085�. �tel-03675828�

HAL Id: tel-03675828

[Link]
Submitted on 25 Jan 2023

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 THÈSE
Pour obtenir le diplôme de doctorat
Spécialité MECANIQUE DES FLUIDES, ENERGETIQUE, THERMIQUE, COMBUSTION,

ACOUSTIQUE

Préparée au sein de l'Université de Rouen Normandie

CFD simulatiοn οf deflagratiοns in οpen cοngested envirοnment

fοr industrial accident applicatiοn
Présentée et soutenue par
CLEANTE LANGREE
Thèse soutenue le 08/12/2021
devant le jury composé de
MAITRE DE CONFERENCES HDR, GROUPE
M. TALIB DBOUK Rapporteur du jury
IMT
DIRECTEUR DE RECHERCHE, ENSMA
M. ARNAUD MURA Rapporteur du jury
POITIERS
INGENIEUR DE RECHERCHE RF, INST. NAT.
M. GUILLAUME LECOCQ Membre du jury
ENVIR. INDUSTRIEL ET RISQUES
PROFESSEUR DES UNIVERSITES, INSA DE
M. RENOU BRUNO Président du jury
ROUEN NORMANDIE
PROFESSEUR DES UNIVERSITES, Université
M. JULIEN REVEILLON Directeur de thèse
de Rouen Normandie
MAITRE DE CONFERENCES HDR,
M. ERWIN FRANQUET Co-directeur de thèse
UNIVERSITE PAU PAYS DE L ADOUR
Thèse dirigée par JULIEN REVEILLON (COMPLEXE DE RECHERCHE
INTERPROFESSIONEL EN AEROTHERMOCHIMIE) et ERWIN FRANQUET (Université
de Rouen Normandie),
ii

PhD, University of Rouen and CORIA - France

This research was done under the supervision of J. Reveillon and E. Franquet with the financial
support of INERIS within a total of 3 years, from 2018 to 2021.

Edition date February 8, 2022

Cléante Langrée University of Rouen - CORIA February 8, 2022

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 éval-
uer 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 per-
met 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.

Cléante Langrée University of Rouen - CORIA February 8, 2022

Abstract - English

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 theo-
retical 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/in-
compressible 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 February 8, 2022

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 [Link] 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

Remerciements

D’abord et avant tout, mes remerciements vont évidemment à l’équipe ayant encadré
cette thèse. Tout d’abord, merci à Guillaume Lecocq, ingénieur recherche supervisant la thèse
à l’INERIS, pour avoir initié ce projet de thèse et avoir tout fait pour que je puisse le mener
à bien. Merci pour sa disponibilité, sa clarté et sa bienveillance à mon égard dès que j’ai
eu besoin de son aide. Enfin merci pour la confiance qu’il m’a accordé en me confiant la
responsabilité de ce travail de thèse.
J’aimerais également remercier Julien Reveillon, de m’avoir accompagné au quotidien
dans ce travail de thèse. Je le remercie de m’avoir aider à organiser mes efforts et d’avoir
tacher de structurer l’organisation de ce travail. Je remercie également Erwin Franquet et
François Xavier Demoulin, pour l’aide et les perspectives nouvelles qu’ils m’ont apporté.
Chacun de ces encadrants m’a aidé et soutenu avec sérieux et bienveillance, et pour cela je
souhaiterais leur exprimer un grand merci.
Mes remerciements vont aussi aux autres membres du jury de thèse : Talib Dbouk,
Arnaud Mura et Bruno Renou. Merci d’avoir évaluer et revu ce travail de thèse. Merci
également pour les nombreuses et pertinentes remarques constructives qu’ils m’ont adressé,
et qui m’ont permis d’améliorer la qualité générale de ce manuscrit.
Ensuite, j’aimerais remercier Emmanuel Leprette et Stéphane Duplantier, respectivement
responsable de l’unité Expérimentation & Modélisation des Explosions et responsable du
pôle Phénomène Dangereux & résistance des Structures. Merci de m’avoir accueilli dans vos
équipes et d’avoir suivi avec intérêt mon avancement durant ces 3 années.
Un grand merci aux équipes de l’INERIS et du CORIA, qui m’ont accueilli en garde
alternée durant ces 3 années de thèse. Tout particulièrement aux ingénieurs des équipes
EMEX et DIEM de l’INERIS, et à l’équipe Spray & Atomisation du CORIA. Qu’il s’agisse de
profiter de la belle forêt de Verneuil le temps d’un footing méridien, ou de quelques sorties au
centre ville de Rouen, ces moments de partage resteront parmi les meilleurs souvenirs de ces
3 années.
Merci également à tous mes collègues doctorants du CORIA, en particulier ceux de
l’équipe Spray & Atomisation : Diego, Victor, Alberto, Leandro, Hector, Emmanuel, Fernando,
Anirudh et Hakim avec qui j’ai eu le plaisir de partager de nombreux "team meetings", et
au moins autant de pauses café. Je n’oublie évidemment pas mes deux collègues de bureau
Hassan et Niccolo. Merci à eux d’avoir indulgé un voisin un peu bruyant et d’avoir participé
aux "white board time" qui faisaient tout le sel de mes journées au laboratoire. Un petit
merci supplémentaire à ceux d’entre eux qui se reconnaitront, d’avoir offert l’asile sportif à un
supporter dépité après l’élimination de son équipe lors de la coupe d’Europe de football 2021.

Cléante Langrée UVCE Simulation February 8, 2022

viii

Enfin, j’aimerais terminer ces remerciements par un énorme merci à ma famille, mes
parents, ma soeur et mes grand-parents pour le soutien et l’amour qu’il m’ont apporté pendant
ces 3 ans ... et à vrai dire, depuis plus longtemps encore. A chaque pas que j’ai fait en 24
ans, je n’ai reçu de leur part que de l’affection, des encouragements et du soutien. Et même
si pour Papy et Mamie je resterai un éternel étudiant, venant leur rendre visite entre deux
semaines de cours, voir de la fierté dans les yeux de ceux par qui j’ai eu le privilège d’être
accompagné jusqu’au jour de ma soutenance est le plus beau diplôme que je puisse obtenir.
J’aimerais finalement adresser tous les remerciements qu’il est possible d’adresser à
Clémence. Vivre avec un doctorant n’est pas chose facile, vivre avec moi encore moins. Merci
d’être restée à mes côtés quand j’avais le plus besoin de toi. Merci de m’avoir supporté (dans
tous les sens du terme) pendant 3 années d’études d’ingénierie et 3 années de doctorat. Et
pour toutes les années qui suivront, d’avance, merci.

Cléante Langrée University of Rouen - CORIA February 8, 2022

Table of Contents

Page

1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Context of this PhD 1
1.2 UVCE 2
1.2.1 Some major UVCEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.2 Usual calculation methods for the explosion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Objective and approach of the study 12

2 Theoretical basis for the numerical simulation of

turbulent combustion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1 Physical description 15
2.1.1 Combustion . . . . . . . . . . . . . . . . . . . . . . . . . . ........................... 15
2.1.2 Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . ........................... 25
2.1.3 Turbulence - Flame interaction . . . . . . . . . . . . . ........................... 26
2.1.4 Unconfined Vapour Cloud Explosion (UVCE) . . . ........................... 28
2.2 Modeling 34
2.2.1 Reactive flow : Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.2.2 Turbulence modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.2.3 Turbulent combustion modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.3 Conclusion 49

3 Numerical tools for UVCE simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

3.1 Introduction 51
3.2 OpenFoam : A brief description 52
3.3 Finite Volume Method 53
3.3.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
3.3.6 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . .................. ....... 61
3.4 Pressure-velocity coupling 63
3.4.1 Semi Implicit Method for Pressure Linked Equations - SIMPLE . . . . . . . . . . . . . . . . . 65
3.4.2 Pressure Implicit with Splitting Operator - PISO . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

Cléante Langrée UVCE Simulation February 8, 2022

x TABLE OF CONTENTS

3.4.3 PIMPLE : PIso + sIMPLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

3.5 XiFoam solver 66
3.5.1 General scheme . . . . . . . . . . . . . . . . . . . . . . . ........................... 66
3.5.2 Continuity equation . . . . . . . . . . . . . . . . . . . . . ........................... 68
3.5.3 Momentum equation . . . . . . . . . . . . . . . . . . . . ........................... 69
3.5.4 Progress variable equation . . . . . . . . . . . . . . . . ........................... 70
3.5.5 Model for mixture inhomogeneity . . . . . . . . . . . ........................... 74
3.5.6 Energy equation . . . . . . . . . . . . . . . . . . . . . . . ........................... 75
3.5.7 Pressure equation . . . . . . . . . . . . . . . . . . . . . . ........................... 76
3.5.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . ........................... 77

4 Flames expanding through obstacles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

4.1 Introduction 79
4.2 Laminar premixed flames 79
4.2.1 1D Configurations . . . . . . . . . . . . . . . . . . . . . . . .......................... 80
4.2.2 Ignition setup . . . . . . . . . . . . . . . . . . . . . . . . . . .......................... 82
4.2.3 Flame velocity . . . . . . . . . . . . . . . . . . . . . . . . . .......................... 82
4.2.4 Pressure and boundary conditions . . . . . . . . . . . . .......................... 86
4.2.5 Expanding spherical flame . . . . . . . . . . . . . . . . . .......................... 87
4.2.6 Configurations . . . . . . . . . . . . . . . . . . . . . . . . . .......................... 88
4.3 Flame expansion 90
4.3.1 Flame structure . . . . . . . . . . . . . . . . . . . . . . . ........... ................ 90
4.3.2 Flame radius . . . . . . . . . . . . . . . . . . . . . . . . . . ........... ................ 90
4.3.3 Mesh convergence . . . . . . . . . . . . . . . . . . . . . . ........... ................ 94
4.3.4 Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . ........... ................ 94
4.3.5 Impact of obstacles . . . . . . . . . . . . . . . . . . . . . ........... ................ 96
4.3.6 Flame structure . . . . . . . . . . . . . . . . . . . . . . . ........... ................ 101
4.3.7 Combustion correlation model . . . . . . . . . . . . . ........... ................ 108
4.4 Conclusion 111

5 Industrial Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113

5.1 Introduction 113
5.2 EXJET experiments 115
5.3 Experimental results 117
5.3.1 Image analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.3.2 Flame velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
5.3.3 Pressure signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
5.4 Numerical Simulation 124
5.4.1 Geometry and meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
5.4.2 Boundary and initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
5.4.3 Global flame propagation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.4.4 Flame velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
5.4.5 Pressure levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
5.5 Conclusion 135

6 Conclusion and perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137

Cléante Langrée University of Rouen - CORIA February 8, 2022

List of Figures

Figure Page

1.1 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. . . . . . . . . . . . 7
1.7 Overpressure against radius abacus from TNT equivalent method [3]. . . . . . . 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 [4]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.1 Fire triangle, or combustion triangle. . . . . . . . . . . . . . . . . . . . . . . . . 16

2.2 Structure of a laminar premixed flame, from [5]. . . . . . . . . . . . . . . . . . . 17
2.3 Flame speed definitions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.4 Flame thickness definitions. From [6]. . . . . . . . . . . . . . . . . . . . . . . . 19
2.5 Flame curvature definitions. From [7]. . . . . . . . . . . . . . . . . . . . . . . . 21
2.6 Curvature effects on flame front. From [6]. . . . . . . . . . . . . . . . . . . . . 22
2.7 Gravity induced instabilities : stable configuration (left) and unstable configu-
ration (right). From [8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.8 Darrieus Landau instabilities. From [9]. . . . . . . . . . . . . . . . . . . . . . . . 23
2.9 Thermo diffusive instabilities. From [10]. . . . . . . . . . . . . . . . . . . . . . . 25
2.10 Kolmogorov’s turbulent energy cascade. . . . . . . . . . . . . . . . . . . . . . . 27
2.11 Borghi diagram, modified by Proust [11]. . . . . . . . . . . . . . . . . . . . . . 28
2.12 Flammability limits. Midbar : stoichiometry. . . . . . . . . . . . . . . . . . . . . 28
2.13 Schematic representation of leak scenarios. . . . . . . . . . . . . . . . . . . . . 29
2.14 Stages of an UVCE from [12]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.15 Turbulence intensity effects on the flame velocity. From [13]. . . . . . . . . . . . 32
2.16 Overpressure wave in a detonation scenario (a) and a deflagration scenario
(b). From [12]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.17 Turbulent spectrum representation of each turbulence modeling method. . . . . 37
2.18 Classification of turbulence modeling methods. Adapted from [14]. . . . . . . . . 37
2.19 Illustration of PDF modeling for turbulent combustion, from [15]. . . . . . . . . 43
2.20 Illustration of a flame surface crossing a control volume [16] . . . . . . . . . . . 45
2.21 wrinkling concept representation. Dashed line : flame front, as captured by
the mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.1 Schematic representation of a control volume. From [17]. . . . . . . . . . . . . 54

3.2 Mesh generated by blockMesh. . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Cléante Langrée UVCE Simulation February 8, 2022

xii LIST OF FIGURES

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. . . . . . . . . . . . . . 57
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. . . . . . . . . . . . . . 58
3.5 Schematic representation of face interpolation. . . . . . . . . . . . . . . . . . . 61

4.1 Schematic representation of the uni dimensional configuration. . . . . . . . . . . 80

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. . . . . . . . . . . . . . . 81
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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.5 Impact of the air mixture definition (molar proportion of oxygen vs nitrogen)
on the velocity bell shape curve. Obtained with Cantera. . . . . . . . . . . . . . 85
4.6 Grid Convergence : impact of the grid step on the premixed flame velocity. . . . 86
4.7 Impact of the waveTransmissive boundary condition on the pressure field. . . 87
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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.10 Average temperature (top) and velocity norm (bottom) profiles over time -
case free-01_00. The colored areas represent the standard deviation of the
corresponding color curve. The vertical red bars indicate the flame position. . . . 92
4.11 Flame radius evolution depending on the selected method to compute flame
front position (rf : burnt/fresh gas equilibrium, re : volume equivalent sphere,
ro : outer radius) - case free-01_00. . . . . . . . . . . . . . . . . . . . . . . . 93
4.12 Flame radius evolution depending on the mesh refinement, from 700k cells
up to 7M cells - case free-01_00. . . . . . . . . . . . . . . . . . . . . . . . . 94
4.13 Turbulent burning velocity St evolution and bending apparition in the obstacle
free environment - case free-01_00. . . . . . . . . . . . . . . . . . . . . . . . 95
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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.16 Azimuthal and elevation cuts of the sta6-0400 flame front, equivalent flame
radius : 2 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
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. 99
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

Cléante Langrée University of Rouen - CORIA February 8, 2022

LIST OF FIGURES xiii

4.21 Time evolution of the PDF angle between the flame normal vector and the
sphere normal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.22 Time evolution of the PDF of the flame radius, sta6-04_00 case. . . . . . . . 103
4.23 Comparison of the PDF of the flame radius : staggered (sta6-04_00) vs in
line obstacles (lin6-05_00) case. . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.24 Comparison of the PDF of the flame surface orientation : staggered (sta6-
04_00) vs in line obstacles lin6-05_00 case. . . . . . . . . . . . . . . . . . . . 105
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. . . . . . . 106
4.26 Azimuthal and elevation cuts of the sta2-02_00 flame front, t = 0.095 s. . . . 107
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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
4.28 Evaluation of various turbulent combustion correlation models. . . . . . . . . . . 109
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 dif-
ferent combustion correlation models. All the flames have the same equivalent
radius : 2 m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

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. . . . . . . . . . . . . . . . . . . . . 113
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. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
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. . . . . . 115
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. . . . . . . . . . 116
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. . . . . . . 118
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. . . . . . . . . . . . 119
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. . . . . . . . . . . . . . . . . . . . . . 121

Cléante Langrée UVCE Simulation February 8, 2022

xiv LIST OF FIGURES

5.8 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). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
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). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
5.10 Arrangement of the crossing between the obstacles. The simplified layout
allows an economy in the number of meshes. . . . . . . . . . . . . . . . . . . . 124
5.11 Paraview screenshot of the two mesh configurations. . . . . . . . . . . . . . . . 126
5.12 Simplification of the ignition process. . . . . . . . . . . . . . . . . . . . . . . . . 126
5.13 Temperature and pressure spatial profiles over the main y axis of the obstacle
module. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
5.14 Comparison of numerical (bottom) and experimental (top) flame front, ex-
perimental time 0.03 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
5.15 Comparison of numerical (bottom) and experimental (top) flame front, ex-
perimental time 0.16 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
5.16 Evolution of numerical and experimental flame front position. The fuchsia
lines have similar slopes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
5.17 Evolution of numerical flame front position for 3 different turbulent flame
speed correlation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
5.18 Numerical (orange) and experimental (blue) pressure signal at 0.28 meters
from the ignition side of the obstacle module. . . . . . . . . . . . . . . . . . . . 133
5.19 Numerical (orange) and experimental (blue) pressure signal at 0.98 meters
from the ignition side of the obstacle module. . . . . . . . . . . . . . . . . . . . 134
5.20 Numerical (orange) and experimental (blue) pressure signal at 1.96 meters
from the ignition side of the obstacle module. . . . . . . . . . . . . . . . . . . . 134

6.1 Turbulent kinetic energy profile along the obstacle module axis. . . . . . . . . . 139

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 [19], Seveso [20], or AZF [21], 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 (Ii ) and the probability of
occurrence (No). Table 1.1 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 – Example of a qualitative scale for the intensity.

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’.

Cléante Langrée UVCE Simulation February 8, 2022

 2 Introduction

Figure 1.1 – Example of a scoring matrix.

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:

• 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’.

1.2 UVCE

An UVCE (Unconfined Vapor Cloud Explosion) is a gas explosion in the open air. In the
case of a flammable gas, such as LPG (Liquefied Petroleum Gas), this explosion produces
thermal and pressure effects. It includes the following steps:

• Release of a fuel into the atmosphere

Cléante Langrée University of Rouen - CORIA February 8, 2022

1.2 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 [12]. 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.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 %

Table 1.2 – Flammability limits of some flammable products encountered in the industrial
processes. Source [22].

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 [23] 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

Cléante Langrée UVCE Simulation February 8, 2022

4 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.

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

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 [24], or jet fires
[25]. 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

1.2 UVCE 5

Figure 1.3 – Characteristics of a shock wave produced by a detonation. From [2].

Overpressure [mBar] Damages on humans Damages on structures

20 Indirect damages Glass/Window destruction
50 Irreversible damages Light damages on buildings
140 Lethal effects Severe damages on buildings
200 Significant lethal effects

Table 1.3 – Overpressure effects thresholds retained by the French legislation for explosion
scenarios.

Cléante Langrée UVCE Simulation February 8, 2022

 6 Introduction

1.2.1 Some major UVCEs

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 [26]

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.

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

Flixborough, 1974 [27]

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

1.2 UVCE 7

Figure 1.5 – Aerial view of the devastation after the Flixborough accident.

Buncefield, 2005 [28]

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 m2 , 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 [29] whereas others prefer the
assumption of a DDT occurrence in a part of the cloud covering dense vegetation [30] [31].

Figure 1.6 – Aerial view of the devastation after the Buncefield accident.

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.

Cléante Langrée UVCE Simulation February 8, 2022

 8 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 [32] and
La Mede in 1992 [33]) 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.

1.2.2 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 [32, 33]. 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 [34, 35]. 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 [36], it will be highly
instrumented whereas the sensors used for large-scale experiments will be less numerous [37].

While gaining knowledge, predictive models were proposed. They can be phenomenolog-
ical 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 MTNT is found by :

Egas
MTNT = a (1.1)
ETNT

with Egas being the potential energy that the flammable mixture can release, ETNT 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.

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

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 Egas 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.

Cléante Langrée UVCE Simulation February 8, 2022

10 Introduction

Congestion Assessment Method [38]

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 propagation - flow of fresh gases ahead of the
flame front - turbulence generation in this zone - flame acceleration’.

The maximum overpressure has to be quantified thanks to a decision tree [39]. 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 [40]. The measured curves were generated
with quasi-hemispherical deflagrations propagating at varying speed.

Baker-Strelhow method [41]

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 [42]. 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 [43] 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.

Multi-Energy Method [4]

It should be kept in mind that an explosion generates high overpressures if the propaga-
tion 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).

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 [4].

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) [44] and George Gabriel Stokes (1819 – 1903) [45]. 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.

Cléante Langrée UVCE Simulation February 8, 2022

 12 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 [46].

In the industrial area, CFD simulation of premixed flames has been used largely for
engine’s combustion chamber studies [47, 48] and explosions [49]. One of the first CFD
modeling explosion was performed in the early 80s [50]. The geometry was a tube of 50 m3
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 [51]. 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 [52, 53]. This is partly explained
by the increase in computational power.

Phenomenological tools can give a result very quickly, simplifying the explosion char-
acterization 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.

1.3 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.

Cléante Langrée UVCE Simulation February 8, 2022

14 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.

2.1 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.

2.1.1 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 :

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
16 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.

Figure 2.1 – Fire triangle, or combustion triangle.

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, auto-
ignition), 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 :
� �� � � �
Y Y Y
CX HY + X + O2 + 3.76N2 → XCO2 + H2 O + 3.76 + X N2 (2.2)
4 2 4

Cléante Langrée University of Rouen - CORIA February 8, 2022

2.1 Physical description 17

Laminar flame

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 [54] 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.

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

In practice, two distinct physical phenomena can drive the diffusion of species in the
flame front :

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
18 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. H2 ) 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 [55]. 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,
Su 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. Fig. 2.3 illustrate new
definitions of "flame speeds".

Figure 2.3 – Flame speed definitions.

Cléante Langrée University of Rouen - CORIA February 8, 2022

2.1 Physical description 19

• The absolute flame speed Sa 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 Sd defined as the speed at which the flame front
propagates with respect to the flow. It is therefore related to the absolute flame speed :

Sd = Sa − u · n (2.3)

The velocity Su is affected by flame front deformations such as stretching (more

detailed therein after in a dedicated paragraph). The notations S u0 will be used to designate
the laminar propagation and burning speed of a plane, adiabatic and non stretched flame.
This characteristic speed Su0 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 [56] 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 [57].

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.4.

Figure 2.4 – Flame thickness definitions. From [6].

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
20 turbulent combustion

A first definition was given by Zeldovich [58] :

λ Dt h
δp = 0 = 0 (2.4)
ρu Cp SL SL

where Dth denotes the thermal diffusivity of fresh gas, such as :

λ
Dth = (2.5)
ρu Cp

Nonetheless, this definition of flame thickness falls short of experimental measurements

and it necessitates knowledge of the thermal properties of the fresh gas Dth and the laminar
burning rate SL0 .

The second most used is the flame front’s thermal thickness, which is defined by
Spalding [59] as :
Tb − Tu
δth = (2.6)
|∇ (T ) |max

Where Tb and Tu are the temperatures of the burnt and unburnt gases, respectively.
Compared to δp , this definition is closer to experimentally available flame thickness measure-
ments. 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,
T − Tu
δtot , as the interval between which the reduced temperature ( ) is between 0.01 and
Tb − Tu
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 [60] :
1 ∂δA
κ= (2.7)
δA ∂t

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 Sf 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 [61] :

κ = κcurvature + κshear = Sf · nC + κshear = Sfn · C + κshear (2.8)

Cléante Langrée University of Rouen - CORIA February 8, 2022

2.1 Physical description 21

Figure 2.5 – Flame curvature definitions. From [7].

where n is the exponent, C is the curvature-related factor and Kshear is the shear-related
factor and are defined as :

� �
1 1
C =∇·n =− + (2.9)
R1 R2

¯ · n = (u · n)C + ∇ · u
Kshear = −n · Ē (2.10)
τ τ

Where 1/R1 and 1/R2 are two characteristic curvature radii, R1 being the largest
curvature radius and R2 being the curvature radius obtained in the normal direction n. Lastly,
¯ stands for the stress tensor, which can be written as 1 �U + U T �. The first term on
Ē
2
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 :

κ = Sfn ∇ · n + (u · n)∇ · n + ∇τ · uτ
= (Sfn − nn )∇ · n + ∇τ · uτ (2.11)
= (Sd )∇ · n + ∇τ · uτ

Candel & Poinsot [62] 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 :

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
22 turbulent combustion

• 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 [6].

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 [63].

Cléante Langrée University of Rouen - CORIA February 8, 2022

2.1 Physical description 23

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

Hydrodynamic instabilities

Darrieus [64] and Landau [65] 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
ρb
speed at which the gases cross the flame front is then multiplied by the ratio, known as
ρu
the expansion factor.

Figure 2.8 – Darrieus Landau instabilities. From [9].

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

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
24 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
r
to the Péclet number, which is the ratio of the flame’s radius to thickness (P e = ). When
δ
the radius is large enough, cell structures appear when the Péclet number exceeds a critical
rcrit
value (P ecrit = ). Bauwens et al. [66] and Kim et al. [67] 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 ecrit ) increases with the
Markstein number (Ma). The same trends were observed in flame speed measurements for
propane/air mixtures by Bauwens et al [66]. 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. [68] 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.1 Physical description 25

Figure 2.9 – Thermo diffusive instabilities. From [10].

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.

2.1.2 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 :

|u|l
Re = (2.12)
ν

Where l and u are the characteristic length and speed of the considered flow, respectively,
and ν being the kinematic viscosity.

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
26 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 [69, 70], 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 lt . The structures in the integral
zone have length scales lt and velocity scales u � (lt ) 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 Ret is introduced :

u � (lt )lt
Ret = (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 u � (r )3
�= r = (2.14)
d u � (r ) r
• The dissipative zone: it has the highest wave numbers and its scale range is limited by
the Kolmogorov scale lK , 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 lK and characteristic velocity uK are calculated as follows :
� � 14
ν3
lK = (2.15)
�
1
uK = (ν�) 4 (2.16)

2.1.3 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

2.1 Physical description 27

Figure 2.10 – Kolmogorov’s turbulent energy cascade.

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
[71], 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 [11] for the situation of industrial UVCE. It is displayed
in the Borghi diagram presented below.

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
28 turbulent combustion

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

2.1.4 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.

H2

CH4

C2 H6
"Fuel"

C3 H8

C4 H10

C6 H12
0 20 40 60 80 100
"Fuel volume fraction"
Non-flammable mixture Flammable mixture

Figure 2.12 – Flammability limits. Midbar : stoichiometry.

Cléante Langrée University of Rouen - CORIA February 8, 2022

2.1 Physical description 29

Step by step description

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

Flammable cloud formation

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 :

• 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.

Figure 2.13 – Schematic representation of leak scenarios.

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 :

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
30 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 :

ST AT
= (2.17)
SU A

With ST and SU being the turbulent and laminar flame speeds respectively and AT 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

2.1 Physical description 31

(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 – Stages of an UVCE from [12].

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
32 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 [72] and Leyer [73] gave a general description of the
drf (τ )
pressure evolution in which ∆P is directly correlated to the flame propagation speed ( )
dt
2
d rf (τ )
and its acceleration ( ):
dt 2

� �� � � �
1 rf (τ ) drf (τ ) r 2 (τ ) d 2 rf (τ )
∆P (r, t) = ρ0 1− 2 + f . (2.18)
α r dt r dt 2

With rf being the flame radius, ρ0 the atmosphere density, α the burnt gases rate of
r
expansion (i.e. ratio of unburnt over burnt density) and τ = t − , with c0 the sound velocity
c0
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 accel-

erating factors can be identified. First the interaction between flame and turbulence is a
well known accelerating factor (fig. 2.15) [55, 74]. 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.

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

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 [75, 76]. 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 [77, 78]. 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.

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
34 turbulent combustion

Figure 2.16 – Overpressure wave in a detonation scenario (a) and a deflagration scenario (b).
From [12].

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.

2.2 Modeling

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

2.2.1 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.

Continuity equation

∂ρ
+ ∇ · (ρU) = 0 (2.19)
∂t

Cléante Langrée University of Rouen - CORIA February 8, 2022

2.2 Modeling 35

Momentum equation

∂ρU
+ ∇ · (ρUU) = −∇p + ∇ · τ (2.20)
∂t

Species mass fraction equation

∂ρYk
+ ∇ · (ρUYk ) = −∇ · J k + ω˙k (2.21)
∂t

Enthalpy equation

∂ρH ∂p � �
+ ∇ · (ρUH) = + ∇ · J h + τ · U + Q̇ (2.22)
∂t ∂t

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 :

2
τ = µ∇U + µ∇U T − µ (∇ · U)δij (2.23)
3

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 Ts :
√
As T
µ= (2.24)
Ts
1+
T

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.

µ
Jk = ∇Yk (2.25)
Sck

� N � �
h µ Pr
J = ∇h + − 1)hk ∇Yk (2.26)
Pr Sck
k=1

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) :

µ
Sck = (2.27)
ρDk

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
36 turbulent combustion

µCp
Pr = (2.28)
λ

Sck λ
Lek = = (2.29)
Pr ρCp Dk

Finally, Q̇ represents the chemical source term, i.e. the heat release by the combustion
reaction. �
Q̇ = − hi ω̇i (2.30)
i

To invert enthalpy in order to retrieve the temperature the janaf model is used [79].
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 :

�� � �
� �
cP = R (a4 T + a3 ) T + a2 T + a1 T + a0 (2.31)

Finally the state equation is used, linking pressure and density :

R
p=ρ T (2.32)
�
M

2.2.2 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

2.2 Modeling 37

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

Figure 2.18 – Classification of turbulence modeling methods. Adapted from [14].

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
38 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, neverthe-
less 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
φ= φ(xi , t) dt (2.33)
T T

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
φ� = = ρ(xi , t)φ(xi , t)dt (2.36)
ρ̄ ρT T

with

� + ��
Φ=Φ (2.37)

ρΦ�� = 0 (2.38)

Cléante Langrée University of Rouen - CORIA February 8, 2022

2.2 Modeling 39

ρφ = ρφ� + ρΦ�� = ρφ� (2.39)

Favre Averaged Navier-Stokes equations

When this operator is applied to our instantaneous Navier Stokes equation (eq. 2.19,
2.20, 2.22, 2.21), the Favre averaged equations are obtained :

∂ρ � �
+ ∇ · ρU� = 0 (2.40)
∂t

∂ρU� � � � �
�
+ ∇ · ρU�U� = −∇ · ρU �� U �� − ∇p + ∇ · τ (2.41)
∂t

∂ρY�k � � � �
�
+ ∇ · ρU�Y�k = ∇ · ρU Yk + ∇ · J k + ω˙k
�� ��
(2.42)
∂t

∂ρH� � � � � � �
� = ∂p + ∇ · ρU
+ ∇ · ρU�H � H + ∇ · J h + τ U + Q̇
�� ��
(2.43)
∂t ∂t

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 �� U �� H ��
k

• 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.

µt � �
��� ��
U Yk = − ∇ · Y�k (2.44)
Sc kt

µt � �
�
U �� ��
H =− �
∇· H (2.45)
Prt

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
40 turbulent combustion

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

� � � �
�
ρU �� U �� = −µ
t
� T − 2 ∇ ∇ · U� ¯Ī + 2 ρk¯Ī
∇U� + (∇U) (2.46)
3 3

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 [80]. 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.

k2
µt = ρCµ (2.47)
�

� � �� � �
∂ρk � =∇· µt
+ ∇ · ρUk µ+ ∇k + Pk − ρ� (2.48)
∂t σk

� � �� � �
∂ρ� � µt � �2
+ ∇ · ρU� = ∇ · µ+ ∇� + C�1 Pk − C�2 ρ Pk (2.49)
∂t σ� k k

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 Pk is also defined as :

�
Pk = −ρ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 [81, 82]. 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 [83] :

a1 k
µt = ρ (2.51)
max(a1 ω, SF2 )

∂ρk � �
� = ∇ · ((µ + σk µt ) ∇k) + Pk − β ∗ ρkω
+ ∇ · ρUk (2.52)
∂t

∂ρω � � 1
�
+ ∇ · ρUω = ∇ · ((µ + σω µt ) ∇ω) − βρω 2 + 2ρ(1 − F1 )σω2 ∇k · ∇ω (2.53)
∂t ω

With F1 and F2 being sensor function, with F1 = 1 in near wall regions, F1 = 0 away
from the surface, F2 = 1 for boundary layer flow and F2 = 0 in free shear layers. This model
also includes several constants with usual values attributed to them : a1 = 0.31, σk = 0.85,
β ∗ = 0.09, β = 0.075, σω = 0.5, σω2 = 0.856. A production term Pk is also defined as :

Pk = mi n(τ ∇ (U) , 10β ∗ kω) (2.54)

2.2.3 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 �� U �� H ��
k

• 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.

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
42 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 [84], 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 [53, 85].

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.

T − Tu
c= (2.55)
Tb − Tu

or

YF − YFu
c= (2.56)
YFb − YF

Cléante Langrée University of Rouen - CORIA February 8, 2022

2.2 Modeling 43

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

∂ρc
+ ∇ · (ρUc) = ∇ · (α∇c) + ω̇c (2.57)
∂t

If eq.2.57 is averaged :

∂ρ�c � � � �
c − ρU�
+ ∇ · ρU�c� = ∇ · α∇� �� ��
c + ω̇c (2.58)
∂t

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).

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

�
ω˙c = ω˙c (c ∗ ) × P (c ∗ , x, t)dc ∗ (2.59)

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
44 turbulent combustion

Turbulent mixing approaches

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 [86, 87], 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 = CEBU ρ c�(1 − c�) (2.60)
k

One can also cite the Bray-Moss-Libby analysis [88], 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.

2 �
ω˙c = c (1 − c�)
ρ� (2.61)
2cm − 1 k

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
[89], 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

2.2 Modeling 45

Figure 2.20 – Illustration of a flame surface crossing a control volume [16]

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 m2 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 Su Σ (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.

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
46 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 :

ω˙c = ρu Su ∇ (c) Ξ (2.64)

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

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 [90] :

∂Ξ
+ Us · ∇Ξ = GΞ − R(Ξ − 1) + (σs − σt )Ξ (2.65a)
∂t
Ξeq − 1
G=R (2.65b)
Ξeq
0.28 Ξ∗eq
R= (2.65c)
τη Ξ∗eq − 1
Ξeq = 1 + 2(1 − b)(Ξ∗eq − 1) (2.65d)
�
u�
∗
Ξeq = 1 + Ξcoeff Rη (2.65e)
Su

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 Su ∇ (b) Ξ (2.66)

Cléante Langrée University of Rouen - CORIA February 8, 2022

2.2 Modeling 47

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

example :

�
u�
Ξ = 1 + 2(1 − b)Ξcoeff Rη (2.67)
Su

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 :

� +∞ � 1
ω̇b
ρu Su = ω̇YF dx = db (2.68)
−∞ 0 �∇b�

In addition, when applied to a turbulent flame, it can be stated that :

� 1
ρu Su Ξ�∇b�
ρu St = db (2.69)
0 �∇b�

� 1
St
= Ξdb (2.70)
Su 0

Now, if Ξ is replaced in 2.70 with the expression in 2.67 :

� 1 �
St u�
= 1 + 2(1 − b)Ξcoeff Rη db (2.71)
Su 0 Su

that is to say :

�
St u�
= 1 + Ξcoeff Rη (2.72)
Su Su

Consequently, with a value of Ξcoeff = 0.6, the turbulent flame speed correlation from
Gulder [91] can be obtained. Logically, other correlations can also be used, such as the Zimont
correlation [92, 93]. In any case, a generic formula for turbulent flame speed velocity can be
written as :

ST u� lt p Cp Cη
= 1 + C( )Cu ( )Cl ( ) Rη (2.73)
Su Su δf pr ef

With C, Cu , Cl , Cp , Cη different correlation constants. Different correlations exist for

this general formula with coefficient adaptations and will be detailed later.

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
48 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 AT . 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 SL in the flamelet regime. The average
turbulent flame front then moves at a turbulent velocity ST equal to the laminar flame velocity
multiplied by the ratio of the total pleated area AT 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", lt , etc.) and the flame (SL0 , δL0 , 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 [94] :

St
= 1 + f (Re, Da, P r ) (2.74)
Su

with, the following definition for the Reynolds, Damköhler and Prandtl numbers

u � lt τt ν
Re = , Da = , P r = (2.75)
ν τf D

Cléante Langrée University of Rouen - CORIA February 8, 2022

 2.3 Conclusion 49

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 :

� � �1/2
St u 1/4
= 1 + CG Ret with CG = 0.62 (2.76)
Su Su

Zimont [92]

Zimont correlation is the turbulent flame speed correlation used in ANSYS-Fluent.

St
= 1 + CZ u � Da1/4 with CZ = 0.52 (2.77)
Su

FLACS [95]

The CFD code FLACS use a specific turbulent flame speed correlation [95].

� �0.716
St u�
= 1 + CF Ret0.196 with CF = 0.96 (2.78)
Su Su

2.3 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.

Cléante Langrée UVCE Simulation February 8, 2022

 Theoretical basis for the numerical simulation of
50 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.

Cléante Langrée University of Rouen - CORIA February 8, 2022

 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 [96], 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 [Link] addition, it is quite difficult
to ensure a mathematical conservation for key quantities [97].
(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 [98]. 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.

3.2 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

 3.3 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
[99]. 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.

3.3 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.

Cléante Langrée UVCE Simulation February 8, 2022

 54 Numerical tools for UVCE simulation

3.3.1 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 Sf be the unit vector orthogonal to the common face Sf 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 .

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

In 3.1, an example of a control volume, or cell, of center P, volume VP , faces fi and

their associated face normal vectors SF i , is presented. The solution of the discrete equations
at the point P, will give the values of any quantity Q (i.e. velocity, mass, energy ...) in the
control volume such as :

�
1
QP = Q = Q(x)dV (3.1)
VP VP

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

3.3 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

56 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

3.3 Finite Volume Method 57

(a) Background mesh. (b) Geometry.

(c) Castellated mesh. (d) Snapped mesh.

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.

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.

Cléante Langrée UVCE Simulation February 8, 2022

 58 Numerical tools for UVCE simulation

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.2 Equations discretization

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

∂ρQ
+ ∇ · (ρUQ) − ∇ · ρΓQ ∇ (Q) = SQ (Q) (3.2)
∂t

It can be integrated over any control volume VP :

� � � �
∂ρQ
dV + ∇ · (ρUQ) dV − ∇ · ρΓQ ∇ (Q) dV = SQ (Q)dV (3.3)
VP ∂t VP VP VP

Then, the Gauss theorem, or divergence theorem, is used to rewrite the equation. The Gauss
theorem can be expressed, for any vector a as :

� �
∇ · (a) dV = [Link] (3.4)
VP δV

with δV the external closed surface of the volume VP , and n the outward pointing normal
vector of the considered surface element. Therefore, eq. 3.3 can be rewritten as :

� � � �
∂ρQ
dV + ρ[Link] − ρΓQ ∇ (Q) .ndS = SQ (Q)dV (3.5)
VP ∂t δV δV VP

Cléante Langrée University of Rouen - CORIA February 8, 2022

 3.3 Finite Volume Method 59

3.3.3 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 :

� �
ρ[Link] = nSf .(ρuQ)f (3.6)
δV f

� �
ρΓQ ∇ (Q) .ndS = nSf .(ρΓQ ∇ (Q))f (3.7)
δV f

Notice that mean values at surface, such as (Qf ), 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 = fx (Q)P + (1 − fx )(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 :

�
SQ (Q)dV = Sc VP + Sp VP QP (3.9)
VP

With Sc , Sp corresponding respectively to the further explicit and implicit part of the
source term. With those approximation our discretized equation can be written :

� � �
∂ρQ
dV + nSf .(ρuQ)f − nSf .(ρΓQ ∇ (Q))f = Sc VP + Sp VP QP (3.10)
VP ∂t
f f

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 [100].

3.3.4 Temporal scheme

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

Cléante Langrée UVCE Simulation February 8, 2022

 60 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 Q0 , and the current time step value
Q itself :

∂Q Q − Q0
= (3.11)
∂t ∆t

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 (Q0 and Q00 , making it a second order "version" of the Euler method.

∂Q 1 3 1
= ( Q − 2Q0 + Q00 ) (3.12)
∂t ∆t 2 2

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.

3.3.5 Spatial schemes

As for the convective term in eq. 3.10, a formulation of the value Qf , 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.

Qf = fx QP + (1 − fx )QN (3.13)

Cléante Langrée University of Rouen - CORIA February 8, 2022

 3.3 Finite Volume Method 61

Figure 3.5 – Schematic representation of face interpolation.

with

fN
fx = (3.14)
PN

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 :

�
QP if F ≥ 0
Qf = (3.15)
QN if F < 0

It is a first order, bounded scheme. It can be extended to second order scheme with a
weighted upwind interpolation for the linearUpwind variant.

3.3.6 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 applica-
tions :

Cléante Langrée UVCE Simulation February 8, 2022

62 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

 3.4 Pressure-velocity coupling 63

waveTransmissive

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 :

�
φP γ
wP = + (3.16)
|Sf | ψP

With φP being the patch face flux, Sf 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 :

∂Q
DDt(W , Q) = + W · ∇ (Q) = 0 (3.17)
∂t

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 :

pP = p0 − 0.5ρ|U|2 (3.18)

nφP
UP = Uc − n([Link] ) + (3.19)
|Sf |

With Up being the velocity at the patch, Uc the velocity in cells adjacent to the patch
and n the patch normal vector.

3.4 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)

Cléante Langrée UVCE Simulation February 8, 2022

64 Numerical tools for UVCE simulation

with M the [n × n] matrix of coefficient for the n cells. For each direction i related to
xi :

     ∂p 
M11 M12 ... M1n (Uxi )1 ( ∂x )1
M21 M22     ∂pi 
 ... M2n  (Uxi )2  ( ∂xi )2 
 ... . =  (3.21)
... ... ...   ...   ... 
Mn1 Mn2 ... Mnn (Uxi )n ∂p
( ∂x i
)n

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 :

∂ρ ∂p
=ψ (3.25)
∂t ∂t

And this expression can be injected in the continuity equation (2.19) :

∂p
ψ + ∇.(ρ(A−1 H(U) − A−1 ∇p)) = 0 (3.26)
∂t

Which give a Poisson equation for the pressure :

∂p
∇.(ρA−1 ∇p) = ψ + ∇.(ρ(A−1 H(U))) (3.27)
∂t

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

3.4.1 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.
• 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).

This loop is often called SIMPLE loop, or outer loop.

3.4.2 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"
• 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).

This loop is often called PISO loop, or inner loop.

3.4.3 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.

Cléante Langrée UVCE Simulation February 8, 2022

 66 Numerical tools for UVCE simulation

3.5 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.

3.5.1 General scheme

XiFoam is a combustion solver based on the b-Xi model [101]. In XiFoam, Navier Stokes
equations are solved (eq. 2.40), (eq. 2.41), (eq. 2.43) in addition to a "b" equations (eq.
2.58), 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.

Cléante Langrée University of Rouen - CORIA February 8, 2022

 3.5 XiFoam solver 67

t=t0

t=t+∆t

rhoEqn.H

UEqn.H, bEqn.H, Eqn.H

pEqn.H

new p guess

p converge ?
no

yes
new U, b, E

Solution converge ?
no

yes

write t data

It’s implemented in OpenFoam in XiFoam with (lst. 3.1) :

Listing 3.1 – XiFoam algorithm

1 while ([Link]())
2 {
3 #include "readTimeControls.H"
4 #include "compressibleCourantNo.H"
5 #include "setDeltaT.H"
6
7 runTime++;
8 Info<< "Time = " << [Link]() << nl << endl;
9
10 #include "rhoEqn.H"
11
12 // --- Pressure -velocity PIMPLE corrector loop
13 while ([Link]())
14 {
15 #include "UEqn.H"
16

Cléante Langrée UVCE Simulation February 8, 2022

 68 Numerical tools for UVCE simulation

17 #include "ftEqn.H"
18 #include "bEqn.H"
19 #include "EauEqn.H"
20 #include "EaEqn.H"
21
22 if (![Link]())
23 {
24 [Link]() == [Link]();
25 }
26
27 // --- Pressure corrector loop
28 while ([Link]())
29 {
30 #include "pEqn.H"
31 }
32
33 if ([Link]())
34 {
35 turbulence ->correct();
36 }
37 }
38
39 rho = [Link]();
40
41 [Link]();
42
43 Info<< "ExecutionTime = " << [Link]() << " s"
44 << " ClockTime = " << [Link]() << " s"
45 << nl << endl;
46 }

3.5.2 Continuity equation

The Favre averaged continuity equation (eq. 2.40) is solved. It’s implemented in XiFoam
with (lst. 3.2) :

∂ρ � �
+ ∇ · ρU� = 0 (3.28)
∂t

Listing 3.2 – Continuity equation - rhoEqn.H

1 fvScalarMatrix rhoEqn
2 (
3 fvm::ddt(rho)
4 + fvc::div(phi)
5 ==
6 fvOptions(rho)
7 );
8
9 [Link](rhoEqn);
10
11 [Link]();
12
13 [Link](rho);

Cléante Langrée University of Rouen - CORIA February 8, 2022

 3.5 XiFoam solver 69

3.5.3 Momentum equation

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

∂ρU� � � � �
�
+ ∇ · ρU�U� = −∇ · ρU �� U �� − ∇p + ∇ · τ + F (3.29)
∂t

The body force source term is neglected.

∂ρU� � � � �
�
+ ∇ · ρU�U� + ∇ · ρU �� U �� − τ = −∇ (p) (3.30)
∂t

The viscous stress tensor τ is defined as :

2 � ij
τ = µ∇U� + µ∇U�T − µ (∇ · U)δ (3.31)
3

And representing the turbulent Reynolds stress with the hypothesis of turbulent viscosity
µt and µef f = µ + µt :

� � � � � �
� T − 2 tr (∇U)I)
−τ = ∇ · µt ∇U� − ∇ · µef f ∇U� − ∇ · µ((∇U) � (3.32)
3

�
Moreover the Reynolds stress tensor is defined as R = U �� U �� . In XiFoam, ρR − τ ,

is implemented via divDevRhoReff in ReynoldsStress.C (lst. 3.3) and TensorI.H (lst.

3.4), both OpenFoam source files.

Listing 3.3 – divDevRhoReff in OpenFoam - ReynoldsStress.C

1 fvc::laplacian
2 (
3 this ->alpha_*rho*this ->nut(),
4 U,
5 "laplacian(nuEff ,U)"
6 )
7 + fvc::div(this ->alpha_*rho*R_)
8 - fvc::div(this ->alpha_*rho*this ->nu()*dev2(T(fvc::grad(U))))
9 - fvm::laplacian(this ->alpha_*rho*this ->nuEff(), U)

Listing 3.4 – dev2 in OpenFoam - TensorI.H

1 inline Tensor <Cmpt> dev2(const Tensor <Cmpt >& t)
2 {
3 return t - SphericalTensor <Cmpt >::twoThirdsI *tr(t);
4 }

And then it’s implemented in XiFoam in UEqn.H (lst. 3.5) :

Cléante Langrée UVCE Simulation February 8, 2022

 70 Numerical tools for UVCE simulation

Listing 3.5 – Momentum equation - UEqn.H

1 [Link](U);
2
3 tmp<fvVectorMatrix > tUEqn
4 (
5 fvm::ddt(rho, U) + fvm::div(phi, U)
6 + [Link](rho, U)
7 + turbulence ->divDevRhoReff(U)
8 ==
9 fvOptions(rho, U)
10 );
11 fvVectorMatrix& UEqn = [Link]();
12
13 [Link]();
14
15 [Link](UEqn);
16
17 if ([Link]())
18 {
19 solve(UEqn == -fvc::grad(p));
20
21 [Link](U);
22 K = 0.5*magSqr(U);
23 }

The [Link](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.

3.5.4 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 � � � �
c ) − ρU�
+ ∇ · ρU�c� = ∇ · α∇ (� �� ��
c + ω̇c (3.34)
∂t

∂ρ�c � �
+ ∇ · ρU�c� = ∇ · (αef f ∇ (�
c )) + ω̇c (3.35)
∂t

A change of variable is then applied :

b =1−c (3.36)

Cléante Langrée University of Rouen - CORIA February 8, 2022

 3.5 XiFoam solver 71

with :

Tb − T
b= (3.37)
Tb − Tu

Therefore :

∂ρ�b � � � � ��
+ ∇ · ρU��
b = ∇ · αef f ∇ �
b + ω̇b (3.38)
∂t

ω̇b is modeled as :

ω̇b = −ρu Su Ξ | ∇�
b| (3.39)

The flame flux term is defined as :

Φst = ρu Su Ξnf (3.40)

Listing 3.6 – Φst calculation - bEqn.H

1 surfaceScalarField phiSt("phiSt", fvc::interpolate(rhou*StCorr*Su*Xi)*nf);

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.

∇�
b
nf = .Sf (3.41)
�
| ∇b |

Listing 3.7 – nf calculation - bEqn.H

1 volVectorField n("n", fvc::grad(b));
2
3 volScalarField mgb(mag(n));
4
5 dimensionedScalar dMgb = 1.0e-3*
6 (b*c*mgb)().weightedAverage(mesh.V())
7 /((b*c)().weightedAverage(mesh.V()) + small)
8 + dimensionedScalar("ddMgb", [Link](), small);
9
10 mgb += dMgb;
11
12 surfaceVectorField SfHat([Link]()/[Link]());
13 surfaceVectorField nfVec(fvc::interpolate(n));
14 nfVec += SfHat*(fvc::snGrad(b) - (SfHat & nfVec));
15 nfVec /= (mag(nfVec) + dMgb);
16 surfaceScalarField nf(([Link]() & nfVec));
17 n /= mgb;

Cléante Langrée UVCE Simulation February 8, 2022

 72 Numerical tools for UVCE simulation

With :

� � � �
| ∇�
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∇ · (3.45)
| ∇�
b| | ∇�
b|

� � � �
∇�b ∇�b
b |= ∇ · �
ρu Su Ξ | ∇� bρu Su Ξ −�
b∇ · ρu Su Ξ (3.46)
| ∇�b| | ∇�b|

� �
ω̇b = ∇ · Φst �
b −�
b∇ · (Φst ) (3.47)

And finally :

∂ρ�b � � � � � � ��
+ ∇ · ρU��
b + ∇ · Φst �
b −�
b∇ · (Φst ) − ∇ · αef f ∇ �
b =0 (3.48)
∂t

It is implemented in XiFoam in bEqn.H (lst. 3.34)

Listing 3.8 – b equation - bEqn.H

1 fvScalarMatrix bEqn
2 (
3 fvm::ddt(rho, b)
4 + mvConvection ->fvmDiv(phi, b)
5 + fvm::div(phiSt , b)
6 - fvm::Sp(fvc::div(phiSt), b)
7 - fvm::laplacian(turbulence ->alphaEff(), b)
8 ==
9 fvOptions(rho, b)
10 );

Cléante Langrée University of Rouen - CORIA February 8, 2022

 3.5 XiFoam solver 73

Su model

In OpenFoam, the laminar flame speed is retrieved via correlation. The Gulders correla-
tions are used :

T α p
Su 0 (T, p, φ) = W × φη × e −ξ(φ−σ) × ( ) × ( )β (3.49)
T0 p0

Where W , η, σ, α, β are the Weller model constants and φ, the equivalence ratio,
Yair
being recovered from ft and ( )sto , the stoichiometric air-fuel ratio.
Yf uel
Listing 3.9 – Flame speed correlation - equivalence ratio - Gulder.C
1 inline Foam::scalar Foam::laminarFlameSpeedModels::Gulders::SuRef
2 (
3 scalar phi
4 ) const
5 {
6 if (phi > small)
7 {
8 return W_*pow(phi, eta_)*exp(-xi_*sqr(phi - 1.075));
9 }
10 else
11 {
12 return 0.0;
13 }
14 }

Listing 3.10 – Flame speed correlation - Pressure/Temperature - Gulder.C

1 inline Foam::scalar Foam::laminarFlameSpeedModels::Gulders::Su0pTphi
2 (
3 scalar p,
4 scalar Tu,
5 scalar phi,
6 scalar Yres
7 ) const
8 {
9 static const scalar Tref = 300.0;
10 static const scalar pRef = 1.013e5;
11
12 return SuRef(phi)*pow((Tu/Tref), alpha_)*pow((p/pRef), beta_)*(1 - f_*Yres);
13 }

Cléante Langrée UVCE Simulation February 8, 2022

 74 Numerical tools for UVCE simulation

3.5.5 Model for mixture inhomogeneity

In addition to the progress variable, a "total fuel mass fraction" ft can be transported.
Note that ft 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).

∂ρft � t ) = ∇.(µ∇ft )
+ ∇.(ρUf (3.50)
∂t

With µ the dynamic viscosity. The effects of equivalence ratio inhomogeneity on the
flame front propagation speed is taken into account by the Su term in the b equation with
the Gulder flame speed correlation (3.49) ([101]).

Yair ft
φ=( )sto (3.51)
Yf uel 1 − ft

Listing 3.11 – total fuel mass fraction equation - ftEqn.H

1 if ([Link]("ft"))
2 {
3 volScalarField& ft = composition.Y("ft");
4
5 fvScalarMatrix ftEqn
6 (
7 fvm::ddt(rho, ft)
8 + mvConvection ->fvmDiv(phi, ft)
9 - fvm::laplacian(turbulence ->alphaEff(), ft)
10 ==
11 fvOptions(rho, ft)
12 );
13
14 [Link](ftEqn);
15
16 [Link]();
17
18 [Link](ft);
19 }

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),
fu (unburnt fuel mass fraction) and fr es the residual fuel at a certain stoichiometric ratio.
They are linked so that :

fu = bft + (1 − b)fr es (3.52)

Yair
ox = 1 − ft − (ft − fu )( )sto (3.53)
Yf uel

pr = 1 − fu − ox (3.54)

Cléante Langrée University of Rouen - CORIA February 8, 2022

 3.5 XiFoam solver 75

Listing 3.12 – XiFoam composition calculation - inhomogeneousMixture.C

1 scalar fu = b*ft + (1.0 - b)*fres(ft, stoicRatio().value());
2 scalar ox = 1 - ft - (ft - fu)*stoicRatio().value();
3 scalar pr = 1 - fu - ox;

With fr es being the residual fuel mass fraction left in the burnt mixture.

1 − ft
fr es = max(ft − ; 0) (3.55)
Yair
( )sto
Yf uel

In summary, mixture composition can be recovered from the two parameters b and ft
as follow (in the case of a stoichiometric mixture) :

b\ft 0 0.055 1
0 Ambiant air Burnt gases Fuel
1 Ambiant air Fresh gases Fuel

Listing 3.13 – Residual burnt fuel mass fraction left calculation - basicCombustionMix-
tureI.H
1 inline Foam::tmp<Foam::volScalarField > Foam::basicCombustionMixture::fres
2 (
3 const volScalarField& ft,
4 const dimensionedScalar& stoicRatio
5 ) const
6 {
7 return max(ft - (scalar(1) - ft)/[Link](), scalar(0));
8 }

3.5.6 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 ht . The exact solved equation is :

∂ρh�t � � ∂ρK� � � � � ��
+ ∇ · ρU�h�t + � − ∂p − ∇ · αef f ∇ · h�t
+ ∇ · ρU�K =0 (3.56)
∂t ∂t ∂t

Listing 3.14 – Unburnt gases energy equation - EauEqn.H

Cléante Langrée UVCE Simulation February 8, 2022

 76 Numerical tools for UVCE simulation

1 fvScalarMatrix heauEqn
2 (
3 fvm::ddt(rho, heau) + mvConvection ->fvmDiv(phi, heau)
4 + (fvc::ddt(rho, K) + fvc::div(phi, K))*rho/[Link]()
5 + (
6 [Link]() == "eau"
7 ? fvc::div
8 (
9 fvc::absolute(phi/fvc::interpolate(rho), U),
10 p,
11 "div(phiv,p)"
12 )*rho/[Link]()
13 : -dpdt*rho/[Link]()
14 )
15 - fvm::laplacian(turbulence ->alphaEff(), heau)
16 ==
17 fvOptions(rho, heau)
18 );

Listing 3.15 – Total energy equation - EaEqn.H

1 fvScalarMatrix EaEqn
2 (
3 fvm::ddt(rho, hea) + mvConvection ->fvmDiv(phi, hea)
4 + fvc::ddt(rho, K) + fvc::div(phi, K)
5 + (
6 [Link]() == "ea"
7 ? fvc::div
8 (
9 fvc::absolute(phi/fvc::interpolate(rho), U),
10 p,
11 "div(phiv,p)"
12 )
13 : -dpdt
14 )
15 - fvm::laplacian(turbulence ->alphaEff(), hea)
16 + fvOptions(rho, hea)
17 );

3.5.7 Pressure equation

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).

Listing 3.16 – Pressure equation - pEqn.H

1 fvScalarMatrix pEqn
2 (
3 fvm::ddt(psi, p)
4 + fvc::div(phiHbyA)
5 - fvm::laplacian(rhorAUf , p)
6 ==
7 fvOptions(psi, p, [Link]())
8 );

Cléante Langrée University of Rouen - CORIA February 8, 2022

 3.5 XiFoam solver 77

3.5.8 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.

Cléante Langrée UVCE Simulation February 8, 2022

78 Numerical tools for UVCE simulation

Cléante Langrée University of Rouen - CORIA February 8, 2022

 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.

4.2 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

Cléante Langrée UVCE Simulation February 8, 2022

 80 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.

4.2.1 1D Configurations

Figure 4.1 – Schematic representation of the uni dimensional configuration.

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.

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

4.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 – 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.

Cléante Langrée UVCE Simulation February 8, 2022

 82 Flames expanding through obstacles

4.2.2 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.

4.2.3 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 :

� � ��
∂T
xf = P os Max | | . (4.1)
∂x
Then, it is possible to estimate
dxf
SL = . (4.2)
dt

Cléante Langrée University of Rouen - CORIA February 8, 2022

4.2 Laminar premixed flames 83

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 – Adiabatic temperature vs equivalence ratio. Comparison of the XiFoam solver
using a progress variable transport equation with Cantera using the GRIMech 3.0 mechanism.

Cléante Langrée UVCE Simulation February 8, 2022

84 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 :

�
Vf = bdv = bAdx , (4.3)
V

with A the constant surface normal to the direction of propagation. The equivalent
distance of the flame front is then introduced to the origin Rf such that the volume Rf · A is
equal to the volume of burned gas Vf . Finally :

�
Rf = bdx , (4.4)
X

and consequently,

dRf
SL = , (4.5)
dt

visible in fig. 4.3, triangle symbols.

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.

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

4.2 Laminar premixed flames 85

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

with the two solvers have been plotted in fig. 4.4. For a lean mixture a good match is observed.
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.

Cléante Langrée UVCE Simulation February 8, 2022

 86 Flames expanding through obstacles

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

4.2.4 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 [102].
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

 4.2 Laminar premixed flames 87

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

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.

4.2.5 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 [103, 104, 105], 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.

Cléante Langrée UVCE Simulation February 8, 2022

 88 Flames expanding through obstacles

The expanding flame is an unsteady flame. Its evolution is the result of different
flame/turbulence interactions over time [106]. 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 [107] (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. [107] 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.

4.2.6 Configurations

(a) - Staggered (b) - Aligned

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.

Cléante Langrée University of Rouen - CORIA February 8, 2022

4.2 Laminar premixed flames 89

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.

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.1.

name Obs St /Sl Coeff u�

01_00 free Gulder 0.62 0.5
01_01 free Gulder 0.62 1
01_02 free Gulder 0.62 1.5
02_00 sta2 Gulder 0.62 0.5
03_00 lin2 Gulder 0.62 0.5
04_00 sta6 Gulder 0.62 0.5
04_04 sta6 Zimont 0.832 0.5
04_05 sta6 Flacs 1.536 0.5
05_00 lin6 Gulder 0.62 0.5

Table 4.1 – 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).

Cléante Langrée UVCE Simulation February 8, 2022

 90 Flames expanding through obstacles

4.3 Flame expansion

4.3.1 Flame structure

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 [108].

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 [Link] 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.

4.3.2 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. [109] proposed to take the radius that corresponds to the equivalent
turbulent flame radius rf where the volume of burned gas vb outside the sphere by this radius
is equal to the volume of fresh gas inside this same sphere. By defining ro , the largest radius
encompassing the maximum position of the burned gases and ri the minimum radius that
does not contain any fresh gases, the average flame radius r f is such that :

� rf � ro
(1 − c) dv = c dv (4.6)
ri rf

Cléante Langrée University of Rouen - CORIA February 8, 2022

4.3 Flame expansion 91

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.

Cléante Langrée UVCE Simulation February 8, 2022

92 Flames expanding through obstacles

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

Cléante Langrée University of Rouen - CORIA February 8, 2022

4.3 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 re = f (Vb ) including the presence of obstacles. The volume
of burnt
� gases Vb is defined by the integral of the progress variable c over the whole domain,
Vb = Ω 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 ro , the ’outer’ radius defined above.

Figure 4.11 – Flame radius evolution depending on the selected method to compute flame
front position (rf : burnt/fresh gas equilibrium, re : volume equivalent sphere, ro : outer
radius) - case free-01_00.

It is interesting to consider the three possible definitions of the flame radius simulta-
neously and to compare their evolution in an unobstructed case. The result can be seen in
Fig. 4.11. Naturally, the outer radius ro that marks the farthest distance from the flame
to the center of the sphere is the largest. The radius rf recommended by Bradley [109]
because of its balance between fresh and burning gases volumes comes next. Finally, the
smallest radius is the equivalent radius re 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.

Cléante Langrée UVCE Simulation February 8, 2022

 94 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.

4.3.3 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.

4.3.4 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 [110] proposed the principle of flame front wrinkling as the main

Cléante Langrée University of Rouen - CORIA February 8, 2022

4.3 Flame expansion 95

Figure 4.13 – Turbulent burning velocity St evolution and bending apparition in the obstacle
free environment - case free-01_00.

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 [111], 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 [112]. 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 :

ρb dre
St = . (4.7)
ρu dt

According to accurate DNS analysis carried out in [113], 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.

Cléante Langrée UVCE Simulation February 8, 2022

 96 Flames expanding through obstacles

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

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.

4.3.5 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

4.3 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.15 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 propa-
gates through the obstacles, which naturally generates a strong deformation of its front along
the azimuthal plane (plane perpendicular to the obstacles - fig. 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 m2 , 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 m2 . 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 m2 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. [111]. On the other hand, the four cases with obstacles

Cléante Langrée UVCE Simulation February 8, 2022

98 Flames expanding through obstacles

(a) - sta2

(b) - sta6

Figure 4.15 – Flame propagation, staggered obstacles (2 and 6 cm radius), equivalent flame
radius : 2 m.

Cléante Langrée University of Rouen - CORIA February 8, 2022

4.3 Flame expansion 99

Figure 4.16 – Azimuthal and elevation cuts of the sta6-0400 flame front, equivalent flame
radius : 2 m.

Figure 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.

Cléante Langrée UVCE Simulation February 8, 2022

100 Flames expanding through obstacles

Figure 4.18 – Flame surface evolution.

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.

Cléante Langrée University of Rouen - CORIA February 8, 2022

 4.3 Flame expansion 101

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.

4.3.6 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.

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 – Time evolution of the PDF of the flame radius.

Cléante Langrée UVCE Simulation February 8, 2022

102 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 25o .
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 60o . 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. [108].

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

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

4.3 Flame expansion 103

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

At time 0.045 s, the first double peak (orange curve) corresponds to the crossing of
the first obstacle row positioned at 0.8 m from the center. Indeed, a part of the flame front
stagnates in front of the obstacle while another part propagates in the gap between the
bars and forms a second front. It is these two fronts that cause the double peak to appear.
After this first row of bars, the propagating flame front shows disturbances and continues
to advance towards the second row of bars. From this level on, the ongoing obstacle bypass
presents three levels of flame front presence that can be observed very clearly on the PDFs
presented in fig. 4.23, blue curve. With a first peak before the obstacle (flame ’blocked’ by
the bar), a second peak behind the obstacle (the two flame fronts bypassing the bars meet)
and finally a last peak corresponding to the front that has not met the previous obstacle and
that has already reached the free space between the bars of the next row.

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).

Cléante Langrée UVCE Simulation February 8, 2022

104 Flames expanding through obstacles

Figure 4.23 – Comparison of the PDF of the flame radius : staggered (sta6-04_00) vs in
line obstacles (lin6-05_00) case.

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 20o /30o ) 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 30o . 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
60o and 70o . 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 0o and 45o , the zone of maximum
probability between 45o and 90o degrees and finally the tail of the PDF beyond 90o .

The 0 − 45o 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

4.3 Flame expansion 105

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

of having a ’lateral’ front is much lower (Fig. 4.24, orange curve). It even shows a plateau
between 10o and 40o . 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 60o ).

Finally, the third zone of the PDF P (β) concerns angles greater than 90o , 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.

Cléante Langrée UVCE Simulation February 8, 2022

106 Flames expanding through obstacles

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.

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 : (1)
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

4.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 – Azimuthal and elevation cuts of the sta2-02_00 flame front, t = 0.095 s.

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 0o and 45o ). 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 90o (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 90o 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.

Cléante Langrée UVCE Simulation February 8, 2022

 108 Flames expanding through obstacles

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.7 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 ([114] 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

4.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 (Cg = 0.62) and the Zimont and FLACS parameters are Cz = 0.96 and
Cf = 1.536, respectively.

Figure 4.28 – Evaluation of various turbulent combustion correlation models.

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.

Cléante Langrée UVCE Simulation February 8, 2022

110 Flames expanding through obstacles

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.

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 45o and 90o , 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 > 90o ) 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.

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

 4.4 Conclusion 111

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.

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.

4.4 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.

Cléante Langrée UVCE Simulation February 8, 2022

112 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.

Cléante Langrée University of Rouen - CORIA February 8, 2022

 Chap. 5 | Industrial Application

5.1 Introduction

The majority of methods to evaluate the consequences of an explosion are based on

experimental results [115]. 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.

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.

Cléante Langrée UVCE Simulation February 8, 2022

114 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.

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.

Cléante Langrée University of Rouen - CORIA February 8, 2022

 5.2 EXJET experiments 115

Figure 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.2 EXJET experiments

The experimental work considered here has been carried out by INERIS and GDF-SUEZ
([18]). 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 con-
figuration 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 confi[Link] 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 [Link] 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.

Cléante Langrée UVCE Simulation February 8, 2022

116 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 experi-
mental 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. Figure
5.2 (blue discs) depicts their exact location.
• 2 high frequency cameras (frequency : 2000 images per second) were used to follow
the flame front propagation.

(a) - Initial picture (b) - Final picture

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.

Cléante Langrée University of Rouen - CORIA February 8, 2022

 5.3 Experimental results 117

5.3 Experimental results

5.3.1 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.

Cléante Langrée UVCE Simulation February 8, 2022

118 Industrial Application

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.

Cléante Langrée University of Rouen - CORIA February 8, 2022

5.3 Experimental results 119

Figure 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.

Cléante Langrée UVCE Simulation February 8, 2022

 120 Industrial Application

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.

5.3.2 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

5.3 Experimental results 121

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.

Cléante Langrée UVCE Simulation February 8, 2022

122 Industrial Application

Figure 5.8 – 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).

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).

Cléante Langrée University of Rouen - CORIA February 8, 2022

 5.3 Experimental results 123

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.

5.3.3 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).

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

Cléante Langrée UVCE Simulation February 8, 2022

 124 Industrial Application

figure 5.9. In ascending order of their position, these are the colors blue, orange and yellow.
Although the signal is slightly noisy, the shape of the profile is easily recognizable with a rise
in pressure upstream of the flame passage then a decrease before a slight depression. The
pressure signals are very similar to each other. It shows an overpressure of about 6mbar and
a depression of 2 mbar . The data are relative to the atmospheric pressure.

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.

(a) - Original disposition (b) - Simplified disposition

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

5.4 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

 5.4 Numerical Simulation 125

5.4.1 Geometry and meshing

Contact between obstacles

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.

5.4.2 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.

Field wall ambiant air

p zeroGradient waveTransmissive
U fixed value : null waveTransmissive
b, k, ω zeroGradient zeroGradient

Table 5.1 – Boundary condition type used for the EXJET configuration simulation.

Cléante Langrée UVCE Simulation February 8, 2022

126 Industrial Application

(a) - Original mesh

(b) - Simplified mesh

Figure 5.11 – Paraview screenshot of the two mesh configurations.

(a) - Ignition : original disposition (b) - Ignition : simplified disposition

Figure 5.12 – Simplification of the ignition process.

Cléante Langrée University of Rouen - CORIA February 8, 2022

 5.4 Numerical Simulation 127

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 ft = 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 ft , 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.

5.4.3 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.

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.

(a) - time 0.025 s (b) - time 0.05 s

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

Cléante Langrée UVCE Simulation February 8, 2022

 128 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.

5.4.4 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

5.4 Numerical Simulation 129

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

Cléante Langrée UVCE Simulation February 8, 2022

130 Industrial Application

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

Cléante Langrée University of Rouen - CORIA February 8, 2022

5.4 Numerical Simulation 131

Figure 5.16 – Evolution of numerical and experimental flame front position. The fuchsia lines
have similar slopes.

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

Cléante Langrée UVCE Simulation February 8, 2022

 132 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.16. The blue
curve represents the Gulder correlation that used by default.

The green curve represents the FLACS correlation, used with the coefficient CF = 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 CZ = 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 [94]:

St
= 1 + f (Re, Da, P r ) , (5.1)
Sl

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.

5.4.5 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

5.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).

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.

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

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 UVCE Simulation February 8, 2022

134 Industrial Application

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

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

Cléante Langrée University of Rouen - CORIA February 8, 2022

 5.5 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.

5.5 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.

Cléante Langrée UVCE Simulation February 8, 2022

136 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.

Cléante Langrée UVCE Simulation February 8, 2022

138 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

 139

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.

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

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, hydro-
gen/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

Cléante Langrée UVCE Simulation February 8, 2022

140 Conclusion and perspectives

Cléante Langrée University of Rouen - CORIA February 8, 2022

Bibliography

[1] S. Trelat. Impact de fortes explosions sur les batiments representatifs d’une installation
industrielle. Phd, 2006.

[2] S. Eveillard. Propagation d’une onde de choc en presence d’une barriere de protection.
Phd, 2013.

[3] TM5-1300. U.s. department of the army, s"structure to resist the effects of acci-
dental explosions", army tm 5-1300, navy navfac p-397, afr 88-22. Technical report,
Washington D.C., Department of the Army, Navy and Air Force, 2008.

[4] A.C. van den Berg. The multi-energy method - a framework for vapour cloud explosion
blast prediction. Technical report, TNO Prins Maurits Laboratory, report PML 1984-
C72, 1984.

[5] T.C. Treurniet, F.T.M. Nieuwstadt, and B.J. Boersma. Direct numerical simulation
of homogeneous turbulence in combination with premixed combustion at low Mach
number modelled by the $G$-equation. Journal of Fluid Mechanics, 565, 10 2006.

[6] A. 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. PhD
thesis, Rouen, INSA, 2016.

[7] C.K. Law. Dynamics of stretched flames. Symposium (International) on Combustion,

22(1):1381–1402, 1989.

[8] C. Endouard. Etude expérimentale de la dynamique des flammes de prémélange

isooctane/air en expansion laminaire et turbulente fortement diluées. Theses, Université
d’Orléans, November 2016.

[9] B. Galmiche. Caractérisation expérimentale des flammes laminaires et turbulentes en

expansion. PhD thesis, Université d’Orléans, 2014.

[10] J. Goulier. Comportements aux limites de flammes de prémélange hydrogène/air/dilu-

ants. Étude de la transition flamme laminaire-flamme turbulente. PhD thesis, Thesis,
University of Orléans, 2015.

[11] C. Proust. Formation, inflammation, combustion des atmosphères explosives (atex) et

effets associés. Mémoire d’HdR présenté à l’Institut National Polytechnique de Lorraine,
12, 2004.

Cléante Langrée UVCE Simulation February 8, 2022

142 BIBLIOGRAPHY

[12] C. Proust J. Daubech, E. Leprette. 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. Technical report, Institut National de l’Environnement Industriel et des Risques,
2016.

[13] R.C. Aldredge, V. Vaezi, and P.D. Ronney. Premixed-flame propagation in turbulent
taylor–couette flow. Combustion and flame, 115(3):395–405, 1998.

[14] P. Sagaut, M. Terracol, and S. Deck. Multiscale and multiresolution approaches in

turbulence-LES, DES and Hybrid RANS/LES Methods: Applications and Guidelines.
World Scientific, 2013.

[15] H Pitsch. Cefrc combustion summer school - modeling turbulent combustion, 2014.

[16] U. Ahmed, N. Chakraborty, and M. Klein. On the stress-strain alignment in premixed

turbulent flames. Scientific reports, 9(1):1–9, 2019.

[17] J. Guerrero. Introduction to Computational Fluid Dynamics: Governing Equations,

Turbulence Modeling Introduction and Finite Volume Discretization Basics. 07 2015.

[18] J. Sail, V. Blanchetiere, B. Geniaut, K. Osman, J. Daubech, D. Jamois, and J. Hebrard.

Review of knowledge and recent works on the influence of initial turbulence in methane
explosion. page 401. GexCon AS, June 2014.

[19] E. Broughton. The Bhopal disaster and its aftermath: a review. Environmental Health,
4(1):6, May 2005.

[20] E. Homberger, G. Reggiani, J. Sambeth, and H.K. Wipf. The seveso accident: its
nature, extent and consequences. The Annals of Occupational Hygiene, 22(4):327–370,
January 1979.

[21] First lessons of the Toulouse ammonium nitrate disaster, 21st September 2001, AZF
plant, France. Journal of Hazardous Materials, 111(1-3):131–138, July 2004.

[22] J-M. Petit and J-L. Poyard. Ed911 - les melanges explosifs - 1. gaz et vapeurs. Technical
report, INRS, 2004.

[23] K. Gugan. Unconfined vapour cloud explosions. Institute of Chemical Engineers.

Reports ISBN 0 85295 114 0 and 0 7114 5102 8, 1979.

[24] S. Jolly, F. Prats, G. Chantelauve, and S. Duplantier. Formalisation du savoir et des

outils dans le domaine des risques majeurs (eat-dra-76) modélisations de feux industriels
- ω 2. Technical report, Institut National de l’Environnement Industriel et des Risques,
2014.

[25] S. Jolly, G. Leroy, B. Truchot, F. Prats, and S. Duplantier. Formalisation du savoir et

des outils dans le domaine des risques majeurs (eat-dra-76) feu torche - ω 8. Technical
report, Institut National de l’Environnement Industriel et des Risques, 2014.

[26] D.S. Burgess and M.G. Zabetakis. Detonation of a flammable cloud following a propane
pipeline break: the december 9, 1970, explosion in port hudson, mo. 1 1973.

[27] J.E.S. Venart. Flixborough: the explosion and its aftermath. Process Safety and
Environmental Protection, 82(2):105 – 127, 2004. Fire and Explosion.

Cléante Langrée University of Rouen - CORIA February 8, 2022

BIBLIOGRAPHY 143

[28] M.J. Johnson. The potential for vapour cloud explosions – lessons from the buncefield
accident. Journal of Loss Prevention in the Process Industries, 23(6):921 – 927, 2010.
Papers Presented at the 2009 International Symposium of the Mary Kay O’Connor
Process Safety Center.
[29] Atkinson G., Hall J., and McGillivray A. Hse research report rr1113 : Review of vapor
cloud explosion incidents. Technical report, Health and Safety Executive, 2017.
[30] S. Davis, J. Pagliaro, D. Botwinick, T. DeBold, K. Wingerden, D. Allason, and D.M.
Johnson. Do not believe the hype: Using case studies and experimental evidence to
show why the HSE is wrong about excluding deflagration-to-detonation transitions.
Process Safety Progress, 38(2):e11998, June 2019.
[31] M. Abdel-jawad and F. Gavelli. The myth of episodic deflagration in large, unconfined
vapor clouds. Process Safety Progress, 40(1), March 2021.
[32] J-F. Lechaudel and Y. Mouilleau. Assessment of an accidental vapour cloud explosion -
a case study : Saint herblain. Proceedings of the 8 th International Loss Prevention
Symposium, pages 333–348, 1995.
[33] P. Michaelis et al. Methods applied to investigate the major uvce that occured in
the total refinery’s fluidcatalytic cracking unit at la mede. Proceedings of the 8 th
International Loss Prevention Symposium, pages 365–376, 1995.
[34] S.G. Davis, P. Hinze, O.R. Hansen, and K. Van Wingerden. Investigation techniques
used to determine the massive vapor cloud explosion at the buncefield fuel depot.
Technical report, GexCon U.S., 2010.
[35] J. Daubech, C. Proust, and G. Lecocq. Propagation of a confined explosion to an
external cloud. Journal of Loss Prevention in the Process Industries, 49:805–813, 2017.
[36] J. Daubech, J. Hebrard, S. Jallais, E. Vyazmina, D. Jamois, and F. Verbecke. Un-ignited
and ignited high pressure hydrogen releases: concentration - turbulence mapping and
overpressure effects. International symposium on hazards, prevention, and mitigation
of industrial explosions, X:93–104, 2014.
[37] Puttock J.S., Cresswell T.M., Marks P.R., Samuels, and Prothero A. Explosion
assessment in confined vented geometries. solvex large-scale explosion tests and scope
model development. project report. Technical report, Health and Safety Executive,
OTO 96 004 (Shell Research Limited, Rep. TRCP 3688R2), 1996.
[38] Developments in the Congestion Assessment Method for the prediction of vapour-cloud
explosions. Loss Prevention and Safety Promotion in the Process Industries, pages
1107–1133, January 2001.
[39] J.S. Puttock. Fuel gas explosions guidelines - the congestion assessment method.
Institution of Chemical Engineers. Symposium Series Num 139. EFCE Event num 548.
EFECEPublication num 113, pages 267–284, 1995.
[40] W.P.M. Mercx. Overall final report of the merge project. Technical report, CEC
contract STEP-CT-0111 (SSMA), 1993.
[41] Q.A. Baker, C.M. Doolittle, G.A. Fitzgerald, and M.J. Tang. Recent developments in
the Baker-Strehlow VCE analysis methodology. Process Safety Progress, 17(4):297–301,
1998.

Cléante Langrée UVCE Simulation February 8, 2022

144 BIBLIOGRAPHY

[42] Strehlow R. Blast wave generated by spherical flames. Combustion and Flame,
35:297–310, 1979.

[43] K. Van Wingerden. Investigation into blast produced by non-steady flames. Technical
report, Prins Maurits Laboratory Report - C- 66, 1984.

[44] C. Navier. Mémoire sur les lois du mouvement des fluides. Mémoires de l’Académie
Royale des Sciences de l’Institut de France, 6(1823):389–440, 1823.

[45] G.G. Stokes. printed in 1847. On the theories of the internal friction of fluids in motion
and of the equilibrium and motion of elastic solids.’Transactions of the Cambridge
Philosophical Society, 8:287–319, 1845.

[46] D. Veynante and L. Vervisch. Turbulent combustion modeling. Progress in Energy and
Combustion Science, page 74, 2002.

[47] S. Richard et al. Towards large eddy simulation of combustion in spark ignition engines.
Proc. Combust. Inst., 31:3059–3066, 2007.

[48] L. Gicquel, G. Staffelbach, and T. Poinsot. Large eddy simulations of gaseous flames
in gas turbine combustion chambers. Progress in Energy and Combustion Science,
38:782–817, 2012.

[49] O.R. Hansen, P. Hinze, D. Engel, and S. Davis. Using computational fluid dynamics
(cfd) for blast wave predictions. Journal of Loss Prevention in the Process Industries,
23:885–906, 2010.

[50] B.H. Hjertager. Computed simulation of turbulent reactive gas dynamics. Modeling,
Identification and Control, 5:211–236, 1985.

[51] B.H. Hjertager, T. Solberg, and K.O. Nymoen. Computer modelling of gas explosion
propagation in offshore modules. Journal of Loss Prevention in the Process Industries,
5:165–174, 1992.

[52] C.R. Bauwens, J. Chaffee, and S.B. Dorofeev. Vented explosion overpressures from
combustion of hydrogen and hydrocarbon mixtures. International Journal of Hydrogen
Energy, 36(3):2329–2336, February 2011.

[53] P. Quilattre. Simulation aux grandes echelles d’explosions en domaine semi-confine.

Phd, 2014.

[54] E. Mallard and H.L. Chatelier. Recherches expérimentales et théoriques sur la combus-
tion des mélanges gazeux explosives. Dunod, H. and Pinat, E., 1883.

[55] T. Poinsot and D. Veynante. Theoretical and numerical combustion. CNRS, Toulouse,
3. ed edition, 2011.

[56] F.A. Williams. Combustion Theory. CRC Press, March 2018.

[57] E. Varea. Experimental analysis of laminar spherically expanding flames. PhD thesis,
03 2013.

[58] Y.B. Zeldovich. The theory of combustion and detonation. In Academy of Sciences,
1944.

Cléante Langrée University of Rouen - CORIA February 8, 2022

BIBLIOGRAPHY 145

[59] D.B. Spalding. Some fundamentals of combustion, volume 2. Butterworths scientific

publications, 1955.

[60] F.A. Williams. A review of some theoretical considerations of turbulent flame structure.
In AGARD Conference Proceeding, 1975, 1975.

[61] M. Matalon. On flame stretch. Combustion science and technology, 31(3-4):169–181,

1983.

[62] S. Candel and T. Poinsot. Flame Stretch and the Balance Equation for the Flame
Area. Combustion Science and Technology, 70(1-3):1–15, March 1990.

[63] P.D. Ronney and H.Y. Wachman. Effect of gravity on laminar premixed gas combustion
I: Flammability limits and burning velocities. Combustion and Flame, 62(2):107–119,
November 1985.

[64] G Darrieus. Propagation d’un front de flamme. La Technique Moderne, 30:18, 1938.

[65] L Landau. On the theory of slow combustion. Acta Physicochimica. URSS, 19:77 –
85, 1944.

[66] C.R. Bauwens, J.M. Bergthorson, and S.B. Dorofeev. Experimental study of spherical-
flame acceleration mechanisms in large-scale propane–air flames. Proceedings of the
Combustion Institute, 35(2):2059–2066, 2015.

[67] W.K. Kim, T. Mogi, and R. Dobashi. Effect of propagation behaviour of expanding
spherical flames on the blast wave generated during unconfined gas explosions. Fuel,
128:396–403, July 2014.

[68] Y.A. Gostintsev, A.G. Istratov, and Y.V. Shulenin. Self-similar propagation of a free
turbulent flame in mixed gas mixtures. Combustion, Explosion and Shock Waves,
24(5):563–569, 1988.

[69] A. Kolmogorov. The local structure of turbulence in incompressible viscous fluid for
very large reynolds numbers. Cr Acad. Sci. URSS, 30:301–305, 1941.

[70] A. Kolmogorov. Dissipation of energy in the locally isotropic turbulence. In Dokl. Akad.
Nauk SSSR A, volume 32, pages 16–18, 1941.

[71] R. Borghi. On the structure and morphology of turbulent premixed flames. In Recent
advances in the Aerospace Sciences, pages 117–138. Springer, 1985.

[72] B. Deshaies and J.C. Leyer. Flow field induced by unconfined spherical accelerating
flames. Combustion and flame, 40:141–153, 1981.

[73] J.C. Leyer. Effets de pression engendres par l’explosion dans l’atmosphere de melanges
gazeux d’hydrocarbures et d’air. 1982.

[74] D. Lee and K.Y. Huh. Validation of analytical expressions for turbulent burning velocity
in stagnating and freely propagating turbulent premixed flames. Combustion and flame,
159(4):1576–1591, 2012.

[75] J. Zeeuwen, C. V. Wingerden, and R. Dauwe. Experimental investigation into the blast
effect produced by unconfined vapour cloud explosions. 1983.

Cléante Langrée UVCE Simulation February 8, 2022

146 BIBLIOGRAPHY

[76] C.J.M. Van Wingerden and J.P. Zeeuwen. Flame propagation in the presence of
repeated obstacles: influence of gas reactivity and degree of confinement. Journal of
hazardous materials, 8(2):139–156, 1983.

[77] J. 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. PhD thesis, Orléans, 2008.

[78] J. Zhou, K. Nishida, T. Yoshizaki, and H. Hiroyasu. Flame propagation characteristics

in a heterogeneous concentration distribution of a fuel-air mixture. SAE transactions,
pages 1440–1460, 1998.

[79] A. Burcat. Thermochemical Data for Combustion Calculations, pages 455–473. Springer
New York, New York, NY, 1984.

[80] B.E. Launder and B.I. Sharma. Application of the energy-dissipation model of turbulence
to the calculation of flow near a spinning disc. Letters in heat and mass transfer,
1(2):131–137, 1974.

[81] D.C. et al. Wilcox. Turbulence modeling for CFD, volume 2. DCW industries La
Canada, CA, 1998.

[82] D.C. Wilcox. Formulation of the kw turbulence model revisited. AIAA journal,
46(11):2823–2838, 2008.

[83] F.R. Menter. Two-equation eddy-viscosity turbulence models for engineering applica-
tions. AIAA Journal, 32(8):1598–1605, August 1994.

[84] Ducros Colin, O., Veynante F., D., and T. Poinsot. A thickened flame model for large
eddy simulations of turbulent premixed combustion. Physics of fluids, 12(7):1843–1863,
2000.

[85] O. Vermorel, P. Quillatre, and T. Poinsot. LES of explosions in venting chamber: A test
case for premixed turbulent combustion models. Combustion and Flame, 183:207–223,
September 2017.

[86] B.F. Magnussen and B.H. Hjertager. On mathematical models of turbulent combus-
tion with special emphasis on soot formation and combustion. In 16 th International
Symposium on Combustion. The Combustion Institute, Pittsburgh, PA, 1976.

[87] B.F. Magnussen. On the structure of turbulence and a generalized eddy dissipation
concept for chemical reaction in turbulent flow. In 9th AIAA Meeting, St. Louis, 1981.

[88] K.N.C. Bray, P.A. Libby, and J.B. Moss. Unified modeling approach for premixed
turbulent combustion—Part I: General formulation. Combustion and Flame, 61(1):87–
102, 1985.

[89] F.A. Williams. Turbulent combustion. In The mathematics of combustion, pages

97–131. SIAM, 1985.

[90] H.G. Weller. The development of a new flame area combustion model using conditional
averaging. Thermo-fluids section report TF, 9307, 1993.

[91] Ö.L. Gülder. Turbulent premixed flame propagation models for different combustion
regimes. Symposium (International) on Combustion, 23(1):743–750, January 1991.

Cléante Langrée University of Rouen - CORIA February 8, 2022

BIBLIOGRAPHY 147

[92] V.L. Zimont. Gas premixed combustion at high turbulence. turbulent flame closure
combustion model. Experimental Thermal and Fluid Science, page 8, 2000.

[93] V.L. Zimont, F. Biagioli, and K. Syed. Modelling turbulent premixed combustion in the
intermediate steady propagation regime. Progress in Computational Fluid Dynamics,
An International Journal, 1(1/2/3):14, 2001.

[94] P. Flohr and H. Pitsch. A turbulent flame speed closure model for LES of industrial
burner flows. In Center for Turbulence Research, Proceedings of the Summer Program,
pages 169–179, 2000.

[95] GEXCON AS. Flacs-cfd v21.1 user’s manual. Technical report, 2021.

[96] J.H. Ferziger and M. Perić. Computational methods for fluid dynamics. Springer, Berlin
; New York, 3rd, rev. ed edition, 2002.

[97] Peter H. Aspects of conservation in finite element flow computations. Computer

Methods in Applied Mechanics and Engineering, 117(3):423–437, 1994.

[98] J.H. Ferziger and M. Perić. Computational Methods for Fluid Dynamics. Springer, 2nd
edition, 1999.

[99] H.G. Weller, Tabor, G., H. Jasak, and C. Fureby. A tensorial approach to computational
continuum mechanics using object-oriented techniques. Computers in Physics, 12(6),
1998.

[100] H. Jasak. Error analysis and estimation for the finite volume method with applications
to fluid flows. Direct, M, 01 1996.

[101] H.G. Weller, G. Tabor, A.D. Gosman, and C. Fureby. Application of a flame-wrinkling
les combustion model to a turbulent mixing layer. Symposium (International) on
Combustion, 27(1):899 – 907, 1998. Twenty-Seventh Sysposium (International) on
Combustion Volume One.

[102] F. Piscaglia, A. Montorfano, and A. Onorati. Development of a non-reflecting boundary

condition for multidimensional nonlinear duct acoustic computation. Journal of Sound
and Vibration, 332(4):922–935, 2013.

[103] V. Di Sarli, A. Di Benedetto, and G. Russo. Using Large Eddy Simulation for under-
standing vented gas explosions in the presence of obstacles. Journal of Hazardous
Materials, 169(1-3):435–442, September 2009.

[104] V. Di Sarli, A. Di Benedetto, and G. Russo. Sub-grid scale combustion models for large
eddy simulation of unsteady premixed flame propagation around obstacles. Journal of
Hazardous Materials, 180(1-3):71–78, August 2010.

[105] A.M. Coates, D.L. Mathias, and B.J. Cantwell. Numerical investigation of the effect
of obstacle shape on deflagration to detonation transition in a hydrogen–air mixture.
Combustion and Flame, 209:278–290, 2019.

[106] T. Abbasi, H.J. Pasman, and S.A. Abbasi. A scheme for the classification of explosions
in the chemical process industry. Journal of Hazardous Materials, 174(1-3):270–280,
February 2010.

Cléante Langrée UVCE Simulation February 8, 2022

148 BIBLIOGRAPHY

[107] D. Bradley, M. Lawes, and M.S. Mansour. The problems of the turbulent burning
velocity. Flow, Turbulence and Combustion, 87:191–204, 2011.

[108] I. Ahmed and N. Swaminathan. Simulation of Spherically Expanding Turbulent Premixed

Flames. Combustion Science and Technology, 185(10), 10 2013.

[109] D. Bradley, M.Z. Haq, R.A. Hicks, T. Kitagawa, M. Lawes, C.G.W. Sheppard, and
R. Woolley. Turbulent burning velocity, burned gas distribution, and associated flame
surface definition. Combustion and Flame, 133(4), 6 2003.

[110] G. Damköhler. Zeilschrift für Elektrochemie und angewandte physikalische Chemie.

Der Einfluß der Turbulenz auf die Flammengeschwindigkeit, 46(11):601–626, 1940.

[111] N. Peters. The turbulent burning velocity for large-scale and small-scale turbulence.
Journal of Fluid Mechanics, 384, 4 1999.

[112] S.A. Filatyev, J.F. Driscoll, C.D. Carter, and J.M. Donbar. Measured properties of
turbulent premixed flames for model assessment, including burning velocities, stretch
rates, and surface densities. Combustion and Flame, 141(1-2), 4 2005.

[113] G. Nivarti and S. Cant. Direct Numerical Simulation of the bending effect in turbulent
premixed flames. Proceedings of the Combustion Institute, 36(2), 2017.

[114] A.E. Dahoe, T. Skjold, D.J.E.M. Roekaerts, H.J. Pasman, R.K. Eckhoff, K. Hanjalic,
and M. Donze. 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. Flow, Turbulence and Combustion, 91(2), 9
2013.

[115] L. Paris. Évaluation des effets d’une explosion de gaz à l’air libre. base documentaire :
TIB157DUO.(ref. article : se5062), 2009.

Cléante Langrée University of Rouen - CORIA February 8, 2022