Figure 6.19: Fuel mass fraction fields at three different times, i.e. (from left to right): 0, 0.08, 0.16, 0.24 and 0.32ms for case 1.
Simulations performed in this section showed that the LB-based numerical schemes con-sidered in the present manuscript can be used to model low Mach number combustion in multi-dimensional configurations, for both premixed and diffusion flames, and can correctly capture flame/flow interaction.
6.3. Concluding remarks
t = 0.276ms
t = 0.456ms
t= 0.684ms
t = 0.78ms
t= 0.168ms
t= 0.3ms
t= 0.48ms
t= 0.66ms
t= 0.084ms
t= 0.18ms
t= 0.372ms
t= 0.54ms
2230
800 T[K]
Figure 6.20: Instantaneous snapshots of temperature field as obtained for different sets of parameters: (from left to right) vortex core radii 3.43×10−4, 1.715×10−4and 8.575×10−5m;
and (from top to bottom) vortex strengths 6.71×10−2, 3.355 ×10−2, 1.6775×10−2 and 8.3875×10−3m2/s.
time of the scheme. Based on these results, the temperature-scaled central moments MRT collision operator is the most efficient way to take into account the full Hermite expansion in EDF and extend the stability domain of the collision operator as it only introduces a factor of 2.5 in the computation time. The low Mach formulation on the other hand, not needing the third-order components of the EDF and relying on a SRT collision operator is observed to match the performances of the bare second-order SRT solver and outperform the com-pressible solver by a factor of 2.7. Throughout the different simulations, it was also observed that using finite differences to solve additional fields in 3-D, i.e. energy and species mass balance, would result in a performance of 5.59 MLUPS per field while the corresponding LB solver would peak at 5.1 MLUPS for the smallest possible stencil, i.e. D3Q7. It must be noted that time spent on evaluating thermo-physical parameters and source terms is not taken into account here. The coupling to REGATH and computation of thermo-physical pa-rameters and source terms, assuming a flow with three species and a single reaction accounts
SRT-EQ2SRT-EQ3SRT-EQ6RR-EQ3RR-EQ6CM-MRT 0
0.1 0.2 0.3 0.4 0.5
MLUPS
SRT-EQ2SRT-EQ3low Mach 0
0.1 0.2 0.3 0.4 0.5 0.6
Figure 6.21: Performances of (left) the different collision models for the compressible formula-tion and (right) comparison of compressible solver to low Mach formulaformula-tion. EQN designates the order of the EDF, CM-MRT refers to the MRT collision operator in temperature-scaled central moment space, RR to the recursive regularized collision operator and low Mach to the low Mach dilatable solver with a second-order EDF, in terms of MLUPS, short for Million Lattice Updates Per Second.
on the evaluation of these parameters using REGATH might be explained, at least in part, by the fact that ALBORZ is a C++ code while REGATH is written in Fortran. Given the different approaches to memory allocation for multi-dimensional arrays in each one of these programming languages, non-negligible resources and time are spent on converting data ar-rays between them. As such, elaboration of a more efficient communication scheme between the two solvers is expected to have non-negligible impact on the performances of the overall code.
In conclusion, while all three approaches considered in this section are able to model low Mach number combustion, the low Mach formulation of section 5.1 coupled to FD solvers (approach I) for the additional fields allows both for larger time-steps and better perfor-mances. Furthermore, the use of FD solvers for energy and species balance equations allows for more flexibility in terms of the order of discretization in both space and time. It must also be noted that while FD solvers based on a first-order time-discretization were slightly faster than their LB counterparts, the latter were second-order accurate under both acoustic and diffusive scaling. To achieve a similar order of accuracy (for example using a second-order RK scheme) the cost would nearly double. The compressible solver on the other hand is a promising formulation for fully compressible flows, and can be further improved through concepts such as flow state-adaptive discretization of the distribution function. Some of these concepts and possible directions are reviewed and studied in the next chapter.
Shifted stencils: A promising
approach to extend the stability domain of the lattice Boltzmann method
Contents
7.1 Adaptive phase-space discretization of distribution function:
brief overview . . . 153 7.2 Galilean invariance of Gauss-Hermite quadratures . . . 154 7.3 Shifted stencils with on-lattice streaming . . . 155 7.4 Final remarks . . . 162
[1] S.A. Hosseini, C. Coreixas, N. Darabiha and D. Th´evenin. Extensive analysis of the lattice Boltzmann method on shifted stencils. Physical Review E, 100(6), p.063301, 2019
7.1 Adaptive phase-space discretization of distribution function: brief overview
Throughout the different sections of the present thesis, it was shown and observed through both theoretical analysis (asymptotic and spectral) and numerical applications that the non-dimensionalization strategy and parameters, also known as reference state, have a huge impact on the performances of the LB-based solver. They impact both stability and accu-racy of the scheme. The error terms appearing at higher-order levels of the CE expansion, of no physical interest, usually tied to higher-order (not supported by the stencil) moments of the EDF are controlled by the scaling of the conserved variables. In chapters 2 and 5 it was clearly shown that deviations of the local flow state from reference temperature θ = 1 tied to δ2x/δt2 and velocity u = 0 resulted in reduced stability domains and larger errors coming from higher-order moments. The same behavior was also observed for advection-diffusion based solvers in chapter 4. While one would expect the use of higher-order quadratures to help broaden the stability domain and reduce errors, their effect remains limited. Further-more, these quadratures come at much larger computational costs (especially in terms of memory consumption and memory access bandwidth). A new concept, first introduced in [224, 225, 226] and later improved in [95], extended the applicability domain of the LBM by using adaptive discrete velocity stencils. While new in the LBM community, it is interesting to note that this approach to extending discrete Boltzmann solvers to high Mach number applications had already been proposed in [11, 12, 13] in the context of kinetic solvers with Eulerian discretization of physical space. In the model proposed in [11, 12, 13] the origin of the discrete velocity set was determined by the local fluid velocity while its unit was tied to the local internal energy (tied to local temperature). This approach is illustrated in Fig. 7.1.
In this figure, the blue dashed vectors represent the discrete velocities used for phase-space discretization. These discrete velocities are composed of the local fluid velocity (red vector) and peculiar velocity components (black vectors), as appearing in the Maxwell-Boltzmann EDF. The peculiar components, are scaled with the local internal energy. The final discrete velocities are therefore direct results of the local fluid state. Another interesting point noted in [12] is that for a stencil based on first neighbors, in the limit of vanishing viscosities, this approach leads to the Euler equations. The benefits of such an adaptive construction of the discrete velocities are limited in the context of an Eulerian discretization of space, as the solver would still be limited by the CFL condition. The LB formulation however, based on an exact integration of the convection term in the Boltzmannn equation could benefit much more from an adaptive reconstruction of the discrete populations, as it would guarantee that numerical scheme characteristics embed the analytical (from the continuous PDE) domain of dependence. That is why Sun et al. were able to model high Mach number flows using this approach and the LBM formalism in [224, 225, 226]. Adaptive discrete velocities were recently updated (with a predictor-corrector formulation) and reformulated in a more con-sistent theoretical framework in [95] as the PonD (Particles on Demand) formulation. In the PonD method, the populations are not guaranteed to fall on a discrete grid-point after the streaming step, as the discrete velocities have a continuous component stemming from
reconstruction step (of populations on the discrete grid-points) using Lagrange polynomials.
An on-grid version of the method (minus the adaptive nature of the velocity shift) was also used in [94]. To clarify the effect of adaptive stencils and their potential in modeling flows
Figure 7.1: The concept of adaptive stencils as proposed in [11, 12, 13].
involving large temperature variations and Mach number, as a proof of concept, we will present a brief numerical and theoretical analysis of the LBM on shifted stencils.
7.2 Galilean invariance of Gauss-Hermite quadratures
To demonstrate the Galilean invariance of the Gauss-Hermite quadrature, one can start with the Maxwell-Boltzmann EDF and write it in a non-dimensional form as:
f(eq)(ξ∗,u∗, θ) = ρ(2πθ)−D/2exp−(ξ∗ −u∗)2
2θ∗ , (7.1)
whereDis the physical dimension of the system, θ∗ = kBcT /m
s whileξ∗ = ξ−Uc
s andu∗ = u−Uc
s . With this change of variable the reference velocity changes from zero toU, while the reference sound speed is set tocs. Proceeding as for the classical stencils to expand the EDF in terms of Hermite polynomials [90, 227, 88]:
f(eq)(ξ∗,u∗, θ∗) = w(ξ∗)
∞
X
n=0
1
n!a(eq)n (u∗, θ∗) :Hn(ξ∗), (7.2) where the weight function w(ξ∗) is defined as:
w(ξ∗) = (2π)−D/2exp
−ξ∗2 2
, (7.3)
one obtains the first few Hermite polynomials as a function of the re-defined velocities:
H0 = 1, (7.4a)
H1,i1 =ξi∗1, (7.4b)
H2,i1i2 =ξi∗
1ξi∗
2 −δi1i2, (7.4c)
H3,i1i2i3 =ξi∗1ξi∗2ξi∗3−
ξ∗i1δi2i3
cyc, (7.4d)
H4,i1i2i3i4 =ξi∗1ξi∗2ξi∗3ξi∗4 + [δi1i2δi3i4]cyc−
ξi∗1ξi∗2δi3i4
cyc (7.4e)
It is interesting to note that the forms and expression of the polynomials are exactly similar to those in Eq. 2.44. Using these polynomials and the EDF, the corresponding coefficients can be computed:
a(eq)0 =ρ, (7.5a)
a(eq)i
1 =ρu∗i
1, (7.5b)
a(eq)i1i2 =ρu∗i1u∗i2 +ρ(θ∗−1)δi1i2, (7.5c) a(eq)i1i2i3 =ρu∗i1u∗i2u∗i3 +ρ(θ∗−1) [ui1δi2i3]cyc, (7.5d) a(eq)i
1i2i3i4 =ρu∗i
1u∗i
2u∗i
3u∗i
4 +ρ(θ∗−1)2[δi1i2δi3i4]cyc+ρ(θ∗−1) u∗i
1u∗i
2δi3i4
cyc. (7.5e) As for the Hermite polynomials it is readily observed that, setting U = 0 and cs to pkBT0/m0, one would recover the classical LBM formulation Hermite coefficients. Discretiz-ing phase-space usDiscretiz-ing the Gauss-Hermite quadrature [227, 88], at order three for example, one obtains exactly the same non-dimensional abscissae c∗α with (in the case of the D2Q9 stencil) c∗α,i ∈ {−√
3,0,√
3}, and corresponding weights, i.e. wα = Πi=x,ywα,i. This clearly establishes the Galilean invariance of the phase-space discretization procedure. To better convey the possibilities of such an approach let us first write down the resulting discrete EDF (expanded up to order two for simplicity) in terms of flow field variables:
fα(eq,2) =wαρ
"
1 + (cα−U)·(u−U)
c2s +[(cα−U)·(u−U)]2
2c4s − (u−U)2 2c2s +
(cα−U)2−D
(θ∗−1) 2c4s
#
. (7.6) SettingU = 0 andcs= δx
δt
√
3 one recovers the classical LBM formulation.