Two models for the computation of laminar flames in dust clouds

TL;DR

This paper compares two models for simulating premixed dust cloud combustion, introducing a phase-field approach for the second model and developing a finite volume scheme to ensure physical bounds.

Contribution

It presents a novel phase-field-like model for dust flame simulation and a finite volume scheme that maintains physical bounds, assessing their consistency.

Findings

Both models are consistent with physical bounds.

The phase-field model effectively tracks flame location.

Finite volume scheme ensures numerical stability.

Abstract

We address two models for the simulation of dust clouds premixed combustion: the first one consists in usual balance equations, to derive the second one, we suppose that the solution takes the form of a travelling combustion wave and track the location of the flame brush by a phase-field-like technique. We build a finite volume fractional step scheme for both models, which respects the natural physical bounds of the unknowns. Then we assess the consistency of both formulations.

Equations

$${\partial }_t\rho +\mathrm{div}(\rho \mathbf{u})=0,$$

$${\nu }_{F}F+{\nu }_{O}O+N\to {\nu }_{P}P+N,$$

$${\partial }_t(\rho y_i)+\mathrm{div}(\rho y_i\mathbf{u})={\dot{\omega }}_i,\text{mbox}for1\le i\le N_s,$$

$${\dot{\omega }}_F=-{\nu }_F{W}_F\dot{\omega },{\dot{\omega }}_O=-{\nu }_O{W}_0\dot{\omega },{\dot{\omega }}_P={\nu }_P{W}_P\dot{\omega }\text{mbox}and{\dot{\omega }}_N=0,$$

$$\sum_{i\in\mathcal{I}}c_{{\rm p},i}\ \Bigl{[}\partial_{t}(\rho y_{i}\theta)+{\rm div}(\rho y_{i}\theta{\boldsymbol{u}})\Bigr{]}-{\rm div}(\lambda{\boldsymbol{\nabla}}\theta)=\dot{\omega}_{\theta},\quad\dot{\omega}_{\theta}=-\sum_{i\in\mathcal{I}}\Delta h_{f,i}^{0}\dot{\omega}_{i}.$$

$$\rho=\varrho\bigl{(}\theta,(y_{i})_{1\leq 1\leq N_{s}}\bigr{)}=\frac{1}{\displaystyle\frac{R\theta}{P_{th}}\sum_{i=O,P,N}\frac{y_{i}}{W_{i}}+\frac{y_{F}}{\rho_{F}}},$$

$$z=\frac{s{y}_F+1-y_O}{1+s},\text{mbox}withs=\frac{{\nu }_O{W}_O}{{\nu }_F{W}_F}.$$

$$(y_N,z,y_F,y_P{)}^{n+1}$$

$$\forall K\in \mathcal{M},$$

$$\displaystyle\hskip 12.91663pt\dfrac{1}{\delta t}\bigl{[}\rho^{n}_{K}(y_{N})_{K}^{n+1}-\rho^{n-1}_{K}(y_{N})_{K}^{n}\bigr{]}+{\rm div}\Bigl{[}\rho^{n}y_{N}^{n+1}{\boldsymbol{u}}^{n}\bigr{]}_{K}=0,$$

$$\displaystyle\hskip 12.91663pt\dfrac{1}{\delta t}\bigl{[}\rho^{n}_{K}z_{K}^{n+1}-\rho^{n-1}_{K}z_{K}^{n}\bigr{]}+{\rm div}\bigl{[}\rho^{n}z^{n+1}{\boldsymbol{u}}^{n}\bigr{]}_{K}=0,$$

$$\displaystyle\hskip 12.91663pt\dfrac{1}{\delta t}\bigl{[}\rho^{n}_{K}(y_{F})_{K}^{n+1}-\rho^{n-1}_{K}(y_{F})_{K}^{n}\bigr{]}+{\rm div}\bigl{[}\rho^{n}y_{F}^{n+1}{\boldsymbol{u}}^{n}\bigr{]}_{K}=(\dot{\omega}_{F})_{K}^{n+1},$$

$$(y_F{)}_K^{n+1}+(y_O{)}_K^{n+1}+(y_N{)}_K^{n+1}+(y_P{)}_K^{n+1}=1.$$

$${\theta }^{n+1}$$

$$\forall K\in \mathcal{M},$$

$$\displaystyle\hskip 12.91663pt\begin{matrix}\displaystyle\sum_{i\in\mathcal{I}}c_{{\rm p},i}\Bigl{[}\dfrac{1}{\delta t}\bigl{[}\rho^{n}_{K}(y_{i})_{K}^{n+1}\theta_{K}^{n+1}-\rho^{n-1}_{K}(y_{i})_{K}^{n}\theta_{K}^{n}\bigr{]}\Bigr{]}+{\rm div}\bigl{[}\rho^{n}y_{i}^{n+1}\theta^{n+1}{\boldsymbol{u}}^{n}\bigr{]}_{K}\hskip 43.05542pt\\ \hfill-{\rm div}(\lambda{\boldsymbol{\nabla}}\theta^{n+1})_{K}=(\dot{\omega}_{\theta})^{n+1}_{K}.\end{matrix}$$

$$\rho_{K}^{n+1}=\varrho\bigl{(}\theta_{K}^{n+1},((y_{i})_{K}^{n+1})_{1\leq 1\leq N_{s}}\bigr{)}$$

$${\mathbf{u}}^{n+1}$$

$$\displaystyle\quad\forall K\in{\mathcal{M}},\qquad\dfrac{1}{\delta t}\bigl{[}\rho_{K}^{n+1}-\rho_{K}^{n}\bigr{]}+{\rm div}\bigl{[}\rho^{n+1}{\boldsymbol{u}}^{n+1}\bigr{]}_{K}=0.$$

$${\partial }_t(\rho G)+\mathrm{div}(\rho G\mathbf{u})+{\rho }_u{u}_f|\nabla G|=0,$$

$$\dot{\omega }=\frac{u_f}{\delta }\eta (y_F,y_O)(G-0.5{)}^{-},\eta (y_F,y_O)=\min (\frac{y_F}{{\nu }_F{W}_F},\frac{y_O}{{\nu }_O{W}_O}),$$

$$G^{n+1}$$

$$\displaystyle\quad\forall K\in{\mathcal{M}},\qquad\dfrac{1}{\delta t}\bigl{[}\rho^{n}_{K}G_{K}^{n+1}-\rho^{n-1}_{K}G_{K}^{n}\bigr{]}+{\rm div}\Bigl{[}\rho^{n}G^{n+1}{\boldsymbol{u}}^{n}\bigr{]}_{K}+\rho_{u}u_{f}|{\boldsymbol{\nabla}}G|^{n+1}_{K}=0.$$

$$|\nabla G|=\frac{\nabla G}{|\nabla G|}\cdot \nabla G,\text{mbox}so|\nabla G{|}_K^{n+1}=({\mathbf{N}}_f^n\cdot \nabla G^{n+1}{)}_K,$$

$$(y_F{)}_0=(y_O{)}_0=0.4,(y_N{)}_0=0.2,(y_P{)}_0=0,{\theta }_0=300^{\circ }K.$$

$${\dot{\omega }}_K=10^4{y}_F{y}_O{e}^{-900/\theta }.$$

$$({\rho }_u-{\rho }_b)u_p={\rho }_u{u}_u-{\rho }_b{u}_b,$$

Full text

1. Introduction

We address in this paper two alternative models dedicated to the simulation of laminar flames in dust suspensions in a gaseous atmosphere, for which a one-dimensional representation, supposing a low Mach number flow, is sufficient. The combustible particules are supposed to be in mechanical and thermal equilibrium with the continuous phase (or, in other words, no drift nor temperature deviation between the gas and solid phases is taken into account). We consider two descriptions of the combustion phenomenon:

• –

the first one is obtained by collecting mass balance for the chemical species, the energy balance and the momentum balance for the mixture; the reaction term $\dot{\omega}$ is expressed by a closure law depending of the temperature, derived on the basis of physical arguments. This model will be refered to in the following as the primitive formulation.

• –

The second one relies on the assumption that the solution consists in a travelling reaction thin interface (the so-called flame front) separating a zone where the combustion is complete (the ”burnt zone”) from a zone where no combustion has yet occured (the ”fresh zone”). This representation offers the possibility to reduce the problem to an explicit tracking of the front location, through the solution of a transport equation for a color function $G$ ($G\in[0,1]$, $G<0.5$ in the burnt zone, $G\geq 0.5$ in the fresh atmosphere); the reaction term is governed by the value of $G$: $\dot{\omega}=0$ if $G\geq 0.5$ and $\dot{\omega}$ is proportional to $1/\tau$ otherwise, where $\tau$ is a time-scale closely correlated to the flame front thickness. In the rest of this paper, we will call this model the flame velocity formulation.

The first option is standard for the computation of laminar flames. Variants of the flame velocity formulation are often chosen to compute turbulent deflagrations in industrial applications [5, 4], as in nuclear safety studies performed at the French Institut de Radioprotection et de Sûreté Nucléaire (IRSN). Indeed, this latter model seems easier to solve, in the sense that stable segregated algorithms may be designed for its purpose; in addition, the flame brush incorporates structures which are very small compared to the system scales, and the flame velocity approach allows an upscaling of this complex physical phenomenon through a single parameter (the turbulent flame velocity) which may be inferred from experimental data.

A finite volume fractional step numerical scheme was developed in [2] for the solution of the system of primitive equations; we shortly describe it here and review its main properties. The aim of this paper is then to assess the accuracy of the switch from the first model to the second one: first, we check that the solution to the primitive formulation, in conditions representative of the target physical reality, is indeed a flame front propagating through the medium; then, we compare such a solution with the one obtained with the corresponding flame velocity model.

2. The primitive formulation

2.1. The governing equations

The flow is supposed to be governed by the balance equations modelling a variable density flow in the asymptotic limit of low Mach number flows, namely the mass balance of the chemical species and of the mixture, the enthalpy balance, and the momentum balance equations. For a one-dimensional flow in such a quasi-incompressible situation, the role played by the mass and momentum balance equations is quite different than in the multi-dimensional case: the velocity may be seen as the solution of the mass balance equation, and the momentum balance yields the dynamic pressure. Since this latter unknown does not appear in the other equations, its computation is of poor interest, and the momentum balance equation may be disregarded.

Except for this aspect, equations in this section are written in the usual multi-dimensional form. The computational domain is denoted by $\Omega$, and its boundary $\partial\Omega$ is supposed to be split into an inflow part $\partial\Omega_{I}$ (where the flow enters the domain, i.e. ${\boldsymbol{u}}\cdot{\boldsymbol{n}}_{\partial\Omega}<0$, with ${\boldsymbol{u}}$ the flow velocity and ${\boldsymbol{n}}_{\partial\Omega}$ the normal vector to $\partial\Omega$ outward $\Omega$) and an outflow one $\partial\Omega_{O}$ (where the flow leaves the domain, i.e. ${\boldsymbol{u}}\cdot{\boldsymbol{n}}_{\partial\Omega}\geq 0$) of positive $(d-1)$-measure. The problem is posed over the time interval $(0,T)$.

Mass balance equations – The mass balance reads:

[TABLE]

where $\rho$ stands for the fluid density. This equation must be complemented by an initial condition and a boundary condition on $\partial\Omega_{I}$ for the density; both functions are obtained by the data of the temperature and flow composition, thanks to the equation of state (see below).

Only four chemical species are supposed to be present in the flow, namely the dust, or fuel (denoted by $F$), the oxydant ($O$), the product ($P$) of the reaction, and a neutral gas ($N$). A one-step irreversible total chemical reaction is considered:

[TABLE]

where $\nu_{F}$, $\nu_{O}$ and $\nu_{P}$ are the molar stoichiometric coefficients of the reaction. Chemical species other than the fuel are supposed to be gases. The system of the mass balance equations for the chemical species reads:

[TABLE]

where $y_{i}$ and $\dot{\omega}_{i}$ stand respectively for the mass fraction and the reaction rate of the species $i$. The number of species is denoted by $N_{s}$, with, by assumption, $N_{s}=4$, and we indifferently use the notation $(y_{i})_{1\leq i\leq N_{s}}$ or $y_{F}$, $y_{0}$, $y_{P}$ and $y_{N}$ for the fuel, oxydant, product and neutral gas mass fractions, respectively. To simplify the exposition, the mass diffusion fluxes have been supposed in the set (2) of equations to vanish. This system must be complemented by initial and a Dirichlet boundary conditions for $(y_{i})_{1\leq i\leq N_{s}}$ on the inflow part of the domain boundary $\partial\Omega_{I}$. The prescribed values of the mass fractions at the initial time and on the inflow boundary lie in the interval $[0,1]$. The reaction rate of each chemical species may be written as:

[TABLE]

where $W_{F}$, $W_{0}$ and $W_{P}$ stand for the molar masses of the fuel, oxydant and product respectively, and $\dot{\omega}$ is a non-negative reaction rate, which is supposed to vanish when either $y_{F}=0$ or $y_{O}=0$. Since $\nu_{F}W_{F}+\nu_{O}W_{0}=\nu_{P}W_{P}$, we have $\dot{\omega}_{F}+\dot{\omega}_{O}+\dot{\omega}_{P}=0$.

Energy balance – In the low Mach number approximation, the total enthalpy balance reads:

[TABLE]

where $\theta$ stands for the temperature, $c_{{\rm p},i}$ for the specific heat of the species $i$ (supposed to be constant), $\Delta h_{f,i}^{0}$ for the formation enthalpy at $0^{\circ}$K and $\lambda$ the thermal conductivity. This equation is complemented by a total flux boundary condition at the inlet boundary, and we suppose that the diffusion flux vanishes at the outlet boundary.

Equation of state – We suppose that the gas phase is a mixture of perfect gases and that the density $\rho_{F}$ of the solid phase (i.e. of the fuel) is constant, so:

[TABLE]

where $R=8.31451\,\mbox{JK}^{-1}\mbox{mol}^{-1}$ stands for the perfect gases constant. Since the computational domain is supposed not to be closed, the so-called thermodynamic pressure $P_{th}$ is constant in time and space, and given by the initial state.

2.2. The numerical scheme

For the solution of the equations of the model, we define the variable $z$ as follows:

[TABLE]

Note that, combining the mass balance equation for the fuel and the oxydant, the variable $z$ satisfies an homogeneous equation; for this reason, we replace the oxydant mass balance equation by the balance equation for $z$ (since, given the values of $z$ and $y_{F}$, we may deduce $y_{O}$).

Let us consider a partition $0=t_{0}<t_{1}<\ldots<t_{N}=T$ of the time interval $(0,T)$, which we suppose uniform. Let $\delta t=t_{n+1}-t_{n}$ for $n=0,1,\ldots,N-1$ be the constant time step. We suppose that the interval $\Omega$ is split into a family of control volumes (sub-intervals of $\Omega$) which realizes a partition of $\Omega$; we denote these control volumes by $(K)_{K\in{\mathcal{M}}}$. The scalar unknowns, i.e. the density, mass fractions and temperature, are associated to the control volumes, and the corresponding unknowns read $\rho_{K}^{n}$, $(y_{i})_{K}^{n}$, $z_{K}^{n}$ and $\theta_{K}^{n}$ for $K\in{\mathcal{M}}$, $0\leq n\leq N$ and $i\in\mathcal{I}$. The velocity is discretized at the faces of the mesh, which we denote by $(\sigma)_{\sigma\in\mathcal{E}}$, so the corresponding unknowns are $u^{n}_{\sigma}$ for $\sigma\in\mathcal{E}$ and $0\leq n\leq N$. We implement a fractional-step algorithm, which consists in four steps, in order to calculate recursively the unknowns $(y_{i})^{n+1}_{i\in\mathcal{I}}$, $z^{n+1}$, $\theta^{n+1}$, $\rho^{n+1}$ and $u^{n+1}$ for $0\leq n<N$:

[TABLE]

The discrete operators appearing in these relations are approximated by finite-volume techniques. Thanks to a careful definition of the convection fluxes, derived to fulfill the conditions introduced in [3] to obtain a maximum-principle-preserving convection operators, the scheme is proven in [2] to preserve the physical bounds of the unknowns: mass fractions in the interval $[0,1]$, positivity of the temperature and so of the density.

3. A model based on an explicit tracking of the flame front

3.1. The governing equations

The physical model addressed in this section is based on an explicit computation of the flame brush location, by a phase-field-like technique. The ”color function” is called $G$ in this context; its transport equation is referred to as the $G$-equation, and reads:

[TABLE]

Initial conditions are $G=0$ at the location where the flame starts and $G=1$ elsewhere. The quantity $\rho_{u}$ is a constant density, which, from a physical point of view, stands for a characteristic value for the unburnt gases density, and $u_{f}$ is the flame brush velocity. The reactive term $\dot{\omega}$ is given by:

[TABLE]

where $\delta$ is a quantity homogeneous to a length scale, which governs the thickness of the reaction zone.

The flame velocity model consists of Equation (6), of the mixture mass balance equation (1), the chemical species mass balance equations (2) (with the modified expression (7) for the chemical reaction term $\dot{\omega}$) and of the energy balance (3). Note that, under some assumptions which are usually not valid in industrial applications, this model may be simplified: for instance, in perfectly premixed situations (i.e. constant in space initial data for the chemical mass fractions and the temperature) and supposing an infinitely fast chemical reaction (i.e., in the present formalism, making $\delta$ tend to zero), the variable $G$ may be indentified to a progress variable and all the other unknowns (i.e. the mass fractions and the temperature) may be deduced from $G$ through an algebraic relation.

3.2. Numerical scheme

The $G$ function is discretized on the primal mesh, so the discrete unknowns are $G^{n}_{K}$, for $K\in{\mathcal{M}}$ and $0\leq n\leq N$. The numerical algorithm differs from the scheme for the primitive formulation, i.e. System (5), by the insertion, as a first step, of a discrete counterpart to Equation (6):

[TABLE]

For the discretization of the last term in this relation, we write:

[TABLE]

where ${\boldsymbol{N}}_{f}$ is an approximation of the advection field ${\boldsymbol{\nabla}}G/|{\boldsymbol{\nabla}}G|$ and we use an upwind finite volume formulation of the transport operator, i.e. the formulation obtained by writing ${\boldsymbol{N}}_{f}\cdot{\boldsymbol{\nabla}}G={\rm div}(G\,{\boldsymbol{N}}_{f})-G\,{\rm div}{\boldsymbol{N}}_{f}$ and using an upwind finite volume (first or second order) discretization of the convection operator. For the present solver, the convection operators are discretized by an explicit MUSCL-like technique [6].

4. Results

Computations presented in this section are performed with MATLAB for the primitive formulation and by the open-source CALIF3S software developped at IRSN [1] for the flame-velocity model.

Data is chosen in order to allow to check the scheme properties (i.e. to avoid unrealistic simplifications, as, for instance, a same specific heat diffusion coefficient for all the chemical species), and to be in the range of practical applications. The mixture is initially at rest, homogeneous and with a uniform temperature:

[TABLE]

In the primitive formulation, the reaction rate follows an Arrhenius law:

[TABLE]

The molar masses of the chemical species are considered to be equal to $20$ g mol*-1* for all the species, so the combustion reaction reads $F+O+N\longrightarrow 2P+N$, and the initial atmosphere composition is stoichiometric. The temperature diffusion coefficient is $\lambda=0.005$, the specific heat coefficients (J Kg*-1* K*-1*) are $c_{{\rm p},N}=3.\,10^{3}$, $c_{{\rm p},F}=1.\,10^{3}$, $c_{{\rm p},O}=2.\,10^{3}$ and $c_{{\rm p},P}=4.\,10^{3}$ and the formation enthalpies (J Kg*-1*) are $\Delta h_{f,N}^{0}=3.\,10^{6}$, $\Delta h_{f,F}^{0}=1.\,10^{6}$, $\Delta h_{f,O}^{0}=-2.\,10^{6}$ and $\Delta h_{f,P}^{0}=-4.\,10^{6}$ (so the reaction is exothermic). The fuel density is equal to $100$  Kg m*-3*. Ignition is obtained in the primitive formulation by making $\dot{\omega}$ depend in a very thin zone on a fictitious elevated temperature, to trigger a reaction at the initial time. In the flame velocity model, $G$ is imposed to zero in the same zone (while $G=1$ elsewhere). Since the inflamation zone is very thin, the consequent initial burst is not too violent.

First of all, we observe that, for the primitive formulation, the solution tends after a transition period to a travelling combustion wave separating a fresh (or unburnt) zone from a burnt zone, where $y_{F}=y_{O}=0$, $y_{P}=0.8$ and the temperature is equal to the adiabatic combustion temperature. By construction of the scheme, the neutral gas mass fraction $y_{N}$ and the reduced variable $z$ are kept constant in time and space and equal to their initial value. Since the profile in the interface does not vary in space and time up to a translation velocity $u_{p}$ (the velocity of the flame brush), we may write the jump conditions for the mixture mass balance equation, to obtain:

[TABLE]

where $\rho_{b}$ and $u_{b}$ (resp. $\rho_{u}$ and $u_{u}$) stand for the constant density and velocity in the burnt (resp. unburnt) zone. Thanks to symmetry conditions (due to the fact that the combustion takes place in an atmosphere initially at rest), $u_{b}=0$ and we deduce from the previous relation that the flame velocity is given by: $u_{f}=u_{p}-u_{u}=u_{u}\ \rho_{b}/(\rho_{u}-\rho_{b}).$ The obtained value is injected in the flame velocity model, and we choose the length $\delta_{f}$ to fit as closely as possible the travelling profiles of the unknowns. Results for the fuel mass fraction and the temperature are given on Figure 1. We observe that, as expected, the thickness of the combustion zone is scaled by $\delta_{f}$ and that a reasonable agreement is obtained with $\delta_{f}=0.1$ mm.