Regularized characteristic boundary conditions for the Lattice-Boltzmann methods at high Reynolds number flows
Gauthier Wissocq Nicolas Gourdain Orestis Malaspinas Alexandre Eyssartier








ISAE,Altran,Cerfacs]Gauthier Wissocq
ISAE]Nicolas Gourdain
UGEN]Orestis Malaspinas
Altran]Alexandre Eyssartier

[ISAE]ISAE, Dpt. of Aerodynamics, Energetics and Propulsion, Toulouse, France
[Altran]Altran, DO ME, Blagnac, France
[Cerfacs]Centre Européen de Recherche et de Formation Avancée en Calcul Scientifique (CERFACS), CFD Team, 42 avenue Gaspard Coriolis, 31057 Toulouse Cedex 01, France
[UGEN]SPC - Centre Universitaire d'Informatique, Université de Genève 7, route de Drize, CH-1227 Switzerland



This paper reports the investigations done to adapt the Characteristic Boundary Conditions (CBC) to the Lattice-Boltzmann formalism for high Reynolds number applications. Three CBC formalisms are implemented and tested in an open source LBM code: the baseline local one-dimension inviscid (BL-LODI) approach, its extension including the effects of the transverse terms (CBC-2D) and a local streamline approach in which the problem is reformulated in the incident wave framework (LS-LODI). Then all implementations of the CBC methods are tested for a variety of test cases, ranging from canonical problems (such as 2D plane and spherical waves and 2D vortices) to a 2D NACA profile at high Reynolds number (Re = 10^5), representative of aeronautic applications. The LS-LODI approach provides the best results for pure acoustics waves (plane and spherical waves). However, it is not well suited to the outflow of a convected vortex for which the CBC-2D associated with a relaxation on density and transverse waves provides the best results. As regards numerical stability, a regularized adaptation is necessary to simulate high Reynolds number flows. The so-called regularized FD (Finite Difference) adaptation, a modified regularized approach where the off-equilibrium part of the stress tensor is computed thanks to a finite difference scheme, is the only tested adaptation that can handle the high Reynolds computation.




Lattice Boltzmann method characteristic boundary conditions LODI high Reynolds number flows












§ INTRODUCTION

A better understanding of turbulent unsteady flows is a necessary step towards a breakthrough in the design of modern aircraft and propulsive systems. Due to the difficulty of predicting turbulence with complex geometry, the flow that develops in these engines remains difficult to predict. At this time, the most popular method to model the effect of turbulence is still the Reynolds Averaged Navier-Stokes (RANS) approach. However there is some evidence that this formalism is not accurate enough, especially when a description of time-dependent turbulent flows is desired (high incidence angle, laminar-to-turbulent transition, etc.) <cit.>. With the increase in computing power, Large Eddy Simulation (LES) applied to the Navier-Stokes equations emerges as a promising technique to improve both knowledge of complex physics and reliability of flow solver predictions <cit.>. It is still the most popular and mature approach to describe the behavior of turbulent flow in complex geometries (e.g. aircraft and gas turbines). However, the resolution of the NS equations requires to add artificial dissipation to ensure numerical stability <cit.>. The consequence is an over-dissipation which affects the flow and limits the capability to transport flow patterns (like turbulence) on a long distance. In some specific cases, like aero-acoustic (far-field noise), NS can thus face some difficulties to predict the flow.

In this context, there is an increasing interest in the fluid dynamics community for emerging methods, based on the Lattice Boltzmann approach <cit.>. The Lattice-Boltzmann Method (LBM) has already demonstrated its potential for complex geometries, thanks to immersed boundary conditions (that allow the use of cartesian grids) and low dissipation properties required for capturing the small acoustic pressure fluctuations <cit.>. LBM also provides the advantage of an easy parallelization, making it well suited for High-Performance Computing <cit.>. However, the most widely used Lattice-Boltzmann models still suffer from weaknesses like a lack of robustness for high Mach number flows (M > 0.4), a limitation to low compressible isothermal flows <cit.> and the use of artificial boundary conditions (Dirichlet/Neumann types can lead to the reflection of outgoing acoustic waves that have a significant influence on the flow field <cit.>). While the use of artificial boundary conditions is also critical for NS methods, it is more problematic for LBM due to the low dissipation of the method.

A potential way to avoid unphysical acoustic reflections at the boundary is to use a "sponge layer" inside the computational domain, on which artificial dissipation (by upwinding) is introduced or physical viscosity is increased (viscosity sponge zones). Acoustic waves (physical or not) are thus damped in such a zone, which allows to eliminate or limit numerical reflections <cit.>. This solution has however important drawbacks. First the calibration of sponge layers is difficult, as a balance must be found between a brutal increase of the viscosity (that will generate acoustic reflections) and a too low dissipation that will not be effective. Such sponge layers also have an impact on the computational cost, since a part of the domain is dedicated to slowly increasing the viscosity. Last, some boundary conditions cannot be treated with a sponge layer, for instance an inlet with turbulence injection.

For NS methods, a successful approach is the use of non-reflective boundary conditions based on a treatment of the characteristic waves of the local flow <cit.>. However, the extension of this approach to LBM is not straightforward given the difficulty to find a bridge between the LBM that describes the world at the mesoscopic level (a population of particles) and the NS world based on a macroscopic description of the flow. But some progress has recently been made on the adaptation of characteristic boundary conditions to the LBM formalism. Izquierdo and Fueyo <cit.> used a pressure antibounceback boundary condition <cit.> adapted to the multiple relaxation time (MRT) collision scheme <cit.> to impose the Dirichlet density and velocity conditions given by the local one-dimensional inviscid (LODI) equations, which provided non-reflective outflow boundary conditions for one-dimensional waves. More recently, Jung et al. <cit.> extended the previous work to include the effects of transverse and viscous terms in the Characteristic Boundary Conditions (CBC) and showed good performance for vortex outflow. Meanwhile, Heubes et al. <cit.> adapted the solution given by a modified-Thompson approach by imposing the corresponding equilibrium populations and Schlaffer <cit.> assessed a modified Zou/He boundary condition <cit.>. 

Still, previous researches are limited to low Reynolds number applications while Latt et al. showed that increasing the Reynolds number can have drastic impact on the numerical stability of the LBM boundary condition <cit.>. Even when the MRT collision is used, as in <cit.>, the numerical stability of the characteristic boundary conditions has not been demonstrated. The aim of this study is thus to develop a numerically stable adaptation of the CBC to the LBM formalism taking advantage of the regularized collision scheme <cit.>, which has proved to be numerically stable at high Reynolds number flows <cit.>, and the corresponding regularized boundary conditions <cit.>. Different kinds of CBC will also be evaluated. This article is structured as follows. The first section describes the LBM framework. Then, the second section presents three kinds of CBCs and three possible adaptations to the LBM formalism for 2D problems in the low compressible isothermal case. In the third section, these models are assessed for simple cases: normal, oblique and spherical waves and a convected vortex at Re=10^3. Finally, the method is assessed on a high Reynolds number application: a NACA0015 airfoil at Re=10^5.



§ NUMERICAL METHOD AND GOVERNING EQUATIONS




 §.§ Lattice Boltzmann framework for isothermal flows

A description of the Lattice-Boltzmann Method can be found in <cit.>. The governing equations describe the evolution of the probability density of finding a set of particles with a given microscopic velocity at a given location: 


    f_i(𝐱+𝐜_𝐢Δ t , t+Δ t) = f_i(𝐱,t) + Ω_i (𝐱,t)

for 0≤ i < q, where 𝐜_𝐢 is a discrete set of q velocities, f_i(𝐱,t) is the discrete single particle distribution function corresponding to 𝐜_𝐢 and Ω_i is an operator representing the internal collisions of pairs of particles.

Macroscopic values such as density, ρ, and the flow velocity, 𝐮, can be deduced from the set of probability density functions f_i(𝐱,t), such as:

    ρ = ∑_i=0^q-1 f_i,      ρ𝐮 =∑_i=0^q-1 f_i 𝐜_𝐢.


Some of the most popular choices for the set of velocities are D2Q9 and D3Q27 lattices, respectively 9 velocities in 2D and 27 velocities in 3D (see Fig. <ref>). For both of these lattices, the sound speed in lattice units (normalized by the ratio between the spatial resolution and the time step Δ x/Δ t) is given by c_s = 1/√(3) <cit.>.







The collision operator Ω_i is usually modelled with the Bhatnagar-Gross-Krook (BGK) approximation <cit.>, which consists in a relaxation, with a relaxation time τ, of every population to the corresponding equilibrium probability density function f_i^(eq):

    Ω_i = -1/τ[f_i(x,t)-f_i^(eq)(x,t)].


The equilibrium distribution function f_i^(eq) is a local function that only depends on density and velocity in the isothermal case. It can be computed thanks to a second order development of the Maxwell-Boltzmann equilibrium function <cit.>:

    f_i^(eq) = w_i ρ[ 1+𝐜_𝐢·𝐮/c_s^2 + ( 𝐜_𝐢·𝐮/2c_s^2)^2- 𝐮^2/2c_s^2],

where w_i are the gaussian weights of the lattice.

A Chapman-Enskog expansion, based on the assumption that f_i is given by the sum of the equilibrium distribution plus a small perturbation f^(1)_i

    f_i=f_i^(eq)+f_i^(1), with  f_i^(1)≪ f_i^(eq),

can be applied to (<ref>) in order to recover the exact Navier-Stokes equation for quasi-incompressible flows in the limit of long-wavelength <cit.>. The pressure is thus given by p=c_s^2 ρ and the kinematic viscosity is linked to the BGK relaxation parameter through

    ν = c_s^2(τ - 1/2).

The Chapman–Enskog expansion also relates the second order tensor Π^(1) defined as

    Π^(1)=∑_i𝐜_𝐢𝐜_𝐢f^(1)_i,

with the strain rate tensor 𝐒=(∇𝐮+(∇𝐮)^T)/2 through the relation

    Π^(1)=-2 c_s^2 ρτ𝐒.

In turn to the leading order f_i^(1) can be approximated by

    f_i^(1)≅w_i/2 c_s^2𝐐_i:Π^(1),

where 𝐐_i≡(𝐜_i𝐜_i-c_s^2𝐈). The colon symbol stands for the double contraction operator and 𝐈 is the identity matrix. A regularization step consisting in the reconstruction of the off-equilibrium parts using (<ref>) and (<ref>) can improve the precision and the numerical stability of the single relaxation time BGK collision <cit.>. This model will be used in the next parts for high Reynolds number flow simulations.



 §.§ The Palabos open-source library

The LBM flow solver used in this work is the Palabos[Copyright 2011-2012 FlowKit Ltd.] open-source library. The Palabos library is a framework for general-purpose CFD with a kernel based on the lattice Boltzmann method. The use of C++ code makes it easy to install and to run on every machine. It is thus possible to set up fluid flow simulations with relative ease and to extend the open-source library with new methods and models, which is of paramount importance for the implementation of new boundary conditions. The numerical scheme is divided in two steps:


  * A collision step where the BGK model is applied:

    f_i(𝐱, t+1/2) = f_i(𝐱,t) + 1/τ[f_i^(eq)(𝐱,t)-f_i(𝐱,t)],

with f_i^(eq) computed using the macroscopic values at time t and f_i can be regularized in order to increase numerical stability for high Reynolds number flows.

  * A streaming step:

    f_i(𝐱+𝐜_𝐢, t+1) = f_i(𝐱, t+1/2).



The streaming step consists in an advection of each discrete population to the neighbor node located in the direction of the corresponding discrete velocity. Since a boundary node has less neighbors than an internal node (less than 9 neighbors in 2D or 27 neighbors in 3D), some populations are missing at the boundary after each iteration. These populations need to be reconstructed, which is the purpose of the implementation of boundary conditions in LBM. Up to now, different methods can be used in Palabos, such as regularized BC <cit.>

 or Zou/He BC <cit.>

to implement open boundaries. However, none of them can be used as they stand for an outflow boundary condition and the use of sponge zones is necessary to avoid non physical reflections. The next sections will aim at developing a more natural boundary condition that minimize acoustic reflections for an outflow boundary type, based on the Characteristic Boundary Conditions (CBC).



§ ADAPTATION OF CHARACTERISTIC BOUNDARY CONDITIONS TO THE LBM FORMALISM


One of the most popular methods in the NS community for subsonic non reflective outflow boundary conditions is the CBC method <cit.>. The adaptation of the CBC to the LBM formalism is presented here for an isothermal flow, in lattice units (normalized by Δ x and Δ t). Acoustic waves thus propagate at the constant lattice sound speed c_s. In the isothermal case, pressure is defined as p=c_s^2 ρ. Three CBC methods will be introduced below: the local one-dimensional inviscid (LODI) approximation, a 2D extension of the LODI approximation including transverse waves and a last method called local-streamline LODI.



 §.§ LODI Approximation (Baseline LODI)


Let us consider a domain outlet located at x=L, as descripted on Fig. <ref>. A diagonalisation of the x-derivative terms in the Navier-Stokes equation allows to define five waves ℒ_i that propagate respectively at velocity u-c_s, u, u, u and u+c_s, where u is the x-component (streamwise) of the non-dimensional macroscopic velocity 𝐮= [u, v, w]. These waves are represented on Fig. <ref> on the inlet (x=0) and outlet (x=L) of a computational domain.







At the outlet (x=L on Fig. <ref>), ℒ_2, ℒ_3, ℒ_4 and ℒ_5 leave the computational domain and are obtained with the general expression of characteristic waves:

    ℒ_2=u( c_s^2 ∂ρ/∂ x - ∂ p/∂ x)=0,


    ℒ_3=u ∂ v/∂ x,


    ℒ_4=u ∂ w/∂ x,


    ℒ_5=(u+c) ( ∂ p/∂ x + ρ c_s ∂ u/∂ x).

Let us notice that ℒ_2, the entropy wave, is null in the isothermal case.The x-derivative terms can be computed using the interior points by one-sided finite difference. The treatment of ℒ_1 is different : since it comes from the outside, it can not be computed using the interior points. The perfectly non-reflecting case is obtained by fixing ℒ_1=0, which ensures eliminating the incoming wave. However, this is known to be unstable because of lack of control of the outlet flow variables. A simple way to ensure well-posedness is to set

    ℒ_1=K_1(p-p_∞),

where K_1=σ (1-M^2)c_s/L, p_∞ is the target pressure at the outlet, σ is a constant, M is the maximum Mach number in the flow and L is a characteristic size of the domain <cit.>.

The time-derivative of the primitive variables can be computed in function of the wave amplitudes by examining a LODI problem :

    ∂ρ/∂ t + 1/c_s^2[ ℒ_2 + 1/2(ℒ_5+ℒ_1) ]=0,


    ∂ p/∂ t + 1/2(ℒ_5+ℒ_1)=0,


    ∂ u/∂ t+1/2ρ c_s(ℒ_5-ℒ_1)=0,


    ∂ v/∂ t + ℒ_3 = 0,


    ∂ w/∂ t+ℒ_4=0.

In the isothermal case, (<ref>) and (<ref>) are equivalent. Finally, with a temporal discretization using an explicit second-order scheme, the physical values that must be imposed at the next time step in order to avoid acoustic reflections can be computed. This implementation will be referred to as the baseline LODI (BL-LODI) in the rest of the paper.



 §.§ LODI approximation including transverse terms


The previous relations are perfectly non-reflecting for the only case of a pure 1D plane wave. For a non-normal wave, the LODI approximation is not verified and a reflected wave, all the more important as the incidence increases, can appear. To take this phenomenon into account, a possible solution is to add the influence of transverse waves in the LODI equations <cit.>:

    ∂ρ/∂ t + 1/2c_s^2(ℒ_5+ℒ_1)=1/2c_s^2(𝒯_5+𝒯_1),


    ∂ u/∂ t+1/2ρ c_s(ℒ_5-ℒ_1)=1/2ρ c_s(𝒯_5-𝒯_1),


    ∂ v/∂ t + ℒ_3 = 𝒯_3,


    ∂ w/∂ t+ℒ_4=𝒯_4.

The transverse waves can be computed as follows:

    𝒯_1=-[ 𝐮_𝐭·∇_tp +p∇_t·𝐮_𝐭 -ρ c_s 𝐮_𝐭·∇_tu ],


    𝒯_3=-[ 𝐮_𝐭·∇_tv + 1/ρ∂ p/∂ y],


    𝒯_4=-[ 𝐮_𝐭·∇_tw + 1/ρ∂ p/∂ z],


    𝒯_5=-[ 𝐮_𝐭·∇_tp +p∇_t·𝐮_𝐭 +ρ c_s 𝐮_𝐭·∇_tu ],

where 𝐮_𝐭=[v,w] and ∇_t=[∂_y, ∂_z].
The second transverse wave 𝒯_2 is not introduced here since it is null in the isothermal case, as well as ℒ_2. The non reflective outflow boundary condition needs now to be set as:

    ℒ_1=K_1(p-p_∞) - K_2(𝒯_1-𝒯_1,exact) + 𝒯_1,

with K_1=σ (1-M^2)c_s/L, K_2 should be equal to the Mach number of the mean flow and 𝒯_1,exact is a desired steady value of 𝒯_1  <cit.>. In the rest of the paper, this method will be named 2D-CBC in the perfectly non-reflecting case (K_1=K_2=0) and 2D-CBC relaxed when a relaxation is done on ℒ_1.



 §.§ Local streamline LODI




Another potential solution is to compute the LODI equation in the local streamline based frame R (Fig. <ref>) <cit.>. In order to compute the new characteristic waves L_i, the non-dimensional velocity vector is projected into the new reference frame with the difficulty to compute x̃-derivative terms from the lattice discretization. A simple approximation is to set 

    ∂ϕ̃/∂x̃ = ∂ϕ̃/∂ x,

and thus to compute it by a first-order upwind scheme using the lattice discretization. This implementation of the CBC condition will be referred to as the local streamline LODI (LS-LODI).



 §.§ Adaptation to the Lattice Boltzmann scheme


The main difficulty in LBM is then to find a set of populations in order to impose the physical values obtained by the CBC theory and the correct associated gradients. The possible adaptations can be divided in two families: those preserving the known particle populations (e.g. Zou/He BC <cit.>) and those replacing all particle populations (e.g. Regularized BC <cit.>). Izquierdo and Fueyo <cit.> and Jung et al. <cit.> chose to modify the missing populations by adapting a pressure anti-bounceback boundary condition, which can be used as far as the MRT collision scheme is adopted <cit.>. Heubes et al. <cit.> decided to impose the CBC physical values thanks to the equilibrium populations, which is known to impose incorrect gradients at the boundary <cit.>. Three possibilities are introduced below and will be further evaluated: a modified Zou/He method and two declinations of the more stable regularized adaptation.






  §.§.§ Adaptation with Zou/He boundary conditions

Let us consider an outflow boundary located at x=L  on a D2Q9 lattice (Fig. <ref>). After streaming, three incoming populations are missing : f_1, f_2 and f_3. 
The zeroth and first-order hydrodynamic moments are:

    ρ = f_1+f_2+f_3+ρ_0+ρ_+,


    ρ u = ρ_+ - (f_1+f_2+f_3),

where ρ_+ = f_5+f_6+f_7 and ρ_0=f_0+f_4+f_8. These equations can be combined, by eliminating the missing populations (f_1+f_2+f_3), to obtain

    ρ = 1/1-u(ρ_0+2ρ_+),

where ρ and u are still non-dimensional. Thus, ρ and u are linked with relation (<ref>), which proves that it is impossible to impose both ρ_b and 𝐮_𝐛 at the boundary without modifying any of the known populations.The same relation can be obtained for every 1D, 2D or 3D lattice as far as there is only one level of velocity. This relation concerns every boundary condition of the first family (those preserving the knwon populations) and will be used later. 

Let us note g_i the corrected populations at the boundary. As proposed in <cit.>, the missing populations can be computed as follows:

    g_1=f_1^(eq) + f_5^(neq) + 1/2(f_4^(neq) - f_8^(neq)),


    g_2=f_2^(eq) + f_6^(neq),


    g_3=f_3^(eq) + f_7^(neq) - 1/2(f_4^(neq) - f_8^(neq)),


while every other populations are kept unchanged: 

    g_i = f_i,        i=0, 4, 5, 6, 7, 8,


with f_i^(neq)=f_i - f_i^(eq) and where the equilibrium populations are computed with the physical values imposed at the boundary, ρ_b and 𝐮_𝐛=(u_b, v_b). 
As shown in <cit.>, corrections (<ref>), (<ref>) and (<ref>) then allow to impose the first-order moment at the boundary ρ_b 𝐮_𝐛. However, density and velocity are still linked with (<ref>). This is not a problem for a dirichlet Zou/He boundary condition where only one physical value (either ρ or 𝐮) is imposed and the other one is computed with (<ref>). On the contrary, for a non reflective outflow, both of them need to be imposed as the result given by the CBC method, so that the only condition set by the user is the value of a characteristic wave. It is then necessary to modify at least one known population. Schlaffer suggests correcting every populations in order to impose the correct density <cit.>. The solution proposed here is to add a correction on the population associated to a null velocity only:

    g_0=f_0 + ρ_b - 1/1-u(ρ_0 + 2ρ_+),

which ensures the value of ρ_b. This choice is motivated by the fact that this added correction will only affect the collision phase and will not be streamed into the computational domain. 

To sum up this method, g_1, g_2, g_3 and g_0 are computed with respectively (<ref>), (<ref>), (<ref>) and (<ref>) while other populations are kept unchanged after streaming:

    g_i=f_i,        i=4, 5, 6, 7, 8.


This method can be easily transposed on every 3D lattice (except for high order lattices <cit.>) by using the general formula of the Zou/He boundary conditions that can be found in <cit.> for the missing populations, and correction (<ref>) to ensure the value of ρ_b and consequently 𝐮_𝐛.



  §.§.§ Adaptation with the regularized method


The Zou/He boundary condition provides the advantage of a very good precision in the definition of the boundary physical values. Unfortunately, this solution suffers from lack of stability for large Reynolds numbers. For example, it has been shown that Zou/He boundary conditions become unstable at Re>100 for a given resolution N=200 nodes per characteristic length in a 2D channel flow <cit.>. 
Another possible adaptation is to use the regularized boundary condition in order to impose the physical values computed by the CBC theory. This solution is yet less accurate but well more stable, as shown by Latt et al. <cit.>.

More details about the regularized method for boundary conditions can be found for example in <cit.>. The purpose
of this section is to explain how this particular boundary condition 
is used to impose ρ_b and 𝐮_b on a flat boundary. 

The leading order of the populations f_i can be expressed (see end of (<ref>) and (<ref>)) as

    f_i=f_i^(eq)(ρ_b,𝐮_𝐛)+f_i^(1)(Π^(1)).

On a boundary node the density ρ_b and velocity 𝐮_b are imposed. Therefore in order 
to be able to use this last equation one needs a way to compute Π^(1). This is achieved
by using the fact that 𝐐_i is a symmetric tensor with respect to i which means that
𝐐_i=𝐐_opp(i), where 

    opp(i)={j|𝐜_i=-𝐜_j}.

At the leading order, this property leads to 

    f^(1)_i=f^(1)_opp(i).

The known f_i^(1) can be straightforwardly computed by the following formula

    f^(1)_i=f_i-f_i^(eq)(ρ_b,𝐮_b).

With the last two equations, the set of f_i^(1) is complete (they are all known) and can be used to compute Π^(1) through (<ref>).
Then using (<ref>) one recomputes regularized f_i^(1) populations and the total populations f_i are computed with the relation (<ref>). This method is valid in both 2D and 3D and will be called Regularized Bounceback (or Regularized BB) adaptation in the next sections. Another possibility is to compute Π^(1) by a second order finite difference scheme thanks to (<ref>) and recompute f^(1) using (<ref>). This method will be called the Regularized FD adaptation.




 §.§ Summary of the method


The non reflecting outflow boundary condition using CBC theory for LBM can be summarized as follows, considering everything is known at (non-dimensional) time t:


  * Computation of the physical values that must be imposed at t+1 to avoid non physical reflections by the CBC theory using either BL-LODI, CBC-2D or LS-LODI method. These values are stored to be used in the last step.
  * Collision step.
  * Streaming step: some populations are missing at the boundary.
  * Correction of the set of populations at the boundary so that the physical values stored in the first step are imposed by using the Zou/He adaptation the so-called regularized Bounceback adaptation or the Regularized FD adapatation. 



§ APPLICATION TO ACADEMIC CASES


In this section, the CBC approach for LBM is assessed on simple 2D cases: a normal plane wave, a plane wave with different incidence angles, a spherical density wave and a convected vortex.



 §.§ 2D normal plane wave

The computational domain, a square of 200×200 cells, is initiated with a gaussian plane wave as follows (in lattice units)

    ρ_0=1+0.1*exp(-(x-x_0)^2/2σ^2),


    u_0=0.1,


    v_0=0.1*exp(-(x-x_0)^2/2R_c^2),

with x_0=110 and 2R_c^2=20 in lattice units (i.e. in number of cells). The Reynolds number, computed with the horizontal non-dimensional initial velocity u_0, the characteristic size of the box in number of voxels and the viscosity in lattice units, is equal to 100.
The boundary conditions are: vertical periodicity, reflecting inlet on the left and perfectly non reflecting outflow on the right. For this pure 1D case, the choice of the CBC condition (baseline, local streamline LODI or LODI with transverse terms) has no effect on the absorption rate. Moreover, all the three adaptations provide the same results at 10^-6 and no stability issues have been encountered with the Zou/He adaptation for this test case. Thus, the results presented below have been obtained with baseline LODI and Zou/He adaptation.
This test case has been chosen so that the reflection rate for every macroscopic value (ρ, u and v) can be computed, since two pressure and axial velocity waves propagates at speed (u-c_s) and (u+c_s) and one transverse velocity wave propagates at speed u. The computed density waves are represented on Fig. <ref> at two different time steps: before reflection of the (u+c_s) wave on the non-reflective outlet and shortly after reflection. 



A very low reflected amplitude can be distinguished. The (u-c_s) wave propagates to the left of the domain without encountering any boundary at the two observed time steps. It is only affected by viscous dissipation and is used as a reference amplitude in the computation of the reflection rates of density and axial velocity. For the transverse velocity wave, one can compute the reflection rate as the ratio between the reflected wave amplitude and the amplitude shortly before reflection. The obtained results are presented on Table <ref>.



As often with CBC, the treatment is more difficult for the pressure wave, but the results obtained here are in good agreement compared to what is found in the literature <cit.>. 

It is noticed that ρ=1 is not correctly recovered at the outlet after reflection, as the boundary has been set as perfectly non-reflecting (ℒ_1=0). In order to impose the correct boundary condition, a relaxation should be implemented, as in eq (<ref>). However, in that case, the reflection rate will increase as shown in <cit.>.



 §.§ 2D plane wave with incidence


At t=0, the computational domain is at rest (ρ_0=1, 𝐮_0=0) except for an oblique line on which the density is set to ρ_0=1.1 in order to generate a plane wave with an incidence α. The reflection coefficient is measured by computing the ratio of maximal amplitudes in density waves only, as this is the most critical hydrodynamic variable. The three implementations are tested for this case: BL-LODI, LS-LODI and CBC-2D in the perfectly non-reflecting case. Again, the results obtained with the Zou/He adaptation only will be presented here, as the same results at 10^-6 have been obtained with the Regularized BB and Regularized DF adaptations.

[][h!]
[width=0.35]ObliqueWaveProblem
(49,68)α
(53,56)M
(63,65)c_s
(65,35)h
(90,75)270 Non-reflective outlet
(5,20)90 Reflective inlet
(40,88)Periodicity
(40,5)Periodicity

Schematic plot of the initialization of the plane wave test case with an incidence angle α. Spherical waves instantaneously appear at the two extremities to distort the plane wave. 


Because of a spurious phenomenon appearing at high incidence angles, only angles below 45 ^∘could be computed. Indeed, contrary to the previous normal wave test case, the incident wave cannot be infinite in its transverse direction. Then, spherical waves appear at each extremity of the oblique line initializing the wave, as shown on Fig. <ref>. Let us imagine a point M located at a distance h of one extremity and moving with the oblique plane wave at velocity c_s. The spherical wave reaches M after a time h/c_s. The point M reaches the non-reflective outlet boundary condition at a time h tan (α)/c_s. The condition for this point of the wave to reach the outlet before being distorted is: 

    h/c_stan(α) < h/c_s⇔tan(α) < 1,

which means than only incidence angles below 45^∘ can be computed with such a test case. This problem is avoided in <cit.> by imposing an 'exact' solution obtained on a larger computational domain until the desired wave reaches the boundary. The exact solution is then switched to the CBC condition. It has not been tested in this paper in order to avoid the eventual acoustic waves generated by a brutal change in the boundary condition.

The reflection coefficient of the density wave with respect to the incidence angle is represented on Fig. <ref> for BL-LODI, LS-LODI and CBC-2D. The results are compared with what is obtained for the same test case with the modified Thompson method with a coefficient γ=3/4 as introduced in <cit.>.



As expected, for the baseline approach (a), the reflection coefficient increases as the incidence angle increases. For the modified Thompson approach, the results are close to what can be found in <cit.>: the reflection rate slightly increases and reaches 5% at 40^∘ of incidence. On the contrary, when the CBC method is extended with the effects of transverse waves as in <cit.>, the reflection rate decreases until 30^∘ and then begins to increase. In the case of a local streamline LODI implementation, the coefficient remains stable at around 2% of reflected wave.



 §.§ 2D Spherical wave


The computational domain, a square of 600×600 cells is initiated with a gaussian density in order to generate a spherical wave, as follows :

    ρ=1+0.1*exp(-((x-x_0)^2+(y-y_0)^2)/2R_c^2),
 
with R_c = 3.2, x_0=520 and y_0=300 nodes.
The boundary conditions are periodic in the vertical direction, a CBC condition is set on the right boundary and a reflective boundary condition on the left (located far enough to avoid its influence on the spherical wave at the studied time steps). Four cases are computed: (a) baseline LODI, (b) CBC with transverse terms, (c) local streamline LODI and (d) reference case where the domain is enlarged in the horizontal direction so that the spherical wave is not affected by the boundary (Fig. <ref>). As for the previous test cases, the LBM adaptation of the CBC condition had no impact on the results.



The local reflection coefficient for such a spherical wave can be computed by the following formula:

    r = R-A_ref/|I|,

where I is the amplitude of the original density wave running towards the boundary condition, R is the amplitude of the density wave after reflection at the outlet and A_ref is the amplitude of the density wave at the same lattice node compared to the reference case (d) (no reflection). 



A map of reflection rates after reflection at the outlet (after 400 iterations) is shown on Fig. <ref>. For the baseline approach, one can see that the reflected wave becomes more important as the local incidence of the spherical waves increases, because the incident wave can be approximated as a local normal plane in the center of the outlet while, in the corners, the incidence becomes important (up to 75^∘). When adding transverse term without any relaxation as in <cit.>, the observation confirms the prospective behavior of the plane waves simulations: even if the reflection rate is consistently reduced for angles below 40^∘, it reaches 20% for the greatest observed angles. On the contrary, the reflection rate remains constant with the local streamline LODI implementation.













 §.§ 2D convected vortex

A 2D vortex is convected from left to right and exits the computational domain, a square of 600×600 grid points. A particular attention must be paid to the initialization of the Lamb-Oseen vortex <cit.> which has to be adapted to the isothermal case in order to avoid spurious waves generated by the adaptation of a wrong initial density. The initial conditions, in lattice units, are imposed as follows:


    u = u_0 - β u_0 (y-y_0)/R_cexp( -r^2/2R_c),


    v = β u_0 (x-x_0)/R_cexp( -r^2/2R_c),


    ρ = [ 1- (β u_0)^2/2C_vexp( -r^2/2) ] ^1/(γ-1),


where u_0=0.1 in lattice units, β=0.5, x_0=y_0=300 nodes (the vortex is initially centered on the box), R_c=20 nodes and r^2=(x-x_0)^2+(y-y_0)^2. With the BGK collision operator with a single relaxation time, the simulated gas has the following constants <cit.>:

    γ=D+2/D,


    C_v=D/2 c_s^2,

where D=2 is the dimension of the problem. The specific heat capacity at constant volume C_v appears instead of C_p in (<ref>) because of an error in the heat flux obtained in the Navier-Stokes equations after the Chapman-Enskog development for athermal lattices <cit.>. The Reynolds number based on u_0 and the size of the computational domain is equal to 1000 and the Regularized BGK scheme is chosen for the collision step. This convected vortex test case is a well known test for boundary conditions as it often reveals spurious distorsions at boundaries due to the presence of transverse terms in the Navier-Stokes equation <cit.>. As previously, top and bottom boundary conditions are periodic, the left condition is a regularized inlet and four different CBC conditions will be evaluated at the right boundary:(a) baseline LODI with K_1=0,

(b) local streamline LODI with K_1=0,

(c) CBC-2D with ℒ_1=𝒯_1,

(d) CBC-2D relaxed including transverse terms with: 
    ℒ_1=σ (1-M^2)c_s^3/R_c(ρ-1) + (1-M)𝒯_1,
 where M=0.2 and σ=0.9.

For low Reynolds numbers (Re < 1000), simulations — not shown here — provided the same results at 10^-6. However, for Re=1000, the Zou/He adaptation was no more stable for this test case, contrary to both regularized methods. The results presented here have been obtained with the Regularized BB adaptation.






Fig.<ref> and <ref> show isovalues of non dimensional longitudinal velocity for the four studied cases. First, it can be noticed that, contrary to what was observed in the previous test cases, the use of local streamline boundary condition does not allow to reduce non physical reflections compared to baseline LODI. A possible explanation would be that, inside the vortex, local streamlines are nearly perpendicular to the direction of propagation of the local wave, whereas they were aligned in the case of a pure acoustic wave. The LODI equations are thus applied in a wrong frame which generates non physical reflections. Results are slightly better for the BL-LODI for which the local frame is the good one for at least local normal waves. The addition of the transverse waves in the CBC-2D boundary allows the vortex to keep a correct shape at least until 1000 iterations. However, once the first half of the vortex has reached the outlet, it is distorted as in the BL-LODI case. The addition of relaxation coefficients allows the vortex to keep its shape at the last iteration, even if it is a bit distorted.
It can be noticed that the configuration appears to be symmetric for boundaries (a) and (c), which is not the case for (b) and (d). Indeed, with the LS-LODI, the streamlines used in the change of frame are not symmetric and in the last case, the density frame appears to be asymmetric which has an influence on ℒ_1 through the relaxed pressure, and thus on the longitudinal velocity fields. 



 §.§ Stability analysis


An analysis of numerical stability of the previous convected vortex test case at different Reynolds numbers and different grid resolutions has been carried out. It has been noticed that the CBC type (BL-LODI, LS-LODI or CBC-2D) had no impact on the stability as far as ℒ_1 is not relaxed: the numerical stability comes from the LBM adaptation itself.

Fig. <ref> compares the stability of the Zou/He adaptation and the Regularized BB adaptation. As predicted from <cit.>, the classical regularized boundary condition is more stable than the Zou/He adaptation. However, a huge resolution is still required in order to reach high Reynolds numbers. On the contrary, the Regularized FD adaptation has shown to be unconditionnally numerically stable: in every configurations of the convected vortex test case, the first numerical instabilities came from the inside of the domain and not the CBC boundary condition.






§ APPLICATION TO A NACA0015 PROFIL AT HIGH REYNOLDS NUMBER

The objective of this section is to demonstrate the robustness of the CBCs in a case relevant for high Reynolds aerodynamics applications. The configuration is a NACA0015 profile, in a 8C×8C domain (with C=1 m the chord of the profile). The Mach number is set to 0.04 and the Reynolds number is set to Re=10^5. The lattice dimension is Δ x=1/400 m, corresponding to a time step Δ t=4.25×10^-6 s. Each simulation is run for 200 000 time steps, but only the last 100 000 time steps are kept for data post-processing. The simulations are performed on a 3,200×3,200 points 2D grid, to ensure numerical stability and a proper resolution of the flow patterns with a LES formalism. The sub grid scale model is the Smagorinsky model, with a constant C_s=0.18, associated with a Regularized BGK collision scheme. In order to generate complex flow patterns, the profile is inclined with an angle of 15^o, compared to the freestream velocity direction which is purely axial. Three CBCs are tested: BL-LODI, LS-LODI and CBC-2D. The CBC solutions are compared with a reference solution obtained on 16 C×8C domain, Fig. <ref>. In this reference case, vortices do not leave the extended computational domain at the observed time steps and the initial acoustic wave is evacuated thanks to a Neumann outlet (populations at the boundary are copied from the neighbor voxel) associated with a 4C-wide viscosity sponge zone inside which the relaxation time is increased up to Δ t with a sinusoidal shape in order to smoothly increase the numerical dissipation. Among CBC computations, only the regularized FD adaptation has been sufficiently robust to achieve the simulations.

At this Reynolds number, the real flow should be 3D in the vicinity of the profile and in the wake, due to the development of turbulent flow patterns. However, the present 2D approach is not able to reproduce such 3D effects. Despite that, the trajectory of the vortices exhibits a chaotic behavior, as observed on the time-averaged flow field in Fig. <ref>. The challenges with the present test case are to ensure that:


  * the initial acoustic wave generated by the presence of the airfoil leaves the domain with minimum reflection,

  * the vortices generated by the boundary layer separation are correctly convected beyond the outlet plane with minimum spurious wave reflection. 




The difficulty for such a test case lies in the chaotic nature of the flow, since the trajectory of the vortices is affected by small perturbations. Because of this phenomenon, the visualization of error fields, computed as the difference between the reference solution and each tested CBC, would not be compellant since vortices are not superposed in each computation. A better choice is to plot the root mean square of the pressure field, defined as p_RMS=√(p'^2)/p, to underline the behavior of each CBC during the whole computation as in Fig. <ref>. 



It can be noticed that the baseline LODI approach (a) allows to evacuate the initial acoustic wave and the convected vortices with minimal reflection. The only differences with the reference p_RMS field are the small perturbations observed at the outlet and a background noise which can be due to the reflection of the initial acoustic wave. The local streamline LODI approach (b) is less efficient: the background noise is more pronounced and the dissymetry observed in the previous academic test cases is still present, which increases pressure fluctuations at the outlet. In the CBC-2D case (c), the map of pressure fluctuations seems more smooth and the vortices does not seem do produce spurious reflections. But the overall prms is increased in the whole domain, which can be due to a drift of the mean pressure because of the absence of relaxation in the boundary condition.

To quantify more the performance of the three CBC approaches, the time-averaged fluctuations of pressure p_RMS and velocity u_RMS, v_RMS are shown in Fig. <ref> at x/C=7.5 (close to the outlet where the CBC is applied). All CBCs show a good ability to predict the mean flow as well as pressure and velocity fluctuations. This test case remains a challenge for CBCs since the pressure and velocity fluctuations outside the wake are four magnitude orders lower than the mean field value. In that regard, all CBC methods predict pressure and velocity fluctuations of the same magnitude order than in the reference case. The CBC that provides the best results in term of accuracy for both pressure and velocity fluctuations is the BL-LODI approach. The LS-LODI gives satisfactory results outside the wake but it overestimates p_RMS by a factor 2 in the wake (at y/C=4) because of the dissymetry in the vortices reflections. The CBC-2D approach overestimates p_RMS outside the wake by a factor 4 but it gives good results in the wake, similar to those obtain with the BL-LODI approach. Similar conclusions can be drawn for axial and transverse velocity fluctuations, except that all methods predict the correct velocity fluctuations in the wake, including the LS-LODI approach.
 





§  SUMMARY AND CONCLUSION


An implementation of a non reflective outflow boundary condition, which does not need additional absorbing layers nor extended domains, has been proposed for a lattice Boltzmann solver. The methods presented here are based on the Characteristic Boundary Conditions (CBC) with the classical LODI approach (BL-LODI), its extension to transverse waves (CBC-2D) and the LODI approach in the local streamline based frame (LS-LODI). Three ways of computing missing populations in order to impose the CBC physical values have been introduced. The first one is based on the classical Zou/He boundary condition, while the other ones are based on the more stable regularized boundary condition : the so-called "Regularized BB" adaptation, where the off-equilibrium part of the stress tensor is evaluated thanks to a bounceback rule, and the "Regularized FD" adaptation, where it is computed thanks to an upwind finite difference scheme. All these methods provided very good results in the test case of a normal plane wave, where computed reflection rates were about 1%. Testcases of a plane wave with incidence and a spherical wave showed that the reflection rate increases with the incidence angle for the BL-LODI adaptation. When adding the effect of transverse terms, the reflected wave is considerably reduced but begins to increase for incidence angles greater than 30^∘. For these pure acoustic waves, the LS-LODI adaptation provided the best results since the reflection rate remained below 5% whatever the incidence angle. However, this method is not adapted for a convected vortex, for which the CBC-2D adaptation with relaxation on density and transverse waves provided the best results and lead to only slight distortions of the velocity fields. As regards the numerical stability of the implemented CBC, the regularized adaptations have shown to be well more stable than the Zou/He one. The regularized FD adaptation associated with a regularized BGK scheme allowed to run the NACA0015 case at high Reynolds number (Re=10^5) in 2D with N=400 cells per chord length, which was not possible with Zou/He or Regularized BB adaptations. Thus, the use of the regularized FD adaptation is well advised for high Reynolds computations.




§ ACKNOWLEDGEMENTS

The authors are grateful to the Calmip computing center of the Federal University of Toulouse (Project account number P1425) for providing all resources that have been used in this work. The authors would also thank J.F. Boussuge from CERFACS for his help on post-processing and the discussion about the method.




















elsarticle-num
10
url<#>1urlprefixURL href#1#2#2 #1#1

Tucker:2014
P. G. Tucker, J. R. DeBonis, Aerodynamics, computers and the environment,
  Philosophical Transactions of the Royal Society A: Mathematical, Physical and
  Engineering Sciences 372 (2022) (2014) 20130331–20130331.

Chen_AnnuRevFluid_30_1998
S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid fows, Annual
  Review of Fluid Mechanics 30 (1) (1998) 329–364.

Succi_2001
S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and
  Beyond, Numerical Mathematics and Scientific Computation, Clarendon Press,
  2001.

Lallemand_PhysRevE_61_2000
P. Lallemand, L.-S. Luo, Theory of the lattice Boltzmann method:
  Dispersion, dissipation, isotropy, Galilean invariance, and stability,
  Phys. Rev. E 61 (2000) 6546–6562.

Buick_EPL_43_1998
J. M. Buick, C. A. Greated, D. M. Campbell, Lattice BGK simulation of
  sound waves, EPL (Europhysics Letters) 43 (3) (1998) 235.

Marie_JCP_228_2009
S. Marié, D. Ricot, P. Sagaut, Comparison between lattice Boltzmann
  method and Navier-Stokes high order schemes for computational
  aeroacoustics, Journal of Computational Physics 228 (4) (2009) 1056 – 1070.

Heuveline_CMAP_58_2009
V. Heuveline, M. J. Krause, J. Latt, Towards a hybrid parallelization of
  lattice Boltzmann methods, Computers & Mathematics with Applications
  58 (5) (2009) 1071 – 1080, mesoscopic Methods in Engineering and Science.

Colonius_AnnuRevFluid_36_2004
T. Colonius, Modeling artificial boundary conditions for compressible flow,
  Annual Review of Fluid Mechanics 36 (1) (2004) 315–345.

Bodony_JCP_212_2006
D. J. Bodony, Analysis of sponge zones for computational fluid mechanics,
  Journal of Computational Physics 212 (2) (2006) 681 – 702.

Israeli_JCP_41_1981
M. Israeli, S. A. Orszag, Approximation of radiation boundary conditions,
  Journal of Computational Physics 41 (1) (1981) 115 – 135.

Poinsot_JCP_101_1992
T. Poinsot, S. Lele, Boundary conditions for direct simulations of compressible
  viscous flows, Journal of Computational Physics 101 (1) (1992) 104 – 129.

Yoo_CTM_9_2005
C. S. Yoo, Y. Wang, A. Trouvé, H. G. Im, Characteristic boundary conditions
  for direct simulations of turbulent counterflow flames, Combustion Theory and
  Modelling 9 (4) (2005) 617–646.

Lodato_JCP_227_2008
G. Lodato, P. Domingo, L. Vervisch, Three-dimensional boundary conditions for
  direct and large-eddy simulation of compressible viscous flows, Journal of
  Computational Physics 227 (10) (2008) 5105 – 5143.

Izquierdo_PhysRevE_78_2008
S. Izquierdo, N. Fueyo, Characteristic nonreflecting boundary conditions for
  open boundaries in lattice Boltzmann methods, Phys. Rev. E 78 (2008)
  046707.

Ginzburg:2008
I. Ginzburg, F. Verhaeghe, D. d'Humieres, Two-relaxation-time lattice boltzmann
  scheme: About parametrization, velocity, pressure and mixed boundary
  conditions, Comm. Comp. Phys. 3 (2008) 427–478.

Dhumieres:1992
D. D'Humieres, Generalized lattice-Boltzmann equations, Progress in
  Astronautics and Aeronautics (159) (1992) 450–458.

Jung:2015
N. Jung, H. W. Seo, C. S. Yoo, Two-dimensional characteristic boundary
  conditions for open boundaries in the lattice Boltzmann methods, Journal of
  Computational Physics 302 (August) (2015) 191–199.

Heubes_JCAM_262_2014
D. Heubes, A. Bartel, M. Ehrhardt, Characteristic boundary conditions in the
  lattice Boltzmann method for fluid and gas dynamics, Journal of
  Computational and Applied Mathematics 262 (2014) 51 – 61, selected Papers
  from NUMDIFF-13.

Schlaffer:2013
M. B. Schlaffer, Non-reflecting Boundary Conditions for the Lattice Boltzmann
  Method, Ph.D. thesis, Technische Universität Münschen (2013).

Zou_PhysFluids_9_1997
Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice
  Boltzmann BGK model, Physics of Fluids 9 (6) (1997) 1591–1598.

Latt_PhysRevE_77_2008
J. Latt, B. Chopard, O. Malaspinas, M. Deville, A. Michler, Straight velocity
  boundaries in the lattice Boltzmann method, Phys. Rev. E 77 (2008) 056703.

Latt:2006
J. Latt, B. Chopard, Lattice Boltzmann method with regularized pre-collision
  distribution functions, Mathematics and Computers in Simulation 72 (2-6)
  (2006) 165–168.
http://arxiv.org/abs/0506157 arXiv:0506157.

Malaspinas:2015
O. Malaspinas, Increasing stability and accuracy of the lattice Boltzmann
  scheme: recursivity and regularization (2015) 1–31http://arxiv.org/abs/1505.06900 arXiv:1505.06900.

Bhatnaghar_PhysRev_94_1954
P. L. Bhatnagar, E. P. Gross, M. Krook, A Model for Collision Processes
  in Gases. I. Small Amplitude Processes in Charged and Neutral
  One-Component Systems, Phys. Rev. 94 (1954) 511–525.

Qian_EPL_17_1992
Y. H. Qian, D. D'Humières, P. Lallemand, Lattice BGK Models for
  Navier-Stokes Equation, EPL (Europhysics Letters) 17 (6) (1992) 479.

Chapman_1952
S. Chapman, T. Cowling, The mathematical theory of non-uniform gases: an
  account of the kinetic theory of viscosity, thermal conduction, and diffusion
  in gases, no. vol. 2, University Press, 1952.

Malaspinas_CompFluids_49_2011
O. Malaspinas, B. Chopard, J. Latt, General regularized boundary condition for
  multi-speed lattice Boltzmann models, Computers & Fluids 49 (1) (2011) 29
  – 35.

Yoo:2007
C. S. Yoo, H. G. Im, Characteristic boundary conditions for simulations of
  compressible reacting flows with multi-dimensional, viscous and reaction
  effects, Combustion Theory and Modelling 11 (2) (2007) 259–286.

Albin_CompFluids_51_2011
E. Albin, Y. D’Angelo, L. Vervisch, Flow streamline based Navier-Stokes
  Characteristic Boundary Conditions: Modeling for transverse and
  corner outflows, Computers & Fluids 51 (1) (2011) 115 – 126.

Philippi2006
P. C. Philippi, L. A. Hegele, L. O. E. dos Santos, R. Surmas, From the
  continuous to the lattice Boltzmann equation: The discretization problem and
  thermal models, Physical Review E 73 (5) (2006) 056702.

Lamb:1932
H. Lamb, Hydrodynamics, 6th Edition, Cambridge University Press, 1932.

Guo:2007
Z. Guo, C. Zheng, B. Shi, T. S. Zhao, Thermal lattice Boltzmann equation for
  low Mach number flows: Decoupling model, Phys. Rev. E 75 (3) (2007) 1–15.

Yoo_CTM_11_2007
C. S. Yoo, H. G. Im, Characteristic boundary conditions for simulations of
  compressible reacting flows with multi-dimensional, viscous and reaction
  effects, Combustion Theory and Modelling 11 (2) (2007) 259–286.




