Numerical simulation of deflagration-to-detonation transition in large confined volumes

Accepted Manuscript Numerical simulation of deflagration-to-detonation transition in large confined volumes Josef Hasslberger , Lorenz R. Boeck , Thomas Sattelmayer PII:

S0950-4230(14)00211-3

DOI:

10.1016/j.jlp.2014.11.018

Reference:

JLPP 2864

To appear in:

Journal of Loss Prevention in the Process Industries

Received Date: 12 September 2014 Revised Date:

21 November 2014

Accepted Date: 26 November 2014

Please cite this article as: Hasslberger, J., Boeck, L.R., Sattelmayer, T., Numerical simulation of deflagration-to-detonation transition in large confined volumes, Journal of Loss Prevention in the Process Industries (2014), doi: 10.1016/j.jlp.2014.11.018. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights

ACCEPTED MANUSCRIPT

EP

TE D

M AN U

SC

RI PT

CFD simulations of hydrogen-air explosions in large confined volumes are presented. Adaptive mesh refinement is a key technique to minimize computational costs. The model is successfully validated by large-scale experiments in the RUT facility. Safety-relevant lean mixtures close to the lower detonation limit are considered. The influence of mixture inhomogeneity on flame acceleration and DDT is discussed.

AC C

    

ACCEPTED MANUSCRIPT

RI PT

Numerical simulation of deflagration-to-detonation transition in large confined volumes Josef Hasslbergera,∗, Lorenz R. Boecka , Thomas Sattelmayera

f¨ ur Thermodynamik, Technische Universit¨ at M¨ unchen, Boltzmannstrasse 15, 85748 Garching, Germany

SC

a Lehrstuhl

Abstract

M AN U

A methodology for the computationally efficient CFD simulation of hydrogenair explosions (including transition to detonation) in large volumes is presented. The model is validated by means of the largest ever conducted indoor DDT experiments in the RUT facility. A combination of models is proposed with a particular focus on the influence of flame-instabilities, especially of thermaldiffusive nature, which are crucial for very lean mixtures. Excellent agreement is achieved in terms of flame acceleration. The quality of DDT predictions

TE D

itself depends on the underlying mechanism. Whereas DDT by shock-focusing is successfully simulated on under-resolved meshes, DDT by local explosions in the vicinity of the turbulent flame brush remains a challenge. Adaptive mesh refinement therefore emerges as a key technique to resolve more of the essential phenomena at reasonable computational costs affordable by industry. Finally,

EP

a generic case demonstrates the influence of mixture inhomogeneity, which can promote flame acceleration and ultimately DDT.

AC C

Keywords: Hydrogen safety, gas explosion, flame acceleration, DDT, RUT facility, under-resolved CFD

∗ Corresponding

author Email address: [email protected] (Josef Hasslberger)

Preprint submitted to Journal of Loss Prevention in the Process IndustriesNovember 19, 2014

ACCEPTED MANUSCRIPT

1. Introduction

RI PT

In case of a core meltdown in water-cooled nuclear power plants, a considerable amount of hydrogen can be released. Together with the ambient air in the

reactor containment, a combustible mixture can be formed in a wide stoichiome5

try spectrum. According to Kuczera [1], several hundred kilograms of hydrogen

accident, eventually leading to severe explosions.

SC

were produced in each of the affected reactor blocks in the Fukushima-Daiichi

Homogeneous mixtures have intensely been studied in the past, experimentally as well as numerically. However, mixtures with concentration gradients represent more realistic scenarios of hydrogen distribution within a reactor con-

M AN U

10

tainment [2]. Due to gravity and the different densities of hydrogen and air (approximately one order of magnitude), hydrogen tends to accumulate below the roof of facilities. Following the scaling law of diffusion, the homogenizing effect of molecular diffusion is comparably slow in large volumes. 15

Once the mixture is ignited (e.g. by an electric spark or electrostatic dis-

TE D

charge), the combustion process usually starts as a slow laminar deflagration, but can quickly transform into a fast turbulent deflagration through several accelerating mechanisms (e.g. instabilities and turbulence generation by obstacles) and finally undergo transition to detonation. Especially the abrupt 20

transition between the two different combustion regimes, deflagration and det-

EP

onation, is of interest since a detonation poses a major threat to the integrity of the containment.

Consequently, the goal in nuclear reactor risk assessment is to develop tools

AC C

for the prediction of the hazardous potential in such accident situations. The

25

analysis chain consists of 1. reactor core behavior (hydrogen release rate), 2. hydrogen dispersion (mixture composition), 3. combustion process (pressure loads), 4. structural response of the containment.

2

ACCEPTED MANUSCRIPT

30

The output given in brackets is the necessary input for the subsequent step. The

RI PT

present work is supposed to make a contribution to the third step. In principle, the methodology is not limited to nuclear safety, but is generally applicable to other fields like process industries.

Regarding the numerical modeling of reactive flows, a large variety of codes 35

exists that are optimized for particular combustion regimes depending on the un-

SC

derlying physics. Whereas deflagrations are subsonic combustion waves driven by diffusive and turbulent exchange of heat and species, detonations are super-

sonic combustion waves mainly driven by gas-dynamic effects and auto-ignition.

40

M AN U

Since all different regimes are relevant to explosive combustion, a robust and efficient model with a wide range of validity has to be established. The fundamentally different mechanisms behind deflagrations and detonations complicate the development of consistent approaches. Therefore, the state-of-the-art in large-scale explosion modeling [3] is to evaluate empirical transition criteria providing information about the worst-case combustion regime that can be ex45

pected. Optimized solvers can then be used for consequence analysis. For the

TE D

transition between slow and fast deflagrations, the so-called σ-criterion σ=

ρu > 3.75 ρb

(1)

has to be satisfied, i.e. a sufficiently high expansion ratio (ratio of densities of reactants ρu and products ρb ) of the hydrogen-air mixture must be given.

50

EP

Contrary to slow deflagrations, the propagation velocity of the flame front with respect to an external observer surpasses the speed of sound of the reactants

AC C

for fast deflagrations. The second criterion, defined by means of the detonation cell width λ,

L > 7λ

(2)

describes a necessary requirement for the onset of detonation from a fast deflagration. Criticism about this approach is related to the characteristic geometric

55

length scale L which is difficult to define for complex facilities as well as the stoichiometry-dependent detonation cell width that is approximated from experiments with some uncertainty. Furthermore, the effect of potential mixture 3

ACCEPTED MANUSCRIPT

inhomogeneity is not incorporated. Such empirical transition criteria are mainly

60

RI PT

deduced from experiments on laboratory scale and do not account for all of the underlying initial and boundary conditions. In contrast, CFD simulations are able to consider the complex three-dimensional geometry and the local fuel-

oxidator ratio of the particular system. A group of combustion and nuclear

experts [3] stated that it would be a significant improvement to simulate flame

65

SC

acceleration, DDT and detonation propagation within one solver framework,

not being reliant on empirical transition criteria and a combination of different codes.

M AN U

Common approaches for DDT or detonation simulations are based on onestep Arrhenius kinetics, require high spatial resolution and are therefore not qualified for real-world scenarios and particularly not for nuclear power plants 70

for which very large domains (approximately 50,000 m3 ) have to be discretized. Estimating the increase in computational power by means of Moore’s law, it is clear that even in the next few decades it will not be possible to resolve all smallscale phenomena for large-scale applications. Thus, there is an urgent need for

75

TE D

efficient simulation techniques that perform well on comparably coarse meshes. Encouraging results on under-resolved meshes have recently been achieved for DDT simulations (Middha & Hansen [4] and Gaathaug et al. [5]) and pure detonation simulations (Zbikowski et al. [6]). Earlier (two-dimensional) deflagration

EP

simulations of the RUT facility (Sect. 3) have been conducted by Breitung & Kotchourko [7], Smiljanovski et al. [8] and Efimenko & Dorofeev [9]. CFD com80

bustion analysis of a full-scale AREVA EPR reactor containment is presented

AC C

by Dimmelmeier et al. [10]. Latter study is, however, limited to deflagrations and not considering DDT events. Since for safety analysis, it is primarily interesting to know key characteristics

such as flame speeds and pressure loads on the containing structure, we abstain

85

from resolving all small-scale mechanisms. The focus of this study is placed on the global propagation behavior rather than the micro-structure of DDT.

4

ACCEPTED MANUSCRIPT

2. Numerical approach

RI PT

2.1. Solver basis The main characteristics of the finite-volume solver, implemented in the 90

open-source CFD package OpenFOAM, are briefly described. The code is writ-

ten in C++ and is able to handle geometrically complex unstructured meshes. Due to the mixed parabolic-hyperbolic nature of the explosion problems stud-

SC

ied (slow quasi-laminar deflagration shortly after weak ignition; fast detonation

propagation after DDT), the numerical methodology initially developed by Et95

tner et al. [11] is built on the unsteady compressible Navier-Stokes equations

M AN U

(in a Favre-averaged sense [12]) and a transport equation for the total internal energy. The system of equations is closed by the ideal gas law. Turbulence is included by the k-ω-SST model [13] in order to combine the advantages of the k-ω model in near-wall regions and the k-ε model in free-stream regions. An explicit 100

and density-based formulation ensures accurate reproduction of gas-dynamic effects, i.e. shocks, which is especially important since auto-ignition phenomena have to be considered. Convective fluxes are computed by the second order

TE D

accurate Harten-Lax-van Leer with contact (HLLC) scheme [14] with linear reconstruction of the left and right cell face values and multi-dimensional gradi105

ent limiters. The time step is adjusted dynamically according to the CourantFriedrichs-Lewy (CFL) condition calculated with the maximum characteristic

AC C

EP

wave velocity (magnitude of convective flow velocity ui plus speed of sound of √ an ideal gas κRT ) in the domain:   p |˜ ui | + κRT˜ ∆t CFL = . (3) ∆x A threshold value of CFL = 0.3 guarantees stability. Temperature-dependent

110

gas properties are obtained from the Chemkin [15] database, molecular transport coefficients are calculated by Sutherland’s formula [16]. 2.2. Combustion model The comparably efficient reaction model consists of a flamelet deflagration model and a two-step auto-ignition model. Coupling is realized on the basis of 5

ACCEPTED MANUSCRIPT

a transport equation ∂ ∂ ∂ (¯ ρc˜) + (¯ ρc˜u ˜j ) − ∂t ∂xj ∂xj

  ∂˜ c 1−e c ρ¯Deff = ρ¯u GΞSL |∇˜ c| + ρ¯ H(e τ − 1) (4) ∂xj ∆t

RI PT

115

for the reaction progress variable c, where c = 0 indicates the unburned mixture and c = 1 denotes the completely burned state. The effective mass diffusivity 1 µT ScT ρ

(5)

SC

Deff = D + DT = D +

is composed of a laminar part D (negligible in highly turbulent flows) as well as a turbulent part DT which can be closed by defining a turbulent Schmidt number ScT = 1.0 and by using the eddy viscosity µT provided by the turbulence model.

M AN U

120

The frequently used gradient formulation defines the first source term for flame propagation by deflagration. Unlike the common approach for large-scale simulations, to directly close the TFC (turbulent flame speed closure) deflagration model via empirical correlations for the turbulent burning velocity ST , an 125

additional transport equation is solved for the sub-grid flame wrinkling factor

ST AT = . SL A⊥

TE D

(following Damk¨ ohler’s ansatz)

Ξ=

(6)

The transport equation’s right hand side

EP

ρΞGΞ − ρ(Ξ − 1)RΞ

(7)

includes the generation rate GΞ =

0.28 τη

(8)

Ξeq Ξeq − 1

(9)

AC C

and the reduction rate

130

RΞ = GΞ

of flame surface area A, respectively. The Kolmogorov time scale of turbulence

is estimated as

r τη =

νu ε

(10)

through dimensional analysis whereas the constant 0.28 was determined by DNS simulations [17]. This approach incorporates non-equilibrium effects (in terms 6

ACCEPTED MANUSCRIPT

of flame wrinkling) which are likely to arise for rapid changes of state as in accelerated flames. Turbulence-chemistry interaction is therefore an inherent

RI PT

135

part of the combustion model itself. Ξ is related to the better-known flame area density Σ by Σ = Ξ|∇˜ c|. The laminar burning velocity SL is a priori

known from a polynomial representation, complemented by a pressure- and temperature-correction. G incorporates quenching effects at high turbulence intensities [18, 19, 20]: G=

SC

140

  1 1   εcr  σ  + erfc − √ ln 2 ε 2 2σ

(11)

M AN U

with erfc indicating the complementary error function and σ representing the standard deviation of the turbulent dissipation rate ε. Modeling by   lT σ = 0.28 ln lη

(12)

uses the integral turbulent length scale lT and the Kolmogorov length scale 

TE D

lη =

ν3 ε

0.25

.

(13)

The critical dissipation rate

2 εcr = 15νgcr

depends on the critical velocity gradient gcr which is approximated by

EP

145

(14)

gcr =

SL2 . au

(15)

Differing from Weller’s original formulation [17], the equilibrium wrinkling

AC C

Ξeq (the flame wrinkling when source and sink term balance) calculation is replaced by the recently published correlation of Dinkelacker et al. [21]. Earlier simulations showed that the original model clearly underestimates the flame

150

acceleration for lean mixtures in RUT experiments [22]. Insufficient modeling of flame-instabilities was identified as the main deficit as only turbulence-induced flame wrinkling was included. As proposed by e.g. Bradley et al. [23], the Lewis number of the premixture Le = 7

a D

(16)

ACCEPTED MANUSCRIPT

can be used to describe the strong influence of flame-instabilities, particularly of thermal-diffusive nature with subsequent secondary (hydrodynamic-) instabili-

RI PT

155

ties, on effective burning velocity for lean mixtures. It is given by the thermal diffusivity a divided by the diffusion coefficient D of the deficient reactant (hy-

drogen in case of lean mixtures). Besides standard turbulence input parameters (rms velocity fluctuation u0 and integral turbulent length scale lT ), Dinkelacker’s model Ξeq

0.46 =1+ ReT 0.25 Le



u0 SL

with the turbulent Reynolds number

p p0

0.2

M AN U

u0 lT , νu

0.3 

SC

160

ReT =

(17)

(18)

additionally contains a dependency on Le and pressure p. The Lewis number is a priori calculated by the chemical package Cantera [24] using a mixture-based approach and ranges approximately between 0.36 (12.5 % of hydrogen) and 0.38 165

(15.5 % of hydrogen) for the simulations performed. In order to handle inho-

TE D

mogeneous mixtures, Le is prescribed in the CFD code as a fuel concentration dependent polynomial. Dinkelacker et al. found that the Lewis number effect is especially important for increased pressure. In the context of explosion problems which are massively accompanied by pressure increase, this seems to be an 170

important finding. Moreover, even at high turbulence intensities, a significant

EP

influence of the Lewis number has been observed for lean premixed combustion. Other ways of instability modeling, e.g. via Markstein number (generally quantifying the effect of stretch rate on burning velocity) based correlations, might

AC C

be a reasonable alternative.

175

Since not only deflagrations have to be described, but also the transition to

detonation, the combustion model is extended for auto-ignition effects via the second source term (similar to approaches by Colin et al. [25] or Gaathaug et al. [5]). To do so, another transport equation   ∂ ∂ ∂ ∂e τ ρ (ρe τ) + (ρe τu ej ) − ρDeff = ∂t ∂xj ∂xj ∂xj tign

8

(19)

ACCEPTED MANUSCRIPT

is introduced in which the normalized quantity τ = t/tign (T, p, fH ) = y/ycrit compares the current simulation time t with the local auto-ignition delay time

RI PT

180

tign (T, p, fH ) depending on temperature T , pressure p and hydrogen mixture

fraction fH . Alternatively, τ can be seen as the ratio between the mass fraction of a fictitious intermediate species y and its critical mass fraction ycrit at which

auto-ignition occurs. The Heaviside function H activates the second source term in Eq. 19 when τ˜ = 1 and, as a result, c˜ = 1 is locally reached in the following

SC

185

time step t + ∆t. Ignition delay times are provided via multi-dimensional tables obtained from zero-dimensional isochoric explosion calculations [24] with the

M AN U

detailed reaction mechanism of O’Conaire et al. [26] during pre-processing. In doing so, the time until the inflection point of the temperature curve is measured. 190

Stiff hydrogen-oxygen chemistry (meaning that the chemical time scales are smaller than the characteristic flow time scales), when using Arrhenius source terms during runtime instead of pre-tabulation, is consequently not a limiting factor for the CFD calculations performed.

In this context, the code is optimized for massive parallelization as combustion usually takes place only in a small portion of the domain in the explosion

TE D

195

problems under investigation. It would cause heavy load imbalances if chemistry was calculated during runtime. Parallelization in OpenFOAM is realized via the domain decomposition technique, communication between processors is

200

EP

implemented by the OpenMPI library.

Apart from initial and boundary conditions, the model requires no tuning (particularly of the reaction rate) and is therefore valid for a wide range of dif-

AC C

ferent configurations. Furthermore, it is not necessary to incorporate empirical combustion regime transition criteria as e.g. in [4] where the likelihood of DDT is evaluated in terms of a parameter proportional to the spatial pressure gradient

205

across the flame front. Both is essential as we want to investigate the model’s predictive capabilities for general large-scale accident scenarios involving hydrogen as the fuel and air as the oxidizer. In case of inhomogeneous mixtures, an

9

ACCEPTED MANUSCRIPT

additional transport equation ∂ f˜H ρ¯Deff ∂xj

! =0

(20)

RI PT

 ∂  ˜ ∂  ˜  ∂ ρ¯fH + ρ¯fH u ˜j − ∂t ∂xj ∂xj

with no source term is solved for the hydrogen mixture fraction. fH can be 210

interpreted as the mass fraction of hydrogen that would be present if no chemical

reaction occurred. Following the safety analysis chain mentioned in Sect. 1,

SC

realistic multi-dimensional mixture fields (e.g. gained from CFD dispersion simulations) can be imported as initial conditions for combustion simulations.

M AN U

2.3. Adaptive mesh refinement

To keep the overall computational costs at a reasonable level, we apply adap-

215

tive mesh refinement (AMR). Manually refining the mesh is not an option due to the highly unsteady nature of explosions. The technique facilitates having high spatial resolution only in regions where small-scale phenomena are crucial while the base mesh can be chosen comparably coarse. Primarily the flame front 220

requires finer resolution as it dominates the propagation behavior. Obtaining

TE D

realistic reaction rates is of highest importance in making accurate predictions of the flame acceleration process in closed systems because of the feedback mechanism via the expansion of combustion products. Since the position of the flame front is known from the reaction progress field, the refinement region is specified as

EP

225

0.001 < c˜ < 0.999,

(21)

implying unrefinement in the completely burned mixture. In this context, it is

AC C

important to note that, because of the gradient formulation of the deflagration source term, the reaction rate does not scale with the thickness of the reaction zone which could be altered by AMR.

230

Additionally, the mesh is refined in the unburned part of the domain Ωu

(˜ c < 0.001) where the normalized velocity gradient |∇|˜ ui || > 0.1. maxΩu (|∇|˜ ui ||)

10

(22)

ACCEPTED MANUSCRIPT

Normalization allows setting reasonable relative threshold values (10 % here)

RI PT

for per se unbounded quantities. The combined refinement criterion includes relevant physics (which stems from or generates velocity gradients) ahead of 235

the flame as is evident from Fig. 3. Enhanced turbulence production in the obstacles’ wake and precursor shock waves – two further phenomena of major importance in turbulent reactive flows with gas-dynamic effects.

SC

The underlying mechanism of the algorithm is isotropic cell division. Accord-

ingly, each mother cell fulfilling the refinement rule is split into eight daughter 240

cells for three-dimensional simulations. In principle, multiple refinement op-

M AN U

erations are possible and are therefore limited by a maximum refinement level related to the affordable computational effort. Good mesh quality is maintained through obeying the 2:1 criterion, i.e. cell faces can have a maximum of four neighbors in three-dimensional simulations. For the present calculations, the al245

gorithm is executed every ten time steps due to the profoundly unsteady nature of the problem. There might be potential to further reduce the computational cost by suspending AMR for a larger number of time steps. Contrary to the

TE D

applied feature-driven strategy, an error-driven strategy is unfeasible due to the lack of analytical solutions as well as residuals in the explicit solver framework.

250

3. Experiment and computational setup

EP

For the validation of the computational model, the largest ever conducted indoor DDT experiments in the RUT facility at Kurchatov Institute, Moscow region [27], have been chosen. So far, the methodology described in Sect. 2 has

AC C

successfully been validated for laboratory-scale DDT experiments undertaken by

255

Vollmer et al. [28]. Doing three-dimensional instead of two-dimensional simulations [11] had a relevant influence on the flame acceleration process [29]. Despite significantly higher computational costs, former strategy is the preferred choice to represent the complex geometry of the RUT facility (Fig. 1) with an approximate inner volume of 424 m3 (exclusive of the obstacles). The first channel

260

(length ≈ 35 m; height 2.3 m; depth 2.5 m) causes effective flame acceleration

11

ACCEPTED MANUSCRIPT

Main direction of flame propagation 2.3 m

RI PT

Ignition point

4m ≈ 10.5 m

SC

≈ 35 m 2.5 m

M AN U

Figure 1: RUT configuration 22 (above) and configuration 16 (below); From left to right: obstructed first channel, canyon, curved second channel

via turbulence generation in the wake of obstacles etc. and is followed by a sudden jump in cross-section. A curved ramp (second channel) succeeds the canyon (length ≈ 10.5 m; height 4 m) after a sudden decrease in cross-section. It is therefore possible to examine various influences on flame acceleration and DDT, such as obstructions, transversal expansion, shock focusing, shock-flame

TE D

265

interaction etc. For better interpretation, beginning and end of the canyon are indicated by straight lines in the subsequent plots describing the flame acceleration etc.

270

EP

Two different cases are investigated regarding fuel concentration and obstacle configuration (Tab. 1 and Fig. 1). Geometrical details (e.g. chamber bevels or obstacle gaps due to constructional reasons), which are expected to be important

AC C

in terms of flame acceleration, are taken from Moser [30] and Breitung et al. [31]. Measurement infrastructure and further (turbulence-inducing) channel installations have to be neglected due to their small dimensions. Depending on

275

the setup, the observed mechanism behind DDT differs. Klein & Rehm [32] generally distinguish between DDT by shock focusing e.g. near corners (strong solution or mode A) and DDT originating from hot spots and ultimately local explosions in the vicinity of the turbulent flame brush (weak solution or mode B).

12

ACCEPTED MANUSCRIPT

Table 1: Investigated RUT configurations

RUT 16

DDT mode

A

B

Hydrogen concentration (vol.)

14 %

12.5 %

Obstacle blockage ratio

60 %

30 %

Obstacle count

6 (5 m spacing)

12 (2.5 m spacing)

Base mesh

1,284,792 cells

1,277,155 cells

SC

RI PT

RUT 22

Latter can be triggered by shock-flame interaction or the explosion of enclosed unburned pockets for instance. While in case 22, DDT happens through shock-

M AN U

280

obstacle interaction in the first channel, DDT is likely to function via mechanism B in the canyon for case 16.

The average cell size of the base mesh, ≈ 7 cm, is apparently much larger than best practice guidelines (if available) would suggest for detailed DDT simu285

lations. It is by far not possible to resolve the half-reaction length or other characteristic scales. Regardless, the method has to prove robust for under-resolved

TE D

meshes in order to be applicable to large-scale scenarios. The significance of sub-grid models for turbulence and flame wrinkling becomes obvious. Much finer meshes are not reasonable, anyway, since the application to realistic scales 290

would then demand access for supercomputers which are not broadly available

EP

to industry.

Initial values of the turbulent quantities k and ω are chosen comparably low, and in any case, are not crucial as the main stage of flame acceleration

AC C

is dominated by the turbulence production behind obstacles. Basically, some

295

initial turbulence must have been present since mixing fans were used to assure a nearly homogeneous mixture at the time of ignition, yielding a mixing accuracy of 0.1 % in terms of the hydrogen content [27]. The turbulence level was not measured however. No-slip adiabatic walls are imposed at all boundaries including the vent opening at the end of the second channel which had no effect

300

on the transition to detonation according to Dorofeev et al. [27]. Heat losses

13

ACCEPTED MANUSCRIPT

2000

Exp. (14 %)

CJ theory (14 %)

Sim. (14 %)

40

Sim. (14 % - No AMR)

1600

30 20

Sim. (14 %)

Sim. (14 % - No AMR) 1200 800

400

10 0

0

0

0.05

0.1

0.15

0

10

20

30

40

50

Distance [m]

SC

Time [s]

RI PT

50

Velocity [m/s]

Distance [m]

60

Figure 2: RUT configuration 22 flame-tip acceleration

M AN U

to the walls are assumed to be negligible for the fast explosion process. Initial temperature and pressure are set to 293 K and 1 atm, respectively. The delicate issue of accurately modeling weak spark ignition is avoided by starting at a predefined flame front, i.e. by patching a sphere of 1 m radius around the 305

ignition point (beginning of the first channel) to the fully burned state.

TE D

4. Results and discussion 4.1. RUT configuration 22

Fig. 2 shows the numerical results together with experimental data from Dorofeev et al. [27]. The time shifting of experimental data (− 0.29 s) enables 310

better comparison and does not distort the flame acceleration behavior which

EP

is of interest here. Excellent agreement can be observed between simulation and experiment. The flame gradually accelerates until DDT occurs at the last

AC C

obstacle. DDT is recognizable by the kink in the flame-tip trajectory in the distance-time diagram as well as the sudden jump in velocity in the velocity-

315

distance diagram. After DDT, the detonation front steadily propagates at the theoretically predicted Chapman-Jouguet (CJ) speed. Latter diagram provides further insight into the effect of obstacles. Instead of continuous acceleration, the flame-tip velocity oscillates which is due to the jet behavior when the flow passes each obstacle.

14

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 3: Level 1 adaptive mesh refinement for regions of high turbulence production (left; 5th obstacle; 100 ms after ignition) and precursor shock waves (right; 6th obstacle; 104 ms

M AN U

after ignition) ahead of the flame; Pressure field p is given in Pa; The white-colored c˜ = 0.5 contour line indicates the flame front.

320

To demonstrate its benefit, AMR is switched off in a second computation of the same case. It can be seen that the flame acceleration in the obstructed channel is more or less the same. However, no DDT takes places without AMR, presumably because of insufficient reproduction of the strength of shocks. Whereas

325

TE D

the deflagration decelerates when passing the canyon (momentum loss due to transversal expansion), the detonation continues to propagate at constant velocity independent of changes in cross-section. AMR was essential in this case to capture the critical conditions for the onset of detonation. The effect of the refinement rule described in Sect. 2.3 is highlighted in Fig. 3. The mesh is

330

EP

locally refined in the vicinity of the flame front (white-colored c˜ = 0.5 contour line), known as the turbulent flame brush. Refinement ahead of the flame includes regions of enhanced turbulence production, e.g. in the wake of obstacles

AC C

(left picture), and precursor shock waves (right picture). The leading pressure wave in the left picture is not yet strong enough to be refined. The burned part of the domain is regarded as minor important and thus, in order to save

335

computational resources, not refined. AMR comes at extra costs of course. Fig. 4 gives an impression of the

relative increase in computational effort which scales with the overall number of cells in the domain. To reflect the dynamic behavior, the current number

15

ACCEPTED MANUSCRIPT

1.6E+06

RI PT

Number of cells [-]

1.7E+06

1.5E+06 1.4E+06

Sim. (14 %)

1.3E+06 1.2E+06

0

10

20

30

SC

Sim. (14 % - No AMR)

40

50

60

M AN U

Distance [m]

Figure 4: Overall number of computational cells evolution

of cells is plotted versus the flame-tip position. Resulting from the refinement 340

rule, the interaction with obstacles in the first channel as well as entering and leaving the canyon manifest in the number of cells. Note the general decrease after transition to stable detonation propagation, for which the mixture ahead

TE D

of the front is predominantly unaffected. A detailed look at DDT is presented in Fig. 5. In the first picture (103 ms), 345

from left-to-right propagating precursor shock waves and the flame are clearly separated, indicating a deflagration. At 104 ms after ignition, the leading shocks

EP

are reflected off the sixth obstacle. Obviously, the reflected shock is not strong enough to cause auto-ignition directly at the obstacle. The local explosion rather seems to originate from the intersection point of shock waves reflected from the obstacle and the bottom of the channel (105 ms). Equally to the

AC C

350

overdriven forward propagating wave, a retonation wave emerges, propagating into the burned part of the domain (106 ms). The overdriven state, characterized by propagation velocities higher than the CJ velocity, can be recognized in the velocity-distance diagram in Fig. 2 as well. At this point in time, the

355

strong forward propagating wave and the flame are already coupled. The overdriven wave quickly decays to the steady CJ solution, transition to detonation

16

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 5: DDT by shock-obstacle interaction in detail (times 103 ms, 104 ms, 105 ms, 106 ms

TE D

107 ms after ignition); Pressure field p is given in Pa (note change of scale); The white-colored c˜ = 0.5 contour line indicates the flame front.

is finished. In the final picture (107 ms), detonation and retonation wave are stably propagating in opposite directions. The mechanism in the simulation

360

[27].

EP

corresponds very well to the DDT mode conclusions drawn by Dorofeev et al.

AC C

4.2. RUT configuration 16 An overview of the results is given in Fig. 6. The time shifting for measured

data is − 0.34 s in this case. In the original publication [27], the mechanism

behind DDT is analyzed by means of pressure recordings and two-dimensional

365

flame shape reconstructions from photodiode measurements. The authors conclude that the onset of detonation did not directly result from the reflection of the leading shock at the canyon end wall but from local explosions in the

17

ACCEPTED MANUSCRIPT

60

2000 Exp. (12.5 %) Sim. (12.5 %) Sim. (13.5 %)

30 20

Sim. (13.5 %)

1200 800

400

10 0

0

0

0.05

0.1

0.15

0.2

0

10

20

30

40

50

Distance [m]

SC

Time [s]

RI PT

40

Sim. (12.5 %) 1600 Velocity [m/s]

Distance [m]

50

Figure 6: RUT configuration 16 flame-tip acceleration

M AN U

close vicinity of the flame brush, i.e. DDT mode B. This particular experiment represents the lowest ever measured detonation limit (12.5 %) of hydrogen-air 370

mixtures. For laboratory-scale facilities, the lower detonation limit usually is several percent of hydrogen higher (minimum of 15 %) [27].

Again, the simulation matches the experiment very well in terms of flame acceleration. However, DDT does not occur for the same hydrogen concentration as in the experiment (12.5 %). Unfortunately, experimental data points for the canyon section (Fig. 6) are missing in the original publication. To en-

TE D

375

force DDT for this geometrical configuration, the molar hydrogen concentration is increased by 1 %. As a consequence, DDT seems to originate from a local explosion nearby the end wall of the canyon - similar to, but not exactly by the

380

EP

same mechanism as in the experiment. The clearly stronger flame acceleration underlines the strong sensitivity on hydrogen concentration close to the lower detonation limit which is well-known from experiments.

AC C

DDT mode B is obviously more ambitious to simulate than DDT mode A.

Klein & Rehm [32] have drawn similar conclusions. The reason is probably connected to several flame-instabilities [33] that are crucial but not adequately

385

treated (modeled or resolved). Particularly the Richtmyer-Meshkov instability could play a key role for DDT mode B via shock-flame interaction [34]. If pressure gradients (due to shocks) and density gradients (due to flames) are not aligned with each other, baroclinic vorticity is generated that wrinkles the flame

18

ACCEPTED MANUSCRIPT

80 Sim. (12.5 %)

Pressure [bar]

40 20 0

0.1

0.15

0.2

0.25

M AN U

Time [s]

SC

RI PT

Sim. (13.5 %)

60

Figure 7: RUT configuration 16 pressure at the canyon rear wall

and suddenly leads to enhanced heat release. Secondary instabilities like the 390

Kelvin-Helmholtz instability follow [35]. Another important mechanism might be the turbulent mixing of hot products and cold reactants, eventually causing explosions of enclosed unburned pockets. Relevant instabilities are usually act-

TE D

ing on very small scales compared to the length scales of the geometry which explains the difficulties in doing under-resolved CFD simulations. Due to the 395

lack of appropriate sub-grid models, there is not really an alternative to trying to resolve such phenomena. AMR is certainly a step in the right direction.

EP

Numerical pressure recordings impressively demonstrate the hazards associated with DDT. The probing point is located at the canyon rear wall, 1.5 m above the canyon bottom, close to the origin of the local explosion in the 13.5 % case. The peak pressure is at about 10 bar for the case without DDT in

AC C

400

the simulation (12.5 %), contrary to almost 80 bar for the DDT case (13.5 %). According to Boeck et al. [36], 10 bar roughly represents a critical pressure level for the onset of detonation. 4.3. Generic RUT configuration with a concentration gradient

405

Mixture inhomogeneity - an issue that is rarely considered in explosion research. A generic study is therefore presented on the influence of concentration 19

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 8: Generic RUT configuration molar hydrogen concentration distribution (120 ms after

M AN U

ignition)

gradients. The geometry is the same as for RUT configuration 22. For the sake of simplicity, we assume a one-dimensional vertical concentration gradient (i.e. perpendicular to the main direction of flame propagation) in this work. Basi410

cally, there is no impediment in importing realistic multi-dimensional mixture fields from CFD dispersion simulations. A linear vertical distribution yields 15.5 % of hydrogen at the top and 10.5 % at the bottom of the obstructed channel.

TE D

Accordingly, the gradient is given by the difference of 5 % along the channel height of 2.3 m. The average concentration in the obstructed channel is ≈ 13 415

% whereas it is ≈ 11.6 % in the whole domain because of the very lean region in the canyon. A snapshot of the hydrogen distribution, 120 ms after ignition,

EP

can be seen in Fig. 8. Large-scale mixing behind obstacles can be identified. It should be kept in mind that the mixture ahead of the flame (entering the picture from the left) is not equal to the unperturbed initial field - contrary to homogeneous mixtures.

AC C

420

Despite the same average fuel concentration in the first channel, the flame-

tip trajectory in Fig. 9 reveals transition to detonation for the gradient case but not for the homogeneous case. It seems like the stronger flame acceleration is responsible. For comparison, the 14 % reference experiment is additionally

425

plotted. Although there is no direct validation in the RUT facility, the result generally agrees with experimental observations of Boeck et al. [37] on laboratory scale. In both studies, the conclusion is that inhomogeneous mixtures 20

ACCEPTED MANUSCRIPT

60

Exp. (14 %) Sim. (13 %)

RI PT

Distance [m]

50

Sim. (13 % - Gradient)

40 30 20

0

0

0.05

0.1

0.15

0.2

M AN U

Time [s]

SC

10

Figure 9: Generic RUT configuration flame-tip acceleration

can promote flame acceleration and ultimately DDT, especially in the very lean regime. 430

The reason for the diverging flame acceleration might be connected to a different evolution of the flame surface area. Since the (laminar) burning velocity

TE D

is higher at the top than at the bottom, the flame is likely to elongate more for the gradient case, thus relatively increasing its surface area. Fig. 10 shows the evolution of the overall macroscopic flame surface area depending on the flame435

tip position. In this context, macroscopic is understood as the c˜ = 0.5 contour

EP

surface which excludes the wrinkling of the flame on sub-grid level. A contourbased piecewise-linear reconstruction algorithm [38] is therefore implemented and evaluated during the runtime of the solver. The higher flame surface area

AC C

for the gradient case consequently implies a higher integral heat release rate 440

which in turn promotes the feedback mechanism behind flame acceleration via the expansion of combustion products.

5. Concluding remarks The methodology presented seems to be able to capture the essential physics behind flame acceleration and DDT in obstructed channels. A decisive improve-

21

ACCEPTED MANUSCRIPT

100

RI PT

80 Area [m2]

60 40 Sim. (13 %)

20

0

10

20

30

40

M AN U

Distance [m]

SC

Sim. (13 % - Gradient) 0

Figure 10: Macroscopic flame surface area evolution

445

ment in terms of flame acceleration was achieved by modifying Weller’s original deflagration model. The recently published, Lewis number based, flame wrinkling correlation of Dinkelacker et al. models the influence of flame-instabilities, especially of thermal-diffusive nature, which are crucial for very lean mixtures

450

TE D

(. 15 % with Lewis numbers far smaller than unity). In contrast to very successful simulations of DDT mode A (DDT by shockfocusing as in RUT configuration 22), DDT mode B (RUT configuration 16) remains a challenge. The reason is probably related to insufficiently modeled or

EP

resolved flame-instabilities like the Richtmyer-Meshkov instability or the explosion of enclosed unburned pockets in the vicinity of the flame front. 455

Adaptive mesh refinement is an appropriate numerical technique to resolve

AC C

more of the underlying phenomena at reasonable computational costs. A combined refinement criterion is applied to dynamically refine the turbulent flame brush as well as relevant physics ahead of the flame (primarily precursor shock waves and regions of enhanced turbulence generation).

460

Since homogeneous mixtures are generally not realistic for accident scenarios,

the effect of mixture inhomogeneity is discussed by means of a generic case. The conclusion that inhomogeneous mixtures can promote flame acceleration and

22

ACCEPTED MANUSCRIPT

ultimately DDT corresponds well to laboratory-scale experiments.

465

RI PT

In the future, we want to combine adaptive mesh refinement with massive parallelization on a high performance cluster in order to reproduce DDT mode B

and to come closer to the measured lower detonation limit of 12.5 % [27]. After

successful validation for the RUT facility, the long-term objective is to analyze

a full-scale Konvoi-type (PWR) reactor containment, similar to the work of

470

SC

Dimmelmeier et al. [10] but with a focus on DDT.

Acknowledgements

M AN U

The presented work is funded by the German Federal Ministry of Economic Affairs and Energy (BMWi) on the basis of a decision by the German Bundestag (project no. 1501338 and 1501425) which is gratefully acknowledged.

References 475

[1] B. Kuczera, Das schwere Tohoku-Seebeben in Japan und die Auswirkun-

TE D

gen auf das Kernkraftwerk Fukushima-Daiichi, International Journal for Nuclear Power 56 (2011) 4/5.

[2] OECD Group of Experts, State-of-the-art report on containment thermalhydraulics and hydrogen distribution, Tech. rep., OECD/NEA/CSNI/R 16 (1999).

EP

480

[3] W. Breitung, C. Chan, S. Dorofeev, A. Eder, B. Gerland, M. Heitsch, R. Klein, A. Malliakos, J. Shepherd, E. Studer, P. Thibault, Flame accel-

AC C

eration and deflagration–to–detonation transition in nuclear safety, Tech. rep., OECD/NEA/CSNI/R 7 (2000).

485

[4] P. Middha, O. R. Hansen, Predicting deflagration to detonation transition in hydrogen explosions, Process Safety Progress 27 (3) (2008) 192–204.

[5] A. V. Gaathaug, K. Vaagsaether, D. Bjerketvedt, Experimental and numerical investigation of DDT in hydrogen–air behind a single obstacle, International Journal of Hydrogen Energy 37 (22) (2012) 17606–17615. 23

ACCEPTED MANUSCRIPT

490

[6] M. Zbikowski, D. Makarov, V. Molkov, Numerical simulations of large–

Hazardous Materials 181 (1-3) (2010) 949–956.

RI PT

scale detonation tests in the RUT facility by the LES model, Journal of

[7] W. Breitung, A. Kotchourko, Numerische Simulation von turbulenten Wasserstoffverbrennungen bei schweren Kernreaktorunf¨allen, Tech. rep., FZK-Nachrichten 28 (1996).

SC

495

[8] V. Smiljanovski, V. Moser, R. Klein, A capturing–tracking hybrid scheme for deflagration discontinuities, Combustion Theory and Modelling 1 (2)

M AN U

(1997) 183–215.

[9] A. A. Efimenko, S. B. Dorofeev, CREBCOM code system for description of 500

gaseous combustion, Journal of Loss Prevention in the Process Industries 14 (6) (2001) 575–581.

[10] H. Dimmelmeier, J. Eyink, M.-A. Movahed, Computational validation of the EPR combustible gas control system, Nuclear Engineering and Design

505

TE D

249 (2012) 118–124.

[11] F. Ettner, K. G. Vollmer, T. Sattelmayer, Numerical simulation of the deflagration-to-detonation transition in inhomogeneous mixtures, Journal of Combustion 2014.

EP

[12] A. Favre, Equations des gaz turbulents compressibles, Journal de M´ecanique 4 (1965) 361–390. [13] F. R. Menter, Two-equation eddy-viscosity turbulence models for engineer-

AC C

510

ing applications, AIAA Journal 32 (1994) 1598–1605.

[14] E. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL-Riemann solver, Shock Waves 4 (1994) 25–34.

[15] R. J. Kee, F. M. Rupley, J. A. Miller, M. E. Coltrin, J. F. Grcar, E. Meeks, 515

H. K. Moffat, A. E. Lutz, G. Dixon-Lewis, M. D. Smooke, J. Warnatz, G. H. Evans, R. S. Larson, R. E. Mitchell, L. R. Petzold, W. C. Reynolds, 24

ACCEPTED MANUSCRIPT

M. Caracotsios, W. E. Stewart, P. Glarborg, C. Wang, O. Adigun, The

RI PT

Chemkin thermodynamic database (2000). [16] R. B. Bird, W. E. Stewart, E. N. Lightfoot, Transport phenomena, Wiley, 520

New York, 2001.

[17] H. Weller, The development of a new flame area combustion model using

SC

conditional averaging, Thermo–Fluids Section Report TF/9307, Imperial College London (1993).

[18] K. Bray, Complex Chemical Reaction Systems, Vol. 47 of Springer Series in Chemical Physics, Springer, 1986, Ch. Methods of Including Realistic

M AN U

525

Chemical Reaction Mechanisms in Turbulent Combustion Models, pp. 356– 375.

[19] V. Zimont, W. Polifke, M. Bettelini, W. Weisenstein, An efficient computational model for premixed turbulent combustion at high Reynolds numbers 530

based on a turbulent flame speed closure, Journal of Engineering for Gas

TE D

Turbines and Power 120 (1998) 526–532.

[20] W. Polifke, P. Flohr, M. Brandt, Modeling of inhomogeneously premixed combustion with an extended TFC model, Journal of Engineering for Gas Turbines and Power 124 (1) (2002) 58–65. [21] F. Dinkelacker, B. Manickam, S. Muppala, Modelling and simulation of lean

EP

535

premixed turbulent methane/hydrogen/air flames with an effective Lewis

AC C

number approach, Combustion and Flame 158 (9) (2011) 1742–1749. [22] J. Hasslberger, L. Boeck, T. Sattelmayer, Numerical simulation of deflagration-to-detonation transition in large confined volumes, in: Tenth

540

International Symposium on Hazards, Prevention, and Mitigation of Industrial Explosions, Bergen, Norway, 2014.

[23] D. Bradley, A. Lau, M. Lawes, Flame stretch rate as a determinant of turbulent burning velocity, Philosophical Transactions of the Royal Society

25

ACCEPTED MANUSCRIPT

of London. Series A: Physical and Engineering Sciences 338 (1650) (1992) 359–387.

RI PT

545

[24] D. Goodwin, N. Malaya, H. Moffat, R. Speth, Cantera: an object-oriented

software toolkit for chemical kinetics, thermodynamics and transport processes (2009).

550

SC

[25] O. Colin, A. Pires da Cruz, S. Jay, Detailed chemistry-based auto-ignition

model including low temperature phenomena applied to 3D engine calculations, Proceedings of the Combustion Institute 30 (2005) 2649–2656.

M AN U

[26] M. O’Conaire, H. J. Curran, J. M. Simmie, W. J. Pitz, C. K. Westbrook, A comprehensive modeling study of hydrogen oxidation, International Journal of Chemical Kinetics 36 (2004) 603–622. 555

[27] S. Dorofeev, V. Sidorov, A. Dvoinishnikov, W. Breitung, Deflagration to detonation transition in large confined volume of lean hydrogen–air mixtures, Combustion and Flame 104 (1-2) (1996) 95 – 110.

TE D

[28] K. Vollmer, F. Ettner, T. Sattelmayer, Deflagration-to-detonation transition in hydrogen-air mixtures with a concentration gradient, Combustion 560

Science and Technology 184 (2012) 1903–1915. [29] J. Hasslberger, F. Ettner, L. Boeck, T. Sattelmayer, 2D and 3D flame sur-

EP

face analysis of flame acceleration and deflagration–to–detonation transition in hydrogen–air mixtures with concentration gradients, in: 24th Inter-

AC C

national Colloquium on the Dynamics of Explosions and Reactive Systems, 565

Taipei, Taiwan, 2013.

[30] V. Moser, Simulation der Explosion magerer Wasserstoff–Luft–Gemische in grossskaligen Geometrien, Ph.D. thesis, Rheinisch–Westf¨alische Technische Hochschule Aachen (1997).

[31] W. Breitung, S. Dorofeev, A. Kotchourko, R. Redlinger, W. Scholtyssek, 570

A. Bentaib, J.-P. L’Heriteau, P. Pailhories, J. Eyink, M. Movahed, K.-G.

26

ACCEPTED MANUSCRIPT

Petzold, M. Heitsch, V. Alekseev, A. Denkevits, M. Kuznetsov, A. Efi-

RI PT

menko, M. Okun, T. Huld, D. Baraldi, Integral large scale experiments on hydrogen combustion for severe accident code validation-HYCOM, Nuclear Engineering and Design 235 (2005) 253 – 270. 575

[32] R. Klein, W. Rehm, Models and criteria for prediction of deflagration-todetonation transition (DDT) in hydrogen-air-steam systems under severe

SC

accident conditions, Forschungszentrum J¨ ulich, Zentralbibliothek, 2000.

[33] G. Ciccarelli, S. Dorofeev, Flame acceleration and transition to detonation

580

M AN U

in ducts, Progress in Energy and Combustion Science 34 (4) (2008) 499– 550.

[34] A. Khokhlov, E. Oran, G. Thomas, Numerical simulation of deflagrationto-detonation transition: the role of shock–flame interactions in turbulent flames, Combustion and Flame 117 (1) (1999) 323–339.

[35] M. Brouillette, The Richtmyer-Meshkov instability, Annual Review of Fluid Mechanics 34 (1) (2002) 445–468.

TE D

585

[36] L. Boeck, S. Msalmi, F. Koehler, J. Hasslberger, T. Sattelmayer, Criteria for DDT in hydrogen-air mixtures with concentration gradients, in: Tenth International Symposium on Hazards, Prevention, and Mitigation of In-

590

EP

dustrial Explosions, Bergen, Norway, 2014. [37] L. Boeck, J. Hasslberger, T. Sattelmayer, Flame acceleration in hydrogen-

AC C

air mixtures with concentration gradients, Submitted to Combustion Science and Technology.

[38] C. Kunkelmann, Numerical modeling and investigation of boiling phenomena, Ph.D. thesis, Technische Universit¨at Darmstadt (2011).

27