Search

We propose a procedure to implement Dirichlet velocity boundary conditions for complex shapes that use data from a single node only, in the context of the lattice Boltzmann method. Two ideas are at the base of this approach. The first is to generalize the geometrical description of boundary conditions combining bounce-back rule with interpolations. The second is to enhance them by limiting the interpolation extension to the proximity of the boundary. Despite its local nature, the resulting method exhibits second-order convergence for the velocity field and shows similar or better accuracy than the well-established Bouzidi's scheme for curved walls [M. Bouzidi, M. Firdaouss, and P. Lallemand, Phys. Fluids 13, 3452 (2001)]PHFLE61070-663110.1063/1.1399290. Among the infinite number of possibilities, we identify several meaningful variants of the method, discerned by their approximation of the second-order nonequilibrium terms and their interpolation coefficients. For each one, we provide two parametrized versions that produce viscosity independent accuracy at steady state. The method proves to be suitable to simulate moving rigid objects or surfaces moving following either the rigid body dynamics or a prescribed kinematic. Also, it applies uniformly and without modifications in the whole domain for any shape, including corners, narrow gaps, or any other singular geometry.

International audience; We introduce two new approaches, called A-LSOB and N-MR, for boundary and interface-conjugate conditions on flat or curved surface shapes in the advection-diffusion lattice Boltzmann method (LBM). The Local Second-Order, single-node A-LSOB enhances the existing Dirichlet and Neumann normal boundary treatments with respect to locality, accuracy, and Peclet parametrization. The normal-multi-reflection (N-MR) improves the directional flux schemes via a local release of their nonphysical tangential constraints. The A-LSOB and N-MR restore all first- and second-order derivatives from the nodal non-equilibrium solution, and they are conditioned to be exact on a piece-wise parabolic profile in a uniform arbitrary-oriented tangential velocity field. Additionally, the most compact and accurate single-node parabolic schemes for diffusion and flow in grid-inclined pipes are introduced. In simulations, the global mass-conservation solvability condition of the steady-state, two-relaxation-time (S-TRT) formulation is adjusted with either (i) a uniform mass-source or (ii) a corrective surface-flux. We conclude that (i) the surface-flux counterbalance is more accurate than the bulk one, (ii) the A-LSOB Dirichlet schemes are more accurate than the directional ones in the high Peclet regime, (iii) the directional Neumann advective-diffusive flux scheme shows the best conservation properties and then the best performance both in the tangential no-slip and interface-perpendicular flow, and (iv) the directional non-equilibrium diffusive flux extrapolation is the least conserving and accurate. The error Peclet dependency, Neumann invariance over an additive constant, and truncation isotropy guide this analysis. Our methodology extends from the d2q9 isotropic S-TRT to 3D anisotropic matrix collisions, Robin boundary condition, and the transient LBM.

This work addresses the Dirichlet boundary condition for momentum in the lattice Boltzmann method (LBM), with focus on the steady-state Stokes flow modelling inside non-trivial shaped ducts. For this task, we revisit a local and highly accurate boundary scheme, called the local second-order boundary (LSOB) method. This work reformulates the LSOB within the two-relaxation-time (TRT) framework, which achieves a more standardized and easy to use algorithm due to the pivotal parametrization TRT properties. The LSOB explicitly reconstructs the unknown boundary populations in the form of a Chapman–Enskog expansion, where not only first- but also second-order momentum derivatives are locally extracted with the TRT symmetry argument, through a simple local linear algebra procedure, with no need to compute their non-local finite-difference approximations. Here, two LSOB strategies are considered to realize the wall boundary condition, the original one called Lwall and a novel one Lnode, which operate with the wall and node variables, roughly speaking. These two approaches are worked out for both plane and curved walls, including the corners. Their performance is assessed against well-established LBM boundary schemes such as the bounce-back, the local second-order accurate CLI scheme and two different parabolic multi-reflection (MR) schemes. They are all evaluated for 3D duct flows with rectangular, triangular, circular and annular cross-sections, mimicking the geometrical challenges of real porous structures. Numerical tests confirm that LSOB competes with the parabolic MR accuracy in this problem class, requiring only a single node to operate. This article is part of the theme issue ‘Fluid dynamics, soft matter and complex systems: recent results and new methods’.

The extended method of moments (EMM) is elaborated in recursive algorithmic form for the prediction of the effective diffusivity, the Taylor dispersion dyadic and the associated longitudinal high-order coefficients in mean-concentration profiles and residence-time distributions. The method applies in any streamwise-periodic stationary d-dimensional velocity field resolved in the piecewise continuous heterogeneous porosity field. It is demonstrated that EMM reduces to the method of moments and the volume-averaging formulation in microscopic velocity field and homogeneous soil, respectively. The EMM simultaneously constructs two systems of moments, the spatial and the temporal, without resorting to solving of the high-order upscaled PDE. At the same time, the EMM is supported with the reconstruction of distribution from its moments, allowing to visualize the deviation from the classical ADE solution. The EMM can be handled by any linear advection-diffusion solver with explicit mass-source and diffusive-flux jump condition on the solid boundary and permeable interface. The prediction of the first four moments is decisive in the optimization of the dispersion, asymmetry, peakedness and heavy-tails of the solute distributions, through an adequate design of the composite materials, wetlands, chemical devices or oil recovery. The symbolic solutions for dispersion, skewness and kurtosis are constructed in basic configurations: diffusion process and Darcy flow through two porous blocks in “series”, straight and radial Poiseuille flow, porous flow governed by the Stokes–Brinkman–Darcy channel equation and a fracture surrounded by penetrable diffusive matrix or embedded in porous flow. We examine the moments dependency upon porosity contrast, aspect ratio, Peclet and Darcy numbers, but also for their response on the effective Brinkman viscosity applied in flow modeling. Two numerical Lattice Boltzmann algorithms, a direct solver of the microscopic ADE in heterogeneous structure and a novel scheme for EMM numerical formulation, are called for validation of the constructed analytical predictions.

International audience; A simple local two-relaxation-time Lattice Boltzmann numerical formulation (TRT-EMM) of the extended method of moments (EMM) is proposed for analysis of the spatial and temporal Taylor dispersion in d-dimensional streamwise-periodic stationary mesoscopic velocity field resolved in a piecewise-continuous porous media. The method provides an effective diffusivity, dispersion, skewness and kurtosis of the mean concentration profile and residence time distribution. The TRT-EMM solves a chain of steady-state heterogeneous advection-diffusion equations with the pre-computed space-variable mass-source and automatically undergoes diffusion-flux jump on the abrupt-porosity streamwise-normal interface. The temporal and spatial systems of moments are computed within the same run; the symmetric dispersion tensor can be restored from independent computations performed for each periodic mean-velocity axis; the numerical algorithm recursively extends for any order moment. We derive an exact form of the bulk equation and implicit closure relations, construct symbolic TRT-EMM solutions and determine specific relation between the equilibrium and the collision degrees of freedom viewing an exact parameterization by the physical non-dimensional numbers in two alternate situations: "parallel" fracture/matrix flow and "perpendicular" Darcy flow through porous blocks in "series". Two-dimensional simulations in linear Brinkman flow around solid and through porous obstacles validate the method in comparison with the two heterogeneous direct LBM-ADE schemes with different anti-numerical-diffusion treatment which are proposed and examined in parallel. On the coarse grid, accuracy of the three moments is essentially determined by the free-tunable collision rate in all schemes, and especially TRT-EMM. However, operated within a single periodic cell, the TRT-EMM is many orders of magnitude faster than the direct solvers, numerical-diffusion free, more robust and much more capable for accuracy improving, high Péclet range and free-parameter influence reduction with the mesh refinement. The method is an efficient predicting tool for the Taylor dispersion, asymmetry and peakedness; moreover, it allows for an optimal analysis between the mutual effect of the flow regime, Péclet number, porosity, permeability and obstruction geometry.

Impact of the unphysical tangential advective-diffusion constraint of the bounce-back (BB) reflection on the impermeable solid surface is examined for the first four moments of concentration. Despite the number of recent improvements for the Neumann condition in the lattice Boltzmann method-advection-diffusion equation, the BB rule remains the only known local mass-conserving no-flux condition suitable for staircase porous geometry. We examine the closure relation of the BB rule in straight channel and cylindrical capillary analytically, and show that it excites the Knudsen-type boundary layers in the nonequilibrium solution for full-weight equilibrium stencil. Although the d2Q5 and d3Q7 coordinate schemes are sufficient for the modeling of isotropic diffusion, the full-weight stencils are appealing for their advanced stability, isotropy, anisotropy and anti-numerical-diffusion ability. The boundary layers are not covered by the Chapman-Enskog expansion around the expected equilibrium, but they accommodate the Chapman-Enskog expansion in the bulk with the closure relation of the bounce-back rule. We show that the induced boundary layers introduce first-order errors in two primary transport properties, namely, mean velocity (first moment) and molecular diffusion coefficient (second moment). As a side effect, the Taylor-dispersion coefficient (second moment), skewness (third moment), and kurtosis (fourth moment) deviate from their physical values and predictions of the fourth-order Chapman-Enskog analysis, even though the kurtosis error in pure diffusion does not depend on grid resolution. In two- and three-dimensional grid-aligned channels and open-tubular conduits, the errors of velocity and diffusion are proportional to the diagonal weight values of the corresponding equilibrium terms. The d2Q5 and d3Q7 schemes do not suffer from this deficiency in grid-aligned geometries but they cannot avoid it if the boundaries are not parallel to the coordinate lines. In order to vanish or attenuate the disparity of the modeled transport coefficients with the equilibrium weights without any modification of the BB rule, we propose to use the two-relaxation-times collision operator with free-tunable product of two eigenfunctions Λ. Two different values Λ_{v} and Λ_{b} are assigned for bulk and boundary nodes, respectively. The rationale behind this is that Λ_{v} is adjustable for stability, accuracy, or other purposes, while the corresponding Λ_{b}(Λ_{v}) controls the primary accommodation effects. Two distinguished but similar functional relations Λ_{b}(Λ_{v}) are constructed analytically: they preserve advection velocity in parabolic profile, exactly in the two-dimensional channel and very accurately in a three-dimensional cylindrical capillary. For any velocity-weight stencil, the (local) double-Λ BB scheme produces quasi-identical solutions with the (nonlocal) specular-forward reflection for first four moments in a channel. In a capillary, this strategy allows for the accurate modeling of the Taylor-dispersion and non-Gaussian effects. As illustrative example, it is shown that in the flow around a circular obstacle, the double-Λ scheme may also vanish the dependency of mean velocity on the velocity weight; the required value for Λ_{b}(Λ_{v}) can be identified in a few bisection iterations in given geometry. A positive solution for Λ_{b}(Λ_{v}) may not exist in pure diffusion, but a sufficiently small value of Λ_{b} significantly reduces the disparity in diffusion coefficient with the mass weight in ducts and in the presence of rectangular obstacles. Although Λ_{b} also controls the effective position of straight or curved boundaries, the double-Λ scheme deals with the lower-order effects. Its idea and construction may help understanding and amelioration of the anomalous, zero- and first-order behavior of the macroscopic solution in the presence of the bulk and boundary or interface discontinuities, commonly found in multiphase flow and heterogeneous transport.

International audience; The effect of the heterogeneity in the soil structure or the nonuniformity of the velocity field on the modeled resident time distribution (RTD) and breakthrough curves is quantified by their moments. While the first moment provides the effective velocity, the second moment is related to the longitudinal dispersion coefficient ( k T )in the developed Taylor regime; the third and fourth moments are characterized by their normalized values skewness (Sk) and kurtosis (Ku), respectively. The purpose of this investigation is to examine the role of the truncation corrections of the numerical scheme in k T, Sk, and Ku because of their interference with the second moment, in the form of the numerical dispersion, and in the higher-order moments, by their definition. Our symbolic procedure is based on the recently proposed extended method of moments (EMM). Originally, the EMM restores any-order physical mo- ments of the RTD or averaged distributions assuming that the solute concentration obeys the advection-diffusion equation in multidimensional steady-state velocity field, in streamwise-periodic heterogeneous structure. In our work, the EMM is generalized to the fourth-order-accurate apparent mass-conservation equation in two- and three-dimensional duct flows. The method looks for the solution of the transport equation as the product of a long harmonic wave and a spatially periodic oscillating component; the moments of the given numerical scheme are derived from a chain of the steady-state fourth-order equations at a single cell. This mathematical technique is exemplified for the truncation terms of the two-relaxation-time lattice Boltzmann scheme, using plug and parabolic flow in straight channel and cylindrical capillary with the d2Q9 and d3Q15 discrete velocity sets as simple but illustrative examples. The derived symbolic dependencies can be readily extended for advection by another, Newtonian or non-Newtonian, flow profile in any-shape open-tabular conduits. It is established that the truncation errors in the three transport coefficients k T , Sk, and Ku decay with the second-order accuracy. While the physical values of the three transport coefficients are set by Peclet number, their truncation corrections additionally depend on the two adjustable relaxation rates and the two adjustable equilibrium weight families which independently determine the convective and diffusion discretization stencils. We identify flow- and dimension-independent optimal strategies for adjustable parameters and confront them to stability requirements. Through specific choices of two relaxation rates and weights, we expect our results be directly applicable to forward-time central differences and leap-frog central-convective Du Fort-Frankel-diffusion schemes. In straight channel, a quasi-exact validation of the truncation predictions through the numerical moments becomes possible thanks to the specular-forward no-flux boundary rule. In the staircase description of a cylindrical capillary, we account for the spurious boundary- layer diffusion and dispersion because of the tangential constraint of the bounce-back no-flux boundary rule.

International audience; An analytical study is devised for the problem of bimodal porous flow across a periodic array of permeable cylindrical inclusions. Such a configuration is particularly relevant for porous media systems of dual granulometry, an idealization often taken, e.g. in the modelling of membranes and fibrous applications. The double-porosity system is governed by the Stokes–Brinkman–Darcy equations, the most general description in this class of flow problems characterized by the permeabilities of the surrounding matrix and inclusions, their porosities and the relative volume fraction. We solve this problem with the Kuwabara cell model and lubrication approach, providing analytical solutions for the system effective permeability in closed analytical form.The ensemble of results demonstrates the self-consistency of the bimodal solutions in eight possible limit configurations and supports the validity of the Beavers–Joseph interface stress jump condition for transmission from the open Stokes flow to low-permeable Darcy region. At the same time, these solutions bring further insight on the relative significance of the governing parameters on the effective permeability, with a focus on the role of the effective viscosity (porosity) distribution. Furthermore, although the cell model is restricted to relatively small volume fractions in openflow, its validity extends in less-permeable background flow inside Brinkman/Brinkman description. In turn, the lubrication approximation remains more adequate in the opposite limit of the dense impermeable inclusions. These conclusions are drawn from comparisons with the numerical solutions obtained with the developed lattice Boltzmann model and the standard finite element method. The two methods principally differ in the treatment of the interface conditions: implicit and explicit, respectively. The purpose of this task is therefore twofold. While the numerical schemes help quantifying the validity limits of the theoretical approach, the analytical solutions offer a non-trivial benchmark for numerical schemes in highly heterogeneous soil.

This paper establishes relations between the stability and the high-order truncated corrections for modeling of the mass conservation equation with the two-relaxation-times (TRT) collision operator. First we propose a simple method to derive the truncation errors from the exact, central-difference type, recurrence equations of the TRT scheme. They also supply its equivalent three-time-level discretization form. Two different relationships of the two relaxation rates nullify the third (advection) and fourth (pure diffusion) truncation errors, for any linear equilibrium and any velocity set. However, the two relaxation times alone cannot remove the leading-order advection-d if fusion error, because of the intrinsic fourth-order numerical diffusion. The truncation analysis is carefully verified for the evolution of concentration waves with the anisotropic diffusion tensors. The anisotropic equilibrium functions are presented in a simple but general form, suitable for the minimal velocity sets and the d2Q9, d3Q13, d3Q15 and d3Q19 velocity sets. All anisotropic schemes are complemented by their exact necessary von Neumann stability conditions and equivalent finite-difference stencils. The sufficient stability conditions are proposed for the most stable (OTRT) family, which enables modeling at any Peclet numbers with the same velocity amplitude. The heuristic stability analysis of the fourth-order truncated corrections extends the optimal stability to larger relationships of the two relaxation rates, in agreement with the exact (one-dimensional) and numerical (multi-dimensional) stability analysis. A special attention is put on the choice of the equilibrium weights. By combining accuracy and stability predictions, several strategies for selecting the relaxation and free-tunable equilibrium parameters are suggested and applied to the evolution of the Gaussian hill.

This work focuses on the numerical solution of the Stokes-Brinkman equation for a voxel-type porous-media grid, resolved by one to eight spacings per permeability contrast of 1 to 10 orders in magnitude. It is first analytically demonstrated that the lattice Boltzmann method (LBM) and the linear-finite-element method (FEM) both suffer from the viscosity correction induced by the linear variation of the resistance with the velocity. This numerical artefact may lead to an apparent negative viscosity in low-permeable blocks, inducing spurious velocity oscillations. The two-relaxation-times (TRT) LBM may control this effect thanks to free-tunable two-rates combination Λ. Moreover, the Brinkman-force-based BF-TRT schemes may maintain the nondimensional Darcy group and produce viscosity-independent permeability provided that the spatial distribution of Λ is fixed independently of the kinematic viscosity. Such a property is lost not only in the BF-BGK scheme but also by "partial bounce-back" TRT gray models, as shown in this work. Further, we propose a consistent and improved IBF-TRT model which vanishes viscosity correction via simple specific adjusting of the viscous-mode relaxation rate to local permeability value. This prevents the model from velocity fluctuations and, in parallel, improves for effective permeability measurements, from porous channel to multidimensions. The framework of our exact analysis employs a symbolic approach developed for both LBM and FEM in single and stratified, unconfined, and bounded channels. It shows that even with similar bulk discretization, BF, IBF, and FEM may manifest quite different velocity profiles on the coarse grids due to their intrinsic contrasts in the setting of interface continuity and no-slip conditions. While FEM enforces them on the grid vertexes, the LBM prescribes them implicitly. We derive effective LBM continuity conditions and show that the heterogeneous viscosity correction impacts them, a property also shared by FEM for shear stress. But, in contrast with FEM, effective velocity conditions in LBM give rise to slip velocity jumps which depend on (i) neighbor permeability values, (ii) resolution, and (iii) control parameter Λ, ranging its reliable values from Poiseuille bounce-back solution in open flow to zero in Darcy's limit. We suggest an "upscaling" algorithm for Λ, from multilayers to multidimensions in random extremely dispersive samples. Finally, on the positive side for LBM besides its overall versatility, the implicit boundary layers allow for smooth accommodation of the flat discontinuous Darcy profiles, quite deficient in FEM.

International audience; This work demonstrates that in advection–diffusion Lattice Boltzmann schemes, the local mass-conserving boundary rules, such as bounce-back and local specular reflection, may modify the transport coefficients predicted by the Chapman–Enskog expansion when they enforce to zero not only the normal, but also the tangential boundary flux. In order to accommodate it to the bulk solution, the system develops a Knudsen-layer correction to the non-equilibrium part of the population solution. Two principal secondary effects- (i) decrease in the diffusion coefficient, and (ii) retardation of the average advection velocity, obtained in a closed analytical form, are proportional, respectively, to freely assigned diagonal weights for equilibrium mass and velocity terms. In addition, due to their transverse velocity gradients, the boundary layers affect the longitudinal diffusion coefficient similarly to Taylor dispersion, as they grow as the square of the Péclet number. These numerical artifacts can be eliminated or reduced by a proper space distribution of the free-tunable collision eigenvalue in two-relaxation-time schemes.

We analyze the intrinsic impact of free-tunable combinations of the relaxation rates controlling viscosity-independent accuracy of the multiple-relaxation-times (MRT) lattice-Boltzmann models. Preserving all MRT degrees of freedom, we formulate the parametrization conditions which enable the MRT schemes to provide viscosity-independent truncation errors for steady state solutions, and support them with the second- and third-order accurate ("linear" and "parabolic", respectively) boundary schemes. The parabolic schemes demonstrate the advanced accuracy with weak dependency on the relaxation rates, as confirmed by the simulations with the D3Q15 model in three regular arrays (SC, BCC, FCC) of touching spheres. Yet, the low-order, bounce-back boundary rule remains appealing for pore-scale simulations where the precise distance to the boundaries is undetermined. However, the effective accuracy of the bounce-back crucially depends on the free-tunable combinations of the relaxation rates. We find that the combinations of the kinematic viscosity rate with the available "ghost" antisymmetric collision mode rates mainly impact the accuracy of the bounce-back scheme. As the first step, we reduce them to the one combination (presented by so-called "magic" parameter ? in the frame of the two-relaxation-times (TRT) model), and study its impact on the accuracy of the drag force/permeability computations with the D3Q19 velocity set in two different, dense, random packings of 8000 spheres each. We also run the simulations in the regular (BCC and FCC) packings of the same porosity for the broad range of the discretization resolutions, ranging from 5 to 750 lattice nodes per sphere diameter. A special attention is given to the discretization procedure resulting in significantly reduced scatter of the data obtained at low resolutions. The results reveal the identical ?-dependency versus the discretization resolution in all four packings, regular and random. While very small ? values overestimate the drag measurements several-fold on the coarse grids, ? 1 may overestimate the permeability at the same extent. In low resolution region we provide practical guidelines, extending previously known solutions for the straight/diagonal Poiseuille flow. Analysis of the high-resolution region reveals the collapse of the solutions obtained with all the considered ? values with the average rate of -1.3, followed by their common, smooth, first-order convergence with the rate of -1.0 as the best, towards the reference solutions provided by the "parabolic" schemes. High-quality power-law fits estimate that the bounce-back would reach their accuracy (obtained at about 200 nodes per sphere) for two-order magnitude higher grid resolution. The parametrization conditions enable the MRT LBE to provide viscosity-independent truncation spatial errors.These schemes are supported with the high-order accurate boundary conditions confirmed in touching arrays of spheres.The effect of free relaxation rates on drag-force computations is examined for bounce-back and high-order conditions.Random and regular packings demonstrate identical numerical errors for various relaxation rates and resolutions.For low resolution we provide practical guidelines while high resolutions revealed first-order convergence to reference drag.

International audience; Using as a benchmark the porous flow in a square array of solid or permeable cylindrical obstacles, we evaluate the numerical performance of the two-relaxation-time lattice Boltzmann method (TRT–LBM) and the linear finite element method (FEM). We analyze the bulk, boundary and interface properties of the Brinkman-based schemes in staircase discretization on the voxel-type grids typical of porous media simulations. The effect of flow regime, grid resolution, and TRT collision degree of freedom is assessed. In coarse meshes, the TRT may outperform the FEM by properly selecting . Further, FEM is more oscillatory, a defect virtually suppressed in TRT with an improved strategy IBF and implicit accommodation of interface/boundary layers.

International audience; We extend lattice Boltzmann (LB) methods to advection and anisotropic-dispersion equations (AADE). LB methods are advocated for the exactness of their conservation laws, the handling of different length and time scales for flow/transport problems, their locality and extreme simplicity. Their extension to anisotropic collision operators (L-model) and anisotropic equilibrium distributions (E-model) allows to apply them to generic diffusion forms. The AADE in a conventional form can be solved by the L-model. Based on a link-type collision operator, the L-model specifies the coefficients of the symmetric diffusion tensor as linear combination of its eigenvalue functions. For any type of collision operator, the E-model constructs the coefficients of the transformed diffusion tensors from linear combinations of the relevant equilibrium projections. The model is able to eliminate the second order tensor of its numerical diffusion. Both models rely on mass conserving equilibrium functions and may enhance the accuracy and stability of the isotropic convection-diffusion LB models.The link basis is introduced as an alternative to a polynomial collision basis. They coincide for one particular eigenvalue configuration,the two-relaxation-time (TRT) collision operator, suitable for both mass and momentum conservation laws. TRT operator is equivalent to the BGK collision in simplicity but the additional collision freedom relates it to multiple-relaxation-times (MRT) models. Optimal convection and optimal diffusion eigenvalue solutions for the TRT E-model allow to remove next order corrections to AADE. Numerical results confirm the Chapman-Enskog and dispersion analysis.

International audience; We address a multi-reflection approach to model Dirichlet and Neumann time-dependent boundary conditions in lattice Boltzmann methods for arbitrarily shaped surfaces. The multi-reflection condition for an incoming population represents a linear combination of the known population solutions. The closure relations are first established for symmetric and anti-symmetric parts of the equilibrium functions, independently of the nature of the problem. The symmetric part is tuned to build second- and third-order accurate Dirichlet boundary conditions for the scalar function specified by the equilibrium distribution. The focus is on two approaches to advection and anisotropic-dispersion equations (AADE): the equilibrium technique when the coefficients of the expanded equilibrium functions match the coefficients of the transformed dispersion tensor, and the eigenvalue technique when the coefficients of the dispersion tensor are built as linear combinations of the eigenvalue functions associated with the link-type collision operator. As a particular local boundary technique, the anti-bounce-back condition is analyzed. The anti-symmetric part of the generic closure relation allows to specify normal flux conditions without inversion of the diffusion tensor. Normal and tangential constraints are derived for bounce-back and specular reflections. The bounce-back closure relation is released from the nonphysical tangential flux restriction at leading orders. Solutions for the Poisson equation and for convection-diffusion equations are presented for isotropic/anisotropic configurations with specified Dirichlet and Neumann boundary conditions.

International audience; This article describes a generalization of the method of moments, called extended method of moments (EMM), for dispersion in periodic structures composed of impermeable or permeable porous inclusions. Prescribing pre-computed steady state velocity field in a single periodic cell, the EMM sequentially solves specific linear stationary advection-diffusion equations and restores any-order moments of the resident time distribution or the averaged concentration distribution. Like the pioneering Brenner's method, the EMM recovers mean seepage velocity and Taylor dispersion coefficient as the first two terms of the perturbative expansion. We consider two types of dispersion: spatial dispersion, i.e., spread of initially narrow pulse of concentration, and temporal dispersion, where different portions of the solute have different residence times inside the system. While the first (mean velocity) and the second (Taylor dispersion coefficient) moments coincide for both problems, the higher moments are different. Our perturbative approach allows to link them through simple analytical expressions. Although the relative importance of the higher moments decays downstream, they manifest the non-Gaussian behaviour of the breakthrough curves, especially if the solute can diffuse into less porous phase. The EMM quantifies two principal effects of bi-modality, as the appearance of sharp peaks and elongated tails of the distributions. In addition, the moments can be used for the numerical reconstruction of the corresponding distribution, avoiding time-consuming computations of solute transition through heterogeneous media. As illustration, solutions for Taylor dispersion, skewness, and kurtosis in Poiseuille flow and open/impermeable stratified systems, both in rectangular and cylindrical channels, power-law duct flows, shallow channels, and Darcy flow in parallel porous layers are obtained in closed analytical form for the entire range of Péclet numbers. The high-order moments and reconstructed profiles are compared to their predictions from the advection-diffusion equation for averaged concentration, based on the same averaged seepage velocity and Taylor dispersion coefficient. In parallel, we construct Lattice-Boltzmann equation (LBE) two-relaxation-times scheme to simulate transport of a passive scalar directly in heterogeneous media specified by discontinuous porosity distribution. We focus our numerical analysis and assessment on (i) truncation corrections, because of their impact on the moments, (ii) stability, since we show that stable Darcy velocity amplitude reduces with the porosity, and (iii) interface accuracy which is found to play the crucial role. The task is twofold: the LBE supports the EMM predictions, while the EMM provides non-trivial benchmarks for the numerical schemes.

International audience; This paper develops a symmetrized framework for the analysis of the anisotropic advection–diffusion Lattice Boltzmann schemes. Two main approaches build the anisotropic diffusion coefficients either from the anisotropic anti-symmetric collision matrix or from the anisotropic symmetric equilibrium distribution. We combine and extend existing approaches for all commonly used velocity sets, prescribe most general equilibrium and build the diffusion and numerical-diffusion forms, then derive and compare solvability conditions, examine available anisotropy and stable velocity magnitudes in the presence of advection. Besides the deterioration of accuracy, the numerical diffusion dictates the stable velocity range. Three techniques are proposed for its elimination: (i) velocity-dependent relaxation entries; (ii) their combination with the coordinate-link equilibrium correction; and (iii) equilibrium correction for all links. Two first techniques are also available for the minimal (coordinate) velocity sets. Even then, the two-relaxation-times model with the isotropic rates often gains in effective stability and accuracy. The key point is that the symmetric collision mode does not modify the modeled diffusion tensor but it controls the effective accuracy and stability, via eigenvalue combinations of the opposite parity eigenmodes. We propose to reduce the eigenvalue spectrum by properly combining different anisotropic collision elements. The stability role of the symmetric, multiple relaxation-times component, is further investigated with the exact von Neumann stability analysis developed in diffusion-dominant limit.

International audience; Despite the growing popularity of Lattice Boltzmann schemes for describing multi-dimensional flow and transport governed by non-linear (anisotropic) advectiondiffusion equations, there are very few analytical results on their stability, even for the isotropic linear equation. In this paper, the optimal two-relaxation-time (OTRT) model is defined, along with necessary and sufficient (easy to use) von Neumann stability conditions for a very general anisotropic advection-diffusion equilibrium, in one to three dimensions, with or without numerical diffusion. Quite remarkably, the OTRT stability bounds are the same for any Peclet number and they are defined by the adjustable equilibrium parameters. Such optimal stability is reached owing to the free (kinetic) relaxation parameter. Furthermore, the sufficient stability bounds tolerate negative equilibrium functions (the distribution divided by the local mass), often labeled as unphysical. We prove that the non-negativity condition is (i) a sufficient stability condition of the TRT model with any eigenvalues for the pure diffusion equation, (ii) a sufficient stability condition of its OTRT and BGK/SRT subclasses, for any linear anisotropic advection-diffusion equation, and (iii) unnecessarily more restrictive for any Peclet number than the optimal sufficient conditions. Adequate choices of the two relaxation rates and the free-tunable equilibrium parameters make the OTRT subclass more efficient than the BGK one, at least in the advection-dominant regime, and allow larger time steps than known criteria of the forward time central finite-difference schemes (FTCS/MFTCS) for both, advection and diffusion dominant regimes.

International audience; We show that a consistent modeling of porous flows needs at least one free collision relaxation rate to avoid a nonlinear dependency of the numerical errors on the viscosity. This condition is necessary to get the viscosity-independent permeability from the Stokes flow and to parametrize properly with nondimensional physical numbers the lattice Boltzmann Brinkman schemes. The two-relaxation-time TRT collision operator controls all coefficients of the higher-order corrections in steady solutions with a specific combination of its two collision rates, a possibility lacking for the Bhatnagar-Gross-Krook BGK-based single-relaxation-time schemes. The analysis is based on exact recurrence equations of the evolution equation and illustrated for the exact solutions of the Brinkman scheme in simply oriented parallel and diagonal channels. The apparent viscosity coefficient of the TRT Stokes-Brinkman scheme in arbitrary flow is only approximated. The compatibility of one-dimensional arbitrarily rotated flows with the nonlinear Navier-Stokes equilibrium is examined. An explicit dependency for all coefficients on the relaxation rates is presented for the infinite steady state Chapman-Enskog expansion.

[Departement_IRSTEA]RE [TR1_IRSTEA]RIE / PHYLEAU; International audience; We develop a two-relaxation-time (TRT) Lattice Boltzmann model for hydrodynamic equations with variable source terms based on equivalent equilibrium functions. A special parametrization of the free relaxation parameter is derived. It controls, in addition to the non-dimensional hydrodynamic numbers, any TRT macroscopic steady solution and governs the spatial discretization of transient flows. In this framework, the multi-reflection approach [16, 18] is generalized and extended for Dirichlet velocity, pressure and mixed pressure/tangential velocity) boundary conditions. We propose second and third-order accurate boundary schemes and adapt them for corners. The boundary schemes are analyzed for exactness of the parametrization, uniqueness of their steady solutions, support of staggered invariants and for the effective accuracy in case of time dependent boundary conditions and transient flow. When the boundary scheme obeys the parametrization properly, the derived permeability values become independent of the selected viscosity for any porous structure and can be computed efficiently. The linear interpolations [5, 46] are improved with respect to this property.

International audience; Research conducted for the last 35 years has shown that subsurface drainage has a significant impact on hydrology and contaminant transport. This can be observed at the field-scale and also at the watershed scale. Impacts are always associated with modifying otherwise natural flow paths. Most computer model representations of drainage have been drawn at the field-scale. These models require relatively precise data that are usually unavailable when simulating hydrology and water quality in large watersheds. We believe that in this case drainage representation should be simplified and yet closely match observations. As a first step towards incorporating drainage systems into large-scale hydrological models, we propose an equivalent representation of drains buried in a soil profile by using a homogeneous anisotropic porous medium without drains. This representation is based on a self-consistent' approach and on geometrical considerations. Simplification is such that calculating the equivalent hydraulic conductivity requires only information on the main length and spacing of the tile drains and not on their precise location. This approach also provides a much simpler discretisation of the domain because of the absence of internal boundary conditions on the drainage pipes. Compared to other methods that have simplified drainage representation in existing watershed models, it requires no parameter fitting. Two alternatives to the method are presented: in the first one, the soil profile equipped with the actual drain pipes is represented by an equivalent, horizontally layered system with no pipes; in the second, the layered system has been replaced with an equivalent homogeneous profile. The efficiency of these approaches was tested against a classical representation of tile drains using the SWMS 3D code, which solves the Richards equation for a typical drained plot configuration. The equivalent-medium approach appears to give satisfying results for global water outflow and mean water table elevation.

[Departement_IRSTEA]RE [TR1_IRSTEA]RIE / PHYLEAU; International audience; We present analytical and Lattice Boltzmann (LB) solutions for steady-state saturated flows in 2D and 3D anisotropic heterogeneous aquifers. The analytical solution is easy to use and extends the known ones for ground-water whirls to more general combinations of the anisotropic properties of two-layered systems. The Bakker and Hemker's multi-layered' semi-analytical solution and the LB results are compared to the analytical solution for a broad range of anisotropic heterogeneous diffusion tensors. The main components of the LB scheme, the eigenvalues of the linear collision operator and/or the equilibrium functions, become discontinuous when the anisotropy changes between the layers. It is shown that the evolution equation of the LB method needs to be modified at the interfaces in order to satisfy the continuity conditions for the diffusion function and/or its tangential derivatives. The existing LB schemes for anisotropic advectiondispersion equations are formulated in a more general framework in which the leading-order interface corrections are constructed and analyzed for linear and highly nonlinear exact solutions. We also present some stability aspects of these schemes, introduce specified normal gradient boundary conditions and discuss the computation of total and local fluxes. The interface analysis developed here applies to generic LB schemes with discontinuous collision operators.

International audience; Irrespective of the nature of the modeled conservation laws, we establish first the microscopic interface continuity conditions for Lattice Boltzmann (LB) multiple-relaxation time, link-wise collision operators with discontinuous components (equilibrium functions and/or relaxation parameters). Effective macroscopic continuity conditions are derived for a planar implicit interface between two immiscible fluids, described by the simple two phase hydrodynamic model, and for an implicit interface boundary between two heterogeneous and anisotropic, variably saturated soils, described by Richard's equation. Comparing the effective macroscopic conditions to the physical ones,we show that the range of the accessible parameters is restricted, e.g. a variation of fluid densities or a heterogeneity of the anisotropic soil properties. When the interface is explicitly tracked, the interface collision components are derived from the leading order continuity conditions. Among particular interface solutions, a harmonic mean value is found to be an exact LB solution, both for the interface kinematic viscosity and for the interface vertical hydraulic conductivity function. We construct simple problems with the explicit and implicit interfaces, matched exactly by the LB hydrodynamic and/or advectiondiffusion schemes with the aid of special solutions for free collision parameters.

International audience; This paper addresses the numerical solution of highly nonlinear parabolic equations with Lattice Boltzmann techniques. They are first developed for generic advection and anisotropic dispersion equations (AADE). Collision configurations handle the anisotropic diffusion forms by using either anisotropic eigenvalue sets or anisotropic equilibrium functions. The coordinate transformation from the orthorhombic (rectangular) discretization grid to the cuboid computational grid is equivalent for the AADE to the anisotropic rescaling of the convection/diffusion terms. The collision components (eigenvalues and/or equilibrium functions) become discontinuous on the boundaries of the computational sub-domains which have different space scaling factors. We focus on the analysis of the boundary continuity conditions by using anisotropic LB techniques. The developed schemes are applied to Richards_ equation for variably saturated flow. The anisotropy of the Richard_s equation originates from distinct soil conductivity values in both the vertical and horizontal directions. The method should on the interface between the heterogeneous layers maintain the continuity of the normal component of the Darcy_s velocity (total flux). Also the method should accommodate steep jumps of the moisture content variable (conserved quantity) resulting from the continuity of the pressure variable, a given non-linear function of the moisture content. The coupling between heterogeneity and the anisotropy is examined by using the distinct space steps in neighboring layers and tested against uniform grid solutions. Different formulations of the Richard_s equation illustrate the construction of distinct diffusion forms and their integral transforms via specification of the equilibrium components. Integral transforms are used to overcome the difficulties coming from the rapid change of the main variables on sharp fronts. The numerical assessment of the stability criteria and the interface boundary conditions extend the analysis of the Lattice Boltzmann schemes to nonlinear problems with discontinuous coefficients.

[Departement_IRSTEA]RE [TR1_IRSTEA]RIE / PHYLEAU; International audience; This paper presents a model to simulate overland flow genesis induced by shallow water table movements in hillslopes. Variably saturated subsurface flows are governed by the Richards equation discretized by continuous finite elements on unstructured meshes. An obstacle-type formulation is used to determine where saturation conditions, and thus seepage face conditions, are met at the ground surface. The impact of hillslope geometry, boundary conditions, and soil hydraulic parameters on model predictions is investigated on two-dimensional test cases at the metric and hectometric scales. The obstacle-type formulation is also compared with a more detailed model coupling subsurface and overland flow, the latter being described by the shallow water equations in the diffusive wave regime.

This article provides a concise exposition of the multiple-relaxation-time lattice Boltzmann equation, with examples of 15-velocity and 19-velocity models in three dimensions. Simulation of a diagonally lid-driven cavity flow in three dimensions at Re = 500 and 2000 is performed. The results clearly demonstrate the superior numerical stability of the multiple-relaxation-time lattice Boltzmann equation over the popular lattice Bhatnagar-Gross-Krook equation.

The filling process of viscoplastic metal alloys and plastics in expanding cavities is modelled using the lattice Boltzmann method in two and three dimensions. These models combine the regularized Bingham model for viscoplastic fluids with a free interface algorithm. The latter is based on a modified immiscible lattice Boltzmann model in which. one species is the fluid and the other one is considered to be a vacuum. The boundary conditions at the curved liquid-vacuum interface are met without any geometrical front reconstruction from a first-order Chapman- Enskog expansion. The numerical results obtained with these models are found in good agreement with available theoretical and numerical analysis.

During heavy rainfall episodes, subsurface flow can saturate the soil in various regions near the surface and, therefore, contribute to the production of overland flow. This paper investigates numerically the dynamics of water tables in partially saturated porous media and their role in the genesis of surface runoff. The water movement in variably saturated soils is modeled with Richards' equation. The water table position being an unknown of the problem, its intersection with the ground surface yields an unsteady obstacle-type problem. The governing equations are discretized by finite elements in space and an implicit Euler scheme in time. At each time step, the approximate solution is obtained using Newton's method embedded into a fixed-point iteration to determine those points lying on the soil surface where artesian conditions are met. Numerical results are presented for a simple hillslope test case. Particular attention is given to the impact of both the soil water retention curve and the unsaturated hydraulic conductivity function on model results, especially near saturation.

A Lattice Boltzmann model with two relaxation times for the 2D/3D advection and anisotropic diffusion equation (AADE) is introduced. The method is applied to Richards' equation for variably saturated flow in isotropic homogeneous media by extending retention curves into the saturated zone in a linear manner. The Darcy velocity is computed locally from the population solution. The method possesses intrinsic mass conservation, it is explicit and especially suitable for parallel computations. Designed for regular grids, the LB approach meets the boundary conditions accurately with an unified “multi-reflexion” technique, introduced to fit pressure head and/or specified flux conditions on static and seepage boundaries. The physical space can assume an uniform rectangular discretization grid which is transformed into the cubic computational grid after proper rescaling [9] of the AADE. The diffusion term is considered in two forms: the conventional one and the transformed one. The integral transformation may avoid problems encountered with the unbounded diffusion coefficients at the residual and saturated limits. An analytical expression for the transformed diffusion function is obtained for the original and modified VGM retention curves [14]. Analytical instationary solutions for constant flux infiltration with non-linear models [2,15] are revised. An exact unstationary solution [1] valid for unsaturated, saturated or variably saturated flow is constructed using the BCM hydraulic conductivity function [3,11], moisture tension being fixed at the surface. Stationary infiltration profiles are generated for the BCM and the VGM conductivity functions. The LB method is validated against these and other reference solutions.

International audience; A generalized lattice Boltzmann model to simulate free-surface is constructed in both two and three dimensions. The proposed model satisfies the interfacial boundary conditions accurately. A distinctive feature of the model is that the collision processes is carried out only on the points occupied partially or fully by the fluid. To maintain a sharp interfacial front, the method includes an anti-diffusion algorithm. The unknown distribution functions at the interfacial region are constructed according to the first-order Chapman-Enskog analysis. The interfacial boundary conditions are satisfied exactly by the coefficients in the Chapman-Enskog expansion. The distribution functions are naturally expressed in the local interfacial coordinates. The macroscopic quantities at the interface are extracted from the least-square solutions of a locally linearized system obtained from the known distribution functions. The proposed method does not require any geometric front construction and is robust for any interfacial topology. Simulation results of realistic filling process are presented: rectangular cavity in two dimensions and Hammer box, Campbell box, Sheffield box, and Motorblock in three dimensions. To enhance the stability at high Reynolds numbers, various upwind-type schemes are developed. Free-slip and no-slip boundary conditions are also discussed.

We present a general framework for several previously introduced boundary conditions for lattice Boltzmann models, such as the bounce-back rule and the linear and quadratic interpolations. The objectives are twofold: first to give theoretical tools to study the existing link-type boundary conditions and their corresponding accuracy; second to design boundary conditions for general flows which are third-order kinetic accurate. Using these new boundary conditions, Couette and Poiseuille flows are exact solutions of the lattice Boltzmann models for a Reynolds number Re=0 (Stokes limit) for arbitrary inclination with the lattice directions. Numerical comparisons are given for Stokes flows in periodic arrays of spheres and cylinders, linear periodic array of cylinders between moving plates, and for Navier-Stokes flows in periodic arrays of cylinders for Re

International audience; We prove for generic steady solutions of the Lattice Boltzmann (LB) models that the variation of the numerical errors is set by specific combinations (called ``magic numbers') of the relaxation rates associated with the symmetric and anti-symmetric collision moments. Given the governing dimensionless physical parameters, such as the Reynolds or Peclet numbers, and the geometry of the computational mesh, the numerical errors remain the same for any change of the transport coefficients only when the ``free' (``kinetic') antisymmetric rates and the boundary rules are chosen properly. The single-relaxation-time (BGK) model has no free collision rate and yields viscosity dependent errors with any boundary scheme for hydrodynamic problems. The simplest and most efficient collision operator for invariant errors is the two-relaxation-times (TRT) model. As an example, this model is able to compute viscosity independent permeabilities for any porous structure. These properties are derived from steady recurrence equations, obtained through linear combinations of the LB evolution equations, in which the equilibrium and non-equilibrium components are directly interconnected via finite-difference link-wise central operators. The explicit dependency of the non-equilibrium solution on the relaxation rates is then obtained. This allows us, first, to confirm the governing role of the ``magic' combinations for steady solutions of the Stokes equation, second, to extend this property to steady solutions of the NavierStokes and anisotropic advectiondiffusion equations, third, to develop a parametrization analysis of the microscopic and macroscopic closure relations prescribed via link-wise boundary schemes.

[Departement_IRSTEA]RE [TR1_IRSTEA]RIE / PHYLEAU; International audience; For simple hydrodynamic solutions, where the pressure and the velocity are polynomial functions of the coordinates, exact microscopic solutions are constructed for the two-relaxation-time (TRT) Lattice Boltzmann model with variable forcing and supported by exact boundary schemes. We show how simple numerical and analytical solutions can be interrelated for Dirichlet velocity, pressure and mixed (pressure/tangential velocity) multi-reflection (MR) type schemes. Special care is taken to adapt themfor corners, to examine the uniqueness of the obtained steady solutions and staggered invariants, to validate their exact parametrization by the non-dimensional hydrodynamic and a kinetic (collision) number. We also present an inlet/outlet constant mass flux condition. We show, both analytically and numerically, that the kinetic boundary schemes may result in the appearance of Knudsen layers which are beyond the methodology of the Chapman-Enskog analysis. Time dependent Dirichlet boundary conditions are investigated for pulsatile flow driven by an oscillating pressure drop or forcing. Analytical approximations are constructed in order to extend the pulsatile solution for compressible regimes.