Hybrid Discretization Methods for Transient Numerical Simulation of Combustion in Porous Media
aa r X i v : . [ m a t h . NA ] A ug Hybrid Discretization Methodsfor Transient Numerical Simulationof Combustion in Porous Media
Peter Knabner ∗ and Gerhard Summ Institute of Applied Mathematics, Martensstraße 3, D-91058 Erlangen, Germany
Abstract.
We present an algorithm for the numerical solution of the equationsgoverning combustion in porous inert media. The discretization of the flow problem isperformed by the mixed finite element method, the transport problems are discretizedby a cell-centered finite volume method. The resulting nonlinear equations arelineararized with Newton’s method, the linearized systems are solved with a multigridalgorithm. Both subsystems are recoupled again in a Picard iteration. Numericalsimulations based on a simplified model show how regions with different porositystabilize the reaction zone inside the porous burner
Submitted to:
Combust. Theory Modelling
1. Introduction
In recent years, the request for low-emission combustion systems has led to thedevelopment of a new burner concept: combustion in porous media. In these burners apremixed gas-air mixture is constrained to flow through and combust within the poresof a porous medium, typically a ceramic foam. Due to the presence of the solid matrix,the temperature in the combustion zone can be controlled such that the emission ofNO x and other pollutants can be lowered. Further advantages of this technique arementioned in [4]. A list of technical applications can be found in [13].Several authors performed numerical simulations to support the development ofporous burners. Most of these simulations are based on mathematical models thatare spatially one-dimensional [8] or stationary in time [4],[12], or one-dimensional andstationary [9]. The employed numerical algorithms are often simple extensions ofexisting codes designed for other applications and use finite difference or finite volumemethods for the discretization of the governing equations. Some of these simulationsconsider the pressure distribution and the velocity field of the flow as given. In [4] and[12] the pressure distribution and the velocity field are computed from the continuityequation and a momentum equation. Then pressure and velocity are coupled using the ∗ [email protected] imulation of Combustion in Porous Media
2. Mathematical model
The simulation of the flow and transport processes inside the porous medium basedon models that resolve the complex pore structure in detail would be computationallytoo expensive. Hence these simulations are based on homogenized or volume-averagedequations, where the fluid and the solid are considered as a pseudohomogeneous medium.We restrict our considerations here to the simplest possible model, which can be usedto describe the flow of the gas in the porous burner and the transport of heat andchemical species. This model is sufficiently realistic to present the essential features ofour algorithm.
The flow in the porous medium is governed by the Darcy–Forchheimer equation µk u + β Fo ρ | u | u + ∇ p = 0 , where | · | denotes the Euclidean norm, the continuity equation φ ∂ρ∂t + div( ρ u ) = 0and the ideal gas law ρ = p WR T .
The unknowns here are the pressure p , the density ρ and the volumetric flow rate u of the gas. Porosity φ , permeability k and Forchheimer coefficient β Fo of the porousmedium, viscosity µ , molecular weight W and temperature T of the gas mixture, andthe universal gas constant R are given. Assuming ρ > S = | p | p and m = | ρ | u these equations can be transformed into( α + β | m | ) m + ∇ S = 0 , (1) φ ∂ t ρ ( S ) + div( m ) = 0 , (2) imulation of Combustion in Porous Media γ := WR T , α := 2 µγ k , β := 2 β Fo γ , and the equation of state ρ = ρ ( S ) is defined by ρ ( S ) := γ S p | S | . The solid matrix of the porous burner typically consists of ceramic foams havinghigh heat transfer rates. Furthermore the specific surface of these foams is very high.Therefore we can assume thermal equilibrium between the gas and the solid. Thetemperature of this pseudo-homogeneous medium is governed by the (effective) heatequation φc p ∂ t ( ρT ) + (1 − φ ) c s ∂ t ( ρ s T ) + div( c p m T − λ eff ∇ T ) = φQ ˙ r + (1 − φ ) F Q , (3)where the effective heat conductivity λ eff := φλ g + (1 − φ ) λ s is a mean value of the heatconductivities λ g of the gas and λ s of the solid. More realistic representations of λ eff ,including effects of radiation and dispersion (see [4]), could also be used. Furthermore c p denotes the heat capacity for constant pressure of the gas, c s the specific heat capacityand ρ s the density of the solid. Q ˙ r is the heat produced by the combustion reaction and F Q the power density of an external heat source needed for ignition.Consider the simplest case only, we restrict the reaction model to an irreversibleone-step reaction mechanism (e.g. methane oxidation CH + 2O −→ CO + 2H O).Then it suffices to consider the conservation equation of the reactant
Rφ∂ t ( ρy ) + div( m y − φD ∇ y ) = − φ ˙ r , (4)where y denotes the mass fraction of the reactant and the sink term − ˙ r models theconsumption of the reactant by the chemical reaction. Following the Arrhenius model˙ r is given by ˙ r = Bρy exp (cid:18) − ER T (cid:19) with the frequency factor B and the activation energy E . Since (1) and (2) are equivalent to a single parabolic equation, we have to prescribeinitial conditions for S and boundary conditions for S or m only. Hence we have toprovide the following initial conditions for the problem consisting of equations (1)–(4): S ( · ,
0) = S in Ω , T ( · ,
0) = T in Ω , y ( · ,
0) = y in Ω , where Ω ⊂ R is the domain corresponding to the porous burner. Then the initialcondition for the mass flux m ( · ,
0) can be calculated from (1).The boundary conditions are listed in Table 1. Note that the mixed boundaryconditions for the heat equation model thermal equilibrium at the inflow boundary andNewton’s law of cooling at the burner wall. The coefficients h b and T b denote the heattransfer coefficient and the ambient temperature, resp. imulation of Combustion in Porous Media Table 1.
Boundary conditions
Flow problem :Inflow boundary Γ I : m · n = m b Outflow boundary Γ O : S = S b Burner wall Γ W : m · n = 0Line of symmetry Γ S : m · n = 0Heat equation :Inflow boundary Γ I : λ eff ∇ T · n = c p ( T − T b ) m b Outflow boundary Γ O : ∇ T · n = 0Burner wall Γ W : λ eff ∇ T · n = h b ( T b − T )Line of symmetry Γ S : ∇ T · n = 0Reactant conservation equation :Inflow boundary Γ I : y = y b Outflow boundary Γ O : ∇ y · n = 0Burner wall Γ W : ∇ y · n = 0Line of symmetry Γ S : ∇ y · n = 0
3. Numerical solution algorithm
Using the Rothe method, we first discretize equations (1)–(4) in time with the implicitEuler method. Hence, in each time step we have to solve the following coupled system:( α + β | m | ) m + ∇ S = 0 , (5) φ ∆ t ( ρ ( S ) − ρ − ) + div( m ) = 0 , (6) φc p ∆ t ( ρ ( S ) T − ρ − T − ) + (1 − φ ) c s ∆ t ( ρ s T − ρ − s T − ) (7)+ div ( c p m T − λ eff ∇ T ) = φQ ˙ r + (1 − φ ) F Q ,φ ∆ t ( ρ ( S ) y − ρ − y − ) + div ( m y − φD ∇ y ) = − φ ˙ r . (8)Here ρ − , T − and y − denote the values of the unknowns from the last time step.The coupling between the flow problem (governed by equations (5) and (6)) andthe transport problems (governed by equations (7) and (8)) is rather weak. Thereforewe can decouple the flow problem and the transport problems in each time step andcompute their solutions ( m , S ) and ( T, y ) alternatingly in a Picard-iteration, as sketchedin Figure 1. This strategy allows us to use different methods for the discretization ofthe flow and the transport problem. imulation of Combustion in Porous Media ✗✖ ✔✕ m , S ❄ transportproblemsolver ✗✖ ✔✕ T , y ✻ Figure 1.
Decoupling of flow and transport problem
The spatial discretization of both problems is based on a partition T h of the domainΩ into triangular elements K . We denote by E h the set of edges of T h . Respectingthe partition of ∂ Ω into Dirichlet boundary Γ D (where S = S b is given) and Neumannboundary Γ N (where m · n = m b is given), E h can be subdivided into three disjointsubsets E h = E Ih ∪ E Dh ∪ E Nh . Here E Ih := { e ∈ E h | e ∂ Ω } is the set of inner edges, E Dh := { e ∈ E h | e ⊂ Γ D } is the set of edges lying on the Dirichlet boundary and E Nh := { e ∈ E h | e ⊂ Γ N } is the set of edges lying on the Neumann boundary. Forevery K ∈ T h we denote by P k ( K ), k ≥
0, the set of polynomials of degree ≤ k on K .Analogously, we define P k ( e ) for k ≥ e ∈ E h .To preserve the form given by (5) and (6), and to get a good approximation forthe mass flux density m , which appears also in equations (7) and (8), we use the mixedfinite element method on Raviart–Thomas elements of lowest order (see e.g. [5]) for thespatial discretization of (5) and (6). Thus m is approximated by m h ∈ V h , where V h := RT (Ω; T h ) := { v h ∈ H (div; Ω) | v h | K ∈ RT ( K ) for all K ∈ T h } and RT ( K ) is defined by RT ( K ) := ( P ( K )) + (cid:18) xy (cid:19) P ( K ) , and S is approximated by S h ∈ Q h , where Q h := (cid:8) q h ∈ L (Ω) | q h | K ∈ P ( K ) for all K ∈ T h (cid:9) . Note that any v h | K ∈ RT ( K ) is uniquely defined by the flux across the edges of Kv e := Z e v h | K · n e ds , e ⊂ ∂K , where n e is an arbitrarily oriented unit normal vector to e ∈ E h . The property v h ∈ H (div; Ω) requires the continuity of these flux values. Using the basis { w e } e ∈E h of V h satisfying Z f w e · n f ds = ( , if e = f , , if e = f , imulation of Combustion in Porous Media v e , e ∈ E h .To take into account the flux boundary conditions on Γ N , we consider the subspace V m b , Γ N h of V h , which is defined by V m b , Γ N h := (cid:26) v h ∈ V h (cid:12)(cid:12)(cid:12) Z e v h · n ds = Z e m b ds for all e ∈ E Nh (cid:27) . If m b ≡ N , we obtain the space V , Γ N h .Then the discrete mixed formulation reads as follows:Find ( m h , S h ) ∈ V m b , Γ N h × Q h such that for every ( v h , q h ) ∈ V , Γ N h × Q h Z Ω ( α + β | m h | ) ( m h · v h ) dx − Z Ω div( v h ) S h dx + Z Γ D S b ( v h · n ) ds = 0 , (9) Z Ω φ ∆ t ρ ( S h ) q h dx + Z Ω div ( m h ) q h dx = Z Ω φ ∆ t ρ − q h dx . (10)Unfortunately, the systems of algebraic equations resulting from the mixed finiteelement discretization are difficult to solve numerically. Therefore we apply the followingimplementational technique, called hybridization to our discrete mixed formulation: Weeliminate the continuity constraints in the definition of V h and enforce the requiredcontinuity instead through additional equations involving Lagrange multipliers definedon the edges e ∈ E h . Thus we replace V h by W h := RT − (Ω; T h ) := n v h ∈ (cid:0) L (Ω) (cid:1) | v h | K ∈ RT ( K ) for all K ∈ T h o , and V m b , Γ N h by the corresponding subspace W m b , Γ N h of W h . In addition, we define thespace of Lagrange multipliers byΛ S b , Γ D h := (cid:26) λ h ∈ L ( E h ) (cid:12)(cid:12)(cid:12) λ h | e ∈ P ( e ) ∀ e ∈ E h , Z e ( λ h − S b ) ds = 0 ∀ e ∈ E Dh (cid:27) , where E h = ∪ e ∈E h e . Then the hybridized mixed formulation reads as:Find ( m h , S h , µ h ) ∈ W m b , Γ N h × Q h × Λ S b , Γ D h , such that for every ( v h , q h , λ h ) ∈ W , Γ N h × Q h × Λ ,∂ Ω h Z Ω ( α + β | m h | ) ( m h · v h ) dx − X K ∈T h Z K div ( v h ) S h dx + X K ∈T h Z ∂K µ h ( v h · n K ) ds = 0 , (11) Z Ω φ ∆ t ρ ( S h ) q h dx + X K ∈T h Z K div ( m h ) q h dx = Z Ω φ ∆ t ρ − q h dx , (12) X K ∈T h Z ∂K λ h ( m h · n K ) ds = 0 . (13)The solutions m h and S h of (11)–(13) coincide with the solutions m h and S h of(9)–(10) (cf. [10]). Therefore we are allowed to use the same notation for them. Wenote that the additionally computed Lagrange multipliers can be used to construct amore accurate approximation for S (see e.g. [5] and the references cited there).Note that every v h ∈ W h is uniquely defined by the degrees of freedom v K,e , whichare given by v K,e = Z e v h | K · n K ds , K ∈ T h , e ⊂ ∂K , imulation of Combustion in Porous Media n K denotes the unit outer normal of K . The corresponding basis vectors aredenoted by w K,e . Thus the unknown functions m h , S h and µ h can be represented by m h = X K ∈T h X e ∈ ∂K m K,e w K,e , S h = X K ∈T h S K χ K and µ h = X e ∈E h µ e χ e , where χ K and χ e denote the characteristic functions of an element K and an edge e ,resp..Employing the basis functions w K, ¯ e for ¯ e / ∈ E Nh and K ∈ T h with ¯ e ⊂ ∂K , χ K for K ∈ T h , and χ e for e ∈ E Ih as test functions, we obtain from (11)–(13) the followingsystem of nonlinear equations. F K, ¯ e := Z K ( α + β | m h | ) ( m h · w K, ¯ e ) d x − S K + µ ¯ e = 0 , ¯ e / ∈ E Nh , K ∈ T h , ¯ e ⊂ ∂K , (14) F K := R K φ γ d x ∆ t S K p | S K | + X e ⊂ ∂K m K,e − R K φ ρ − d x ∆ t = 0 , K ∈ T h , (15) F ¯ e := m K, ¯ e + m K ′ , ¯ e = 0 , e ∈ E Ih . (16)For every K ∈ T h the basis functions w K,e ( e ⊂ ∂K ) of W h vanish in Ω \ K . Thereforethe degrees of freedom m K,e for e ⊂ ∂K affect only m h | K = P e ∈ ∂K m K,e w K,e =: m K .Hence, for every K ∈ T h , the degrees of freedom m K,e for e ⊂ ∂K and S K appear onlyin equations (14) and (15) belonging to the element K . This fact can be exploited tocompensate the main drawback of hybridization, the introduction of additional degreesof freedom: For every K ∈ T h , we can eliminate the approximate mass flux values m K,e and the approximate element values S K from the global system by solving locally thesystems consisting of (14) and (15).Owing to the nonlinearity of (14) and (15) it is not possible to find a closed formsolution for m K,e ( e ⊂ ∂K ) and S K . Nevertheless, using monotonicity arguments itcan be shown that the local subsystems are uniquely solvable and their Jacobians areinvertible (see [10, Section 4.3]). Thus it is possible to compute m K,e ( e ⊂ ∂K ) and S K by means of Newton’s method during the assembling procedure. The elimination of m K,e ( e ⊂ ∂K ) and S K introduces nonlinearity into the originally linear equations (16).Using Newton’s method to solve these equations we have to compute the Jacobian ofthe remaining global system. By the implicit function theorem this requires again theinvertibility of the Jacobians of the local subsystems.The linearized global systems are solved by a multigrid method. Chen [6] showedthat – at least for linear elliptic problems – the equations after elimination of fluxand element variables correspond to the equations resulting from a nonconformingdiscretization using the Crouzeix–Raviart ansatz space. Therefore we can employ theintergrid transfer operators developed for the Crouzeix–Raviart ansatz space [3] toconstruct a multigrid algorithm. Using the equation of state ρ = ρ ( S ), we need the values of the temperature T andauxiliary variable S to compute the density ρ . Therefore we choose the same ansatz imulation of Combustion in Porous Media T h as for S h , i.e., T h := P K ∈T h T K χ K ∈ Q h . Furthermore ρ , T and y appearin the reaction rate ˙ r . Hence we approximate the mass fraction y of the reactant by y h := P K ∈T h y K χ K ∈ Q h , too. Since the fluxes corresponding to T or y do not enterfurther equations, we do not need to compute them explicitly using the mixed finiteelement method. Instead, we employ a cell-centered finite volume scheme for the spatialdiscretization of (7) and (8). Note that this finite volume scheme can be derived fromthe mixed finite element method (see [1]).Integrating (7) and (8) over a cell K ∈ T h , applying the divergence theorem andreplacing the continuous unknowns T and y by their discrete approximations T h and y h yields the following equations for all K ∈ T h : R K φc p d x ∆ t ( ρ K T K − ρ − K T − K ) + R K (1 − φ ) c s ρ s d x ∆ t T K − R K (1 − φ ) c s ρ − s d x ∆ t T − K (17)+ X e ⊂ ∂K Z e ( c p m h T h − λ eff ∇ T h ) · n K dσ = Z K φQ ˙ r d x + Z K (1 − φ ) F Q d x , R K φ d x ∆ t ( ρ K y K − ρ − K y − K ) + X e ⊂ ∂K Z e ( m h y h − φD ∇ y h ) · n K dσ = − Z K φ ˙ r d x . (18)To obtain algebraic equations, we have to provide quadrature formulas for theintegrals appearing in (17) and (18). The terms R K φ d x and R K (1 − φ ) F Q d x areevaluated by means of the midpoint rule, i.e., Z K φ d x ≈ | K | φ ( b K ) and Z K (1 − φ ) F Q d x ≈ | K | (1 − φ ( b K )) F Q ( b K ) , where b K is the center of gravity of K . The remaining integrals over K depend onone or several unknowns. Of course, we use the corresponding element value for theseunknowns. Taking into account the Arrhenius law and the temperature dependenceof coefficient functions c p , c s and ρ s , we obtain the following approximations for theseintegrals: Z K φ ˙ r d x = Z K φBρ h y h exp (cid:18) − ER T h (cid:19) d x ≈ | K | φ ( b K ) B ρ K y K exp (cid:18) − ER T K (cid:19) , Z K φc p d x ≈ | K | φ ( b K ) c p ( T K ) , Z K (1 − φ ) c s ρ ( − ) s d x ≈ | K | (1 − φ ( b K )) c s ( T K ) ρ s ( T ( − ) K ) . The approximation of integrals over the edges e ∈ E h is much more difficult. Usingthe generalizing notation z h ∈ { y h , T h } , we have to find approximations of Z e ( c z m h z h − D z ∇ z h ) · n K dσ , where z h ∈ Q h is piecewise constant. In particular, z h is not continuous across interioredges. Consequently, z e := z h | e is not uniquely defined and ∇ z h | e is not defined at all.For an element K ∈ T h and an interior edge e ∈ E Ih let Nb( K, e ) ∈ T h be the element,which shares the edge e with K . imulation of Combustion in Porous Media e ∈ E Ih . In order to discretize theconvective flux in stable way, we use an upwind scheme (see e.g. [7, Section 7]), i.e., Z e c z ( m h · n K ) z h dσ ≈ m + K,e c z | K z K + m − K,e c z | Nb(
K,e ) z Nb(
K,e ) . Here m + K,e := max { m K,e , } and m − K,e := min { m K,e , } .The diffusive flux is approximated by some form of difference quotient. Since theelements K may be of different size and shape and the value of the diffusion coefficient D z may vary from element to element, some points have to be respected to obtain astable discretization. Similar to the computation of the effective heat conductivity oflayered materials, we have to employ some kind of geometric mean for the diffusioncoefficient. The derivation of the cell-centered finite volume method from the mixedfinite element method yields the following discretization (cf. [7, Section 9]): Z e D z ∇ z h · n K dσ ≈ | e | z Nb(
K,e ) − z K d e ( D z ) , where d e ( D z ) = (cid:18) d K,e D z | K + d Nb(
K,e ) ,e D z | Nb(
K,e ) (cid:19) . and d K,e denotes the distance between the center of the circumcircle of K and the edge e . Since the center of the circumcircle of a triangle with an obtuse angle lies outside ofthe triangle, such obtuse-angled triangles have to be avoided in the triangulation T h .Finally, we describe the approximation of the fluxes across boundary edges. In thecase of Dirichlet boundary conditions ( z = z b on e ) we employ the approximation Z e ( c z m h z h − D z ∇ z h ) · n K dσ ≈ c z ( z K ) m K,e z b,e − | e | D z ( b K ) z b,e − z K d K,e
Here z b,e is the mean value of z b on e , approximated by z b,e := 1 | e | Z e z b dσ ≈ ( z b ( n ( e )) + z b ( n ( e ))) / , where n i ( e ), i = 1 ,
2, denote the vertices of the edge e .In the case of Neumann boundary conditions ∇ z · n = q b is given. To evaluate theconvective flux, we choose z e ≈ z K . Hence we obtain Z e ( c z m h z h − D z ∇ z h ) · n K dσ ≈ c z ( z K ) m K,e z K − D z ( b K ) Z e q b dσ ≈ c z ( z K ) m K,e z K − | e | D z ( b K ) ( q b ( n ( e )) + q b ( n ( e ))) / . Finally, we consider the case of mixed boundary conditions, given in the form z + 1 σ ( D z ∇ z · n ) = g ⇐⇒ D z ∇ z · n = σ ( g − z ) . In this case we employ the approximation Z e ( c z m h z h − D z ∇ z h ) · n K dσ ≈ c z ( z K ) m K,e z K − | e | σ e ( g e − z K ) . imulation of Combustion in Porous Media σ e and g e are defined like z b,e above. Note that the cell-centered finite volumemethod can be extended to more general elliptic operators, including discontinuousmatrix diffusion coefficients (see [7, Section 11]).After approximating the integrals in (17) and (18) as above, we arrive at a coupledsystem of nonlinear equations. Again, we propose to use Newton’s method for thelinearization of these equations and to use a multigrid method for the solution of thelinearized equations. For the computations presented below, we used the trivial injectionoperator I tk and its transpose as intergrid transfer operators and obtained satisfactoryconvergence rates. As there may be cases where the V -cycle does not converge for thesetrivial intergrid transfer operators, we mention also the weighted interpolation operator I wk , which has been proposed in [11].
4. Stabilization of the reaction zone
In order to demonstrate the facilities of the discretization methods presented above,we implemented them in the framework of the software toolbox ug [2]. By meansof this algorithm, we studied, how regions with varying porosity can help to controlthe position of the reaction zone inside the porous burner. To this end we considerthe simplified model sketched in Figure 2 of a porous burner prototype developed byTrimis and Durst [13]. Experiments showed that the reaction zone stabilizes at the (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) (cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1)(cid:1) ΓΓ C1 C2 Γ I Γ O Γ Γ S Γφ l φ h Figure 2.
Prototype of a porous burner interface between the preheating region with smaller pores (or lower porosity φ l ) andthe combustion region with larger pores or (or higher porosity φ h ). It was suggestedthat this stabilization effect appears when the pore diameter in the preheating region islower than a quenching diameter observed in tube experiments. But the pore diameteris a microscopic parameter, which does not appear in the macroscopic (averaged)equations (1)–(4). Therefore the numerical simulations presented in [4] implementthis stabilization effect by permitting chemical reactions only in the combustion region.Our simulations show, that this stabilization effect can be explained with the aid ofmacroscopic parameters only. Of course, these simulations can not yield quantitativelyrealistic results, since they are based on the simplified model consisting of equations(1)–(4). Nevertheless we expect to obtain qualitatively correct results concerning basicfeatures of porous burners like the stabilization effect of regions with varying porosity. imulation of Combustion in Porous Media The values of the coefficients used in our simulations are listed in Table 2. For theForchheimer parameter β Fo , we use Ergun’s relation β Fo = c F / √ k , where c F is a constantand k the permeability of the porous medium. Table 2.
Values of coefficient functions
Gas mixture Solid Reaction rate µ . · − Pa s c F . R .
314 J / mol K D . · − kg / m s ρ s / m B . · / s c p / kg K c s
765 J / kg K E / mol λ g .
049 W / K m λ s
36 W / K m Q . · J / kg W .
028 kg / molIn order to study the influence of varying porosity (and permeability) on thestability and position of the reaction zone, we compare the results of numericalsimulations for the following choices of φ and k :a) Different values of φ l and φ h : φ ( x ) = ( . x < . , . x > .
08 and k ( x ) = ( . · − m for x < . , . · − m for x > . . (19)b) Constant low value φ ≡ φ l : φ ( x ) = 0 . x and k ( x ) = 1 . · − m for all x . (20)c) Constant high value φ ≡ φ h : φ ( x ) = 0 . x and k ( x ) = 1 . · − m for all x . (21)We start the simulation with constant initial values S ≡ , y ≡ . T ≡
298 K . Changing the flux boundary condition m ( x , t ) · n = m b ( x , t ) at the inflow boundary Γ I linearly from m b ( · , t ) = 0 kg / m s for t <
50 to m b ( · , t ) = − . / m s for t >
60 aflow field is generated. The homogeneous flux boundary conditions at the burner wallΓ W = Γ C1 ∪ Γ C2 and at the symmetry line Γ S , and the Dirichlet boundary condition S ( x , t ) = S b ( x , t ) ≡ S are kept fixed. At the same time ( t ∈ [50 , y b = 0 to y b = 0 .
05. In order to start the reaction a pointheat source (modelling a glow igniter) of power density 100 000 W / m is located at x = (0 . , .
07) for t < T b , which appears in the mixed imulation of Combustion in Porous Media T b = 298 K. For t >
150 theheat transfer coefficient at the burner wall Γ W is set to h b = 1500 W / K m at Γ C1 and h b = 100 W / K m at Γ C2 . After t = 150 the boundary conditions are kept fixed, the power density of the heatsource is 0. The simulation is continued, until a numerically steady state is reached.This state depends strongly on the choice for the porosity and permeability values. φ l and φ h according to (19) Figure 3 shows the calculated distributions of the temperature T and the mass fractionof the reactant y for t = 150, t = 300, t = 500, t = 1000 and t = 5000. . . . . . . . . . . . . . . . . . . . . . . . . . . (a) (b) Figure 3.
Temperature (a) and mass fraction of reactant (b) for varying porosity
After switching off the external heat source, the reaction zone, where the values of y decrease from 0 .
05 to 0 .
0, widens in direction of the symmetry line, until nearly straightreaction front is generated. In the following, the reaction zone moves towards the inflow imulation of Combustion in Porous Media x < .
08) this movementslows down more and more. Finally (for t > λ eff is larger than in the regionwith high porosity. Therefore the heat released by the reaction is conducted fastertowards the cooled burner wall. This effect is enforced by the fact that the heat transfercoefficient at Γ C1 (no isolation) is much higher than at Γ C2 . Furthermore the reactionrate is proportional to the porosity φ . Thus there is a stronger heat production bychemical reaction in the region with high porosity. b) Results for fixed low value φ ≡ φ l according to (20) Figures 4 shows the calculated distributions of the temperature T and the mass fractionof the reactant y for t = 150, t = 200, t = 300, t = 500 and t = 1000. . . . . . . . . . . . . (a) (b) Figure 4.
Temperature (a) and mass fraction of reactant (b) for low porosity
After switching off the external heat source and starting the heat loss across the imulation of Combustion in Porous Media T is equal tothe ambient temperature. c) Results for fixed high value φ ≡ φ h according to (21) Figure 5 shows the calculated distributions of the temperature T and the mass fractionof the reactant y for t = 100, t = 200, t = 300, t = 500 and t = 1000. . . . . . . . . . . . . . . . . (a) (b) Figure 5.
Temperature (a) and mass fraction of reactant (b) for low porosity
In the beginning there is a behavior similar to case a). But now, the intensity ofthe reaction in the preheating region is not reduced due to lower porosity. Likewise, theeffective heat conductivity is not increased. Hence the movement of the reaction zonetowards the inflow boundary is not retarded strong enough. Finally, the reaction zone imulation of Combustion in Porous Media
5. Conclusion
We proposed discretization and solution methods for the equations governingcombustion in porous inert media. Decoupling the flow problem from the transportproblems, we can use different methods for the discretization of these two problems.The discretization of the flow problem is performed by the mixed finite element method.Thereby we obtain a good approximation for the mass flux, which governs the convectivetransport in the transport equations. The transport problems are discretized by a cell-centered finite volume method. The resulting nonlinear systems of equations for bothproblems are lineararized with Newton’s method, the linearized systems are solved witha multigrid algorithm. Finally both subsystems are recoupled again in a Picard iteration.Although this approach is presented on the basis of a strongly simplified model, it can beextended easily to more realistic situations. Numerical simulations with this simplifiedmodel showed that the stabilization of the reaction zone inside the porous burner canbe explained using only macroscopic quantities like porosity.
References [1] Baranger J, Maitre J F and Oudin F 1996 Connection between finite volume and mixed finiteelement methods
RAIRO Math. Mod. Num. Anal. Comput. Visual. Sci. SIAM J. Numer. Anal. Combustion and Flame
Mixed and Hybrid Finite Element Methods (New York: Springer)[6] Chen Z 1996 Equivalence between and multigrid algorithms for nonconforming and mixed methodsfor second-order elliptic problems
East-West J. Numer. Math. Finite Volume Methods in Ciarlet P G and Lions J L,eds.
Handbook of Numerical Analysis, Vol. VII (Amsterdam: North Holland)[8] Hanamura K and Echigo R 1991 An analysis of flame stabilization mechanism in radiation burners
W¨arme- und Stoff¨ubertragung Int. J. Heat Mass Transfer Math. Comput. (summitted)[11] Kwak D Y, Kwon H J and Lee S 1999 Multigrid algorithm for cell centered finite difference ontriangular meshes
Appl. Math. Comput.
J. Heat Transfer
Combust.Sci. Tech.121