Regularized characteristic boundary conditions for the Lattice-Boltzmann methods at high Reynolds number flows

TL;DR

This paper adapts and tests various characteristic boundary conditions within the Lattice-Boltzmann framework for high Reynolds number flows, identifying the most effective methods and necessary regularization for stability.

Contribution

It introduces and compares three CBC formalisms in LBM for high Reynolds number flows, proposing a regularized finite difference approach for numerical stability.

Findings

LS-LODI best for acoustic waves

CBC-2D excels with vortices

Regularized FD approach ensures stability at high Reynolds numbers

Abstract

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.

Figures

Tables

Table 1. Table 1: Reflection rates for the 2D normal plane wave.

Variable	$\rho$	$u$	$v$
Reflection rate	1.2 %	1.1 %	%

Equations

$$f_i(\mathbf{x}+{\mathbf{c}}_{\mathbf{i}}\Delta t,t+\Delta t)=f_i(\mathbf{x},t)+{\Omega }_i(\mathbf{x},t)$$

$$\rho =\sum _{i=0}^{q-1}{f}_i,\rho \mathbf{u}=\sum _{i=0}^{q-1}{f}_i{\mathbf{c}}_{\mathbf{i}}.$$

$${\Omega }_i=-\frac{1}{\tau }\left[f_i(x,t)-f_i^{(eq)}(x,t)\right].$$

$$f_i^{(eq)}=w_i\rho \left[1+\frac{{\mathbf{c}}_{\mathbf{i}}\cdot \mathbf{u}}{c_s^2}+{\left(\frac{{\mathbf{c}}_{\mathbf{i}}\cdot \mathbf{u}}{2c_s^2}\right)}^2-\frac{{\mathbf{u}}^2}{2c_s^2}\right],$$

$$f_i=f_i^{(eq)}+f_i^{(1)},with{f}_i^{(1)}\ll f_i^{(eq)},$$

$$\nu =c_s^2(\tau -\frac{1}{2}).$$

$${\mathbf{\Pi }}^{(1)}=\sum _i{\mathbf{c}}_{\mathbf{i}}{\mathbf{c}}_{\mathbf{i}}{f}_i^{(1)},$$

$${\mathbf{\Pi }}^{(1)}=-2c_s^2\rho \tau \mathbf{S}.$$

$$f_i^{(1)}\cong \frac{w_i}{2c_s^2}{\mathbf{Q}}_i:{\mathbf{\Pi }}^{(1)},$$

$$f_{i}(\mathbf{x},t+\frac{1}{2})=f_{i}(\mathbf{x},t)+\frac{1}{\tau}\left[f_{i}^{(eq)}(\mathbf{x},t)-f_{i}(\mathbf{x},t)\right],\$$

$$f_i(\mathbf{x}+{\mathbf{c}}_{\mathbf{i}},t+1)=f_i(\mathbf{x},t+\frac{1}{2}).$$

$${\mathcal{L}}_2=u\left(c_s^2\frac{\partial \rho }{\partial x}-\frac{\partial p}{\partial x}\right)=0,$$

$${\mathcal{L}}_3=u\frac{\partial v}{\partial x},$$

$${\mathcal{L}}_4=u\frac{\partial w}{\partial x},$$

$${\mathcal{L}}_5=(u+c)\left(\frac{\partial p}{\partial x}+\rho c_s\frac{\partial u}{\partial x}\right).$$

$${\mathcal{L}}_1=K_1(p-p_{\infty }),$$

$$\frac{\partial \rho }{\partial t}+\frac{1}{c_s^2}\left[{\mathcal{L}}_2+\frac{1}{2}({\mathcal{L}}_5+{\mathcal{L}}_1)\right]=0,$$

$$\frac{\partial p}{\partial t}+\frac{1}{2}({\mathcal{L}}_5+{\mathcal{L}}_1)=0,$$

$$\frac{\partial u}{\partial t}+\frac{1}{2\rho c_s}({\mathcal{L}}_5-{\mathcal{L}}_1)=0,$$

$$\frac{\partial v}{\partial t}+{\mathcal{L}}_3=0,$$

$$\frac{\partial w}{\partial t}+{\mathcal{L}}_4=0.$$

$$\frac{\partial \rho }{\partial t}+\frac{1}{2c_s^2}({\mathcal{L}}_5+{\mathcal{L}}_1)=\frac{1}{2c_s^2}({\mathcal{T}}_5+{\mathcal{T}}_1),$$

$$\frac{\partial u}{\partial t}+\frac{1}{2\rho c_s}({\mathcal{L}}_5-{\mathcal{L}}_1)=\frac{1}{2\rho c_s}({\mathcal{T}}_5-{\mathcal{T}}_1),$$

$$\frac{\partial v}{\partial t}+{\mathcal{L}}_3={\mathcal{T}}_3,$$

$$\frac{\partial w}{\partial t}+{\mathcal{L}}_4={\mathcal{T}}_4.$$

$${\mathcal{T}}_1=-\left[{\mathbf{u}}_{\mathbf{t}}\cdot {\nabla }_{\mathbf{t}}p+p{\nabla }_{\mathbf{t}}\cdot {\mathbf{u}}_{\mathbf{t}}-\rho c_s{\mathbf{u}}_{\mathbf{t}}\cdot {\nabla }_{\mathbf{t}}u\right],$$

$${\mathcal{T}}_3=-\left[{\mathbf{u}}_{\mathbf{t}}\cdot {\nabla }_{\mathbf{t}}v+\frac{1}{\rho }\frac{\partial p}{\partial y}\right],$$

$${\mathcal{T}}_4=-\left[{\mathbf{u}}_{\mathbf{t}}\cdot {\nabla }_{\mathbf{t}}w+\frac{1}{\rho }\frac{\partial p}{\partial z}\right],$$

$${\mathcal{T}}_5=-\left[{\mathbf{u}}_{\mathbf{t}}\cdot {\nabla }_{\mathbf{t}}p+p{\nabla }_{\mathbf{t}}\cdot {\mathbf{u}}_{\mathbf{t}}+\rho c_s{\mathbf{u}}_{\mathbf{t}}\cdot {\nabla }_{\mathbf{t}}u\right],$$

$${\mathcal{L}}_1=K_1(p-p_{\infty })-K_2({\mathcal{T}}_1-{\mathcal{T}}_{1,exact})+{\mathcal{T}}_1,$$

Full text

††journal: Journal of Computational Physics

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.) [1]. 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 [1]. 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 [1]. 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 [2, 3, 4]. 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 [5, 6]. LBM also provides the advantage of an easy parallelization, making it well suited for High-Performance Computing [7]. 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 [2] 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 [8, 9, 10]). 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 [9]. 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 [11, 12, 13]. 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 [14] used a pressure antibounceback boundary condition [15] adapted to the multiple relaxation time (MRT) collision scheme [16] 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. [17] 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. [18] adapted the solution given by a modified-Thompson approach by imposing the corresponding equilibrium populations and Schlaffer [19] assessed a modified Zou/He boundary condition [20].

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 [21]. Even when the MRT collision is used, as in [14, 15], 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 [22], which has proved to be numerically stable at high Reynolds number flows [23], and the corresponding regularized boundary conditions [21]. 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}$.

1 Numerical method and governing equations

Lattice Boltzmann framework for isothermal flows

A description of the Lattice-Boltzmann Method can be found in [2, 3, 4]. 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:

[TABLE]

for $0\leq i<q$, where $\mathbf{c_{i}}$ is a discrete set of $q$ velocities, $f_{i}(\mathbf{x},t)$ is the discrete single particle distribution function corresponding to $\mathbf{c_{i}}$ and $\Omega_{i}$ is an operator representing the internal collisions of pairs of particles.

Macroscopic values such as density, $\rho$, and the flow velocity, $\mathbf{u}$, can be deduced from the set of probability density functions $f_{i}(\mathbf{x},t)$, such as:

[TABLE]

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. 1). For both of these lattices, the sound speed in lattice units (normalized by the ratio between the spatial resolution and the time step $\Delta x/\Delta t$) is given by $c_{s}=1/\sqrt{3}$ [3].

The collision operator $\Omega_{i}$ is usually modelled with the Bhatnagar-Gross-Krook (BGK) approximation [24], which consists in a relaxation, with a relaxation time $\tau$, of every population to the corresponding equilibrium probability density function $f_{i}^{(eq)}$:

[TABLE]

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 [25]:

[TABLE]

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}$

[TABLE]

can be applied to (1) in order to recover the exact Navier-Stokes equation for quasi-incompressible flows in the limit of long-wavelength [26]. The pressure is thus given by $p=c_{s}^{2}\rho$ and the kinematic viscosity is linked to the BGK relaxation parameter through

[TABLE]

The Chapman–Enskog expansion also relates the second order tensor $\bm{\Pi^{(1)}}$ defined as

[TABLE]

with the strain rate tensor $\mathbf{S}=(\bm{\nabla}\mathbf{u}+(\bm{\nabla}\mathbf{u})^{T})/2$ through the relation

[TABLE]

In turn to the leading order $f_{i}^{(1)}$ can be approximated by

[TABLE]

where $\mathbf{Q}_{i}\equiv\left(\mathbf{c}_{i}\mathbf{c}_{i}-c_{s}^{2}\mathbf{I}\right)$. The colon symbol stands for the double contraction operator and $\mathbf{I}$ is the identity matrix.

A regularization step consisting in the reconstruction of the off-equilibrium parts using (9) and (7) can improve the precision and the numerical stability of the single relaxation time BGK collision [22, 23]. 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 Palabos111Copyright 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:

[TABLE]

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.

• 2.

A streaming step:

[TABLE]

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 [27] or Zou/He BC [20] 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).

2 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 [11]. The adaptation of the CBC to the LBM formalism is presented here for an isothermal flow, in lattice units (normalized by $\Delta x$ and $\Delta 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}\rho$. 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. 2. A diagonalisation of the x-derivative terms in the Navier-Stokes equation allows to define five waves $\mathcal{L}_{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 $\mathbf{u}=[u,v,w]$. These waves are represented on Fig. 2 on the inlet ($x=0$) and outlet ($x=L$) of a computational domain.

At the outlet ($x=L$ on Fig. 2), $\mathcal{L}_{2}$, $\mathcal{L}_{3}$, $\mathcal{L}_{4}$ and $\mathcal{L}_{5}$ leave the computational domain and are obtained with the general expression of characteristic waves:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Let us notice that $\mathcal{L}_{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 $\mathcal{L}_{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 $\mathcal{L}_{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

[TABLE]

where $K_{1}=\sigma(1-M^{2})c_{s}/L$, $p_{\infty}$ is the target pressure at the outlet, $\sigma$ is a constant, $M$ is the maximum Mach number in the flow and $L$ is a characteristic size of the domain [11].

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the isothermal case, (17) and (18) 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 [28, 17]:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The transverse waves can be computed as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\mathbf{u_{t}}=[v,w]$ and $\bm{\nabla_{t}}=[\partial_{y},\partial_{z}]$. The second transverse wave $\mathcal{T}_{2}$ is not introduced here since it is null in the isothermal case, as well as $\mathcal{L}_{2}$. The non reflective outflow boundary condition needs now to be set as:

[TABLE]

with $K_{1}=\sigma(1-M^{2})c_{s}/L$, $K_{2}$ should be equal to the Mach number of the mean flow and $\mathcal{T}_{1,exact}$ is a desired steady value of $\mathcal{T}_{1}$  [28]. 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 $\mathcal{L}_{1}$.

Local streamline LODI

Another potential solution is to compute the LODI equation in the local streamline based frame $\tilde{\mathcal{R}}$ (Fig. 3) [29]. In order to compute the new characteristic waves $\tilde{\mathcal{L}}_{i}$, the non-dimensional velocity vector is projected into the new reference frame with the difficulty to compute $\tilde{x}$-derivative terms from the lattice discretization. A simple approximation is to set

[TABLE]

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 [20]) and those replacing all particle populations (e.g. Regularized BC [27]). Izquierdo and Fueyo [14] and Jung et al. [17] 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 [15]. Heubes et al. [18] decided to impose the CBC physical values thanks to the equilibrium populations, which is known to impose incorrect gradients at the boundary [21]. 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. 4). After streaming, three incoming populations are missing : $f_{1}$, $f_{2}$ and $f_{3}$. The zeroth and first-order hydrodynamic moments are:

[TABLE]

[TABLE]

where $\rho_{+}=f_{5}+f_{6}+f_{7}$ and $\rho_{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

[TABLE]

where $\rho$ and $u$ are still non-dimensional. Thus, $\rho$ and $u$ are linked with relation (34), which proves that it is impossible to impose both $\rho_{b}$ and $\mathbf{u_{b}}$ 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 [20], the missing populations can be computed as follows:

[TABLE]

[TABLE]

[TABLE]

while every other populations are kept unchanged:

[TABLE]

with $f_{i}^{(neq)}=f_{i}-f_{i}^{(eq)}$ and where the equilibrium populations are computed with the physical values imposed at the boundary, $\rho_{b}$ and $\mathbf{u_{b}}=(u_{b},v_{b})$. As shown in [20], corrections (35), (36) and (37) then allow to impose the first-order moment at the boundary $\rho_{b}\mathbf{u_{b}}$. However, density and velocity are still linked with (34). This is not a problem for a dirichlet Zou/He boundary condition where only one physical value (either $\rho$ or $\mathbf{u}$) is imposed and the other one is computed with (34). 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 [19]. The solution proposed here is to add a correction on the population associated to a null velocity only:

[TABLE]

which ensures the value of $\rho_{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 (35), (36), (37) and (39) while other populations are kept unchanged after streaming:

[TABLE]

This method can be easily transposed on every 3D lattice (except for high order lattices [30]) by using the general formula of the Zou/He boundary conditions that can be found in [21] for the missing populations, and correction (39) to ensure the value of $\rho_{b}$ and consequently $\mathbf{u_{b}}$.

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 [21]. 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. [21].

More details about the regularized method for boundary conditions can be found for example in [21]. The purpose of this section is to explain how this particular boundary condition is used to impose $\rho_{b}$ and $\mathbf{u}_{b}$ on a flat boundary.

The leading order of the populations $f_{i}$ can be expressed (see end of (5) and (9)) as

[TABLE]

On a boundary node the density $\rho_{b}$ and velocity $\mathbf{u}_{b}$ are imposed. Therefore in order to be able to use this last equation one needs a way to compute $\bm{\Pi^{(1)}}$. This is achieved by using the fact that $\mathbf{Q}_{i}$ is a symmetric tensor with respect to $i$ which means that $\mathbf{Q}_{i}=\mathbf{Q}_{opp(i)}$, where

[TABLE]

At the leading order, this property leads to

[TABLE]

The known $f_{i}^{(1)}$ can be straightforwardly computed by the following formula

[TABLE]

With the last two equations, the set of $f_{i}^{(1)}$ is complete (they are all known) and can be used to compute $\bm{\Pi^{(1)}}$ through (7). Then using (9) one recomputes regularized $f_{i}^{(1)}$ populations and the total populations $f_{i}$ are computed with the relation (41). 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 $\bm{\Pi^{(1)}}$ by a second order finite difference scheme thanks to (8) and recompute $f^{(1)}$ using (9). 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. 2. 2.

Collision step. 3. 3.

Streaming step: some populations are missing at the boundary. 4. 4.

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.

3 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\times 200$ cells, is initiated with a gaussian plane wave as follows (in lattice units)

[TABLE]

[TABLE]

[TABLE]

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 ($\rho$, $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. 5 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 1.

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 [14, 18].

It is noticed that $\rho=1$ is not correctly recovered at the outlet after reflection, as the boundary has been set as perfectly non-reflecting ($\mathcal{L}_{1}=0$). In order to impose the correct boundary condition, a relaxation should be implemented, as in eq (16). However, in that case, the reflection rate will increase as shown in [11].

2D plane wave with incidence

At $t=0$, the computational domain is at rest ($\rho_{0}=1$, $\mathbf{u}_{0}=\mathbf{0}$) except for an oblique line on which the density is set to $\rho_{0}=1.1$ in order to generate a plane wave with an incidence $\alpha$. 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.

Because of a spurious phenomenon appearing at high incidence angles, only angles below $45^{\circ}$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. 6. 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(\alpha)/{c_{s}}$. The condition for this point of the wave to reach the outlet before being distorted is:

[TABLE]

which means than only incidence angles below $45^{\circ}$ can be computed with such a test case. This problem is avoided in [18] 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. 7 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 $\gamma=3/4$ as introduced in [18].

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 [18]: the reflection rate slightly increases and reaches $5\%$ at $40^{\circ}$ of incidence. On the contrary, when the CBC method is extended with the effects of transverse waves as in [28, 17], the reflection rate decreases until $30^{\circ}$ 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\times 600$ cells is initiated with a gaussian density in order to generate a spherical wave, as follows :

[TABLE]

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. 8). 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:

[TABLE]

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. 9. 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^{\circ}$). When adding transverse term without any relaxation as in [28, 17], the observation confirms the prospective behavior of the plane waves simulations: even if the reflection rate is consistently reduced for angles below $40^{\circ}$, 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\times 600$ grid points. A particular attention must be paid to the initialization of the Lamb-Oseen vortex [31] 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:

[TABLE]

[TABLE]

[TABLE]

where $u_{0}=0.1$ in lattice units, $\beta=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 [32]:

[TABLE]

[TABLE]

where $D=2$ is the dimension of the problem. The specific heat capacity at constant volume $C_{v}$ appears instead of $C_{p}$ in (49) because of an error in the heat flux obtained in the Navier-Stokes equations after the Chapman-Enskog development for athermal lattices [32]. 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 [33]. 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 $\mathcal{L}_{1}=\mathcal{T}_{1}$,

(d) CBC-2D relaxed including transverse terms with:

[TABLE]

where $M=0.2$ and $\sigma=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.10 and 11 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 $\mathcal{L}_{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 $\mathcal{L}_{1}$ is not relaxed: the numerical stability comes from the LBM adaptation itself.

Fig. 12 compares the stability of the Zou/He adaptation and the Regularized BB adaptation. As predicted from [21], 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.

4 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 8$C\times 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 $\Delta x=1/400$ m, corresponding to a time step $\Delta t=4.25\times 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\times 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\times 8C$ domain, Fig. 13. 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 $\Delta 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. 13. 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,

• 2.

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}=\sqrt{\overline{p^{\prime 2}}}/{\overline{p}}$, to underline the behavior of each CBC during the whole computation as in Fig. 14.

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. 15 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.

5 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^{\circ}$. 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.