An Adaptive Moving Finite Element Method for Steady Low Mach Number Compressible Combustion Problems

TL;DR

This paper presents an r-adaptive moving mesh finite element method tailored for efficiently solving steady low Mach number compressible combustion problems, especially capturing thin boundary layers with high resolution.

Contribution

It introduces a moving mesh PDE approach based on variational principles to optimize mesh distribution without changing grid topology, improving computational efficiency.

Findings

Effective in resolving thin boundary layers in combustion flows

Reduces computational cost compared to fixed mesh methods

Maintains data structure and load balancing advantages

Abstract

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.

Figures

Tables

Table 1. Table 1: Values for J ​ ( 𝒖 h ) 𝐽 subscript 𝒖 ℎ J({\boldsymbol{u}}_{h}) defined in ( 49 ) for uniform meshes with different numbers of points (noP) and absolute values of the numerical errors J ​ ( 𝒆 ) = J ​ ( 𝒖 − 𝒖 h ) 𝐽 𝒆 𝐽 𝒖 subscript 𝒖 ℎ J({\boldsymbol{e}})=J({\boldsymbol{u}}-{\boldsymbol{u}}_{h}) derived from an adaptive reference solution with about 45000 45000 45000 grid points.

noP	$J\left({𝒖}_h\right)$	$\left|J\left(𝒆\right)\right|$
1105	0.03331313	1.912892 $\times {10}^{-4}$
2835	0.03347317	1.570370 $\times {10}^{-4}$
4257	0.03348598	1.844070 $\times {10}^{-5}$
adaptive	0.03350442	-

Table 2. Table 2: Results for uniform and moving meshes with the quantity of interest 𝝍 = D 1 ​ w O 3 , h 𝝍 superscript 𝐷 1 subscript 𝑤 subscript O 3 ℎ {\boldsymbol{\psi}}=D^{1}w_{\text{O}_{3},h} and parameters M s = 8 subscript 𝑀 𝑠 8 M_{s}=8 and α = 80 𝛼 80 \alpha=80 for the moving mesh design. The reference value is J ​ ( 𝒖 ) ≈ 0.03350442 𝐽 𝒖 0.03350442 J({\boldsymbol{u}})\approx 0.03350442 .

noP	$J\left({𝒖}_h\right)$, uniform	$J\left({𝒖}_h\right)$, adaptive	$\left|J\left(𝒆\right)\right|$, uniform	$\left|J\left(𝒆\right)\right|$, adaptive
1105	0.03331313	0.03356365	2.912892 $\times {10}^{-4}$	5.918997 $\times {10}^{-5}$
2835	0.03347317	0.03344742	1.570370 $\times {10}^{-4}$	5.703389 $\times {10}^{-5}$
4257	0.03348598	0.03344386	1.844070 $\times {10}^{-5}$	6.059835 $\times {10}^{-5}$

Table 3. Table 3: Computation results for moving meshes with two monitor function based on the gradient D 1 ​ w O 3 , h superscript 𝐷 1 subscript 𝑤 subscript O 3 ℎ D^{1}w_{\text{O}_{3,h}} (gradient) and on the interpolation error indicator D 2 ​ w O 3 , h superscript 𝐷 2 subscript 𝑤 subscript O 3 ℎ D^{2}w_{\text{O}_{3,h}} (interr) with M s = 8 subscript 𝑀 𝑠 8 M_{s}=8 , α = 80 𝛼 80 \alpha=80 . The reference value is J ​ ( 𝒖 ) ≈ 0.03350442 𝐽 𝒖 0.03350442 J({\boldsymbol{u}})\approx 0.03350442 .

noP	$J\left({𝒖}_h\right)$, gradient	$J\left({𝒖}_h\right)$, interr	$\left|J\left(𝒆\right)\right|$, gradient	$\left|J\left(𝒆\right)\right|$, interr
2835	0.03344742	0.03344717	5.703389 $\times {10}^{-5}$	3.272824 $\times {10}^{-5}$
4257	0.03344386	0.03504592	6.059835 $\times {10}^{-5}$	1.347509 $\times {10}^{-7}$

Equations

$$\nabla \cdot (\rho \mathbf{v})$$

$$\rho (\mathbf{v}\cdot \nabla )\mathbf{v}+\nabla p+\nabla \cdot \mathbf{\tau }$$

$$\rho \mathbf{v}\cdot \nabla w_i+\nabla \cdot (\rho D_i\nabla w_i)$$

$$c_p\rho \mathbf{v}\cdot \nabla T-\nabla \cdot (\lambda \nabla T)-\mathbf{v}\cdot \nabla p+\mathbf{\tau }:\nabla \mathbf{v}$$

$$p$$

$$\mathbf{\tau }=-\mu \left(\nabla \mathbf{v}+(\nabla \mathbf{v}{)}^T-\frac{2}{3}(\nabla \cdot \mathbf{v})\text{I}\right)$$

$$p(\mathbf{x})=P_{th}+p_{hyd}(\mathbf{x}),$$

$$\nabla \cdot \mathbf{v}$$

$$\rho (\mathbf{v}\cdot \nabla )\mathbf{v}+\nabla p_{hyd}-\nabla \cdot (\mu \nabla \mathbf{v})$$

$$\rho \mathbf{v}\cdot \nabla w_i+\nabla \cdot (\rho D_i\nabla w_i)$$

$$c_p\rho \mathbf{v}\cdot \nabla T-\nabla \cdot (\lambda \nabla T)$$

$$\rho$$

$${\nu }_0=\lambda ,{\nu }_i=\rho D_i,i=1,\dots ,n_s,$$

$$\mathbf{\beta }=\rho \mathbf{v},\mathbf{m}={\overline{M}}^{-1}\nabla \overline{M},l=(\rho T{)}^{-1},$$

$$f_0(\mathbf{w})=-\sum _{i=1}^{n_s}{h}_i{M}_i{\dot{w}}_i,f_i(\mathbf{w})=M_i{\dot{w}}_i,i=1,\dots ,n_s,$$

$$\nabla \cdot \mathbf{v}-l\mathbf{\beta }\cdot \nabla w_0+\mathbf{v}\cdot \mathbf{m}$$

$$(\mathbf{\beta }\cdot \nabla )\mathbf{v}-\nabla \cdot (\mu \nabla \mathbf{v})+\nabla p_{hyd}$$

$$\bar{\mathbf{\beta }}\cdot \nabla w_i-\nabla \cdot ({\nu }_i\nabla w_i)$$

$$(f,g):={\int }_{\Omega }f(\mathbf{x})g(\mathbf{x})d\mathbf{x},\Vert f\Vert :={\left({\int }_{\Omega }|f(\mathbf{x}){|}^{2}d\mathbf{x}\right)}^{1/2},$$

$$p\in Q=L^2(\Omega )/\mathbb{R},\mathbf{v}\in H\subset {\left(H^1(\Omega )\right)}^2,\mathbf{w}\in R\subset {\left(H^1(\Omega )\right)}^{n_s+1}.$$

$$\mathbf{u}\in V:A(\mathbf{u},\mathbf{\varphi })=0\text{for all }\mathbf{\varphi }\in V$$

$$A(\mathbf{u},\mathbf{\varphi }):=$$

$$((\mathbf{\beta }\cdot \nabla )\mathbf{v},\mathbf{\chi })+(\mu \nabla \mathbf{v},\nabla \mathbf{\chi })-(p_{hyd},\nabla \cdot \mathbf{\chi })-(\rho \mathbf{g},\mathbf{\chi })+$$

$$\sum _{i=0}^{n_s}\{(\bar{\mathbf{\beta }}\cdot \nabla w_i,{\pi }_i)+({\nu }_i\nabla w_i,\nabla {\pi }_i)-(f_i(\mathbf{w}),{\pi }_i)\}.$$

$$S_h:=\{f\in C^0(\overline{\Omega }):f_{|K}\in {\mathbb{P}}_1,K\in {\mathcal{T}}_h\},$$

$$Q_h:=S_h/\mathbb{R},H_h:={\left(S_h\right)}^2,R_h:={\left(S_h\right)}^{n_s+1}.$$

$${\mathbf{u}}_h\in V_h:A({\mathbf{u}}_h,\mathbf{\varphi })+S_h({\mathbf{u}}_h,\mathbf{\varphi })=0\text{for all }\mathbf{\varphi }\in V_h$$

$$S_h({\mathbf{u}}_h,\mathbf{\varphi }):=c_h({\mathbf{u}}_h,\theta )+m_h({\mathbf{u}}_h,\mathbf{\chi })+\sum _{i=0}^{n_s}{t}_{i,h}({\mathbf{u}}_h,{\pi }_i),$$

$$c_h({\mathbf{u}}_h,\theta )$$

$$m_h({\mathbf{u}}_h,\mathbf{\chi })$$

$$t_{i,h}({\mathbf{u}}_h,{\pi }_i)$$

$${\alpha }_K$$

$${\beta }_K$$

$${\beta }_{K,i}$$

$${\mathbf{u}}_h(t)\in V_h:(P{\partial }_t{\mathbf{u}}_h(t),\mathbf{\varphi })+A({\mathbf{u}}_h(t),\mathbf{\varphi })+S_h({\mathbf{u}}_h(t),\mathbf{\varphi })=0\text{for all }\mathbf{\varphi }\in V_h$$

$$I[\mathbf{\xi }]=\frac{1}{2}{\int }_{\Omega }\left(\sum _i\nabla {\xi }_i^T{\mathbf{G}}^{-1}\nabla {\xi }_i\right)d\mathbf{x},$$

$$\frac{\partial \mathbf{\xi }}{\partial t}=-\frac{B}{\tau }\frac{\delta I}{\delta \mathbf{\xi }}=\frac{B}{\tau }\nabla \cdot ({\mathbf{G}}^{-1}\nabla \mathbf{\xi }).$$

$$\tau \frac{\partial \mathbf{x}}{\partial t}=B\left(\sum _{i,j}{A}_{ij}\frac{{\partial }^2\mathbf{x}}{\partial {\xi }_i\partial {\xi }_j}-\sum _i{b}_i\frac{\partial \mathbf{x}}{\partial {\xi }_i}\right),$$

$$A_{ij}:=\nabla {\xi }_i^T{\mathbf{G}}^{-1}\nabla {\xi }_j,b_i:=\sum _j\nabla {\xi }_i^T\frac{\partial {\mathbf{G}}^{-1}}{\partial {\xi }_j}\nabla {\xi }_j.$$

Full text

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

1 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 [2]. Two of the authors, Braack and Lang, have also been contributing to the development of AMR methods [1, 5, 10, 16]. 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 [13] 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 [21].

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.

2 Theoretical Fundamentals of Reactive Flows

2.1 Models for stationary chemically reacting flows

We denote the flow velocities by ${\boldsymbol{v}}$, the pressure by $p$, the temperature by $T$ and the density by $\rho$. Furthermore, ${\boldsymbol{w}}=(w_{1},\ldots,w_{n_{s}})^{T}$ is defined as the vector of mass fractions $w_{i}$ of the $i$-th species, where $i=1,\ldots,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:

[TABLE]

where ${\boldsymbol{g}}$ 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, $\dot{w}_{i}$ its molar production rate, and $D_{i}$ its mass diffusion coefficient. The viscous stress tensor ${\boldsymbol{\tau}}$ is given by

[TABLE]

with $\mu$ being the dynamic viscosity of the fluid. Compressible flow equations admit material and acoustic waves that propagate at velocity ${\boldsymbol{v}}$ and speed of sound $c$, respectively. In most practical combustion systems, the flow is in the low Mach number regime, i.e., $M=|{\boldsymbol{v}}|/c\ll 1$. This disparity in scales can be exploited to compute combustion problems much more efficiently. A rigorous low Mach number asymptotic analysis [22] for the behaviour $M\rightarrow 0$ shows that the total pressure can be decomposed as

[TABLE]

where $P_{th}$ is the constant thermodynamic pressure and $p_{hyd}({\boldsymbol{x}})$ 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 [4]

[TABLE]

Here, we have differentiated the new equation of state in (12) to derive a constraint on the velocity in (8). 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 $\sum_{i=1}^{n_{s}}w_{i}=1$ to keep mass conservation and the computation of diffusion velocities consistent [26]. Since only laminar flames are considered, the influence of the gravitation in (9) is also often neglected. We set $P_{th}=101\,325\,$\mathrm{P}\mathrm{a}\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}$$ in our applications.

2.2 Finite Element Discretization

To simplify the notation and write the system (8)-(12) in a more compact form, we extend the vector of mass fractions by $w_{0}:=T$ and introduce

[TABLE]

for the diffusion coefficients of all components of ${\boldsymbol{w}}=(w_{0},\ldots,w_{n_{s}})^{T}$. With the new quantities

[TABLE]

and the convention that $\bar{{\boldsymbol{\beta}}}=c_{p}{\boldsymbol{\beta}}$ for the temperature equation ($i=0$) and $\bar{{\boldsymbol{\beta}}}={\boldsymbol{\beta}}$ otherwise, the stationary low Mach number system can be rewritten in the condensed form

[TABLE]

For the presentation of the finite element approximation of this system, we follow the approach presented in [4, Sec. 2], see also [3].

Let $\Omega\in\mathbb{R}^{2}$ be our bounded computational domain with polygonal boundary $\partial\Omega$. Introducing the inner product and norm

[TABLE]

the Lebesgue space of all square-integrable functions on $\Omega$, $L^{2}(\Omega)$, is defined by all functions $f({\boldsymbol{x}})$ with $\|f\|<\infty$. The functions from $L^{2}(\Omega)$ with square-integrable generalized first-order derivatives form the Sobolev space $H^{1}(\Omega)$. We will use certain subspaces of these two spaces to set up a variational weak formulation for (16)-(18), namely

[TABLE]

Note that the pressure is determined only modulo constants, which is expressed by the special construction of the space $Q$. The solution triple ${\boldsymbol{u}}:=(p,{\boldsymbol{v}},{\boldsymbol{w}})$ is now an element of the space $V:=Q\times H\times R$. Testing the equations (16)-(18) with a function ${\boldsymbol{\phi}}=(\theta,{\boldsymbol{\chi}},{\boldsymbol{\pi}})\in V$, integrating over $\Omega$ and applying integration by parts for the diffusive terms and the pressure gradient yields the variational nonlinear equations

[TABLE]

with the semi-linear form

[TABLE]

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 $\overline{\Omega}$ into a regular partition of triangles ${\cal T}_{h}=\{K\}$ and define the conforming finite element space of continuous piecewise linear functions by

[TABLE]

where $\mathbb{P}_{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 ${\boldsymbol{u}}_{h}=(p_{h},{\boldsymbol{v}}_{h},{\boldsymbol{w}}_{h})$ are determined in the finite element spaces

[TABLE]

The main observation is now that a simple replacement of $V$ by $V_{h}$ in (21) 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

[TABLE]

with the stabilization term $S_{h}$ defined by

[TABLE]

where

[TABLE]

Here, $(\cdot,\cdot)_{K}$ stands for the inner product over the triangle $K$. The density $\rho_{h}$ is a function of ${\boldsymbol{w}}_{h}$ and determined from the equation of state in (12). The mesh-dependent constants $\alpha_{k}$, $\beta_{K}$, and $\beta_{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 [17] and set in each element $K\in{\cal T}_{h}$,

[TABLE]

where $V$ denotes a global reference velocity, $h^{\sharp}_{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 ${\boldsymbol{v}}$.

The highly nonlinear system (25) 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 (25), resulting in

[TABLE]

with a projection matrix $P=\text{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 ${\boldsymbol{u}}_{h}(t)$ is not important, we solve (33) with the linearly implicit Rosenbrock method Ros3pl [20] 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 (33).

3 The Moving Mesh Method

3.1 Basic Formulation of the Moving Mesh PDE

In what follows, we will adopt the time-dependent moving mesh method described in [12, 13] to stationary combustion problems. The mesh equation is formulated in terms of a coordinate transformation between the original physical domain $\Omega$ and a computational domain $\Omega_{C}$ which has the same topology. We denote the corresponding coordinates by ${\boldsymbol{x}}=(x_{1},\dots,x_{n})^{T}$ and ${\boldsymbol{\xi}}=(\xi_{1},\dots,\xi_{n})^{T}$. The time-depending mapping ${\boldsymbol{x}}({\boldsymbol{\xi}},t):\Omega_{C}\rightarrow\Omega$ is defined by the minimizer of the quadratic functional

[TABLE]

where $\boldsymbol{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

[TABLE]

Here, $\tau>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 ${\boldsymbol{x}}({\boldsymbol{\xi}},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:

[TABLE]

where

[TABLE]

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 [11] to obtain a well spatially balanced MMPDE is to scale the terms on the right-hand side of (36) as follows:

[TABLE]

The main advantage of this approach is that the coefficients in the right-hand side of (36) become $\mathcal{O}(1)$. Then, the corresponding equations are ${\boldsymbol{x}}$- and ${\boldsymbol{G}}$-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 $\tau$ in (36) to adjust the time scale of the mesh movement. In general, a smaller $\tau$ results in a faster mesh adaptation with respect to changes in the monitor function ${\boldsymbol{G}}$, while a larger $\tau$ produces slower movement in time [7, 8]. 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 $\tau=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 [11].

3.2 Construction of the Monitor Function

A general approach for constructing monitor functions ${\boldsymbol{G}}$ in two dimensions is based on its eigendecomposition,

[TABLE]

where $\boldsymbol{\mathit{v}}_{1}$ and $\boldsymbol{\mathit{v}}_{2}$ are mutually orthogonal normalized eigenvectors of ${\boldsymbol{G}}$. The eigenvectors of ${\boldsymbol{G}}$ 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 ${\boldsymbol{\psi}}({\boldsymbol{u}}_{h}({\boldsymbol{x}}))$ computed from the stationary numerical solution ${\boldsymbol{u}}_{h}({\boldsymbol{x}}))$, a class of monitor functions can be constructed by harmonic mappings [6]. This reads

[TABLE]

where $\alpha$ denotes a user-defined intensity parameter. Typical examples for the function ${\boldsymbol{\psi}}$ 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., ${\boldsymbol{\psi}}=\nabla u_{h}$, then, by using (40), the intensity of the adaptation determined by $\lambda_{1}$ along the normalized gradient direction $\boldsymbol{\mathit{v}}_{1}$ is much stronger than those determined by $\lambda_{2}$ in the tangential direction $\boldsymbol{\mathit{v}}_{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 ${\boldsymbol{\psi}}$ 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 [7].

The monitor function ${\boldsymbol{G}}$ 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 ${\boldsymbol{x}}_{p}$ be a mesh point in $\Omega$ and ${\boldsymbol{\xi}}_{p}$ the corresponding mesh point in $\Omega_{c}$. Then the following smoothing algorithm has proven to work quite satisfactorily in practice [11, 19]:

[TABLE]

with $C({\boldsymbol{\xi}}_{p})\subset\Omega_{C}$ is the union of neighbouring grid cells having ${\boldsymbol{\xi}}_{p}$ as their common vertex. Here, the initial monitor function ${\boldsymbol{G}}({\boldsymbol{x}}_{p})$ is directly obtained from (39) and (40). $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 $\alpha$ and $M_{s}$ are problem-dependent.

4 Applications

In this section, we present two reactive flow problems which have an increasing degree of complexity. Both problems have been extensively studied in [3, 4]. 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 [14]. The influence of the gravitation is neglected, i.e., we set ${\boldsymbol{g}}=0$ in our applications.

4.1 A Two-Dimensional Example for Ozone Decomposition

We use the mechanism for ozone decomposition reaction that consists of $6$ elementary reactions proposed in [29]:

[TABLE]

Here, M denotes an arbitrary third body, i.e., it can be one of three considered species: Oxygen atoms O, Oxygen molecules $\text{O}_{2}$ or Ozone $\text{O}_{3}$. Further details of the reaction mechanism based on the Arrhenius law can be found in [4, Tab. A.1].

4.1.1 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 $\Omega:=]0,0.02[\times]0,0.005[$ and $\partial\Omega=\Gamma_{\text{in}}\cup\Gamma_{\text{wall}}\cup\Gamma_{\text{out}}$. The boundary conditions are

[TABLE]

The initial mass fractions of species are $\omega_{\text{O}_{3}}=0.2$, $\omega_{\text{O}_{2}}=0.8$ and the initial temperature is

[TABLE]

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_{c}L/\mu\approx 62$, where $v_{c}$ is the mean value of the velocity at the inlet, $L$ stands for the length of the geometry, and ${\mu}$ 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 $\omega_{\text{O}_{3}}$ [4],

[TABLE]

The profiles of the numerical approximation of $w_{\text{O}_{3}}$ and its gradient are shown in Fig. 1. We will use $J({\boldsymbol{u}}_{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 [9, 18] and provide also values for regular meshes with different numbers of uniform grid points. The results are collected in Tab. 1.

4.1.2 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 [24]. 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 $\omega_{\text{O}_{3}}$ illustrated in Fig. 1, 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\in{\cal T}_{h}$. In order to improve the approximation property, we apply a standard gradient recovery operator $D^{1}:S_{h}\rightarrow(S_{h})^{2}$ as introduced in [30], and choose in our first experiment

[TABLE]

as eigenvector ${\boldsymbol{v}}_{1}$ of the monitor function ${\boldsymbol{G}}$ to enforce a mesh adaption in this direction. From Fig. 2, 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. 3. However, Tab. 2 shows that the accuracy of $J({\boldsymbol{u}}_{h})$ does not always improve with increasing grid points as one would expect. A similar behaviour has been observed in [11].

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_{\text{O}_{3},h}$, the accuracy of the other components, which also influence the balance of $\text{O}_{3}$, might be lower due to certain local expansions of mesh lines, resulting in a worse accuracy for $J({\boldsymbol{u}}_{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. 3 illustrates an estimated error distribution which is based on a hierarchical error estimator proposed in [18]. Obviously, the error distribution differs significantly from the density function illustrated in the left graph in Fig. 3, 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.

4.1.3 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}(\Omega)$ which consists of all functions from $H^{1}(\Omega)$ with square-integrable generalized second-order derivatives. In what follows, we will denote by $|\cdot|_{H^{i}},i=1,2$, the usual semi-norm in $H^{i}$ defined by derivatives of $i$-th order only. Let $u\in H^{2}(\Omega)\cap C^{0}(\bar{\Omega})$, for the moment, be a scalar solution, $u_{h}$ its linear finite element approximation, and $\mathnormal{\Pi}_{1}u\in S_{h}$ its linear interpolant defined by $u({\boldsymbol{x}}_{i})=\mathnormal{\Pi}_{1}u({\boldsymbol{x}}_{i})$ for all mesh points ${\boldsymbol{x}}_{i}$. Then a standard argument for the interpolation error

[TABLE]

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

[TABLE]

where $h:=\max_{K\in{\cal T}_{h}}h_{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^{2}u_{h}\in(S_{h})^{2}$ for $(\partial_{x_{1}x_{1}}u,\partial_{x_{2}x_{2}}u)^{T}$ is derived from a quadratic interpolation of $u_{h}$ [7, Sec. 4.2]. Then, we set for our application,

[TABLE]

In Fig. 4, the profile of $w_{\text{O}_{3},h}$ and its approximate interpolation error indicator $|D^{2}w_{\text{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. 3, but is much cheaper. From Tab. 3, we see that the results obtained by using the new monitor function are much better than those based on the gradient of $w_{\text{O}_{3},h}$. From a comparison of both methods shown in Fig 5, 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.

4.2 A 2D Example for Methane Combustion

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

[TABLE]

The reaction mechanism with $15$ species and $84$ elementary reactions is taken from [27]. We will consider a lamella burner which is modified from a prototype constructed by JUNKERS Bosch Thermotechnik [23]. The geometry for the simulation is given in Fig. 6.

4.2.1 Specification of the Simulation

A stoichiometric mixture of methane (CH4) and air (O${}_{2}/$N2) 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 $\Omega$, see Fig. 6. 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 $\Omega$, symmetric boundary conditions are used for all the variables. For more details, we refer to [3, 4, 28].

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.

4.2.2 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 ${\boldsymbol{\hat{u}}}_{h}(t):={\boldsymbol{\hat{u}}}_{h}({\boldsymbol{x}}({\boldsymbol{\xi}},t),t)={\boldsymbol{\hat{u}}}_{h}({\boldsymbol{\xi}},t)$. Under the mapping ${\boldsymbol{x}}({\boldsymbol{\xi}},t)$, we then transform (33) into a system involving the computational coordinates ${\boldsymbol{\xi}}$, i.e.,

[TABLE]

Here, $\hat{\nabla}$ denotes the gradient operator with respect to ${\boldsymbol{\xi}}$ and $J=\partial{\boldsymbol{x}}/\partial{\boldsymbol{\xi}}$ is the Jacobian of the mapping ${\boldsymbol{x}}({\boldsymbol{\xi}},t)$. The operators $\hat{A}$ and $\hat{S}_{h}$ are obtained from $A$ and $S_{h}$ by using the identity $\nabla=J^{-T}\hat{\nabla}$ and replacing ${\boldsymbol{u}}_{h}$ by ${\boldsymbol{\hat{u}}}_{h}$. The additional term $(\partial_{t}{\boldsymbol{x}}\cdot J^{-T}\hat{\nabla})\,{\boldsymbol{\hat{u}}}_{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 $\Omega_{C}$. Equation (55) is simultaneously solved with the MMPDE defined in (36) in the course of the pseudo time-stepping method. The location of the mesh points described by the mapping ${\boldsymbol{x}}({\boldsymbol{\xi}},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 [25]. 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. 7. Moreover, both are very sensible with respect to insufficient mesh resolution [28]. Hence, $w_{\text{HCO}}$ is a good candidate for the mesh design. We set ${\boldsymbol{\psi}}=D^{2}w_{\text{HCO},h}$, which is now updated in each time step. In Fig. 8, we show the mass fraction of CH4, 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. 9 and Fig. 10 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 [15]. Thus, the numerical investigation of the dynamic behaviour of this radical can reveal more information for further research.

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

6 Acknowledgements

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