Next Article in Journal
Validation of the LOGOS Software Package Methods for the Numerical Simulation of Cavitational Flows
Previous Article in Journal
Flow Boiling Heat Transfer Performance and Boiling Phenomena on Various Straight Fin Configurations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Analysis of Hierarchical Grid Refinement Techniques for the Lattice Boltzmann Method by Numerical Experiments

1
Department of Mechanical and Process Engineering, University of Kaiserslautern-Landau, 67663 Kaiserslautern, Germany
2
Department of Mechanical and Process Engineering, Offenburg University, 77652 Offenburg, Germany
*
Author to whom correspondence should be addressed.
Fluids 2023, 8(3), 103; https://doi.org/10.3390/fluids8030103
Submission received: 21 February 2023 / Revised: 7 March 2023 / Accepted: 16 March 2023 / Published: 21 March 2023
(This article belongs to the Topic Computational Fluid Dynamics (CFD) and Its Applications)

Abstract

:
Over the last few decades, several grid coupling techniques for hierarchically refined Cartesian grids have been developed to provide the possibility of varying mesh resolution in lattice Boltzmann methods. The proposed schemes can be roughly categorized based on the individual grid transition interface layout they are adapted to, namely cell-vertex or cell-centered approaches, as well as a combination of both. It stands to reason that the specific properties of each of these grid-coupling algorithms influence the stability and accuracy of the numerical scheme. Consequently, this naturally leads to a curiosity regarding the extent to which this is the case. The present study compares three established grid-coupling techniques regarding their stability ranges by conducting a series of numerical experiments for a square duct flow, including various collision models. Furthermore the hybrid-recursive regularized collision model, originally introduced for cell-vertex algorithms with co-located coarse and fine grid nodes, has been adapted to cell-centered and combined methods.

1. Introduction

Due to the associated significant increase in computational efficiency, local grid refinement has become an indispensable part of modern computational fluid dynamics methods (CFD), allowing for varying spatial resolution throughout the simulation domain. Regions with, e.g., strong gradients of the flow field or flow structures of interest are thereby discretized with a larger number of grid points, whereas a coarser grid size may be used where rather moderate flow field variations are present.
With the lattice Boltzmann Method (LBM) gaining evermore popularity among CFD practitioners over the last few decades, several grid refinement techniques tailored to the method have been developed in order to satisfy the aforementioned need for efficient usage of computational resources. A widely recognized category of such techniques takes advantage of the fact that the LBM was originally developed for Cartesian grids by hierarchically subdividing the grid cells in a quad-(2D) or octree (3D)-like manner, maintaining length and time scale ratios of two between neighboring grid levels.
A further conceptual categorization of methods belonging to this so-called multi-level approach can be made according to the specific arrangement of coarse and fine grid nodes relative to one another as depicted schematically in Figure 1 for two dimensions. Whereas cell-vertex algorithms are characterized by grid nodes residing in cell corners and thus partially co-located coarse and fine nodes along grid transition interfaces (see Figure 1a) [1,2,3], cell-centered approaches apply a volumetric description without the possibility of co-location (see Figure 1b) [4,5,6]. Another type of somewhat mixed arrangement arises when both grids are shifted relative to each other with coarse nodes sitting in the coarse cell corners and fine nodes residing in their respective cell centers (see Figure 1c) [7,8]. These three grid layout types will herein be referred to as cell-vertex, cell-centered and combined approaches, and involve individual coupling procedures to enable necessary information exchange between grid levels. Details on the algorithmic peculiarities of each coupling mechanism will be given in Section 3.
It is well known that the LBM, as an explicit finite difference method, has certain stability limits, which was shown in [9] by means of linear stability analysis. This instability issue is usually amplified by the use of the multi-level approach due to the introduction of disturbances at grid transition interfaces [10]. These disturbances may manifest in the form of spurious acoustic and vorticity waves, corrupting the solution accuracy and stability of multi-level LBM [11], which is of utmost importance for the realization of reliable high Reynolds number simulations.
It seems reasonable to assume the extent of occurring disturbances and thus the stability limit and accuracy of a multi-level scheme to be influenced by the specific coupling mechanisms of each grid arrangement type. While in [12] it is mentioned that the cell-vertex layout qualitatively emits less spurious acoustic noise compared to the cell-centered layout, to the best of the authors’ knowledge no further quantitative stability and accuracy comparisons between the different types are present in the literature. The subject of this paper is therefore an investigation of the aforementioned three grid layouts in terms of their stability and accuracy through a series of numerical experiments on a D3Q19 lattice for a force-driven flow through a square duct. Due to their significant impact on the stability of the LBM, the simulations were carried out employing four different collision operators:
  • Bhatnagar–Gross–Krook collision model (BGK) [13]
  • Multiple relaxation time collision model (MRT) based on the work of d’Humiéres et al. [14]
  • Recursive regularized collision model (RR) as proposed by Malaspinas [15]
  • Hybrid-recursive regularized collision model (HRR) as first described by Jacob et al. [16].
The HRR collision operator was introduced for and applied to cell-vertex grids and originally involves finite difference approximations of strain-rates in order to hybridize the calculation of the stress tensor needed for the reconstruction of second-order non-equilibrium moments. Due to a lack of neighbors on the same level, this requires some adjustments to the multi-level collide-and-stream algorithm in grid transition interfaces if homogeneous use of second-order central differences throughout the whole calculation domain is to be made. Astoul et al. [12] relied on a fictitious transport step to synchronize in time the interface nodes participating in the central finite difference computation in this case. In our method, no such fictitious transport step is necessary to allow a continuous calculation using central differences. Furthermore, the HRR algorithm was adapted in a consistent way to cell-centered and combined grid layouts in the course of this work.
The paper is organized as follows. After an introduction of the theoretical basics regarding the applied collision operators, physical models and boundary conditions in Section 2, a detailed review of the algorithmic procedures of each of the investigated grid layouts, extended by the specifics of the HRR algorithm with regards to the cell-vertex grid layout follows in Section 3. Subsequently, the adaptation of the HRR algorithm to the cell-centered and combined grid layouts is addressed and a consistent solution is presented in Section 4. Finally, the square duct test case and the investigated variants of the numerical experiments are presented and discussed together with the results in Section 5.

2. Lattice Boltzmann Method

The lattice Boltzmann method, which historically evolved from lattice gas automata [17], but mathematically can be derived directly from the Boltzmann equation, represents an explicit system of Equation (1) for the spatiotemporal evolution of a discrete set of distribution functions f i describing phase-space densities of fictitious particles. The equations are typically solved on a Cartesian grid, with a velocity space discretization such that the fictitious particles travel with a velocity ξ i from one grid node to selected neighboring nodes during one time step Δ t . The left-hand side of the lattice Boltzmann Equation (LBE) (1) reflects this node-to-node particle transport, while the right-hand side describes local particle redistribution by scattering processes through the collision operator Ω i x , t :
f i x + ξ i Δ t , t + Δ t f i x , t = Ω i x , t .
Following Qian [18], the discretized phase-spaces or lattices are labeled according to a DdQq nomenclature where d represents the number of spatial dimensions and q reflects the finite number of discrete particle velocities that result from the restriction of velocity space. Figure 2 shows two common velocity sets, namely D2Q9 and D3Q19, the latter being the basis for all numerical experiments discussed in Section 5, whereas the former is used for its simplicity in order to review and clarify concepts regarding the different grid layouts in Section 3. With the molecular velocity ξ ˜ = Δ x Δ t , particle velocity vectors for the D3Q19 lattice are given by
ξ i = 0 , 0 , 0 i = 18 ± ξ ˜ , 0 , 0 , 0 , ± ξ ˜ , 0 , 0 , 0 , ± ξ ˜ i = 0 , , 5 ± ξ ˜ , ± ξ ˜ , 0 , ± ξ ˜ , 0 , ± ξ ˜ , 0 , ± ξ ˜ , ± ξ ˜ i = 6 , , 17 .
The connection between the molecular velocity ξ ˜ , grid spacing Δ x and time-step size Δ t reflects the aforementioned particle transport from one node to its neighbors along discrete velocity directions during each time step.
In the continuous picture, macroscopic physical quantities can be extracted from the information contained within the mesoscopic scale as so-called moments of the velocity distribution function f x , ξ , t . Defined as integrals of f weighted with some function of ξ over infinite velocity space, these moments are evaluated as weighted sums using Gauss–Hermite quadrature with discrete velocity vectors as abscissae in discretized phase space [19]:
M x , t = φ ξ f x , ξ , t d ξ = i = 0 q 1 φ ξ i f i
f i = w i f x , ξ i , t
The quadrature weights w i in Equation (4) emerge from specific conditions imposed upon the equivalence between the moments of the discrete distribution functions f i and the ones based on its continuous counterpart, as expressed in Equation (3) [20]. The quadrature order, i.e., the degree of precision, is defined by its ability to exactly abide Equation (3) up to a certain order of moments [21]. From this, the weight factors for the D3Q19 lattice result to
w i = 1 3 , i = 18 1 18 , i = 0 , , 5 1 36 , i = 6 , , 17
and the relationship between the molecular velocity ξ ˜ and the isothermal speed of sound c s = p / ρ = R T , also referred to as lattice constant, is then given by ξ ˜ = 3 c s . The quadrature order of the D3Q19 lattice leads to a O Ma 3 error and violation of Galilean invariance in the viscous stress tensor [22], which is usually negligible for small Mach numbers or in the incompressible limit.
Regarding the right-hand side of Equation (1), a large number of different collision operators have been developed [13,14,15,16,23,24,25], some of which are used in the present work and will now be briefly discussed.

2.1. Bhatnagar–Gross–Krook Collision Model

Given that the difficulties in solving the Boltzmann equation were caused mainly by the complexity of Boltzmann’s original collision integral, Bhatnagar, Gross and Krook [13] proposed a considerable simplification by approximating the collision process through a relaxation of the velocity distribution function f x , t towards the Maxwell–Boltzmann equilibrium distribution f e q ρ , u , ξ using a single relaxation rate 1 τ , with τ being the so-called collision time:
Ω BGK = 1 τ f x , t f e q ρ , u , ξ .
Since all involved physical processes are equilibrated with the same relaxation rate, the classical BGK approximation is also termed single relaxation time model (SRT) [20]. The equilibrium distribution function as a function of fluid density ρ , local mean or macroscopic velocity u and particle velocity ξ reads in D spatial dimensions [26]
f e q ρ , u , ξ = ρ 2 π c s 2 D 2 exp ξ u 2 2 c s 2 ,
where the temperature dependency is dropped, indicated by the isothermal speed of sound c s , since only isothermal cases are subject of the present work.
Adapting the BGK approximation to the lattice Boltzmann Equation (1), the discrete counterpart to Ω BGK in Equation (6) becomes
Ω i BGK = Δ t τ + 1 2 Δ t f i f i ( 0 ) = c s 2 Δ t ν + 1 2 c s 2 Δ t f i f i ( 0 ) = ω ˜ f i n e q .
Here, the continuous velocity distribution function f and the Maxwell–Boltzmann equilibrium f e q are replaced by their respective discrete equivalent f i and f i ( 0 ) . Compared to Equation (8), the term 1 2 Δ t is added to the relaxation factor combined within dimensionless collision frequency ω ˜ . This additional factor is also known as discrete lattice effect [27] and represents numerical viscosity, absorbed into the physical model to produce correct macroscopic transport coefficients [28].
In order to fulfill the macroscopic conservation laws, it is not necessary to consider the complete mesoscopic Maxwell–Boltzmann distribution function. Utilizing similarities between the generating function of Hermite polynomials ω (see definition in Equation (10)) and the Maxwell–Boltzmann equilibrium f e q a simplified equilibrium in terms of a truncated series expansion in Hermite polynomials is thus used [29], allowing for a reduction in numerical effort [26].
Although the first three terms of the Hermite expansion are sufficient for Navier–Stokes level hydrodynamics [26], some authors claim positive effects on numerical stability and accuracy by including higher-order expansion terms [14]. In order to remove spurious couplings among third-order Hermite moments [30], we rely on orthogonal third-order Hermite polynomials in the expansion of f i ( 0 ) as in [16]. Finally, the discrete equilibrium functions used for all single relaxation time models in the present paper are given by
f i ( 0 ) = w i [ n = 0 2 1 c s 2 n n ! H i ( n ) : A 0 ( n ) + 1 2 c s 6 ( H i , x x y ( 3 ) + H i , y z z ( 3 ) ) A 0 , x x y ( 3 ) + A 0 , y z z ( 3 ) + 1 2 c s 6 ( H i , x x z ( 3 ) + H i , y y z ( 3 ) ) A 0 , x x z ( 3 ) + A 0 , y y z ( 3 ) + 1 2 c s 6 ( H i , x y y ( 3 ) + H i , x z z ( 3 ) ) A 0 , x y y ( 3 ) + A 0 , x z z ( 3 ) + 1 6 c s 6 ( H i , x x y ( 3 ) H i , y z z ( 3 ) ) A 0 , x x y ( 3 ) A 0 , y z z ( 3 ) + 1 6 c s 6 ( H i , x x z ( 3 ) H i , y y z ( 3 ) ) A 0 , x x z ( 3 ) A 0 , y y z ( 3 ) + 1 6 c s 6 ( H i , x y y ( 3 ) H i , x z z ( 3 ) ) A 0 , x y y ( 3 ) A 0 , x z z ( 3 ) ] .
Due to its insufficient quadrature order, isotropic Hermite tensors of the form H i , α α α ( 3 ) or H i , α β γ ( 3 ) have to be excluded for the D3Q19 lattice [31,32].
Hermite equilibrium expansion coefficients A 0 ( n ) in Equation (9) are obtained through a projection of the Maxwell–Boltzmann equilibrium distribution onto a n th -order Hermite polynomial basis in D = 3 spatial dimensions, defined as [12]:
H i ( n ) = c s 2 n ω ξ i ξ n ω ξ i with ω ξ = 1 2 π c s 2 D 2 exp | | ξ ˜ | | 2 2 2 c s 2 .
These expansion coefficients are related to or directly coincide with equilibrium moments and are given by [32]
A 0 ( 0 ) = ρ , A 0 , α ( 1 ) = ρ u α , A 0 , α β ( 2 ) = ρ u α u β + ρ c s 2 δ α β , A 0 , α β γ ( 3 ) = ρ u α u β u γ δ α β where δ α β = 1 , α = β 0 , α β
It should be mentioned that the equilibrium distribution function is sometimes expanded as a series in macroscopic velocity u rather than Hermite polynomials. This as Mach number expansion known procedure [19] yields identical results up to second-order, but differs from the Hermite approach for higher-order terms due to emerging spurious couplings between lower- and higher-order moments [26]. For this reason, the Hermite series expansion is preferred.

2.2. Multiple Relaxation Time Collision Model

Although a significant simplification of the collision operator is achieved by applying the single relaxation time approach described in the previous section, its use corresponds to a relaxation of all relevant physical processes at one and the same rate, namely the collision frequency ω ˜ , towards their respective equilibrium state. In reality, however, these processes unfold on different time scales, yielding various macroscopic transport coefficients. The idea of the multiple relaxation time collision model is therefore to perform collision in moment space, rather than velocity space, individually equilibrating moments of the distribution functions.
For the D3Q19 lattice, multiple relaxation time collision was first proposed by d’Humières et al. [14]. The MRT collision operator to be inserted into the LBE (1), which is then sometimes referred to as generalized LBE [33], reads as
Ω i MRT = j = 0 q 1 M 1 S ^ M i j f j 0 f j .
Here, distribution functions are first mapped onto discrete moment space via the linear transformation matrix M . In case of the D3Q19 lattice, the diagonal matrix S ^ contains 19 relaxation parameters S ^ = diag s 0 , s 1 , , s 18 by means of which moment equilibration is performed. Retransformation back to velocity space is achieved by finally applying the inverse matrix M 1 to the relaxed moments. The transformation matrix M holds 19 orthogonal base vectors φ i , found through Gram–Schmidt orthogonalization procedure, as rows [14], see Equation (A1). Through scalar multiplication with the vector of distribution functions f , 19 linear independent moments are constructed [20]
(13) m = φ 0 · f , φ 1 · f , , φ 18 · f T = M f m = ρ , e , ε , j x , q x , j y , q y , j z , q z , 3 p x x , 3 π x x , p w w , π x x , p x y , p y z , p x z , m x , m y , m z T (14) = m 0 , m 1 , m 2 , m 3 , m 4 , m 5 , m 6 , m 7 , m 8 , m 9 , m 10 , m 11 , m 12 , m 13 , m 14 , m 15 , m 16 , m 17 , m 18 T .
For the most part, the resulting moments in Equation (14) allow a physical interpretation and are categorized based on whether they represent conserved physical quantities, also referred to as hydrodynamic moments or non-conserved physical quantities called kinetic moments [14]. Higher-order moments devoid of any physical interpretation are also termed ghost moments [34].
As mentioned above, the moments contained within the moment vector m are relaxed to their corresponding equilibrium state m ( 0 ) = M f ( 0 ) , where f ( 0 ) holds equilibrium distributions f i ( 0 ) expanded up to second order, that means identical to Equation (9) with only the first term inside the square brackets remaining. Since f ( 0 ) only depends on ρ and u , all equilibrium moments are also describable solely in terms of these quantities, see Equation (A2).
Regarding individual relaxation parameters within S ^ , all parameters except those related to kinematic viscosity ν through the third equal sign in Equation (8), i.e., s i = ω ˜ , i 9 , 11 , 13 , 14 , 15 , can be freely chosen from the interval 0 , 2 . Furthermore, since m i = m i ( 0 ) , i 0 , 3 , 5 , 7 , relaxation parameters for collision invariants fluid density ρ and momentum j α = ρ u α are set to 0 for simplicity. The optimal set of arbitrary parameters in terms of stability can be determined for a particular case through a von Neumann analysis [35], resulting in the following choice for a channel flow: s 1 = 1.19 , s 2 = s 10 = s 12 = 1.4 , s 4 = s 6 = s 8 = 1.2 , s 16 = s 17 = s 18 = 1.98 .
Another possibility, commonly known as regularized scheme [36], would be to specify free relaxation factors as s i = 1 , i 1 , 2 , 4 , 6 , 8 , 10 , 12 , 16 , 17 , 18 , which corresponds to a direct relaxation of the associated moments towards their equilibrium state. As this also affects higher-order non-hydrodynamic moments by effectively dampening their respective modes, stability is positively influenced [37]. These two specific choices of free parameter sets, in this paper refered to as MRT-NEU and MRT-REG, respectively, are used for all simulations with the multiple relaxation time model in Section 5. An alternative formulation of the regularized approach in velocity space, as well as various enhancements will be discussed in more detail in the next section.

2.3. Regularized Collision Models

From a mathematical perspective, the bridge between the mesoscopic scale of the particle distribution function and macroscopic Naiver-Stokes level hydrodynamics is built using Chapman-Enskog analysis [38]. In case of the lattice Boltzmann equation, discrete distribution functions f i are expanded as a perturbation series in Knudsen number around their local equilibrium state f i ( 0 ) :
f i = f i ( 0 ) + ϵ f i ( 1 ) + ϵ 2 f i ( 2 ) + ϵ 3 f i ( 3 ) + = f i ( 0 ) + f i n e q ,
where all deviations from equilibrium are gathered in the non-equilibrium distributions f i n e q . The smallness parameter ϵ expresses each terms relative order in Knudsen number. By Taylor expanding the LBE (1) and applying a similar multiple scale expansion to the emerging time derivatives together with the ansatz that only the two lowest orders in Knudsen number are necessary for an accurate description of macroscopic fluid behavior, the Navier–Stokes equations are recovered [26].
The core idea behind regularized collision operators is to abandon contributions of higher-order terms O ϵ n , n > 1 in (15) and approximate the complete non-equilibrium distribution functions by
f i n e q = f i f i ( 0 ) f i ( 1 ) ,
with an explicit reconstruction of the first-order non-equilibrium part f i ( 1 ) . Applying this regularization procedure to the pre-collision distribution functions together with a single relaxation time collision operator [39,40], results in the following form of the LBE:
f i x + ξ i Δ t , t + Δ t f i REG = Ω i REG (17) f i REG = f i 0 + f i 1 (18) Ω i REG = ω ˜ f i 1 (19) f i x + ξ i Δ t , t + Δ t = f i 0 + 1 ω ˜ f i 1
Similar to the formulation of the discrete equilibrium functions f i ( 0 ) in Equation (9), the regularized first-order non-equilibrium part f i ( 1 ) is reconstructed through Hermite series expansion
f i ( 1 ) = w i n = 1 2 1 c s 2 n n ! H i ( n ) : A 1 ( n ) .
Non-equilibrium expansion coefficients A 1 ( n ) are obtained by projection of the non-equilibrium functions onto the nth-order Hermite tensor [31]
A 1 , α ( 1 ) = i H i , α ( 1 ) f i f i ( 0 ) = i ξ i , α f i n e q
A 1 , α β ( 2 ) = i H i , α β ( 2 ) f i f i ( 0 ) = i ξ i , α ξ i , β δ α β c s 2 f i n e q ,
Notice that the Hermite series expansion of f i ( 1 ) is often presented as beginning from second-order [12,31,40,41], since A 1 ( 0 ) = 0 and A 1 ( 1 ) = 0 by definition. However, Li et al. [42] pointed out that first-order terms must be included in the reconstruction, if semi-implicit forcing schemes with so-called half-force correction, as proposed by Guo et al. [27,43], are deployed, because then only the zeroth-order coefficient vanishes, while A 1 ( 1 ) 0 . The importance of this was demonstrated by means of a body force-driven 2D Poiseuille flow, showing significant differences between simulation and analytical solution. In our study, such forcing scheme is applied to the force-driven square duct flow in Section 5 and thus first-order terms are taken into account. The treatment of forces will be covered in more detail in Section 2.4.
As mentioned in [30], this method can be reinterpreted in terms of a MRT collision model with discrete moment space built upon a Hermite tensor product polynomial basis with non-hydrodynamic moments directly equilibrated by setting their corresponding relaxation parameters to 1.
Malaspinas proposed in [15] to further improve the stability of regularized schemes by including higher-order Hermite expansion terms in the reconstruction of f i ( 1 ) . Using Chapman–Enskog analysis, it was proven that non-equilibrium moments of order higher than two can be obtained recursively from lower order moments, which allows to reconstruct f i ( 1 ) up to any order given the second order expansion coefficients and the macroscopic velocity u . Analogous to the equilibrium functions in Equation (9), we rely on orthogonal third-order non-equilibrium moments in order to deal with spurious couplings among them [16,30], resulting in
f i ( 1 ) = w i [ n = 1 2 1 c s 2 n n ! H i ( n ) : A 1 ( n ) + 1 2 c s 6 ( H i , x x y ( 3 ) + H i , y z z ( 3 ) ) A 1 , x x y ( 3 ) + A 1 , y z z ( 3 ) + 1 2 c s 6 ( H i , x x z ( 3 ) + H i , y y z ( 3 ) ) A 1 , x x z ( 3 ) + A 1 , y y z ( 3 ) + 1 2 c s 6 ( H i , x y y ( 3 ) + H i , x z z ( 3 ) ) A 1 , x y y ( 3 ) + A 1 , x z z ( 3 ) + 1 6 c s 6 ( H i , x x y ( 3 ) H i , y z z ( 3 ) ) A 1 , x x y ( 3 ) A 1 , y z z ( 3 ) + 1 6 c s 6 ( H i , x x z ( 3 ) H i , y y z ( 3 ) ) A 1 , x x z ( 3 ) A 1 , y y z ( 3 ) + 1 6 c s 6 ( H i , x y y ( 3 ) H i , x z z ( 3 ) ) A 1 , x y y ( 3 ) A 1 , x z z ( 3 ) ] .
Again, isotropic tensors H i , α α α ( 3 ) and H i , α β γ ( 3 ) must be ruled out for D3Q19. Due to aforementioned recurrence properties, this method is known as recursive-regularization (RR). Third-order moments and Hermite tensors in Equation (23) are given by
A 1 , α α β ( 3 ) = 2 u α A 1 , α β ( 2 ) + u β A 1 , α α ( 2 ) ,
H i , α α β ( 3 ) = ξ i , α ξ i , α c s 2 ξ i , β
A further improvement regarding stability was achieved with the hybrid-recursive regularized collision operator (HRR), introduced by Jacob et al. [16]. In this method the calculation of second-order non-equilibrium expansion coefficients A 1 ( 2 ) (see Equation (22)), which relate to the viscous stress tensor through [31]
A 1 , α β ( 2 ) = 2 ρ c s 2 Δ t ω ˜ S α β ,
is hybridized before being utilized in the reconstruction of f i ( 1 ) in Equation (23):
A 1 ( 2 ) = σ A 1 ( 2 ) , PR + 1 σ A 1 ( 2 ) , FD where [ 0 σ 1 ] .
A 1 ( 2 ) , PR represents projection-based moments computed by Equation (22), whereas the second term A 1 ( 2 ) , FD is determined by approximating strain-rates S α β = 1 2 β u α + α u β contained within the stress tensor using central finite differences, yielding:
A 1 , α β ( 2 ) , FD = ρ c s 2 Δ t ω ˜ u α x + e β Δ x u α x e β Δ x 2 Δ x + u β x + e α Δ x u β x e α Δ x 2 Δ x .
The recursive-regularized model is recovered by setting the hybridization factor σ in Equation (27) to 1, whereof σ = 0 leads to a complete finite difference reconstruction of the viscous stress tensor. By choosing σ 0 , 1 , compared to the RR model, a significant stability increase is achieved by injection of numerical hyper-viscosity due to the usage of finite differences, strongly damping non-hydrodynamic mode contributions [16,44,45], which negatively affect stability especially for non-uniform grids [12,41].
It was noted before that the insufficient quadrature order of the D3Q19 lattice leads to a viscous stress tensor exhibiting an O Ma 3 error term and suffering from a violation of Galilean invariance. To deal with this issue, Feng et al. [46] proposed a cubic Mach error correction term for thermal LBM, allowing for a great improvement of the HRR collision model in the high-subsonic regime [47]. Details on the calculation of the cubic Mach correction term are summarized in Appendix B. The regularized LBE including the correction term ψ i becomes:
f i x + ξ i Δ t , t + Δ t = f i 0 + 1 ω ˜ f i 1 + Δ t 2 ψ i .
Moreover, the projection-based second-order non-equilibrium moments in Equation (27) need to be modified accordingly
A 1 ( 2 ) , PR = i H i ( 2 ) f i f i ( 0 ) + Δ t 2 ψ i .
Even though the numerical experiments presented in Section 5 were all conducted for a comparatively small Mach number, HRR simulations have been exemplarily performed and evaluated, both, with and without cubic Mach correction. Corrected HRR collision will be referred to as HRR ψ .

2.4. Treatment of Forces

The test case investigated in terms of its numerical accuracy and stability range in Section 5, consists of a body force-driven square duct flow. Since the simulations have been carried out for various collision operators, herein, the adopted force model is the scheme proposed by Guo, adjusted to the respective collision operator type in use. First, the original version [27] for the standard BGK model will be addressed.
The presence of forces in the LBE is taken into account by an additional source term [27]. The right-hand side of the LBE (1) with standard BGK collision operator then becomes:
Ω i BGK = ω ˜ f i f i 0 + Δ t F i
Following the categorization in [42], Guo’s forcing is classified as semi-implicit scheme, with the collision and force terms both being integrated in time by second-order accurate trapezoidal rule [48] during spatiotemporal discretization of the discrete-velocity Boltzmann equation. The forcing term F i is written as part of the collision operator and derived by means of Hermite expansion:
F i = w i ω ˜ ν c s 2 Δ t ξ i v c s 2 + ξ i · v c s 4 ξ i · f .
Here f = ρ a symbolizes the force per unit volume, with a being the acceleration. Notice also the influence of the discrete lattice effect on the forcing term. An important implication of this model is the so-called half-force correction of fluid momentum
ρ v = i = 0 q 1 ξ i f i + 1 2 f Δ t ,
expressing its average value within one discrete time step Δ t due to the impact of force f , yielding Guo’s forcing scheme more robust [31] and accurate compared to other existing schemes [42]. The corrected velocity v is applied to the reconstruction of the equilibrium functions f i ( 0 ) in Equation (9) and the forcing term F i in Equation (32).
Guo’s forcing scheme adapted to the MRT approach presented in Section 2.2 gives the following expression for Ω i MRT including the forcing term F ^ i [49]:
Ω i MRT = j = 0 q 1 M 1 S ^ M i j f j 0 f j + Δ t F ^ i
Since the discrete lattice effect in Equation (32) depends on the relaxation parameter ω ˜ , it is obvious that a mapping to moment space has to take place in the MRT adaption of the forcing scheme. F ^ i thus represents the i th component of the source term vector f ^ I R 19 [20]:
f ^ = M 1 I 1 2 S ^ M f ¯ ,
with vector f ^ consisting of discrete force components F ¯ i excluding the discrete lattice effect:
F ¯ i = w i ξ i v c s 2 + ξ i · v c s 4 ξ i · f .
Regarding the recursive-regularized and hybrid recursive-regularized collision models including Guo’s forcing scheme, we implemented the procedure proposed by Yoo et al. [31]. The pre-collision distribution functions in regularized models with discrete forcing should be defined to not include the half-force corrected velocity, similar to the standard BGK forcing scheme. Following this line of thought, the final form of the regularized LBE with forcing term and cubic Mach correction results to
(37) f i REG = f i 0 v Δ t 2 f + f i 1 = f i 0 + f i 1 Δ t 2 F i ψ i (38) Ω i REG = ω ˜ f i 1 + Δ t F i (39) f i x + ξ i Δ t , t + Δ t = f i 0 + 1 ω ˜ f i 1 + Δ t 2 F i + ψ i .

3. Hierarchical Grid Refinement Techniques

The stream and collide algorithm of the classical LBM relies on its use of rectangular and equidistant grids. To preserve this beneficial algorithm while simultaneously coping with increasing demands on computational efficiency and flexibility, several hierarchical refinement techniques for Cartesian grids are available in the literature. Despite rendering the grid globally irregular or non-uniform, they still maintain uniformity locally at each grid level. Characteristic of these schemes is the fact that particle populations at the various grid levels live on different spatial and temporal scales [5]. As addressed in Section 1, among these, three common methods, referred to as cell-vertex, cell-centered and combined approach within the scope of this work, represent the core of the numerical experiments later on. The specificities inherent to each of these methods will now be introduced.
In order to keep the following algorithm descriptions for all grid layouts vivid and clear, we will restrict ourselves to one additional level of refinement, thus two neighboring grid levels with spacings connected through Δ x c = 2 Δ x f . Here, superscripts c and f indicate the coarse and fine grid levels, respectively. To keep the molecular velocity ξ ˜ and as a consequence the isothermal speed of sound c s constant throughout the whole domain, convective scaling is used to link time steps between levels as Δ t c = 2 Δ t f . The fine grid therefore performs twice as many nested time steps as the coarse grid.

3.1. Cell-Vertex Grid Layout

Pioneering work regarding grid refinement in LBM has been achieved among many others by Filippova and Hänel [1], Lin and Lei [50], Dupuis and Chopard [2] and Yu et al. [51], viewing the numerical grid as conceptually based on a cell-vertex approach (abbr. cv), with the particle populations f i at all grid levels residing in cell corners, resulting in partially co-located coarse and fine nodes along grid transition interfaces, as depicted in Figure 3.
Grid interfaces extend over an overlap width of one coarse cell [51] and are characterized by interface nodes at the coarse grid Fluids 08 00103 i001 and the fine grid Fluids 08 00103 i002, hanging middle nodes Fluids 08 00103 i003 lacking neighbors on the same level in at least one direction and are therefore missing the corresponding populations after streaming. The missing information is recovered using co-located grid nodes, referred to as partner nodes Fluids 08 00103 i004 and Fluids 08 00103 i005, belonging to the other grid, under the condition that a smooth transition of hydrodynamic quantities such as density ρ , fluid momentum ρ u and the viscous stress tensor σ α β is ensured across the interface.
Equilibrium distribution functions f i ( 0 ) solely depend on ρ and u and can thus be simply transferred between interface nodes and their overlaying partners, due to grid level independency of ρ and u :
ρ = i = 0 q 1 f i ( 0 ) , c = i = 0 q 1 f i ( 0 ) , f
ρ u = i = 0 q 1 ξ i f i ( 0 ) , c = i = 0 q 1 ξ i f i ( 0 ) , f .
However, since the viscous stress tensor is related to the non-equilibrium part of f i , which itself involves the discrete lattice effect as mentioned in Section 2.3 (see Equation (26)), demanding equality of σ α β between grid levels leads to the subsequent correlations between coarse and fine non-equilibrium distribution functions:
σ α β c = ! σ α β f
ω ˜ c ν c s 2 Δ t c i = 0 q 1 ξ i ξ i f i n e q , c = ω ˜ f ν c s 2 Δ t f i = 0 q 1 ξ i ξ i f i n e q , f
ω ˜ c Δ t c f i n e q , c = ω ˜ f Δ t f f i n e q , f ,
with dimensionless collision frequencies on both levels denoted as ω ˜ c and ω ˜ f . Thus, when transferring populations from one grid level to another at interfaces, also known as grid coupling, non-equilibrium parts need to be rescaled accordingly, resulting in the following coarse-to-fine and fine-to-coarse grid coupling relationships:
c f : f i f = f i ( 0 ) + ω ˜ c 2 ω ˜ f f i n e q , c
f c : f i c = f i ( 0 ) + 2 ω ˜ f ω ˜ c f i n e q , f .
This single relaxation time grid-coupling procedure is inserted into the core algorithm after streaming, but before the collision step, in accordance with [2] and applied to all collision models in the numerical experiments for the cell-vertex approach, regardless of whether the collision process was performed in velocity or moment space.
Fine hanging nodes without a coarse partner will be referred to as middle nodes Fluids 08 00103 i003, if they are located on the edges of the coarse interface cell and center nodes Fluids 08 00103 i006, if they reside in a cell face center, as depicted in Figure 4 for three spatial dimensions. Missing populations at these nodes are determined through cubic p x = α = 0 3 a α x α [51] and bi-cubic p x , y = α = 0 3 β = 0 3 a α β x α y β [52] polynomial interpolation, respectively, using information of neighboring fine nodes with coarse partners. With stencil numeration as illustrated on the right side of Figure 4, interpolation formulas for middle and center nodes become
Fluids 08 00103 i007
Fluids 08 00103 i008
Another important aspect is a consequence of the nested time-stepping described above. Missing populations at both grid interfaces need to be provided before every respective collision step by applying grid coupling Equations (45) and (46). Nested time-stepping by definition involves stages during which coarse and fine grids are not synchronized in time, specifically at t + Δ t f no information on the coarse grid is available, since the coarse iteration advances at double the pace in the case of two levels. Therefore, coarse populations in Equation (45) are interpolated linearly in time during asynchronous iteration [53], before being transferred to the fine level:
Fluids 08 00103 i009
For scale-resolving turbulent flows, it is necessary to consider that part of the scales resolved at the fine level cannot be represented on the coarse grid due to the reduced spatial and temporal resolution of the latter. A naive transfer of quantities from a high-resolution region to one of reduced resolution would necessarily violate the Nyquist–Shannon criterion [54]. On this account, a restriction of fine scales during the f c (fine-to-coarse) transfer can be found in the literature [3,55,56], consisting of an additional stabilizing filtering of fine non-equilibrium distributions f i n e q , f in Equation (46) before scaling. Lagrava et al. [3] proposed the use of a simple averaging of the non-equilibrium parts over all lattice neighbors in this regard. For the D3Q19 lattice used in this paper, this results in:
f ¯ i n e q , f , LAG x , t = 1 19 j = 0 18 f i n e q , f x + ξ j Δ t f , t
Alternatively, Touil et al. [56] suggested a weighted average, depending on the lattice direction of the corresponding distributions:
f ¯ i n e q , f , TOU x , t = 1 7 f i n e q , f x , t + 1 14 j = 0 5 f i n e q , f x + ξ j Δ t f , t + 1 28 j = 6 17 f i n e q , f x + ξ j Δ t f , t .
As recommended in [3] all populations f i , not only the missing ones, are reconstructed during grid coupling at both coarse and fine interfaces. Even though this restriction operation is in principle advisable for turbulent flow regimes, it is still applied to laminar-regime test cases for the square-duct in Section 5 to investigate its effects on numerical stability.
With all the aforementioned steps combined, the multi-level collide-and-stream algorithm for the cell-vertex grid layout is summerized in Algorithm 1.
Algorithm 1: Cell-Vertex
Reference state  t c = t f = t : Distribution functions f i at all nodes are valid.
1:
Collision of all coarse and fine nodes.
2:
Streaming step of all coarse and fine nodes. Populations on the coarse grid reach t c = t + Δ t c , fine populations reach t f = t + 1 2 Δ t c .
3:
c f coupling: Reconstruction of f i f (45) and transfer onto Fluids 08 00103 i002 with temporal interpolation (49) of f i ( 0 ) and f i n e q , c as well as scaling of f i n e q , c (44).
4:
Spatial interpolation of f i f at Fluids 08 00103 i003 (47) and Fluids 08 00103 i006 (48). Asynchronous iteration is now complete.
5:
Collision of all fine nodes.
6:
Streaming step of all fine nodes. Populations on the fine grid reach t f = t + Δ t c .
7:
f c coupling: Reconstruction of f i c (45) and transfer onto Fluids 08 00103 i001 with optional filtering (50,51) and scaling of f i n e q , f (44).
8:
c f coupling: Reconstruction of f i f (45) and transfer onto Fluids 08 00103 i002 with scaling of f i n e q , c (44).
9:
Spatial interpolation of f i f at Fluids 08 00103 i003 and Fluids 08 00103 i006. Synchronous iteration is now complete.
10:
Return to Algorithm 1 and repeat until the end of the simulation.
This sequential process reproduces the algorithm if the collision models BGK, RR and MRT are used. When employing the HRR collision operator, some adjustments have to be made to the algorithm, which will be discussed in the following.

3.2. Peculiarities of the Hybrid-Recursive Regularized Collision Model

As already noted in Section 2.3, the use of the hybrid-recursive regularized collision model typically requires an approximation of the strain-rate tensor via finite differences, preferably of central type, due to its beneficial locality and accuracy properties. Given the abscence of same-level neighbors at interface nodes, approximated strain-rates need to be obtained from partner nodes (cf. Figure 3). The requirement of central differences, however, cannot be satisfied straight-forward during c f coupling at fine interface nodes Fluids 08 00103 i002, since coarse partners Fluids 08 00103 i005 have at least one direct neighbor Fluids 08 00103 i001 belonging to the coarse interface, i.e., without valid population state after advection. Thus, velocity information needed for a central difference approximation of S α β FD at Fluids 08 00103 i005, can not be provided by these neighbor nodes.
During synchronous iteration, when both the coarse and fine grid are in phase at t + Δ t c , this missing velocity information at the coarse interface nodes can be directly substituted with help of their fine partners by u (Fluids 08 00103 i004, t + Δtc).
Considering asynchronous iteration, when the fine grid reaches t + 1 2 Δ t c and thus grids are out of phase, S α β FD (Fluids 08 00103 i005, t + Δtf) needs to be specified and transferred to the fine interface nodes Fluids 08 00103 i002. Astoul et al. [12] suggested temporal interpolation of the strain-rate tensor, including a fictitious streaming step towards fine partner nodes Fluids 08 00103 i004 to synchronize them with regular coarse nodes Fluids 08 00103 i010, which are also utilized in the central difference stencil and temporally evolve directly from t to t + Δ t c . After the fictitious streaming step, a valid velocity vector can be extracted from the updated fine partners.
We propose a somewhat different approach. Temporal interpolation is performed at regular coarse nodes Fluids 08 00103 i010, in order to obtain velocity values at time t + Δ t f during the finite difference approximation of the strain-rate tensor at Fluids 08 00103 i005. At neighbors corresponding to coarse interface nodes Fluids 08 00103 i001 with incomplete population state after streaming, fine partner nodes Fluids 08 00103 i004 with populations already residing at t + Δ t f are used instead. The advantage of this proposal is, that no fictitious streaming is necessary anymore to temporally synchronize involved nodes.
For example, with the interface depicted in Figure 3 orthogonal to the x-direction, strain-rate S x x FD at coarse partners Fluids 08 00103 i005 would then be given by:
Fluids 08 00103 i011
where u (Fluids 08 00103 i010, t + Δtf) is computed through linear temporal interpolation
Fluids 08 00103 i012
For the coarse interface Fluids 08 00103 i001, no such problem arises during f c coupling, since the underlying fine partners are always surrounded by neighbors with valid states.
With these enhancements, the HRR adaption to the multi-level cell-vertex algorithm is summarized in Algorithm 2.
Algorithm 2: Cell-Vertex HRR
Reference state  t c = t f = t : Distribution functions f i at all nodes are valid.
1:
Collision of all coarse and fine nodes with hybrid-recursive regularized collision operator (27).
2:
Streaming step of all coarse and fine nodes. Populations on the coarse grid reach t c = t + Δ t c , fine populations reach t f = t + 1 2 Δ t c .
3:
c f coupling: Finite difference approximation of strain-rate-tensor S α β at coarse partners Fluids 08 00103 i005 according to (52) including temporal velocity interpolation (53), transfer onto Fluids 08 00103 i002. Reconstruction of f i f (45) and transfer onto Fluids 08 00103 i002 with temporal interpolation (49) of f i ( 0 ) and f i n e q , c as well as scaling of f i n e q , c (44).
4:
Spatial interpolation of f i f and S α β at Fluids 08 00103 i003 (47) and Fluids 08 00103 i006 (48). Asynchronous iteration is now complete.
5:
Collision of all fine nodes with hybrid-recursive regularized collision operator (27).
6:
Streaming step of all fine nodes. Populations on the fine grid reach t f = t + Δ t c .
7:
f c coupling: Finite difference approximation of strain-rate-tensor S α β at fine partners Fluids 08 00103 i004, transfer onto Fluids 08 00103 i001. Reconstruction of f i c (45) and transfer onto Fluids 08 00103 i001 with optional filtering (50,51) and scaling of f i n e q , f (44).
8:
c f coupling: Finite difference approximation of strain-rate-tensor S α β at coarse partners Fluids 08 00103 i005, transfer onto Fluids 08 00103 i002. Reconstruction of f i f (45) and transfer onto Fluids 08 00103 i002 with scaling of f i n e q , c (44).
9:
Spatial interpolation of f i f and S α β at Fluids 08 00103 i003 and Fluids 08 00103 i006. Synchronous iteration is now complete.
10:
Return to Algorithm 1 and repeat until the end of the simulation.

3.3. Cell-Centered Grid Layout

An alternative scheme for grid coupling in LBM with hierarchically refined grids, relies on a cell-centered (abbr. cc) or volume-based arrangement of grid nodes [4,5]. Although continuity of macroscopic physical quantities ( ρ , ρ u , σ α β ) is ensured across grid transitions by rescaling non-equilibrium distribution functions f i n e q during the coupling steps in the cell-vertex algorithm previously outlined in Section 3.1; mass and momentum conservation may be violated, due to the choice of grid layout in connection with the involved spatial and temporal interpolation methods [3,57]. Resorting to a volumetric description, with grid nodes at all levels located in their respective cell centers and particles considered as mass rather than densities, being transferred between grids of various resolution, guarantees mass- and momentum conservation on a global scale in a natural way [5].
Grid coupling in the volumetric approach is characterized by two specific steps, as vividly illustrated using an exemplary interface region in Figure 5. During c f coupling, post-collision particles traveling towards the fine grid and originating from the coarse cell-centered interface nodes Fluids 08 00103 i013, are redistributed among fine interface nodes Fluids 08 00103 i014 contained within the same coarse parent cell. This procedure will be referred to as uniform explosion, since particle mass is homogeneously distributed, without any type of spatial or temporal interpolation applied. Since two rows of fine interface nodes are supplied with coarse post-collision states, they remain valid for two consecutive fine time steps. Furthermore, according to [5] no rescaling of non-equilibrium distributions is performed, in contrast to the cell-vertex method in Algorithm 1. To keep the algorithmic description consistent between grid layouts, mathematical notation of grid coupling steps in the cell-centered scheme will be expressed in terms of particle densities, rather than masses. Explosion then implies equality between particle densities, and coarse provider and fine receivers. With post-collision distribution functions symbolized by f i , the explosion step is simply given as:
c f : f i f x f , t = f i c x c , t .
Here, x f and x c represent position vectors of the involved fine and coarse nodes.
As for f c coupling, particle masses traveling towards the coarse grid are gathered at fine interface grid nodes Fluids 08 00103 i014 inside the parent cell after two consecutive fine cycles and transferred to the coarse interface node Fluids 08 00103 i013 residing in the cell center. In this paper, fine interface nodes with at least one regular fine neighbor Fluids 08 00103 i015 will be referred to as first-layer nodes, whereas second-layer nodes means fine interface node devoid of regular fine neighbors. Particle collision does not take place at first- or second-layer nodes, so that collision has to be performed at the coarse receiver node after f c coupling. In terms of particle densities, this as coalescence known procedure implies missing coarse densities to be obtained by averaging fine particle densities contained within the coarse cell, resulting in the following expression for D spatial dimensions:
Fluids 08 00103 i016
Here again, no rescaling of transferred non-equilibrium distribution functions is performed during the coupling step in accordance with [5]. In contrast to the direct, i.e., without filtering, f c communication step (cf. Equation (46)) in the cell-vertex approach in Algorithm 1, coalescence functions as implicit filtering of fine distributions, as a result of spatial equalization of grids in the cell-centered layout [58].
These main differences in comparison to the grid transition algorithm for the cell-vertex layout, consisting in the absence of any higher-order spatiotemporal interpolation or explicit rescaling of non-equilibrium distribution functions during information exchange in the coupling step and implicit filtering during f c communication, renders the computational implementation of the cell-centered approach rather simple. Still, numerical errors depending on various factors, such as viscosity, body-force, grid spacing and time step are introduced by the method and can be comparatively large, especially for grid interfaces oriented perpendicular to the flow direction, as was shown by Rohde et al. [5]. Moreover, Rohde et al. assume a scaling of the non-equilibrium part to be implicitly present in their method.
An improved accuracy of the coarse-to-fine coupling in the cell-centered approach can be achieved by using a linear distribution of coarse post-collision states during explosion, as proposed by Chen et al. [59]:
f i f x f , t = f i c x c , t + x f x c · F i x c , t ξ i ξ i · F i x c , t | ξ i | 2 ,
where the vector function F i (not to be confused with the force term F i in Section 2.4) represents gradients of coarse post-collision states f i c along directions α parallel to the interface, so that both, x c + e α Δ x c and x c e α Δ x c point to locations of coarse interface nodes in the immediate vicinity of Fluids 08 00103 i013 x c , t (cf. Figure 5):
F i α x c , t = f i c x c + e α Δ x c , t f i c x c e α Δ x c , t 2 Δ x c .
Gradients in directions orthogonal to the interface are disregarded by setting:
F i α x c , t = 0 .
This procedure is referred to as linear explosion in the course of this work. Both types of explosion (54) and (56) are included in the numerical experiments in Section 5. The resulting algorithm for the cell-centered grid layout is summarized in Algorithm 3.
Algorithm 3: Cell-Centered
Reference state  t c = t f = t : Distribution functions f i at all nodes are valid.
1:
Collision of all coarse nodes Fluids 08 00103 i013 Fluids 08 00103 i010 and regular fine nodes Fluids 08 00103 i015 (see Figure 5). Fine interface nodes Fluids 08 00103 i014 do not take part in any collision step in the cell-centered algorithm.
2:
c f coupling (explosion): Redistribution of coarse post-collision populations pointing towards the fine grid among fine interface nodes Fluids 08 00103 i014 according to Equation (54) or (56). Neither rescaling nor temporal interpolation of distribution functions is performed.
3:
Streaming step of all coarse and fine nodes. Coarse populations reach t c = t + Δ t c , populations on the fine grid reach t f = t + 1 2 Δ t c . Asynchronous iteration is now complete.
4:
Collision at regular fine nodes Fluids 08 00103 i015. Since two lines of fine interface nodes are updated during explosion with coarse post-collision states, the corresponding states remain valid for two consecutive fine streaming steps.
5:
Streaming step of all fine nodes except fine second-layer interface nodes, since these nodes lack valid states at this point. Populations on the fine grid reach t f = t + Δ t c .
6:
f c coupling (coalescence): Averaging of populations belonging to fine interface nodes Fluids 08 00103 i014 pointing towards the coarse grid and redistribution to coarse interface node Fluids 08 00103 i013 residing in corresponding cell center according to Equation (55). Synchronous iteration is now complete.
7:
Return to Algorithm 1 and repeat until the end of the simulation.
Qi et al. [60] improved the accuracy during the volumetric coarse-to-fine coupling by applying a compact gradient-based interpolation method based on a second-order polynomial expression for the macroscopic velocity very similar to the one used for the combined grid layout in this paper. In contrast to the uniform and linear explosion in Equations (54) and (56), their method relies on a grid interface overlap of two coarse cells and a preceding execution of coalescence to ensure that all coarse source nodes are valid. Another difference is that non-equilibrium distribution functions are scaled according to Equation (44) during both coupling steps. This scheme was not implemented as part of the present study.

3.4. Combined Grid Layout

Another variant of grid arrangement can be constructed as a combination of both aforementioned layouts (cf. Figure 6), with fine interface nodes Fluids 08 00103 i014 residing in their respective cell centers being enclosed by regular coarse nodes Fluids 08 00103 i010 in cell corners, referred to as combined method in the course of this work (abbr. cm). Similarly, coarse interface nodes Fluids 08 00103 i013 are surrounded by regular fine cell-centered nodes Fluids 08 00103 i015 [7,8,58,61], as depicted in Figure 6 in two dimensions. Unlike in the cell-vertex and cell-centered layouts, here, the grid transition interface has an overlapping width of two coarse cells. Distribution functions at interface nodes are then obtained through a decomposition into equilibrium and non-equilibrium parts, utilizing gradient-based compact interpolation schemes to preserve second-order accuracy of the lattice Boltzmann method. Equilibrium distribution functions are calculated by applying Equation (9) with interpolated values of density and velocity.
There are two aspects that have to be taken into account here. First, in order to recover correct hydrodynamic behavior, second-order derivatives of the momentum field in the Navier–Stokes equations need to be sustained by the grid coupling method. Secondly, since the LBM used here operates in the weakly compressible and isothermal limit and only the first-order pressure gradient is present in the Navier–Stokes equations, it is sufficient to interpolate the density field in a linear manner [7], that means trilinearly in our case.
Regarding the first point, we use a three-dimensional version of the compact gradient-based velocity interpolation proposed in [62], resulting in the following third-order polynomial expression for Cartesian velocity vector components u α I :
u α I x , y , z = a α , 000 + a α , 100 x + a α , 010 y + a α , 001 z + a α , 200 x 2 + a α , 110 x y + a α , 101 x z + a α , 020 y 2 + a α , 011 y z + a α , 002 z 2 + a α , 111 x y z ,
with α = x , y , z and where superscript I indicates an interpolated value. The 33 interpolation coefficients a α , i j k are obtained utilizing velocity gradient information contained locally within strain-rates, for details see Appendix C. Furthermore, besides the density field, non-equilibrium distribution functions are also interpolated trilinearly. On this matter, our method is similar to the one proposed by Qi et al. [60], albeit they rely on an even more compact stencil of only four source nodes for their gradient-based second-order polynomial interpolation in three dimensions. This leads to the following reconstruction of distribution functions at interface nodes during the grid coupling procedures:
c f : f i f = f i ( 0 ) ρ I , u I + ω ˜ c 2 ω ˜ f f i n e q , c , I
f c : f i c = f i ( 0 ) ρ I , u I + 2 ω ˜ f ω ˜ c f i n e q , f , I ,
corresponding to the pre-collision state, so that collision on both grids has to be performed afterwards. Similar to the cell-vertex approach, non-equilibrium functions are rescaled before being transferred to the receiving nodes, however, as mentioned in Section 3.3, filtering of f i n e q , f during f c communication is achieved implicitly as a result of the spatial offset between grids, as in the cell-centered method. Furthermore, temporal interpolation during c f coupling is redundant, since two layers of fine interface nodes are updated and thus remain valid for two consecutive fine cycles. The combined scheme merges both previous approaches together with a beneficial compact interpolation method [58].
The necessary strain-rates in Equations (A9) to (A20) can be obtained as either second-order non-equilibrium moments combining Equations (22) and (26) or alternatively from finite differences. We rely on the first option.
The collide-and-stream algorithm for the combined scheme is summarized in Algorithm 4.
Algorithm 4: Combined
Reference state  t c = t f = t : Distribution functions f i at all nodes are valid.
1:
Collision of all coarse and fine nodes.
2:
Streaming step of all coarse and fine nodes. Populations on the coarse grid reach t c = t + Δ t c , fine populations reach t f = t + 1 2 Δ t c . No grid coupling is performed in the asynchronous iteration, hence asynchronous iteration is now complete.
3:
Collision of all fine nodes except second-layer fine interface nodes, since these nodes lack valid states at this point.
4:
Streaming step of all fine nodes except second-layer fine interface nodes. Populations on the fine grid reach t f = t + Δ t c .
5:
f c coupling according to (61), reconstruction of f i f at Fluids 08 00103 i014 (Figure 6) including trilinear interpolation of ρ and compact gradient-based interpolation of u (59) for the calculation of equilibrium distributions, trilinear interpolation and scaling of f i n e q , c .
6:
c f coupling according to (60), reconstruction of f i c at Fluids 08 00103 i013 (Figure 6) including trilinear interpolation of ρ and compact gradient-based interpolation of u (59) for the calculation of equilibrium distributions, trilinear interpolation and scaling of f i n e q , f . Since two lines of fine interface nodes are updated during this procedure, the corresponding states remain valid for two consecutive fine time steps. Synchronous iteration is now complete.
7:
Return to Algorithm 1 and repeat until the end of the simulation.

4. Adaption of the Hybrid-Recursive Regularized Collision Model to Cell-Centered and Combined Grid Layouts

As mentioned in Section 3.2, a homogenous approximation of strain-rates through second-order central finite differences within the hybrid-recursive regularized collision model, requires some particularities to be considered during the collide-and-stream process in the cell-vertex approach. Similarly, some adjustments need to be made for the cell-centered and combined method, in order to ensure a consistent use of central differences across the grid transition interface.
In the cell-centered procedure summarized in Algorithm 3, collision in the grid transition region is performed exclusively at coarse interface nodes Fluids 08 00103 i013 (cf. Figure 5). Coarse post-collision populations f i c traveling towards the fine grid are distributed among fine interface nodes by either uniform or linear explosion, cf. Equation (54) and (56), respectively. However, in order to perform HRR collision beforehand, strain-rates at Fluids 08 00103 i013 need to be approximated, preferably via central finite differences. Hence, absent coarse neighbor nodes must be replaced. Taking an exemplary interface portion as depicted in Figure 7, we substitute the missing coarse neighbor cell in the fine grid region Fluids 08 00103 i017, using a fictitious coalescence at regular fine nodes covering its image. Fine pre-collision distribution functions in all directions are averaged according to Equation (55) and thus provided for a calculation of the macroscopic velocity within the strain-rate computation at Fluids 08 00103 i013. A similar approach was adopted by Premnath et al. [63] in order to obtain test-filtered quantities across grid transition interfaces during dynamic subgrid scale modeling.
Together with this, two further modifications regarding the cell-centered algorithm previously discussed, need to be considered. First, collision at coarse nodes with subsequent explosion, i.e., part of step 1 and step 2 of the cell-centered algorithm, need to occur in advance to the collision of regular fine nodes Fluids 08 00103 i015. Otherwise, invalid population states of first-layer interface nodes Fluids 08 00103 i014 would be employed in the central difference approximation of the strain-rate tensor S FD at neighboring regular fine nodes. Secondly, in contrast to step 2 in Algorithm 1, a complete set of coarse post-collision distribution functions is exploded onto the first-layer interface nodes. These additional adjustments keep the HRR procedure consistent with the rest of the cell-centered approach and are used for the HRR test simulations later on. Steps 3 to 7 of the cell-centered algorithm remain unchanged.
For the combined method, similar adaptions, consistent with the remainder of the combined algorithm in Section 3.4, are proposed. Here, collision is performed at coarse as well as first- and second-layer fine interface nodes, so that missing velocity values during central difference calculation need to be reconstructed on both sides of the interface by appropriate measures. On the fine grid, only second-layer interface nodes lack necessary neighbors. For this purpose, compact velocity interpolation according to Equation (59) is applied on positions indicated by coarse Fluids 08 00103 i017 and fine Fluids 08 00103 i018 ghost nodes in Figure 8, to achieve a central finite difference stencil relative to second-layer and coarse interface nodes. This modification is included into step 1 of the combined Algorithm 4.
There is an additional aspect that must be taken into account here. Right after the asynchronous iteration is completed in step 2 of Algorithm 4, second-layer interface nodes are excluded from any further collision and propagation, due to incomplete population state. Nevertheless, valid velocity values have to be provided by these nodes in order to compute S FD at first-layer interface nodes during HRR collision throughout synchronous iteration. Since populations on the fine grid still live in t + Δ t f at this point, unlike the coarse interpolation sources, which progress from t directly to t + Δ t c , we combine the compact spatial interpolation with a linear temporal interpolation similar to Equation (49) to obtain u (Fluids 08 00103 i014, t + Δtf) at second-layer interface nodes. Steps 3 to 7 of Algorithm 4 remain unchanged.
In Section 2.3, cubic Mach correction for the D3Q19 lattice was mentioned. According to Equation (A5), spatial derivatives of deviation terms Ψ are approximated via central finite differences, similar to strain-rates in Equation (28). All the algorithmic adaptions discussed previously in this section apply in the same manner to the discretized gradients of the deviation terms to obtain correction terms ψ i . With these consistent modifications, the HRR/ HRR ψ procedure for the combined and cell-centered layouts is complete.

5. Numerical Experiments and Results

5.1. Test Case Setup and Evaluation Strategy

Numerical experiments for the stability evaluation of the grid transition algorithms and collision operators discussed in the previous chapters are performed using a force-driven square duct flow test case, schematically depicted in Figure 9. Within subcritical Reynolds number regime, the analytical solution for the laminar velocity field is given by [64]:
u y , z = 16 h 2 ν π 3 a x i = 1 , 3 , 5 , . . . 1 i 1 2 1 cosh i π z 2 h cosh i π 2 cos i π y 2 h i 3
Here, ν is the kinematic viscosity and a x represents the uniform flow acceleration acting in positive x-direction, i.e., parallel to the duct center axis.
The simulations are performed in a periodic, cubic domain with characteristic width 2 h , see Figure 9. A single level reference case (abbr. sl) with grid spacing 10 Δ x c = 2 h is used for comparison. In the multi-level case, the cubic domain is discretized using two levels of grid resolution Δ x c = 2 Δ x f , with coarse cells located in the duct core and six refinement layers of fine cells at the surrounding duct walls, resulting in approximately 10 4 cells in total. This rather coarse grid resolution was chosen to keep the time frame for the numerical study within acceptable bounds. The grid transition interface with an overlap width of Δ x c for cell-vertex and cell-centered cases and an overlap width of 2 Δ x c for the combined grid layout, is oriented parallel to the flow direction througout the complete domain. No-slip walls are implemented utilizing simple half-way bounce-back scheme, while periodic conditions enclose the domain at boundaries orthogonal to the flow direction. All simulations have been performed with the in-house lattice Boltzmann package SAM-Lattice.
As was shown by means of linear stability analysis for various collision models, the lattice Boltzmann method exhibits instabilities for high Reynolds number flows [9,10], which materialize in the form of physically meaningless negative values for the distribution functions f i . As can be seen directly from the connection between dimensionless collision frequency and kinematic viscosity in Equation (8), a Reynolds number tending towards infinity is associated with ω ˜ 2 . In order to examine the influence of the various grid transition schemes and collision models on accuracy and stability ranges for the square-duct flow, numerical experiments are performed. Through successive variation of the flow acceleration a x , Reynolds number and consequently the collision frequency ω ˜ c are modified together with simultaneous verification of adherence to the following constraints, defining three distinct evaluation categories:
  • Category 1: Simulation remains stable, i.e., f i > 0 , and within the subcritical Reynolds number regime Re b , h < Re b , h , crit 1673 [65] 2060 [66]. Furthermore, iterative convergence u x , t + Δ t u x , t 1 × 10 15 m   s 1 is reached with an averaged relative velocity error of δ u ¯ = 1 N n = 0 N 1 | u x n , t u x n , t analytical | u x n , t analytical < 5 % .
  • Category 2: Simulation remains stable, i.e., f i > 0 , but either averaged relative error increases to δ u ¯ 5 % in subcritical regime or Re b , h > Re b , h , crit .
  • Category 3: Simulation is unstable and diverges, i.e., f i < 0 is detected.
Here, the Reynolds number Re b , h = U b D h / ν is based on the bulk velocity
U b = a x h 2 3 ν 1 192 π 5 i = 1 , 3 , 5 , . . . tanh i π / 2 i 5
and hydraulic equivalent diameter D h = 2 h for the square duct flow. The simulation parameters are summarized in Table 1. Affiliation of a particular test case to one of these categories is expressed in terms of dimensionless collision frequency ω ˜ 1 / 2 c , where superscript c indicates a coarse grid value as before and subscript 1 or 2 represents the respective category. Regarding the first category, the procedure stops once the relative error increases to δ u ¯ 5 % in the subcritical regime and the condition for iterative convergence can no longer be met or Re b , h > Re b , h , crit and thus a comparison with the analytical solution Equation (62) is not permissible anymore. Consequently, in the remainder of this paper ω ˜ 1 c represents the maximum value of ω ˜ c that just about satisfies these critera. As for the second category, ω ˜ c variation terminates as soon as negative distribution functions f i < 0 emerge. Similar to ω ˜ 1 c for category 1, the value ω ˜ 2 c reflects the maximum dimensionless collision frequency for which this scenario has just not yet occurred up to a terminal value of ω ˜ c = 1.99999 , combined with the additional requirement that this must be fulfilled for at least 3 s of physical flow time. It should be mentioned that the relatively coarse grid is not suitable for high-accuracy simulations of turbulent flows in the sense of insufficient resolution. However, the focus here is on finding the limits of the method or seeing how far it can be pushed for a given configuration, as stability of the numerical scheme is of high significance for application to flows of engineering interest.
The choice of ω ˜ c , i.e., dimensionless collision frequency on the coarse grid instead of ω ˜ f for the evaluation is arbitrary. To remain in the incompressible flow regime, the Mach number is kept constantly low at Ma = 0.1 .
An overview of case variants examined in the numerical study is given in Table 2. To emphasize the influence of both, the collision model and the grid transition scheme on the accuracy and stability range, each of the collision models reviewed in Section 2 are tested together with the schemes addressed in Section 3 and Section 4 as well as on an equidistant, i.e., single level grid (abbr. sl). For the cv algorithm with BGK collision, two additional calculations are performed including restriction operations as expressed in Equations (50) and (51). Furthermore, the HRR model is tested additionally with cubic Mach number correction, referred to as HRR ψ . The hybridization factor σ in Equation (27) is set to 0.98 for all HRR cases [12,16,41], that means 2 % of the stress tensor calculation in the reconstruction of non-equilibrium distributions is achieved through centered finite differences, while remaining 98 % stem from non-equilibrium moment computation. Regarding the tests of cc algorithm with BGK and HRR models, explosion is executed with both, uniform and linearly interpolated coarse post-collision particle distribution, as expressed in Equations (54) and (56), respectively. All remaining cc simulations are performed with uniform explosion only. Together with the procedure for the variation of ω ˜ c described above, this results in a total number of simulations of approximately 500.

5.2. Verification

The conducted numerical experiments are verified by comparison with the analytical solution. Figure 10 shows normalized velocity profiles ( u ^ analytical = max u y , z in Equation (62)) along y direction normalized with half channel height h, for all HRR cases with σ = 0.98 and ω ˜ c = 1.94990 , corresponding to a bulk Reynolds number of Re h , b = 67.53 . All simulation results agree very well with the solution in Equation (62), supporting the correct numerical implementation of the adaptation of the HRR collision operator to the cell-centered and combined types of grid transition interface layout, as discussed in Section 4. Velocity profiles for all other tested models exhibit similar agreement with the analytical solution and are shown in Figure A2.
Numerical accuracy of the various grid transition schemes is further examined for the HRR collision operator by means of a grid convergence study, summarized in Figure 11 in a double logarithmic manner. Results of the global absolute error in velocity averaged over all N nodes, defined as
Δ u RMS = 1 N n = 0 N 1 u x n , t u analytical x n , t 2
are given together with a slope of 1 and 2 for three successively refined grids, starting from the coarsest grid, which serves as the basis for stability experiments. Grid resolution is expressed in terms of number of respective fine level grid cells per duct height. The time step is determined by diffusive scaling Δ t Δ x 2 in order to ensure O Δ x 2 convergence of the compressibility error, see Table 3.
The validity and correctness of the presented HRR adaption for the cc and cm algorithms is again confirmed. The cc algorithm with linear explosion (Fluids 08 00103 i019) shows an average slope of 2.03 , while the cm HRR adaption (Fluids 08 00103 i020) results in an average convergence rate of 2.30 with a slight tendency of hyper convergence between the coarsest and medium grids, preserving the second order spatial accuracy of the LB scheme in a consistent manner. Relying on uniform particle redistribution during cc explosion (Fluids 08 00103 i021), expectedly degrades the methods accuracy to first order with an average slope of 1.11 [59]. Choosing a relatively low Mach number of Ma = 0.1 no difference between cv HRR including (Fluids 08 00103 i022) and excluding (Fluids 08 00103 i023) cubic correction terms ψ i could be observed here, with both variants exhibiting an average slope of 2.04 . The single level method (Fluids 08 00103 i024) yields an average convergence rate of 2.00 .

5.3. Accuracy and Stability Range Comparison

As mentioned, accuracy and stability range are considered in the context of the categories defined in Section 5.1 for the test cases summarized in Table 2. The decisive factor is the respective maximum value of ω ˜ c that just fulfills the conditions imposed on the category. Thus, for a concrete test case, ω ˜ 1 c describes the maximum value of ω ˜ c with a precision of five significant digits, whereby the simulation satisfies the constraints u x , t + Δ t u x , t 1 × 10 15 m   s 1 and δ u ¯ < 5 % within a laminar flow regime. Similarly for ω ˜ 2 c , the particular test case under consideration is at the stability limit, so it does not quite fall into the third category, since f i < 0 does not occur.
Figure 12 shows the maximum values ω ˜ 1 c and ω ˜ 2 c defined in this way for the BGK collision model for different grid layouts. It can be seen that the cell-vertex grid without additional restriction of the fine distributions denoted as cv-0 has the lowest values in both categories, with ω ˜ 1 c = 1.98100 ( Re h , b = 180.92 ) in category 1 and a clear distance to the other approaches in category 2 with ω ˜ 2 c = 1.98312 ( Re h , b = 203.85 ).
The cell-centered approach cc-5 with linear interpolation of the coarse post-collision distributions in the explosion step shows the highest stability of ω ˜ 2 c = 1.99961 ( Re h , b = 9256.56 ) for the BGK collision model, closely followed by cc-0, i.e., uniform explosion, with ω ˜ 2 c = 1.99940 ( Re h , b = 6183.25 ). Even though in the literature [3,56] physical justification for the f c restriction step originally consisted in filtering out non-resolvable fine scales in turbulent flows, accuracy and stability range for the cv algorithm can be increased in our test case within laminar flow regime. Regarding category 1, the accuracy range can be expanded by including either filtering operation from Equations (50) or (51) into step 7 of Algorithm 1 to ω ˜ 1 c = 1.99046 ( Re h , b = 362.10 ) and ω ˜ 1 c = 1.99060 ( Re h , b = 367.53 ), respectively, exceeding the corresponding values of cc and cm. Stability range is increased to ω ˜ 2 c = 1.99463 ( Re h , b = 644.72 ) with both filters, which is in the immediate vicinity of ω ˜ 2 c for the cm algorithm. The onset of instability is connected to the occurrence of spurious oscillations in the transverse velocity components v and w, the amplitude of which increases with time in the case of the cell-vertex algorithm without restriction (cv-0). This is shown in Figure 13 along normalized y direction for three exemplary chosen successive time steps and in comparison with the corresponding velocity values for cc-0, cm-0 and cv-5 (filtering according to Equation (50)) at ω ˜ c = 1.98317 ( Re h , b = 205.04 ), i.e., beyond cv-0’s stability limit. The spurious oscillation in transversal velocity w is depicted only in time step t = 8.5   s since its amplitude is several orders of magnitude below the v-oscillation until that point and then increases rapidly.
When filtering is activated in the cv algorithm during the fine-to-coarse coupling step, as in cv-5, a significant suppression of these spurious oscillations comes into effect, yielding the simulation more stable. Due to the spatial separation of the fine and coarse interface nodes and the resulting implicit filtering of the distributions during the respective fine-to-coarse coupling (cf. Equations (55) and (61)), this suppression of spurious oscillations occurs in equal form in the cc and cm layout. Differences in the curve characteristics at t = 2.3   s between cv-5 and cm-0 on the one hand and cc-0 on the other, are possibly associated with missing scaling of f i n e q in the volumetric description.
As for the accuracy and stability range maps of the MRT collision model, both sets of relaxation parameters, the one obtained by von Neumann analysis and the one corresponding to a regularization step, are compared in Figure 14. For all grid layouts except cc, MRT-REG (cases labeled -1) exhibits a slightly superior accuracy as well as stability range in the numerical experiments compared to MRT-NEU (cases labeled -2), while both MRT methods clearly surpass the BGK model in terms of stability. Even though MRT-NEU is explicitly calibrated to yield higher stability, regularizing the scheme by direct equilibriation of higher-order non-hydrodynamic moments as in MRT-REG seems to allow for further moderate stability gain. Between all grid transition types, the volumetric description again performs the strongest in terms of stability and also accuracy range with both MRT-REG and MRT-NEU exhibiting negative distribution functions at our terminal value ω ˜ 2 c = 1.99999 ( Re h , b = 3.470 × 10 5 ), followed this time by the cell-vertex algorithm. Differences to the cm algorithm are, however, of minor extent. Since stabilizing fine-to-coarse restriction is not active in our cv MRT test cases (cf. case overview in Table 2), the stability increase with MRT-NEU and MRT-REG is probably associated to two aspects: An increase in bulk viscosity and non-hydrodynamic mode filtering, weakening any destructive and parasitic contributions that may occur, materializing in the form of spurious oscillations as seen with the BGK model. Any harmful artifact is thus attenuated before it can reach the grid transition interface to be amplified.
Figure 15 summarizes the values of ω ˜ 1 c and ω ˜ 2 c of the RR and HRR collision models for all grid layouts. Here again, differences between the cv algorithm on one side and cc and cm on the other in terms of stability range are overcome by the positive effects of the collision operator: The cell-centered algorithm leads in terms of accuracy range for both, the RR and HRR model and in terms of stability range for the RR model, followed again by the cell-vertex and combined layouts. It should be briefly mentioned here again, that the cc tests for all advanced collision operators except HRR have been contucted exclusively with uniform explosion (cf. Equation (54)), that means with a degraded first order accuracy as can be seen in the grid convergence study for the HRR model in Figure 11. Nonetheless, the differences in terms of stability range compared to linear explosion in Equation (56) are negligable, as discussed above and depicted in Figure 12 for BGK collision. Similarly, no difference has been observed in terms of stability range for HRR with linear explosion and only minor differences are present in the accuracy range map. For lucidity reasons, test case cc-6 is not depicted in Figure 15 and an extension of the statements made here regarding cc’s stability range with uniform explosion can be attributed in a straightforward way to its linear enhancement.
The recursive-regularization step aims at filtering out non-hydrodynamic modes by explicit reconstruction of f i ( 1 ) according to Equation (23). Although the basic idea is the same as in the MRT-REG model discussed earlier, differences in the stability ranges can be observed when both methods are compared to each other, more precisely, the RR model does exhibit slightly reduced stability. For example, cv MRT-REG remains in category 2 (cf. Section 5.1) up to ω ˜ 2 c = 1.99970 ( Re h , b = 1.176 × 10 4 ), while cv RR remains stable only up to ω ˜ 2 c = 1.99883 ( Re h , b = 2.956 × 10 3 ), even though a third-order Hermite expansion of the equilibrium functions is used as compared to a second-order one in MRT model. An explanation for this can be found in the fact that the MRT model utilized here relies on a raw moment space built using the Gram–Schmidt orthogonalization procedure based on an unweighted scalar product [30], while the populations in the RR framework are derived based on Gauss–Hermite quadrature. The specific choice of relaxation parameters in the MRT-REG model leads to an increased bulk viscosity compared to the RR model, contributing to the formers stability.
For the refined cases the highest values in terms of stability and accuracy range can be found for the Hybrid-Recursive Regularized collision model, regardless of the type of grid transition in use. All conducted HRR simulations remained in category 2 in our study, meaning no negative distribution functions have been detected up to the terminal value, thus yielding ω ˜ 2 , HRR c > 1.99999 ( Re h , b = 3.470 × 10 5 ) in all cases. Besides the non-hydrodynamic mode filtering properties of the HRR model [12], its accompanying drastic stability increase is related to the injection of numerical hyper-viscosity through the usage of a hybridized computation of the stress tensor during the reconstruction of f i ( 1 ) , as expressed in Equation (27) [16,31]. For the single level reference case however, the finite difference approximation of strain-rates within the stress tensor together with the rather coarse grid resolution leads to an early increase of the relative arithmetic mean error δ u ¯ beyond the 5 % limit of category 1, as can be seen in Figure 15 and Figure 16a. The direct adjustment of the amount of numerical viscosity ν ¯ num through the hybridization factor σ is depicted in Figure 17 with ω ˜ c = 1.94990 , corrsponding to a bulk Reynolds number of Re h , b = 67.53 . By inserting the numerical velocity into the analytical solution of the square duct flow (62), solving for kinematic viscosity ν and averaging over all nodes of the computational domain, the relative amount of total viscosity ν ¯ total / ν = ν ¯ num + ν / ν is estimated for all collision models and grid transition algorithms. Interestingly, the RR and MRT-REG model lead to similar amounts of total viscosity for all grids, with the MRT-REG slightly above RR, besides showing differences in the stability ranges.
In addition to the values for the HRR hybridization factor used in the stability study, i.e., σ = 1 (corresponding to the RR model, cf. Equation (27)) and σ = 0.98 , the value of σ = 0.99 is included. The linear trend of numerically induced hyper-viscosity can be clearly seen by comparing the RR ( σ = 1 ) and two HRR cases with σ = 0.99 and σ = 0.98 , i.e., 0 % , 1 % and 2 % of stress tensor approximation through central finite differences, respectively. Whether the distinct stability gain with the HRR model could have been achieved with even higher values of hybridization factor, i.e., less added dissipation, will be left open for future work. Comparing the average levels of added numerical viscosity between individual grid transition schemes in Figure 17 may give an explanation for the improved stability properties of the cc algorithm. The increased amount of dissipation introduced into the flow field relative to the other two approaches could be explained by the implicit filtering of the complete distribution function f i during coalescence step 6 of Algorithm 3, as opposed to filtering only the non-equilibrium part f i n e q in step 7 of Algorithm 1 and step 5 of Algorithm 4. Lagrava et al. [3] mentioned that they observed a strong dissipation by application of their averaging procedure in Equation (50) to f i or the macroscopic quantities ρ and u , which lead to the choice of filtering only f i n e q instead. Even though ρ is similarly interpolated trilinearly in our cm approach for the reconstruction of equilibrium functions f i ( 0 ) on interface nodes, u is obtained through a third-order polynomial expression. Moreover, f i n e q instead of f i is interpolated trilinearly here. However, further research needs to be conducted on this regard.
While moderate differences between the grid transition algorithms are evident in terms of accuracy and stability range in Figure 12, Figure 13, Figure 14 and Figure 15 with the volumetric approach showing the highest values for ω ˜ 2 c , the influence of the collision operator is considerably greater in this respect. This is made even more apparent by a rearrangement of the range maps juxtaposing different collision operators for each individual grid transition algorithm as it is shown in Figure 16 and Figure 18. Regardless of the specific interface type in use, the highest values of ω ˜ 2 c is achieved exclusively with the HRR collision model, extending beyond the scope of ω ˜ 2 c > 1.99999 ( Re h , b = 3.470 × 10 5 ) of our numerical experiments, closely followed by the two MRT models and then RR, with MRT-REG slightly above MRT-NEU and lastly BGK.

6. Conclusions and Future Work

In this paper, we compared three established lattice Boltzmann grid transition algorithms for hierarchically refined grids in terms of stability and accuracy by conducting a series of numerical experiments for a square duct flow test case. Depending on the specific arrangement of the nodes on the individual grid levels relative to each other, a distinction can be made between the cell-vertex (cv), cell-centered (cc), and a combined approach (cm). In order to simultaneously assess the influence of the collision operator, the classical BGK model was included in the study in addition to three advanced collision models, namely a Multi Relaxation Time model (MRT) with two sets of relaxation times and the Recursive Regularized (RR) and Hybrid-Recursive Regularized (HRR) models. HRR collision involves hybridization of the stress tensor calculation in the Hermite series reconstruction of first-order non-equilibrium functions f i ( 1 ) through the usage of central finite differences. Since this is not straightforward and the three grid types differ in the communication between fine and coarse neighbor nodes in the grid transition region, an adaptation of the HRR procedure, which was originally presented for cv type grids, to cc and cm layouts was introduced and validated by comparison with the analytical solution for the duct flow and a grid convergence study. The adapted procedures demonstrated excellent agreement between numerical and analytical results and conserved second-order spatial accuracy of the LB scheme.
A total of three categories were introduced for an assessment of the accuracy and stability range. The first two categories differ in whether the relative deviation of the calculated duct flow to its analytical solution exceeds a fixed barrier of five percent, defining the accuracy range, while the transition from the second to the third category requires the occurrence of unphysical negative values of distribution functions f i , defining the stability range. For each case in the study, place value within the category was evaluated by finding the maximum value of the dimensionless coarse grid collision frequency ω ˜ 1 c and ω ˜ 2 c that just meets the requirements/avoids of the respective category 1 or 2, respectively.
Amongst all tests, cv with BGK collision expectedly showed the weakest stability. It turned out that the onset of instability is linked to the occurence of spurious transverse velocity oscillations that occur probably through the interactions of non-hydrodynamic modes with the grid transition interface. By including a filtering operation during the fine-to-coarse communication in the cv algorithm, the accuracy and stability range of cv BGK were enhanced, reaching corresponding values of its cm counterpart, where this filtering is granted per se by the interpolations necessary during coupling steps. However, further work is needed regarding the relation between the onset of instability, the occurence of spurious oscillations and how they can be partially suppressed by a restriction of fine populations when transferring them to the coarse grid.
With advanced collision operators, i.e., all except BGK, weaknesses of the cv grid transition scheme regarding stability and accuracy ranges could be overcome without the use of an explicit filtering step, extending both ranges towards the immediate vicinity of the cm scheme, both methods, however, with a moderate distance to the cc algorithm. A hint on the latter’s superior stability range was given by inserting the numerically obtained flow field into the analytical solution for the square duct flow and solving for the viscosity in order to quantify the relative amount of total viscosity of the system, averaged over the whole domain ν ¯ total / ν . The cc algorithm exhibits an increased amount of numerical dissipation, possibly linked to the implicit filtering of the complete distribution function during coalescence, instead of only the non-equilibrium part as is done in the cv filtering step and trilinear interpolation of the cm scheme. A deeper look into this phenomenon will be left open for future work.
The highest increase in stability was found independently of the specific grid layout, by use of the HRR collision operator. Even with a hybridization factor of σ = 0.98 , i.e., two percent of stress tensor reconstruction through central finite differences, an increase in stability range beyond the scope of our terminal ω ˜ 2 c > 1.99999 value, meaning no negative distribution functions were detected with this collision operator in our study up to a bulk Reynolds number of Re h , b = 3.470 × 10 5 . This leads to the final conclusion that, although peculiar features of the different grid transition schemes can influence the stability to a certain extent, the effect of the collision model is far more crucial. Similar conclusions regarding the impact of the collision model have been drawn by Astoul et al. [41]. The possibility of adjusting the amount of stabilizing numerical hyper-viscosity injection by the hybridization in the HRR model, was shown by evaluating ν ¯ total / ν as described above, for three distinct σ values. The influence of various hybridization parameter values on the stability range of the different grid transition schemes will be examined in the future.
This paper represents the first part of an investigation to better understand the relationships between different grid refinement algorithms in the lattice Boltzmann context and their peculiar properties. In a forthcoming study, this will be extended to turbulent flows and grid transitions of varying orientation relative to the flow. Moreover, an analysis regarding the properties of the various refinement schemes in aeroacoustic scenarios is planned.

Author Contributions

Conceptualization, methodology, validation, formal analysis, A.S. (Alexander Schukmann) and A.S. (Andreas Schneider); writing—Original draft preparation, visualization, A.S. (Alexander Schukmann); writing—Review and editing, A.S. (Andreas Schneider); investigation, A.S. (Alexander Schukmann) and V.H.; data curation, V.H.; supervision, project administration, M.B. All authors have read and agreed to the published version of the manuscript.

Funding

The article processing charge was funded by the Baden-Württemberg Ministry of Science, Research and Culture and Offenburg University in the funding programme Open Access Publishing.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data can be provided by the authors upon request.

Acknowledgments

All simulations in this study were carried out on the high performance cluster Elwetritsch of the University of Kaiserslautern-Landau which is part of the Alliance for High Performance Computing Rheinland-Pfalz (AHRP). We kindly acknowledge the support of the computing center.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
CFDComputational fluid dynamics
LBMLattice Boltzmann method
BGKBhatnagar–Gross–Krook
MRTMultiple relaxation time
REGRegularized model
NEUvon Neumann analysis
RRRecursive regularized
HRRHybrid-recursive regularized
HRR ψ Hybrid-recursive regularized model with cubic correction terms
LBELattice Boltzmann equation
SRTSingle relaxation time
RMSRoot mean square
cvCell-vertex
ccCell-centered
cmCombined
slSingle level

Appendix A. Multiple Relaxation Time Collision Model

Transformation matrix M
( 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ) 1 ( 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 ) ξ ˜ 2 ( 1 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 ) ξ ˜ 4 ( 0 1 1 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 ) ξ ˜ ( 0 2 2 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 ) ξ ˜ 3 ( 0 0 0 1 1 0 0 1 1 1 1 0 0 0 0 1 1 1 1 ) ξ ˜ ( 0 0 0 2 2 0 0 1 1 1 1 0 0 0 0 1 1 1 1 ) ξ ˜ 3 ( 0 0 0 0 0 1 1 0 0 0 0 1 1 1 1 1 1 1 1 ) ξ ˜ ( 0 0 0 0 0 2 2 0 0 0 0 1 1 1 1 1 1 1 1 ) ξ ˜ 3 ( 0 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 ) ξ ˜ 2 ( 0 2 2 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 ) ξ ˜ 4 ( 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 ) ξ ˜ 2 ( 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 ) ξ ˜ 4 ( 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 ) ξ ˜ 2 ( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 ) ξ ˜ 2 ( 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 ) ξ ˜ 2 ( 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 ) ξ ˜ 3 ( 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 ) ξ ˜ 3 ( 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 ) ξ ˜ 3
Equilibrium moments for the MRT collision model
m 0 0 = ρ m 1 0 = ρ u 2 m 3 0 = ρ u x m 5 0 = ρ u y m 7 0 = ρ u z m 9 0 = ρ 2 u x 2 u y 2 u z 2 m 11 0 = ρ u y 2 u z 2 m 13 0 = ρ u x u y m 14 0 = ρ u y u z m 15 0 = ρ u x u z m 2 0 = m 4 0 = m 6 0 = m 8 0 = m 10 0 = m 12 0 = m 16 0 = m 17 0 = m 18 0 = 0

Appendix B. Cubic Mach Correction For D3Q19 Lattice

To correct the O Ma 3 error in the stress tensor, resulting from insufficient quadrature order of the D3Q19 lattice, Feng et al. [46] proposed cubic Mach correction terms of the following form:
ψ i = w i 2 c s 4 [ H i , x x ( 2 ) x Ψ x x x + y Ψ x x y + z Ψ x x z + H i , y y ( 2 ) x Ψ x y y + y Ψ y y y + z Ψ y y z + H i , z z ( 2 ) x Ψ x z z + y Ψ y z z + z Ψ z z z + 2 H i , x y ( 2 ) x Ψ x x y + y Ψ x y y + z Ψ x y z + 2 H i , x z ( 2 ) x Ψ x x z + y Ψ x y z + z Ψ x z z + 2 H i , y z ( 2 ) x Ψ x y z + y Ψ y y z + z Ψ y z z ] ,
with deviation terms in the athermal model employed here, given as:
Ψ α α α = ρ u α 3 Ψ α α β = 1 2 ρ u β u γ 2 Ψ α β γ = ρ u α u β u γ .
Spatial derivatives are approximated using second-order centered finite differences [41]:
α Ψ α β γ Ψ α β γ x + e α Δ x Ψ α β γ x e α Δ x 2 Δ x .

Appendix C. Gradient-Based Velocity Interpolation Coefficients

In addition to the eight known base point velocities u B i , i = 0 7 , further constraints are needed in order to solve for the 33 unknown interpolation coefficients in Equation (59). For this purpose, we average strain-rates on each cell face and utilize these averaged values in second-order centered finite difference approximations of strain-rate gradients in the cell center. For example, given a cell nomenclature as depicted in Figure A1, x-direction gradients of strain-rates are discretized as:
Figure A1. Base point nomenclature for compact gradient-based velocity interpolation.
Figure A1. Base point nomenclature for compact gradient-based velocity interpolation.
Fluids 08 00103 g0a1
x S ¯ α β | 0 , 0 , 0 S ¯ α β Δ x 2 , 0 , 0 S ¯ α β + Δ x 2 , 0 , 0 Δ x ,
where S ¯ α β are cell-face-averaged strain rates, evaluated using local strain rates S α β , B i from second-order non-equilibrium moments computed at base points defining respective face corners:
S ¯ α β Δ x 2 , 0 , 0 = 1 4 S α β , B 0 + S α β , B 3 + S α β , B 4 + S α β , B 7
S ¯ α β + Δ x 2 , 0 , 0 = 1 4 S α β , B 1 + S α β , B 2 + S α β , B 5 + S α β , B 6 .
With this, a sufficient number of contraints is established and the compact gradient-based interpolation coefficients can be determined by solving the corresponding matrix-vector equation. The coefficients can be summarized as follows.
Zeroth-order x-direction coefficient:
a x , 000 = 1 8 u x , B 0 + u x , B 1 + u x , B 2 + u x , B 3 + u x , B 4 + u x , B 5 + u x , B 6 + u x , B 7 + 1 16 u y , B 0 u y , B 1 + u y , B 2 u y , B 3 + u y , B 4 u y , B 5 + u y , B 6 u y , B 7 + 1 16 u z , B 0 u z , B 1 u z , B 2 + u z , B 3 u z , B 4 + u z , B 5 + u z , B 6 u z , B 7 + Δ x 32 S x x , B 0 S x x , B 1 S x x , B 2 + S x x , B 3 + S x x , B 4 S x x , B 5 S x x , B 6 + S x x , B 7 + Δ x 16 S x y , B 0 + S x y , B 1 S x y , B 2 S x y , B 3 + S x y , B 4 + S x y , B 5 S x y , B 6 S x y , B 7 + Δ x 16 S x z , B 0 + S x z , B 1 + S x z , B 2 + S x z , B 3 S x z , B 4 S x z , B 5 S x z , B 6 S x z , B 7
First-order x-direction coefficient:
a x , 100 = 1 4 Δ x u x , B 0 u x , B 1 u x , B 2 + u x , B 3 + u x , B 4 u x , B 5 u x , B 6 + u x , B 7 a x , 010 = 1 4 Δ x u x , B 0 + u x , B 1 u x , B 2 u x , B 3 + u x , B 4 + u x , B 5 u x , B 6 u x , B 7 a x , 001 = 1 4 Δ x u x , B 0 + u x , B 1 + u x , B 2 + u x , B 3 u x , B 4 u x , B 5 u x , B 6 u x , B 7
Second-order x-direction coefficient:
a x , 200 = 1 8 Δ x S x x , B 0 S x x , B 1 S x x , B 2 + S x x , B 3 + S x x , B 4 S x x , B 5 S x x , B 6 + S x x , B 7 a x , 020 = 1 4 Δ x 2 u y , B 0 u y , B 1 + u y , B 2 u y , B 3 + u y , B 4 u y , B 5 + u y , B 6 u y , B 7 1 4 Δ x S x y , B 0 + S x y , B 1 S x y , B 2 S x y , B 3 + S x y , B 4 + S x y , B 5 S x y , B 6 S x y , B 7 a x , 002 = 1 4 Δ x 2 u z , B 0 u z , B 1 u z , B 2 + u z , B 3 u z , B 4 + u z , B 5 + u z , B 6 u z , B 7 1 4 Δ x S x z , B 0 + S x z , B 1 + S x z , B 2 + S x z , B 3 S x z , B 4 S x z , B 5 S x z , B 6 S x z , B 7 a x , 110 = 1 2 Δ x 2 u x , B 0 u x , B 1 + u x , B 2 u x , B 3 + u x , B 4 u x , B 5 + u x , B 6 u x , B 7 a x , 101 = 1 2 Δ x 2 u x , B 0 u x , B 1 u x , B 2 + u x , B 3 u x , B 4 + u x , B 5 + u x , B 6 u x , B 7 a x , 011 = 1 2 Δ x 2 u x , B 0 + u x , B 1 u x , B 2 u x , B 3 u x , B 4 u x , B 5 + u x , B 6 + u x , B 7
Third-order x-direction coefficient:
a x , 111 = 1 Δ x 3 u x , B 0 u x , B 1 + u x , B 2 u x , B 3 u x , B 4 + u x , B 5 u x , B 6 + u x , B 7
Zeroth-order y-direction coefficient:
a y , 000 = 1 8 u y , B 0 + u y , B 1 + u y , B 2 + u y , B 3 + u y , B 4 + u y , B 5 + u y , B 6 + u y , B 7 + 1 16 u x , B 0 u x , B 1 + u x , B 2 u x , B 3 + u x , B 4 u x , B 5 + u x , B 6 u x , B 7 + 1 16 u z , B 0 + u z , B 1 u z , B 2 u z , B 3 u z , B 4 u z , B 5 + u z , B 6 + u z , B 7 + Δ x 32 S y y , B 0 + S y y , B 1 S y y , B 2 S y y , B 3 + S y y , B 4 + S y y , B 5 S y y , B 6 S y y , B 7 + Δ x 16 S x y , B 0 S x y , B 1 S x y , B 2 + S x y , B 3 + S x y , B 4 S x y , B 5 S x y , B 6 + S x y , B 7 + Δ x 16 S y z , B 0 + S y z , B 1 + S y z , B 2 + S y z , B 3 S y z , B 4 S y z , B 5 S y z , B 6 S y z , B 7
First-order y-direction coefficient:
a y , 100 = 1 4 Δ x u y , B 0 u y , B 1 u y , B 2 + u y , B 3 + u y , B 4 u y , B 5 u y , B 6 + u y , B 7 a y , 010 = 1 4 Δ x u y , B 0 + u y , B 1 u y , B 2 u y , B 3 + u y , B 4 + u y , B 5 u y , B 6 u y , B 7 a y , 001 = 1 4 Δ x u y , B 0 + u y , B 1 + u y , B 2 + u y , B 3 u y , B 4 u y , B 5 u y , B 6 u y , B 7
Second-order y-direction coefficient:
a y , 200 = 1 4 Δ x 2 u x , B 0 u x , B 1 + u x , B 2 u x , B 3 + u x , B 4 u x , B 5 + u x , B 6 u x , B 7 1 4 Δ x S x y , B 0 S x y , B 1 S x y , B 2 + S x y , B 3 + S x y , B 4 S x y , B 5 S x y , B 6 + S x y , B 7 a y , 020 = 1 8 Δ x S y y , B 0 + S y y , B 1 S y y , B 2 S y y , B 3 + S y y , B 4 + S y y , B 5 S y y , B 6 S y y , B 7 a y , 002 = 1 4 Δ x 2 u z , B 0 + u z , B 1 u z , B 2 u z , B 3 u z , B 4 u z , B 5 + u z , B 6 + u z , B 7 1 4 Δ x S y z , B 0 + S y z , B 1 + S y z , B 2 + S y z , B 3 S y z , B 4 S y z , B 5 S y z , B 6 S y z , B 7 a y , 110 = 1 2 Δ x 2 u y , B 0 u y , B 1 + u y , B 2 u y , B 3 + u y , B 4 u y , B 5 + u y , B 6 u y , B 7 a y , 101 = 1 2 Δ x 2 u y , B 0 u y , B 1 u y , B 2 + u y , B 3 u y , B 4 + u y , B 5 + u y , B 6 u y , B 7 a y , 011 = 1 2 Δ x 2 u y , B 0 + u y , B 1 u y , B 2 u y , B 3 u y , B 4 u y , B 5 + u y , B 6 + u y , B 7
Third-order y-direction coefficient:
a y , 111 = 1 Δ x 3 u y , B 0 u y , B 1 + u y , B 2 u y , B 3 u y , B 4 + u y , B 5 u y , B 6 + u y , B 7
Zeroth-order z-direction coefficient:
a z , 000 = 1 8 u z , B 0 + u z , B 1 + u z , B 2 + u z , B 3 + u z , B 4 + u z , B 5 + u z , B 6 + u z , B 7 + 1 16 u x , B 0 u x , B 1 u x , B 2 + u x , B 3 u x , B 4 + u x , B 5 + u x , B 6 u x , B 7 + 1 16 u y , B 0 + u y , B 1 u y , B 2 u y , B 3 u y , B 4 u y , B 5 + u y , B 6 + u y , B 7 + Δ x 32 S z z , B 0 + S z z , B 1 + S z z , B 2 + S z z , B 3 S z z , B 4 S z z , B 5 S z z , B 6 S z z , B 7 + Δ x 16 S x z , B 0 S x z , B 1 S x z , B 2 + S x z , B 3 + S x z , B 4 S x z , B 5 S x z , B 6 + S x z , B 7 + Δ x 16 S y z , B 0 + S y z , B 1 S y z , B 2 S y z , B 3 + S y z , B 4 + S y z , B 5 S y z , B 6 S y z , B 7
First-order y-direction coefficient:
a z , 100 = 1 4 Δ x u z , B 0 u z , B 1 u z , B 2 + u z , B 3 + u z , B 4 u z , B 5 u z , B 6 + u z , B 7 a z , 010 = 1 4 Δ x u z , B 0 + u z , B 1 u z , B 2 u z , B 3 + u z , B 4 + u z , B 5 u z , B 6 u z , B 7 a z , 001 = 1 4 Δ x u z , B 0 + u z , B 1 + u z , B 2 + u z , B 3 u z , B 4 u z , B 5 u z , B 6 u z , B 7
Second-order z-direction coefficient:
a z , 200 = 1 4 Δ x 2 u x , B 0 u x , B 1 u x , B 2 + u x , B 3 u x , B 4 + u x , B 5 + u x , B 6 u x , B 7 1 4 Δ x S x z , B 0 S x z , B 1 S x z , B 2 + S x z , B 3 + S x z , B 4 S x z , B 5 S x z , B 6 + S x z , B 7 a z , 020 = 1 4 Δ x 2 u y , B 0 + u y , B 1 u y , B 2 u y , B 3 u y , B 4 u y , B 5 + u y , B 6 + u y , B 7 1 4 Δ x S y z , B 0 + S y z , B 1 S y z , B 2 S y z , B 3 + S y z , B 4 + S y z , B 5 S y z , B 6 S y z , B 7 a z , 002 = 1 8 Δ x S z z , B 0 + S z z , B 1 + S z z , B 2 + S z z , B 3 S z z , B 4 S z z , B 5 S z z , B 6 S z z , B 7 a z , 110 = 1 2 Δ x 2 u z , B 0 u z , B 1 + u z , B 2 u z , B 3 + u z , B 4 u z , B 5 + u z , B 6 u z , B 7 a z , 101 = 1 2 Δ x 2 u z , B 0 u z , B 1 u z , B 2 + u z , B 3 u z , B 4 + u z , B 5 + u z , B 6 u z , B 7 a z , 011 = 1 2 Δ x 2 u z , B 0 + u z , B 1 u z , B 2 u z , B 3 u z , B 4 u z , B 5 + u z , B 6 + u z , B 7
Third-order z-direction coefficient:
a z , 111 = 1 Δ x 3 u z , B 0 u z , B 1 + u z , B 2 u z , B 3 u z , B 4 + u z , B 5 u z , B 6 + u z , B 7

Appendix D. Normalized Square Duct Velocity Profiles

Figure A2. Normalized square duct velocity profiles for various tested models with ω ˜ c = 1.94990 ( Re h , b = 67.53 ).
Figure A2. Normalized square duct velocity profiles for various tested models with ω ˜ c = 1.94990 ( Re h , b = 67.53 ).
Fluids 08 00103 g0a2aFluids 08 00103 g0a2b

References

  1. Filippova, O.; Hänel, D. Grid Refinement for Lattice-BGK Models. J. Comput. Phys. 1998, 147, 219–228. [Google Scholar] [CrossRef]
  2. Dupuis, A.; Chopard, B. Theory and applications of an alternative lattice Boltzmann grid refinement algorithm. Phys. Rev. E 2003, 67, 066707. [Google Scholar] [CrossRef]
  3. Lagrava, D.; Malaspinas, O.; Latt, J.; Chopard, B. Advances in multi-domain lattice Boltzmann grid refinement. J. Comput. Phys. 2012, 231, 4808–4822. [Google Scholar] [CrossRef] [Green Version]
  4. Chen, H.; Teixeira, C.; Molvig, K. Realization of fluid boundary conditions via discrete Boltzmann dynamics. Int. J. Mod. Phys. C 1998, 9, 1281–1292. [Google Scholar] [CrossRef]
  5. Rohde, M.; Kandhai, D.; Derksen, J.J.; van den Akker, H.E.A. A generic, mass conservative local grid refinement technique for lattice-Boltzmann schemes. Int. J. Numer. Methods Fluids 2006, 51, 439–468. [Google Scholar] [CrossRef]
  6. Freitas, R.K.; Meinke, M.; Schröder, W. Lattice-Boltzmann Turbulence simulation via the Lattice-Boltzmann method on hierarchically refined meshes. In Proceedings of the European Conference on Computational Fluid Dynamics, Egmond aan Zee, The Netherlands, 5–8 September 2006; p. 12. [Google Scholar]
  7. Geier, M.; Greiner, A.; Korvink, J.G. Bubble functions for the lattice Boltzmann method and their application to grid refinement. Eur. Phys. J. Spec. Top. 2009, 171, 173–179. [Google Scholar] [CrossRef]
  8. Schönherr, M.; Kucher, K.; Geier, M.; Stiebler, M.; Freudiger, S.; Krafczyk, M. Multi-thread implementations of the lattice Boltzmann method on non-uniform grids for CPUs and GPUs. Comput. Math. Appl. 2011, 61, 3730–3743. [Google Scholar] [CrossRef] [Green Version]
  9. Sterling, J.D.; Chen, S. Stability Analysis of Lattice Boltzmann Methods. J. Comput. Phys. 1996, 123, 196–206. [Google Scholar] [CrossRef] [Green Version]
  10. Ricot, D.; Marié, S.; Sagaut, P.; Bailly, C. Lattice Boltzmann method with selective viscosity filter. J. Comput. Phys. 2009, 228, 4478–4490. [Google Scholar] [CrossRef]
  11. Gendre, F.; Ricot, D.; Fritz, G.; Sagaut, P. Grid refinement for aeroacoustics in the lattice Boltzmann method: A directional splitting approach. Phys. Rev. E 2017, 96, 023311. [Google Scholar] [CrossRef]
  12. Astoul, T.; Wissocq, G.; Boussuge, J.F.; Sengissen, A.; Sagaut, P. Analysis and reduction of spurious noise generated at grid refinement interfaces with the lattice Boltzmann method. J. Comput. Phys. 2020, 418, 109645. [Google Scholar] [CrossRef]
  13. Bhatnagar, P.L.; Gross, E.P.; Krook, M. A Model For Collison Processes in Gases. Phys. Rev. 1954, 94, 511–525. [Google Scholar] [CrossRef]
  14. Coveney, P.V.; Succi, S.; d’Humieres, D.; Ginzburg, I.; Krafczyk, M.; Lallemand, P.; Luo, L.S. Multiple-relaxation-time lattice Boltzmann models in three dimensions. Philos. Trans. R. Soc. London. Ser. Math. Phys. Eng. 2002, 360, 437–451. [Google Scholar]
  15. Malaspinas, O. Increasing stability and accuracy of the lattice Boltzmann scheme: Recursivity and regularization. arXiv 2015, arXiv:1505.06900. [Google Scholar]
  16. Jacob, J.; Malaspinas, O.; Sagaut, P. A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation. J. Turbul. 2018, 19, 1051–1076. [Google Scholar] [CrossRef] [Green Version]
  17. Wolf-Gladrow, D.A. Lattice-Gas Cellular Automata and Lattice Boltzmann Models: An Introduction; Springer: Berlin/Heidelberg, Germany, 2000; Volume 1725. [Google Scholar]
  18. Qian, Y.; d’Humieres, D.; Lallemand, P. Lattice BGK Models for Navier–Stokes Equation. Europhys. Lett. 1992, 17, 479–484. [Google Scholar] [CrossRef]
  19. He, X.; Luo, L.S. A priori derivation of the lattice Boltzmann equation. Phys. Rev. E 1997, 55, R6333–R6336. [Google Scholar] [CrossRef] [Green Version]
  20. Schneider, A. A Consistent Large Eddy Approach for Lattice Boltzmann Methods and its Application to Complex Flows. Ph.D. Thesis, Technical University of Kaiserslautern, Kaiserslautern, Germany, 2015. [Google Scholar]
  21. Chen, H.; Shan, X. Fundamental conditions for N-th order accurate lattice Boltzmann models. Phys. D Nonlinear Phenom. 2008, 237, 2003–2008. [Google Scholar] [CrossRef] [Green Version]
  22. Dellar, P.J. Incompressible limits of lattice Boltzmann equations using multiple relaxation times. J. Comput. Phys. 2014, 190, 351–370. [Google Scholar] [CrossRef] [Green Version]
  23. Ansumali, S.; Karlin, I.V. Single relaxation time model for entropic lattice Boltzmann methods. Phys. Rev. E 2002, 65, 056312. [Google Scholar] [CrossRef] [Green Version]
  24. Geier, M.C. Ab Initio Derivation of the Cascaded Lattice Boltzmann Automaton. Ph.D. Thesis, University of Freiburg, Freiburg, Germany, 2006. [Google Scholar]
  25. Geier, M.; Schönherr, M.; Pasquali, A.; Krafczyk, M. The cumulant lattice Boltzmann equation in threedimensions: Theory and validation. Comput. Math. Appl. 2015, 70, 507–547. [Google Scholar] [CrossRef] [Green Version]
  26. Krüger, T.; Kusumaatmaja, H.; Kuzmin, A.; Shardt, O.; Silva, G.; Viggen, E.M. The Lattice Boltzmann Method: Principles and Practice; Graduate Texts in Physics; Springer: Berlin/Heidelberg, Germany, 2017. [Google Scholar]
  27. Guo, Z.; Zheng, C.; Shi, B. Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys. Rev. E 2002, 65, 046308. [Google Scholar] [CrossRef]
  28. Hübner, T. A Monolithicm, Off-Lattice Approach to the Discrete Boltzmann Equation with Fast and Accurate Numerical Methods. Ph.D. Thesis, Technical Unversity of Dortmund, Dortmund, Germany, 2011. [Google Scholar]
  29. Shan, X.; Yuan, X.F.; Chen, H. Kinetic theory representation of hydrodynamics: A way beyond the Navier–Stokes equation. J. Fluid Mech. 2006, 550, 413–441. [Google Scholar] [CrossRef]
  30. Coreixas, C.; Chopard, B.; Latt, J. Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations. Phys. Rev. E 2019, 100, 033305. [Google Scholar] [CrossRef] [Green Version]
  31. Yoo, H.; Bahlali, M.L.; Favier, J.; Sagaut, P. A hybrid recursive regularized lattice Boltzmann model with overset grids for rotating geometries. Phys. Fluids 2021, 33, 057113. [Google Scholar] [CrossRef]
  32. Coreixas, C. High-Order Extension of the Recursive Regularized Lattice Boltzmann Method. Ph.D. Thesis, Université de Toulouse, Toulouse, France, 2018. [Google Scholar]
  33. Guo, Z.; Shu, C. Lattice Boltzmann Method and Its Application to Engineering; World Scientific: Singapore, 2013. [Google Scholar]
  34. Dellar, P.J. Nonhydrodynamic modes and a priori construction of shallow water lattice Boltzmann equations. Phys. Rev. E 2002, 65, 036309. [Google Scholar] [CrossRef] [Green Version]
  35. Lallemand, P.; Luo, L.S. Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance and stability. Phys. Rev. E 2000, 61, 6546–6562. [Google Scholar] [CrossRef] [Green Version]
  36. Brownlee, R.A.; Levesley, J.; Packwood, D.; Gorban, A.N. Add-ons for Lattice Boltzmann Methods: Regularization, Filtering and Limiters. Prog. Comput. Phys. 2013, 3, 31–52. [Google Scholar]
  37. Marié, S.; Ricot, D.; Sagaut, P. Comparison between lattice Boltzmann method and Navier–Stokes high order schemes for computational aeroacoustics. J. Comput. Phys. 2009, 228, 1056–1070. [Google Scholar] [CrossRef]
  38. Chapman, S.; Cowling, T.G. The Mathematical Theory of Non-Uniform Gases; Cambridge University Press: Cambridge, UK, 1953. [Google Scholar]
  39. Lätt, J.; Chopard, B. Lattice Boltzmann method with regularized pre-collision distribution functions. Math. Comput. Simul. 2006, 72, 165–168, Discrete Simulation of Fluid Dynamics in Complex Systems. [Google Scholar] [CrossRef]
  40. Zhang, R.; Shan, X.; Chen, H. Efficient kinetic method for fluid simulation beyond the Navier–Stokes equation. Phys. Rev. E 2006, 74, 046703. [Google Scholar] [CrossRef] [Green Version]
  41. Astoul, T.; Wissocq, G.; Boussuge, J.F.; Sengissen, A.; Sagaut, P. Lattice Boltzmann method for computational aeroacoustics on non-uniform meshes: A direct grid coupling approach. J. Comput. Phys. 2021, 447, 110667. [Google Scholar] [CrossRef]
  42. Li, Z.; Cao, W.; Touzé, D.L. On the coupling of a direct-forcing immersed boundary method and the regularized lattice Boltzmann method for fluid–structure interaction. Comput. Fluids 2019, 190, 470–484. [Google Scholar] [CrossRef]
  43. Guo, Z.; Zheng, C.; Shi, B.; Zhao, T.S. Thermal lattice Boltzmann equation for low Mach number flows: Decoupling model. Phys. Rev. E 2007, 75, 036704. [Google Scholar] [CrossRef]
  44. Renard, F.; Wissocq, G.; Boussuge, J.; Sagaut, P. A linear stability analysis of compressible hybrid lattice Boltzmann methods. arXiv 2020, arXiv:2006.08477. [Google Scholar] [CrossRef]
  45. Wissocq, G.; Coreixas, C.; Boussuge, J.F. Linear stability and isotropy properties of athermal regularized lattice Boltzmann methods. Phys. Rev. E 2020, 102, 053305. [Google Scholar] [CrossRef]
  46. Feng, Y.; Boivin, P.; Jacob, J.; Sagaut, P. Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows. J. Comput. Phys. 2019, 394, 82–99. [Google Scholar] [CrossRef]
  47. Astoul, T. Towards improved lattice Boltzmann aeroacoustic simulations with non-uniform grids: Application to landing gears noise prediction. Ph.D. Thesis, Aix-Marseille Université, Marseille, France, 2021. [Google Scholar]
  48. He, X.; Chen, S.; Doolen, D. A novel thermal model for the lattice Boltzmann method in incompressible limit. J. Comput. Phys. 1998, 146, 282–300. [Google Scholar] [CrossRef]
  49. Guo, Z.; Zheng, C.; Shi, B. Lattice Boltzmann equation with multiple effective relaxation times for gaseous microscale flow. Phys. Rev. E 2008, 77, 036707. [Google Scholar] [CrossRef]
  50. Lin, C.L.; Lai, Y.G. Lattice Boltzmann method on composite grids. Phys. Rev. E 2000, 62, 2219–2225. [Google Scholar] [CrossRef]
  51. Yu, D.; Mei, R.; Shyy, W. A multi-block lattice Boltzmann method for viscous fluid flows. Numer. Methods Fluids 2002, 39, 99–120. [Google Scholar] [CrossRef]
  52. Freudiger, S. Entwicklung eines parallelen, adaptiven, komponentenbasierten Strömungskerns für hierarchische Gitter auf Basis des Lattice-Boltzmann-Verfahrens. Ph.D. Thesis, Technische Universität Braunschweig, Braunschweig, Germany, 2009. [Google Scholar]
  53. Crouse, B. Lattice-Boltzmann Strömungssimulation auf Baumdatenstrukturen. Ph.D. Thesis, Technische Universität München, München, Germany, 2003. [Google Scholar]
  54. Shannon, C.E. Communication in the presence of noise. Proc. IRE 1949, 37, 10–21. [Google Scholar] [CrossRef]
  55. Quéméré, P.; Sagaut, P.; Couailler, C. A new multi-domain/multi-resolution method for large-eddy simulation. J. Numer. Methods Fluids 2001, 36, 391–416. [Google Scholar] [CrossRef]
  56. Touil, H.; Ricot, D.; Leveque, E. Direct and large-eddy simulation of turbulent flows on composite multi-resolution grids by the lattice Boltzmann method. J. Comput. Phys. 2013, 256, 220–233. [Google Scholar] [CrossRef] [Green Version]
  57. Guzik, S.M.; Weisgraber, T.H.; Colella, P.; Alder, B.J. Interpolation Methods and the Accuracy of Lattice-Boltzmann Mesh Refinement. J. Comput. Phys. 2012, 259, 461–487. [Google Scholar] [CrossRef]
  58. Kutscher, K.; Geier, M.; Krafczyk, M. Multiscale Simulation of Turbulent Flow interacting with Porous Media based on a massively parallel implementation of the Cumulant Lattice Boltzmann Method. Comput. Fluids 2018, 193, 103733. [Google Scholar] [CrossRef]
  59. Chen, H.; Filippova, O.; Hoch, J.; Molvig, K.; Shock, R.; Teixeira, C.; Zhang, R. Grid refinement in lattice Boltzmann methods based on volumetric formulation. Physica A 2006, 362, 158–167. [Google Scholar] [CrossRef]
  60. Qi, J.; Klimach, H.; Roller, S. Implementation of the compact interpolation within the octree based Lattice Boltzmann solver Musubi. Comput. Math. Appl. 2019, 78, 1131–1141. [Google Scholar] [CrossRef] [Green Version]
  61. Feuchter, C. Direct aeroacoustic simulation with a cumulant Lattice-Boltzmann model. Comput. Fluids 2021, 224, 104970. [Google Scholar] [CrossRef]
  62. Tölke, J.; Krafczyk, M. Second order interpolation of the flow field in the lattice Boltzmann method. Comput. Math. Appl. 2009, 58, 898–902. [Google Scholar] [CrossRef] [Green Version]
  63. Premnath, K.N.; Pattison, M.J.; Banerjee, S. Dynamic subgrid scale modeling of turbulent flows using lattice-Boltzmann method. Phys. A 2009, 388, 2640–2658. [Google Scholar]
  64. White, F.M. Viscous Fluid Flow; McGraw-Hill Series in Aeronautical and Aerospace Engineering; McGraw-Hill: New York, NY, USA, 1974. [Google Scholar]
  65. Tosun, I.; Uner, D.; Ozgen, C. Critical Reynolds number for Newtonian flow in rectangular ducts. Ind. Eng. Chem. Res. 1988, 27, 1955–1957. [Google Scholar]
  66. Hanks, R.W.; Ruo, H.C. Laminar-turbulent transition in ducts of rectangular cross section. Ind. Eng. Chem. Fundam. 1966, 5. [Google Scholar] [CrossRef]
Figure 1. Grid transition interface layouts.
Figure 1. Grid transition interface layouts.
Fluids 08 00103 g001
Figure 2. Discrete velocity sets for D2Q9 and D3Q19 lattices.
Figure 2. Discrete velocity sets for D2Q9 and D3Q19 lattices.
Fluids 08 00103 g002
Figure 3. Grid transition interface for two-dimensional cell-vertex layout, Fluids 08 00103 i015: regular fine node, Fluids 08 00103 i014: fine interface node with co-located coarse partner, Fluids 08 00103 i003: hanging middle node, Fluids 08 00103 i010: regular coarse node, Fluids 08 00103 i013: coarse interface node with co-located fine partner.
Figure 3. Grid transition interface for two-dimensional cell-vertex layout, Fluids 08 00103 i015: regular fine node, Fluids 08 00103 i014: fine interface node with co-located coarse partner, Fluids 08 00103 i003: hanging middle node, Fluids 08 00103 i010: regular coarse node, Fluids 08 00103 i013: coarse interface node with co-located fine partner.
Fluids 08 00103 g003
Figure 4. (Left): Grid transition interface for three-dimensional cell-vertex layout, Fluids 08 00103 i015: regular fine node, Fluids 08 00103 i014: co-located fine and coarse nodes at fine interface, Fluids 08 00103 i003: hanging middle node, Fluids 08 00103 i006: hanging center node, Fluids 08 00103 i010: regular coarse node, Fluids 08 00103 i013: co-located fine and coarse nodes at coarse interface. (Right): Stencil for hanging node interpolation.
Figure 4. (Left): Grid transition interface for three-dimensional cell-vertex layout, Fluids 08 00103 i015: regular fine node, Fluids 08 00103 i014: co-located fine and coarse nodes at fine interface, Fluids 08 00103 i003: hanging middle node, Fluids 08 00103 i006: hanging center node, Fluids 08 00103 i010: regular coarse node, Fluids 08 00103 i013: co-located fine and coarse nodes at coarse interface. (Right): Stencil for hanging node interpolation.
Fluids 08 00103 g004
Figure 5. Grid transition interface for two-dimensional cell-centered layout, Fluids 08 00103 i015: Regular fine node, Fluids 08 00103 i014: Fine interface node, Fluids 08 00103 i010: Regular coarse node, Fluids 08 00103 i013: Coarse interface node, inward pointing arrows { ξ 2 , ξ 5 , ξ 6 } represent directions of pre-collision populations averaged and assigned to Fluids 08 00103 i013 during coalescence, outward pointing arrows { ξ 0 , ξ 4 , ξ 7 } represent directions of post-collision populations distributed among Fluids 08 00103 i014 during explosion. Solid arrows show populations provided, dashed arrows populations received.
Figure 5. Grid transition interface for two-dimensional cell-centered layout, Fluids 08 00103 i015: Regular fine node, Fluids 08 00103 i014: Fine interface node, Fluids 08 00103 i010: Regular coarse node, Fluids 08 00103 i013: Coarse interface node, inward pointing arrows { ξ 2 , ξ 5 , ξ 6 } represent directions of pre-collision populations averaged and assigned to Fluids 08 00103 i013 during coalescence, outward pointing arrows { ξ 0 , ξ 4 , ξ 7 } represent directions of post-collision populations distributed among Fluids 08 00103 i014 during explosion. Solid arrows show populations provided, dashed arrows populations received.
Fluids 08 00103 g005
Figure 6. Grid transition interface for two-dimensional combined grid layout, Fluids 08 00103 i015: Regular fine node, Fluids 08 00103 i014: Fine interface node, Fluids 08 00103 i010: Regular coarse node, Fluids 08 00103 i013: Coarse interface node, light-red area indicates validity range of interpolation functions, inward pointing arrows represent pre-collision distribution functions reconstructed on both grids during c f and f c coupling according to Equations (60) and (61), respectively.
Figure 6. Grid transition interface for two-dimensional combined grid layout, Fluids 08 00103 i015: Regular fine node, Fluids 08 00103 i014: Fine interface node, Fluids 08 00103 i010: Regular coarse node, Fluids 08 00103 i013: Coarse interface node, light-red area indicates validity range of interpolation functions, inward pointing arrows represent pre-collision distribution functions reconstructed on both grids during c f and f c coupling according to Equations (60) and (61), respectively.
Fluids 08 00103 g006
Figure 7. Grid transition interface for two-dimensional cell-centered layout including coarse ghost nodes Fluids 08 00103 i017 needed for central-difference approximation of velocity gradients during HRR collision at coarse interface nodes Fluids 08 00103 i013, light-red area indicates coarse ghost cell from which fine pre-collision populations, represented as arrows, are gathered and averaged during fictitious coalescence according to Equation (55).
Figure 7. Grid transition interface for two-dimensional cell-centered layout including coarse ghost nodes Fluids 08 00103 i017 needed for central-difference approximation of velocity gradients during HRR collision at coarse interface nodes Fluids 08 00103 i013, light-red area indicates coarse ghost cell from which fine pre-collision populations, represented as arrows, are gathered and averaged during fictitious coalescence according to Equation (55).
Fluids 08 00103 g007
Figure 8. Grid transition interface for two-dimensional combined grid layout including coarse Fluids 08 00103 i017 and fine ghost nodes Fluids 08 00103 i018 needed for central-difference approximation of velocity gradients during HRR collision at coarse Fluids 08 00103 i013 and fine interface nodes Fluids 08 00103 i014, respectively, light-red area indicates validity range of velocity interpolation function (59).
Figure 8. Grid transition interface for two-dimensional combined grid layout including coarse Fluids 08 00103 i017 and fine ghost nodes Fluids 08 00103 i018 needed for central-difference approximation of velocity gradients during HRR collision at coarse Fluids 08 00103 i013 and fine interface nodes Fluids 08 00103 i014, respectively, light-red area indicates validity range of velocity interpolation function (59).
Fluids 08 00103 g008
Figure 9. Schematic representation of the calculation domain for the square duct test case.
Figure 9. Schematic representation of the calculation domain for the square duct test case.
Fluids 08 00103 g009
Figure 10. Normalized square duct velocity profiles for HRR models with σ = 0.98 and ω ˜ c = 1.94990 ( Re h , b = 67.53 ).
Figure 10. Normalized square duct velocity profiles for HRR models with σ = 0.98 and ω ˜ c = 1.94990 ( Re h , b = 67.53 ).
Fluids 08 00103 g010
Figure 11. Grid convergence for HRR collision model and all grid layouts with ω ˜ c = 1.94990 ( Re h , b = 67.53 ): Sl-4 (Fluids 08 00103 i024), cv-4 (Fluids 08 00103 i022), cv-7 (Fluids 08 00103 i023), cc-4 (Fluids 08 00103 i021), cc-6 (Fluids 08 00103 i019), cm-4 (Fluids 08 00103 i020), slope -1 (—), slope -2 (- -).
Figure 11. Grid convergence for HRR collision model and all grid layouts with ω ˜ c = 1.94990 ( Re h , b = 67.53 ): Sl-4 (Fluids 08 00103 i024), cv-4 (Fluids 08 00103 i022), cv-7 (Fluids 08 00103 i023), cc-4 (Fluids 08 00103 i021), cc-6 (Fluids 08 00103 i019), cm-4 (Fluids 08 00103 i020), slope -1 (—), slope -2 (- -).
Fluids 08 00103 g011
Figure 12. BGK accuracy and stability range comparison for different grid layouts and cv f c filters.
Figure 12. BGK accuracy and stability range comparison for different grid layouts and cv f c filters.
Fluids 08 00103 g012
Figure 13. Spurious transverse velocity oscillations with temporally increasing amplitude for cv-0 with ω ˜ 2 c = 1.98317 ( Re h , b = 205.04 ).
Figure 13. Spurious transverse velocity oscillations with temporally increasing amplitude for cv-0 with ω ˜ 2 c = 1.98317 ( Re h , b = 205.04 ).
Fluids 08 00103 g013
Figure 14. MRT accuracy and stability range comparison for different grid layouts.
Figure 14. MRT accuracy and stability range comparison for different grid layouts.
Fluids 08 00103 g014
Figure 15. RR/HRR accuracy and stability range comparison for different grid layouts.
Figure 15. RR/HRR accuracy and stability range comparison for different grid layouts.
Fluids 08 00103 g015
Figure 16. Accuracy and stability range comparison for all collision models with single level and cell-vertex grid.
Figure 16. Accuracy and stability range comparison for all collision models with single level and cell-vertex grid.
Fluids 08 00103 g016
Figure 17. Comparison of ν ¯ total / ν for all collision models and grid layouts with ω ˜ c = 1.94990 ( Re h , b = 67.53 ).
Figure 17. Comparison of ν ¯ total / ν for all collision models and grid layouts with ω ˜ c = 1.94990 ( Re h , b = 67.53 ).
Fluids 08 00103 g017
Figure 18. Accuracy and stability range comparison for all collision models with cell-centered and combined grid.
Figure 18. Accuracy and stability range comparison for all collision models with cell-centered and combined grid.
Fluids 08 00103 g018
Table 1. Summary of simulation parameters. Top half: Constant parameters, bottom half: Changing parameters at lower and upper limit of investigated Reynolds number range as well as critical flow regime.
Table 1. Summary of simulation parameters. Top half: Constant parameters, bottom half: Changing parameters at lower and upper limit of investigated Reynolds number range as well as critical flow regime.
Grid resolution 2 h / Δ x c 10Spacing Δ x c in m 2 × 10 4
Density ρ in k g   m 3 998.2 Kin. viscosity ν in m 2   s 1 1 × 10 6
Lattice Mach number Ma0.1Refined coarse cells at walls3
Re b , h , min Re b , h , crit Re b , h , max
Flow acceleration a x in m   s 2 2.075 × 10 1 5.140 6.329 1.066 × 10 3
Time step Δ t c in s 1.713 × 10 4 1.449 × 10 5 1.177 × 10 5 3.333 × 10 8
Collision frequency ω ˜ c 1.94990 1.99566 1.99647 1.99999
Bulk Reynolds number Re b , h 6.753 × 10 1 1.673 × 10 3 2.060 × 10 3 3.470 × 10 5
Table 2. Overview of variants examined in the numerical study, sl: Single level, cv: Cell-vertex, cc: Cell-centered, cm: Combined, all cv and cc cases above the light grey midline are computed without f c restriction and with uniform explosion, respectively.
Table 2. Overview of variants examined in the numerical study, sl: Single level, cv: Cell-vertex, cc: Cell-centered, cm: Combined, all cv and cc cases above the light grey midline are computed without f c restriction and with uniform explosion, respectively.
Collision Model Ω slcvcccm
BGKsl-0cv-0cc-0cm-0
MRT-REGsl-1cv-1cc-1cm-1
MRT-NEUsl-2cv-2cc-2cm-2
RRsl-3cv-3cc-3cm-3
HRRsl-4cv-4cc-4cm-4
BGK, f c filter Equation (50)cv-5
BGK, f c filter Equation (51)cv-6
BGK, linear explosion Equation (56)cc-5
HRR, linear explosion Equation (56)cc-6
HRR ψ cv-7
Table 3. Summary of grid parameters for convergence study with ω ˜ c = 1.94990 ( Re h , b = 67.53 ).
Table 3. Summary of grid parameters for convergence study with ω ˜ c = 1.94990 ( Re h , b = 67.53 ).
Grid resolution 2 h / Δ x f 204080
Spacing Δ x f in m 1 × 10 4 5 × 10 5 2.5 × 10 5
Time step Δ t f in s 4.28 × 10 5 1.07 × 10 5 5.35 × 10 6
Lattice Mach number Ma0.10.050.025
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Schukmann, A.; Schneider, A.; Haas, V.; Böhle, M. Analysis of Hierarchical Grid Refinement Techniques for the Lattice Boltzmann Method by Numerical Experiments. Fluids 2023, 8, 103. https://doi.org/10.3390/fluids8030103

AMA Style

Schukmann A, Schneider A, Haas V, Böhle M. Analysis of Hierarchical Grid Refinement Techniques for the Lattice Boltzmann Method by Numerical Experiments. Fluids. 2023; 8(3):103. https://doi.org/10.3390/fluids8030103

Chicago/Turabian Style

Schukmann, Alexander, Andreas Schneider, Viktor Haas, and Martin Böhle. 2023. "Analysis of Hierarchical Grid Refinement Techniques for the Lattice Boltzmann Method by Numerical Experiments" Fluids 8, no. 3: 103. https://doi.org/10.3390/fluids8030103

APA Style

Schukmann, A., Schneider, A., Haas, V., & Böhle, M. (2023). Analysis of Hierarchical Grid Refinement Techniques for the Lattice Boltzmann Method by Numerical Experiments. Fluids, 8(3), 103. https://doi.org/10.3390/fluids8030103

Article Metrics

Back to TopTop