- Email: [email protected]
From mesoscopic models to continuum mechanics: Newtonian and non-newtonian fluids
Accepted Manuscript
From Mesoscopic Models to Continuum Mechanics: Newtonian and Non-Newtonian Fluids R.R. Huilgol, G.H.R. Kefayati PII: DOI: Reference:
S0377-0257(16)30003-9 10.1016/j.jnnfm.2016.03.002 JNNFM 3762
To appear in:
Journal of Non-Newtonian Fluid Mechanics
Received date: Revised date: Accepted date:
1 September 2015 2 March 2016 3 March 2016
Please cite this article as: R.R. Huilgol, G.H.R. Kefayati, From Mesoscopic Models to Continuum Mechanics: Newtonian and Non-Newtonian Fluids, Journal of Non-Newtonian Fluid Mechanics (2016), doi: 10.1016/j.jnnfm.2016.03.002
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Highlights
AC
CE
PT
ED
M
AN US
CR IP T
• Critical examination of the BGK approximation to derive constitutive equations for Newtonian fluids • Limitations due to grid dependence of viscosity in Newtonian and nonNewtonian fluids • A new evolution equation for the lattice distribution function • Derivation of the equations of continuum mechanics for Newtonian and nonNewtonian fluids from the new model
1
ACCEPTED MANUSCRIPT
From Mesoscopic Models to Continuum Mechanics: Newtonian and Non-Newtonian Fluids R. R. Huilgol ∗ and G. H. R. Kefayati
CR IP T
School of Computer Science, Engineering and Mathematics, Flinders University of South Australia, GPO Box 2100, Adelaide, SA 5001, Australia
Abstract
PT
ED
M
AN US
A review of the BGK approximation to obtain the equations of motion for an incompressible fluid is presented and its drawbacks are revealed. In order to overcome these inherent problems, new models for the particle distribution functions are needed. Using the Finite Difference Lattice Boltzmann Method (FDLBM) due to Fu and So [1] and the Thermal Difference Discrete Flux Method (TDDFM) proposed by Fu et al [2], it is shown that the newer distribution functions lead to the mass conservation equation, the equations of motion and the energy balance equation for incompressible fluids in two dimensions, employing the D2Q9 lattice as the model. This derivation is extended to compressible fluids as well. Next, using the D3Q15 lattice as an example, the three dimensional equations of continuum mechanics are derived. Since no restrictions are placed on the constitutive equations, the theoretical development applies to all fluids, whether they be Newtonian, or power law fluids, or viscoelastic and viscoplastic fluids. Finally, some comments are offered regarding the numerical scheme to calculate the particle distribution functions to determine the velocity and temperature fields.
AC
CE
Key words: Lattice Boltzmann equation, BGK approximation, particle distribution function, continuum mechanics
1
Introduction
From the Lattice Boltzmann equation, it is possible to derive the continuity equation and Cauchy’s equations of motion for a compressible medium, when ∗ Corresponding author. Fax: +61 8 8201 2904. Email address: [email protected] (R. R. Huilgol).
Preprint submitted to Elsevier
22 March 2016
ACCEPTED MANUSCRIPT one uses the Bhatnagar-Gross-Krook (BGK) approximation. From this, one can obtain equations relevant to incompressible fluids. However, these require that the pressure be proportional to the density and the viscosity be dependent on the collision relaxation time [3]; see Section 2 below. Clearly, these restrictions on the pressure and the viscosity are unacceptable in modelling the flows of non-Newtonian, incompressible fluids.
M
AN US
CR IP T
In order to overcome these inherent problems, new models for the particle distribution functions are needed. In the Finite Difference Lattice Boltzmann Method (FDLBM) due to Fu and So [1], the particle distribution function leads to the conservation of mass and the equations of motion applicable to incompressible fluids, when the flows are assumed to occur in a two dimensional setting underpinned by a D2Q9 lattice. Our derivation of these results is succinct and is more transparent, for it uses using vector analysis and linear algebra. In addition, the energy equation is also obtained from Thermal Difference Discrete Flux Method (TDDFM) proposed by Fu et al [2]; once again, simple results from vector analysis and linear algebra are employed. The important point to note is that the previous restrictions on the pressure and the viscosity are eliminated in these derivations, meaning that one is free to choose a constitutive equation. That is, one can model a Newtonian fluid, or power law fluids, or viscoelastic and viscoplastic fluids. Moreover, we point out in sub-section 3.1 that it is quite easy to incorporate Dirichlet type boundary conditions into the numerical scheme to determine the particle distribution functions for the velocity and temperature fields.
ED
In Section 4, the method is extended to the flows of all fluids in three dimensions, using the D3Q15 lattice as an example. Once again, the particle distribution functions are such that every type of fluid, compressible or incompressible, can be employed.
BGK Approximation to Continuum Mechanics
AC
2
CE
PT
Finally, in Section 5, some remarks are offered on the numerical scheme employed to determine the particle distribution functions for the velocity and the temperature fields. The CFL condition for the stability of the numerical scheme is also derived.
Beginning in 1986, Lattice Boltzmann equation (LBE) models evolved from their Boolean counterparts, viz., the lattice-gas-automata. The theoretical framework of the LBE was underpinned by the Chapman-Enskog analysis of the LGA models (Frisch et al [4,5]; Wolfram [6]). That is, the statistical mechanics of the LGA played a crucial role in these developments. A decade later, He and Luo [7] proved that the LBE is a specific discretised form of 3
ACCEPTED MANUSCRIPT the continuous Boltzmann equation using the Bhatnagar-Gross-Krook (BGK) approximation [8]. As an application, He and Luo [7] derived the Lattice Boltzmann BGK (LBGK) equation for the D2Q9 lattice model. It is worth noting that Welander [9] published an approximation to the Boltzmann equation at the same time as BGK and applied it to a rarefied monatomic gas exchanging heat with an adjacent wall. In several aspects, his treatment of this problem is exemplary and anticipates by several decades the subsequent Lattice formulation. However, we shall use the more commonly used acronym LBGK equation in the sequel.
CR IP T
From the LBGK equation, Cauchy’s equations of motion for compressible fluids can be derived under the assumption that the Mach number is small and that the density variation is small. Attempts to derive the equations for incompressible fluids from the foregoing end in an physically unacceptable requirements, viz., that the pressure be proportional to the density and that the viscosity depends on the collision relaxation time [3].
AN US
To render these remarks transparent, let f = f (x, ξ, t) be the probability of finding a particle with the velocity ξ near the point x at time t. The Boltzmann equation for this distribution function, in the absence of external forces, is given by
ED
M
1 ∂f + ∇f · ξ + f = N (f ), (2.1) ∂t λ where ξ is the microscopic velocity, λ is the relaxation time due to collision and N (f ) is a measure of the net number of molecules which disappear from the region due to inter-molecular collisions; see Welander [9], for example. This is not applied to the problem at hand, for it depends upon an approximation to this measure of collisions; see (2.3) below.
PT
The macroscopic variables are the density of mass ρ, the velocity field u and the absolute temperature T. These are related to the distribution function f through the following integrals in momentum space:
CE
Z
ρ=
f dξ,
ρu =
Z
ξf dξ,
1Z ρε = |ξ − u|2 f dξ, 2
(2.2)
AC
where ε = D0 RT /2, with D0 being the number of degrees of freedom of a particle and R is the ideal gas constant. Note that we have used dξ = dξ1 · · · dξD , where D ∈ [2, 3] stands for the dimension of the physical space. In the BGK approximation, the collision operator is assumed to be given by N (f ) = g(x, ξ, t)/λ, where the collision function g(x, ξ, t) has the form: "
#
ρ |ξ − u|2 g(x, ξ, t) = exp − . (2πRT )D/2 2RT 4
(2.3)
ACCEPTED MANUSCRIPT Next, using (2.3) and integrating (2.1) formally over a time interval 4t with et /λ as the integrating factor, one obtains [7]: f (x + ξ4t, ξ, t + 4t) = e−4t/λ f (x, ξ, t) Z 4t 1 0 et /λ g(x + ξt0 , ξ, t + t0 ) dt0 . + e−4t/λ λ 0
(2.4)
Assuming that the time interval 4t is small enough and linearising in terms of this time interval, one has:
CR IP T
1 f (x + ξ4t, ξ, t + 4t) − f (x, ξ, t) = − [f (x, ξ, t) − g(x, ξ, t)], τ
(2.5)
where τ = λ/4t is the non-dimensional relaxation time. Next, suppose that the collision function g(x, ξ, t) can be expanded as a Taylor series in u, retaining terms up to order |u|2 . Identifying this as the equilibrium distribution function f (0) , the following can be derived from (2.3): "
AN US
f
(0)
#
|u|2 ρ ξ · u (ξ · u)2 2 exp(−|ξ| /2RT ) + − = 1 + . (2πRT )D/2 RT 2(RT )2 2RT
(2.6)
To derive Cauchy’s equations of motion for a continuous medium, exact evaluation of the following momentum integrals m (0)
|ξ| f
dξ =
Z
exp(−|ξ|2 /2RT )ψ(ξ) dξ,
M
I=
Z
m = 0, 1, 2,
(2.7)
PT
ED
is required. Using a Cartesian coordinate system and noting that the physical dimension of the space D = 2, one can express the microscopic velocity as ξ = (ξx , ξy ). Using this, ψ(ξ) is expressed as a polynomial of the form ψ(ξ) = ξxr ξys .
CE
The integrals in (2.7) can be evaluated through a Gaussian-type quadrature [10] and lead to X I= Wα exp(−|ξ α |2 /2RT )ψ(|ξ α |), (2.8) α
AC
where Wα are the weights and ξ α are the discrete velocities of the quadrature. In the D2Q9 model, there are nine discrete velocities: {ξ α , α = 0, 1, · · · , 8}. The configuration space is discretised into a square lattice with a lattice con√ stant 4x = 3RT 4t. In fact, in an isothermal problem, the temperature T has no physical significance √ and thus, one can choose 4x as a fundamental quantity instead, leading to 3RT = c = 4x/4t, where the lattice speed c is related to the speed of sound cs through c2 = 3c2s . Thus, RT = c2s = c2 /3. The nine lattice points are located at the centre of the lattice, at the four midpoints of the edges of the square and at the four vertices. One identifies 5
ACCEPTED MANUSCRIPT α0 = (0, 0) as the centre of the lattice, {α1 , α3 , α5 , α7 } as the midpoints of the four sides: (±1, 0), (0, ±1), and {α2 , α4 , α6 α8 } = (±1, ±1) as the vertices. The lattice velocity vector is given by ξ 0 = 0. Next, the velocities for α = 1, 3, 5, 7, have the form c(±1, 0), c(0, ±1), and for α = 2, 4, 6, 8, they have the form c(±1, ±1). In other words, except at the centre, at each lattice point αj , j = 1, · · · , 8, the velocity vector points away from the lattice in the direction from α0 to αj .
where
wα =
4/9,
α = 0,
1/9, α = 1, 3, 5, 7,
1/36, α = 2, 4, 6, 8.
CR IP T
The weights Wα in (2.8) can now be computed through quadrature and one obtains [7]: Wα = 2πRT exp(|ξ α |2 /2RT )wα , (2.9)
(2.10)
AN US
The equilibrium distribution function for the D2Q9 model is given by fα(0) = ρwα + ρsα (u(x, t)), where
with
(0, 0),
(2.12)
α = 0,
c(cos θα , sin θα ) α = 1, 3, 5, 7, √ c 2(cos θα , sin θα ), α = 2, 4, 6, 8.
(2.13)
ED
ξα =
#
3(ξ α · u) 9(ξ α · u)2 3|u|2 + − , c2 c4 2c2
M
sα (u) = wα
"
(2.11)
PT
Here, the angles θα are defined through θα = (α − 1)π/4, α = 1, · · · , 8.
CE
Introducing a Mach number M = |u|/cs , one can see that the equilibrium distribution function can be written in the form: fα(0) = fα0 + Mfα1 + M2 fα2 ,
(2.14)
AC
where
fα0 = ρwα ,
fα1 = ρwα .
(2.15)
Assuming that the distribution function fα can be found from the equation corresponding to (2.5), the momentum space integrals in (2.2) lead to the following results for the macroscopic variables: ρ=
X α
fα ,
ρu =
X
fα ξ α ,
α
6
ρε =
1X fα |ξ α − u|2 . 2 α
(2.16)
ACCEPTED MANUSCRIPT When multi-scaling using the Mach number is employed, these equations lead to Cauchy’s equations for a compressible continuous medium of the following form: ∂ρ + ∇ · (ρu) = 0, ∂t
(2.17)
∂(ρu) + ∇ · (ρu ⊗ u) = −∇p + ν[∇2 (ρu) + ∇(∇ · (ρu))]. ∂t
(2.18)
AN US
CR IP T
While the first equation delivers the correct form of the equation for the conservation of mass, the second leads to the equations of motion for an artificial compressible fluid because the pressure p = ρc2s has a thermodynamic form only, with cs being the speed of sound. Of course, this form of the pressure is not surprising due to the absence of any external forces in the original formulation (2.1). However, the expression for the kinematic viscosity is unusual for it is given by ν = (2τ − 1)c2 4t/6, and depends on the non-dimensional relaxation time τ and the time step 4t.
M
Turning to incompressible fluids, we note that the pressure p must be independent of the density. To derive such a model, Guo et al [3] found that a new type of distribution function gα (x, t) is required. The important point is that this new distribution function is patterned along the lines of the BGK model; however, there is no statistical mechanical basis for this formulation. To understand this, one notes that the new equilibrium distribution function gα(0) (x, t)is defined through α = 0,
λ(p/c2 ) + sα (u)
α = 1, 3, 5, 7,
γ(p/c2 ) + sα (u),
α = 2, 4, 6, 8,
−4σ(p/c2 ) + sα (u),
ED
gα(0) =
(2.19)
PT
where σ, λ and γ are parameters satisfying
CE
λ + γ = σ,
λ + 2γ = 1/2.
(2.20)
This new distribution function satisfies the following conservation laws:
AC
8 X
gα =
α=0
8 X
α=0
ceα gα =
8 X
α=0 8 X
gα(0) ,
(2.21)
ceα gα(0) .
(2.22)
α=0
The evolution equation for the system is given by 1 gα (x + ceα 4t, t + 4t) − gα (x, t) = − [gα (x, t) − gα(0) (x, t)]. τ 7
(2.23)
ACCEPTED MANUSCRIPT The velocity and pressure of the flow are given by
u=
8 X
ceα gα ,
(2.24)
α=0 " 2
#
8 c X gα + s0 (u) . p= 4σ α=1
(2.25)
∇ · u = 0,
CR IP T
From the chosen distribution function, a multi-scale expansion as described in full in the Appendix of the paper by Guo et al [3], leads to the following set of equations applicable to an incompressible fluid: (2.26)
∂u + ∇ · (u ⊗ u) = −∇p + ν∇2 u. ∂t
(2.27)
ν=
AN US
While p can now be arbitrary, the kinematic viscosity ν is still relaxation time and grid-dependent, for it is given by (2τ − 1) (4x)2 · . 6 4t
(2.28)
CE
PT
ED
M
In this context, it is necessary to recall that in addition to the work cited earlier [4,5], attempts to derive the Navier-Stokes equations for compressible or incompressible fluids have been made [11–16]. In addition, thermal stability of the hydrodynamic equations has also been investigated [17]. It is important to note that hydrodynamic equations deal, at most, with a restricted class of fluids, viz., the ideal gas or fluids with constant viscosity. Secondly, to obtain the equations for a Newtonian fluid, second and higher order lattice Boltzmann schemes have to employed [15,16]. The derivations given here in Section 3 are based on the results in [1,2] and are applicable to all fluids, compressible or incompressible. Secondly, the equations of continuum mechanics appear at the expansion of the particle distribution functions at the first order, as described in full in (3.3) - (3.50) below.
AC
Turing our attention to non-Newtonian fluids, we note that the relationship between the relaxation time and the viscosity means that the relaxation time has to be prescribed over a specific range only; hence, the viscosity function is restricted and cannot describe a wide range of non-Newtonian fluid models. This is clear from the previous studies employing LBM for non-Newtonian fluids; for example, see the work on power-law fluids [18–25] and Bingham fluids [26]. Even though a second order LBM scheme has been proposed in [19] to solve problems in power-law fluids, the kinematic viscosity is still related to the relaxation time. Hence, a different approach to deriving the equations 8
ACCEPTED MANUSCRIPT of continuum mechanics applicable to all fluids is needed, and we turn to this next.
3
Thermal Difference Discrete Flux Method to Continuum Mechanics
CR IP T
First of all, the lattice counterpart of (2.1) with the assumption that the collision term N (f ) = f eq /λ can be written as follows: 1 ∂fα + ξ α · ∇ x fα = − (fα − fαeq ), ∂t ετ
(3.1)
where ε is a small parameter to be prescribed when numerical simulations are considered. To proceed, one assumes that fα has the following expansion:
AN US
fα = fαeq + εfα(1) + ε2 fα(2) + O(ε3 ).
(3.2)
The novelty of the approach by Fu and So [1] lies in expanding the equilibrium lattice function fαeq as a quadratic in the particle velocity ξ α : fαeq = Aα + ξ α · Bα + (ξ α ⊗ ξ α ) : Cα ,
(3.3)
8 X
fαeq = ρ,
α=0
fαeq ξ α = ρu,
(3.4)
u = ui + vj,
PT
8 X
ED
M
where Bα is a vector and Cα is a 2 × 2 symmetric matrix. Before defining these entities, the following relations must hold (cf. (2.16)):
(3.5)
α=0
fαeq ξ α ⊗ ξ α = M,
CE
8 X
α=0
8 X
fα(n) = 0,
n ≥ 1,
(3.7)
fα(n) ξ α = 0,
n ≥ 1.
(3.8)
AC
α=0
8 X
α=0
(3.6)
In (3.6), M has the matrix form
2 ρu + p − τxx
M=
ρuv − τxy
9
ρuv − τxy . ρv 2 + p − τyy
(3.9)
ACCEPTED MANUSCRIPT In the above set, ρ is the density, u and v are the components of the velocity field u in the x and y directions respectively, τxx , τxy = τyx , τyy are the stresses which can be defined through any relevant constitutive relation. For a D2Q9 lattice, in (3.3), there are nine coefficients: Aα , Bα , and Cα . However, when α = 0, only the coefficient A0 is required. Regarding the others, one makes the assumption that the coefficients having the same energy shell of the lattice velocities are equal. Thus, Aα = A1 , α = 1, 3, 5, 7; Aα = A2 , α = 2, 4, 6, 8.
AN US
CR IP T
Since B0 and C0 do not affect the value of f0eq in (3.3), they can be put equal to zero. By the just mentioned assumption about certain coefficients being equal, we see that Bα = B1 , α = 1, 3, 5, 7; Bα = B2 , α = 2, 4, 6, 8. Similarly, Cα = C1 , α = 1, 3, 5, 7; Cα = C2 , α = 2, 4, 6, 8. Keeping in mind that the matrices are symmetric, there are thirteen independent quantities: three scalars in A0 , A1 , A2 ; four components of the vectors B1 , B2 ; three components each in C1 , C2 . Looking at (3.4) - (3.6), it is obvious that the number of available constraints is six only. Thus, there is a lot of latitude in choosing the variables. Following Fu and So [1], we shall assume that A0 = ρ −
2p ρ|u|2 τxx + τyy − 2 + , σ2 σ σ2 ρu , 2σ 2
M
and
B1 =
and
C11 =
PT
C11 0 C1 = , 0 C22
1 (p + ρu2 − τxx ), 2σ 4
(3.10)
(3.11)
C22 =
1 (p + ρv 2 − τyy ), 2σ 4 (3.12)
0 C12 C2 = , C21 0
CE
B2 = 0.
ED
Next,
A1 = A2 = 0,
C12 = C21 =
1 (ρuv − τxy ). 8σ 4
(3.13)
AC
Here, σ is a constant to be determined keeping in mind that the vectors ξ α have been redefined so that in (2.13), c = σ, the lattice speed. Since the value of σ affects numerical stability, its choice depends on the problem at hand. For instance, in their solution to the Stokes Second problem, Fu and So [1] recommend varying it at every iteration step. This matter is briefly addressed in Section 5 below. Before deriving the macroscopic equations of continuum mechanics for incompressible fluids, we note that in several fluids, such as viscoplastic fluids, the 10
ACCEPTED MANUSCRIPT pressure p has to be defined uniquely. This requires that the trace of the extra stress tensor should be zero. For a proof, see [27], [28] and [29]. Thus, in (3.10), τxx + τyy = 0, and, in these models, the coefficient A0 has the simple form: ρ|u|2 2p A0 = ρ − 2 − 2 . σ σ
(3.14)
CR IP T
We shall now demonstrate that the macroscopic equations for an incompressible continuous medium can be derived from (3.1) - (3.9). To make the calculations more transparent, it is helpful to note that
AN US
ρ ρ u + σ 2 C11 , f5eq = − u + σ 2 C11 , f0eq = A0 , f1eq = 2σ 2σ ρ ρ eq eq 2 2 f3 = v + σ C22 , f7 = − v + σ C22 , 2σ 2σ f2eq = f6eq = C12 , f4eq = f8eq = −C12 .
(3.15)
(3.16)
(3.17)
Hence, Eqs. (3.4) - (3.9) follow from the substitution of (3.10) - (3.13) into (3.4) - (3.6), which is as expected. Secondly, substituting the expression for fα in (3.2) into (3.1), we find that (3.18)
M
∂fαeq 1 + ξ α · ∇x fαeq = − fα(1) + O(ε). ∂t τ
ED
Since the velocity vectors ξ α do not depend on the spatial coordinates, one sees that ∇ · (fαeq ξ α ) = ξ α · ∇x fαeq . (3.19) Hence, (3.18) can be rewritten as
PT
∂fαeq 1 + ∇ · (fαeq ξ α ) = − fα(1) + O(ε). ∂t τ
(3.20)
CE
Summing the equation above, we obtain
AC
∂ ∂t
8 X
α=0
fαeq
!
+∇·
8 X
α=0
fαeq ξ α
!
8 1X =− fα(1) + O(ε). τ α=0
(3.21)
Appealing to (3.4), (3.5), (3.7) and recalling that ρ is a constant, the above equation reduces to ρ(∇ · u) = 0 + O(ε), (3.22) or the continuity equation for an incompressible medium is satisfied. Next,
∇ · (fαeq ξ α ⊗ ξ α ) = (ξ α · ∇x fαeq )ξ α . 11
(3.23)
ACCEPTED MANUSCRIPT Hence, multiplying (3.18) throughout by ξ α , summing over α, and appealing to (3.5) - (3.8), one obtains ρ
∂u + ∇ · M = 0 + O(ε). ∂t
(3.24)
Rearranging, one obtains the equations of motion for an incompressible fluid: ρa + ∇p − ∇ · τ = 0 + O(ε),
a=
∂u + (u · ∇)u. ∂t
(3.25)
AN US
CR IP T
Thus, the set of equations (3.1) - (3.13) deliver Cauchy’s equations of motion for an incompressible medium. It is also clear that there is no relationship between the density ρ and the pressure p, which is as it should be in an incompressible fluid. Moreover, one can specify the material through a chosen constitutive equation for the extra stress tensor τ , whether the fluid be Newtonian, or non-Newtonian. At this juncture, it has to be noted that in the equations of motion (3.25), the body force term is missing. It is possible to do so by altering the evolution equation, suggested in the work by Fu et al [2]. If the body force is given by ρb, it has to be incorporated into the evolution equation (3.1) in such a way that (3.4) - (3.9) are unaffected. Since we are dealing with a D2Q9 lattice, the following modification is made: (3.26)
M
1 ∂fα + ξ α · ∇ x fα − F α = − (fα − fαeq ). ∂t ετ The new set of functions Fα must be such that
ED
8 X
Fα = 0,
(3.27)
α=0
CE
PT
which guarantees that the mass conservation equation (3.22) is unchanged. Next, one requires that 8 X
Fα ξ α = ρb,
(3.28)
α=0
AC
so that the momentum equation (3.25) has the form ρa + ∇p − ∇ · τ − ρb = 0 + O(ε).
(3.29)
Thus, one choice for the set of Fα is: 1 1 ρb · ξ 1 , F3 = 2 ρb · ξ 3 , 2 2σ 2σ 1 1 F5 = 2 ρb · ξ 5 , F7 = 2 ρb · ξ 7 , 2σ 2σ Fα = 0, α = 0, 2, 4, 6, 8. F1 =
12
(3.30) (3.31) (3.32)
ACCEPTED MANUSCRIPT To verify (3.27) - (3.28), one observes that F1 = −F5 , F3 = −F7 . The Energy Equation
In order to obtain the energy equation, an internal energy distribution function gα is introduced and it is assumed to satisfy an evolution equation similar to that for fα . Thus,
CR IP T
∂gα 1 + ξ α · ∇x gα − Gα = − (gα − gαeq ). ∂t ετ
(3.33)
Here, gαeq has a linear form in the lattice velocity vector as follows:
AN US
gαeq = Dα + ξ α · Eα . And,
gα = gαeq + εgα(1) + ε2 gα(2) + O(ε3 ), with the requirement that
gα(n) = 0,
M
8 X
α=0
n ≥ 1.
(3.34)
(3.35)
(3.36)
ED
The energy equation applicable to an incompressible continuous media is given by
PT
ρ
de 1 = τ : A1 − ∇ · q + ρr, dt 2
(3.37)
CE
where e is the internal energy, A1 is the first Rivlin-Ericksen tensor [30], q is the heat efflux vector and r is the external supply. The derivation of this equation can be found in standard books on continuum mechanics and rheology; for example, see Tanner [31].
AC
In order to derive the above equation from the internal energy distribution function, it is assumed that et is the total energy given by the sum of the internal and kinetic energies, i.e., 1 et = e + |u|2 . 2
Next, the following consistency relations must hold: 13
(3.38)
ACCEPTED MANUSCRIPT 8 X
gαeq = ρet ,
(3.39)
α=0
!
8 X
1 gαeq ξ α = p + ρe + ρ|u|2 u − τ u + q, 2 α=0 8 X
α=0
(3.40)
Gα = ρb · u − ρr.
(3.41)
D0 = ρet ,
D1 = 0,
D2 = 0.
CR IP T
One way of satisfying the above is to assume, as before, that the scalars Dα are such that Dα = D1 , α = 1, 3, 5, 7, and Dα = D2 , α = 2, 4, 6, 8, and set (3.42)
In addition, it is assumed that the vectors Eα are defined through E0 = 0, Eα = E1 , α = 1, 3, 5, 7; Eα = E2 , α = 2, 4, 6, 8, where !
Finally, the scalars Gα must satisfy
1 1 ρb · (ξ α ⊗ ξ α )u − 2 ρr(ξ α ⊗ ξ α ) : 1, 2 2σ 4σ Gα = 0, α = 0, 2, 4, 6, 8.
(3.43)
α = 1, 3, 5, 7, (3.44)
M
Gα =
E2 = 0.
AN US
1 1 E1 = 2 p + ρe + ρ|u|2 u − τ u + q, 2σ 2
(3.45)
ED
Letting b = bx i + by j, it is easy to verify that b · (ξ 1 ⊗ ξ 1 )u = b · (ξ 5 ⊗ ξ 5 )u = σ 2 bx u, b · (ξ 3 ⊗ ξ 3 )u = b · (ξ 7 ⊗ ξ 7 )u = σ 2 by v, (ξ α ⊗ ξ α ) : 1 = σ 2 , α = 1, 3, 5, 7,
PT
(3.46) (3.47) (3.48)
CE
which means that (3.41) is satisfied.
AC
Finally, summing Eq. (3.33) over all α, and using (3.34) - (3.36), one finds that de ∂u ·u+∇· ρ +ρ dt ∂t
8 X
gαeq ξ α
α=0
!
− ρb · u = 0 + O(ε).
(3.49)
Next, !
!
1 ∇ · p + ρe + ρ|u|2 u = ∇p + ρ∇e + ρ(u · ∇)u · u, 2 14
(3.50)
ACCEPTED MANUSCRIPT because ∇ · u = 0. Finally, using indicial notation for convenience and noting the symmetry of the stress tensor, 1 (τij uj ),i = τij,i uj + τij ui,j = τji,j ui + τij (A1 ),ij . 2
(3.51)
Thus, 1 ∇ · (ττ u) = (∇ · τ ) · u + τ : A1 . 2 Combining all of the above, we obtain
CR IP T
!
(3.52)
de 1 ρ + ρa + ∇p − ∇ · τ − ρb · u − τ : A1 + ∇ · q − ρr = 0 + O(ε), (3.53) dt 2
Compressible Fluids
AN US
which reduces to (3.37) when one appeals to the equations of motion (3.29).
The continuity, the equations of motion and the energy equations for compressible fluids can be derived quite easily. First of all, Eqs. (3.4) - (3.8) remain unchanged, while the matrix M in Eq. (3.9) has the form:
M
2 ρu − Txx ρuv − Txy
(3.54)
ED
M= , ρuv − Txy ρv 2 − Tyy
PT
where Txx , Txy , Tyy are the components of the total stress tensor. A glance at Eqs. (3.10) - (3.13) shows that A1 , A2 , B1 , B2 are unchanged, while one needs the following: ρ|u|2 Txx + Tyy + , σ2 σ2
(3.55)
CE
A0 = ρ −
1 1 (ρu2 − Txx ), C22 = 4 (ρv 2 − Tyy ), 4 2σ 2σ 1 C12 = C21 = 4 (ρuv − Txy ). 8σ
AC
C11 =
(3.56) (3.57)
Next, the energy equation (3.37) is now given by ρ
de 1 = T : A1 − ∇ · q + ρr. dt 2
Thus, one requires the following changes be made: 15
(3.58)
ACCEPTED MANUSCRIPT 8 X
!
1 gαeq ξ α = ρe + ρ|u|2 − T u + q, 2 α=0
(3.59)
!
1 1 E1 = 2 ρe + ρ|u|2 − T u + q. 2σ 2
(3.60)
Boundary conditions
Three Dimensional Equations for Continua
PT
4
ED
M
AN US
CR IP T
One of the main advantages of the current approach is that boundary conditions can be incorporated in a manner similar to macroscopic methods, in contrast with other methods utilised for solving LBM equations. The latter employ complicated special relationships for the discrete particle distribution function (fα ) and the internal energy distribution function (gα ) for each kind of boundary conditions and problems [32,11]. For example, methods such as on-grid and mid-grid bounce back are used when the velocity is zero on the boundary; when the boundary is in motion, bounce-back is used along with a set of linear equations to determine the boundary values fα . In the method used here, the boundary conditions of fα and gα can be obtained directly from the macroscopic values on the boundaries due to the relationships of the macroscopic values with fα as seen from (3.3) - (3.7); in the case of gα , see (3.34), (3.38) - (3.41). As a result, in this method, boundary conditions, especially the Dirichlet conditions, can be included in various problems similar to macroscopic methods and no special equations for fα and gα are needed to incorporate the boundary conditions.
AC
CE
In order to derive the equations of continuum mechanics applicable to motions in three dimensions, one can employ a D3Q15 lattice [2]. Of course, any spacefilling, symmetric lattice can be used for the derivation of the equations of continuum mechanics, since the method relies on the expansion of the particle distribution function as a polynomial in the particle velocity; see (4.2) below. In the case of the D3Q15 lattice, there are 15 nodes and one lies at the centre of a cube with a side of 2 units each; six are chosen from the midpoints of the faces of the cube and the remaining eight are the vertices. In the Table, the positions of the nodes and the corresponding microscopic velocities ξ α , divided by σ for convenience, are given. In order to obtain the continuity equation and the equations of motion, we observe that Eqs. (3.2 ) and (3.3) hold and the evolution equation for the 16
ACCEPTED MANUSCRIPT particle distribution function is once again given by (cf. (3.26)):
1 ∂fα + ξ α · ∇ x fα − F α = − (fα − fαeq ), ∂t ετ
(4.1)
CR IP T
with the following quadratic form:
fαeq = Aα + ξ α · Bα + (ξ α ⊗ ξ α ) : Cα ,
Node
Location
ξ α /σ
α0
(0, 0, 0)
0
α1
(1, 0, 0)
i
α2
(0, 1, 0)
j
α3
M
AN US
Table of Nodes and Microscopic Velocities
(-1, 0, 0)
−i
(0, -1, 0)
−j
ED
α4
(0, 0, 1)
k
α6
(0, 0, -1)
−k
α7
(1, 1, 1)
i+j+k
α8
(-1, 1, 1)
−i + j + k
α9
(-1, -1, 1)
−i − j + k
α10
(1, -1, 1)
i−j+k
α11
(1, 1, -1)
i+j−k
α12
(-1 , 1, -1)
−i + j − k
α13
(-1, -1, -1) −i − j − k
α14
(1, -1, -1)
AC
CE
PT
α5
i−j−k
Next, Eqs. (3.3) - (3.9) are modified as follows: 17
(4.2)
ACCEPTED MANUSCRIPT 14 X
fαeq = ρ,
(4.3)
α=0 14 X
fαeq ξ α = ρu,
u = ui + vj + wk,
(4.4)
α=0
α=0
fαeq ξ α ⊗ ξ α = M, 14 X
fα(n) = 0,
n ≥ 1,
fα(n) ξ α = 0,
n ≥ 1.
α=0 14 X
(4.5)
α=0
(4.6)
CR IP T
14 X
(4.7)
In (4.5), M has the matrix form
2 ρu − Txx ρuv − Txy ρuw − Txz
2 . M= ρuv − Txy ρv − Tyy ρvw − Tyz
(4.8)
AN US
ρuw − Txz ρvw − Tyz ρw2 − Tzz
In the above set, ρ is the density, u, v and w are the components of the velocity field u in the x, y and z directions respectively, Txx , Txy , Txz , Tyy , Tyz , Tzz are the stresses which can be defined through any relevant constitutive relation for the continuous medium. The body force distribution functions now satisfy
M
14 X
Fα = 0,
(4.9)
Fα ξ α = ρb.
(4.10)
ED
α=0
and
14 X
α=0
F1 =
CE
F0 = 0,
PT
As before, one choice for the set of Fα is: 1 ρb · ξ 1 = −F3 , 2σ 2
F2 =
(4.12)
AC
1 ρb · ξ 2 = −F4 , 2σ 2 1 F5 = 2 ρb · ξ 5 = −F6 , 2σ Fα = 0, α = 7, · · · , 14.
(4.11)
(4.13) (4.14)
Finally, in (4.2), one has to specify the scalars Aα , the vectors Bα and the matrices Cα . We assume that A1 = Aα , α = 1, · · · 6; A2 = Aα = 7, · · · 14, and define ρ|u|2 Txx + Tyy + Tzz A0 = ρ − 2 + , A1 = A2 = 0. (4.15) σ σ2 18
ACCEPTED MANUSCRIPT Next, B0 = 0, and B1 = Bα , α = 1, · · · 6; B2 = Bα = 7, · · · 14, and B1 =
ρu , 2σ 2
B2 = 0.
(4.16)
Next, once again, C0 = 0, and C1 = Cα , α = 1, · · · 6; C2 = Cα = 7, · · · 14, and C11 C1 = 0
C22 =
0
0
, C22 0 0 0 C33
1 (ρv 2 − Tyy ), 2σ 4
C11 =
1 (ρu2 − Txx ), 4 2σ
C33 =
1 (ρw2 − Tzz ). 2σ 4
0 C12 C13
C2 = , C21 0 C23 C31 C32 0
C23 = C32 =
AN US
Finally,
C12 = C21 =
1 (ρvw − Tyz ), 16σ 4
(4.17)
CR IP T
1 (ρuv − Txy ), 16σ 4
C13 = C31 =
1 (ρxz − Txz ). 16σ 4
(4.18)
(4.19)
(4.20)
ED
M
Repeating the calculations employed in Section 3, one can now derive the continuity equation and the equations of motion for a continuous medium, valid in three dimensions. In order to obtain the energy equation, one begins with
PT
1 ∂gα + ξ α · ∇x gα − Gα = − (gα − gαeq ), ∂t ετ
(4.21)
CE
where gαeq has a polynomial expansion: gαeq = Dα + ξ α · Eα .
(4.22)
And,
AC
gα = gαeq + εgα(1) + ε2 gα(2) + O(ε3 ), with the requirement that 14 X
gα(n) = 0,
α=0
n ≥ 1.
(4.23)
(4.24)
The energy equation applicable to a continuous media is given by ρ
de 1 = T : A1 − ∇ · q + ρr. dt 2 19
(4.25)
ACCEPTED MANUSCRIPT In order to derive the above equation from the internal energy distribution function, it is assumed that et is the total energy given by the sum of the internal and kinetic energies, i.e., 1 et = e + |u|2 . 2
(4.26)
Next, the following consistency relations must hold:
gαeq = ρet ,
(4.27)
α=0 14 X
gαeq ξ α
α=0 14 X
α=0
CR IP T
14 X
!
1 = ρe + ρ|u|2 u − Tu + q, 2
(4.28)
Gα = ρb · u − ρr.
(4.29)
D0 = ρet ,
AN US
One way of satisfying the above is to assume, as before, that Dα = D1 , α = 1, · · · , 6, and Dα = D2 , α = 7, . . . 14, and set D1 = 0,
D2 = 0.
(4.30)
In addition, it is assumed that E0 = 0, Eα = E1 , α = 1, · · · , 6; Eα = E2 , α = 7, · · · , 14, where
M
!
ED
1 1 E1 = 2 ρe + ρ|u|2 u − Tu + q, 2σ 2
E2 = 0.
(4.31)
Finally, G0 = 0, and
α = 1, · · · , 6, (4.32) (4.33)
CE
PT
1 1 ρb · (ξ α ⊗ ξ α )u − 2 ρr(ξ α ⊗ ξ α ) : 1, 2 2σ 4σ Gα = 0, α = 7, · · · , 14.
Gα =
AC
Once again, repeating the calculations in Section 3, the energy equation (4.25) can be derived.
5
Concluding Remarks
Using an evolution equation for the modified particle distribution functions, equations of continuum mechanics in two and three dimensions have been derived. It is obvious that the choice of the constitutive equation for the incompressible or compressible fluid is unrestricted. Focussing attention on 20
ACCEPTED MANUSCRIPT
CR IP T
two-dimensional flows for convenience, it can be seen from the definition of the matrix M in (3.9) that the extra stresses τxx , τxy , τyy , in an incompressible fluid, can be chosen at will. That is, the fluid model can be Newtonian or non-Newtonian. Similar observations apply to the matrix M in (4.8). The derivations given here, based on that in [1,2], are applicable to all fluids, compressible or incompressible. As far as incompressible fluids are concerned, it has to be emphasised again that in deriving the continuity equation (3.22) and the momentum equations (3.25) or (3.29), the density ρ and the pressure p are independent of each other. Secondly, the equations of continuum mechanics appear at the expansion of the particle distribution functions at the first order, as described in full in (3.3) - (3.53) and (4.1) - (4.33).
AN US
Next, some remarks regarding the numerical scheme are in order. As is well known, there are two main categories of flows in fluid mechanics. One is where the pressure gradient is specified, which occurs in the case of the flow in a pipe. The second is where the pressure field has to be determined as a part of the solution; consider the flow in a lid-driven cavity as an example. Both types of problems have to be solvable with the chosen numerical scheme for the determination of the distribution functions fα and gα .
M
For instance, consider the evolution equation for fα in (3.26). This is solved by using the splitting method of Toro [33] which is explained in detail in [1,2]. Briefly, the evolution equation is split into two parts. The first one is called the streaming section, equivalent to solving the homogeneous equation:
ED
∂fα + ξ α · ∇x fα − Fα = 0. ∂t
(5.1)
AC
CE
PT
In order to solve this equation, an initial value for fαeq is required. To determine this, one begins with an assumed velocity distribution in the domain, say u0 , from which one determines the initial stress field through the designated constitutive equation. If the pressure field is prescribed, this can be used to calculate the matrix M from (3.9) and the initial value of fαeq from (3.3). If the pressure has to be determined as a part of the solution, one assumes that, in addition to u0 , a value of fαeq is given. From (3.9), the pressure p can be determined from the trace of the matrix M in (3.9). Next, using the initial value of fαeq , determined from the given initial macroscopic variables as in (3.3) - (3.9), the intermediate value of fαI is calculated by solving the homogeneous equation; for instance, the Lax-Wendroff scheme can be employed. Using this, the intermediate macroscopic value of the velocity field uI is determined from (3.5). The second one is called the collision section, equivalent to solving the time21
ACCEPTED MANUSCRIPT dependent equation: ∂fα 1 = − (fα (x, t) − fαeq (x, t)). ∂t ετ
(5.2)
Using the previously determined fαI as the initial condition, the collision equation is now solved; for instance, the Euler method can be used.
CR IP T
To reiterate, considering two-dimensional flows for convenience, we note that the matrix M in (3.9) can be used in two different ways. If the pressure p is specified and the initial velocity vector u0 is known, one uses the relevant constitutive equation in this matrix to find the initial value of fαeq . On the other hand, if the pressure field has to be determined, the trace of the matrix can be employed; see Eqs. (18a) - (18c) in [1]. Clearly, the constant density assumption is met, for it is inherent in the Eqs. (3.3) - (3.13). Moreover, at each iteration, a new value of the lattice speed σ is chosen employing (5.8) or (5.9) below.
AN US
Next, we shall discuss the stability of the numerical scheme. Multiplying fαeq with |ξξ α |2 /2 and the utilisation of (3.15) - (3.17) leads to 8 X
1 eq τxx + τyy 1 fα |ξξ α |2 = p + ρ |u|2 − . 2 2 α=0 2
(5.3)
M
Next, it is easy to verify from (3.30) - (3.32) that 8 X
ED
α=0
Fα |ξξ α |2 = 0.
(5.4)
(5.5)
PT
Hence, Eq. (5.1) becomes
CE
∂ 1 τxx + τyy σ2 p + ρ|u|2 − + ρ(∇ · u) = O (ε) . ∂t 2 2 2
AC
The Courant-Friedrichs-Lewy (CFL) condition states that ([34] and [35])
K=
u∆t v∆t + ≤ 1. ∆x ∆y
(5.6)
This can be used in (5.5) and we obtain "
#
2p − τxx − τyy |u| + + σ 2 K = O (ε) . ρ 2
22
(5.7)
ACCEPTED MANUSCRIPT Thus, the lattice speed σ must satisfy
σ=
v u u τxx Kc t
+ τyy − 2p − |u|2 , ρ
1 Kc = √ ≥ 1. K
(5.8)
If the pressure p has to be uniquely defined, one requires that τxx + τyy = 0; see (3.14). Thus, in these fluid models, (5.8) reduces to v u u −2p Kc t ρ
− |u|
,
CR IP T
1 Kc = √ ≥ 1 (5.9) K As a result, the value σ is modified and changes in each iteration as defined through (5.9). σ=
2
ED
M
AN US
Finally, in the Finite Difference Lattice Boltzmann Method (FDLBM) adopted here, the iteration and recovery of the pressure field is similar to the SIMPLE method of Patankar and Spalding [36,37]. As is well known, the SIMPLE method is a guess-and-correct procedure for the calculation of the pressure field. In each iteration, the velocity field is obtained from the first guessed pressure field. Next, using the corrected velocity field, it is possible to find the corrected pressure and this process continues till a very small or zero mass residual is obtained, since the zero mass residual demonstrates that the divergence of the velocity vector field is zero. In FDLBM, the criteria, which is the mass residual in the SIMPLE method, is the difference between the sum of the distribution functions and the fixed density; see Eq. (3.4). Thus, the correction of the pressure field and the subsequent correction of the velocity field continues till a small or zero difference exists between this sum and the density.
CE
PT
For an in-depth discussion of these matters, see [1,2]. For an application to the lid-driven, mixed convection flow of a Bingham fluid in a cavity, see [38].
AC
Acknowledgement: The research of G. H. R. Kefayati is supported by an International Postgraduate Research Scholarship from the Flinders University of South Australia. We would also like to thank the three reviewers for their extremely helpful comments which have been useful in revising the paper.
23
ACCEPTED MANUSCRIPT References [1] S. C. Fu, R. M. C. So, Modeled Lattice Boltzmann Equation and the ConstantDensity Assumption. AIAA Journal, 47 (2009), 3038–3042. [2] S. C. Fu, R. M. C. So, R. M. C. Leung, Linearized-Boltzmann-type-equationbased finite difference method for thermal incompressible flow. Computers & Fluids, 6 (2012), 67–80.
CR IP T
[3] Z. Guo, B. Shi, N Wang, Lattice BGK model for incompressible Navier-Stokes equation. J. Comp. Phys., 165 (2000), 288–306.
[4] U. Frisch, B. Hasslacher, Y. Pomeau, Lattice gas automata for the Navier-Stokes equation. Phys. Rev. Lett. 56 (1986), 1505–1508.
[5] U. Frisch, D. d’Hum`eres, B. Hasslacher, P. Lallemand, Y. Pomeau, J-P. Rivet, Lattice gas hydrodynamics in two and three dimensions. Complex Syst. 1(1987), 649–707.
AN US
[6] S. Wolfram, Cellular automaton fluids 1: Basic theory. J. Stat. Phys. 45 (1986), 471–526 .
[7] X. He, L-S. Luo, A priori derivation of the lattice Boltzmann equation. Phys. Rev. E, 55 (1997), R6333–R6336.
M
[8] P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. Ser. 2, 94 (1954), 511–525.
ED
[9] P. Welander, On the temperature jump in a rarefied gas. Ark. Fysik, 7 (1954), 507–553. [10] P. J. Davis, P. Rabinowitz, Methods of Numerical Integration. Academic Press, New York, 1984.
PT
[11] Z. Guo, B. Shi, C. Zheng, A coupled lattice BGK model for the Boussinesq equations. Int. J. Numer. Methods Fluids, 39 (2002), 325–342.
CE
[12] M. Junk, W-A Yong, Rigorous Navier-Stokes limit of the lattice Boltzmann equation. Asymptotic Analysis, 35 (2003), 165–186.
AC
[13] M. Junk, A. Klar, L-S. Luo, Asymptotic analysis of the lattice Boltzmann equation. J. Comp. Phys. 210 (2005), 676–704. [14] F. Dubois, Equivalent partial differential equations of a lattice Boltzmann scheme. Comp. Math. Appl. 55 (2008), 1441–1449. [15] F. Dubois, Third order equivalent equation of lattice Boltzmann scheme. Discrete Cont. Dyn. Syst. 23 (2009), 221–248. [16] F. Dubois, P. Lallemand, Towards higher order lattice Boltzmann scheme. J. Stat. Mech. doi:10.1088/1742-5468/2009/06/P06006.
24
ACCEPTED MANUSCRIPT [17] P. Lallemand, L-S. Luo, Theory of the lattice Boltzmann method: Acoustic and thermal properties in two and three dimensions. Phys. Rev E, 68 (2003), 036706. [18] S. Gabbanelli, G. Drazer, J. Koplik, Lattice Boltzmann method for nonNewtonian (powerlaw) fluids, Phys. Rev E 72 (2005) 046312. [19] J. Boyd, J. Buick, S. Green, A second-order accurate lattice Boltzmann nonNewtonian flow model, Journal of Physics A: Mathematical and General, 39 (2006) 14241–14247.
CR IP T
[20] S.P. Sullivan, L.F. Gladden, M.L. Johns, Simulation of powerlaw fluid flow through porous media using lattice Boltzmann techniques, J. Non-Newt. Fluid Mech. 133 (2006) 91–98.
[21] S.P. Sullivan, A.J. Sederman, M.L. Johns, L.F. Gladden, Verification of shear thinning LB simulations in complex geometries, J. Non-Newt. Fluid Mech. 143 (2007) 59–63.
AN US
[22] J. Psihogios, M.E. Kainourgiakis, A.G. Yiotis, A.Th. Papaioannou, A.K. Stubos, A lattice Boltzmann study of non-Newtonian flow in digitally reconstructed porous domains, Trans. Porous Media 70 (2007) 279–292.
[23] M. Yoshino, Y. Hotta, T. Hirozane, M. Endo, A numerical method for incompressible non-Newtonian fluid flows based on the lattice Boltzmann method, J. Non-Newt. Fluid Mech., 147 (2007) 69–78.
M
[24] J.M. Buick, Lattice Boltzmann simulation of power-law fluid flow in the mixing section of a single-screw extruder, Chem. Eng. Sci., 64 (2009) 52–58.
ED
[25] Amir Nejat, Vahid Abdollahi, Koohyar Vahidkhah, Lattice Boltzmann simulation of non Newtonian flows past confined cylinders, J. Non-Newt. Fluid Mech. 166 (2011) 689–697.
PT
[26] C.H. Wang, J.R. Ho, Lattice Boltzmann modeling of Bingham plastics, Physica A 387 (2008) 4740–. [27] K. R. Rajagopal, A. R. Srinivasa, On the nature of constraints for continua undergoing dissipative processes. Proc. R. Soc. A461 (2005), 2785–2795.
CE
[28] R. R. Huilgol, On the definition of pressure in rheology. Rheol. Bull., 78(2) (2009), 12, 14–15, 29.
AC
[29] R. R. Huilgol, Fluid Mechanics of Viscoplasticity. Springer, Berlin Heidelberg, 2015. [30] R. S. Rivlin, J. L. Ericksen, Stress-deformation relations for isotropic materials. J. Rational Mech. Anal. 4 (1955), 323–425. [31] R. I. Tanner, Engineering Rheology, 2nd Ed. Oxford University Press, Oxford, 2000. [32] Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Phys. Fluids, 9 (1997), 1591–1598.
25
ACCEPTED MANUSCRIPT [33] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics: a practical introduction, 2nd Ed. Springer-Verlag, New York, 1999. [34] J. Blazek, Computational Fluid Dynamics: Principles and Applications. Elsevier, Amsterdam, 2001. See pp. 347–350. [35] T. Cebeci, J. P. Shao, F. Kafyeke, E. Laurendeau, Computational Fluid Dynamics for Engineers. Springer-Verlag, New York, 2005. See pp. 311–320.
CR IP T
[36] S. V. Patankar, D. B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int. J. Heat Mass Transfer 15 (1972) 1787–1806. [37] S. V. Patankar, A calculation procedure for two-dimensional elliptic situations. Numerical Heat Transfer 4 (1981) 409–425.
AC
CE
PT
ED
M
AN US
[38] G. H. R. Kefayati, R. R. Huilgol, Lattice Boltzmann Method for the simulation of mixed convection flow of a Bingham fluid in a lid-driven cavity. J. Non-Newt. Fluid Mech. (submitted).
26
From Mesoscopic Models to Continuum Mechanics: Newtonian and Non-Newtonian Fluids R.R. Huilgol, G.H.R. Kefayati PII: DOI: Reference:
S0377-0257(16)30003-9 10.1016/j.jnnfm.2016.03.002 JNNFM 3762
To appear in:
Journal of Non-Newtonian Fluid Mechanics
Received date: Revised date: Accepted date:
1 September 2015 2 March 2016 3 March 2016
Please cite this article as: R.R. Huilgol, G.H.R. Kefayati, From Mesoscopic Models to Continuum Mechanics: Newtonian and Non-Newtonian Fluids, Journal of Non-Newtonian Fluid Mechanics (2016), doi: 10.1016/j.jnnfm.2016.03.002
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT Highlights
AC
CE
PT
ED
M
AN US
CR IP T
• Critical examination of the BGK approximation to derive constitutive equations for Newtonian fluids • Limitations due to grid dependence of viscosity in Newtonian and nonNewtonian fluids • A new evolution equation for the lattice distribution function • Derivation of the equations of continuum mechanics for Newtonian and nonNewtonian fluids from the new model
1
ACCEPTED MANUSCRIPT
From Mesoscopic Models to Continuum Mechanics: Newtonian and Non-Newtonian Fluids R. R. Huilgol ∗ and G. H. R. Kefayati
CR IP T
School of Computer Science, Engineering and Mathematics, Flinders University of South Australia, GPO Box 2100, Adelaide, SA 5001, Australia
Abstract
PT
ED
M
AN US
A review of the BGK approximation to obtain the equations of motion for an incompressible fluid is presented and its drawbacks are revealed. In order to overcome these inherent problems, new models for the particle distribution functions are needed. Using the Finite Difference Lattice Boltzmann Method (FDLBM) due to Fu and So [1] and the Thermal Difference Discrete Flux Method (TDDFM) proposed by Fu et al [2], it is shown that the newer distribution functions lead to the mass conservation equation, the equations of motion and the energy balance equation for incompressible fluids in two dimensions, employing the D2Q9 lattice as the model. This derivation is extended to compressible fluids as well. Next, using the D3Q15 lattice as an example, the three dimensional equations of continuum mechanics are derived. Since no restrictions are placed on the constitutive equations, the theoretical development applies to all fluids, whether they be Newtonian, or power law fluids, or viscoelastic and viscoplastic fluids. Finally, some comments are offered regarding the numerical scheme to calculate the particle distribution functions to determine the velocity and temperature fields.
AC
CE
Key words: Lattice Boltzmann equation, BGK approximation, particle distribution function, continuum mechanics
1
Introduction
From the Lattice Boltzmann equation, it is possible to derive the continuity equation and Cauchy’s equations of motion for a compressible medium, when ∗ Corresponding author. Fax: +61 8 8201 2904. Email address: [email protected] (R. R. Huilgol).
Preprint submitted to Elsevier
22 March 2016
ACCEPTED MANUSCRIPT one uses the Bhatnagar-Gross-Krook (BGK) approximation. From this, one can obtain equations relevant to incompressible fluids. However, these require that the pressure be proportional to the density and the viscosity be dependent on the collision relaxation time [3]; see Section 2 below. Clearly, these restrictions on the pressure and the viscosity are unacceptable in modelling the flows of non-Newtonian, incompressible fluids.
M
AN US
CR IP T
In order to overcome these inherent problems, new models for the particle distribution functions are needed. In the Finite Difference Lattice Boltzmann Method (FDLBM) due to Fu and So [1], the particle distribution function leads to the conservation of mass and the equations of motion applicable to incompressible fluids, when the flows are assumed to occur in a two dimensional setting underpinned by a D2Q9 lattice. Our derivation of these results is succinct and is more transparent, for it uses using vector analysis and linear algebra. In addition, the energy equation is also obtained from Thermal Difference Discrete Flux Method (TDDFM) proposed by Fu et al [2]; once again, simple results from vector analysis and linear algebra are employed. The important point to note is that the previous restrictions on the pressure and the viscosity are eliminated in these derivations, meaning that one is free to choose a constitutive equation. That is, one can model a Newtonian fluid, or power law fluids, or viscoelastic and viscoplastic fluids. Moreover, we point out in sub-section 3.1 that it is quite easy to incorporate Dirichlet type boundary conditions into the numerical scheme to determine the particle distribution functions for the velocity and temperature fields.
ED
In Section 4, the method is extended to the flows of all fluids in three dimensions, using the D3Q15 lattice as an example. Once again, the particle distribution functions are such that every type of fluid, compressible or incompressible, can be employed.
BGK Approximation to Continuum Mechanics
AC
2
CE
PT
Finally, in Section 5, some remarks are offered on the numerical scheme employed to determine the particle distribution functions for the velocity and the temperature fields. The CFL condition for the stability of the numerical scheme is also derived.
Beginning in 1986, Lattice Boltzmann equation (LBE) models evolved from their Boolean counterparts, viz., the lattice-gas-automata. The theoretical framework of the LBE was underpinned by the Chapman-Enskog analysis of the LGA models (Frisch et al [4,5]; Wolfram [6]). That is, the statistical mechanics of the LGA played a crucial role in these developments. A decade later, He and Luo [7] proved that the LBE is a specific discretised form of 3
ACCEPTED MANUSCRIPT the continuous Boltzmann equation using the Bhatnagar-Gross-Krook (BGK) approximation [8]. As an application, He and Luo [7] derived the Lattice Boltzmann BGK (LBGK) equation for the D2Q9 lattice model. It is worth noting that Welander [9] published an approximation to the Boltzmann equation at the same time as BGK and applied it to a rarefied monatomic gas exchanging heat with an adjacent wall. In several aspects, his treatment of this problem is exemplary and anticipates by several decades the subsequent Lattice formulation. However, we shall use the more commonly used acronym LBGK equation in the sequel.
CR IP T
From the LBGK equation, Cauchy’s equations of motion for compressible fluids can be derived under the assumption that the Mach number is small and that the density variation is small. Attempts to derive the equations for incompressible fluids from the foregoing end in an physically unacceptable requirements, viz., that the pressure be proportional to the density and that the viscosity depends on the collision relaxation time [3].
AN US
To render these remarks transparent, let f = f (x, ξ, t) be the probability of finding a particle with the velocity ξ near the point x at time t. The Boltzmann equation for this distribution function, in the absence of external forces, is given by
ED
M
1 ∂f + ∇f · ξ + f = N (f ), (2.1) ∂t λ where ξ is the microscopic velocity, λ is the relaxation time due to collision and N (f ) is a measure of the net number of molecules which disappear from the region due to inter-molecular collisions; see Welander [9], for example. This is not applied to the problem at hand, for it depends upon an approximation to this measure of collisions; see (2.3) below.
PT
The macroscopic variables are the density of mass ρ, the velocity field u and the absolute temperature T. These are related to the distribution function f through the following integrals in momentum space:
CE
Z
ρ=
f dξ,
ρu =
Z
ξf dξ,
1Z ρε = |ξ − u|2 f dξ, 2
(2.2)
AC
where ε = D0 RT /2, with D0 being the number of degrees of freedom of a particle and R is the ideal gas constant. Note that we have used dξ = dξ1 · · · dξD , where D ∈ [2, 3] stands for the dimension of the physical space. In the BGK approximation, the collision operator is assumed to be given by N (f ) = g(x, ξ, t)/λ, where the collision function g(x, ξ, t) has the form: "
#
ρ |ξ − u|2 g(x, ξ, t) = exp − . (2πRT )D/2 2RT 4
(2.3)
ACCEPTED MANUSCRIPT Next, using (2.3) and integrating (2.1) formally over a time interval 4t with et /λ as the integrating factor, one obtains [7]: f (x + ξ4t, ξ, t + 4t) = e−4t/λ f (x, ξ, t) Z 4t 1 0 et /λ g(x + ξt0 , ξ, t + t0 ) dt0 . + e−4t/λ λ 0
(2.4)
Assuming that the time interval 4t is small enough and linearising in terms of this time interval, one has:
CR IP T
1 f (x + ξ4t, ξ, t + 4t) − f (x, ξ, t) = − [f (x, ξ, t) − g(x, ξ, t)], τ
(2.5)
where τ = λ/4t is the non-dimensional relaxation time. Next, suppose that the collision function g(x, ξ, t) can be expanded as a Taylor series in u, retaining terms up to order |u|2 . Identifying this as the equilibrium distribution function f (0) , the following can be derived from (2.3): "
AN US
f
(0)
#
|u|2 ρ ξ · u (ξ · u)2 2 exp(−|ξ| /2RT ) + − = 1 + . (2πRT )D/2 RT 2(RT )2 2RT
(2.6)
To derive Cauchy’s equations of motion for a continuous medium, exact evaluation of the following momentum integrals m (0)
|ξ| f
dξ =
Z
exp(−|ξ|2 /2RT )ψ(ξ) dξ,
M
I=
Z
m = 0, 1, 2,
(2.7)
PT
ED
is required. Using a Cartesian coordinate system and noting that the physical dimension of the space D = 2, one can express the microscopic velocity as ξ = (ξx , ξy ). Using this, ψ(ξ) is expressed as a polynomial of the form ψ(ξ) = ξxr ξys .
CE
The integrals in (2.7) can be evaluated through a Gaussian-type quadrature [10] and lead to X I= Wα exp(−|ξ α |2 /2RT )ψ(|ξ α |), (2.8) α
AC
where Wα are the weights and ξ α are the discrete velocities of the quadrature. In the D2Q9 model, there are nine discrete velocities: {ξ α , α = 0, 1, · · · , 8}. The configuration space is discretised into a square lattice with a lattice con√ stant 4x = 3RT 4t. In fact, in an isothermal problem, the temperature T has no physical significance √ and thus, one can choose 4x as a fundamental quantity instead, leading to 3RT = c = 4x/4t, where the lattice speed c is related to the speed of sound cs through c2 = 3c2s . Thus, RT = c2s = c2 /3. The nine lattice points are located at the centre of the lattice, at the four midpoints of the edges of the square and at the four vertices. One identifies 5
ACCEPTED MANUSCRIPT α0 = (0, 0) as the centre of the lattice, {α1 , α3 , α5 , α7 } as the midpoints of the four sides: (±1, 0), (0, ±1), and {α2 , α4 , α6 α8 } = (±1, ±1) as the vertices. The lattice velocity vector is given by ξ 0 = 0. Next, the velocities for α = 1, 3, 5, 7, have the form c(±1, 0), c(0, ±1), and for α = 2, 4, 6, 8, they have the form c(±1, ±1). In other words, except at the centre, at each lattice point αj , j = 1, · · · , 8, the velocity vector points away from the lattice in the direction from α0 to αj .
where
wα =
4/9,
α = 0,
1/9, α = 1, 3, 5, 7,
1/36, α = 2, 4, 6, 8.
CR IP T
The weights Wα in (2.8) can now be computed through quadrature and one obtains [7]: Wα = 2πRT exp(|ξ α |2 /2RT )wα , (2.9)
(2.10)
AN US
The equilibrium distribution function for the D2Q9 model is given by fα(0) = ρwα + ρsα (u(x, t)), where
with
(0, 0),
(2.12)
α = 0,
c(cos θα , sin θα ) α = 1, 3, 5, 7, √ c 2(cos θα , sin θα ), α = 2, 4, 6, 8.
(2.13)
ED
ξα =
#
3(ξ α · u) 9(ξ α · u)2 3|u|2 + − , c2 c4 2c2
M
sα (u) = wα
"
(2.11)
PT
Here, the angles θα are defined through θα = (α − 1)π/4, α = 1, · · · , 8.
CE
Introducing a Mach number M = |u|/cs , one can see that the equilibrium distribution function can be written in the form: fα(0) = fα0 + Mfα1 + M2 fα2 ,
(2.14)
AC
where
fα0 = ρwα ,
fα1 = ρwα .
(2.15)
Assuming that the distribution function fα can be found from the equation corresponding to (2.5), the momentum space integrals in (2.2) lead to the following results for the macroscopic variables: ρ=
X α
fα ,
ρu =
X
fα ξ α ,
α
6
ρε =
1X fα |ξ α − u|2 . 2 α
(2.16)
ACCEPTED MANUSCRIPT When multi-scaling using the Mach number is employed, these equations lead to Cauchy’s equations for a compressible continuous medium of the following form: ∂ρ + ∇ · (ρu) = 0, ∂t
(2.17)
∂(ρu) + ∇ · (ρu ⊗ u) = −∇p + ν[∇2 (ρu) + ∇(∇ · (ρu))]. ∂t
(2.18)
AN US
CR IP T
While the first equation delivers the correct form of the equation for the conservation of mass, the second leads to the equations of motion for an artificial compressible fluid because the pressure p = ρc2s has a thermodynamic form only, with cs being the speed of sound. Of course, this form of the pressure is not surprising due to the absence of any external forces in the original formulation (2.1). However, the expression for the kinematic viscosity is unusual for it is given by ν = (2τ − 1)c2 4t/6, and depends on the non-dimensional relaxation time τ and the time step 4t.
M
Turning to incompressible fluids, we note that the pressure p must be independent of the density. To derive such a model, Guo et al [3] found that a new type of distribution function gα (x, t) is required. The important point is that this new distribution function is patterned along the lines of the BGK model; however, there is no statistical mechanical basis for this formulation. To understand this, one notes that the new equilibrium distribution function gα(0) (x, t)is defined through α = 0,
λ(p/c2 ) + sα (u)
α = 1, 3, 5, 7,
γ(p/c2 ) + sα (u),
α = 2, 4, 6, 8,
−4σ(p/c2 ) + sα (u),
ED
gα(0) =
(2.19)
PT
where σ, λ and γ are parameters satisfying
CE
λ + γ = σ,
λ + 2γ = 1/2.
(2.20)
This new distribution function satisfies the following conservation laws:
AC
8 X
gα =
α=0
8 X
α=0
ceα gα =
8 X
α=0 8 X
gα(0) ,
(2.21)
ceα gα(0) .
(2.22)
α=0
The evolution equation for the system is given by 1 gα (x + ceα 4t, t + 4t) − gα (x, t) = − [gα (x, t) − gα(0) (x, t)]. τ 7
(2.23)
ACCEPTED MANUSCRIPT The velocity and pressure of the flow are given by
u=
8 X
ceα gα ,
(2.24)
α=0 " 2
#
8 c X gα + s0 (u) . p= 4σ α=1
(2.25)
∇ · u = 0,
CR IP T
From the chosen distribution function, a multi-scale expansion as described in full in the Appendix of the paper by Guo et al [3], leads to the following set of equations applicable to an incompressible fluid: (2.26)
∂u + ∇ · (u ⊗ u) = −∇p + ν∇2 u. ∂t
(2.27)
ν=
AN US
While p can now be arbitrary, the kinematic viscosity ν is still relaxation time and grid-dependent, for it is given by (2τ − 1) (4x)2 · . 6 4t
(2.28)
CE
PT
ED
M
In this context, it is necessary to recall that in addition to the work cited earlier [4,5], attempts to derive the Navier-Stokes equations for compressible or incompressible fluids have been made [11–16]. In addition, thermal stability of the hydrodynamic equations has also been investigated [17]. It is important to note that hydrodynamic equations deal, at most, with a restricted class of fluids, viz., the ideal gas or fluids with constant viscosity. Secondly, to obtain the equations for a Newtonian fluid, second and higher order lattice Boltzmann schemes have to employed [15,16]. The derivations given here in Section 3 are based on the results in [1,2] and are applicable to all fluids, compressible or incompressible. Secondly, the equations of continuum mechanics appear at the expansion of the particle distribution functions at the first order, as described in full in (3.3) - (3.50) below.
AC
Turing our attention to non-Newtonian fluids, we note that the relationship between the relaxation time and the viscosity means that the relaxation time has to be prescribed over a specific range only; hence, the viscosity function is restricted and cannot describe a wide range of non-Newtonian fluid models. This is clear from the previous studies employing LBM for non-Newtonian fluids; for example, see the work on power-law fluids [18–25] and Bingham fluids [26]. Even though a second order LBM scheme has been proposed in [19] to solve problems in power-law fluids, the kinematic viscosity is still related to the relaxation time. Hence, a different approach to deriving the equations 8
ACCEPTED MANUSCRIPT of continuum mechanics applicable to all fluids is needed, and we turn to this next.
3
Thermal Difference Discrete Flux Method to Continuum Mechanics
CR IP T
First of all, the lattice counterpart of (2.1) with the assumption that the collision term N (f ) = f eq /λ can be written as follows: 1 ∂fα + ξ α · ∇ x fα = − (fα − fαeq ), ∂t ετ
(3.1)
where ε is a small parameter to be prescribed when numerical simulations are considered. To proceed, one assumes that fα has the following expansion:
AN US
fα = fαeq + εfα(1) + ε2 fα(2) + O(ε3 ).
(3.2)
The novelty of the approach by Fu and So [1] lies in expanding the equilibrium lattice function fαeq as a quadratic in the particle velocity ξ α : fαeq = Aα + ξ α · Bα + (ξ α ⊗ ξ α ) : Cα ,
(3.3)
8 X
fαeq = ρ,
α=0
fαeq ξ α = ρu,
(3.4)
u = ui + vj,
PT
8 X
ED
M
where Bα is a vector and Cα is a 2 × 2 symmetric matrix. Before defining these entities, the following relations must hold (cf. (2.16)):
(3.5)
α=0
fαeq ξ α ⊗ ξ α = M,
CE
8 X
α=0
8 X
fα(n) = 0,
n ≥ 1,
(3.7)
fα(n) ξ α = 0,
n ≥ 1.
(3.8)
AC
α=0
8 X
α=0
(3.6)
In (3.6), M has the matrix form
2 ρu + p − τxx
M=
ρuv − τxy
9
ρuv − τxy . ρv 2 + p − τyy
(3.9)
ACCEPTED MANUSCRIPT In the above set, ρ is the density, u and v are the components of the velocity field u in the x and y directions respectively, τxx , τxy = τyx , τyy are the stresses which can be defined through any relevant constitutive relation. For a D2Q9 lattice, in (3.3), there are nine coefficients: Aα , Bα , and Cα . However, when α = 0, only the coefficient A0 is required. Regarding the others, one makes the assumption that the coefficients having the same energy shell of the lattice velocities are equal. Thus, Aα = A1 , α = 1, 3, 5, 7; Aα = A2 , α = 2, 4, 6, 8.
AN US
CR IP T
Since B0 and C0 do not affect the value of f0eq in (3.3), they can be put equal to zero. By the just mentioned assumption about certain coefficients being equal, we see that Bα = B1 , α = 1, 3, 5, 7; Bα = B2 , α = 2, 4, 6, 8. Similarly, Cα = C1 , α = 1, 3, 5, 7; Cα = C2 , α = 2, 4, 6, 8. Keeping in mind that the matrices are symmetric, there are thirteen independent quantities: three scalars in A0 , A1 , A2 ; four components of the vectors B1 , B2 ; three components each in C1 , C2 . Looking at (3.4) - (3.6), it is obvious that the number of available constraints is six only. Thus, there is a lot of latitude in choosing the variables. Following Fu and So [1], we shall assume that A0 = ρ −
2p ρ|u|2 τxx + τyy − 2 + , σ2 σ σ2 ρu , 2σ 2
M
and
B1 =
and
C11 =
PT
C11 0 C1 = , 0 C22
1 (p + ρu2 − τxx ), 2σ 4
(3.10)
(3.11)
C22 =
1 (p + ρv 2 − τyy ), 2σ 4 (3.12)
0 C12 C2 = , C21 0
CE
B2 = 0.
ED
Next,
A1 = A2 = 0,
C12 = C21 =
1 (ρuv − τxy ). 8σ 4
(3.13)
AC
Here, σ is a constant to be determined keeping in mind that the vectors ξ α have been redefined so that in (2.13), c = σ, the lattice speed. Since the value of σ affects numerical stability, its choice depends on the problem at hand. For instance, in their solution to the Stokes Second problem, Fu and So [1] recommend varying it at every iteration step. This matter is briefly addressed in Section 5 below. Before deriving the macroscopic equations of continuum mechanics for incompressible fluids, we note that in several fluids, such as viscoplastic fluids, the 10
ACCEPTED MANUSCRIPT pressure p has to be defined uniquely. This requires that the trace of the extra stress tensor should be zero. For a proof, see [27], [28] and [29]. Thus, in (3.10), τxx + τyy = 0, and, in these models, the coefficient A0 has the simple form: ρ|u|2 2p A0 = ρ − 2 − 2 . σ σ
(3.14)
CR IP T
We shall now demonstrate that the macroscopic equations for an incompressible continuous medium can be derived from (3.1) - (3.9). To make the calculations more transparent, it is helpful to note that
AN US
ρ ρ u + σ 2 C11 , f5eq = − u + σ 2 C11 , f0eq = A0 , f1eq = 2σ 2σ ρ ρ eq eq 2 2 f3 = v + σ C22 , f7 = − v + σ C22 , 2σ 2σ f2eq = f6eq = C12 , f4eq = f8eq = −C12 .
(3.15)
(3.16)
(3.17)
Hence, Eqs. (3.4) - (3.9) follow from the substitution of (3.10) - (3.13) into (3.4) - (3.6), which is as expected. Secondly, substituting the expression for fα in (3.2) into (3.1), we find that (3.18)
M
∂fαeq 1 + ξ α · ∇x fαeq = − fα(1) + O(ε). ∂t τ
ED
Since the velocity vectors ξ α do not depend on the spatial coordinates, one sees that ∇ · (fαeq ξ α ) = ξ α · ∇x fαeq . (3.19) Hence, (3.18) can be rewritten as
PT
∂fαeq 1 + ∇ · (fαeq ξ α ) = − fα(1) + O(ε). ∂t τ
(3.20)
CE
Summing the equation above, we obtain
AC
∂ ∂t
8 X
α=0
fαeq
!
+∇·
8 X
α=0
fαeq ξ α
!
8 1X =− fα(1) + O(ε). τ α=0
(3.21)
Appealing to (3.4), (3.5), (3.7) and recalling that ρ is a constant, the above equation reduces to ρ(∇ · u) = 0 + O(ε), (3.22) or the continuity equation for an incompressible medium is satisfied. Next,
∇ · (fαeq ξ α ⊗ ξ α ) = (ξ α · ∇x fαeq )ξ α . 11
(3.23)
ACCEPTED MANUSCRIPT Hence, multiplying (3.18) throughout by ξ α , summing over α, and appealing to (3.5) - (3.8), one obtains ρ
∂u + ∇ · M = 0 + O(ε). ∂t
(3.24)
Rearranging, one obtains the equations of motion for an incompressible fluid: ρa + ∇p − ∇ · τ = 0 + O(ε),
a=
∂u + (u · ∇)u. ∂t
(3.25)
AN US
CR IP T
Thus, the set of equations (3.1) - (3.13) deliver Cauchy’s equations of motion for an incompressible medium. It is also clear that there is no relationship between the density ρ and the pressure p, which is as it should be in an incompressible fluid. Moreover, one can specify the material through a chosen constitutive equation for the extra stress tensor τ , whether the fluid be Newtonian, or non-Newtonian. At this juncture, it has to be noted that in the equations of motion (3.25), the body force term is missing. It is possible to do so by altering the evolution equation, suggested in the work by Fu et al [2]. If the body force is given by ρb, it has to be incorporated into the evolution equation (3.1) in such a way that (3.4) - (3.9) are unaffected. Since we are dealing with a D2Q9 lattice, the following modification is made: (3.26)
M
1 ∂fα + ξ α · ∇ x fα − F α = − (fα − fαeq ). ∂t ετ The new set of functions Fα must be such that
ED
8 X
Fα = 0,
(3.27)
α=0
CE
PT
which guarantees that the mass conservation equation (3.22) is unchanged. Next, one requires that 8 X
Fα ξ α = ρb,
(3.28)
α=0
AC
so that the momentum equation (3.25) has the form ρa + ∇p − ∇ · τ − ρb = 0 + O(ε).
(3.29)
Thus, one choice for the set of Fα is: 1 1 ρb · ξ 1 , F3 = 2 ρb · ξ 3 , 2 2σ 2σ 1 1 F5 = 2 ρb · ξ 5 , F7 = 2 ρb · ξ 7 , 2σ 2σ Fα = 0, α = 0, 2, 4, 6, 8. F1 =
12
(3.30) (3.31) (3.32)
ACCEPTED MANUSCRIPT To verify (3.27) - (3.28), one observes that F1 = −F5 , F3 = −F7 . The Energy Equation
In order to obtain the energy equation, an internal energy distribution function gα is introduced and it is assumed to satisfy an evolution equation similar to that for fα . Thus,
CR IP T
∂gα 1 + ξ α · ∇x gα − Gα = − (gα − gαeq ). ∂t ετ
(3.33)
Here, gαeq has a linear form in the lattice velocity vector as follows:
AN US
gαeq = Dα + ξ α · Eα . And,
gα = gαeq + εgα(1) + ε2 gα(2) + O(ε3 ), with the requirement that
gα(n) = 0,
M
8 X
α=0
n ≥ 1.
(3.34)
(3.35)
(3.36)
ED
The energy equation applicable to an incompressible continuous media is given by
PT
ρ
de 1 = τ : A1 − ∇ · q + ρr, dt 2
(3.37)
CE
where e is the internal energy, A1 is the first Rivlin-Ericksen tensor [30], q is the heat efflux vector and r is the external supply. The derivation of this equation can be found in standard books on continuum mechanics and rheology; for example, see Tanner [31].
AC
In order to derive the above equation from the internal energy distribution function, it is assumed that et is the total energy given by the sum of the internal and kinetic energies, i.e., 1 et = e + |u|2 . 2
Next, the following consistency relations must hold: 13
(3.38)
ACCEPTED MANUSCRIPT 8 X
gαeq = ρet ,
(3.39)
α=0
!
8 X
1 gαeq ξ α = p + ρe + ρ|u|2 u − τ u + q, 2 α=0 8 X
α=0
(3.40)
Gα = ρb · u − ρr.
(3.41)
D0 = ρet ,
D1 = 0,
D2 = 0.
CR IP T
One way of satisfying the above is to assume, as before, that the scalars Dα are such that Dα = D1 , α = 1, 3, 5, 7, and Dα = D2 , α = 2, 4, 6, 8, and set (3.42)
In addition, it is assumed that the vectors Eα are defined through E0 = 0, Eα = E1 , α = 1, 3, 5, 7; Eα = E2 , α = 2, 4, 6, 8, where !
Finally, the scalars Gα must satisfy
1 1 ρb · (ξ α ⊗ ξ α )u − 2 ρr(ξ α ⊗ ξ α ) : 1, 2 2σ 4σ Gα = 0, α = 0, 2, 4, 6, 8.
(3.43)
α = 1, 3, 5, 7, (3.44)
M
Gα =
E2 = 0.
AN US
1 1 E1 = 2 p + ρe + ρ|u|2 u − τ u + q, 2σ 2
(3.45)
ED
Letting b = bx i + by j, it is easy to verify that b · (ξ 1 ⊗ ξ 1 )u = b · (ξ 5 ⊗ ξ 5 )u = σ 2 bx u, b · (ξ 3 ⊗ ξ 3 )u = b · (ξ 7 ⊗ ξ 7 )u = σ 2 by v, (ξ α ⊗ ξ α ) : 1 = σ 2 , α = 1, 3, 5, 7,
PT
(3.46) (3.47) (3.48)
CE
which means that (3.41) is satisfied.
AC
Finally, summing Eq. (3.33) over all α, and using (3.34) - (3.36), one finds that de ∂u ·u+∇· ρ +ρ dt ∂t
8 X
gαeq ξ α
α=0
!
− ρb · u = 0 + O(ε).
(3.49)
Next, !
!
1 ∇ · p + ρe + ρ|u|2 u = ∇p + ρ∇e + ρ(u · ∇)u · u, 2 14
(3.50)
ACCEPTED MANUSCRIPT because ∇ · u = 0. Finally, using indicial notation for convenience and noting the symmetry of the stress tensor, 1 (τij uj ),i = τij,i uj + τij ui,j = τji,j ui + τij (A1 ),ij . 2
(3.51)
Thus, 1 ∇ · (ττ u) = (∇ · τ ) · u + τ : A1 . 2 Combining all of the above, we obtain
CR IP T
!
(3.52)
de 1 ρ + ρa + ∇p − ∇ · τ − ρb · u − τ : A1 + ∇ · q − ρr = 0 + O(ε), (3.53) dt 2
Compressible Fluids
AN US
which reduces to (3.37) when one appeals to the equations of motion (3.29).
The continuity, the equations of motion and the energy equations for compressible fluids can be derived quite easily. First of all, Eqs. (3.4) - (3.8) remain unchanged, while the matrix M in Eq. (3.9) has the form:
M
2 ρu − Txx ρuv − Txy
(3.54)
ED
M= , ρuv − Txy ρv 2 − Tyy
PT
where Txx , Txy , Tyy are the components of the total stress tensor. A glance at Eqs. (3.10) - (3.13) shows that A1 , A2 , B1 , B2 are unchanged, while one needs the following: ρ|u|2 Txx + Tyy + , σ2 σ2
(3.55)
CE
A0 = ρ −
1 1 (ρu2 − Txx ), C22 = 4 (ρv 2 − Tyy ), 4 2σ 2σ 1 C12 = C21 = 4 (ρuv − Txy ). 8σ
AC
C11 =
(3.56) (3.57)
Next, the energy equation (3.37) is now given by ρ
de 1 = T : A1 − ∇ · q + ρr. dt 2
Thus, one requires the following changes be made: 15
(3.58)
ACCEPTED MANUSCRIPT 8 X
!
1 gαeq ξ α = ρe + ρ|u|2 − T u + q, 2 α=0
(3.59)
!
1 1 E1 = 2 ρe + ρ|u|2 − T u + q. 2σ 2
(3.60)
Boundary conditions
Three Dimensional Equations for Continua
PT
4
ED
M
AN US
CR IP T
One of the main advantages of the current approach is that boundary conditions can be incorporated in a manner similar to macroscopic methods, in contrast with other methods utilised for solving LBM equations. The latter employ complicated special relationships for the discrete particle distribution function (fα ) and the internal energy distribution function (gα ) for each kind of boundary conditions and problems [32,11]. For example, methods such as on-grid and mid-grid bounce back are used when the velocity is zero on the boundary; when the boundary is in motion, bounce-back is used along with a set of linear equations to determine the boundary values fα . In the method used here, the boundary conditions of fα and gα can be obtained directly from the macroscopic values on the boundaries due to the relationships of the macroscopic values with fα as seen from (3.3) - (3.7); in the case of gα , see (3.34), (3.38) - (3.41). As a result, in this method, boundary conditions, especially the Dirichlet conditions, can be included in various problems similar to macroscopic methods and no special equations for fα and gα are needed to incorporate the boundary conditions.
AC
CE
In order to derive the equations of continuum mechanics applicable to motions in three dimensions, one can employ a D3Q15 lattice [2]. Of course, any spacefilling, symmetric lattice can be used for the derivation of the equations of continuum mechanics, since the method relies on the expansion of the particle distribution function as a polynomial in the particle velocity; see (4.2) below. In the case of the D3Q15 lattice, there are 15 nodes and one lies at the centre of a cube with a side of 2 units each; six are chosen from the midpoints of the faces of the cube and the remaining eight are the vertices. In the Table, the positions of the nodes and the corresponding microscopic velocities ξ α , divided by σ for convenience, are given. In order to obtain the continuity equation and the equations of motion, we observe that Eqs. (3.2 ) and (3.3) hold and the evolution equation for the 16
ACCEPTED MANUSCRIPT particle distribution function is once again given by (cf. (3.26)):
1 ∂fα + ξ α · ∇ x fα − F α = − (fα − fαeq ), ∂t ετ
(4.1)
CR IP T
with the following quadratic form:
fαeq = Aα + ξ α · Bα + (ξ α ⊗ ξ α ) : Cα ,
Node
Location
ξ α /σ
α0
(0, 0, 0)
0
α1
(1, 0, 0)
i
α2
(0, 1, 0)
j
α3
M
AN US
Table of Nodes and Microscopic Velocities
(-1, 0, 0)
−i
(0, -1, 0)
−j
ED
α4
(0, 0, 1)
k
α6
(0, 0, -1)
−k
α7
(1, 1, 1)
i+j+k
α8
(-1, 1, 1)
−i + j + k
α9
(-1, -1, 1)
−i − j + k
α10
(1, -1, 1)
i−j+k
α11
(1, 1, -1)
i+j−k
α12
(-1 , 1, -1)
−i + j − k
α13
(-1, -1, -1) −i − j − k
α14
(1, -1, -1)
AC
CE
PT
α5
i−j−k
Next, Eqs. (3.3) - (3.9) are modified as follows: 17
(4.2)
ACCEPTED MANUSCRIPT 14 X
fαeq = ρ,
(4.3)
α=0 14 X
fαeq ξ α = ρu,
u = ui + vj + wk,
(4.4)
α=0
α=0
fαeq ξ α ⊗ ξ α = M, 14 X
fα(n) = 0,
n ≥ 1,
fα(n) ξ α = 0,
n ≥ 1.
α=0 14 X
(4.5)
α=0
(4.6)
CR IP T
14 X
(4.7)
In (4.5), M has the matrix form
2 ρu − Txx ρuv − Txy ρuw − Txz
2 . M= ρuv − Txy ρv − Tyy ρvw − Tyz
(4.8)
AN US
ρuw − Txz ρvw − Tyz ρw2 − Tzz
In the above set, ρ is the density, u, v and w are the components of the velocity field u in the x, y and z directions respectively, Txx , Txy , Txz , Tyy , Tyz , Tzz are the stresses which can be defined through any relevant constitutive relation for the continuous medium. The body force distribution functions now satisfy
M
14 X
Fα = 0,
(4.9)
Fα ξ α = ρb.
(4.10)
ED
α=0
and
14 X
α=0
F1 =
CE
F0 = 0,
PT
As before, one choice for the set of Fα is: 1 ρb · ξ 1 = −F3 , 2σ 2
F2 =
(4.12)
AC
1 ρb · ξ 2 = −F4 , 2σ 2 1 F5 = 2 ρb · ξ 5 = −F6 , 2σ Fα = 0, α = 7, · · · , 14.
(4.11)
(4.13) (4.14)
Finally, in (4.2), one has to specify the scalars Aα , the vectors Bα and the matrices Cα . We assume that A1 = Aα , α = 1, · · · 6; A2 = Aα = 7, · · · 14, and define ρ|u|2 Txx + Tyy + Tzz A0 = ρ − 2 + , A1 = A2 = 0. (4.15) σ σ2 18
ACCEPTED MANUSCRIPT Next, B0 = 0, and B1 = Bα , α = 1, · · · 6; B2 = Bα = 7, · · · 14, and B1 =
ρu , 2σ 2
B2 = 0.
(4.16)
Next, once again, C0 = 0, and C1 = Cα , α = 1, · · · 6; C2 = Cα = 7, · · · 14, and C11 C1 = 0
C22 =
0
0
, C22 0 0 0 C33
1 (ρv 2 − Tyy ), 2σ 4
C11 =
1 (ρu2 − Txx ), 4 2σ
C33 =
1 (ρw2 − Tzz ). 2σ 4
0 C12 C13
C2 = , C21 0 C23 C31 C32 0
C23 = C32 =
AN US
Finally,
C12 = C21 =
1 (ρvw − Tyz ), 16σ 4
(4.17)
CR IP T
1 (ρuv − Txy ), 16σ 4
C13 = C31 =
1 (ρxz − Txz ). 16σ 4
(4.18)
(4.19)
(4.20)
ED
M
Repeating the calculations employed in Section 3, one can now derive the continuity equation and the equations of motion for a continuous medium, valid in three dimensions. In order to obtain the energy equation, one begins with
PT
1 ∂gα + ξ α · ∇x gα − Gα = − (gα − gαeq ), ∂t ετ
(4.21)
CE
where gαeq has a polynomial expansion: gαeq = Dα + ξ α · Eα .
(4.22)
And,
AC
gα = gαeq + εgα(1) + ε2 gα(2) + O(ε3 ), with the requirement that 14 X
gα(n) = 0,
α=0
n ≥ 1.
(4.23)
(4.24)
The energy equation applicable to a continuous media is given by ρ
de 1 = T : A1 − ∇ · q + ρr. dt 2 19
(4.25)
ACCEPTED MANUSCRIPT In order to derive the above equation from the internal energy distribution function, it is assumed that et is the total energy given by the sum of the internal and kinetic energies, i.e., 1 et = e + |u|2 . 2
(4.26)
Next, the following consistency relations must hold:
gαeq = ρet ,
(4.27)
α=0 14 X
gαeq ξ α
α=0 14 X
α=0
CR IP T
14 X
!
1 = ρe + ρ|u|2 u − Tu + q, 2
(4.28)
Gα = ρb · u − ρr.
(4.29)
D0 = ρet ,
AN US
One way of satisfying the above is to assume, as before, that Dα = D1 , α = 1, · · · , 6, and Dα = D2 , α = 7, . . . 14, and set D1 = 0,
D2 = 0.
(4.30)
In addition, it is assumed that E0 = 0, Eα = E1 , α = 1, · · · , 6; Eα = E2 , α = 7, · · · , 14, where
M
!
ED
1 1 E1 = 2 ρe + ρ|u|2 u − Tu + q, 2σ 2
E2 = 0.
(4.31)
Finally, G0 = 0, and
α = 1, · · · , 6, (4.32) (4.33)
CE
PT
1 1 ρb · (ξ α ⊗ ξ α )u − 2 ρr(ξ α ⊗ ξ α ) : 1, 2 2σ 4σ Gα = 0, α = 7, · · · , 14.
Gα =
AC
Once again, repeating the calculations in Section 3, the energy equation (4.25) can be derived.
5
Concluding Remarks
Using an evolution equation for the modified particle distribution functions, equations of continuum mechanics in two and three dimensions have been derived. It is obvious that the choice of the constitutive equation for the incompressible or compressible fluid is unrestricted. Focussing attention on 20
ACCEPTED MANUSCRIPT
CR IP T
two-dimensional flows for convenience, it can be seen from the definition of the matrix M in (3.9) that the extra stresses τxx , τxy , τyy , in an incompressible fluid, can be chosen at will. That is, the fluid model can be Newtonian or non-Newtonian. Similar observations apply to the matrix M in (4.8). The derivations given here, based on that in [1,2], are applicable to all fluids, compressible or incompressible. As far as incompressible fluids are concerned, it has to be emphasised again that in deriving the continuity equation (3.22) and the momentum equations (3.25) or (3.29), the density ρ and the pressure p are independent of each other. Secondly, the equations of continuum mechanics appear at the expansion of the particle distribution functions at the first order, as described in full in (3.3) - (3.53) and (4.1) - (4.33).
AN US
Next, some remarks regarding the numerical scheme are in order. As is well known, there are two main categories of flows in fluid mechanics. One is where the pressure gradient is specified, which occurs in the case of the flow in a pipe. The second is where the pressure field has to be determined as a part of the solution; consider the flow in a lid-driven cavity as an example. Both types of problems have to be solvable with the chosen numerical scheme for the determination of the distribution functions fα and gα .
M
For instance, consider the evolution equation for fα in (3.26). This is solved by using the splitting method of Toro [33] which is explained in detail in [1,2]. Briefly, the evolution equation is split into two parts. The first one is called the streaming section, equivalent to solving the homogeneous equation:
ED
∂fα + ξ α · ∇x fα − Fα = 0. ∂t
(5.1)
AC
CE
PT
In order to solve this equation, an initial value for fαeq is required. To determine this, one begins with an assumed velocity distribution in the domain, say u0 , from which one determines the initial stress field through the designated constitutive equation. If the pressure field is prescribed, this can be used to calculate the matrix M from (3.9) and the initial value of fαeq from (3.3). If the pressure has to be determined as a part of the solution, one assumes that, in addition to u0 , a value of fαeq is given. From (3.9), the pressure p can be determined from the trace of the matrix M in (3.9). Next, using the initial value of fαeq , determined from the given initial macroscopic variables as in (3.3) - (3.9), the intermediate value of fαI is calculated by solving the homogeneous equation; for instance, the Lax-Wendroff scheme can be employed. Using this, the intermediate macroscopic value of the velocity field uI is determined from (3.5). The second one is called the collision section, equivalent to solving the time21
ACCEPTED MANUSCRIPT dependent equation: ∂fα 1 = − (fα (x, t) − fαeq (x, t)). ∂t ετ
(5.2)
Using the previously determined fαI as the initial condition, the collision equation is now solved; for instance, the Euler method can be used.
CR IP T
To reiterate, considering two-dimensional flows for convenience, we note that the matrix M in (3.9) can be used in two different ways. If the pressure p is specified and the initial velocity vector u0 is known, one uses the relevant constitutive equation in this matrix to find the initial value of fαeq . On the other hand, if the pressure field has to be determined, the trace of the matrix can be employed; see Eqs. (18a) - (18c) in [1]. Clearly, the constant density assumption is met, for it is inherent in the Eqs. (3.3) - (3.13). Moreover, at each iteration, a new value of the lattice speed σ is chosen employing (5.8) or (5.9) below.
AN US
Next, we shall discuss the stability of the numerical scheme. Multiplying fαeq with |ξξ α |2 /2 and the utilisation of (3.15) - (3.17) leads to 8 X
1 eq τxx + τyy 1 fα |ξξ α |2 = p + ρ |u|2 − . 2 2 α=0 2
(5.3)
M
Next, it is easy to verify from (3.30) - (3.32) that 8 X
ED
α=0
Fα |ξξ α |2 = 0.
(5.4)
(5.5)
PT
Hence, Eq. (5.1) becomes
CE
∂ 1 τxx + τyy σ2 p + ρ|u|2 − + ρ(∇ · u) = O (ε) . ∂t 2 2 2
AC
The Courant-Friedrichs-Lewy (CFL) condition states that ([34] and [35])
K=
u∆t v∆t + ≤ 1. ∆x ∆y
(5.6)
This can be used in (5.5) and we obtain "
#
2p − τxx − τyy |u| + + σ 2 K = O (ε) . ρ 2
22
(5.7)
ACCEPTED MANUSCRIPT Thus, the lattice speed σ must satisfy
σ=
v u u τxx Kc t
+ τyy − 2p − |u|2 , ρ
1 Kc = √ ≥ 1. K
(5.8)
If the pressure p has to be uniquely defined, one requires that τxx + τyy = 0; see (3.14). Thus, in these fluid models, (5.8) reduces to v u u −2p Kc t ρ
− |u|
,
CR IP T
1 Kc = √ ≥ 1 (5.9) K As a result, the value σ is modified and changes in each iteration as defined through (5.9). σ=
2
ED
M
AN US
Finally, in the Finite Difference Lattice Boltzmann Method (FDLBM) adopted here, the iteration and recovery of the pressure field is similar to the SIMPLE method of Patankar and Spalding [36,37]. As is well known, the SIMPLE method is a guess-and-correct procedure for the calculation of the pressure field. In each iteration, the velocity field is obtained from the first guessed pressure field. Next, using the corrected velocity field, it is possible to find the corrected pressure and this process continues till a very small or zero mass residual is obtained, since the zero mass residual demonstrates that the divergence of the velocity vector field is zero. In FDLBM, the criteria, which is the mass residual in the SIMPLE method, is the difference between the sum of the distribution functions and the fixed density; see Eq. (3.4). Thus, the correction of the pressure field and the subsequent correction of the velocity field continues till a small or zero difference exists between this sum and the density.
CE
PT
For an in-depth discussion of these matters, see [1,2]. For an application to the lid-driven, mixed convection flow of a Bingham fluid in a cavity, see [38].
AC
Acknowledgement: The research of G. H. R. Kefayati is supported by an International Postgraduate Research Scholarship from the Flinders University of South Australia. We would also like to thank the three reviewers for their extremely helpful comments which have been useful in revising the paper.
23
ACCEPTED MANUSCRIPT References [1] S. C. Fu, R. M. C. So, Modeled Lattice Boltzmann Equation and the ConstantDensity Assumption. AIAA Journal, 47 (2009), 3038–3042. [2] S. C. Fu, R. M. C. So, R. M. C. Leung, Linearized-Boltzmann-type-equationbased finite difference method for thermal incompressible flow. Computers & Fluids, 6 (2012), 67–80.
CR IP T
[3] Z. Guo, B. Shi, N Wang, Lattice BGK model for incompressible Navier-Stokes equation. J. Comp. Phys., 165 (2000), 288–306.
[4] U. Frisch, B. Hasslacher, Y. Pomeau, Lattice gas automata for the Navier-Stokes equation. Phys. Rev. Lett. 56 (1986), 1505–1508.
[5] U. Frisch, D. d’Hum`eres, B. Hasslacher, P. Lallemand, Y. Pomeau, J-P. Rivet, Lattice gas hydrodynamics in two and three dimensions. Complex Syst. 1(1987), 649–707.
AN US
[6] S. Wolfram, Cellular automaton fluids 1: Basic theory. J. Stat. Phys. 45 (1986), 471–526 .
[7] X. He, L-S. Luo, A priori derivation of the lattice Boltzmann equation. Phys. Rev. E, 55 (1997), R6333–R6336.
M
[8] P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys. Rev. Ser. 2, 94 (1954), 511–525.
ED
[9] P. Welander, On the temperature jump in a rarefied gas. Ark. Fysik, 7 (1954), 507–553. [10] P. J. Davis, P. Rabinowitz, Methods of Numerical Integration. Academic Press, New York, 1984.
PT
[11] Z. Guo, B. Shi, C. Zheng, A coupled lattice BGK model for the Boussinesq equations. Int. J. Numer. Methods Fluids, 39 (2002), 325–342.
CE
[12] M. Junk, W-A Yong, Rigorous Navier-Stokes limit of the lattice Boltzmann equation. Asymptotic Analysis, 35 (2003), 165–186.
AC
[13] M. Junk, A. Klar, L-S. Luo, Asymptotic analysis of the lattice Boltzmann equation. J. Comp. Phys. 210 (2005), 676–704. [14] F. Dubois, Equivalent partial differential equations of a lattice Boltzmann scheme. Comp. Math. Appl. 55 (2008), 1441–1449. [15] F. Dubois, Third order equivalent equation of lattice Boltzmann scheme. Discrete Cont. Dyn. Syst. 23 (2009), 221–248. [16] F. Dubois, P. Lallemand, Towards higher order lattice Boltzmann scheme. J. Stat. Mech. doi:10.1088/1742-5468/2009/06/P06006.
24
ACCEPTED MANUSCRIPT [17] P. Lallemand, L-S. Luo, Theory of the lattice Boltzmann method: Acoustic and thermal properties in two and three dimensions. Phys. Rev E, 68 (2003), 036706. [18] S. Gabbanelli, G. Drazer, J. Koplik, Lattice Boltzmann method for nonNewtonian (powerlaw) fluids, Phys. Rev E 72 (2005) 046312. [19] J. Boyd, J. Buick, S. Green, A second-order accurate lattice Boltzmann nonNewtonian flow model, Journal of Physics A: Mathematical and General, 39 (2006) 14241–14247.
CR IP T
[20] S.P. Sullivan, L.F. Gladden, M.L. Johns, Simulation of powerlaw fluid flow through porous media using lattice Boltzmann techniques, J. Non-Newt. Fluid Mech. 133 (2006) 91–98.
[21] S.P. Sullivan, A.J. Sederman, M.L. Johns, L.F. Gladden, Verification of shear thinning LB simulations in complex geometries, J. Non-Newt. Fluid Mech. 143 (2007) 59–63.
AN US
[22] J. Psihogios, M.E. Kainourgiakis, A.G. Yiotis, A.Th. Papaioannou, A.K. Stubos, A lattice Boltzmann study of non-Newtonian flow in digitally reconstructed porous domains, Trans. Porous Media 70 (2007) 279–292.
[23] M. Yoshino, Y. Hotta, T. Hirozane, M. Endo, A numerical method for incompressible non-Newtonian fluid flows based on the lattice Boltzmann method, J. Non-Newt. Fluid Mech., 147 (2007) 69–78.
M
[24] J.M. Buick, Lattice Boltzmann simulation of power-law fluid flow in the mixing section of a single-screw extruder, Chem. Eng. Sci., 64 (2009) 52–58.
ED
[25] Amir Nejat, Vahid Abdollahi, Koohyar Vahidkhah, Lattice Boltzmann simulation of non Newtonian flows past confined cylinders, J. Non-Newt. Fluid Mech. 166 (2011) 689–697.
PT
[26] C.H. Wang, J.R. Ho, Lattice Boltzmann modeling of Bingham plastics, Physica A 387 (2008) 4740–. [27] K. R. Rajagopal, A. R. Srinivasa, On the nature of constraints for continua undergoing dissipative processes. Proc. R. Soc. A461 (2005), 2785–2795.
CE
[28] R. R. Huilgol, On the definition of pressure in rheology. Rheol. Bull., 78(2) (2009), 12, 14–15, 29.
AC
[29] R. R. Huilgol, Fluid Mechanics of Viscoplasticity. Springer, Berlin Heidelberg, 2015. [30] R. S. Rivlin, J. L. Ericksen, Stress-deformation relations for isotropic materials. J. Rational Mech. Anal. 4 (1955), 323–425. [31] R. I. Tanner, Engineering Rheology, 2nd Ed. Oxford University Press, Oxford, 2000. [32] Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model. Phys. Fluids, 9 (1997), 1591–1598.
25
ACCEPTED MANUSCRIPT [33] E. F. Toro, Riemann solvers and numerical methods for fluid dynamics: a practical introduction, 2nd Ed. Springer-Verlag, New York, 1999. [34] J. Blazek, Computational Fluid Dynamics: Principles and Applications. Elsevier, Amsterdam, 2001. See pp. 347–350. [35] T. Cebeci, J. P. Shao, F. Kafyeke, E. Laurendeau, Computational Fluid Dynamics for Engineers. Springer-Verlag, New York, 2005. See pp. 311–320.
CR IP T
[36] S. V. Patankar, D. B. Spalding, A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int. J. Heat Mass Transfer 15 (1972) 1787–1806. [37] S. V. Patankar, A calculation procedure for two-dimensional elliptic situations. Numerical Heat Transfer 4 (1981) 409–425.
AC
CE
PT
ED
M
AN US
[38] G. H. R. Kefayati, R. R. Huilgol, Lattice Boltzmann Method for the simulation of mixed convection flow of a Bingham fluid in a lid-driven cavity. J. Non-Newt. Fluid Mech. (submitted).
26