A minimally-dissipative low-Mach number solver for complex reacting flows in OpenFOAM

TL;DR

This paper introduces an open-source LES solver for low-Mach number turbulent combustion in OpenFOAM, emphasizing minimal kinetic energy dissipation for improved accuracy in complex geometries.

Contribution

It develops a minimally-dissipative, collocated-mesh LES solver tailored for unstructured grids in OpenFOAM, enhancing simulation fidelity for complex reacting flows.

Findings

Validated on manufactured solutions for accuracy.

Demonstrated on canonical and turbulent sooting flame cases.

Achieved minimal kinetic energy dissipation in simulations.

Abstract

Large eddy simulation (LES) has become the de-facto computational tool for modeling complex reacting flows, especially in gas turbine applications. However, readily usable general-purpose LES codes for complex geometries are typically academic or proprietary/commercial in nature. The objective of this work is to develop and disseminate an open source LES tool for low-Mach number turbulent combustion using the OpenFOAM framework. In particular, a collocated-mesh approach suited for unstructured grid formulation is provided. Unlike other fluid dynamics models, LES accuracy is intricately linked to so-called primary and secondary conservation properties of the numerical discretization schemes. This implies that although the solver only evolves equations for mass, momentum, and energy, the implied discrete equation for kinetic energy (square of velocity) should be minimally-dissipative. Here, a specific spatial and temporal discretization is imposed such that this kinetic energy dissipation is minimized. The method is demonstrated using manufactured solutions approach on regular and skewed meshes, a canonical flow problem, and a turbulent sooting flame in a complex domain relevant to gas turbines applications.

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{u})=0,$$

$$\frac{\partial \rho \mathbf{u}}{\partial t}+\nabla \cdot (\rho \mathbf{u}\mathbf{u})=-\nabla p+\nabla \cdot \overline{\mathbf{\sigma }},$$

$$\frac{\partial \rho \varphi }{\partial t}+\nabla \cdot (\rho \varphi \mathbf{u})=\nabla \cdot (D\nabla \varphi )+\dot{\omega },$$

$$\frac{\rho {\mathbf{u}}^{\ast }-\rho {\mathbf{u}}^n}{\Delta t}+\mathcal{C}=\mathcal{D},$$

$$\rho^{\prime}\approx\Big{(}\frac{\partial\rho}{\partial p}\Big{)}_{T}p^{\prime},$$

$$V_p\frac{{\rho }^{n+1}-{\rho }^n}{\Delta t}+\sum _f{\mathbf{S}}_{\mathbf{f}}\cdot (\rho {\mathbf{u}}_f^{\ast }+{\mathbf{\varphi }}_{\mathbf{f}}^{\prime })=0,$$

$${\mathbf{\varphi }}_{\mathbf{f}}^{\prime }=({\rho }^{n+1}{\mathbf{u}}^{\prime }{)}_f+({\mathbf{u}}_m^{\ast }{\rho }^{\prime }{)}_f,$$

$${\nabla }^2{p}_m=\frac{1}{\Delta t}(\frac{\frac{{\rho }^{n+1}}{p_{m-1}}{p}_m-\frac{{\rho }^n}{p^n}{p}^n}{\Delta t}+\nabla \cdot (\rho {\mathbf{u}}^{\ast })).$$

$${\nabla }^{2}p=\frac{1}{\Delta t}(\frac{{\rho }_m-{\rho }^{n-1}}{\Delta t}+\nabla \cdot (\rho {\mathbf{u}}^{\ast })).$$

$${\rho }_f{{\mathbf{u}}_{\mathbf{f}}}^{n+1}\cdot {\mathbf{S}}_{\mathbf{f}}=2{\rho }_f{{\mathbf{u}}_{\mathbf{f}}}^n\cdot {\mathbf{S}}_{\mathbf{f}}-{\rho }_f{{\mathbf{u}}_{\mathbf{f}}}^{n-1}\cdot {\mathbf{S}}_{\mathbf{f}}.$$

$$\sum _{f,in}{\rho }_f{\mathbf{u}}_{\mathbf{f}}\cdot {\mathbf{S}}_{\mathbf{f}}-\sum _{f,out}{\rho }_f{\mathbf{u}}_{\mathbf{f}}\cdot {\mathbf{S}}_{\mathbf{f}}+\sum _{cell}\frac{\delta \rho }{\delta t}\delta V_{cell}=0,$$

$$mas{s}_{err}(x,t)=V_p\frac{{\rho }^{n+1}-{\rho }^n}{\Delta t}+\sum _f{\rho }_f{\mathbf{u}}_{\mathbf{f}}(x,t)\cdot {\mathbf{S}}_{\mathbf{f}}.$$

$$\frac{\partial Z_{mix}}{\partial t}+\nabla \cdot (\mathbf{u}{Z}_{mix})=\nabla \cdot (D\nabla Z_{mix}).$$

$$\frac{\partial \rho u^2/2}{\partial t}+\nabla \cdot (\rho \mathbf{u}{u}^2/2)=-\nabla \cdot (p\mathbf{u})+p\nabla \cdot \mathbf{u}+\mathbf{u}\cdot (\nabla \cdot \overline{\sigma }).$$

$$\mathbf{u}\cdot \{\frac{\partial \rho \mathbf{u}}{\partial t}\}+\mathbf{u}\cdot \{\nabla \cdot (\rho \mathbf{u}\mathbf{u}))\}=\frac{\partial \rho u^2/2}{\partial t}+\nabla \cdot (\rho \mathbf{u}{u}^2/2).$$

$${\overline{\sqrt{{\overline{\rho }}^t}}}^t\frac{{\delta }_1\sqrt{{\overline{\rho }}^t}\mathbf{u}}{{\delta }_{1}t}+\frac{1}{2}(\frac{{\delta }_1{\overline{{\mathbf{\varphi }}_f}}^t{\overline{\hat{\mathbf{u}}}}^x}{{\delta }_{1}x}+{\overline{{\overline{{\mathbf{\varphi }}_f}}^t\frac{{\delta }_1\hat{\mathbf{u}}}{{\delta }_{1}x}}}^x),$$

$$\hat{\mathbf{u}}=\frac{{\overline{\sqrt{{\overline{\rho }}^t}u}}^t}{{\overline{\sqrt{{\overline{\rho }}^t}}}^t},$$

$$\sum _f{\overline{{\mathbf{\varphi }}_{\mathbf{f}}}}^t\cdot {\mathbf{S}}_{\mathbf{f}}+\frac{{\delta }_1{\overline{\rho }}^t}{{\delta }_{1}t}=0,$$

$$\frac{1}{2}(\sqrt{{\overline{\rho }}^t}\mathbf{u}{)}^2.$$

$${\dot{\omega }}^{n+1}=\frac{\dot{\omega }({\varphi }^{n+3/2})+\dot{\omega }({\varphi }^{n+1/2})}{2}.$$

$$\frac{{\rho }^{n+3/2}-{\rho }^{n+1/2}}{\Delta t}+\nabla \cdot ({{\mathbf{\varphi }}_{\mathbf{f}}}^{\ast })=\nabla \cdot (a^{-1}.{\overline{\rho }}^t\nabla p^{n+1/2}),$$

$$u_x(x,y)=-cos(\pi x)sin(\pi y),$$

$$u_y(x,y)=sin(\pi x)cos(\pi y),$$

$$p(x,y)=-\frac{1}{4}(cos(2\pi x)+cos(2\pi y)).$$

$$x^{\prime }=x+0.2sin(\pi y),$$

$$y^{\prime }=y+0.2sin(\pi x).$$

$$u_x(x,y,t)=\frac{{\rho }_1-{\rho }_0}{\rho (x,y,t)}(\frac{-\omega }{4k})cos(\pi kx)sin(\pi ky)sin(\pi \omega t),$$

$$u_y(x,y,t)=\frac{{\rho }_1-{\rho }_0}{\rho (x,y,t)}(\frac{-\omega }{4k})sin(\pi kx)cos(\pi ky)sin(\pi \omega t).$$

$$\frac{{\rho }^{n+1}{\mathbf{u}}^g-{\rho }^n{\mathbf{u}}^n}{\Delta t}+{\mathcal{C}}^g=-\nabla p^n.$$

$$\frac{{\rho }^{n+1}{\mathbf{u}}^{\ast }-{\rho }^n{\mathbf{u}}^n}{\Delta t}+{\mathcal{C}}^g=0.$$

Full text

Large eddy simulation (LES) has become the de-facto computational tool for modeling complex reacting flows, especially in gas turbine applications. However, readily usable general-purpose LES codes for complex geometries are typically academic or proprietary/commercial in nature. The objective of this work is to develop and disseminate an open source LES tool for low-Mach number turbulent combustion using the OpenFOAM framework. In particular, a collocated-mesh approach suited for unstructured grid formulation is provided. Unlike other fluid dynamics models, LES accuracy is intricately linked to so-called primary and secondary conservation properties of the numerical discretization schemes. This implies that although the solver only evolves equations for mass, momentum, and energy, the implied discrete equation for kinetic energy (square of velocity) should be minimally-dissipative. Here, a specific spatial and temporal discretization is imposed such that this kinetic energy dissipation is minimized. The method is demonstrated using manufactured solutions approach on regular and skewed meshes, a canonical flow problem, and a turbulent sooting flame in a complex domain relevant to gas turbines applications.

1 Introduction

The use of large eddy simulation in complex turbulent flows has increased substantially in the last decade. This advance is being driven by the rapid growth in computational power, as well as advances in numerical algorithms for such complex flows. This has led to the development of a number of LES solvers that are routinely deployed in industrial applications [1, 2, 3, 4, 5]. Nevertheless, many of these solvers remain proprietary and not open to the research community. This deficiency is particularly challenging for LES due to the sensitivity of LES results to numerical accuracy. For instance, when head-to-head comparison of LES models are made, such studies rely on different LES solvers and numerical methods. Consequently, the conclusions are highly sensitive to these numerical details, and may not even be relevant to model comparison [6]. Furthermore, many combustion models and subfilter closures are developed for canonical flow problems and are not exercised in full-scale geometries or application-relevant flow problems. This sparsity in full-scale validation is at least partially due to the lack of an easy approach to porting models to complex flow problems. As LES models mature, there is a clear need for a robust open source platform to demonstrate their performance for practical configurations. In this context, OpenFOAM [4] is a prime candidate for such a framework. Based on a field operation approach [7], OpenFOAM provides a convenient code-base for numerically solving partial differential equations. Moreover, OpenFOAM has developed a broad community of developers and users, who have added valuable tools and methods to the base solver [8]. This robust development ecosystem has been leveraged extensively in the turbulent flow and turbulent combustion research and application communities [9, 10, 11, 12]. The current work intends to upgrade the implementation of LES algorithms that are indispensable in solving variable density reacting flows of interest to gas turbine applications. With the will of gradually moving combustion research applications to OpenFOAM, key issues regarding LES numerical accuracy are addressed.

Many gas turbine applications, especially related to the combustor section, operate in the low-Mach number regime, characterized by velocities smaller than 0.3 Ma, where the Mach number is defined based on local fluid properties. In this regime, the acoustic component of the Navier-Stokes equations is decoupled from the basic flow physics. Consequently, by reformulating the governing equations for this low-Ma regime, it is possible to go beyond the CFL restriction imposed by the acoustic wave speed (i.e., local sound speed). In many applications, this can amount to an order of magnitude increase in timestep used. Although OpenFOAM has been used in combustion applications, many of these are related to the compressible flow regime (see, for instance [10]). In the base distribution of OpenFOAM, the variable density solver the closest to a low-Mach number solver uses an all-Mach approach, where some level of compressibility-related coupling of the governing equations is retained even in the limit of zero Mach number [13]. Hence, developing a robust low-Ma solver for variable density flows will be of interest to combustion applications.

Another issue of importance to LES is the coupling between numerical discretization and modeling errors, as discussed in [14, 15, 16, 17, 18, 19]. Briefly, there are two types of numerical errors. First, the spatial discretization of derivative operators for fields containing high wave number components can be highly erroneous [20]. In practical LES, features comparable to the filter size can thus be contaminated by discretization errors [21, 14, 18]. In turbulent combustion applications, where small-scale models are critical for capturing the mixing and reactions processes, this discretization error represents a major concern. The most comprehensive solution is to use an explicit filtering technique, where the small-scales are removed through a filter during the simulation in order to prevent deposition of energy at such scales [22, 14, 15, 23]. However, extending these techniques to variable density flows has not been achieved yet [24]. Several approaches [19, 25] have been proposed to mitigate these numerical effects. The second issue concerns the so-called conservation of secondary quantities. Finite-volume approaches discretely conserve primary quantities such as mass and momentum. With LES, the accurate representation of the turbulent energy spectrum is important the validity of the modeling assumptions. Therefore, minimizing numerical dissipation of kinetic energy is a key point of the LES solvers. In low-Mach number incompressible flows, and in the absence of viscous dissipation, kinetic energy should be exactly conserved. In this case, the spatial and temporal discretization should ensure that such conservation or, at the least, a minimization of dissipation is achieved [13]. In this regard, there has been considerable progress in the design of energy-conserving numerical schemes [26, 27, 28, 29, 30]. Again, a comprehensive implementation of such tools in the OpenFOAM framework is not available.

With this introduction, the objectives of this work are as follows: 1) Analyze the variable density solvers in OpenFOAM, and implement a consistent low-Ma solver that preserves spatial and temporal accuracy in the limit of zero Mach number, and 2) evaluate the energy conserving properties of existing OpenFOAM solvers, and implement a minimally-dissipative approach. In Sec. 2, a variable density low-Mach number solver is designed and its implementation in OpenFOAM is described. In particular, it is stressed that the baseline implementation of variable density flow solver in OpenFOAM was previously only allowed for compressible flow cases. The solver is then tested using the method of manufactured solutions. Sec. 3 provides the necessary theoretical background for the design of energy conservative solvers. In Sec. 4, the implementation of a minimally-dissipative solver is described for the OpenFOAM framework. Finally in Sec. 5, a set of verification and validation cases are used to demonstrate the capabilities of the new solver. Sec. 3, energy conservation properties of OpenFOAM discretization operators are studied, followed by the implementation of a fully-conservative scheme.

2 Low-Mach number solvers in OpenFOAM

The governing equations of fluid flow of interest here are written as

[TABLE]

[TABLE]

where $\rho$ is the flow density, $\boldsymbol{u}$ is the local gas phase velocity vector, $p$ is the mechanical pressure and $\overline{\boldsymbol{\sigma}}$ is the viscous stress tensor.

When considering chemical reactions, additional equations that describe transport of chemical species and a formulation that couples any heat addition to the density term need to be included. A generic scalar transport equation that may be used for this purpose is given by:

[TABLE]

where $\phi$ denotes the transported scalar and $D$ denotes mass diffusivity and $\dot{\omega}$ is some volumetric source term.

To solve this system of equations, a numerical approach is used. Three different types of flow solvers are defined here for the sake of future discussion. An incompressible solver is defined as a solver that does not take into account any density change. In other words, the density field is treated as a constant value. A compressible solver is defined as a solver that takes into account the dynamic coupling between the pressure and density fields. In other words, this solver allows mechanical energy (kinetic energy noted KE) to be converted to thermal energy (through pressure). A low-Mach number solver is defined as a solver that does not couple density changes with instantaneous pressure changes. In particular, the pressure field is split into a thermodynamic pressure and a mechanical pressure, with the former held constant while the latter is allowed to vary through velocity changes in the flow. The thermodynamic pressure is used in the equation of state, while the mechanical pressure appears in the momentum transport equation (Eq. 2). Density variations in low-Mach number solvers occur through heat addition or removal, for instance due to chemical reactions. Such low-Ma number flows and associated solvers are the focus of this work.

2.1 Hybrid solvers in OpenFOAM for variable density flows

The base distribution of OpenFOAM [4] contains a suite of solvers targeting reacting flows that involve density changes. Although they nominally fall under the compressible solver definition provided above, some aspects of the low-Mach representation is also included. The main choice is as follows: all variable density solvers use the equation of state to reduce the number of partial differential equations needed (for the thermochemical part) to two of the three variables (pressure, temperature, density). Assuming that an energy equation that provides temperature is available, the choice is then to solve for density and obtain pressure from equation of state or vice-versa. These two types of compressible solvers are: the density-based solvers where the density is transported using the continuity equation, and the pressure-based solvers where the thermodynamical pressure is computed using a pressure correction method. In the pressure-based solvers handling variable density available in OpenFOAM, a generic pressure correction procedure [13, Ch. 10.2] is used. This approach is easily applicable even with unstructured and complex grids. The pressure correction equation contains an incompressible divergence term (correcting the mass fluxes) and a compressible convective term (correcting the density). Each one of these term becomes dominant when the flow is largely incompressible or compressible, respectively. This pressure-correction strategy is referred to as a hybrid approach in the remainder of the paper. This type of pressure correction has been used to simulate combustion cases at high and low-Mach numbers [31, 9].

Since the pressure correction has an impact on the density field, this procedure involving momentum, scalar transport and pressure correction can be repeated in order to fully couple all transported variables. These iterations are called outer-iterations as opposed to iterations which denote the time advancement of the variables. Within each outer-iteration (denoted by $(.)_{m}$), the momentum field is first advanced from the old timestep (denoted by $(.)^{n}$ below) to an intermediate fractional time (denoted by $(.)^{*}$ below) without taking into account the pressure gradient.

[TABLE]

where $\boldsymbol{\mathcal{C}}$ denotes the convective term, $\boldsymbol{\mathcal{D}}$ the diffusion term and $\Delta t$ is the flow timestep.

The time at which $\boldsymbol{\mathcal{C}}$ and $\boldsymbol{\mathcal{D}}$ is evaluated in Eq. 4 is not specified since the following discussion is about completing the time derivative term with a pressure correction. The field $\boldsymbol{u}^{*}$ does not necessarily respect the integral continuity equation. It is therefore corrected with a pressure gradient which is tailored to enforce mass conservation in the domain. Because this correction is done through pressure, the thermodynamical state also depends on this correction procedure. This pressure correction equation is therefore solved to determine the correction to apply to the mass fluxes as obtained from the intermediate field $\boldsymbol{u}^{*}$.

In pressure-based compressible solvers, the pressure correction equation is formulated by replacing density changes by pressure changes under the assumption that temperature is held constant:

[TABLE]

where $\rho^{\prime}$ denotes a density change and $p^{\prime}$ denotes a pressure change. In the notations adopted in OpenFOAM 4.1, the multiplicative factor on the right-hand-side (RHS) of Eq. 5 is called psi. The derivation of this method is covered in detail in [13, Ch. 10.2]. Similar to all predictor corrector methods, the pressure field is used as a corrector which adjusts the mass fluxes in order to ensure conservation of mass at the new timestep (denoted by $(.)^{n+1}$). Unlike a low-Mach implementation where the pressure gradient only corrects the flow velocity, the pressure correction also acts on the density in pressure-based compressible solvers. The goal of the procedure is to find the mass flux correction $\boldsymbol{\phi^{\prime}_{f}}$ such that

[TABLE]

leading to

[TABLE]

where $V_{p}$ is the cell volume, $(.)_{f}$ denotes a variable defined at a cell face, $\boldsymbol{S_{f}}$ denotes the cell face normal vector which norm is the cell face surface area. The first term on the RHS of Eq. 7 comes from the additional pressure gradient due to the computed pressure correction. The second term of the RHS comes from the correction of density which goes to zero in a low-Mach number case [13, Ch. 10.2]. Whenever the hybrid solver is compared to a low-Mach number solver, only the first RHS term is therefore included. Injecting Eq. 7 into Eq. 6, and using the compressibility factor to transform the derivation of density into pressure, the pressure equation can be formulated for the hybrid solver. Only the time derivative term is discretized in Eq. 8 since it is where the numerical differences with a strict low-Mach number solver become apparent.

[TABLE]

The above equation requires a solution to the elliptic Poisson equation system at each outer-iteration step, which is typically computationally expensive. Hence, the cost of the simulation is tied to the number of outer-iterations used in each timestep.

2.2 Strict low-Mach number solvers for variable density flows

The hybrid approach used in the pressure-based solver of OpenFOAM offers flexibility and allows to deal with flows spanning a large range of Mach numbers. This method is based on the assumption that the compressible and the incompressible terms of the pressure correction equation will be effectively negligible in the right context. However, even in the limit of a low-Mach number flow, some compressibility is retained, and the solver cannot be considered as a strictly acoustically decoupled formulation.

In the context of direct and large eddy simulations, several low-Mach formulations have been developed in the past [32, 33, 34, 35], which are generally semi-implicit in the temporal direction. The objective of this work is to develop a similar approach for implementation in OpenFOAM. In these methods, the solver is segregated in that the different partial differential equations are solved sequentially within each outer-iteration, and multiple outer-iterations are used to couple the scalar and the velocity fields [32]. Within each outer-iteration, the scalar fields are first advanced. Second, the impact of the scalar fields on the flow variables (density or viscosity for example) are then computed and these variables are updated. Third, the fractional timestep method is used to first advance the momentum equation without the pressure gradient. For an implicit Euler method, the momentum equation discretized for time only takes the form shown in Eq. 4.

Similar to the pressure based solver presented in Sec. 2.1, the momentum equation is updated using the pressure gradient obtained from a pressure-correction equation. As opposed to the hybrid solver, which links the pressure to the thermodynamical variables in the pressure equation, the sole role of pressure in a low-Mach number solver is to enforce mass conservation. Hence, this pressure variable is referred to as the mechanical pressure, which is different from the thermodynamical pressure that determined the thermophysical quantities such as density and temperature. For a low-Mach number solver, the pressure equation is derived from the continuity equation and takes the form

[TABLE]

Finally, the velocity field is corrected using the newly computed pressure gradient. At the start of the first outer-iteration, the most updated fields available still correspond to the previous timestep. Any term that is required to be defined at the next timestep is estimated using an explicit extrapolation in time. Typically, the mass flux at the face needs to be obtained in this way for the convective terms.

[TABLE]

A similar requirement holds for density in variable density cases in order to accurately compute the time derivative terms.

2.3 Implementation of a low-Mach number solver for combustion applications

In this section, the details pertaining to the combustion model that is necessary for implementing the low-Mach number solver in OpenFOAM are discussed.

2.3.1 Flow-combustion coupling through density field

Here, it is important to consider the nature of combustion models for low-Mach number flows. While different combustion models are available [36, 37, 38, 39, 40], the coupling procedure between reacting scalars and the density/momentum equation set is approximately similar. In general, a set of scalars is used to obtain the gas-phase density, viscosity and individual species diffusivity. When the scalars are species mass fractions, the local thermodynamical properties of the mixture can be directly extracted. When the scalars are not species mass fractions, a mapping between these scalars and the thermodynamical properties is first established and used throughout the domain. In some implementations [41, 42], an energy equation is solved, for which the source term due to chemical reactions is provided based on the scalar fields. The density field is then obtained from the energy field, thermodynamic pressure, and scalar fields through a non-linear inversion process.

For variable density unsteady flows (called compressible in OpenFOAM 4.1), the existing turbulence models in OpenFOAM take as an input a correlation between density and pressure. In order to decouple the variations in thermophysical quantities (density, diffusivity, viscosity) from pressure fluctuations, it is necessary to formulate a turbulence model which does not take the thermodynamical model as an input. Instead it takes the thermophysical variables obtained from the combustion model. For this purpose, a new turbulence model compSansThermo has been developed. The full source code can be obtained as a supplementary material of this paper.

2.3.2 Pressure solver

As mentioned in Sec. 2.1, the pressure correction method used in pressure-based solvers is written to generically handle compressible or incompressible flows. Although similar to Eq. 9, the hybrid pressure equation given by Eq. 2.1 will always result in a different correction for the mass fluxes from the correction that would be provided by low-Mach number solvers. In the hybrid formulation available in the base OpenFOAM distribution, the pressure equation reduces to the low-Mach number pressure equation if and only if the pressure does not change between two subsequent outer-iterations. In practice this is never exactly satisfied, and this results in numerical compressibility effects, even for zero Mach number flows (Sec. 2.4).

The pressure solver ensures discrete mass conservation by re-arranging velocities throughout the domain. A pre-requisite for this solver to function is that the initial discretely non-conserving field satisfies global mass conservation. In other words, the integral of the continuity equation over the computational domain should be satisfied. Here, this is achieved by first altering the outflow surface fluxes to match the sum of density changes in each control volume as well as the fluxes across the boundary surfaces. The global continuity equation can then be written as

[TABLE]

where $\frac{\delta}{\delta t}$ denotes some time discretization that is not of interest for now.

In order to use the low-Mach number solver in OpenFOAM, the boundary conditions of the pressure fields need to be adequately chosen. Neumann boundary conditions should be applied where velocity must not be corrected, typically at the inlet and at the walls. At the outlet of the domain, the user has the choice between Neumann and Dirichlet boundary conditions. With a Neumann boundary conditions, the outlet velocity adjusted to ensure global mass conservation through the function adjustPhi in OpenFOAM 4.1 is preserved. In the original function adjustPhi of OpenFOAM 4.1, the last term of the left-hand side of Eq. 11 is missing which is expected since incompressible solvers were not designed to handle variations of density. With a Dirichlet boundary condition, the outlet velocity is not preserved anymore and the global mass conservation is only enforced up to the addition of the mass conservation errors in each cell of the domain (this is the telescopic property of mass conservation). Although the Dirichlet boundary conditions can lead to non-exact global mass conservation, it was observed to speed up the calculations, especially in open geometries.

2.4 Low-Mach number solver verification and comparison with baseline hybrid solver

The method of manufactured solutions (MMS) is used to verify the solvers and analyze their performance. The one-dimensional variable density flow problem of Shunn et al. [43] is studied. Since many of the changes implemented for the low-Mach number solver are done for the mass fluxes and the pressure correction, this test case clearly highlights the difference with the hybrid procedure. The 1D case is computed with the following parameters (see [43] for details of the flow): $\rho_{0}=10;\rho_{1}=1;\mu=0.03;k_{1}=4;k_{2}=2;w_{0}=50$. The density field is obtained analytically from the value of the transported scalar field $\phi$. Note that no equation of state link the pressure and the density. The analytical solution is shown in Fig. 1.

The grid sizes are varied from $64$ to $512$ cells for a 2 m domain. The timestep is held constant to $6.25$ $\mu$s. These cases are used to compare how the original mixed compressible and incompressible approach (referred to as hybrid in the following) performs in OpenFOAM compares to the strict low-Mach approach. The results shown in Fig. 2 as the L2-error to the analytical solution integrated over the domain and plotted against time.

The hybrid pressure correction consistently leads to an increased L2-error for the considered test case, often exceeding by at least one order magnitude. It is likely that the error comes from the approximation of the time derivative of density by the left hand side (LHS) of Eq. 8.

Finally, the impact of the hybrid formulation on mass conservation is investigated. In Fig. 3, the mass balance at every cell is plotted for the hybrid and the low-Mach number solver. The mass balance is computed using mass fluxes obtained from the velocity interpolations at each face (Eq. 12) as follows:

[TABLE]

While the low-Mach number solver ensures mass conservation up to machine precision, the hybrid solver generates an error proportional to $\frac{\rho^{n+1}-\rho^{n}}{\Delta t}-\frac{\frac{p_{m}}{p_{m-1}}\rho^{n+1}-\rho^{n}}{\Delta t}$. This difference numerically comes from the difference between the pressure at two successive outer-iterations. This is particularly prominent at locations where the density varies the most, which in this problem is near the left boundary of the 1D domain.

3 Non-dissipative solvers

3.1 Motivation for reducing numerical dissipation in OpenFOAM

When simulating reacting flows using large eddy simulation (LES), there has to be a special focus on numerical errors due to the strong coupling between discretization schemes and the notion of filtering in this technique [16, 18, 19, 44]. Unlike in Reynolds-averaged Navier Stokes (RANS) formulations, the simulated field is not smooth with respect to the mesh spacing. More specifically, numerical dispersion errors arising from Taylor series based finite-difference/finite-volume methods introduce scale-dependent errors, with the highest numerical errors occurring at scales close to the filter size [18]. In order to assess and maintain accuracy of solvers, the concept of secondary conservation has been considered. Here, errors in quantities that are not directly solved for, but obtained from output fields, are minimized by appropriate choice of numerical schemes. In the context of turbulent flows, one such key quantity is kinetic energy (KE), which dictates the level of resolved turbulence captured by the numerical approach.

To illustrate the numerical dissipation problem, an constant density bluff-body flow simulation is carried out using two codes: the pisoFoam solver used with linear interpolations in OpenFOAM, and the NGA solver used with second order schemes [33]. The pisoFoam solver does not include variable density i.e. it is written using an incompressible pressure correction equation and is therefore a low-Mach number solver. Thus, the differences observed between both codes are not due to the numerical compressibility effect noted in Sec. 2. The flow configuration is illustrated in Fig. 4. When used for reacting flows, fuel is injected through the central pipe while air-coflow is injected through the outer annulus. The bluff body acts as a flame stabilization device because the flow recirculation that it creates. In the inert configuration constitutes an interesting test case for assessing the capability to capture recirculating flows.

A structured cylindrical grid is used in NGA with a resolution of 192 $\times$ 92 $\times$ 32 in the axial, radial and azimuthal direction. The same grid is used in OpenFOAM except near the centerline. The mesh used in OpenFOAM is illustrated in Figure 5. While the axial location of the cell exactly matches (Figure 5), the radial location is slightly shifted for the OpenFOAM case because of the particular treatment at the centerline. These differences are deemed minor and the grids are considered identical for following discussion.

The spatial and temporal schemes used in NGA and OpenFOAM are of the same order. However, as opposed to OpenFOAM, NGA uses a staggered arrangement of variables [33] which ensures exact kinetic energy conservation when the mesh is uniform. In OpenFOAM, numerical dissipation of KE occurs because of the collocated arrangement of variables and because of the numerical discretization of the momentum equation. The latter will be addressed in Sec. 4. The baseline OpenFOAM solver is non-kinetic energy conservative (non-KEcons solver).

To illustrate the difference in the mixing process, a mixture fraction variable [45] is transported using the following equation.

[TABLE]

Due to the absence of source terms and the uniform mass diffusivity, the mixture fraction is a conserved variable.

At the fuel inlet boundary, a mixture fraction $Z_{mix}$ is set to 1 and the axial velocity $U$ to 54.8 m/s. At the coflow boundary, the mixture fraction is set to 0 and the axial velocity is set to 20 m/s. The timestep is held constant at 5e-7 s and the turbulence model is disabled in order to remove any differences which would stem from the subgrid scale (SGS) model. The dynamic viscosity and the mass diffusivity are held constant to a molecular value set to $1\times 10^{-5}$ kg.m*-1*.s*-1*. The jet diameter is set to 3.6 mm.

Figure  6 shows the comparison of the mean mixture fraction fields and the root mean square (RMS) axial velocity, which is a direct measurement of the resolved turbulent kinetic energy. Although the mean structures are recovered, it is apparent that the RMS velocity is overall underpredicted compared to NGA, implying that the turbulent fluctuations are numerically damped more in OpenFOAM. The mean mixture also shows similar pattern as NGA, but the jet seems to mix at a slower rate.

3.2 Theoretical background

3.2.1 Continuous and discrete operations

While kinetic energy and mass conservation are uniquely defined in the context of continuous fields, such conservation is valid only in a specified sense in discrete representations. In order to prescribe the notion of conservation, a description of the discrete operation is necessary.

The ultimate goal of a numerical simulation is to be able to approximate the continuous functions $\boldsymbol{u}(\boldsymbol{x},t)$ and $p(\boldsymbol{x},t)$ which satisfy Eq. 1 and Eq. 2. The continuous solutions also satisfy other equations obtained from the manipulation of the mass and momentum equation. For example, applying the curl operator to the momentum equation results in the vorticity conservation equation. In this work the dot product of the momentum equation with the velocity is considered. It results in a continuous sense in the KE transport equation written below

[TABLE]

The discrete fields obtained numerically satisfy the volume integrated momentum conservation equation and mass conservation equation, according to the numerical procedure outlined in Sec. 2.2. The velocity at the cell centers and the pressure at the cell centers are defined as fields respecting the momentum conservation equation and the mass conservation equation.

Using a similar strategy as the one used for the continuous fields, the discrete KE equation could be derived from the discrete momentum equation. It can be found that the continuous operations (like summation by parts) used to derive Eq. 14 are not necessarily valid discretely and lead to extra-terms which act as sink terms or source terms in each cell of the domain [13, Ch. 10.2]. Two observations can be made: 1) satisfying the discrete momentum equation does not ensure that the KE equation is satisfied 2) the fact that velocity and pressure were solely defined as satisfying momentum and mass conservation does not leave room for introducing another conservation equation. Therefore, the conservation of other quantities like KE must come from an appropriate choice of the discretization schemes, and should be indirectly enforced.

Let $K$ be one form of KE for which a transport equation can be deduced from the already enforced momentum and mass conservation equations. An energy-conservative scheme guarantees that for each internal face (as opposed to boundary face) of the domain, the amount of $K$ being shared from a computational cell to a neighbour is the exact opposite of the one being shared from this neighbour to this cell.

As opposed to this definition, with non-energy-conservative schemes, one can find an internal face for which the property cited previously is not ensured. Two types of non-energy-conservative schemes can be then defined. A dissipative scheme is a scheme which removes KE from a control volume. An unstable scheme is a scheme which creates KE in a control volume. Therefore, non-energy-conservative schemes do not necessarily lead to more numerical stability through dissipation. Demonstrating that a scheme always dissipate KE (upwind schemes for example or certain discrete pressure gradient formulations [29]) is a considerable achievement in ensuring numerical stability.

The ability to conserve quadratic properties does not imply better accuracy in predicting primary quantities (velocity, pressure etc.). In fact, it will be seen in Sec. 4 that some first order interpolation schemes are better suited than second order centered schemes to ensure KE conservation. Hence, a trade-off between accuracy and dissipation needs to be made. Prior implementations have used both high-order spatial schemes and lower order energy conservative spatial schemes [46, 47, 48]. For example, the compressible AVBP code from the CERFACS research group [49, 46] has second and third order schemes capabilities, but high-order schemes are usually preferred [50]. For low-Mach number solvers, high-order capabilities are sparse (for example YALES2 solver from the CORIA research group [47, 51]) and also exist alongside energy conservative schemes (See [52] for the YALES2 solver, and [29, 48] for the CDP solver from the CTR research group).

With this background, the goal now is to develop a solver in the OpenFOAM framework that is able to minimize KE losses in the sense described above.

4 Minimal KE dissipation in collocated mesh formulations

The OpenFOAM solver, similar to be many unstructured mesh formulations, involves collocated variables. Here, all variables (pressure, velocity, density and scalar) are located at the center of a control volume. The alternative approach is the staggered grid technique [53], where the velocity components are located on the surfaces of the control volume. In collocated mesh solvers, complete KE conservation is not feasible due to irreducible interpolation errors (described below and in [29, 54]). Consequently, the objective here is to implement schemes that minimize energy dissipation.

Schemes that conserve KE are designed such that by taking the dot product of the discrete momentum equation, a discrete KE equation can be obtained with minimal artificial source terms. Here the schemes derived in [26] are chosen for the time derivative and the convective term of the momentum equation. The justification of this choice is provided in Sec. 4.1. As for the pressure gradient term, the scheme introduced in [30] for collocated arrangement of variables is used. Finally the pressure Laplacian is treated following [29].

4.1 Convective term and time-derivative term

The energy conservation statement for the time-derivative and convective term can be written in the continuous form as:

[TABLE]

The goal then is to obtain a discrete approximation that is consistent with this continuous representation. In this regard, skew-symmetric formulations have been widely used [16, 55]. Achieving discrete energy conservation for constant density flows was done using first interpolations on non-uniform grids (interpolation denoted as midpoint in the following) [56, 29, 30]. A trade-off is therefore made on accuracy to attain energy conservation. In the case variable density flows, the derivations referenced above are not applicable anymore because they require the velocity field to be solenoidal. As a result, extra-terms proportional to the density change appear in the KE equation [32]. Recently, the skew-symmetric schemes were again used to ensure KE conservation of variable density flows in collocated formulations [26]. This formulation removed the requirement of the velocity field to be solenoidal. It should be noted that this formulation still relies on midpoint interpolations of the velocity from the cell centers to the cell faces.

The skew-symmetric form of the momentum convection term used for the derivations of energy conservative properties is written as follows:

[TABLE]

where

[TABLE]

with $\frac{\delta_{1}}{\delta_{1}x}$ the second order finite volume approximation of gradients, $\frac{\delta_{1}}{\delta_{1}t}$ the Euler scheme, $\overline{(.)}^{t}$ the midpoint time interpolation and $\overline{(.)}^{x}$ the midpoint spatial interpolation as defined in [26].

In this approach, a divergence form equivalent to the skew-symmetric form can also be derived which makes this scheme well suited for finite volume methods. The equivalence between the divergence and the skew symmetric form relies on the enforcement of the discrete continuity equation given by:

[TABLE]

which is used to formulate the pressure-correction Poisson system (see C).

As a side note, in many combustion applications, the presence of liquid (e.g., spray combustion) or solid (e.g., soot formation) phases lead to source/sink terms in the continuity equation representing exchange of mass between phases. It is shown in B that the above formulation remains energy-conservative for such modified continuity equations as well.

The form of KE conserved by the schemes is written as

[TABLE]

As stated in Sec. 3.2, this is a non-trivial form of KE which needs to be precisely defined when deriving or describing a solver that conserves secondary properties.

4.2 Pressure term

In the context of minimizing energy dissipation for variable density flows, the pressure gradient term of the momentum equation has not received particular attention. The approach used here relies on the derivations done for constant density flows [30, 29]. Therefore these schemes ensure that far from density variations in the domain, the KE dissipation is minimized. As noted elsewhere [29, 30, 54], collocated formulation always introduce an extra-term in the KE budget. The only way of dealing with this term is to make it as small as possible. A least square formulation has been proposed by Mahesh [30] in the computation of the pressure gradient and the application of the pressure gradient to the cell center. However it resulted in instabilities in the solver [29, 30]. Other methods have been proposed recently but require intrusive modifications to the solver [57]. Here, the approach of Ham and Iaccarino [29] is employed. The discretization of the pressure laplacian is tailored to ensure that the energy dissipation scales as $\mathcal{O}(\Delta t\Delta x^{2})$ for skewed meshes. This derivation remains valid in the context of variable density flows. In order to apply the aforementioned methods, the pressure gradient term should be discretized using a midpoint scheme.

4.3 Scalar transport equation schemes

Similar to KE for velocity, a quadratic conservation property can be defined for scalars based on the square of the scalar variables. However, scalar equations do not have an accompanying pressure equation. It is known that nonlinear convection terms can lead to the formation of dispersive waves for any discretization scheme of order more than 1. Consequently, the use of midpoint scheme, which will reduce to a second-order central scheme, will cause the generation of location maxima and minima outside the bounds of the scalar [58]. In general, some form of numerical dissipation is necessary to ensure stability of the solution. While several techniques have been formulated [59, 60, 61, 62, 58], none of these schemes will preserve the quadratic properties associated with scalars.

Here, the focus is on the temporal scheme used for the discretization of Eq. 3. In the solver description, the momentum equation advances the velocity field from the time $n$ to the time $n+1$. The solver procedure involves outer-iterations to couple the velocity and the scalars like mixture fraction or density. In order to make the procedure more implicit, it is useful to use the newly estimated velocity in the formulation of the scalar transport equation. Therefore the scalars are fully transported with the most updated velocity field. Using a centered time-scheme, the scalars are advanced from the time $n+1/2$ to the time $n+3/2$ as shown in Fig. 7. Since the velocity is known at the time $n+1$, this allows to use a central time scheme for the scalar transport equation. The solver is said to be time-staggered.

In variable density cases, a time-staggered solver allows the use of a second order scheme for time variation of density field. Without time-staggering, the density would be evaluated only at time $n+1$ and $n$. The pressure-Poisson equation is derived using the $\nabla\cdot(\boldsymbol{u}^{n+1})$ which involves $\frac{\delta_{1}\rho^{n+1}}{\delta_{1}t}$ (See Sec. 4.4.3 for a complete derivation). In a non time-staggered scheme $\frac{\delta_{1}\rho^{n+1}}{\delta_{1}t}=\frac{\rho^{n+1}-\rho^{n}}{\Delta t}$ which is a 1st order approximation in time. In the time-staggered version $\frac{\delta_{1}\rho^{n+1}}{\delta_{1}t}=\frac{\rho^{n+3/2}-\rho^{n+1/2}}{\Delta t}$, which is a $2^{nd}$ order approximation in time.

Similar to the other terms of the scalar transport equation, the scalar source term $\dot{\omega}$ is written using a time interpolation. Since the source term is a function of scalar values, it is written as:

[TABLE]

4.4 Implementation in OpenFOAM

4.4.1 Spatial schemes

As highlighted in Sec. 4.1 and Sec. 4.2, midpoint schemes (Gauss midPoint in OpenFOAM) are better suited than second order central schemes (Gauss linear in OpenFOAM) for the convective term of the momentum equation and the pressure gradient. In addition, in order for the above derivation to be valid, the mass flux used in the convective momentum term should satisfy the continuity equation in an integral sense.

The Laplacian computation of the pressure follows the approach of Ham and Iaccarino [29]. To take into account non-orthogonality of the computational mesh, the cell-face pressure gradient is obtained using a combination of cell-center pressure values, and cell center pressure gradient values. This family of laplacian schemes is the corrected laplacian schemes in OpenFOAM. The ratio of orthogonal to non-orthogonal contribution can be adjusted using the limited parameter [7, Ch. 3.3.1.3]. Since the non-orthogonal part of the laplacian is treated explicitly, the pressure equation is repeated nNonOrthogonalCorrectors$+1$ times to account for this explicit part.

The viscosity terms in the momentum equations dissipate KE. Hence, their discretization is governed primarily by truncation errors in spatial and temporal directions. To ensure temporal consistency, these terms are evaluated at $n+1/2$ in the time-staggered scheme, using a midpoint interpolation in time.

In OpenFOAM syntax, the final momentum equation takes the form.

rhoTimeInterp = 0.5*(rho+rho.oldTime()); rhoTimeInterp.oldTime() = 0.5*(rho.oldTime()+rho.oldTime().oldTime());

sqrt_rhoTimeInterp = sqrt(rhoTimeInterp); sqrt_rhoTimeInterp.oldTime() = sqrt(rhoTimeInterp.oldTime()); sqrt_rhoTimeInterp_timeInterp = 0.5*(sqrt_rhoTimeInterp+sqrt_rhoTimeInterp.oldTime());

fvVectorMatrix UEqn ( fvm::ddt(rhoTimeInterp,U)

• (sqrt_rhoTimeInterp/sqrt_rhoTimeInterp_timeInterp)* fvm::div((phi+phi.oldTime())/4.0, U)
• (sqrt_rhoTimeInterp.oldTime()/sqrt_rhoTimeInterp_timeInterp)* fvc::div((phi+phi.oldTime())/4.0, U.oldTime())

• fvm::laplacian(0.5*turbulence->muEff(),U)
• fvc::laplacian(0.5*turbulence->muEff()*U.oldTime())
• fvc::div(turbulence->muEff()0.5dev2(Foam::T(fvc::grad(U))))
• fvc::div(turbulence->muEff()0.5dev2(Foam::T(fvc::grad(U.oldTime())))) );

4.4.2 Temporal schemes

Regarding energy dissipation, there exist a constraint on the convective term: conservation is fully achieved only when the velocity convected is the same as the velocity updated. This can be satisfied using some semi-implicit treatment in the algorithm (through outer-iterations). Given that this treatment is correctly implemented, the solver structure discussed above is equally valid for implicit or explicit Euler timestepping.

The temporal treatment of the pressure formulation needs to be further discussed, since this can affect the computational expense. More precisely, the choice is between the PISO procedure, which is already implemented in OpenFOAM, and the fractional timestep approach that is used by most research codes [33, 35][63, Ch. 3.2.2.1]. Numerical differences between PISO and the fractional timestep method are briefly described in A.

The PISO procedure involves twice the number of pressure-Poisson solutions as the fractional timestep method, and is therefore twice as expensive. However, the PISO procedure allows the momentum equations to be treated fully implicitly, which allows the use of larger timesteps. It has been argued [13] that the PISO procedure was tailored for steady-state problems. This has to do with the better stability properties of the PISO algorithm. When large CFL numbers can be used (for instance, in time-marching to steady-state), the PISO procedure becomes attractive. However, for transient problems where time accuracy is important, a large timestep is not useful and the PISO procedure has no added value compared to the fractional timestep method. Furthermore, for reacting cases, outer-iterations that couple scalar equations are used and are likely to be redundant with the pressure iterations involved in the PISO algorithm. Since their impact on the KE conservation property is not significant (Sec. 5.1), the final application will determine the approach needed. The PISO implementation is described here and can be easily extended to fractional timestep using details in A and C.

4.4.3 Pressure Poisson equation

The final piece of this low-Mach number solver is the solution of the pressure-Poisson equation. Based on the description in Sec. 4.2, the pressure equation has to be written such that the discrete continuity equation is enforced after the correction step. The derivation of the Poisson equation is detailed in C. The consistent pressure equation can be written as

[TABLE]

where $a^{-1}$ is the component at the cell considered of the inverse of the split operator $[A]$. In the OpenFOAM framework it translates to

phi = fvc::interpolate(rhoTimeInterp*U) & mesh.Sf();

fvScalarMatrix pEqn ( fvm::laplacian(rhoTimeInterp*rAU,p)

• fvc::ddt(rho)
• fvc::div(phi) );

5 Verification of the minimally KE dissipative solver implementation

The solver verification is conducted in various configurations with increasing complexity. First, the minimal dissipation property of the solver is assessed using an inert periodic configuration. Using the upgraded solver, the inert bluff-body flame configuration which was presented in Sec. 3.1 as a motivation for improving the KE conservation property of the solver is revisited. This test is used to assess the conservation properties in flows with unsteady inflow, non-periodic boundary conditions, and complex flow patterns with a range of length scales. Then a variable density case using periodic boundaries is studied. Finally a variable density case in a complex geometry using multiphase physics is studied.

5.1 Constant density Taylor-Green (TG) Vortex Case

The steady Taylor-Green vortex problem is used to test KE conservation properties of the new solver. This solver is called KEcons solver as opposed to the baseline non-KEcons solver available in OpenFOAM. The initial flow field consists in a series of vortices inside a periodic domain. The viscosity is set to zero so that any dissipation is purely numerical. Two cases are investigated, one with a uniform and isotropic grid, and another with a skewed mesh. The parameters chosen for the size and timestepping of the problem can be found in [29]. The flow described by the following set of equations:

[TABLE]

[TABLE]

[TABLE]

for $x\in[-1;1]$ and $y\in[-1;1]$. For the first test case, the grid is regular and structured. Similar to [29], a second case is investigated where the grid has undergone a skewing operation defined by

[TABLE]

[TABLE]

The skewed grid is shown in Fig. 8.

Since the solution is time-invariant, KE should be preserved. However, errors in discretization will lead to change in KE. The proposed method leads to decay of KE, albeit at a slow rate, compared to existing OpenFOAM formulations without the pressure correction and interpolation schemes. The rate of decay is presented against results from Ham et al. [29], where a similar collocated solver was used for the same test case.

The structured grid solution shows little difference between the conservative and non-conservative formulations (Fig. 9). This is expected, since energy loss is primarily associated with the interpolation schemes of the convection term and the pressure laplacian formulation. Uniform structured grids minimize such interpolation errors, especially when the flow is well-resolved using small computational cells. In this case, linear interpolation adequately captures spatial variation. However, skewed mesh calculations show larger difference between conservative and non-conservative solvers (Fig. 9). It is seen that the OpenFOAM formulation produces results consistent with the collocated arrangement solver of Ham et al. [29]. It is also seen that non-conservative approaches such as that implemented in the OpenFOAM icoFoam solver with no pressure correction schemes, lead to much higher KE loss.

This test case is also used to assess the influence of the timestepping technique (PISO or fractional timestep. See A for details about these methods). Figure  9 shows that both techniques lead to the same amount of numerical dissipation.

For a constant density case, the definition of $\boldsymbol{\hat{u}}$ in Sec. 17 reduces to $\frac{\boldsymbol{u}^{n}+\boldsymbol{u}^{n+1}}{2}$. The definition of KE at time $n$ becomes $\frac{\boldsymbol{u}^{n}\cdot\boldsymbol{u}^{n}}{2}$. In this case, the budget of the discrete KE equation integrated over the simulation domain can be used to analyze the contribution of the individual terms. The contribution of the convective term to the energy dissipation budget is plotted in Fig. 10. This KE error is called the skew-symmetric error. It is first observed that in both the PISO and the fractional timestep methods, the contribution of the convective term to the total energy dissipation is negligible and is therefore dominated by the pressure term $\sum_{P}V_{P}(\hat{u}_{i}\cdot\frac{\delta_{1}P}{\delta_{1}x_{i}})$. The time derivative term $(\hat{u}_{i}\cdot\frac{\delta_{1}u_{i}}{\delta_{1}t}-\frac{u_{i}^{2,n+1}/2-u_{i}^{2,n}/2}{\Delta t})$ is exactly zero by construction. Here, $V_{P}$ is the local cell volume and $P$ is the identifier of the cell. To have complete cancellation of the KE dissipation due to convection, it is essential to ensure that mass is exactly conserved to a high precision in the entire domain. In the skewed case, using the same tolerances for the pressure equation, the mass conservation error was higher.

The TG test case is also solved using an unstructured grid, which is representative of the practical configurations that will be simulated using this solver. To ensure direct comparison, the number of triangles in the tetrahedral case is set equal to the number of hexahedrons on the boundaries. The results from these simulations are shown in Fig. 11. It is observed that KE dissipation is considerably larger than in the structured and skewed mesh cases (Fig. 9). In fact, the rate of dissipation increases with time initially, but reaches a constant value after this initial transient. This leads to a concave shaped dissipation plot. These results were also observed elsewhere [64], and suggests accumulation of numerical error. It is interesting to note that the dissipation of energy does not improve with a reduction in timestep, implying that the dissipation error obtained from the interpolation of the pressure gradient at the cell faces to the cell centers is $\mathcal{O}(\frac{1}{\Delta t})$. However, the results are improved compared to the original icoFoam solver.

5.2 Non-periodic boundary verification

The construction KEcons solver involved several modifications which led to improved energy conservation properties of the OpenFOAM solver. One could ask which one these modifications is the most critical. In Fig. 12, the baseline icoFoam solver using pressure gradient corrections is compared to the KEcons solver and the baseline icoFoam solver. It appears that most of the improvements between the described formulation and the baseline OpenFOAM solver are due to the pressure scheme. In order to observe the improvements due to the formulation of the convective scheme as well as the time-staggering, a different test problem is studied. Similar to the discussion in Sec. 3.1, the new formulation is compared to the pisoFoam formulation using linear schemes and the corrected scheme for pressure. Therefore, the changes observed for the velocity field will reflect the influence of the convective term in the new momentum equation formulation. The changes observed in the mixture fraction field will reflect the improvements due also to the convective term formulation and time-staggering implementation.

The mean mixture fraction field and the RMS axial velocity are plotted against the NGA results in Fig. 13. The modified low-Ma solver shows visibly improved results in terms of the core jet length, as compared to the non-conservative solver (Fig. 6). This is consistent with the improvements seen in the TG vortex cases (Sec. 5.1). The radial profiles of these statistics are plotted in Fig. 14. Overall, the new solver produces results closer to the structured grid NGA solver, with the maximum improvement near the centerline.

5.3 Variable Density TG Vortex Case

The TG case can be modified to include density changes by including an additional transport equation. This manufactured method solution (MMS) technique corresponds to problem 3 in [43]. The relevant parameters (see [43] for details) are: $\rho_{0}=5;\rho_{1}=1;k=2;\omega=2;u_{F}=0;v_{F}=0;\mu=0.001$. Unlike the steady Taylor-Green vortex case in Sec. 5.1, the flow field oscillates at the same rate as the density field. Similar to the case studied in Sec. 2.4, the density field is directly obtained from a transported scalar field (analytical solution available in [43]), for which a source term is provided. No additional source terms are used for the momentum or mass conservation equations. The analytical velocity field is now given by the following set of equations:

[TABLE]

and

[TABLE]

For the solver to be accurate, it is important to capture the coupling between the density changes and the velocity field. Since the momentum equation does not have a direct forcing term, its temporal variations are a result of the density changes alone. For these tests, the maximum convective CFL number is held constant at 0.15 in a fashion similar to [43]. An instantaneous density field is shown in Fig. 15. It was first found that the variable density solver using the hybrid formulation of the pressure correction and the non-KEcons solver were prone to instabilities.

Using the new low-Mach number solver, second order accuracy in space could be obtained as shown in Fig. 16. This is consistent with tests of similar collocated numerical scheme presented in [43]. The same configuration is simulated using a 2D grid composed of triangular cells. Similar to the constant density TG case, this tetrahedral mesh based simulation is more applicable to practical OpenFOAM computations. The spatial convergence remains second order (Fig. 17), which is encouraging for practical application.

The results are displayed next to the ones generated by Shunn et al. [43] who used non-zero convective velocity parameters $u_{F}$ and $v_{F}$. The interpolation methods are also likely to be different since Shunn et al. [43] use a node-based code, meaning that variables are stored at the nodes rather than the cell centers. The discretization of the scalar source terms also plays a large role in the solver accuracy especially in MMS procedure. Depending on the time at which the source terms have been computed ($n$ or $n+1/2$), significant accuracy differences can be observed. One should therefore not infer any conclusions from the magnitude of errors between both solvers but rather focus on the rate of convergence for the error, which appears to be 2nd order in space.

Since the test case contains viscosity, and non-zero dilatation, KE will not be conserved. Instead, the solver should minimize the difference between the analytical and computed KE values. In order to verify KE conservation, the dot product of the time and convective terms of the momentum equation with $\boldsymbol{\hat{u}}$ is computed, and is effectively $\frac{\overline{\rho}^{t,n+1}(u^{2})^{n+1}-\overline{\rho}^{t,n}(u^{2})^{n}}{\Delta t}$. Any error comes from an inaccurate estimation of $u^{n+1}$ (explicit term of the momentum convection) or from a lack of primary mass conservation in the continuity equation. In Fig. 18, the contribution of this error is plotted along with the total KE error for the 1024 cells case using a timestep of $6.25ms$. This error is called the skew-symmetric error, similarly to Sec. 5.3. It can be observed that most of the KE error does not come from the skew-symmetric error but from other components of the equations solved. This gives confidence in the ability of the new methodology to ensure energy conservation for the time derivative and the momentum convection terms.

5.4 Complex geometry with density model and non-mass conservative flow

In order to demonstrate the solver performance in a practical geometry, combustion and soot formation in a pressured model aircraft combustor is simulated. The flow configuration is based on the experimental set up at DLR [65]. Here, the new solver is used to simulate combustor operation at 3 bar pressure.

The combustor is similar to the one used by Geigle et al. [66] and Koo et al. [67]. The combustor geometry is shown in Fig. 19 with typical streamline in the combustor from inlet to outlet. It was designed for operation at 10 kW/bar power and installed with large optical access for simultaneous laser acquisition of velocity, temperature, species mass fractions and soot volume fraction. It has a cross-sectional area of 68 $\times$ 68 mm2, with a height of 120 mm. The inflow consists of three concentric nozzles: two room temperature air inlets with swirling velocity and 60 annular straight channel fuel ($C_{2}H_{4}$) inlet in between the two air flows with a size of 0.5 $\times$ 0.4 mm2 each. A single constricted exit of diameter 40 mm removes the combustion products. At 80 mm height, four additional air ducts of 5 mm diameter inject secondary air into the combustor radially, meeting at the combustor central axis, forming a stagnation point. Due to the high pressure, combined with locally rich fuel-air conditions, significant soot formation is observed. Estimation of soot profiles is the main target of these computations. Due to the extreme sensitivity of soot formation to local thermochemical conditions, as well as the trajectory of the fluid particles within the combustor, the use of temporally accurate LES solvers become important. In this study, the specific case of 3 bar operation with 460.3 slpm of primary air, 39.3 slpm of ethylene fuel and 187.4 slpm of secondary oxidation air is considered. The global equivalence ratio is approximately 0.86, which is considerably below the sooting limit. Hence, any soot particle observed is generated due to local inefficiencies in mixing and oxidation of the fuel.

The computational domain consists of approximately 12 million tetrahedral cells, with grid refinement applied near the inlets and the near-wall region (Fig. 20). In LES, it is necessary to resolve the large-scale structures, and a metric for such refinement is the Pope criterion [68]. A fractional energy $M$, defined as the ratio betwen the sub-filter KE and the total KE is used to determine the resolution adequacy. Since the sub-filter KE is not directly available from the resolved fields, a model is used to estimate this quantity. The cell-size is refined until this ratio is below 0.30 everywhere in the domain. Note that even with this requirement, there might be instantaneous $M$ values that are above the cut-off threshold. Figure 20 shows the $M$ field plotted for the final grid used in this study.

The low-Mach number LES equations are solved along with a flamelet-based model for turbulent combustion, and moments-based model for soot particle evolution [69]. This requires additional transport equations for mixture fraction, progress variable, enthalpy, and soot moments. details of the model can be found in [67, 70].

The LES governing equations are obtained by Favre-filtering the momentum and scalars including soot moments. Unclosed terms in the subfilter flux are closed with the dynamic SGS model [71]. Filtered chemistry reaction terms for the enthalpy equation are closed using the presumed-PDF approach described in Mueller et al. [69]. For these tests, the timestep is held constant at 5 $\times$ 10*-7* s. The LES simulation was performed on 2048 cores with time-averaged data taken after 10 flow-through times, totaling approximately 80 wall-clock hours.

Figure 21 shows time-averaged two-dimensional images of axial velocity, tangential velocity and mixture fraction, obtained using the minimally-dissipative solver. It can be seen that the flow structure consists of a large inner recirculation zone (marked by negative axial velocity component), as well as a smaller but persistent outer recirculation zone. The size and strength of these recirculation zones determine the level of mixing, and the tendency of the combustor to form soot particles.

Figure 22 shows time-averaged velocity statistics compared against experimental data. The time-averaged axial velocity is well-captured by the solver, including the presence of the high-velocity regions along the fuel injection path. Further, the low-velocity recirculation region is also captured accurately. Similarly, the mean tangential velocity shows very good agreement with the experimental data. In particular, the swirl component that determines the flame stabilization mechanism is well represented by the LES solver.

Perhaps, of more interest is the root-mean square (RMS) velocity, which is directly tied to the ability of the LES solver to maintain turbulent fluctuations from dissipating through numerical inaccuracies. Figure 23 shows that at all axial locations considered, the RMS velocity of the both the axial and tangential components are well captured. This provides confidence that the method is not dissipative even for complex geometries, and the turbulent structures remain accessible on unstructured grids as well.

6 Conclusion

A minimally-dissipative low-Mach number solver has been implemented in OpenFOAM. It was established that the default hybrid pressure correction approach leads to significant errors in mass conservation in a discrete sense, and was numerically unstable for certain flow conditions. The skew-symmetric formulation of the low-Mach number solver ensured discrete conservation of mass, while reducing dissipation of quadratic quantities. For variable density and reacting flows problems, a formulation for coupling scalar transport equations was also implemented. Verification cases demonstrated that the new solver significantly reduces KE dissipation even for highly skewed meshes. Further, the method is shown to be second-order accurate for variable density flow problems in structured and unstructured mesh cases.

Since LES computations are highly sensitive to numerical errors, ensuring such primary and secondary conservation properties will help isolate modeling and numerical errors, leading to more useful model validation exercises. The algorithm developed will benefit from formulations for a suitable pressure gradient treatment on collocated meshes. In conclusion, the current study establishes a reliable and open source solver for LES with applications in complex turbulent reacting flows.

Appendix A Distinction between PISO and fractional timestep method

The fractional timestep [34] and PISO [72] algorithm follow the same idea, which relies on enforcing mass conservation with a pressure correction term. Below, the distinction between the two approaches is described.

The fractional timestep method is the most intuitive method but the least-implicit approach, in the sense that not all terms of the momentum equation are updated using the same velocity vector. First, the velocity at the new timestep is guessed (and is denoted by $\boldsymbol{u}^{g}$) using the most updated estimate for the pressure term. Then a discrete momentum equation is written using $\boldsymbol{u}^{g}$ everywhere except in the time derivative term. The viscous term is dropped without changing the purpose of the discussion.

[TABLE]

The convection and dissipation terms are written implicitly, leading to an equation for $\boldsymbol{u}^{g}$, which is an estimated velocity field. Using $\boldsymbol{u}^{g}$ in the convection term, the fractional timestep procedure can be progressed as:

[TABLE]

This equation is used to obtain $\boldsymbol{u}^{*}$ and subsequently, a pressure equation is formulated as

[TABLE]

The above equation can be used to solver for pressure, as well as the velocity field $\boldsymbol{u}^{n+1}$.

The PISO procedure follows the same idea but splits the derivative operator between a part treated implicitly and a part treated explicitly. The first step of the method is exactly similar to Eq. 29. The following correction step is formulated differently. Each operator (convection, diffusion, time) is treated as a sparse square matrix of size equal to the number of cells. This matrix is split between components that will be treated implicitly ($A$) and explicitly ($H$).

The pressure correction procedure is initialized as

[TABLE]

where $A_{time}$ denotes the implicit part of the time derivative operator, $A_{conv}$ and $H_{conv}$ denote respectively the implicit and explicit part of the convection operator, and $[\boldsymbol{u}]$ denotes the full velocity field.

The fractional velocity $\boldsymbol{u}^{*}$ is then corrected into a mass conservative velocity $\boldsymbol{\overline{u}}$ using a pressure correction $[\nabla p_{1}]$ field.

[TABLE]

In order to achieve second order accuracy in time, another layer of correction is needed. Therefore next, $\boldsymbol{\overline{u}}$ is used as a guess in the exact same way $\boldsymbol{u}^{g}$ was used.

[TABLE]

[TABLE]

Comparing Eq. 35 and Eq. 31, it is clear that the momentum equations are advanced differently in the PISO method and the fractional timestep method. The difference lies in the implicit treatment of some terms and the additional pressure correction layer in the PISO method.

Appendix B Application to non-mass-conserving systems

Many reacting flow applications involve exchange of mass between different phases, especially in the context of dispersed-phase flows [73, 74, 75]. In such cases, the continuity equation contains a source term that represents the mass added or removed from the continuous gas-phase. Here, the applicability of the low-Mach number minimally-dissipative solver to such systems is investigated.

Assuming that mass addition only affects the continuity equation, the modified mass conservation equation can be written as:

[TABLE]

where $\dot{S}$ denotes the mass exchange source term. For the sake of simplicity, pressure and viscosity are not included without any loss of generality. The continuous momentum equation can be written as

[TABLE]

Multiplying by $\boldsymbol{u}$, the KE equation can be written as

[TABLE]

Based on the continuity equation, the last term is non-zero and equal to $\dot{S}$. Following Morinishi [26], the momentum equation is still written in discrete conservative form as [26, Eq.128] and becomes:

[TABLE]

The skewed-symmetric form varies from the divergence form, only by the last term on the RHS. Its contribution to the conservation of KE is then similar to the continuous equation and is a fortiori proportional to the evaporation rate.

Appendix C Pressure Poisson equation

The pressure Poisson equation formulation is derived here for both the fractional timestep and PISO schemes.

In the fractional timestep method, the momentum equation in discrete form is given by:

[TABLE]

, while the continuity equation is written as:

[TABLE]

Applying the divergence operator to the above equation leads to

[TABLE]

In the PISO formulation, the derivation is not as straightforward since it involves the splitting operators (see A). In what follows, the brackets $[.]$ denote fields as opposed to cell values. Each pressure correction results in the same following relation:

[TABLE]

Multiplying by $\rho^{n+1}$ leads to

[TABLE]

Applying the divergence operator and rearranging results in

[TABLE]

Recall that $\rho$ is defined at timesteps $n+\frac{3}{2}$ and $n+\frac{1}{2}$. Hence $\overline{\rho}^{t}$ is defined at timestep $n+1$. Noting that $\overline{\rho}^{t}U^{*}=\phi^{*}$, Eq. 21 is obtained.