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

1.5 MB·

by  Gauthier Wissocq

Regularized characteristic boundary conditions for the Lattice-Boltzmann
methods at high Reynolds number flows
GauthierWissocqa,b,c,NicolasGourdaina,OrestisMalaspinasd,AlexandreEyssartierb
aISAE,Dpt.ofAerodynamics,EnergeticsandPropulsion,Toulouse,France
bAltran,DOME,Blagnac,France
cCentreEurope´endeRechercheetdeFormationAvance´eenCalculScientifique(CERFACS),CFDTeam,42avenueGaspardCoriolis,31057
ToulouseCedex01,France
dSPC-CentreUniversitaired’Informatique,Universite´deGene`ve7,routedeDrize,CH-1227Switzerland
7
1
0
2

n Abstract
a
J This paper reports the investigations done to adapt the Characteristic Boundary Conditions (CBC) to the Lattice-

6 BoltzmannformalismforhighReynoldsnumberapplications. ThreeCBCformalismsareimplementedandtestedin
2 anopensourceLBMcode: thebaselinelocalone-dimensioninviscid(BL-LODI)approach, itsextensionincluding
theeffectsofthetransverseterms(CBC-2D)andalocalstreamlineapproachinwhichtheproblemisreformulatedin

]
n theincidentwaveframework(LS-LODI).ThenallimplementationsoftheCBCmethodsaretestedforavarietyoftest
y cases,rangingfromcanonicalproblems(suchas2Dplaneandsphericalwavesand2Dvortices)toa2DNACAprofile
d athighReynoldsnumber(Re=105),representativeofaeronauticapplications. TheLS-LODIapproachprovidesthe
-
bestresultsforpureacousticswaves(planeandsphericalwaves). However, itisnotwellsuitedtotheoutflowofa
u
l convected vortex for which the CBC-2D associated with a relaxation on density and transverse waves provides the
f
. bestresults. Asregardsnumericalstability,aregularizedadaptationisnecessarytosimulatehighReynoldsnumber
cs flows. The so-called regularized FD (Finite Difference) adaptation, a modified regularized approach where the off-
i equilibriumpartofthestresstensoriscomputedthankstoafinitedifferencescheme,istheonlytestedadaptationthat
s
y canhandlethehighReynoldscomputation.
h
p Keywords: LatticeBoltzmannmethod,characteristicboundaryconditions,LODI,highReynoldsnumberflows
[

1
Introduction complex geometries (e.g. aircraft and gas turbines).
v
4 However, the resolution of the NS equations requires
Abetterunderstandingofturbulentunsteadyflowsis
3 to add artificial dissipation to ensure numerical stabil-
7 a necessary step towards a breakthrough in the design ity [1]. The consequence is an over-dissipation which
7 of modern aircraft and propulsive systems. Due to the affects the flow and limits the capability to transport
0 difficultyofpredictingturbulencewithcomplexgeom-
flow patterns (like turbulence) on a long distance. In
.
1 etry,theflowthatdevelopsintheseenginesremainsdif- somespecificcases,likeaero-acoustic(far-fieldnoise),
0 ficulttopredict. Atthistime,themostpopularmethod NScanthusfacesomedifficultiestopredicttheflow.
7 to model the effect of turbulence is still the Reynolds
1
Averaged Navier-Stokes (RANS) approach. However In this context, there is an increasing interest in
:
v there is some evidence that this formalism is not ac- the fluid dynamics community for emerging methods,
i
X curate enough, especially when a description of time- basedontheLatticeBoltzmannapproach[2,3,4]. The
dependentturbulentflowsisdesired(highincidencean- Lattice-BoltzmannMethod(LBM)hasalreadydemon-
r
a gle, laminar-to-turbulent transition, etc.) [1]. With the strated its potential for complex geometries, thanks
increase in computing power, Large Eddy Simulation to immersed boundary conditions (that allow the use
(LES) applied to the Navier-Stokes equations emerges of cartesian grids) and low dissipation properties re-
as a promising technique to improve both knowledge quired for capturing the small acoustic pressure fluc-
of complex physics and reliability of flow solver pre- tuations [5, 6]. LBM also provides the advantage of
dictions [1]. It is still the most popular and mature aneasyparallelization, makingitwellsuitedforHigh-
approach to describe the behavior of turbulent flow in PerformanceComputing[7]. However,themostwidely
PreprintsubmittedtoJournalofComputationalPhysics January27,2017
used Lattice-Boltzmann models still suffer from weak- arycondition[20].
nesses like a lack of robustness for high Mach number Still,previousresearchesarelimitedtolowReynolds
flows (M > 0.4), a limitation to low compressible number applications while Latt et al. showed that in-
isothermal flows [2] and the use of artificial boundary creasing the Reynolds number can have drastic impact
conditions(Dirichlet/Neumanntypescanleadtothere- on the numerical stability of the LBM boundary con-
flection of outgoing acoustic waves that have a signif- dition [21]. Even when the MRT collision is used, as
icant influence on the flow field [8, 9, 10]). While the in [14,15], thenumerical stability ofthe characteristic
use of artificial boundary conditions is also critical for boundary conditions has not been demonstrated. The
NSmethods,itismoreproblematicforLBMduetothe aimofthisstudyisthustodevelopanumericallystable
lowdissipationofthemethod. adaptationoftheCBCtotheLBMformalismtakingad-
A potential way to avoid unphysical acoustic reflec- vantageoftheregularizedcollisionscheme[22],which
tions at the boundary is to use a ”sponge layer” inside has proved to be numerically stable at high Reynolds
the computational domain, on which artificial dissipa- number flows [23], and the corresponding regularized
tion(byupwinding)isintroducedorphysicalviscosity boundaryconditions[21]. DifferentkindsofCBCwill
is increased (viscosity sponge zones). Acoustic waves also be evaluated. This article is structured as follows.
(physicalornot)arethusdampedinsuchazone,which ThefirstsectiondescribestheLBMframework. Then,
allows to eliminate or limit numerical reflections [9]. the second section presents three kinds of CBCs and
This solution has however important drawbacks. First three possible adaptations to the LBM formalism for
thecalibrationofspongelayersisdifficult,asabalance 2D problems in the low compressible isothermal case.
mustbefoundbetweenabrutalincreaseoftheviscos- Inthethirdsection, thesemodelsareassessedforsim-
ity(thatwillgenerateacousticreflections)andatoolow ple cases: normal, oblique and spherical waves and a
dissipationthatwillnotbeeffective. Suchspongelay- convected vortex at Re = 103. Finally, the method
ersalsohaveanimpactonthecomputationalcost,since is assessed on a high Reynolds number application: a
a part of the domain is dedicated to slowly increasing NACA0015airfoilatRe=105.
the viscosity. Last, some boundary conditions cannot
betreatedwithaspongelayer,forinstanceaninletwith
1. Numericalmethodandgoverningequations
turbulenceinjection.
For NS methods, a successful approach is the use LatticeBoltzmannframeworkforisothermalflows
ofnon-reflectiveboundaryconditionsbasedonatreat- A description of the Lattice-Boltzmann Method can
ment of the characteristic waves of the local flow befoundin[2,3,4]. Thegoverningequationsdescribe
[11, 12, 13]. However, the extension of this approach theevolutionoftheprobabilitydensityoffindingaset
to LBM is not straightforward given the difficulty to ofparticleswithagivenmicroscopicvelocityatagiven
find a bridge between the LBM that describes the location:
world at the mesoscopic level (a population of parti-
cles) and the NS world based on a macroscopic de- fi(x+ci∆t,t+∆t)= fi(x,t)+Ωi(x,t) (1)
scription of the flow. But some progress has recently
for0 ≤ i < q, wherec isadiscretesetofqvelocities,
been made on the adaptation of characteristic bound- i
f(x,t)isthediscretesingleparticledistributionfunction
ary conditions to the LBM formalism. Izquierdo and i
corresponding to c and Ω is an operator representing
Fueyo [14] used a pressure antibounceback boundary i i
theinternalcollisionsofpairsofparticles.
condition [15] adapted to the multiple relaxation time
Macroscopic values such as density, ρ, and the flow
(MRT) collision scheme [16] to impose the Dirichlet
velocity, u, can be deduced from the set of probability
densityandvelocityconditionsgivenbythelocalone-
densityfunctions f(x,t),suchas:
dimensionalinviscid(LODI)equations,whichprovided i
non-reflective outflow boundary conditions for one- (cid:88)q−1 (cid:88)q−1
dimensionalwaves. Morerecently,Jungetal.[17]ex- ρ= fi, ρu= fici. (2)
tendedthepreviousworktoincludetheeffectsoftrans- i=0 i=0
verseandviscoustermsintheCharacteristicBoundary Someofthemostpopularchoicesforthesetofveloc-
Conditions (CBC) and showed good performance for itiesareD2Q9andD3Q27lattices,respectively9veloc-
vortexoutflow. Meanwhile, Heubesetal.[18]adapted itiesin2Dand27velocitiesin3D(seeFig.1). Forboth
the solution given by a modified-Thompson approach of these lattices, the sound speed in lattice units (nor-
byimposingthecorrespondingequilibriumpopulations malizedbytheratiobetweenthespatial√resolutionand
and Schlaffer [19] assessed a modified Zou/He bound- thetimestep∆x/∆t)isgivenbyc =1/ 3[3].
s
2
c12 c22 c24 withthestrainratetensorS = (∇u+(∇u)T)/2through
c5 c15 therelation
c1 c8 c7 c13 c6 c21 c23 c17 Π(1) =−2c2sρτS. (8)
c2 c0 c6 c7 c1 c16 cc319 c14 c20 Inturntotheleadingorder fi(1)canbeapproximatedby
y c3x c4 c5 c4 c10 c2 c8 c18c26 y x fi(1) (cid:27) 2wci2sQi :Π(1), (9)
z (cid:16) (cid:17)
c11 c9 c25 where Q ≡ cc −c2I . The colon symbol stands for
i i i s
thedoublecontractionoperatorandIistheidentityma-
Figure1: SchematicplotofvelocitydirectionsoftheD2Q9model trix.
(left)andtheD3Q27model(right)
Aregularizationstepconsistinginthereconstructionof
theoff-equilibriumpartsusing(9)and(7)canimprove
the precision and the numerical stability of the single
relaxationtimeBGKcollision[22,23]. Thismodelwill
The collision operator Ω is usually modelled with beusedinthenextpartsforhighReynoldsnumberflow
i
theBhatnagar-Gross-Krook(BGK)approximation[24], simulations.
which consists in a relaxation, with a relaxation time
τ,ofeverypopulationtothecorrespondingequilibrium ThePalabosopen-sourcelibrary
probabilitydensityfunction f(eq): The LBM flow solver used in this work is the Pal-
i
abos1 open-source library. The Palabos library is
1(cid:104) (cid:105)
Ω =− f(x,t)− f(eq)(x,t) . (3) a framework for general-purpose CFD with a kernel
i τ i i
based on the lattice Boltzmann method. The use of
The equilibrium distribution function f(eq) is a local C++ code makes it easy to install and to run on ev-
i ery machine. It is thus possible to set up fluid flow
functionthatonlydependsondensityandvelocityinthe
simulations with relative ease and to extend the open-
isothermalcase. Itcanbecomputedthankstoasecond
sourcelibrarywithnew methodsandmodels, whichis
order development of the Maxwell-Boltzmann equilib-
ofparamountimportancefortheimplementationofnew
riumfunction[25]:
boundaryconditions. Thenumericalschemeisdivided
fi(eq) =wiρ1+ cic·2u +(cid:32)c2ic·2u(cid:33)2− 2uc22, (4) intwosteps:
s s s • AcollisionstepwheretheBGKmodelisapplied:
wherew arethegaussianweightsofthelattice.
i 1 1(cid:104) (cid:105)
AChapman-Enskogexpansion,basedontheassump- fi(x,t+2)= fi(x,t)+τ fi(eq)(x,t)− fi(x,t) , (10)
tionthat f isgivenbythesumoftheequilibriumdistri-
i
butionplusasmallperturbation fi(1) with f(eq) computedusingthemacroscopicvalues
i
at time t and f can be regularized in order to in-
f = f(eq)+ f(1), with f(1) (cid:28) f(eq), (5) i
i i i i i creasenumericalstabilityforhighReynoldsnum-
berflows.
can be applied to (1) in order to recover the exact
Navier-Stokes equation for quasi-incompressible flows
• Astreamingstep:
in the limit of long-wavelength [26]. The pressure is
thus given by p = c2ρ and the kinematic viscosity is 1
s f(x+c,t+1)= f(x,t+ ). (11)
linkedtotheBGKrelaxationparameterthrough i i i 2
ν=c2(τ− 1). (6) The streaming step consists in an advection of each
s 2 discrete population to the neighbor node located in the
directionofthecorrespondingdiscretevelocity. Sincea
The Chapman–Enskog expansion also relates the sec-
ondordertensorΠ(1)definedas boundarynodehaslessneighborsthananinternalnode
(cid:88)
Π(1) = cc f(1), (7)
i i i 1Copyright2011-2012FlowKitLtd.
i
3
(less than 9 neighbors in 2D or 27 neighbors in 3D),
y
somepopulationsaremissingattheboundaryaftereach
iteration. These populations need to be reconstructed,
which is the purpose of the implementation of bound-
ary conditions in LBM. Up to now, different methods
canbeusedinPalabos, suchasregularizedBC[27]or (u+c ) (u+c )
5 s 5 s
Zou/HeBC[20]toimplementopenboundaries. How- L L
ever,noneofthemcanbeusedastheystandforanout- 2(u) 2(u)
L L
flow boundary condition and the use of sponge zones (u) (u)
3 3
L L
isnecessarytoavoidnonphysicalreflections. Thenext (u) (u)
4 4
sectionswillaimatdevelopingamorenaturalboundary L L
conditionthatminimizeacousticreflectionsforanout- (u c ) (u c )
1 s 1 s
L � L �
flowboundarytype,basedontheCharacteristicBound-
0 L x
aryConditions(CBC).
2. Adaptation of Characteristic Boundary Condi- z
tionstotheLBMformalism
One of the most popular methods in the NS com- Figure 2: Waves leaving and entering the computational domain
throughaninletplane(x=0)andanoutletplane(x= L)forasub-
munity for subsonic non reflective outflow boundary
sonicflow[11]inlatticeunits.
conditions is the CBC method [11]. The adaptation of
the CBC to the LBM formalism is presented here for
∂v
an isothermal flow, in lattice units (normalized by ∆x L =u , (13)
3 ∂x
and∆t). Acousticwavesthuspropagateattheconstant
latticesoundspeedcs. Intheisothermalcase,pressure L =u∂w, (14)
is defined as p = c2ρ. Three CBC methods will be 4 ∂x
s
introduced below: the local one-dimensional inviscid (cid:32) (cid:33)
∂p ∂u
(LODI) approximation, a 2D extension of the LODI L5 =(u+c) ∂x +ρcs∂x . (15)
approximation including transverse waves and a last
methodcalledlocal-streamlineLODI. Let us notice that L , the entropy wave, is null in the
2
isothermalcase.Thex-derivativetermscanbecomputed
using the interior points by one-sided finite difference.
The treatment of L is different : since it comes from
LODIApproximation(BaselineLODI) 1
the outside, it can not be computed using the interior
Let us consider a domain outlet located at x = L,
points. Theperfectlynon-reflectingcaseisobtainedby
as descripted on Fig. 2. A diagonalisation of the x- fixingL = 0,whichensureseliminatingtheincoming
1
derivative terms in the Navier-Stokes equation allows
wave. However, this is known to be unstable because
to define five waves L that propagate respectively at
i oflackofcontroloftheoutletflowvariables. Asimple
velocity u − c , u, u, u and u + c , where u is the x-
s s waytoensurewell-posednessistoset
component(streamwise)ofthenon-dimensionalmacro-
scopic velocity u = [u,v,w]. These waves are repre- L1 = K1(p−p∞), (16)
sentedonFig.2ontheinlet(x = 0)andoutlet(x = L) where K = σ(1− M2)c /L, p is the target pressure
1 s ∞
ofacomputationaldomain. attheoutlet,σisaconstant, M isthemaximumMach
numberintheflowand L isacharacteristicsizeofthe
domain[11].
At the outlet (x = L on Fig. 2), L , L , L and L
2 3 4 5 Thetime-derivativeoftheprimitivevariablescanbe
leave the computational domain and are obtained with
computedinfunctionofthewaveamplitudesbyexam-
thegeneralexpressionofcharacteristicwaves:
iningaLODIproblem:
(cid:32) (cid:33) (cid:34) (cid:35)
∂ρ ∂p ∂ρ 1 1
L =u c2 − =0, (12) + L + (L +L ) =0, (17)
2 s∂x ∂x ∂t c2 2 2 5 1
s
4
∂p + 1(L +L )=0, (18) the isothermal case, as well as L2. The non reflective
∂t 2 5 1 outflowboundaryconditionneedsnowtobesetas:
∂∂ut + 2ρ1c (L5−L1)=0, (19) L1 = K1(p−p∞)−K2(T1−T1,exact)+T1, (30)
s
with K = σ(1− M2)c /L, K should be equal to the
1 s 2
∂v +L =0, (20) MachnumberofthemeanflowandT1,exact isadesired
∂t 3 steady value of T [28]. In the rest of the paper, this
1
method will be named 2D-CBC in the perfectly non-
∂w
+L =0. (21) reflecting case (K = K = 0) and 2D-CBC relaxed
∂t 4 1 2
whenarelaxationisdoneonL .
1
Intheisothermalcase,(17)and(18)areequivalent. Fi-
nally, with a temporal discretization using an explicit LocalstreamlineLODI
second-order scheme, the physical values that must be
imposedatthenexttimestepinordertoavoidacoustic (a)geometry-basedframeR
reflectionscanbecomputed. Thisimplementationwill
be referred to as the baseline LODI (BL-LODI) in the Computational
Outside
restofthepaper. domain
LODIapproximationincludingtransverseterms
Thepreviousrelationsareperfectlynon-reflectingfor v u
theonlycaseofapure1Dplanewave.Foranon-normal
(u c )
wave,theLODIapproximationisnotverifiedandare- L1 � s
u (u+c )
flected wave, all the more important as the incidence L5 s
increases, can appear. To take this phenomenon into
account, a possible solution is to add the influence of
y
transversewavesintheLODIequations[28,17]:
x
∂ρ 1 1 R
∂t + 2c2(L5+L1)= 2c2(T5+T1), (22) (b)localstreamline-basedframeR˜
s s
Computational
∂u 1 1 Outside
+ (L −L )= (T −T ), (23) domain
∂t 2ρc 5 1 2ρc 5 1
s s
∂v
∂t +L3 =T3, (24) v L5
u
e
∂w
+L =T . (25)
∂t 4 4 u
L1
Thetransversewavescanbecomputedasfollows:
e
T =−(cid:2)u ·∇ p+p∇ ·u −ρc u ·∇ u(cid:3), (26) y˜ x˜
1 t t t t s t t
R
(cid:34) (cid:35)
1∂p
T =− u ·∇ v+ , (27) e
3 t t ρ ∂y
Figure3:Schematicofcharacteristicwaveprojectionsintwocoordi-
(cid:34) (cid:35)
1∂p natesystems: a)geometrybasedframeR. b)localstreamlinebased
T4 =− ut·∇tw+ ρ ∂z , (28) frameR˜.
T =−(cid:2)u ·∇ p+p∇ ·u +ρc u ·∇ u(cid:3), (29) Another potential solution is to compute the LODI
5 t t t t s t t
equation in the local streamline based frame R˜
whereu = [v,w]and∇ = [∂ ,∂ ]. Thesecondtrans- (Fig. 3) [29]. In order to compute the new character-
t t y z
verse wave T is not introduced here since it is null in istic waves L˜ , the non-dimensional velocity vector is
2 i
5
projected into the new reference frame with the diffi- incoming populations are missing : f , f and f . The
1 2 3
cultytocomputex˜-derivativetermsfromthelatticedis- zerothandfirst-orderhydrodynamicmomentsare:
cretization. Asimpleapproximationistoset
ρ= f1+ f2+ f3+ρ0+ρ+, (32)
∂φ˜ ∂φ˜
∂x˜ = ∂x, (31) ρu=ρ+−(f1+ f2+ f3), (33)
and thus to compute it by a first-order upwind scheme where ρ+ = f5 + f6 + f7 and ρ0 = f0 + f4 + f8. These
using the lattice discretization. This implementation equationscanbecombined,byeliminatingthemissing
of the CBC condition will be referred to as the local populations(f + f + f ),toobtain
1 2 3
streamlineLODI(LS-LODI).
1
ρ= 1−u(ρ0+2ρ+), (34)
AdaptationtotheLatticeBoltzmannscheme
where ρ and u are still non-dimensional. Thus, ρ and
The main difficulty in LBM is then to find a set of
u are linked with relation (34), which proves that it is
populations in order to impose the physical values ob-
impossible to impose both ρ and u at the boundary
b b
tained by the CBC theory and the correct associated
without modifying any of the known populations.The
gradients. The possible adaptations can be divided in
same relation can be obtained for every 1D, 2D or 3D
twofamilies: thosepreservingtheknownparticlepop-
latticeasfarasthereisonlyonelevelofvelocity. This
ulations (e.g. Zou/He BC [20]) and those replacing
relation concerns every boundary condition of the first
all particle populations (e.g. Regularized BC [27]).
family (those preserving the knwon populations) and
Izquierdo and Fueyo [14] and Jung et al. [17] chose
willbeusedlater.
to modify the missing populations by adapting a pres-
sureanti-bouncebackboundarycondition,whichcanbe
Letusnoteg thecorrectedpopulationsatthebound-
i
usedasfarastheMRTcollisionschemeisadopted[15].
ary. As proposed in [20], the missing populations can
Heubes et al. [18] decided to impose the CBC physi-
becomputedasfollows:
calvaluesthankstotheequilibriumpopulations,which
is known to impose incorrect gradients at the bound- 1
g = f(eq)+ f(neq)+ (f(neq)− f(neq)), (35)
ary [21]. Three possibilities are introduced below and 1 1 5 2 4 8
will be further evaluated: a modified Zou/He method
and two declinations of the more stable regularized g2 = f2(eq)+ f6(neq), (36)
adaptation.
1
g = f(eq)+ f(neq)− (f(neq)− f(neq)), (37)
3 3 7 2 4 8
Computational
Outside
domain whileeveryotherpopulationsarekeptunchanged:
f f g = f, i=0,4,5,6,7,8, (38)
f 8 7 i i
1
f0 with fi(neq) = fi − fi(eq) and where the equilibrium
f f populations are computed with the physical values im-
2 6
y posedattheboundary, ρ andu = (u ,v ). Asshown
b b b b
in[20],corrections(35),(36)and(37)thenallowtoim-
x f3 f4 f5 posethefirst-ordermomentattheboundaryρbub.How-
ever,densityandvelocityarestilllinkedwith(34).This
isnotaproblemforadirichletZou/Heboundarycondi-
Figure4: Schematicplotofthepopulationsetattheboundaryofa tionwhereonlyonephysicalvalue(eitherρoru)isim-
D2Q9model.Afterthestreamphase, f1, f2and f3aremissing. posedandtheotheroneiscomputedwith(34). Onthe
contrary,foranonreflectiveoutflow,bothofthemneed
to be imposed as the result given by the CBC method,
so that the only condition set by the user is the value
AdaptationwithZou/Heboundaryconditions ofacharacteristicwave. Itisthennecessarytomodify
Let us consider an outflow boundary located at x = at least one known population. Schlaffer suggests cor-
L on a D2Q9 lattice (Fig. 4). After streaming, three rectingeverypopulationsinordertoimposethecorrect
6
density[19].Thesolutionproposedhereistoaddacor- Attheleadingorder,thispropertyleadsto
rection on the population associated to a null velocity
only: f(1) = f(1) . (43)
i opp(i)
1
g0 = f0+ρb− 1−u(ρ0+2ρ+), (39) The known fi(1) can be straightforwardly computed by
thefollowingformula
whichensuresthevalueofρ . Thischoiceismotivated
b
by the fact that this added correction will only affect f(1) = f − f(eq)(ρ ,u ). (44)
i i i b b
the collision phase and will not be streamed into the
computationaldomain. With the last two equations, the set of f(1) is complete
i
(they are all known) and can be used to compute Π(1)
To sum up this method, g , g , g and g are com- through(7). Thenusing(9)onerecomputesregularized
1 2 3 0
puted with respectively (35), (36), (37) and (39) while f(1) populations and the total populations f are com-
i i
otherpopulationsarekeptunchangedafterstreaming: putedwiththerelation(41).Thismethodisvalidinboth
2Dand3DandwillbecalledRegularizedBounceback
g = f, i=4,5,6,7,8. (40)
i i (orRegularizedBB)adaptationinthenextsections.
This method can be easily transposed on every 3D Another possibility is to compute Π(1) by a second or-
lattice(exceptforhighorderlattices[30])byusingthe derfinitedifferenceschemethanksto(8)andrecompute
generalformulaoftheZou/Heboundaryconditionsthat f(1) using(9). ThismethodwillbecalledtheRegular-
can be found in [21] for the missing populations, and izedFDadaptation.
correction (39) to ensure the value of ρ and conse-
b
quentlyub. Summaryofthemethod
Adaptationwiththeregularizedmethod Thenonreflectingoutflowboundaryconditionusing
The Zou/He boundary condition provides the ad- CBC theory for LBM can be summarized as follows,
considering everything is known at (non-dimensional)
vantage of a very good precision in the definition
timet:
of the boundary physical values. Unfortunately, this
solutionsuffersfromlackofstabilityforlargeReynolds
1. Computation of the physical values that must be
numbers. For example, it has been shown that Zou/He imposed at t+1 to avoid non physical reflections
boundary conditions become unstable at Re > 100 for
by the CBC theory using either BL-LODI, CBC-
a given resolution N = 200 nodes per characteristic
2D or LS-LODI method. These values are stored
length in a 2D channel flow [21]. Another possible
tobeusedinthelaststep.
adaptationistousetheregularizedboundarycondition
inordertoimposethephysicalvaluescomputedbythe
2. Collisionstep.
CBCtheory. Thissolutionisyetlessaccuratebutwell
morestable,asshownbyLattetal. [21].
3. Streaming step: some populations are missing at
theboundary.
Moredetailsabouttheregularizedmethodforbound-
ary conditions can be found for example in [21]. The
purposeofthissectionistoexplainhowthisparticular 4. Correction of the set of populations at the bound-
boundaryconditionisusedtoimposeρ andu onaflat ary so that the physical values stored in the first
b b
boundary. step are imposed by using the Zou/He adaptation
The leading order of the populations f can be ex- the so-called regularized Bounceback adaptation
i
pressed(seeendof(5)and(9))as ortheRegularizedFDadapatation.
f = f(eq)(ρ ,u )+ f(1)(Π(1)). (41)
i i b b i
Onaboundarynodethedensityρ andvelocityu are
b b 3. Applicationtoacademiccases
imposed. Therefore in order to be able to use this last
equation one needs a way to compute Π(1). This is
In this section, the CBC approach for LBM is as-
achievedbyusingthefactthatQ isasymmetrictensor
i sessed on simple 2D cases: a normal plane wave, a
withrespecttoiwhichmeansthatQ =Q ,where
i opp(i) planewave withdifferentincidence angles, a spherical
opp(i)={j|c =−c }. (42) densitywaveandaconvectedvortex.
i j
7
2Dnormalplanewave
The computational domain, a square of 200 × 200
cells,isinitiatedwithagaussianplanewaveasfollows
(inlatticeunits)
(cid:32) (cid:33)
(x−x )2
ρ =1+0.1∗exp − 0 ,
0 2σ2
u =0.1,
0
(cid:32) (cid:33)
(x−x )2
v =0.1∗exp − 0 ,
0 2R2
c
withx =110and2R2 =20inlatticeunits(i.e. innum-
0 c
berofcells). TheReynoldsnumber,computedwiththe
horizontalnon-dimensionalinitialvelocityu ,thechar-
0
acteristic size of the box in number of voxels and the
viscosityinlatticeunits,isequalto100. Theboundary
conditions are: vertical periodicity, reflecting inlet on
theleftandperfectlynonreflectingoutflowontheright.
Figure5:Densitywavesbeforeandafterreflectionofthe(u+cs)wave
Forthispure1Dcase,thechoiceoftheCBCcondition ontherightboundary(non-reflectingoutflow).
(baseline, local streamline LODI or LODI with trans-
verseterms)hasnoeffectontheabsorptionrate. More-
over,allthethreeadaptationsprovidethesameresultsat As often with CBC, the treatment is more difficult
10−6andnostabilityissueshavebeenencounteredwith forthepressurewave, buttheresultsobtainedhereare
the Zou/He adaptation for this test case. Thus, the re- in good agreement compared to what is found in the
sultspresentedbelowhavebeenobtainedwithbaseline literature[14,18].
LODI and Zou/He adaptation. This test case has been It is noticed that ρ = 1 is not correctly recovered at
chosensothatthereflectionrateforeverymacroscopic theoutletafterreflection, astheboundaryhasbeenset
value(ρ,uandv)canbecomputed,sincetwopressure asperfectlynon-reflecting(L1 =0). Inordertoimpose
andaxialvelocitywavespropagatesatspeed(u−c )and the correct boundary condition, a relaxation should be
s
(u+c )andonetransversevelocitywavepropagatesat implemented, as in eq (16). However, in that case, the
s
speed u. The computed density waves are represented reflectionratewillincreaseasshownin[11].
onFig.5attwodifferenttimesteps:beforereflectionof
the(u+cs)waveonthenon-reflectiveoutletandshortly 2Dplanewavewithincidence
afterreflection.
Att=0,thecomputationaldomainisatrest(ρ =1,
Averylowreflectedamplitudecanbedistinguished. 0
The (u−c ) wave propagates to the left of the domain u0 = 0) except for an oblique line on which the den-
s sity is set to ρ = 1.1 in order to generate a plane
withoutencounteringanyboundaryatthetwoobserved 0
timesteps. Itisonlyaffectedbyviscousdissipationand wave with an incidence α. The reflection coefficient
is measured by computing the ratio of maximal ampli-
is used as a reference amplitude in the computation of
tudes in density waves only, as this is the most criti-
thereflectionratesofdensityandaxialvelocity. Forthe
calhydrodynamicvariable. Thethreeimplementations
transverse velocity wave, one can compute the reflec-
aretestedforthiscase: BL-LODI,LS-LODIandCBC-
tionrateastheratiobetweenthereflectedwaveampli-
2D in the perfectly non-reflecting case. Again, the re-
tude and the amplitude shortly before reflection. The
sults obtained with the Zou/He adaptation only will be
obtainedresultsarepresentedonTable1.
presented here, as the same results at 10−6 have been
obtainedwiththeRegularizedBBandRegularizedDF
Variable ρ u v
Reflectionrate 1.2% 1.1% <10−5% adaptations.
Becauseofaspuriousphenomenonappearingathigh
Table1:Reflectionratesforthe2Dnormalplanewave. incidence angles, only angles below 45◦could be com-
puted.Indeed,contrarytothepreviousnormalwavetest
8
Periodicity Figure6:
Schematic
plot of the
α N initialization
c o of the plane
s n
M -re wave test case
fl with an inci-
einlet ective dαe.nceSphearnigclael
eflectiv h outlet wsatpaapnveteaasrneoautsltiyhne-
R
two extremi-
ties to distort
the plane
Periodicity
wave.
case, the incident wave cannot be infinite in its trans- Figure7: Reflectioncoefficient(%)withrespecttotheincidenceof
verse direction. Then, spherical waves appear at each the plane wave: (a) baseline LODI, (b) local streamline LODI, (c)
CBC-2D[17]and(d)Mod-Thompson3/4[18].
extremity of the oblique line initializing the wave, as
shownonFig.6. Letusimagineapoint M locatedata
distancehofoneextremityandmovingwiththeoblique LODIimplementation,thecoefficientremainsstableat
planewaveatvelocitycs.ThesphericalwavereachesM around2%ofreflectedwave.
afteratimeh/c .ThepointMreachesthenon-reflective
s
outlet boundary condition at a time htan(α)/cs. The 2DSphericalwave
condition for this point of the wave to reach the outlet The computational domain, a square of 600 × 600
beforebeingdistortedis:
cellsisinitiatedwithagaussiandensityinordertogen-
h h erateasphericalwave,asfollows:
tan(α)< ⇔tan(α)<1, (45)
cs cs (cid:32) ((x−x )2+(y−y )2)(cid:33)
which means than only incidence angles below 45◦ ρ=1+0.1∗exp − 0 0 ,
2R2
c
can be computed with such a test case. This problem
is avoided in [18] by imposing an ’exact’ solution with R = 3.2, x = 520 and y = 300 nodes. The
c 0 0
obtained on a larger computational domain until the boundary conditions are periodic in the vertical direc-
desiredwavereachestheboundary. Theexactsolution tion, a CBC condition is set on the right boundary and
is then switched to the CBC condition. It has not a reflective boundary condition on the left (located far
been tested in this paper in order to avoid the eventual enough to avoid its influence on the spherical wave at
acoustic waves generated by a brutal change in the the studied time steps). Four cases are computed: (a)
boundarycondition. baselineLODI,(b)CBCwithtransverseterms,(c)local
streamline LODI and (d) reference case where the do-
The reflection coefficient of the density wave with main is enlarged in the horizontal direction so that the
respect to the incidence angle is represented on Fig. 7 sphericalwaveisnotaffectedbytheboundary(Fig.8).
for BL-LODI, LS-LODI and CBC-2D. The results are As for the previous test cases, the LBM adaptation of
compared with what is obtained for the same test case theCBCconditionhadnoimpactontheresults.
withthemodifiedThompsonmethodwithacoefficient The local reflection coefficient for such a spherical
γ=3/4asintroducedin[18]. wavecanbecomputedbythefollowingformula:
As expected, for the baseline approach (a), the re-
R−A
flection coefficient increases as the incidence angle in- r= ref, (46)
|I|
creases. For the modified Thompson approach, the re-
sults are close to what can be found in [18]: the re- where I is the amplitude of the original density wave
flection rate slightly increases and reaches 5% at 40◦ running towards the boundary condition, R is the am-
of incidence. On the contrary, when the CBC method plitudeofthedensitywaveafterreflectionattheoutlet
is extended with the effects of transverse waves as andA istheamplitudeofthedensitywaveatthesame
ref
in [28, 17], the reflection rate decreases until 30◦ and lattice node compared to the reference case (d) (no re-
thenbeginstoincrease. Inthecaseofalocalstreamline flection).
9
1200 A map of reflection rates after reflection at the out-
let (after 400 iterations) is shown on Fig. 9. For the
baseline approach, one can see that the reflected wave
300
becomes more important as the local incidence of the
spherical waves increases, because the incident wave
600
y canbeapproximatedasalocalnormalplaneinthecen-
teroftheoutletwhile,inthecorners,theincidencebe-
80
x comes important (up to 75◦). When adding transverse
term without any relaxation as in [28, 17], the obser-
600 vation confirms the prospective behavior of the plane
wavessimulations: evenifthereflectionrateisconsis-
Figure8: Schematicplotofthecomputationaldomain-solidline:
tentlyreducedforanglesbelow40◦,itreaches20%for
CBCcases(a),(b)and(c),dottedline:referencecase(d).Dimensions
aregiveninnumberofcells. the greatest observed angles. On the contrary, the re-
flectionrateremainsconstantwiththelocalstreamline
(a)
LODIimplementation.
2Dconvectedvortex
A2Dvortexisconvectedfromlefttorightandexits
the computational domain, a square of 600×600 grid
points. Aparticularattentionmustbepaidtotheinitial-
izationoftheLamb-Oseenvortex[31]whichhastobe
adapted to the isothermal case in order to avoid spuri-
ous waves generated by the adaptation of a wrong ini-
tial density. The initial conditions, in lattice units, are
(b) imposedasfollows:
(cid:32) (cid:33)
(y−y ) r2
u=u −βu 0 exp − , (47)
0 0 R 2R
c c
(cid:32) (cid:33)
(x−x ) r2
v=βu 0 exp − , (48)
0 R 2R
c c
(cid:34) (βu )2 (cid:32) r2(cid:33)(cid:35)1/(γ−1)
ρ= 1− 0 exp − , (49)
2C 2
v
(c) where u = 0.1 in lattice units, β = 0.5, x = y =
0 0 0
300 nodes (the vortex is initially centered on the box),
R = 20nodesandr2 = (x−x )2+(y−y )2. Withthe
c 0 0
BGK collision operator with a single relaxation time,
thesimulatedgashasthefollowingconstants[32]:
D+2
γ= , (50)
D
D
C = c2, (51)
v 2 s
where D = 2 is the dimension of the problem. The
specific heat capacity at constant volume C appears
Figure9:Densitywaveafterreflectionofasphericalwaveattheoutlet v
(400iterations):(a)baselineLODI,(b)CBC-2D[28,17]and(c)local instead of Cp in (49) because of an error in the
streamlineCBC. heat flux obtained in the Navier-Stokes equations af-
tertheChapman-Enskogdevelopmentforathermallat-
tices [32]. The Reynolds number based on u and the
0
10