Studies of accurate multi-component lattice Boltzmann models on benchmark cases required for engineering applications

0.28 MB·

by  Hiroshi Otomo

Studies of accurate multi-component lattice Boltzmann models on
benchmark cases required for engineering applications
Hiroshi Otomo, Hongli Fan, Yong Li, Marco Dressler, Ilya Staroselsky, Raoyang Zhang, Hudong Chen
Exa Corporation, 55 Network Drive,Burlington, Massachusetts 01803, USA
6
1
0
2

n
a Abstract
J
WepresentrecentdevelopmentsinlatticeBoltzmannmodelingformulti-componentflows,implemented
6
on the platform of a general purpose, arbitrary geometry solver PowerFLOW. Presented benchmark

] cases demonstrate the method’s accuracy and robustness necessary for handling real world engineering
n
applicationsatpracticalresolutionandcomputationalcost. Thekeyrequirementsforsuchapproachare
y
d that the relevant physical properties and flow characteristics do not strongly depend on numerics. In
- particular, the strength of surface tension obtained using our new approach is independent of viscosity
u
and resolution, while the spurious currents are significantly suppressed. Using a much improved surface
l
f wetting model, undesirable numerical artifacts including thin film and artificial droplet movement on
.
s inclined wall are significantly reduced.
c
i
s Keywords: Multi-component flow, Lattice Boltzmann
y 2010 MSC: 00-01,99-00
h
p
[
1. Introduction tionality are specified, and an improved scheme is

1 proposed and tested via the simulation of a two-
v Recently,therehasbeenincreasedinterestinen-
dimensional droplet. In Sec. 4, cases with wall
7 gineeringapplicationsofmulti-componentflowsim-
boundaries are discussed and typical issues asso-
3
ulation with the lattice Boltzmann (LB) method,
2 ciated with boundary models are pointed out. A
because of its advantages for complex geometry
1 new boundary model is tested on some benchmark
0 and turnaround time efficiency. The lattice Boltz- cases. We discuss results in Sec. 5. In this paper,
. mannmethod(LBM)isbasedonthekinetictheory
1 all physical quantities are written in lattice units,
whichallowstoconstructphysicalmodelsfrommi-
0 and the discrete lattice time and space increments
6 croscopic as well as from macroscopic viewpoints. are ∆x=∆t=1.
1 While multi-component LBM’s have shown
v: promising results on a large number of academic
2. Multi-component lattice Bolzmann
i cases, numerical accuracy and stability still repre-
X method
sent challenges under extreme conditions such as
r
coarse resolution and low viscosity. Actually, it is
a The commonly used LB equation for multi-
these conditions that are likely to be encountered
component flow can be written as:
in many engineering applications. In this paper,
issues with currently existing schemes and models fα(x+c ∆t,t+∆t)−fα(x,t)=C α+Fα, (1)
i i i i i
are pointed out, and an improved LB scheme is
tested. where α stands for different components (species),
The paper is organizedas follows. In Sec. 2, the c is the discrete velocity and Fα is inter-
i i
standard LBM for multi-component flow is briefly componentinteractionforce[1]. TheD3Q19lattice
reviewed. InSec. 3,issuesrelatedtothebasicfunc- model[2]is usedsothat irangesfrom1to 19,and
Preprint submitted to Computational Science January 7, 2016
C α is the particle collision operator. The simplest force representation at the first order in resolu-
i
andcommonly used one is the BGK collisionoper- tion/time step, but different at the second and
ator [2–5] with a single relaxation time τ for the higher orders. The high order difference does have
α
α-species: significant influence on simulation quality. In this
work we use the forcing term described in [7].
1
C α =− (fα−feq,α). (2) Theresultingfluidvelocityu isthevelocityav-
i τ i i F
α eraged over pre- and post- collision steps,
The equilibrium state feq,α with the third order
i
expansion is defined as: uF =u+g∆t/2 (8)
2 g = gαρα/ρ. (9)
c ·u (c ·u)
feq,α(ρα,u)=ραw [1+ i + i Xα
i i T0 2T02
u2 (c ·u)3 c ·u where gα is the acceleration of the component α
−2T0 + i6T03 − 2iT02 u2]. (3) derivFeαd,βfr/oρmα.tThehiisntqeuracnotmitpyounenitsfhoernceceFfoαr,tβh:cgaαlle=d
β F
Here T0 =1/3 is the lattice temperature, wi is the sPimply velocity.
isotropic weight in D3Q19, ρα is the density of the
component α, and u is the mixture flow velocity:
3. High accuracy bulk solver
ρα = fα (4)
i
Xi
Engineering applications usually require simula-
ρ= ρα = fiα (5) tionsinvolvingvariousmaterialpropertiesandflow
Xα Xα Xi scenarios. Due to the jump of physical character-
u= α ici·fiα. (6) istics at the interface between components, accu-
P Pρ rate representation and simulation of these inter-
faces represents a significant difficulty. There is a
There exist several models that introduce local in-
consensusthatnumericalstabilityandaccuracyre-
teractionsbetweencomponentsthatareresponsible
maintwomajorchallengesindevelopmentofmulti-
for separationbetween the components [1, 6]. One
phase/multi-componentLBflowsolvers. Toensure
of the most commonly used ones is the Shan-Chen
numerical stability, the viscosities cannot be too
potential force:
small,andalsotheviscosityratiobetweendifferent
Fα,β(x)=Gα,βρα(x) w c ρβ(x+c ∆t). components cannot be too large. Numerical arti-
i i i
facts including spurious current could often con-
Xi
(7) taminate flow physics near the interface region. It
Here,thematrixGα,β definesparameterswhichde- becomesevenmorechallengingwhenthesolid-fluid
termine the strength of interaction between com- interaction, i.e. surface wetting, is also considered.
ponents. If Gα,α = 0, the interaction forces only A new LB algorithm for the multi-component flow
exist between different components. The equa- used in this work improves these numerical issues.
tion of state for each component is that of ideal Evenwhenthe interfaceisstatic,numericalarti-
gas. If Gα,α is nonzero, in addition to the inter- facts could provide a source of artificial velocity,
action forces between different components, there which is called spurious velocity (cf. [8]). The
is also a repulsive force within the α-component. propertreatmentofthesephenomenaisrecognized
As a result, the α-component acquires the equa- as one of the key requirements for accuracy and
tion of state of a non-ideal gas and phase transi- stability of the multi-component flow modeling. In
tion within that component becomes possible. In previous studies [9, 10], it is pointed out that the
this paper, phase transitions of single components spurious current is associated with the insufficient
are neglected and Gα,α =0. isotropy of the of the numerical system caused by
There areseveralwaysto apply the forcingterm the discretization.
Fα. The existing approaches have the same body Instead of BGK, we use here a regularized filter
i
2
collision operator [11]: the viscosity and droplet size. Achieving such in-
dependence is an important first step towardssim-
fα(x+c ∆t,t+∆t)=feq,α(ρ ,u)
i i i α ulating complex practical problems.
1 As mentioned above the spurious current prob-
+ 1− fneq,α+Fα. (10)
(cid:18) τ (cid:19) i i lemisbelievedtobecausedbyinsufficientisotropy
α
of discrete schemes [9, 10]. In Fig. 2, maximum
Here τ is the relaxation time of the fluid compo-
α spurious currents are plotted in terms of τ2 and R.
nent α that is related to the kinematic viscosity of
In the left figure, τ2 is varied while τ1 is fixed cor-
that component ν [12]. fneq,α is the regularized
α i respondingtotheinitialR=48. Itisseenthatthe
non-equilibrium distribution function,
spurious current of the modified scheme is lower
fneq,α =Φα :Πα, (11) than the original one for all cases. Furthermore,
i with the modified version the spurious current de-
where Φ is a regularized filter collision operator pendence upon τ and R is much reduced. As a
based on Hermite polynomials [11, 13, 14] and Πα result,onecanestimatethespuriousvelocityquan-
is the non-equilibrium momentum flux tensor for titativelyevenbeforesimulation,evaluateitseffect
different components. Fα is the interaction body on the main flow, and reduce numerical artifacts.
i
force. The general idea and relevant algorithm de- In Fig. 3, the distributions of the velocity field
tails of the regularization can be found in [11, 13– and ρ2 are presented. Here the initial R is 48 and
16]. Here we would like to emphasize that this the relaxation times are τ1 = 0.525 and τ2 = 3.0.
filter collision operator keeps the nonequilibrium The results demonstrate that the new scheme sig-
information of moments up to the desired order, nificantly reduces the spurious current while pre-
for example the 2nd order for the momentum flux serving the density profile and the interface thick-
and the 3rd order for the energy flux, and removes ness.
other higher order nonequilibrium moments in the
Hermite space. Such a filtering procedure could 4. Wall boundary condition
substantially reduce unphysical noise and numer-
ical artifacts and improve numerical stability and For developed methodology of multi-component
accuracy. flow prediction, accurate handling of complex ge-
Asafirsttestofthisapproach,atwo-dimensional ometryis required. Moreover,mostusagecasesfor
static droplet is simulated with the variable ini- multi-componentflowrequireaccuratetreatmentof
tial droplet radius, R = {4,8,16,24,32,48}, and solidwallpropertiessuchasno-slipandwettability.
relaxation time, τα = {0.525,0.55,1.0,1.5,3.0}for InthestandardLBM,theno-slipconditionisre-
eachcomponent. Thesimulationdomainsizeisfive alized with pointwise particles’bounce back model
times the droplet radius and the initial density for [17]. In cases with coarse resolution or low vis-
each component is 0.22. After a steady state is cosity, accuracy may deteriorate. When the geom-
reached, the droplet radius is measured by fitting etry or the local shear velocity is under-resolved,
the hyperbolic tangentcurveto the densityprofile. theshearstressalongthe wallcannotbe estimated
In Fig.1, the pressure differences across the accurately. Numerical smearing could also easily
droplet interface, dP, are plotted with respect to contaminateflowfieldinthe nearwallregionwhen
the inverse droplet radius 1/R, using four sets of the physical viscosity is low. As a result, the nu-
τα combinationswith the maximumviscosityratio merically simulated absolute permeability in the
of 100. The subscripts 1 and 2 for τ denote quan- porousmediamayhavesignificantviscositydepen-
tities inside and outside the droplet, respectively. dence even if the Reynolds number is low enough
Results for all τ options are fitted by a line. Ac- to satisfy the Darcy’s law approximately [18]. In
cording to the Young-Laplace law, the previous studies [19, 20], the multiple relax-
ation time (MRT) scheme along with the interpo-
σ
dP = , (12) lated and multi-reflected boundary conditions im-
R
proved this issue. However those schemes are diffi-
the slope of the fitted line is the numerically culttogeneralizeandimplementforcomplexprob-
achievedsurfacetensionσ, whichisindependentof lems. Both of these types of boundary conditions
3
0.0035 0.0035
0.003 0.003
0.0025 0.0025
0.002 0.002
P P
d d
0.0015 0.0015
0 0.0.000015 ττ11 ==ττ1103 ==..50312 τ..5002 =τττ2220 ===.53312...5000 0 0.0.000015 ττ11 ==ττ1103 ==..50132 ..τ0052 =τττ2220 ===.51332...0050
Best fitting line Best fitting line
0 0
0 0.04 0.08 0.12 0.16 0 0.04 0.08 0.12 0.16
1/R 1/R
Figure 1: Pressuredifference across the droplet interface, dP, as a function of the inverse droplet radius, 1/R, using the
originalBGKschemes (left) andnew schemes (right)withfourcombinations ofrelaxation times. Thesubscripts 1and2
ofdenote τ insideandoutsideofthedroplet,respectively. Alineisfittedbasedonresultswithfourcombinationsofτ.
s current 000...000111246 ττ11 = =ττττ011110 .====5.5132132..5..00500 ( (((((PPPOOOrrrerrreeiiisgggsseiiieennnnnnaaattt)lll))))) s current 000...000345 ττττ1111 ====ττ110303 ==....5050112,2,. . 05τ05τ22,,,, ==ττττ122200= ===.. 55131322....005005 ((((((PPPOOOrrrrrreeeiiigggsssiiieeennnnnnaaatttlll))))))
u u
o 0.01 o
uri uri 0.02
p 0.008 p
s s
x x
Ma 0.006 Ma 0.01
0.004
0
0.5 1 2 3 4 8 12 16 20 24
τ R
2
Figure2: Maximumspuriousvelocitywiththeoriginalandmodifiedschemeasafunctionofτ2 (left)andR(right)using
variousviscositycombinations. Intheleftfigure,theinitialRis48.
4
Figure 3: Color contours of ρ2 and the velocity field using the original scheme (right) and the modified scheme (left).
τ1=0.525andτ2=3.0. InitialRis48.
do not conserve the local mass and require com- 4.1. Modified wall models
plicated local numerical interpolations in the near
In our work,the bounce back model is an exten-
wall regions. Thus a more robust, accurate, and
sionofthevolumetricboundaryconditionproposed
simple model is desired.
by Chen et al in 1998 [21–24], which has been ex-
tensivelystudiedforarbitrarygeometry. Themain
Thewettabilityconditionisoftenmodeledbyan features of this model are:
interaction force between the fluid and the wall. A
wall potential ρα is assignedto eachwall to enable 1. Boundary surfaces are discretized into piece-
s
wise linear surface facets in two dimensions
interactionsbetweenthesolidsurfaceandfluidpar-
and triangular polygons in three dimensions;
ticles using a point-wise concept along the lines of
Eq. (7). The correspondence between ρα and the 2. During the fluid dynamics calculation, the
s
contact angle can be defined by simulating some facets/polygons gather incoming particles
test cases, after which it can be used for general from neighboring cells in a volumetric way;
cases. Although this wettability model works well 3. Atthe wall,the directionofthe incoming par-
for certain cases [12], it could sometimes generate ticleisflippedandtheoutgoingparticleiscon-
an artificial thin film along a wetting wall. This structed;
film originates from slugs or droplets on the wall, 4. The outgoing particles are scattered back to
andthisnumericalartifactisobviouslyundesirable neighbor cells in the similar volumetric way.
because large amounts of mass may artificially dif-
fuse and escape from the inlet or outlet boundary. More details can be found in [21], in particular a
Another well known undesirable issue is the artifi- proof that conservation laws are obeyed locally as
cial movement of static droplets on inclined walls. well as globally.
Because of insufficient isotropy of discrete numer- In addition to the correction of surface scatter-
ical schemes, even if the droplet is not subject by ing described in [22], a hybrid solid wall boundary
anydrivingforce,thedropletmayeasilydescendor conditionforthenoslipwallisproposedhereinor-
even climb the slope. To address these issues, the dertofurtherreducenumericalsmearingforcoarse
improved wall boundary conditions are developed resolution simulations. The distribution function
and validated. of outgoing particles fout,α is a combination of
i
5
bouncedbackparticlesandtheMaxwellianequilib- long as the flow Reynolds number is small. In our
riumparticledistributioninaccordancetoviscosity simulations, the absolute permeability is evaluated
values. as K = φνU/g, where φ is the porosity and U is
In our wettability model, the interaction force the spatial averagedvelocity.
Fα,β is extended from the inter-component force A streamline inthe case ofτ =1 is presentedon
w
formgivenby Eq. (7) using the volumetric bound- the left side of Fig. 4. The color represents the ve-
ary scheme as, locitymagnituderangingfrom0to0.0025. Around
the regionwherethe spheresaretouching,the flow
Fαw,β(x)=Gρα(x) wiciρ′β(x+ci∆t), (13) pathisofanarrowconcaveshapeandthereforethe
Xi flowisfast. Insucharegion,theresolutiontendsto
be relatively coarse and numerical smearing could
′
where ρ is constructured in a volumetric way so
β be quite pronounced when fluid viscosity is low.
that∂ρ /∂n=0[21]. Itisworthpointingoutthat
β On the right side of Fig. 4, the simulated abso-
such a volumetric wettability scheme has sufficient
lutepermeabilitywiththeoriginalbouncebackand
isotropy in complex geometry.
the currently presented models are shown, along
with the analytical solution obtained through the
4.2. Results
analytical dimensionless drag force from [19, 20,
Themodifiedmodelsarevalidatedinordertoes- 25]. In all simulations the Reynolds number UL/ν
tablishthattheissuesspecifiedaboveareimproved. is less than 0.4. The absolute permeability simu-
First, the kinematic viscosity ν dependence of the lated using the original model shows a measurable
absolute permeability is investigatedby simulating dependence upon viscosity, consistent with previ-
the simple cubic (SC) array of spheres. Second, ous studies by other researchers[18]. On the other
the thin film along the wall is tested on the case hand,themodifiedmodelreportedinthisworksig-
of a static slug between plates. Third, the artifi- nificantlyreducesthisdependenceevenatverylow
cial movement of a droplet on an inclined wall is viscosities. Simulation results agree well with the
studied. analytical solution. Thus the new models can be
applied for the complex geometry including sharp
4.2.1. The SC array of spheres convexandconcaveshapes,whilemaintaininghigh
The single-component gravity driven flow accuracy.
through the SC array of spheres is simulated.
The simulation domain is set as a cube with the 4.2.2. Static slug between flat plates
edge length L = 34, containing a sphere in the Atwo-dimensionalstaticslugbetweenflatplates
center. Periodic boundaries are assigned to each is simulatedunder variouswetting conditions. The
pair of faced cube surfaces. The relative volume channel height is set at 32 and τ for both compo-
fraction χ = (c/c )1/3 = 1 and the spheres are nentsis1. Thesecondcomponentismainlylocated
max
touching each other, where c is the ratio of the in the middle section of the channel. The wall po-
solid volume to the cube volume and c is its tential ρα is set as ρ1 = 0, ρ2 = 0.088 so that
max s s s
maximum value, π/6 in the SC array. The choices the contact angle is roughly 40 degrees. In Fig.
of τ are {0.505,0.51,0.55,0.6,1.0,1.5,2.5}. The 5, the color contours show the second component
gravity g is in the perpendicular direction to a density distribution with the originaland modified
cube surface and g = 1.e − 4 when ν = 0.1666. wettability models, respectively. It is seenthat the
The value of g for the other values of ν is adjusted modified scheme reduces the thin film artifact sig-
so that g/ν2 is kept constant. Such a choice keeps nificantly without changing the contact angle.
the Reynolds number constant for the case of the In Fig. 6, the contact angle and the thin film
Hagen-Poiseuille flow when the channel height is density are shown as a function of wall potential.
1 2 1
fixed. UnderthelowReynoldsnumberassumption, The wall potentials ρ and ρ are varied as ρ =
s s s
the Darcy’s law is obeyed and the flow is within −ρ0·α·Θ(−α)andρ2s =ρ0·α·Θ(α). Theparameter
the Stokes flow regime. Therefore for a given α is actually a dimensionless wall potential, Θ is
geometry, the absolute permeability K should be the Heaviside step function, and ρ0 = 0.22. The
a constant independent of physical properties as contact angle was measured by fitting a circle to
6
4
Present
Original
y 3.6 Analytical
bilit
a 3.2
e
m
er 2.8
p
e
ut 2.4
ol
s
Ab 2
1.6
0.01 0.1 1
ν
Figure4: Geometryandflow stream linescoloredbyvelocitymagnitude rangingfrom0to0.0025(left), andtheabsolute
permeabilityasafunctionofthekinematicviscosityν withthepresentandoriginalmodel(right). Theanalyticalsolution
forpermeabilityisshownbythedotted line.
the interface. Thin film density is detected in the whiletheonewiththeoriginalmodelismoving. In
first lattice cell of the simulation domain’s edge. Table 1, results for various τ combinations and in-
Observe from Fig. 6 that the our modification of clinationanglesarepresented. Thedropletvelocity
the wallmodels does notchangethe contactangle, is calculated by measuring the position change of
as intended, but at the same time reduces the thin the droplet center mass between 8.e+4 and 8.5e+4
filmandwettabilitydependenceofthethinfilmfor time steps. Therefore the velocity less than 2.e-
all cases. 4 cannot be measured accurately and is set to 0.
Utilizing the ratio of the droplet volume to total
4.2.3. A droplet on the inclined wall volume, the mass diffusion can be evaluated and
A two-dimensional droplet that is not subject thequantitativecomparisonbetweenthetwomod-
to any explicit driving force on an inclined wall is els is presented. The droplet volume is detected
simulated. The droplet is composed by the sec- by adding volume where the second component
ond component, surrounded by the first compo- density is more than 0.18. The results show that
nent. The channel height is 16 and the wall po- the modified model improves the artificial droplet
tential is set similarly to the slug case above. The movement and mass diffusion on the inclined wall
periodic boundaries are enforced on the top/left for all tested inclinations.
andbottom/rightedges. Theinclinationanglesare
{10,30,50,70}degrees. In the caseswith the same
inclination angle, exactly the same initial droplet 5. Summary
mass is set. In Fig. 7 and Fig. 8, the results ob-
tained with the original and modified models are An enhanced multi-component LB flow solver is
comparedusing the inclinationangles of 10 and 30 presentedwhoseaccuracyandstabilityaredemon-
degrees, for τ1 = 0.55 and τ2 = 1.5. Fig. 7 shows strated on a set of difficult benchmark cases. The
that the droplet is less diffusive with the modified algorithm and numerical scheme are generalized
models, which is likely due to the reduction of the for practical applications. The surface tension
thin film effect that has been shown for the previ- achieved in the simulation is independent of fluid
ously discussed cases. In Fig. 8 the droplet simu- viscosity and resolution. The spurious current is
lated with the modified model is static, as desired, significantly reduced. The new surface wetting
7
Figure5: Densityofthesecondcomponent,ρ2,withtheoriginal(left)andmodifiedwettabilitymodel(right). Thechannel
heightis32andτ ofbothcomponents are1.
180 0.024
Present Present
160 Original Original
e) 0.02
Degre 112400 nsity 0.016
e ( 100 de
gl m 0.012
act an 6800 hin fil 0.008
nt 40 T
o 0.004
C 20
0 0
-0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4
α α
Figure 6: The static contact angle and the thin film density as a function of α with original and modified models. The
wallpotentials ρ1s andρ2s arevariedfollowingρ1s =−ρ0·α·Θ(−α)andρ2s =ρ0·α·Θ(α)whereΘistheHeavisidestep
functionandρ0 is0.22. τ ofbothcomponents are1.0.
Table1: Dropletvelocityandtheratioofdropletvolumetototal volumeontheinclinedwallwithrespecttocomponent
τ andtheinclinationanglewiththeoriginalmodel(left)andthemodifiedmodel(right).
τ1 τ2 Inclination Droplet Droplet volume τ1 τ2 Inclination Droplet Droplet volume
angle velocity /total volume angle velocity /total volume
0.55 1.5 10 0 0.077 0.55 1.5 10 0 0.25
0.55 1.5 30 6.e-4 0.25 0.55 1.5 30 0 0.35
0.55 1.5 50 5.e-4 0.21 0.55 1.5 50 0 0.33
0.55 1.5 70 3.e-4 0.27 0.55 1.5 70 0 0.40
1.0 1.0 10 0 0.078 1.0 1.0 10 0 0.24
1.0 1.0 30 4.e-4 0.23 1.0 1.0 30 0 0.31
1.0 1.0 50 0 0.20 1.0 1.0 50 0 0.28
1.0 1.0 70 3.e-4 0.26 1.0 1.0 70 0 0.37
1.5 0.55 10 0 0.057 1.5 0.55 10 0 0.24
1.5 0.55 30 3.e-3 0.21 1.5 0.55 30 0 0.26
1.5 0.55 50 5.e-4 0.19 1.5 0.55 50 0 0.26
1.5 0.55 70 4.e-4 0.24 1.5 0.55 70 0 0.34
8
Figure7: Thedensityofthesecondcomponentatdifferenttimestepswiththeoriginal(left)andmodified(right)models.
The density range is from 0 to 0.22. The angle of the inclined wall is 10 degrees; τ1 = 0.55 and τ2 = 1.5. No explicit
drivingforceisapplied. Theperiodicboundariesareenforcedbetween thetop/leftandbottom/rightedges.
Figure8: Thedensityofthesecondcomponentatdifferenttimestepswiththeoriginal(left)andmodified(right)models.
Thedensityrangeisfrom0to0.22. Theangleoftheinclinedwallis30degrees;τ1=0.55andτ2=1.5. Nodrivingforce
isapplied. Theperiodicboundariesareenforcedbetween thetop/leftandbottom/right edges.
scheme for complex geometry improves the near [5] H. Chen, S. Chen, W.H. Matthaeus, Recov-
wall algorithm isotropy and the overall quality of ery of the Navier-Stokes equations using a
numerical solution in the near wall region, in par- lattice-gas Boltzmann method, Phys Rev A
ticular it reduces the unphysical surface thin film 45(1992)R5339
and mass diffusion. The model enables simulation
of multi-component flows in an extended viscos- [6] X. Shan, H. Chen, Simulation of non-ideal
ity range and in complex geometry with much im- gases and liquid-gas phase transitions by
proved accuracy, stability, and robustness. lattice Boltzmann equation, Phys Rev E
49(1994)2941
References
[7] Q. Li, K.H. Luo, X.J. Li, Forcing scheme in
[1] X. Shan, H. Chen, Lattice Boltzmann model pseudopotential lattice Boltzmann model for
for simulating flows with multiple phases and multiphase flows, Phys Rev E 86(2012)016709
components, Phys Rev E 47(1993)1815
[8] K. Connington, T. Lee, A review of spurious
[2] Y. Qian, D. d’Humi´eres, P. Lallemand, Lat- currents in the lattice Boltzmann method for
tice BGK models for Navier-stokes equation, multiphase flows, Journal of mechanical sci-
Europhys.Lett 17(1992)479 ence and technology 26(2012)3857-3863
[3] P.L. Bhatnagar, E.P. Gross, M. Krook, A
[9] X. Shan, Analysis and reduction of the
model for collisions in gases I. Small ampli-
spurious current in a class of multiphase
tude processes in charged and neutral one-
lattice Boltzmann models, Phys Rev E
component systems, Phys Rev 94(1954)511-
77(2006)047701
525
[4] S.Chen,H.Chen,D.Martnez,W.Matthaeus, [10] M. Sbragaglia,R. Benzi, L. Biferale, S. Succi,
Lattice Boltzmann model for simulation K. Sugiyama, F. Toschi, Generalized lattice
of magnetohydrodynamics, Phys Rev Lett Boltzmannmethodwithmultirangepseudopo-
67(1991)3776 tential, Phys Rev E 75(2007)026702
9
[11] R. Zhang, X. Shan, H. Chen, Efficient kinetic [22] Y. Li, R. Zhang, R. Shock, H. Chen, Pre-
methodforfluidsimulationbeyondtheNavier- diction of vortex shedding from a circular
Stokes equation, Phys. Rev. E, 74, (2006) cylinderusingavolumetricLattice-Boltzmann
046703 boundary approach, Eur.Phys.J.Special Top-
ics 171(2009)91-97
[12] H. Otomo, H. Fan, R. Hazlett, Y. Li, I.
Staroselsky,R.Zhang,H. Chen, Simulationof [23] Y.Li,R.shock,R.Zhang,H.Chen,Numerical
residual oil displacement in a sinusoidal chan- studyofflowpastanimpulsivelystartedcylin-
nel with the lattice Boltzmann method, C R der by the lattice-Boltzmann method, J.Fluid
Mecanique 343(2015)559-570 Mech 519(2004)273-300
[13] X. Shan, X. Yuan, H. Chen, Kinetic the-
[24] H. Fan, R. Zhang, H. Chen, Extended volu-
ory representation of hydrodynamics: a way
metric scheme for lattice Boltzmann models,
beyond the Navier-Stokes equation, J.Fluid
Phys Rev E 73(2006)066708
Mech 550(2006)413-441
[25] A.S. Sangani and A. Acrivos, Slow flow
[14] H.Chen,X.Shan,Fundamentalconditionsfor
throughaperiodicarrayofspheres,IntJMul-
N-th-orderaccuratelatticeBoltzmannmodels,
tiphase Flow 8(1982)343-360
Physica D 237(2008)2003-2008
[15] H. Chen, R. Zhang, I. Staroselsky, M. Jhon,
Recovery of full rotational invariance in lat-
ticeBoltzmannformulationsforhighKnudsen
number flows, Physica A 362 (1) (2006), 125
[16] J. Latt, B. Chopard, Lattice Boltzmann
methodwithregularizedpre-collisiondistribu-
tionfunctions,Math.Comput.Simulat.72(2-
6) (2006), 165
[17] S. Chen, G. Doolen, Lattice Boltzmann
method for Fluid Flows, Annual Rev. Fluid
Mechs., 30 (1998), 329-364
[18] H. Liu, Q. Kang, C.R. Leonardi, B.D. Jones,
S. Schmieschek, A. Narva´ez, J.R. Willianms,
A.J. Valocchi, J. Harting, Multiphase lattice
Boltzmann simulations for porous media ap-
plications, arXiv:1404.7523
[19] C. Pan, L. Luo, C.T. Miller, An evaluation of
latticeBoltzmannschemesforporousmedium
flow simulation, Comput. Fluids 35(2006)898-
909
[20] E. Fattahi, C. Waluga, B. Wohlmuth, U.
Rude, M. Manhart, R. Helmig, Pore-scale
lattice Boltzmann simulation of laminar
and turbulent flow through a sphere pack,
arXiv:1508.02960
[21] H. Chen, C. Teixeria, K. Molving, Realiza-
tion of fluid boundary conditions via dis-
crete Boltzmann dynamics, Int J Mod Phys
C 9(1998)1281
10