Model formalism

The idea behind the zero-Mach limit formulation of the compressible fluid system of equations consists of operating a multi-scale perturbation analysis as a function of the Mach number.

This in turn decouples the flow from acoustics and can allow for much larger time-steps.

Considering a low Mach number flow, the different independent variables governing the fluid flow are expanded in terms of the Mach number (a very small parameter) [202]:

p=p⁽⁰⁾+ Ma p⁽¹⁾+ Ma²p⁽²⁾+O Ma³

, (5.1a)

u=u⁽⁰⁾+ Mau⁽¹⁾+O Ma²

, (5.1b)

ρ=ρ⁽⁰⁾+ Maρ⁽¹⁾+O Ma²

, (5.1c)

T =T⁽⁰⁾+ Ma T⁽¹⁾+O Ma²

. (5.1d)

Upon introduction of these expansions into the PDEs one finds the following equations at order zero:

∇p⁽⁰⁾ = 0, (5.2a)

∂tρ⁽⁰⁾+∇·ρ⁽⁰⁾u⁽⁰⁾ = 0, (5.2b)

∂_tρ⁽⁰⁾u⁽⁰⁾+∇·ρ⁽⁰⁾u⁽⁰⁾⊗u⁽⁰⁾+∇p⁽¹⁾

−∇·µ

∇⊗u⁽⁰⁾+∇⊗u^(0)T −2

3∇u⁽⁰⁾I

= 0 (5.2c)

p⁽⁰⁾ =ρ⁽⁰⁾rT¯ ⁽⁰⁾, (5.2d)

where p⁽⁰⁾ is the thermodynamic pressure, uniformly distributed in space, and p⁽¹⁾ the hy-drodynamic pressure admitting fluctuations in time and space. In practice, the acoustics are filtered by operating a Helmholtz-Hodge decomposition of the velocity vector field into a solenoidal and irrotational field. Reconstruction of purely hydrodynamic effects through a divergence-free velocity field allows to filter out acoustic effects while the second component

is determined through the EoS and the uniform thermodynamic pressure. The aim of the next section is to propose a LB-based model capable to represent this low Mach number formulation.

5.1.1.1 Decoupling of density: introducing hydrodynamic pressure

Over the past score years, a number of attempts have been made at developing models for combustion simulation based (partially or entirely) on LB solvers. A number of these models based their developments on low Mach (decoupled) formulations [37, 38, 43]. However none of these proposed models have gone beyond simple steady configurations. The present part of the thesis will introduce a pressure/density decoupling approach initially developed for simulation of multi-phase flows with large density ratios, see for instance [203, 204]. This formulation was later extended and used (under different assumptions) for 2-D thermo-compressible flows in [48, 205, 206]. The main idea behind the development of a LB-based solver for thermo-compressible flows, following the previously introduced philosophy of low Mach formulations is to decouple the flow density from velocity and hydrodynamic pressure space fluctuations. To achieve this decoupling, and given that the thermodynamic pressure tied to density through the EoS, equal to p = ρc²_s in the LBM, is assumed to be uniform in space, the first step is to introduce a hydrodynamic pressure independent from the local density, at the Euler level. This is achieved by introducing a body force defined as [44]:

F =∇ρc²_s−∇p_h+F_b, (5.3)

where Fb designates other external body forces such as gravity. While this body force modifies the pressure term at the Euler level, it does not eliminate the velocity-dependence of the local fluid density; Furthermore, since the hydrodynamic pressure does not appear in the Boltzmann equation (except in the introduced body force), it needs closure. To recover the hydrodynamic pressure as the zeroth-order moment and completely decouple the flow density from hydrodynamics the following new distribution function is introduced:

g_α⁰ =c²_sfα+wα ph−ρc²_s

. (5.4)

It can readily be observed that the zeroth-order moment of this new distribution function is p_h. Introducing this new distribution back into Boltzmann’s equation the following contin-uous time-evolution PDE is obtained:

∂_tg_α⁰ +c_α·∇g_α⁰ = 1 τ

g_α^(eq)

0

−g_α⁰

+w_α(∂_tp_h+c_α·∇p_h)

−wαc²_s(∂tρ+cα·∇ρ)−F ·∇cαfα, (5.5) where ∂_tρ+c_α·∇ρ can be re-written using the continuity equation as:

∂_tρ+c_α·∇ρ= (c_α−u)·∇ρ−ρ∇·u. (5.6)

Furthermore, the last term on the RHS of Eq. 5.5 can be approximated as [204]:

F ·∇_cαf_α =−F · c_α−u

ρc²_s f_α^(eq), (5.7)

and terms of the form∂_tp_h andu^j·∇p_h, withj ≥1, can be neglected as they are third-order in non-dimensional velocity. Eventually, before discretization in space and time the following PDE is recovered:

∂_tg_α⁰ +c_α·∇g_α⁰ = 1 τ

g_α^(eq)

0

−g_α⁰

+c²_s fα^(eq)

ρ −w_α

!

(c_α−u)·∇ρ+w_αc²_sρ∇·u +F_b·(c_α−u)fα^(eq)

ρ . (5.8) 5.1.1.2 Introducing thermo-compressibility: Evaluation of velocity divergence The previously introduced formalism decoupled the density from the flow field, i.e. velocity and hydrodynamic pressure field. However, it also eliminated the implicit equation of state in the original LB formulation. Furthermore, the zeroth-order moment of the distribution function being the hydrodynamic pressure, the continuity equation is not imposed in the model anymore. Following the original low Mach formulation, the density is now computed using the thermodynamic pressure, local temperature and EoS as:

ρ= p_th

¯

rT. (5.9)

Furthermore, for the model to satisfy the continuity equation, it is used to evaluate the velocity divergence appearing in Eq. 5.8 as:

∇·u=−∂_tρ+u·∇ρ

ρ , (5.10)

and further developed using the EoS as:

∇·u =−∂_tp_th+u·∇p_th

p_th − ∂_tT¹ +u·∇_T¹

1 T

− ∂_t¹_r¯+u·∇¹_r¯

1

¯ r

. (5.11)

The space-derivative of p_th cancels out as it is assumed to be uniform. Furthermore the second term on the LHS can be re-written as:

− ∂_tT¹ +u·∇_T¹

1 T

= ∂_tT +u·∇T

T , (5.12)

while the third one is readily developed as:

− ∂_t¹_r¯+u·∇¹_r¯

1

¯ r

=

Nsp

X

k=1

M¯ Mk

(∂_tY_k+u·∇Y_k), (5.13) where we have used ¯r =RPNsp

k=1Y_k/M_k. Putting all these terms back together one gets:

∇·u=−∂tpth

p_th + ∂tT +u·∇T

T +

Nsp

X

k=1

M¯

M_k(∂_tY_k+u·∇Y_k), (5.14) for a multi-species flow. The previous equation along with Eqs. 5.8 and 5.10 and:

X

α

g_α =p_h, (5.15a)

X

α

c_αg_α =ρc²_su, (5.15b)

constitute the main equations for the LB-based low Mach model for thermo-compressible flows.

Using the CE formalism introduced in subsection 3.1.1, it can readily be observed that at order ε¹ the following PDE is recovered for the hydrodynamic pressure:

1

ρc²_s∂_t⁽¹⁾p_h+∇⁽¹⁾·u=−∂_t⁽¹⁾pth

p_th +∂_t⁽¹⁾T +u·∇⁽¹⁾T

T +

Nsp

X

k=1

M¯ M_k

∂_t⁽¹⁾Y_k+u·∇⁽¹⁾Y_k

. (5.16) While this LB-based low Mach formulation does not exactly impose the intended velocity divergence, the error term _ρc¹2

s∂_t⁽¹⁾p_h is of order O( ^u3

(δx/δt)³) and therefore negligible for small non-dimensional velocities. This argument is further comforted by the restrictive CFL condi-tion on different collision operators (especially SRT) in the limit of vanishing non-dimensional viscosities (relaxation coefficients).

5.1.1.3 Space and time-discretized equations

Given that the space discretization process relying on integration along characteristic lines has been thoroughly detailed in subsection 2.4.2 for the classical LB formulation, re-deriving the discrete equations for the previously-introduced model would be redundant. As such only the final equations are given here. The discrete collision-streaming equation for the new distribution function is:

g_α⁰ (x+c_αδ_t, t+δ_t)−g⁰_α(x, t) = δ_t τ

g^(eq)

0

α(x, t)−g_α⁰ (x, t) +

1− δ_t

2τ

w_αc²_sρ∇·u

+

1− δ_t 2τ

(c_α−u) fα^(eq)

ρ −w_α

!

c²_s∇ρ+F_α,b, (5.17)

where we have dropped the overbar of the post-discretization distribution function for the sake of simplicity andF_α,b is the contribution of the external body forces, e.g. gravity, which can be evaluated through anyone of the available LB forcing schemes. It is interesting to note that apart from the presence of species, the model as used here is different from those proposed in [48, 205] as it includes a factor 1−2τ^δt. The absence of this factor, also needed for the correction term of the compressible scheme in the next section, leads to serious stability issues. The moments of the distribution function are defined as:

X

α

g_α =p_h− c²_sδ_t

2 (ρ∇·u+u·∇ρ), (5.18a)

X

α

c_αg_α =ρc²_su− c²_sδ_t

2 F_b. (5.18b)

It is also worth mentioning that this formulation where the temperature is defined through the EoS, which would need the local temperature and composition, can not be used in combi-nation with solvers for the conservative form of the energy and species mass balance equations with explicit coupling. Given that in the conservative form extensive forms of variables are transported, to get to the intensive parameters one needs the local density. The local den-sity on the other hand in this low Mach formulation is unknown and determined through the EoS and local intensive variables, i.e. temperature and mass fractions. Therefore, to only correct way to couple these two solvers is implicit, and involves an iterative solver. As such for the remainder of this study, the low Mach formulation is always used in combination with a finite-difference solver for the non-conservative form of the energy and species balance equations. Furthermore, given that in this model the thermodynamic pressure is supposed to be uniform in space, for energy balance, the sensible enthalpy formulation is used and the space gradient of pressure is set to zero.

Im Dokument Development of a lattice Boltzmann-based numerical method for the simulation od reacting flows (Seite 130-134)