# Verification of a low Mach variable-density Navier-Stokes solver for turbulent combustion

R. Mullyadzhanov, E. Palkin, B. Nićeno, L. Vervisch, K. Hanjalić

► To cite this version:

R. Mullyadzhanov, E. Palkin, B. Nićeno, L. Vervisch, K. Hanjalić. Verification of a low Mach variable-density Navier-Stokes solver for turbulent combustion. Journal of Physics: Conference Series, 2016, 754 (6), pp.062005. 10.1088/1742-6596/754/6/062005. hal-01611262

HAL Id: hal-01611262

https://hal.science/hal-01611262

Submitted on 29 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
PAPER • OPEN ACCESS

Verification of a low Mach variable-density Navier-
Stokes solver for turbulent combustion

To cite this article: R Mullyadzhanov et al 2016 *J. Phys.: Conf. Ser.* **754** 062005

View the [article online](#) for updates and enhancements.

You may also like

- Analysis of the process of loosening the rocks with different strength properties using the undercutting bolts
M Siegmund and J Jonak

- The Students' Errors in Answering Geometric Tests with Newman Procedures
Thamrin Tayeb, Andi Dian Angriani, Sri Ratu Humaerah et al.

- Evaluating math test of students through different item functioning: Who's taken benefit and harmed?
S Solihatun, I B Rangka, J Siregar et al.
# Verification of a low Mach variable-density Navier-Stokes solver for turbulent combustion

R Mullyadzhanov¹,², E Palkin¹,², B Ničeno³, L Vervisch⁴, K Hanjalić²,⁵

¹ Institute of Thermophysics SB RAS, Novosibirsk, Russia

² Novosibirsk State University, Novosibirsk, Russia

³ Paul Scherrer Institute, Villigen, Switzerland

⁴ CORIA-CNRS and INSA Rouen, Normandie Université, Rouen, France

⁵ Delft University of Technology, Delft, The Netherlands

E-mail: rustammul@gmail.com

**Abstract.** We describe the low Mach variable-density Navier-Stokes numerical iterative solution procedure implemented in the finite-volume unstructured T-FlowS code. As the test cases we use a number of analytic manufactured solutions and Rayleigh-Taylor instability problem from the literature for algorithm verification purposes. The tests show that the code is second-order accurate in agreement with the spatial discretization scheme. We outline the recent combustion ADEF model implemented in the program.

## 1. Introduction

During the last few decades rapid development of the computational technologies stimulated the research using computational codes for fundamental and applied studies. Increasing complexity of codes makes the verification process essential for research purposes [1]. In the present work we assess the accuracy of the implemented algorithm in the variable-density Navier-Stokes solver related to the combustion problems where quite a broad range of problems can serve as a benchmark [2–9]. As the test cases we use the analytic manufactured solutions derived by Shunn et al [10], the Rayleigh-Taylor instability problem employed by Desjardins et al [11] and describe the combustion modeling framework [12] used in Large-eddy simulations of the Cambridge stratified burner [13].

## 2. Governing equations and iterative procedure

At low Mach numbers variable-density reacting flows can be described by the following conservation equations for mass, momentum and scalars together with a suitable equation of state:

$$ \frac{\partial \rho}{\partial t} + \frac{\partial \rho u_j}{\partial x_j} = \dot{Q}_{\rho}, \qquad (1) $$

$$ \frac{\partial \rho u_i}{\partial t} + \frac{\partial \rho u_i u_j}{\partial x_j} = -\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j} (2\mu S_{ij}) + \dot{Q}_{u_i}, \qquad (2) $$

<figure><img src="image_2.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
$$ \frac{\partial \rho \phi_k}{\partial t} + \frac{\partial \rho \phi_k u_j}{\partial x_j} = \frac{\partial}{\partial x_j} \left( \rho \alpha_k \frac{\partial \phi_k}{\partial x_j} \right) + \dot{Q}_{\phi_k}, \quad (3) $$

$$ \rho = f(\phi_1, \phi_2, \dots, \phi_k, \dots), \quad (4) $$

where

$$ S_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) - \frac{1}{3} \delta_{ij} \frac{\partial u_l}{\partial x_l} \quad (5) $$

is the rate-of-strain tensor. The variables $\rho$, $u_i$, $p$ and $\phi_k$ denote the density, velocity vector components, pressure and additional scalar fields. The dynamic viscosity and kinematic diffusivity coefficients are denoted by $\mu$ and $\alpha_k$. The source terms $\dot{Q}_{\rho}$, $\dot{Q}_{u_i}$ and $\dot{Q}_{\phi_k}$ in Eqs. (1)-(3) are used to construct the analytic manufactured solutions or denote the subgrid-scale terms appearing in the framework of Large-eddy simulations (LES). In the present work we consider only one additional scalar $\phi$, thus, further we omit the subscript index 'k'.

The following iteration procedure similar to the one described by [14] is implemented into the finite-volume T-FlowS code [15] employing unstructured cell-centered collocated grids and featuring second-order accuracy in time and space. Below the superscript $n$ and $n-1$ refers to solutions that are known from two previous time levels, the subscript 0 indicates the initial guess or the current provisional value of some quantity.

**Step 0:** The solution from the previous time step $n$ for momentum, pressure and scalar are used as the initial best guess:

$$ (\rho u_i)_0^{n+1} = (\rho u_i)_n, \quad p_0^{n+1} = p^n, \quad \phi_0^{n+1} = \phi^n. \quad (6) $$

The predictor $\rho_0^{n+1}$ for the density is calculated from the continuity equation, Eq. (1), integrated over a control volume (CV):

$$ \int_{CV} \frac{\partial \rho}{\partial t} dV = \frac{\Delta V}{\Delta t} \left( \frac{3}{2} \rho_0^{n+1} - 2\rho^n - \frac{1}{2} \rho^{n-1} \right) = - \int_{CV} \frac{\partial \rho u_j}{\partial x_j} dV = - \int_{faces} \rho u_j n_j dS, \quad (7) $$

where the Gauss theorem is used to go from the volume to the surface integration. The time derivative is approximated with a three-point backward-difference scheme (second order).

**Step 1:** The scalar equation is advanced so that a better estimate for the density can be obtained early in the iteration process. Time advancing Eq. (3) yields $(\rho\phi)^{n+1}$, from which a provisional estimate for $\phi_0^{n+1}$ is obtained using the current density predictor:

$$ \phi_0^{n+1} = (\rho\phi)^{n+1} / \rho_0^{n+1}. \quad (8) $$

**Step 2:** Update the density from the equation of state using the provisional scalar values:

$$ \rho^{n+1} = f(\phi_0^{n+1}). \quad (9) $$

**Step 3:** Update the scalar based on the new density:

$$ \phi^{n+1} = (\rho\phi)^{n+1} / \rho^{n+1}. \quad (10) $$

**Step 4:** Advance in time the momentum equations, Eq. (2), to obtain provisional values for the momentum components, $g_i = \rho u_i$. Compute provisional velocity field:

$$ u_{i,0}^{n+1} = g_{i,0}^{n+1} / \rho^{n+1}. \quad (11) $$
Step 5: A Poisson equation is solved to satisfy the continuity equation adjusting the values of momentum and pressure. We consider additive corrections, $\delta g_i$ and $\delta p$, to the momentum and pressure such that

$$g_i^{n+1} = g_{i,0}^{n+1} + \delta g_i^{n+1}, \quad p^{n+1} = p_0^{n+1} + \delta p^{n+1}, \qquad (12)$$

We then substitute these expressions into the momentum equation. Following Pierce [14] we group all other terms into a generic right-hand-side term, $R_i$,

$$\frac{\partial g_i^{n+1}}{\partial t} = \frac{\partial(g_{i,0}^{n+1} + \delta g_i^{n+1})}{\partial t} = -\frac{\partial(p_0^{n+1} + \delta p^{n+1})}{\partial x_i} + R_i. \qquad (13)$$

Provided that the following equation has been satisfied on the previous step

$$\frac{\partial g_{i,0}^{n+1}}{\partial t} = -\frac{\partial p_0^{n+1}}{\partial x_i} + R_i, \qquad (14)$$

we have to satisfy the rest of the equation

$$\frac{\partial \delta g_i^{n+1}}{\partial t} = -\frac{\partial \delta p^{n+1}}{\partial x_i}. \qquad (15)$$

The time derivative is approximated using the three-point backward-difference to be consistent with Eq. (7). Taking the divergence of this expression and integrating over a control volume, we arrive to

$$\frac{1}{\Delta t} \int_{CV} \frac{\partial}{\partial x_i} \left( \frac{3}{2} \delta g_i^{n+1} - 2 \delta g_i^n + \frac{1}{2} \delta g_i^{n-1} \right) dV = - \int_{CV} \frac{\partial^2 \delta p^{n+1}}{\partial x_i^2} dV. \qquad (16)$$

Note that $\delta g_i^n = \delta g_i^{n-1} = 0$ provided the algorithm has satisfied both the momentum and continuity equations on previous time steps. The derivative of the momentum correction can be expressed as the residual of the continuity equation:

$$\frac{\partial \delta g_i^{n+1}}{\partial x_i} = \frac{\partial(g_i^{n+1} - g_{i,0}^{n+1})}{\partial x_i} = -\frac{\partial \rho^{n+1}}{\partial t} - \frac{\partial g_{i,0}^{n+1}}{\partial x_i} = -\delta(\text{cont.})^{n+1}. \qquad (17)$$

Going from volume to surface integrals, we obtain

$$\frac{3}{2\Delta t} \int_{CV} \delta(\text{cont.})^{n+1} dV = \int_{\text{faces}} \frac{\partial \delta p^{n+1}}{\partial x_i} n_i dS \qquad (18)$$

which is solved for the pressure correction derivatives in the cell centers.

Step 6: The momentum components, velocity and pressure at the cell centers are updated:

$$g_i^{n+1} = g_{i,0}^{n+1} - \Delta t \frac{\partial \delta p^{n+1}}{\partial x_i}, \quad u_i^{n+1} = g_i^{n+1} / \rho^{n+1}, \qquad (19)$$

$$p^{n+1} = p_0^{n+1} + \delta p^{n+1}. \qquad (20)$$

In the next section we describe a set of test cases to verify the described algorithm.
<table><caption>**Table 1.** Parameters of problems introduced by Shunn et al [10]</caption><thead><tr><th>Problem 1</th><th>Value</th><th>Problem 2</th><th>Value</th><th>Problem 3</th><th>Value</th></tr></thead><tbody><tr><td>ρ<sub>0</sub></td><td>20</td><td>ρ<sub>0</sub></td><td>20</td><td>ρ<sub>0</sub></td><td>5</td></tr><tr><td>ρ<sub>1</sub></td><td>1</td><td>ρ<sub>1</sub></td><td>1</td><td>ρ<sub>1</sub></td><td>1</td></tr><tr><td>k<sub>1</sub></td><td>4</td><td>u<sub>F</sub></td><td>1</td><td>k = ω</td><td>2</td></tr><tr><td>k<sub>2</sub></td><td>2</td><td>v<sub>F</sub></td><td>0.5</td><td>u<sub>F</sub> = v<sub>F</sub></td><td>0.5</td></tr><tr><td>w<sub>0</sub></td><td>5</td><td>a</td><td>0.2</td><td>ρα<sub>φ</sub> = μ</td><td>0.001</td></tr><tr><td>ρα<sub>φ</sub> = μ</td><td>0.03</td><td>b</td><td>20</td><td></td><td></td></tr><tr><td></td><td></td><td>k</td><td>4π</td><td></td><td></td></tr><tr><td></td><td></td><td>ω</td><td>1.5</td><td></td><td></td></tr><tr><td></td><td></td><td>ρα<sub>φ</sub> = μ</td><td>0.001</td><td></td><td></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Profiles of $\phi$, $\rho$ and $u$ for Problem 1 and different time instants (increment $\Delta t = 0.2$).</figcaption></figure>

### 3. Laminar test cases

#### 3.1. Problem 1: one-dimensional diffusion

We consider a one-dimensional problem described by the following relations:

$$ \phi(x,t) = \frac{\exp(-k_1 t) - \cosh(w_0 x \exp(-k_2 t))}{\exp(-k_1 t)(1 - \frac{\rho_0}{\rho_1}) - \cosh(w_0 x \exp(-k_2 t))}, \quad (21) $$

$$ \rho(x,t) = \left( \frac{\phi(x,t)}{\rho_1} + \frac{1 - \phi(x,t)}{\rho_2} \right)^{-1}, \quad (22) $$

$$ u(x,t) = 2k_2 \exp(-k_1 t) \frac{\Delta \rho}{\rho(x,t)} \left( \frac{\hat{u}x}{\hat{u}^2 + 1} + \frac{\left(\frac{k_1}{k_2} - 1\right)\left(\arctan \hat{u} - \frac{\pi}{4}\right)}{w_0 \exp(-k_2 t)} \right), \quad (23) $$

where $\hat{u} = \exp(w_0 x \exp(-k_2 t))$ and $\Delta\rho = \rho_0 - \rho_1$, $w_0$, $k_1$ and $k_2$ are constants (see Table 1). The corresponding source terms in Eqs. (1)-(3) are explicitly derived by Shunn et al [10] (see their Appendices). The computational domain for this problem is $0 \le x \le 2$ and $0 \le t \le 1$ and the time step is $\Delta t = 0.00125$. The number of uniformly distributed cells is $n_x = 64$, 128 and 256. Figure 1 shows the profiles of $\phi$, $\rho$ and $u$ for different time instants. The rate of convergence of the algorithm is assessed by calculating the $L_2$-error for different grid refinement levels, where

$$ L_2 = \sqrt{\frac{\sum_{i=1}^{N} (\Phi_{ex} - \Phi_i)^2}{\sum_{i=1}^{N} \Phi_{ex}^2}} \quad (24) $$

and $\Phi$ is some test function. The subscript 'ex' denotes the exact solution while 'i' corresponds to the value in the i-th control volume with the total number N. Figure 2 shows the $L_2$-error
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Log-log plot of $L_2$-error for $\phi$, $\rho$ and $u$ for Problem 1 and $t = 1$. The dashed line shows the slope corresponding to the second-order accurate discretization.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Contours of $\rho(t, x, y)$ for $t = 0$ (top), $t = 0.5$ (middle), $t = 1$ (bottom). Red color corresponds to $\rho = 1$ and blue color to $\rho = 20$.</figcaption></figure>

for $\phi$, $\rho$ and $u$ at the moment $t = 1$. Note that the decrease of the error when $n_x$ is increased confirms the second-order accuracy of the algorithm.

### 3.2. Problem 2: two-dimensional advection and diffusion

The second problem describes the advection of a diffusing two-dimensional ‘flame’ front (see Fig. 3):

$$ \phi(x, y, t) = \frac{1 + \tanh(b\hat{x} \exp(-\omega t))}{(1 + \frac{\rho_0}{\rho_1}) + (1 - \frac{\rho_0}{\rho_1}) \tanh(b\hat{x} \exp(-\omega t))}, \quad (25) $$

$$ \rho(x, y, t) = \left( \frac{\phi(x, y, t)}{\rho_1} + \frac{1 - \phi(x, y, t)}{\rho_2} \right)^{-1}, \quad (26) $$
<figure><img src="image_6.png" /><figcaption>Figure 4. L₂-error for u, v, p, φ and ρ for Problem 2 and t = 1.</figcaption></figure>

$$u(x, y, t) = \frac{\Delta \rho}{\rho} \left( \omega \hat{x} - \frac{\omega \hat{x} - u_F}{\exp(2b\hat{x}) \exp(-\omega t)) + 1} - \frac{\omega \log(\exp(2b\hat{x}) \exp(-\omega t)) + 1}{2b \exp(-\omega t)} \right), \quad (27)$$

$$v(x, y, t) = v_F, \quad p(x, y, t) = 0, \qquad (28)$$

where $\hat{x}(x, y, t) = u_F t - x + a \cos(k(v_F t - y))$, $\Delta \rho = \rho_0 - \rho_1$ and $a, b, k, \omega, u_F$ and $v_F$ are parameters (see Table 1). The computational domain for this problem is $-1 \le x \le 3$, $-1/2 \le y \le 1/2$ and $0 \le t \le 1$. Four meshes are considered with the number of uniform cells in each direction $n_x \times n_y = 200 \times 50, 400 \times 100, 800 \times 200$ and $1600 \times 400$. The time step is $\Delta t = 0.00125$. Figure 4 shows the spatial convergence of the L₂-error with the grid refinement approximately indicating the second-order accuracy of the code.

### 3.3. Problem 3: two-dimensional oscillating density field

Further we consider a time-periodic solution of the following form (see Fig. 5):

$$\phi(x, y, t) = \frac{1 + \sin(\pi k \hat{x}) \sin(\pi k \hat{y}) \cos(\pi \omega t)}{\left(1 + \frac{\rho_0}{\rho_1}\right) + \left(1 - \frac{\rho_0}{\rho_1}\right) \sin(\pi k \hat{x}) \sin(\pi k \hat{y}) \cos(\pi \omega t)}, \quad (29)$$

$$\rho(x, y, t) = \left( \frac{\phi(x, y, t)}{\rho_1} + \frac{1 - \phi(x, y, t)}{\rho_2} \right)^{-1}, \qquad (30)$$

$$u(x, y, t) = \frac{\rho_1 - \rho_0}{\rho(x, y, t)} \left( \frac{-\omega}{4k} \right) \cos(\pi k \hat{x}) \sin(\pi k \hat{y}) \sin(\pi \omega t), \quad (31)$$

$$v(x, y, t) = \frac{\rho_1 - \rho_0}{\rho(x, y, t)} \left( \frac{-\omega}{4k} \right) \sin(\pi k \hat{x}) \cos(\pi k \hat{y}) \sin(\pi \omega t), \quad (32)$$

$$p(x, y, t) = \frac{1}{2} \rho(x, y, t) u(x, y, t) v(x, y, t), \qquad (33)$$

where $\hat{x} = x - u_F t$ and $\hat{y} = y - v_F t$. The computational domain for this problem is $-1 \le x \le 1$, $-1 \le y \le 1$ and $0 \le t \le 1$. Structured hexahedral and unstructured triangular grids have been used for calculations with four levels of refinement, i.e. $n_x \times n_y = 32^2, 64^2, 128^2$ and $256^2$. The time step decreases from $\Delta t = 0.025$ for the coarsest mesh to $\Delta t = 0.003125$ for the finest. Figure 6 shows the L₂-error convergence for hexahedral and triangular grids in comparison with the results of Shunn et al [10]. While unstructured meshes demonstrate similar absolute level of the L₂-error, the hexahedral grids with the present algorithm bring lower error compared to the data from the literature.
<figure><img src="image_7.png" /><figcaption>Figure 5. Colorplots of $\phi$, $\rho$ and $p$ and streamlines at $t = 0$ for Problem 3.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Left: uniform hexahedral mesh. Right: unstructured triangular mesh. Solid lines correspond to the present simulations while dashed lines denote results of Shunn et al [10].</figcaption></figure>

### 3.4. Problem 4: Rayleigh-Taylor instability

Following Desjardins et al [11] we consider the two-dimensional Rayleigh-Taylor instability problem. A rectangular domain of the size $[-0.5, 0.5] \times [-0.5, 0.5]$ filled with two miscible fluids separated by a horizontal perturbed interface. The heavy fluid with $\rho_0 = 1$ is above the light fluid with $\rho_1 = 0.1$. The exact location of the interface is given by

$$y_{int}(x) = -\gamma \sum_{k=1}^{8} \cos(\omega_k \pi x), \quad (34)$$

where $\gamma = 0.001$ is the perturbations amplitude and $\omega_k = 4, 14, 23, 28, 33, 42, 51, 59$. The initial scalar field is set by the relation:

$$\phi(x, y, t = 0) = \frac{1}{2} \left[ 1 + \tanh \left( \frac{y_{int}(x) - y}{2\delta} \right) \right], \quad (35)$$

with $\delta = 0.002$. The density obeys the same equation of state used in previous problems:

$$\rho(x, y, t) = \left( \frac{\phi(x, y, t)}{\rho_0} + \frac{1 - \phi(x, y, t)}{\rho_1} \right)^{-1}. \quad (36)$$

The two fluids have identical kinematic viscosity $\nu = \mu/\rho = 0.001$ and kinematic diffusivity $\alpha_{\phi} = 0.0005$. The gravity body force is added to the momentum equation with the gravity
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Contour plot of $\rho$ at $t = 0.1$ (left), $t = 0.4$ (middle) and $t = 0.75$ (right) computed on a 128 $\times$ 128 grid. The dashed lines correspond to $\rho = 0.2, 0.5$ and 0.8 levels. Red color corresponds to the heavy fluid ($\rho = 1$) and blue color to the light one ($\rho = 0.1$).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** The plot of $\rho$ at $t = 0.75$ and $y = 0.2$ (the line is shown in Fig. 7, right) for different mesh resolutions. The data from the literature [11] is also shown.</figcaption></figure>

acceleration constant set to $g = 9$ so that the Reynolds number is $Re = \sqrt{gL_y}L_x/\nu = 3000$. Simulations have been performed on a uniform hexahedral mesh with $n_x \times n_y = 128^2$ and $512^2$. The time step size is $\Delta t = 0.001$ for the coarse mesh and $\Delta t = 0.00025$ for the fine mesh. The comparison of our results with the data of Desjardins et al [11] shows good agreement although not excellent. In our case coarse and fine simulations robustly capture the same peaks in the profile of $\rho(x, y = 0.2, t = 0.75)$ while in case of Desjardins et al the results are more different with the mesh refinement. We imposed no-slip conditions on top and bottom wall and periodic conditions in $x$-direction for the side walls while Desjardins et al did not mention the boundary conditions. This issue could be a possible source of deviations.

## 4. Turbulent test case

In this section we introduce a recent approach to the Large-eddy simulation of premixed turbulent combustion called the ADEF model (Approximate Deconvolution and Explicit flame Filtering) for subgrid-scale modeling of scalar fields. First, let us define the spatial filtering
procedure denoted by an overbar symbol of some test function $\phi$:

$$ \bar{\phi}(x, t) = \int_{-\infty}^{+\infty} \phi(x', t)G(x - x')dx', \qquad (37) $$

where G is the filter kernel. The density weighted filtering (similar to Favre time-averaging) is introduced as follows

$$ \bar{\rho}\tilde{\phi}(x, t) = \bar{\rho}\phi(x, t) = \int_{-\infty}^{+\infty} \rho(x', t)\phi(x', t)G(x - x')dx', \qquad (38) $$

Applying the filtering operation to the system of Eqs. (1)-(4) we obtain

$$ \frac{\partial \bar{\rho}}{\partial t} + \frac{\partial \bar{\rho} \tilde{u}_j}{\partial x_j} = 0, \qquad (39) $$

$$ \frac{\partial \bar{\rho} \tilde{u}_i}{\partial t} + \frac{\partial \bar{\rho} \tilde{u}_i \tilde{u}_j}{\partial x_j} = -\frac{\partial \bar{\rho}}{\partial x_i} + \frac{\partial}{\partial x_j} \left( 2\mu(\tilde{\phi})\tilde{S}_{ij} \right) - \frac{\partial}{\partial x_j} \left( \tau_{ij}^{c,u} - \tau_{ij}^{d,u} \right), \qquad (40) $$

$$ \frac{\partial \bar{\rho} \tilde{\phi}}{\partial t} + \frac{\partial \bar{\rho} \tilde{\phi} \tilde{u}_j}{\partial x_j} = \dot{\omega}(\tilde{\phi}) + \frac{\partial}{\partial x_j} \left( \bar{\rho} \alpha(\tilde{\phi}) \frac{\partial \tilde{\phi}}{\partial x_j} \right) - \frac{\partial}{\partial x_j} \left( \tau_j^{c,\phi} - \tau_j^{d,\phi} \right), \qquad (41) $$

$$ \bar{\rho} = f(\tilde{\phi}), \qquad (42) $$

where

$$ \tau_{ij}^{c,u} = \overline{\rho u_i u_j} - \bar{\rho} \tilde{u}_i \tilde{u}_j, \qquad \tau_{ij}^{d,u} = \overline{2\mu(\tilde{\phi})S_{ij}} - 2\mu(\tilde{\phi})\tilde{S}_{ij}, \qquad (43) $$

$$ \tau_j^{c,\phi} = \overline{\rho\phi u_j} - \bar{\rho}\tilde{\phi}\tilde{u}_j, \qquad \tau_j^{d,\phi} = \overline{\rho\alpha(\tilde{\phi})\frac{\partial\tilde{\phi}}{\partial x_j}} - \bar{\rho}\alpha(\tilde{\phi})\frac{\partial\tilde{\phi}}{\partial x_j}, \qquad (44) $$

The subgrid-scale terms, Eq. (43), in the momentum equation are expressed with the dynamic Smagorinsky model:

$$ \tau_{ij}^{c,u} - \tau_{ij}^{d,u} = -2\mu_t \tilde{S}_{ij}, \qquad (45) $$

where $\mu_t$ is the time and space-dependent parameter which is calculated using a standard dynamic routine [16]. The unclosed $\tau_j^{c,\phi}$ term in Eq. (41) can also be modeled using the Boussinesq hypothesis:

$$ \tau_j^{c,\phi} = -\bar{\rho}\alpha_t \frac{\partial \tilde{\phi}}{\partial x_j}, \qquad (46) $$

At the same time according to the thickened flame approach [17]:

$$ \overline{\rho\alpha(\tilde{\phi})\frac{\partial\tilde{\phi}}{\partial x_j}} = \mathcal{F}\bar{\rho}\alpha(\tilde{\phi})\frac{\partial\tilde{\phi}}{\partial x_j}, \qquad (47) $$

where

$$ \mathcal{F} = \frac{\dot{\omega}(\tilde{\phi})}{\dot{\omega}(\tilde{\phi})} \qquad (48) $$
In the end, the equation for the filtered scalar field is the following:

$$ \frac{\partial \bar{\rho} \tilde{\phi}}{\partial t} + \frac{\partial \bar{\rho} \tilde{\phi} \tilde{u}_j}{\partial x_j} = \overline{\dot{\omega}(\phi)} + \frac{\partial}{\partial x_j} \left[ \bar{\rho} (\mathcal{F}\alpha + \alpha_t) \frac{\partial \tilde{\phi}}{\partial x_j} \right]. \quad (49) $$

Momentum and scalar equations contain unclosed terms which can be expressed if the inverse filtering operation is defined [12], for example, the chemical source term is

$$ \overline{\dot{\omega}(\phi)} = \overline{\dot{\omega}(\tilde{\mathcal{L}}^{-1}[\tilde{\phi}])}, \quad (50) $$

where $\tilde{\mathcal{L}}^{-1}$ denotes the inverse of the filtering operation. For a Gaussian filter the inverse operator can be approximately described as (see [18])

$$ \phi(x, t) = \mathcal{L}^{-1}[\bar{\phi}(x, t)] = \bar{\phi}(x, t) - \frac{\Delta^2}{24} \frac{\partial^2 \bar{\phi}}{\partial x_k^2}, \quad (51) $$

where $\Delta$ is the local filter width. The above described model further will be applied to the laboratory stratified burner [13] for the validation purposes.

## 5. Conclusions

In this work we described the numerical algorithm and assessed the accuracy of the low Mach variable-density Navier-Stokes solver on a number of analytic manufactured solutions and Rayleigh-Taylor instability problem. It was shown that the code reproduces second order of accuracy expected from the program. After the verification step we will validate the recent Approximate Deconvolution and Explicit flame Filtering combustion model on experimental data from a number of laboratory burners.

## Acknowledgments

This work is funded by RFBR grants No. 15-08-06166 and No. 14-08-01233a. The computational resources are provided by Siberian Supercomputer Center SB RAS (Novosibirsk) and Joint Supercomputer Center RAS (Moscow).

## References

[1] Oberkampf WL, Roy CJ 2010 Verification and Validation in Scientific Computing (Cambridge: Cambridge University Press)

[2] Shunn L, Ham F 2007 Method of manufactured solutions applied to variable-density flow solvers Annu. Res. Briefs, CTR 155-168

[3] Knikker R 2011 A comparative study of high-order variable-property segregated algorithms for unsteady low Mach number flows Int. J. Numer. Meth. Fluids 66 403-427

[4] Vedovoto JM, da Silveira Neto A, Mura A, da Silva LF 2011 Application of the method of manufactured solutions to the verification of a pressure-based finite-volume numerical scheme Comp. Fluids 51 (1) 85-99

[5] Ferrer PJ, Buttay R, Lehnasch G, Mura A 2014 A detailed verification procedure for compressible reactive multicomponent Navier-Stokes solvers Comp. Fluids 89 88-110

[6] Waltz J, Canfield TR, Morgan NR, Risinger LD, Wohlbier JG 2014 Manufactured solutions for the three-dimensional Euler equations with relevance to Inertial Confinement Fusion J. Comp Phys. 267 196-209

[7] Tyliszczak A 2014 Projection method for high-order compact schemes for low Mach number flows in enclosures Int. J. Numer. Meth. Fluids 24 (5) 1141-1174

[8] Tyliszczak A 2016 High-order compact difference algorithm on half-staggered meshes for low Mach number flows Comp. Fluids 127 131-145

[9] Choudhary A, Roy CJ, Dietiker JF, Shahnam M, Garg R, Musser J 2016 Code verification for multiphase flows using the method of manufactured solutions Comp. Fluids 80 150-163

[10] Shunn L, Ham F, Moin P 2012 Verification of variable-density flow solvers using manufactured solutions J. Comp. Phys. 231 (9) 3801-3827
[11] Desjardins O, Blanquart G, Balarac G, Pitsch H 2008 High order conservative finite difference scheme for variable density low Mach number turbulent flows *J. Comp. Phys.* **227** (15) 7125–7159

[12] Domingo P, Vervisch L 2015 Large Eddy Simulation of premixed turbulent combustion using approximate deconvolution and explicit flame filtering *Proc. Comb. Inst.* **35** (2) 1349–1357

[13] Zhou R, Balusamy S, Sweeney MS, Barlow RS, Hochgreb S 2013 Flow field measurements of a series of turbulent premixed and stratified methane/air flames *Comb. Flame* **160** (10) 2017–2028

[14] Pierce CD 2001 *Progress-variable approach for large-eddy simulation of turbulent combustion* (PhD, Stanford University)

[15] Ničeno B, Hanjalić K. 2005 Unstructured large-eddy-and conjugate heat transfer simulations of wall-bounded flows in *Modeling and Simulation of Turbulent Heat Transfer (Developments in Heat Transfer)* (UK: WIT)

[16] Moin P, Squires K, Cabot W, Lee S 1991 A dynamic subgridscale model for compressible turbulence and scalar transport *Phys. Fluids A* **3** (11) 2746–57

[17] Colin O, Ducros F, Veynante D, Poinsot T 2000 A thickened flame model for large eddy simulations of turbulent premixed combustion *Phys. Fluids* **12** (7) 1843-63

[18] Katopodes FV, Street RL, Ferziger JH 2000 A theory for the subfilter-scale model in large-eddy simulation *Tech. Report, Stanford K1*