An Adaptive Moving Finite Element Method for Steady Low Mach Number Compressible Combustion Problems
Zhen Sun Malte Braack Jens Lang

Atom-specific persistent homology and its application to protein flexibility analysis 
    David Bramer^1 and
Guo-Wei Wei^1,2,3 [ Address correspondences  to Guo-Wei Wei. E-mail:wei@math.msu.edu]


^1 Department of Mathematics 

Michigan State University, MI 48824, USA

^2  Department of Biochemistry and Molecular Biology

Michigan State University, MI 48824, USA 

^3 Department of Electrical and Computer Engineering 

Michigan State University, MI 48824, USA 


    Received / Accepted
===========================================================================================================================================================================================================================================================================================================================================================================================




This work surveys an r-adaptive moving mesh finite element method for the numerical solution of premixed laminar flame problems. Since the model of chemically reacting
flow involves many different modes with diverse length scales, the computation of such a problem is often extremely time-consuming. Importantly, to capture the significant characteristics of the flame structure when using detailed chemistry, a much more stringent requirement on the spatial resolution of the interior layers of some intermediate species is necessary. Here, we propose a moving mesh method in which the mesh is obtained from the solution of so-called moving mesh partial differential equations. Such equations result from the variational formulation of a minimization problem for a given target functional that characterizes the inherent difficulty in the numerical approximation of the underlying physical equations. Adaptive mesh movement has emerged as an area of intense research in mesh adaptation in the last decade. With this approach points are only allowed to be shifted in space leaving the topology of the grid unchanged. In contrast to methods with local refinement, data structure hence is unchanged and load balancing is not an issue as grid points remain on the processor where they are. We will demonstrate the high potential of moving mesh methods for effectively optimizing the distribution of grid points to reach the required resolution for chemically reacting flows with extremely thin boundary layers.


Keywords: low Mach number combustion; adaptive moving meshes;
stabilized finite elements; Rosenbrock time integrators



§ INTRODUCTION

The numerical simulation of chemically reacting flows with detailed
chemistry is still a challenging task due to the presence of
a large range of multi-scale aspects. High spatial resolution is
needed in the vicinity of the flame zone, while simultaneously
coarser meshes can be used in regions with relatively large flow
structures, e.g., downstream of the flame. Adaptive mesh refinement (AMR)
based on dynamically refining or coarsening the computational mesh
to match local features has proven to be a very efficient strategy to
perform a multi-scale combustion simulation. An overview of the basic
design concepts used to develop block-structured AMR algorithms
for the low Mach number model has been given by Bell and Day <cit.>. Two of the authors, Braack and Lang, have also been contributing to the development of AMR methods <cit.>.
However, since reacting flow applications with detailed
chemistry and transport models already tend to be significantly complex, an AMR implementation on parallel computers might quickly become prohibitively intricate.
In such a situation, adaptive moving mesh methods are considered to be an attractive
alternative since the mesh topology and hence the underlying data structures
are left unchanged. Such methods have become available through the implementation of variational principles that define a moving mesh as the solution of so-called moving mesh partial differential equations (MMPDE) introduced by Huang and
Russell, see <cit.> for a thorough overview.
Recently, adaptive moving meshes (AMM) have been successfully applied
to large eddy simulations of complex turbulent flows by Liersch, Frankenbach,
Fröhlich, and Lang <cit.>.

In this paper, we investigate the application of the AMM technique to steady low
Mach number compressible combustion problems in two spatial dimensions. Special emphasize
is put on the design of the monitor function that constitutes the heart of the MMPDE and links
the movement of the grid points to the physical solution. We will discuss two approaches that
are based on the gradient and the curvature of selected chemical components, respectively. The latter one is motivated by a standard error estimate of the interpolation error. The steady combustion model
is discretized by stabilized linear finite elements that use pressure stabilization and
streamline diffusion to enhance the main diagonal of the resulting algebraic system. A
pseudo linearly implicit time-stepping scheme is applied to solve the highly nonlinear algebraic
systems. Two examples are considered to demonstrate the potential of our AMM method: a benchmark three-component ozone decomposition and a methane flame in the practical application of a
prototype lamelle burner from the Bosch company.



§ THEORETICAL FUNDAMENTALS OF REACTIVE FLOWS




 §.§ Models for stationary chemically reacting flows

We denote the flow velocities by , the pressure by p, the temperature by T and the density by ρ. Furthermore, =(w_1,…,w_n_s)^T is defined as the
vector of mass fractions w_i of the i-th species, where i=1,…,n_s. The basic equations for a stationary reactive viscous flow express the conservation and balance laws for the total mass, momentum, energy and each species, respectively, accomplished by the equation of state:


    ∇· (ρ)    =    0, 
    ρ ( ·∇ )  + ∇ p + ∇·   =   ρ, 
    ρ·∇ w_i + ∇· ( ρ D_i ∇ w_i )    =    M_i  _i(T, ω),    i=1,…, n_s,
    
    c_p ρ·∇ T - ∇· (λ∇ T)
    - ·∇ p +  : ∇   =    -∑_i=1^n_s h_i M_i
    _i (T, ),
    
    p    =   ρ RT/,   = ( ∑_i=1^n_sw_i/M_i)^-1,


where  is the gravitational force and c_p is the heat capacity of the mixture at constant pressure. For each species, M_i is the molecular weight, h_i its (mass-)specific enthalpy, _i its molar production rate, and D_i its mass
diffusion coefficient. The viscous stress tensor  is given by


    = - μ(∇ + (∇)^T
    - 2/3(∇·) I)


with μ being the dynamic viscosity of the fluid.
Compressible flow equations admit material and acoustic waves that propagate
at velocity  and speed of sound c, respectively. In most practical
combustion systems, the flow is in the low Mach number regime, i.e., M=||/c≪ 1.
This disparity in scales can be exploited to compute combustion problems much
more efficiently. A rigorous low Mach number asymptotic analysis <cit.>
for the behaviour M→ 0 shows that the total pressure can be decomposed as


    p()  = P_th + p_hyd (),


where P_th is the constant thermodynamic pressure and
p_hyd() is the hydrodynamic pressure. With this decomposition, P_th defines
the thermodynamic state and p_hyd plays the role of a Lagrange multiplier to constrain
the flow so that the thermodynamic pressure is equilibrated everywhere in space. Hence,
compressibility effects due to chemical heat release and other thermal processes are
retained while acoustic wave propagation is entirely eliminated.
Eventually, the simplified low Mach number approximation for steady compressible reactive flows becomes <cit.>


    ∇·   =   1/T·∇ T -
    1/·∇, 
    ρ ( ·∇ )  + ∇ p_hyd -
    ∇· ( μ∇ )    =   ρ, 
    ρ·∇ w_i + ∇· ( ρ D_i ∇ w_i )    =    M_i  _i(T,),    i=1,…, n_s,
    
    c_p ρ·∇ T - ∇· (λ∇ T)
       =    -∑_i=1^n_s h_i M_i
    _i (T, ),
    ρ   =   P_th/RT,   = ( ∑_i=1^n_sw_i/M_i)^-1.


Here, we have differentiated the new equation of state in
(<ref>) to derive a
constraint on the velocity in (<ref>). As
usual, parts of the stress tensor have been absorbed into p_hyd.
We will only consider n_s-1 species
equations together with the requirement ∑_i=1^n_sw_i=1 to
keep mass conservation and the computation of
diffusion velocities consistent <cit.>. Since
only laminar flames are considered, the influence of the gravitation
in (<ref>) is also often neglected. We set
P_th=101 325 Pascal in our applications.



 §.§ Finite Element Discretization

To simplify the notation and write the system
(<ref>)-(<ref>) in a more compact form,
we extend the vector of mass fractions by w_0:=T and introduce


    ν_0 = λ,   ν_i = ρ D_i,  i=1,…,n_s,


for the diffusion coefficients of all components of =(w_0,…,w_n_s)^T.
With the new quantities


    = ρ,  =^-1∇,   l=(ρ T)^-1,
        f_0()= -∑_i=1^n_s h_i M_i_i,  
    f_i()=M_i _i, i=1,…,n_s,


and the convention that =c_p for the temperature equation
(i=0) and = otherwise,
the stationary low Mach number system can be rewritten in the condensed form


    ∇· - l ·∇ w_0 + ·    =  0, 
    
    ( ·∇ )  - ∇· (μ∇ ) + ∇ p_hyd    = ρ, 
    ·∇ w_i - ∇· ( ν_i ∇ w_i )     =  f_i(), 
    i=0,…,n_s.


For the presentation of the finite element approximation of this system, we follow
the approach presented in <cit.>, see also <cit.>.

Let Ω∈^2 be our bounded computational domain with polygonal boundary
∂Ω. Introducing the inner product and norm


    (f,g) := ∫_Ω f()g() d,   ‖ f ‖ := ( ∫_Ω| f() |^2  d)^1/2,


the Lebesgue space of all square-integrable functions on Ω, L^2(Ω),
is defined by all functions f() with ‖ f ‖ < ∞. The functions
from L^2(Ω) with square-integrable generalized first-order derivatives form
the Sobolev space H^1(Ω). We will use certain subspaces of these two spaces to
set up a variational weak formulation for (<ref>)-(<ref>),
namely


    p∈ Q=L^2(Ω)/,  ∈ H⊂( H^1(Ω))^2,
      ∈ R⊂(H^1(Ω))^n_s+1.


Note that the pressure is determined only modulo constants, which is expressed
by the special construction of the space Q. The solution triple :=(p,,)
is now an element of the space V:=Q× H× R. Testing the equations
(<ref>)-(<ref>) with a function
=(θ,,)∈ V, integrating over Ω and applying integration
by parts for the diffusive terms and the pressure gradient yields the variational
nonlinear equations


    ∈ V:    A(,) = 0   for all ∈ V


with the semi-linear form


    A(,) :=      ( ∇·, θ ) - (l ·∇ w_0, θ )
    + ( ·, θ )  + 
         (( ·∇ ) ,  ) + ( μ∇, ∇ )
    - ( p_hyd, ∇· ) - ( ρ,  )  + 
        ∑_i=0^n_s{ ( ·∇ w_i, π_i ) +
    (ν_i ∇ w_i, ∇π_i ) - ( f_i(),π_i) }.


Here, for simplicity, we have used the free-stream outflow condition and homogeneous
Dirichlet conditions. Other boundary conditions can be handled by natural modifications.

We will now describe our finite element approximation. First, we decompose Ω
into a regular partition of triangles _h={K} and
define the conforming finite element space of continuous piecewise linear functions by


    S_h := { f ∈ C^0(Ω): f_|K∈_1, K∈_h },


where _1 is the space of all polynomials of degree not larger than one.
Then, the infinite dimensional space V is approximated by a finite dimensional space
V_h=(Q_h,H_h,R_h) and the corresponding discrete approximations _h=(p_h,_h,_h) are
determined in the finite element spaces


    Q_h := S_h/,    H_h := (S_h)^2,   R_h := (S_h)^n_s+1.


The main observation is now that a simple replacement of V by V_h in (<ref>)
does not give a stable discretization. Instead, we use pressure stabilization and
streamline diffusion to enhance the main diagonal of the resulting algebraic system.
This gives the finite element approximation


    _h∈ V_h:    A(_h,) + S_h(_h,)= 0   for all ∈ V_h


with the stabilization term S_h defined by


    S_h(_h,) := c_h(_h,θ) + m_h(_h,) + ∑_i=0^n_s t_i,h(_h,π_i),


where


    c_h(_h,θ)     := ∑_K∈_hα_K (
    ( ·∇ )  + ∇ p_hyd - ρ, ∇θ ))_K,
    
    m_h(_h,)      := ∑_K∈_hδ_K (
    ( ·∇ )  + ∇ p_hyd - ρ, ·∇ ))_K,
    
    t_i,h(_h,π_i)     := ∑_K∈_hδ_K,i(
    ( ·∇ ) w_h,i - f_i(_h), ·∇π_i ))_K, i=0,…,n_s.


Here, (·,·)_K stands for the inner product over the triangle K.
The density ρ_h is a function of _h and determined from the equation of state in
(<ref>). The mesh-dependent constants α_k, β_K, and
β_K,i must be chosen carefully in order to locally reflect convection-dominated as well
as diffusion-dominated flows in an appropriate manner. We follow <cit.> and set in each element
K∈_h,


    α_K     = h^♯_K/2VRe/√(1+(Re)^2),  
    Re := h^♯_K V/μ,
    β_K     = h_K/2||Re/√(1+(Re)^2),  
    Re := h_K ||/μ,
    β_K,i    = h_K/2||Re/√(1+(Re)^2),  
    Re := h_K ||/ν_i,   i=0,…,n_s,


where V denotes a global reference velocity, h^♯_K is the diameter of the
two-dimensional ball having the same area as K, and h_K ist the length of the
element K in the direction of the local velocity .

The highly nonlinear system (<ref>) is often hard to solve without
the knowledge of a good initial guess. A common technique is the use of a homotopy
approach in order to stabilize Newton's method. This can be best realized by a pseudo
implicit time-stepping scheme. We formally add time derivatives to equation
(<ref>), resulting in


    _h(t)∈ V_h:    (P∂_t _h(t),) + A(_h(t),) + S_h(_h(t),)
    = 0   for all ∈ V_h


with a projection matrix P=diag(0,1,1). In this way, the
hydrodynamic pressure p_hyd is
still determined by a stationary equation and therefore the characteristic property
of the low-Mach number model is preserved. Since an accurate resolution of _h(t) is not
important, we solve (<ref>) with the linearly implicit Rosenbrock
method Ros3pl <cit.> employing variable step sizes and a low tolerance
until a stationary solution is obtained. This method has very good stability properties and
is suitable for the solution of differential-algebraic equations like (<ref>).




§ THE MOVING MESH METHOD



 §.§ Basic Formulation of the Moving Mesh PDE

In what follows, we will adopt the time-dependent moving mesh method described
in <cit.> to stationary combustion problems.
The mesh equation is formulated in terms of a coordinate transformation between
the original physical domain Ω and a computational domain Ω_C which
has the same topology. We denote the corresponding coordinates by
=(x_1,…,x_n)^T and =(ξ_1,…,ξ_n)^T. The time-depending mapping
(,t):Ω_C→Ω is defined by the minimizer of
the quadratic functional


    I[] = 1/2∫_Ω( ∑_i∇ξ_i^T
    ^-1∇ξ_i ) d ,


where G stands for the so-called monitor function that links
to the physical solution and is used to control the mesh density. Using a
variational approach, the minimizer is approximated by the smooth solution
of a modified gradient flow equation which reads


    ∂/∂ t = -B/τδ I/δ=
    B/τ∇· ( ^-1∇ ).


Here, τ>0 is a user specified parameter that controls the response time
of the mesh movement and B is a balance function to control the spatial movement
of the mesh points.
In order to compute the mapping (,t) directly, we interchange the roles
of dependent and independent variables. The final form of the moving mesh PDE (MMPDE)
can be expressed in the following form:


    τ∂ x/∂ t  = B ( ∑_i,j A_ij∂^2  x/∂ξ_i ∂ξ_j
    - ∑_i b_i ∂ x/∂ξ_i),


where


    A_ij := ∇ξ_i^T  G^-1∇ξ_j,         
    b_i:= ∑_j  ∇ξ_i^T ∂ G^-1/∂ξ_j∇ξ_j .


The parameter B should be determined such that all mesh points move with a
uniform time scale. This allows an easier numerical integration of the MMPDE.
The most convenient approach <cit.> to obtain a well spatially balanced
MMPDE is to scale the terms on the right-hand side of (<ref>) as follows:


    B = 1/√(∑_i (A_ii^2 +b_i^2)).


The main advantage of this approach is that the coefficients in the
right-hand side of (<ref>) become 𝒪(1). Then, the corresponding
equations are - and -scaling invariant, i.e., the MMPDE will not change if
we rescale the physical domain and the monitor function.  This allows to use the
positive parameter τ in (<ref>) to adjust the time scale of the
mesh movement. In general, a smaller τ results in a faster mesh adaptation
with respect to changes in the monitor function , while a larger τ produces
slower movement in time <cit.>. Since we
only focus on the design of a static quasi-optimal mesh to improve the approximation
quality of the stationary solution, it is sufficient to set τ=1 and to apply a
pseudo time-stepping method. For a complete specification of the coordinate transformation,
we also need to supply the MMPDE with appropriate boundary conditions. In our applications,
we always fix the mesh points on the boundary. For more advanced strategies,
we refer to <cit.>.



 §.§ Construction of the Monitor Function

A general approach for constructing monitor functions  in two
dimensions is based on its eigendecomposition,


    = λ_1_1_1^T + λ_2_2_2^T,


where _1 and _2 are mutually orthogonal normalized
eigenvectors of . The eigenvectors of  determine the
directions of the mesh adaptation, while the associated eigenvalues dictate
the intensity of the concentration of the mesh in these directions.

Given a vector-values quantity of interest (_h()) computed from
the stationary numerical solution _h()), a class of monitor functions
can be constructed by harmonic mappings <cit.>. This reads


    [              _1 = /‖‖,   _2 = _1^,; λ_1 = √(1+α‖‖^2),
        λ_2 = 1/λ_1, ]


where α denotes a user-defined intensity parameter. Typical examples for
the function  are error indicators, curvature or gradient data taken from
certain solution components. Heuristically, the numerical error is larger in regions
where the solution changes dramatically. For example, suppose a numerical scalar solution
u_h exhibits a steep front and the solution gradient is chosen to be one of the
eigenvectors, i.e., =∇ u_h, then, by using (<ref>),
the intensity of the adaptation determined by λ_1 along the normalized gradient
direction _1 is much stronger than those determined by λ_2 in
the tangential direction _2. Thus, it is expected that coordinate expansion
and compression will mainly occur in the gradient direction. However, it should be
pointed out that using the solution gradient as  may not always be the best
option for some problems. This topic will be discussed in our
applications. Similar conclusion can already be found in <cit.>.

The monitor function  is usually constructed by non-smooth operations applied
to the numerical solutions. This can result in very stiff MMPDEs. A usual remedy
is to smooth the monitor function. Let _p be a mesh point in Ω and _p
the corresponding mesh point in Ω_c. Then the following smoothing algorithm has
proven to work quite satisfactorily in practice <cit.>:


    [                          ^0(_p) :=(_p)                for all grid points _p,; ^m+1(_p) := ∫_C(_p)^m(_p()) d/| C(_p)|                for all grid points _p,;                                                               m=0,1,…, M_s-1, ]


with C(_p) ⊂Ω_C is the union of neighbouring grid cells having
_p as their common vertex. Here, the initial monitor function (_p)
is directly obtained from (<ref>) and (<ref>).
M_s is a user-specified parameter standing for the number of the smoothing cycles to
be performed. A smaller M_s gives a more accurate description of the profiles of the
monitor function, since it characterizes the exact features of the solutions. However,
this usually makes the MMPDE harder to solve.  Generally, the value of the intensity
parameter α and M_s are problem-dependent.



§ APPLICATIONS

In this section, we present two reactive flow problems which have an increasing degree of complexity. Both problems have been extensively studied in <cit.>. The first example is the ozone decomposition flame. The chemical processes are modelled with 3 species and 6 elementary chemical reactions. The motivation for this application is to demonstrate and compare the numerical results for various monitor functions. Benefits and shortcomings of the moving mesh method will also be discussed. The second example describes a methane flame in a complex geometry. The underlying chemical model contains 15 species and 84 elementary reactions. Some intermediate species have extremely thin flame layers that need to be resolved sufficiently. Furthermore, the whole process possesses both extremely fast and slow motions, implying a strict requirement on the time resolution.
All transport coefficients for viscosity, thermal conductivity and diffusion are evaluated
from kinetic models collected in the data bases of the Sandia National Laboratories <cit.>. The influence of the gravitation is neglected, i.e., we
set =0 in our applications.



 §.§ A Two-Dimensional Example for Ozone Decomposition

We use the mechanism for ozone decomposition reaction  that consists of 6 elementary reactions proposed in <cit.>:


    O_3 +  M   12⇌   O_2 + O + M
    O + O + M   34⇌   O_2  + M
    O+ O_3     56⇌   O_2 + O_2


Here, M denotes an arbitrary third body, i.e., it can be one of three considered species: Oxygen atoms O, Oxygen molecules O_2 or Ozone O_3.
Further details of the reaction mechanism based on the Arrhenius law can be found in
<cit.>.



  §.§.§ Specification of the Simulation

The geometry of the simulation is supposed to consist of two flat plates (regarded as rigid walls) with an inflow at the left side, comprising a cold mixture of ozone and oxygen molecules, and a free-stream outflow at the right side. The computational domain is defined as Ω:=]0, 0.02[ × ]0, 0.005[ and ∂Ω=Γ_in∪Γ_wall∪Γ_out. The boundary conditions are


    inflow on Γ_in:       =_in,    T=T_in,  =_in, 
    rigid walls on Γ_wall:       = 0,    T=T_wall,  ∇· = 0,  
    outflow on Γ_out:      
    -μ∇· + p  = 0,  ∇ T· = 0,
    ∇· = 0.


The initial mass fractions of species are ω_O_3=0.2, ω_O_2=0.8
and the initial temperature is


    T(x) = 298 + (800-298) e^-10^5(x_1-0.005)^2,


which is kept fixed at the inflow and at the rigid wall as Dirichlet boundary conditions.
The velocity profile on the inflow boundary is parabolic with a maximum velocity of 0.25  m/s.  The Reynolds number is approximately Re=v_cL/μ≈ 62, where v_c is the mean value of the velocity at the inlet, L stands for the length of the geometry, and μ is the viscosity of the mixture at the inflow.

A good indicator of the correct location of the flame front is the
mean value of ω_O_3 <cit.>,


    J()= 1/|Ω|∫_Ω w_O_3 d.


The profiles of the numerical approximation of w_O_3 and its
gradient are shown in Fig. <ref>.
We will use J(_h) to demonstrate the performance of our moving
mesh approach. To compare the accuracies for different meshes, we first
compute a reference solution with an adaptive, locally refined mesh with
about 45000 grid points using the Kardos software
<cit.> and provide also values for
regular meshes with different numbers of uniform grid points. The
results are collected in Tab. <ref>.







  §.§.§ A First Monitor Function

A central issue when applying a moving mesh method is to define a proper monitor function to control the distribution of the grid points. In most cases, the monitor function should ensure that the grid points are concentrated in regions where the physical solutions require a higher resolution. So the grid points on a uniform mesh should be relocated so that the requirement on the number of the mesh cells necessary to appropriately resolve a flame layer can be satisfied <cit.>.  Common practice when constructing a monitor function is to employ some heuristic choices. For problems where the solutions change dramatically within a small region, such as the flame layer of the ozone molecules ω_O_3 illustrated in Fig. <ref>, there are usually larger numerical errors in regions with large gradients. Thus the gradients of certain solution components often represent a good error indicator. The gradient of a piecewise linear function is a constant vector in each
element K∈_h. In order to improve the approximation property, we apply a standard
gradient recovery operator D^1: S_h→ (S_h)^2 as introduced in <cit.>, and choose in our first experiment


    : = D^1 w_O_3,h


as eigenvector _1 of the monitor function  to enforce a mesh adaption
in this direction. From Fig. <ref>, we can make the following observations: The contour of the flame layer becomes more distinct and can be better approximated through adjusting the locations of the grid points close to the sharp flame front.  The mesh cells are compressed to different extents along the convective direction of the flow so that an obvious anisotropic behaviour is visible.  Moreover, the mesh points are correctly concentrated in the area where the gradient is large. This is exactly the location where the mesh density takes its peak values, as illustrated in Fig. <ref>. However,
Tab. <ref> shows that the accuracy of J(_h) does not always improve with increasing grid points as one would expect. A similar behaviour has been observed in <cit.>.







In the following, we would like to discuss two reasons for this behaviour. First,
while the grid points are concentrated in regions with large gradients of
w_O_3,h, the accuracy of the other components, which also influence the
balance of O_3, might be lower due to certain local expansions of mesh lines,
resulting in a worse accuracy for J(_h). Second, although regions of large solution gradients can usually reflect areas with large numerical errors, this is not
always true for the exact error distribution. Because of the coupling of all components, using the gradient of a single component may lead to an over-concentration of mesh points in areas that have already been sufficiently resolved. The right graph in
Fig. <ref> illustrates an estimated error distribution which is based on a hierarchical error estimator proposed in <cit.>. Obviously, the error distribution differs significantly from the density function illustrated in the left graph in Fig. <ref>, which governs the concentration of the grid points in the MMPDE. The error distribution has a double-layer structure, where the numerical errors are mainly accumulated in regions with large flame curvature.  More specifically, at the downstream edge, the numerical errors are larger than those on the upstream side. The downstream edge of the flame layer is located in a transitional reaction zone, where the dynamics of the species is governed by both convection and source terms. There the intermediate species are also frequently produced and consumed. Thus, a higher resolution in this region is preferable.
Next we will construct a monitor function which is based on local error information.



  §.§.§ A Second Monitor Function

Now, we will introduce a second monitor function which is based on an approximate
interpolation error indicator. First, we recall the Sobolev space H^2(Ω) which
consists of all functions from H^1(Ω) with square-integrable generalized second-order derivatives. In what follows, we will denote by |·|_H^i, i=1,2, the usual semi-norm in
H^i defined by derivatives of i-th order only.
Let u∈ H^2(Ω)∩ C^0(Ω̅), for the moment, be a scalar solution,
u_h its linear finite element approximation, and Π_1u∈ S_h
its linear interpolant defined by u(_i)=Π_1u(_i) for all
mesh points  _i. Then a standard argument for the interpolation error


    | u - Π_1 u |_H^1(Ω)  =
    ( ∑_K ∈_h| u - Π_1 u |_H^1(K))^1/2≤  C ∑_K ∈𝒯_h h_K | u |_H^2(K),


gives the a priori error estimate for the linear finite element solution,


    ‖ u - u_h ‖_H^1(Ω) ≤  C h | u |_H^2(Ω),


where h:=max_K∈_hh_K. This, together with the observations made in
the previous section, motivates to use second derivatives of u_h to control
the local density of the moving mesh. An approximation D^2u_h∈ (S_h)^2 for
(∂_x_1x_1u,∂_x_2x_2u)^T is derived from a quadratic interpolation
of u_h <cit.>. Then, we set for our application,


    : = D^2 w_O_3,h






In Fig. <ref>, the profile of w_O_3,h and its
approximate interpolation error indicator | D^2 w_O_3,h|
are shown. We observe that this error indicators gives nearly as good information
on the local errors as the more sophisticated error estimator based on a hierarchical
basis, compare Fig. <ref>, but is much cheaper. From Tab. <ref>, we see that the results obtained by using the new
monitor function are much better than those based on the gradient of w_O_3,h.
From a comparison of both methods shown in Fig <ref>,
it becomes visible that with the new monitor function grid points are concentrated in
a larger area due to the higher second derivatives and hence numerical errors.






 §.§ A 2D Example for Methane Combustion

In this section, we consider a methane combustion problem described by
the following global reaction


    CH_4 + 2  O_2  →CO_2  + 2   H_2 O.


The reaction mechanism with 15 species and 84 elementary reactions is
taken from <cit.>. We will consider a lamella burner which is modified from a
prototype constructed by JUNKERS Bosch Thermotechnik <cit.>.
The geometry for the simulation is given in Fig. <ref>.



  §.§.§ Specification of the Simulation

A stoichiometric mixture of methane (CH_4) and air (O_2 /N_2) flows from the bottom of
the burner through a sample of slots which have a uniform width interval of 2 mm and three different heights varying from 15 mm to 11 mm. All slots have the same width of 1.5 mm. The inflow velocity is v_2=0.27 m/s. The solution is assumed to be spatially periodic. Thus it is sufficient to restrict the computational domain to three slots which define Ω, see
Fig. <ref>. The lamellae can be considered as obstacles. Dirichlet boundary conditions are used for the temperature, no-slip conditions for the velocities and Neumann boundary conditions for the species on the wall of the lamellae. To specify the temperature on the wall of the slots, we use three piecewise linear functions varying from 298 K to 393 K, 453 K, 463 K, respectively. On the cut boundary of Ω, symmetric boundary conditions are used for all the variables. For more details, we
refer to <cit.>.

This problem is characterized by the interaction of different physical processes and
largely separated scales. A simulation on uniform meshes can be very prohibitive.
We expect that moving mesh strategies designed with an appropriate monitor function will improve the numerical approximation even on relatively coarse meshes.





  §.§.§ Application of Moving Meshes

Since the underlying combustion mechanism is relatively complicated, it is
not possible to calculate a stable stationary solution on quasi-uniform
meshes with around 40000 grid points. Higher resolution is needed near the
flame front in the reaction zone and in the pre-heating zone, where the chemical
and convection-diffusion processes are strongly coupled. These regions
are unknown a priori. Therefore, we directly couple the physical PDE with the
MMPDE by an Arbitrary Lagrangian-Eulerian approach.

We first define _h(t):=_h((,t),t)=_h(,t). Under the mapping
(,t), we then transform (<ref>) into
a system involving the computational coordinates , i.e.,


    [                                       _h(t)∈ V_h:         ( P ( ∂_t _h(t) -
    (∂_t· J^-T∇̂ ) _h(t)),);                                                   +  Â(_h(t),) + Ŝ_h(_h(t),) = 0     for all ∈ V_h. ]


Here, ∇̂ denotes the gradient operator with respect to  and
J=∂/∂ is the Jacobian of the mapping (,t).
The operators Â and Ŝ_h are obtained from A and S_h by
using the identity ∇ = J^-T∇̂ and replacing _h by
_h.
The additional term (∂_t· J^-T∇̂) _h on the left
side can be viewed as a correction for the convective effects of the mesh motion.
Note that the linear finite element space V_h is now related to a
time-independent regular
partition of triangles of the computational domain Ω_C. Equation
(<ref>) is simultaneously solved with the MMPDE defined
in (<ref>) in the course of the pseudo time-stepping method. The location
of the mesh points described by the mapping (,t) is
immediately adapted to the evolving flame structure until the stationary state is
reached. Adjusting the mesh points in this way yields a stable numerical solution,
even on relatively coarse meshes. In our calculation, we start with an isotropic
quasi-uniform mesh with 29864 grid points delivered from the 2D-mesh
generator Triangle <cit.>. We first compute the
temperature-dependent flow field without combustion and start then our mesh moving
approach including all reaction terms.

For the choice of the monitor function, we study a one-dimensional simplification
of the chemical reaction process. A corresponding simulation shows that the formyl
radical HCO and its production rate have an extremely thin flame layer, see
Fig. <ref>. Moreover, both are very sensible with respect to
insufficient mesh resolution <cit.>. Hence, w_HCO is a good candidate
for the mesh design. We set = D^2 w_HCO,h, which is now updated
in each time step. In Fig. <ref>, we show the mass fraction of
CH_4, H, and HCO at the steady state. The thin reaction zone of HCO is
clearly visible. A closer look at the top of the lamellae presented
in Fig. <ref> and
Fig. <ref> reveals that
the mesh cells are strongly compressed there and aligned with respect
to the flame structure. Their areas are correspondingly reduced.





We would like to mention that the detection of the formyl radical
HCO is also of great interest since it can provide information
about the local heat release rate
which is a key parameter in the understanding of combustion processes.
Due to the low signal level, some optical diagnostic techniques, such
as HCO planar laser-induced fluorescence (PLIF), are not capable for
visualization and quantitative measurements
<cit.>. Thus, the numerical investigation
of the dynamic behaviour of this radical can reveal
more information for further research.







§ CONCLUSIONS

In this paper, we have presented a two-dimensional adaptive moving mesh
method based on
combustion-specific design criteria in order to improve the numerical
resolution of steady state premixed flames in the low Mach number regime.
We have discussed the key ideas needed to design moving meshes which
allow a relocation of mesh points without changing their connectivity.
An iterative and simultaneous moving mesh strategy that balances solution
gradients or interpolation errors over the whole spatial mesh have been
applied to two combustion problems in 2D: a ozone decomposition and
a methane flame for a lamella burner. In both simulations, a
curvature-based adaptation of the mesh points yielded an improved flame
resolution and a stable calculation of the steady state. Adaptive moving
meshes led to an enhanced simulation capability that has made it possible
to simulate realistic flames without using more sophisticated local mesh
refinement methods.



§ ACKNOWLEDGEMENTS

Z. Sun and J. Lang are supported by the Deutsche Forschungsgemeinschaft within
the Graduate School of Excellence Energy Science and Engineering (DFG GSC1070).

plain


