Asymptotic Analysis of the Lattice Boltzmann Method for Generalized - - PowerPoint PPT Presentation

✬ ✫ ✩ ✪

Asymptotic Analysis of the Lattice Boltzmann Method for Generalized Newtonian Fluid Flows

Wen-An Yong

Tsinghua University

(joint with) Zaibao Yang Sino-German Symposium on Advanced Numerical Methods for Compressible Fluid Mechanics and Related Problems

Beijing, May 2014

✬ ✫ ✩ ✪ Outline (7 parts)

• Generalized Newtonian Fluids
• Lattice Boltzmann Method
• Asymptotic Expansion
• Theorem
• Corollaries
• Observation
• Conclusions

✬ ✫ ✩ ✪

• 1. Generalized Newtonian Fluids

Newtonian fluids: air, water (in ordinary conditions) non-Newtonian fluids: shampoo, toothpaste, silly putty, blood, polymers, foams, ceramics, cement slurry, paint, food like oil, vinegar, ketchup, mayonnaise, yogurt, etc. The mechanical behaviors of non-Newtonian fluids are quite different from those of the Newtonian ones. For example, the stress tensor of a non-Newtonian fluid depends on the rate-of-deformation tensor NOT in a simple linear fashion.

✬ ✫ ✩ ✪ Generalized Newtonian fluids are among the most common non-Newtonian fluids. The macroscopic motion of a generalized Newtonian fluid is governed by the Navier-Stokes equations ∇ · v = 0, (mass) vt + v · ∇v + ∇p = ∇ · T + external force, (momentum). Here v = v(x, t): velocity, p = p(x, t): pressure, T = µ(˙ γ)S: stress tensor,

✬ ✫ ✩ ✪ where the viscosity µ = µ(˙ γ) is not a constant (Newtonian), but a non-negative function of the shear rate ˙ γ = √ tr(S2)/2 with S the rate-of-strain tensor or rate-of-deformation tensor S = ∇v + (∇v)t, where the superscript t indicates the transpose.

Ref.: R. B. Bird & R. C. Amstrong & O. Hassager, Dynamics of polymeric liquids, Vol. 1: fluid mechanics (2nd Ed.), John Wiley & Sons, Inc., York, 1987 (Chapter 4).

✬ ✫ ✩ ✪ Examples: The power-law model: µ(˙ γ) = µp ˙ γn−1. (n < 1: shear-thinning or pseudo-plastic fluid, n > 1: shear-thickening

• r dilatant fluid, n = 1: classical Newtonian fluid.)

The Bingham plastic fluid: µ(˙ γ) =    µp + τ0

˙ γ

if |T| ≥ τ0 ∞ if |T| < τ0. Here µp is the plastic viscosity, τ0 is the yield stress, and |T| is the magnitude of the stress tensor T.

✬ ✫ ✩ ✪

• 2. Lattice Boltzmann Method (LBM)

The general form is fi(x + cih, t + δt) − fi(x, t) = Ωi(f1, f2, · · · , fN) for i = 0, 1, 2, · · · , N. Here N is a given integer, fi = fi(x, t) is the i-th density distribution function of particles at the space-time point (x, t), ci is the i-th given velocity, h is the lattice spacing, δt is the time step, and Ωi = Ωi(f1, f2, · · · , fN) is the i-th given collision term.

✬ ✫ ✩ ✪ Motivated by the diffusive-scaling analysis (De Masi & Esposito & Lebowitz, 1989; Sone, 2002; Inamuro & Yoshino & Ogino, 1997; Junk & Y., 2003; Junk & Klar & Luo, 2005), we take δt = h2 in what follows, implying that the non-dimensional speed h/δt = 1/h → ∞ as h → 0. This indicates that we are using transport (advection) processes (LBM) to approximate the diffusion ones (NS equations).

✬ ✫ ✩ ✪ For simplicity, we only consider the D2Q9 lattice: 2D problem, N = 8 and ci =        (0, 0) for i = 0 (1, 0), (0, 1), (−1, 0), (0, −1) for i = 1, 2, 3, 4 (1, 1), (−1, 1), (−1, −1), (1, −1) for i = 5, 6, 7, 8. For i ∈ {0, 1, 2, · · · , 8}, we define ¯ i =        for i = 0 3, 4, 1, 2 for i = 1, 2, 3, 4 7, 8, 5, 6 for i = 5, 6, 7, 8. It is clear that ci = −c¯

• i. In this sense, ci is said to be odd.

✬ ✫ ✩ ✪ The first collision term is (Aharonov & Rothman, 1993; Gabbanelli et

• al. 2005; Psihogios et al. 2007; · · ·)

Ωi = 1 τ(S/h) ( f eq

i (ρ, v) − fi

) . Here the relaxation time is taken as τ(S) = 3µ(˙ γ) + 1 2, while the equilibrium distribution is quite standard: f eq

i

= f eq

i (ρ, v) = wi[ρ + 3ci · v + 9

2(ci · v)2 − 3 2v · v]

✬ ✫ ✩ ✪ (· indicates the inner product) with ρ =

8

∑

i=0

fi, v =

8

∑

i=0

cifi and weighted coefficients wi = 1 36        16 for i = 0 4 for i = 1, 2, 3, 4 1 for i = 5, 6, 7, 8. Note that wi = w¯

i, that is, wi is even.

✬ ✫ ✩ ✪ For convenience, we decompose the equilibrium distribution above as f eq

i

= fiL(ρ, v) + fiQ(v, v) with fiL(ρ, v) := wi(ρ + 3ci · v), fiQ(u, v) := wi[ 9

2(ci · u)(ci · v) − 3 2u · v].

Since ci is odd and wi is even, it is easy to see that ∑

i fiL(ρ, v) ≡ ρ,

∑

i fiQ(u, v) ≡ 0,

∑

i cifiL(ρ, v) ≡ v,

∑

i cifiQ(u, v) ≡ 0.

✬ ✫ ✩ ✪ The second collision term (M. Yoshino & Y. Hotta & T. Hirozane &

• M. Endo, 2007; C.H. Wang & J.R. Ho, 2008)

Ωi = f eq

i (ρ, v; S, h) − fi

τ with a constant relaxation time τ. The equilibrium distribution is f eq

i (ρ, v; S, h) = fiL(ρ, v) + fiQ(v, v) + wihA(S

h )S : ct

ici

with scalar function A(S) = 3 2(τ − 1 2) − 9 2µ(˙ γ). Here S : ct

ici is the standard contraction of two symmetric tensors S

and ct

ici.

✬ ✫ ✩ ✪

• 3. Asymptotic Expansion

Motivated by G. Strang (1964), we notice that the LB solution fi = fi(x, t; h) depends on the lattice spacing h which is small. Thus, we seek a power-series expansion fi(x, t; h) ∼ ∑

n≥0

hnf (n)

i

(x, t). Different from the Chapman-Enskog expansion!

✬ ✫ ✩ ✪ Referring to the expansion above, we introduce ρ(n) :=

8

∑

i=0

f (n)

i

, v(n) :=

8

∑

i=0

cif (n)

i

, S(n) := (∇v(n))+(∇v(n))T and ρh ∼ ∑

n≥0

hnρ(n), vh ∼ ∑

n≥0

hnv(n), Sh ∼ ∑

n≥0

hnS(n).

✬ ✫ ✩ ✪ Take f (0)

i

= wi. It follows clearly from the even/odd properties of wi and ci that ρ(0) = 1, v(0) = 0, S(0) = 0. By using the Taylor expansion, we can write the the left-hand side of the LBM as a power series fi(x + cih, t + h2) − fi(x, t) = ∑

n≥2

hn

n−1

∑

s=1

Di,sf (n−s)

i

(x, t) with Di,s = ∑

l+2m=s 1 l!m!∂m t (ci · ∇)l a differential operator.

✬ ✫ ✩ ✪ For the first collision, we expand τ(Sh/h) = τ(S(1) + hS(2) + · · ·) ∼ ∑

n≥0 hnF (n).

It is not difficult to see that F (n) is determined with S(l) for l = 1, 2, · · · , n + 1. In particular, F (0) = τ(S(1)). Thus, we may rewrite the LBM as ∑

r≥0 hrF (r) ∑ n≥2 hn ∑n−1 s=1 Di,sf (n−s) i

= − ∑

n≥0 hn[f (n) i

− fiL(ρ(n), v(n))] + ∑

n≥1 hn ∑ p+q=n fiQ(v(p), v(q)).

✬ ✫ ✩ ✪ By equating the coefficient of hk in the two sides of the last equation and using v(0) = 0, we obtain h0 : f (0)

i

= fiL(ρ(0), v(0)) = fiL(1, 0), h1 : f (1)

i

= fiL(ρ(1), v(1)), h2 : τ(S(1))Di,1f (1)

i

+ f (2)

i

= fiL(ρ(2), v(2)) + fiQ(v(1), v(1)), hk : τ(S(1)) ∑k−1

s=1 Di,sf (k−s) i

+ +F (1) ∑k−2

s=1 Di,sf (k−1−s) i

+ · · · + F (k−2)Di,1f (1)

i

+ f (k)

i

= fiL(ρ(k), v(k)) + ∑

p+q=k fiQ(v(p), v(q))

for k ≥ 3.

✬ ✫ ✩ ✪ With this hierarchy of equations, the expansion coefficient f (k)

i

can be uniquely determined in terms of (ρ(l), v(l)) for l = 1, 2, · · · , k. Moreover, (ρ(l), v(l)) can be inductively obtained by solving a hierarchy

• f quasilinear or linear PDEs.

In particular, ρ(1) ≡ 0 and (ρ(2), v(1)) satisfies the Navier-Stokes equations ∇ · v(1) = 0

∂v(1) ∂t

+ ∇ ρ(2)

3

+ v(1) · ∇v(1) = ∇ · [ 1

3

( τ(S(1)) − 1

2

) S(1)]. Namely, v(1) and ρ(2)

3

are the respective velocity and pressure of the generalized Newtonian fluid, for the relaxation time is taken as τ(S) = 3µ(˙ γ) + 1

2.

✬ ✫ ✩ ✪

• 4. Theorem

Assume ρ(2k+1) |t=0= 0 and v(2k) |t=0= 0 for k = 0, 1, 2, · · ·, and the viscosity µ = µ(˙ γ) is smooth and satisfies that ˙ γµ(˙ γ) is increasing with respect to ˙ γ ≥ 0. Then, for periodic boundary-value problems, the expansion coefficients possess the following nice property f (k)

i

= (−1)kf (k)

¯ i

for all i and all k. Remark: Thanks to √ T : T = µ(˙ γ) √ S : S = √ 2˙ γµ(˙ γ), the structural condition is equivalent to that the magnitude of the stress tensor T increases with increasing the shear rate ˙ γ.

✬ ✫ ✩ ✪ Moreover, (ρ(k+2), v(k+1)) with k ≥ 1 can be consecutively obtained by solving the following linear PDEs ∇ · v(k+1) = Gk,

∂v(k+1) ∂t

+ ∇ ρ(k+2)

3

+ v(1) · ∇v(k+1) + v(k+1) · ∇v(1) =

1 3∇ · [(τ(S(1)) − 1 2)S(k+1) + F (k)S(1)] + ˜

Gk, where Gk and ˜ Gk depend only on (ρ(l), v(l)) and their derivatives with l ≤ k. This theorem is valid also for the second collision mechanism.

✬ ✫ ✩ ✪ Besides those in Junk & Yong (2003) and Junk & Klar & Luo (2005), the new crucial points of the proof contain

• Induction on k.
• To show that F (j) = 0 with j odd.
• The structural condition is used to prove a uniqueness conclusion,

that is, v(k) = 0 with k even.

✬ ✫ ✩ ✪

• 5. Corollaries

(1). From the theorem above, we see that the density and velocity moments of the LB solution fi = fi(x, t; h) can be expanded as ρh := ∑

i fi ∼ ∑ n≥0 hnρ(n) = 1 + h2ρ(2) + h4ρ(4) + h6ρ(6) + · · · ,

vh := ∑

i cifi ∼ ∑ n≥0 hnv(n) = hv(1) + h3v(3) + h5v(5) + · · · .

✬ ✫ ✩ ✪ Thus we have

vh h − v(1) = h2v(3) + h4v(5) + · · · , ρh−1 h2

− ρ(2) = h2ρ(4) + h4ρ(6) + · · · . These suggest that the rescaled LB moments vh

h and ρh−1 3h2

should be taken as approximations of the velocity v(1) and pressure ρ(2)/3 of the generalized Newtonian fluids with second-order accuracy in space and first-order accuracy in time. They also provide a basis for using the Richardson extrapolation to improve the accuracy up to higher orders.

✬ ✫ ✩ ✪ (2). Following Yong and Luo (2012), we can deduce from the theorem that τ(S/h)S = 3 ∑

i(f eq i

− fi)ci ⊗ ci h + O(h2). This relation hints an alternative way to implement the LB method with the first collision mechanism, instead of computing the shear-rate tensor S = ∇v + (∇v)T directly from v with finite difference or other

• methods. Indeed, we could use the relation

Σ = 3 ∑8

i=1(f eq i

− fi)ci ⊗ ci h2τ(Σ) to compute Σ := S/h and thereby τ(Σ) by iteration (Boyd et al. 2006, Tang et al. 2009, Chai et al. 2011).

✬ ✫ ✩ ✪ (3). Initial condition fi(0, x) = f eq

i (1, hv0) + 3h2wi{p0 − τ(S0)ci · ∇(v0 · ci)}.

We will have second-order accuracy for the velocity and first-order for the pressure! (4). To have second-order accuracy for the pressure, the initial condition can also be constructed. The key is solving a Poisson equation for u(3) in terms of v0. Details can be found in

Ref.: Juntao Huang & Hao Wu & Y., On initial conditions for the lattice Boltzmann method, preprint (2013).

✬ ✫ ✩ ✪

• 6. Observation

The structural condition that ˙ γµ(˙ γ) is increasing with respect to ˙ γ ≥ 0 can be easily verified for the following 7 models. The power-law model: µ(˙ γ) = µp ˙ γn−1. Here µp is the flow consistency coefficient and n is the power-law index of fluid. The case n < 1 corresponds to a shear-thinning or pseudo-plastic fluid, whereas n > 1 corresponds to a shear-thickening

• r dilatant fluid, and n = 1 reduces to the classical Newtonian fluid.

✬ ✫ ✩ ✪ The Carreau model: µ(˙ γ) = µ∞ + (µ0 − µ∞)[1 + (λ˙ γ)2](n−1)/2 for 0 < n ≤ 1. Here µ0 is the zero-shear-rate viscosity(˙ γ → 0), µ∞ is the infinity-shear-rate viscosity (˙ γ → ∞), and λ is a time constant. Notice that µ0 > µ∞ for shear-shinning fluids. We compute d d˙ γ [˙ γµ(˙ γ)] = µ∞ + (µ0 − µ∞)[1 + λ2 ˙ γ2]

n−3 2 (1 + nλ2 ˙

γ2) > 0 for 0 < n ≤ 1 and µ0 > µ∞.

✬ ✫ ✩ ✪ The Carreau-Yasuda model: µ(˙ γ) = µ∞ + (µ0 − µ∞)[1 + (λ˙ γ)a](n−1)/a, where the parameters have the same meaning as above and the new parameter a is an extra material constant. The Cross model (for shear-thinning): µ(˙ γ) = µ∞ + µ0 − µ∞ 1 + (λ˙ γ)(1−n) . These models can also describe shear-thickening fluids, where n > 1 and the meanings of µ0, µ∞ exchange as follows lim

˙ γ→0 µ(˙

γ) = µ∞, lim

˙ γ→∞ µ(˙

γ) = µ0.

✬ ✫ ✩ ✪ A smoothed version of the Bingham model (T.C. Papanastasiou, 1987): µ(˙ γ) = τ0 ˙ γ (1 − e−m ˙

γ) + ηp.

Here τ0 is the yield stress, m is the stress growth exponent (regularization parameter), and ηp is the plastic viscosity. A smoothed version of the Casson model (P. Neofytou & D. Drikakis, 2003): µ(˙ γ) = [√τ0 ˙ γ (1 − e−√m ˙

γ) + √ηp

]2 . Here the parameters τ0, m and ηp are same as in the last model.

✬ ✫ ✩ ✪ The Powell-Eyring model: µ(˙ γ) = µ∞ + (µ0 − µ∞)sinh−1(λ˙ γ) λ˙ γ , where µ0, µ∞ and λ are material constants. More models?

✬ ✫ ✩ ✪

• 7. Conclusions
• We present a detailed asymptotic analysis of two LB-BGK models

for generalized Newtonian fluids. We do NOT use the Chapman-Enskog expansion, leading to the compressible NS

• equations. Ours gives the incompressible equations directly.
• Our analysis exposes certain important features of the LB

solutions: f (k)

i

= (−1)kf (k)

¯ i

. This has several useful corollaries including a basis for using the Richardson extrapolation to improve the accuracy, an iteration method to compute the the rate-of-deformation tensor, etc.

✬ ✫ ✩ ✪

• We introduce the lattice spacing h into the collision terms.
• We observe that, for many widely used generalized Newtonian

models, the magnitude of the stress tensor increases with increasing the shear rate.

• Our analysis can be directly extended to MRT models, to 3-D

problems, and to rectangular lattices.

✬ ✫ ✩ ✪

THANK You for Your Attention!