Steady-state two-relaxation-time lattice Boltzmann formulation for transport and flow, closed with the compact multi-reflection boundary and interface-conjugate schemes
Introduction
The LBM, Lattice Boltzmann method [57], [58], destines for the intensive computations of the single/multi-phase flow and transport in a complex geometry, e.g. [10], [53], [76], [77], [98], [2], [92], [74], [87], [25], [4], [119], due to the extreme simplicity of its explicit-time marching bulk algorithm, local mass-conserving boundary handling and implicit-interface tracking. However, the porous media simulations are often intended for stationary velocity and phase field distributions, which are attainable after a very long number of computational steps. Hence, several strategies have been proposed to improve this issue. On the one side, the iterative momentum relaxation techniques [71], [68], [69], that adjust the driven body forcing until its complete cancellation by friction, enables the time steps reduction by half using the single-relaxation-times BGK model [97]. On the other side, much more radically, Verberg & Ladd [110] reduce the BGK Stokes-flow model with to linear algebraic system with respect to the macroscopic variables alone, pressure and momentum, and then gain one to two orders of magnitude in step reduction with the bi-conjugate matrix inversion in a random arrays of spheres, for about longer single-node update. The subsequent LBM “matrix” formulations concentrate their efforts on the efficiency of the linear and non-linear solvers, e.g. [108], [109], [88], [80]; a recent work [84] reports an enhanced efficiency of the non-linear, implicit, finite-volume, steady-state BGK algorithm for high Reynolds number simulations. Alternatively, the “preconditioned” equilibrium formulations [54], [66], [96], [89] aim at introducing of small adjustable prefactors in front of the time-derivatives in the modelled macroscopic equations. An independent scaling of the mass and diffusive-flux variables also accelerates the advection-diffusion schemes [33]; however their gain in time steps is counterbalanced by a reduction in stable velocity amplitude [51]. In flow schemes, the number of steps to steady-state depends on the lattice value of the kinematic viscosity but its optimal choice is not obvious and depends upon the convergence criteria [106], [70], irregularity of structure [106], the penetrability of porous inclusions [44] and the image size [87]. Finally, the “squared sound velocity” , relating pressure to density, is free tunable in linear flow computations: the number of steps and the convergence dependency upon decreases as increases [106] and diminishes, in accord with the “preconditioned” techniques.
It should be understood that the free tuning of the lattice transport coefficients, like kinematic viscosity, is only admissible when the numerical scheme respects the dimensionless group of the considered problem, basically formed from the Darcy, Péclet, Reynolds, Froude numbers, the transport coefficient contrasts and aspect ratios. The two-relaxation-times operator [34], [37], [38] accomplishes this requirement by keeping the specific freely-adjustable combination of its two collision rates at a fixed value, when the transport coefficient varies linearly with one of the two, [31], [65]. The control has been proven [65] through the generic form of the steady-state recurrence solutions, also there extended to more complicated operators. Moreover, the multiple-relaxation-times collisions [64] may achieve the same objective by keeping fixed all combinations of the symmetric and anti-symmetric mode rates, [38], [65], [70]; this property is not available either with the or regularized operators with fixed kinetic “ghost” rates [58], [32], [16], [86], where depends, respectively, quadratically and linearly on the transport coefficient. On the other side, it has been shown [37], [51] that the steady-state implicit-interface tracking can be identically reformulated through the anti-bounce back [34], [72] and the bounce-back, and interface-conjugate directional relations. The steady-state solutions of these boundary rules share the bulk parametrization [65]; consequently, the same values locate exactly the mid-grid straight/diagonal boundary and interface on the piece-wise parabolic, velocity and diffusion scalar, profiles, [28], [29], [37], [18], [49]. In this context, several numerical studies have performed consistent measurements of the transport coefficients at fixed , and examined their impact, including the prediction of the permeability and drag in homogeneous [106], [70], [8] or heterogeneous [44], [100] structures, effective permeability in non-Newtonian fluids [107], effective diffusivity [26], [11] and thermal conductivity of the composite materials [19]. On the theoretical perspective, impacts the truncation [42], [24], the numerical and Taylor dispersion [46], the high-order moments [48], the asymptotic convergence rate: quite uniformly in regular and random sphere packings [70], and stability [41], [75], [42].
The von Neumann stability analysis [41] reveals one singular choice: , called the “optimal” subclass, where the non-equilibrium becomes fully expressed through the directional central-difference operators of the equilibrium alone. The advanced stability is confirmed numerically through the Taylor-dispersion modelling in a broad Péclet range [46], but also in high Reynolds number dipole-wall collisions [90], [91]. The also distinguishes itself for an accurate modeling of the Taylor dispersion [46], effective diffusivity [11] and liquid/solid flux exchange [63], [120], [85]. The steady-state recurrence equations [65], [40] also allow to find exact solutions of the scheme in grid-aligned geometry, e.g. discrete-exponential, heterogeneous Brinkman channel flow, [44], [47]. The initial suspicion that the exhibits spurious Knudsen-type accommodation [17] in the grid-inclined channels [30] was later algebraically expressed [65] and quantified (i) with the Navier–Stokes steady-state channel flow [39]; and (ii), for tangential bounce-back moments impact in the transient Gaussian evolution along the flat wall [45], [49]. The underlying discrete-exponential mechanism is beyond the scope of the Chapman-Enskog or perturbative analysis, but its presence is not negligible, as it degrades the accuracy to the first order in inclined [30], [99], [102] and spherical [70] geometries. The steady-state linear model excites the whole, equilibrium and non-equilibrium, perturbation even on diagonally-oriented implicit interface, when the advection and/or diffusion equilibrium weights are freely rotated with respect to the lattice, and the scalar field evolves along all interface-cut lattice directions, [52].
The symbolic solution of the recurrence equations becomes tedious in grid-inclined slabs and motivates our interest for a robust steady-state analysis tool, which is able to handle any boundary and interface modifications. This seemingly difficult task turns out to be relatively simple and perhaps of much broader, numerical interest. The task of this work is hence threefold: (i) to formulate the system of the steady-state bulk equations; (ii) to extend and recast the multi-reflection boundary rules, and (iii), to adapt them with the directional interface-conjugate equations. Every couple of bulk equations are simply the aforementioned and , linear directional combinations of the equilibrium (macroscopic) and non-equilibrium variables, expressed with the half discrete velocity set owing to the symmetry argument; their form is the same for any lattice; it allows for non-uniform or discontinuous external sources and collision rates. The directional equations are complemented with the local mass and, in flow problem, momentum conservation steady-state solvability conditions. The implicit-interface tracking is matched automatically; otherwise, a couple of two directional equations are to be replaced by two interface-conjugate conditions, formulated in terms of the bulk unknowns. We use this opportunity to systematize and further generalize the Dirichlet approach for velocity [31], [38] and scalar fields [35], [38], [81], but also to extend the for particular pressure-stress, advection-diffusive flux, diffusive flux and Robin conditions.
By construction, the rule is equivalent to directional Taylor expansion along the cut-link, from the grid boundary node to an arbitrarily-shaped surface. Hereafter, we call “linear” those boundary techniques which only account for the first-order gradients in the reconstruction of the wall-incoming populations, e.g. with finite-differences [105], [20], interpolations [21], [9], [116], [16], three-population “linear” multi-reflection [35], [38], [81], or locally, through the non-equilibrium projections [32], [67], [86], [117]. Although the linear rules neglect the “kinetic” non-equilibrium component, their effective accuracy depends on the “ghost” collision rates, very noticeably in slow (porous) flow [31], [39] and diffusion-dominate problems, at least. The parabolic rules [31], [35], [38] account for the second-order derivatives through the third-order accurate Chapman-Enskog population expansion [28], [31], [35] or the recurrence solution form [38]. They attain then a quasi-analytical quality of velocity field and notably smooth pressure fluctuations on modest grids, e.g. in low Reynolds number flow around circular and spherical surface [31], [94], [70]; the parabolic rules also surpasses the partially-saturated or immersed boundary approaches in porous media simulations [94], [16], [95], [14]. The is further adapted for moving bodies [31], free-interface pressure-stress condition [7], spatial resistance variation in Brinkman porous models [102] and the finite Knudsen number slip-flows [103], [104]. Our particular emphasis will be put on the parabolic rules because they are able to locate exactly grid-rotated parabolic profiles, unreachable with the , and linear rules. These solutions were originally achieved with the local but parabolic, boundary concept LSOB [30]; its linear counterpart [117] is recently revived, and the parabolic LSOB is now extended [101] for smooth shaped-surface and acute angles. Finally, in this topic, it is worth mentioning that the LSOB and the moment-based boundary approach [93], [6] developed [90], [73], [91] for on-grid straight boundaries, are more suitable than the for independent prescription of the normal and tangential stress conditions [32].
Beyond of the formal accuracy, our main concern is the exact steady-state parametrization of the boundary rules by the governing physical numbers as, otherwise, the numerical dimensionless solutions and their relative errors will depend on the lattice values assigned for physical coefficients, such as the kinematic viscosity or molecular diffusion. The sufficient parametrization conditions [38] are expressed through the coefficients, and we extend them for all type schemes developed in the present work. In particular, we show that the linear Dirichlet scalar family [81], [82], [83], referred to as hereafter, and its local equivalent [118], are not parametrized by the Péclet number, except for their particular coefficient sets. In the past, linear Dirichlet flow schemes [9], [116] have been improved for parametrization within the infinite linear family [38], [39]. In the same fashion, we extend the to two infinite parametrized families, called and [simultaneously extending the pressure-linear scheme [38], [39]], and examine their respective accuracy for diffusion-and advection-dominant regimes. The compact steady-state formulation is recasted through the bulk unknowns alone; it is “one-node” for the linear schemes and “two-node” for the parabolic ones; hence, the computationally-efficient directional structure [87], [4] remains preserved. We also postulate and satisfy the “portability” conditions on the parabolic coefficients, from the to the principal, flow and , collisions, and outline their heuristic, adjustable stability bounds.
The interface-conjugate is built for two linear problems exemplified in arbitrarily-rotated slabs: the well known benchmark of the two-phase Stokes flow with the flat interface, e.g. [29], [37], [56], [122], [79], and the with the discontinuous mass-source and diffusion coefficient due to heterogeneous porosity [111], [50], [51], which is typical boundary problem found in the heat and mass transfer applications, e.g. [112], [115], [82], [55], [19], [78], [61], [60]. Our transient interface-conjugate update [] is inspired by the “decoupled interface treatment” [82], [55]; the steady-state interface-conjugate equations [] then replace every couple of bulk equations on interface cut-links. The is able to model exactly the piece-wise parabolic rotated velocity profile thanks to parabolic, Dirichlet velocity and pressure-stress, rules. In turn, the exact piece-wise parabolic solutions are closed by a broad panel of the interface-jump and boundary conditions, in the presence of grid-rotated interface-parallel advective velocity. A particular emphasis is placed on the “opposite” problem [50], [51] of a constant, interface-perpendicular Darcy velocity through the heterogeneous blocks in series; this tough benchmark has not been yet examined with the interface-conjugate, or sometimes is forbidden with it, e.g. [59].
Finally, we tackle another challenging problem. Basically, the explicit interface-conjugate algorithms maintain the flux balance only on the straight midway interface, e.g. [36], [82], while the grid-shifted boundary and interface, but also inlet/outlet conditions, e.g. [87], most often do not guarantee it exactly. Although the mass deficit can be recompensed, e.g. locally through the immobile population, it has been demonstrated [30] that exact boundary schemes for grid-inclined parabolic flow conserve the outgoing mass over the whole periodic segment; a harmful effect of the artificial mass conservation has been numerically demonstrated for flow field [1], [113] and scalar solutions with the Neumann flux schemes [23]. Typically, when the Dirichlet velocity condition does not preserve the global mass, the stationary momentum solution copes with the constant-rate mass flux over the computational domain; in other words, the immobile population is not at equilibrium and hence, an additive pressure constant will continuously vary in time [31]. We will propose and approve a simple, uniform modification of the local mass-conservation solvability condition able to mimics “quasi-stationary” solutions.
The paper is organized as follows. Section 2 summarizes the transient bulk update with external sources and introduces its steady-state counterpart. Section 3 analyses the and , defines the “linear” and “parabolic” closure relations, recalls the multi-reflection concept, recasts it in steady-state compact form and builds the sufficient parametrization conditions for Dirichlet, flux and shear-stress schemes. Section 4 extends the Dirichlet velocity class , constructs and applies the interface-conjugate and schemes in a two-phase stratified grid-rotated flow. Section 5 extends the “linear” and “parabolic” Dirichlet-scalar families and builds the particular, directional advective-diffusive flux, diffusive-flux and Robin boundary schemes. Section 6 extends the and for the problem with discontinuous collision rates, mass-sources and interface jumps. Section 7 examines the existence and accuracy of the stationary solutions in the interface-parallel and interface-perpendicular flow. Section 8 concludes the paper. Appendix A provides details on the construction of the schemes and tabulates their coefficients. Appendix B outlines the steady-state numerical algorithm with the boundary and interface-conjugate closure.
Section snippets
Transient (standard) and steady-state algorithms
The is operated on the cuboid computational grid within the penetrable -dimensional domain composed from grid points; a regular mesh-size and time step are set equal to one lattice unit and one iteration, respectively; the characteristic size is . The neighboring nodes are interconnected by the discrete velocity set : it consists of zero-vector and vectors ; the one-half of the non-zero velocity vectors is numbered with the positive numbers , their
Transient and steady-state closure relations
The population update in Eq. (3) is closed by the boundary rules for the populations incoming from the outside of the computational domain .
Definition I Any node with at least one wall-cut velocity such that , , is referred to as a boundary node; the set of the wall-cut velocities is denoted by .
Consider a boundary node and a wall-cut velocity , . The transient update with Eq. (3) computes with the boundary rule. In turn,
Flow schemes
The “linear” and “parabolic”, Stokes/Navier–Stokes flow families for the Dirichlet velocity condition have been introduced [31], generalized [38], [39] and extended for the Brinkman flow [102]. Namely, an infinite three-population “linear” family is correctly parametrized [38], [39], [70], and it includes the [9] and [116] schemes improved for this property with the help of the suitable post-collision correction . The member with is called ; is most
The Dirichlet scalar and Neumann flux schemes
We first extend and parametrize the Dirichlet scalar families [34], [38], [81]: the “linear” and “parabolic” is constructed in Sections 5.1 and 5.2, respectively. Section 5.3 re-builds the Neumann flux scheme [81] and its novel extension . Section 5.4 introduces the diffusion-flux Neumann families . Section 5.5 combines the and for mixed (Robin) family.
Interface-conjugate schemes
We construct the transient [] and steady-state [] interface-conjugate schemes for the heterogeneous advection-diffusion equation with the piece-wise continuous diffusion coefficients and sources. The interface-continuity or jump conditions, prescribed for the scalar and flux variables, are respectively based on the and schemes. They extend the “decoupled interface treatment” [82], [55] under specific conditions imposed on the tangential advection-diffusive flux. In
Heterogeneous in the interface-parallel and interface-normal flow
We consider steady-state solutions to Eq. (93) in the inclined stratified channel , , using the rotated frame from Eq. (69). The geometry is set periodic along , the interface is parallel with the axis and placed arbitrarily with respect to the grid at . Eq. (93) is modelled in the form:
The mass source is set velocity-dependent:
Concluding remarks
We proposed the steady-state formulation sustained with the compact boundary and interface closure equations; the method is expressed through the equilibrium and post-collision non-equilibrium components in the half discrete velocity space. The underlying matrix structure is mainly directional; in-node variables are only inter-connected through the conservation solvability constraints. When the stationary solution exists, it is the same with the transient and steady-state approaches,
Conflict of interest
The authors declare no conflict of interest.
Declaration of Competing Interest
The authors report no declarations of interest.
Acknowledgments
The author thanks Dr. Goncalo Silva for his critical reading of the manuscript.
References (116)
- et al.
A scalable multiphysics algorithm for massively parallel direct numerical simulations of electrophoretic motion
J. Comput. Sci.
(2018) - et al.
Boundary conditions for free interfaces with the lattice Boltzmann method
J. Comput. Phys.
(2015) - et al.
Drag correlation for dilute and moderately dense fluid-particle systems using the lattice Boltzmann method
Int. J. Multiph. Flow
(2015) - et al.
A comparative study on the lattice Boltzmann models for predicting effective diffusivity of porous media
Int. J. Heat Mass Transf.
(2016) - et al.
A Knudsen layer theory
Physica D
(1991) - et al.
Application and accuracy issues of TRT lattice Boltzmann method for solving elliptic PDEs commonly encountered in heat transfer and fluid flow problems
Int. J. Therm. Sci.
(2016) - et al.
Grad's approximation for moving and stationary walls in entropic lattice Boltzmann siimulations
J. Comput. Phys.
(2015) - et al.
Lattice Boltzmann simulation of gas-particle glow in filters
Comput. Fluids
(1997) - et al.
Parametrization of the cumulant lattice Boltzmann method for fourth order accurate diffusion part I, Derivation and validation
J. Comput. Phys.
(2017) - et al.
Parametrization of the cumulant lattice Boltzmann method for fourth order accurate diffusion part II, Application to flow around a sphere at drag crisis. M. Geier
J. Comput. Phys.
(2017)
Lattice Boltzmann model for free-surface flow and its application to filling process in casting
J. Comput. Phys.
Variably saturated flow described with the anisotropic Lattice Boltzmann methods
J. Comput. Fluids
Equilibrium-type and Link-type lattice Boltzmann models for generic advection and anisotropic-dispersion equation
Adv. Water Res.
Generic boundary conditions for Lattice Boltzmann models and their application to advection and anisotropic-dispersion equations
Adv. Water. Res.
Lattice Boltzmann and analytical modeling of flow processes in anisotropic and heterogeneous stratified aquifers
Adv. Water Resour.
Multiple anisotropic collisions for advection-diffusion lattice Boltzmann schemes
Adv. Water Res.
Local boundary reflections in lattice Boltzmann schemes: spurious boundary layers and their impact on the velocity, diffusion and dispersion
Comptes Rendus Mecanique
Truncation effect on Taylor-Aris dispersion in lattice Boltzmann schemes: accuracy towards stability
J. Comput. Phys.
Comment on “An improved gray lattice Boltzmann model for simulating fluid flow in multi-scale porous media”: intrinsic links between LBE Brinkman schemes
Adv. Water Resour.
Determination of the diffusivity, dispersion, skewness and kurtosis in heterogeneous porous flow. Part I: Analytical solutions with the extended method of moments
Adv. Water. Res.
Determination of the diffusivity, dispersion, skewness and kurtosis in heterogeneous porous flow. Part II: Lattice Boltzmann schemes with implicit interface
Adv. Water Res.
Lattice Boltzmann method for conjugate heat and mass transfer with interfacial jump conditions
Int. J. Heat Mass Transf.
Second-order curved interface treatments of the lattice Boltzmann method for convection-diffusion equations with conjugate interfacial conditions
Comput. Fluids
Second-order curved boundary treatments of the lattice Boltzmann method for convection-diffusion equations
J. Comput. Phys.
Phase interface effects in the total enthalpy-based lattice Boltzmann model for solid-liquid phase change
J. Comput. Phys.
Viscosity independent numerical errors for Lattice Boltzmann models: from recurrence equations to “magic” collision numbers
Comput. Math. Appl.
Implementation aspects of 3d lattice-bgk: boundaries, accuracy and a new fast relaxation method
J. Comput. Phys.
Iterative momentum relaxation for fast lattice-Boltzmann simulations
Future Gener. Comput. Syst.
A role of the kinetic parameter on the stability of two-relaxation-times advection-diffusion lattice Boltzmann scheme
Comput. Math. Appl.
Numerical evaluation of two recoloring operators for an immiscible two-phase flow lattice Boltzmann model
Appl. Math. Model.
Boundary conditions for thermal lattice Boltzmann equation method
J. Comput. Phys.
Lattice Boltzmann models for the convection-diffusion equation: D2Q5 vs D2Q9
Int. J. Heat Mass Transf.
An optimal two-relaxation-time lattice Boltzmann equation for solid-liquid phase change: the elimination of unphysical numerical diffusion
Int. J. Therm. Sci.
General regularized boundary condition for multi-speed lattice Boltzmann models
Comput. Fluids
A prospect for computing in porous materials research: very large fluid flow simulations
J. Comput. Sci.
Multigrid solution of the steady-state lattice Boltzmann equation
Comput. Fluids
Assessing moment-based boundary conditions for the lattice Boltzmann equation: a study of dipole-wall collisions
Comput. Fluids
Parallel performance of an IB-LBM suspension simulation framework
J. Comput. Sci.
An evaluation of lattice Boltzmann schemes for porous media simulations
J. Comput. Fluids
Implementation issues and benchmarking of lattice Boltzmann method for moving rigid particle simulations in a viscous flow
Comput. Math. Appl.
Steady state convergence acceleration of the generalized lattice Boltzmann equation with forcing term through preconditioning
J. Comput. Phys.
Error analysis and correction for Lattice Boltzmann simulated flow conductance in cappillaries of different shapes and alignements
J. Comput. Phys.
The permeability and quality of velocity field in a square array of solid and permeable cylindrical obstacles with the TRT-LBM and FEM Brinkman-schemes
CR Mechanique
Low- and high-order accurate boundary conditions: from Stokes to Darcy porous flow modeled with standard and improved Brinkman lattice Boltzmann schemes
J. Comput. Phys.
Second order accurate Lattice Boltzmann flow simulations in reconstructed porous media
Int. J. Comput. Fluid Dyn.
Lattice Boltzmann method for complex flows
Annu. Rev. Fluid Mech.
A Lattice Boltzmann Model for Diffusion of Binary Gas Mixtures. PhD Thesis
Momentum transfer of a Boltzmann-lattice fluid with boundaries
Phys. Fluids
The permeability of random medium: comparison of simulations with theory
Phys. Fluids A.
A multiple-relaxation-time lattice Boltzmann model for general anonlinear anistropic convection-diffusion equations
J. Sci. Comput.
Cited by (14)
Two pressure boundary conditions for multi-component multiphase flow simulations using the pseudo-potential lattice Boltzmann model
2022, Computers and FluidsCitation Excerpt :For such conditions, popular single-phase BC schemes are, for example, the bounce-back approach [27], the Non-Equilibrium Extrapolation Method (NEEM) [28], the Non-Equilibrium Bounce Back (NEBB) Method (sometimes referred to as the Zou–He Method [29]), and the anti-bounce-back (ABB) approach [30,31]. The ABB schemes in the last two works can prescribe pressure and velocity conditions at any distance to a straight boundary, and they have been recently extended to pressure boundary condition at a curved boundary [32]. Due to its simplicity and generalizability, the NEBB method has been applied to diverse pressure boundary condition simulations.
Lattice Boltzmann simulation of complex thermal flows via a simplified immersed boundary method
2022, Journal of Computational ScienceCitation Excerpt :In previous studies, the IB was frequently combined with some traditional NS-CD solvers [21], among which are the finite difference, finite element and finite volume schemes. Since the base of NS-CD equations, i.e., the continuous medium hypothesis would collapse in some scenarios, one can start from the kinetic theory [22] for interpreting the phenomenon of fluid motion as well as heat transfer. In this framework, the flow system is assumed to consist of some fictitious particles, which execute two procedures of collision and streaming continuously.
An exact-interface-fitted mesh generator and linearity-preserving finite volume scheme for anisotropic elliptic interface problems
2022, Journal of Computational PhysicsA conservative and stable explicit finite difference scheme for the diffusion equation
2021, Journal of Computational ScienceCitation Excerpt :It is also suitable for use in real-world problems because its implementation is simple and shows reasonable approximations. Ginzburg [22] proposes a simple and uniform modification of the local mass-conservation solvability condition to preserve mass when using the ADE scheme. Although various studies related to the ADE scheme have been performed because it is both simple and stable, in general, the ADE scheme does not conserve a mass, which is a disadvantage.
Spurious interface and boundary behaviour beyond physical solutions in lattice Boltzmann schemes
2021, Journal of Computational PhysicsCitation Excerpt :A poor accuracy because of an implicit, straight interface tangential-coupling has been also encountered with the discontinuous anisotropic weights, and improved for their leading-order mismatch [25]. Our recent work [35] confirms that the second-order interface-conjugate extension of this last approach vanishes the exponential accommodations of the arbitrary-placed rotated interface in the parabolic profiles. It would be interesting to examine and extend from this perspective the known interface, boundary and grid-refining approaches.