Asymptotic expansion for nonlinear LBM

3.1.1 Asymptotic expansion for nonlinear LBM

In order to removethe zeroth order of the residue (3.12), we have to find a functionf⁽⁰⁾ which satisfies

A(f^eq,(0)(f⁽⁰⁾)−f⁽⁰⁾) = 0. (3.20) We can easily verify that equation (3.20) and

f⁽⁰⁾ =f^eq,(0)(f⁽⁰⁾) (3.21)

are equivalent. In fact, any solution of (3.21) satisfies the problem (3.20). On the other hand, multiplying (3.20) withA^† from the left and usingA^†A=P as well as the equality (3.19), we obtain equation (3.21). Notice that

f⁽⁰⁾ =F^L(ρ⁰,u⁰) +F^Q(u⁰,u⁰) (3.22) with anyρ⁰ and u⁰ is a solution of equation (3.21). In view ofρ₀=h1,f⁽⁰⁾i= ρ⁰,u₀ =h1,Vf⁽⁰⁾i=u⁰, the functionsρ⁰,u⁰ are actually the momentsρ₀,u₀ of f⁽⁰⁾. We conclude that defining f⁽⁰⁾ according to (3.22) together with any arbitrarily prescribedρ⁰ =ρ₀ andu⁰=u₀ removes the zeroth order residue in (3.12).

Next, to removethe first order term in residue (3.12), the equation A

f^eq,(1)(f⁽⁰⁾,f⁽¹⁾)−f⁽¹⁾

−ξ⁽¹⁾(f⁽⁰⁾) = 0 (3.23) must be fulfilled by some f⁽¹⁾. We first prove that a necessary and sufficient condition to solve (3.23) is

h1,ξ⁽¹⁾(f⁽⁰⁾)i= 0, h1,Vξ⁽¹⁾(f⁽⁰⁾)i= 0. (3.24) Since A is symmetric and has a kernel with basis {1, v₁, . . . ,v_d}, it follows from (3.23) that

h1,ξ⁽¹⁾(f⁽⁰⁾)i= 0, hvα,ξ⁽¹⁾(f⁽⁰⁾)i= 0, (3.25) i.e.,ξ⁽¹⁾(f⁽⁰⁾) must be orthogonal to the kernel ofA. Moreover, (3.25) is equiv-alent to (3.24) by applying (2.15), hence (3.24) is a necessary condition to solve (3.23).

On the other hand, if (3.24) holds, it follows that

Pξ⁽¹⁾(f⁽⁰⁾) =ξ⁽¹⁾(f⁽⁰⁾), Qξ⁽¹⁾(f⁽⁰⁾) = 0, (3.26) so that any solution of

f⁽¹⁾ =f^eq,(1)(f⁽⁰⁾,f⁽¹⁾)−A^†ξ⁽¹⁾(f⁽⁰⁾) (3.27) satisfies equations (3.23). This can be easily verified by inserting (3.27) into (3.23) and applying (3.26) as well asAA^† =P. Conversely, solutions of (3.23)

satisfy (3.27) too, which is a consequence of (3.19) as well as A^†A = P and (3.26).

Based on (3.27), we can define f⁽⁰⁾ and f⁽¹⁾. To see this, we first notice that (3.24) is equivalent to partial differential equations for u₀ and ρ0 by taking (3.21) andξ⁽¹⁾(f⁽⁰⁾,f⁽¹⁾) = (V· ∇)f⁽⁰⁾ into account. They are

h1,ξ⁽¹⁾(f⁽⁰⁾)i=∇ ·u₀= 0,

h1,Vξ⁽¹⁾(f⁽⁰⁾)i=c²_s∇ρ₀+ (u₀· ∇)u₀ = 0. (3.28) Assuming we have momentsu₀andρ₀satisfying (3.28), and arbitrary functions u₁ and ρ₁, thenf⁽¹⁾ defined by

f⁽¹⁾ =F^L(ρ₁,u₁) +F^Q(u₀,u₁)−A^†ξ⁽¹⁾(f⁽⁰⁾) (3.29) satisfies (3.23). This is checked by direct computation using (3.7) and (3.8).

Hence (3.24) is also a sufficient condition to findf⁽¹⁾.

We stress that the necessary and sufficient condition to solve (3.23) for f⁽¹⁾ gives a constraint on the lower order moments. More preciselyf⁽⁰⁾ is restricted by equation (3.22) with the moments u₀ and ρ₀ governed by (3.28). In the following, we choose a particular solution of (3.28),

u₀ = 0, ρ₀ ≡ρ¯₀, ρ¯₀ is a constant, (3.30) which turns out to be compatible with the regime in which lattice Boltzmann algorithms are usually applied. In fact, a relation to the incompressible Navier-Stokes equation is only obtained in this regime. Hence

f⁽⁰⁾ = ¯ρ₀f^∗ (3.31)

is constant intandx, so that the derivative terms ofu₀,ρ₀ andf⁽⁰⁾ disappear in the residue (3.12) and the functions ξ^(k) are actually independent of f⁽⁰⁾. For f⁽¹⁾ defined in (3.29) which removes the 1st order residue in (3.12), it has the consequence

f⁽¹⁾=f^eq,(1) =f^L(1). (3.32)

Next, we remove the second order term in residue (3.12). f⁽²⁾ is required to satisfy

A

f^eq,(2)(f⁽⁰⁾,f⁽¹⁾,f⁽²⁾)−f⁽²⁾

−ξ⁽²⁾(f⁽⁰⁾,f⁽¹⁾) = 0, (3.33) ξ⁽²⁾(f⁽⁰⁾,f⁽¹⁾) = (V· ∇)f⁽¹⁾. (3.34) Proceeding as above, we find that a sufficient and necessary condition to solve (3.33) for f⁽²⁾ is

h1,ξ⁽²⁾(f⁽¹⁾)i= 0, h1,Vξ⁽²⁾(f⁽¹⁾)i= 0, (3.35) which are actually equations for the moments ρ₁ and u₁,

∇·u₁= 0, ∇ρ₁= 0. (3.36)

3.1.1 Asymptotic expansion for nonlinear LBM 27 At this point, we restrict our considerations to expansions with a constant ρ₁ = ¯ρ₁. It turns out, that these expansions are still flexible enough to describe the behavior of the lattice Boltzmann solutions to high orders. (In contrast, settingu₁ to be zero would generally stall the expansion at order 3).

In addition, under the condition (3.35), equation (3.33) is equivalent to f⁽²⁾ =f^eq,(2)(f⁽⁰⁾,f⁽¹⁾,f⁽²⁾)−A^†ξ⁽²⁾(f⁽⁰⁾,f⁽¹⁾). (3.37) Prescribing arbitrary values of ρ₂ and u₂, and restricting ρ₁ = ¯ρ₁ and u₁ following the equations (3.36), we can define

f⁽²⁾=F^L(ρ₂,u₂) +F^Q(u₁,u₁)−A^†ξ⁽²⁾(f⁽⁰⁾,f⁽¹⁾). (3.38) which satisfies (3.37). As a result, the second order term in the residue (3.12) is removed.

Continuing to removethe third order term in the residue (3.12), we seek f⁽³⁾ to satisfy

A(f^eq,(3)(f⁽⁰⁾, . . . ,f⁽³⁾)−f⁽³⁾)−ξ⁽³⁾(f⁽⁰⁾,f⁽¹⁾,f⁽²⁾) = 0, ξ⁽³⁾(f⁽⁰⁾,f⁽¹⁾,f⁽²⁾) = (V· ∇)f⁽²⁾+1

2(V· ∇)²f⁽¹⁾+∂tf⁽¹⁾−g. (3.39) Again we find a sufficient and necessary condition to solve (3.39)

h1,ξ⁽³⁾(f⁽¹⁾,f⁽²⁾)i= 0, h1,Vξ⁽³⁾(f⁽¹⁾,f⁽²⁾)i= 0. (3.40) We first note that

h1,(V· ∇)²f⁽¹⁾i=∇ ⊗ ∇:p⁽¹⁾ =c²_s∇²ρ₁ = 0,

wherep⁽¹⁾=h1,V⊗Vf⁽¹⁾i, and the equation (3.36) forρ₁ is applied. Due to the symmetry ofVand f^∗

h1,Vf^∗i=0, h1,gi=c⁻²_s G· h1,Vf^∗i= 0, so that the actual form of the first equation in (3.40) is

∇·u₂= 0. (3.41)

Let us see the second equation in (3.40) which gives

∂tu₁+∇·p⁽²⁾+1

2h1,V(V· ∇)²f⁽¹⁾i=G. (3.42) with p⁽²⁾ = h1,V⊗Vf⁽²⁾i. Applying (2.28) to the second order derivative term, we obtain

h1,V(V· ∇)²f⁽¹⁾i=c²_s∇²u₁. (3.43) Proceeding further we need to know p⁽²⁾ which is based on the second order coefficient

f⁽²⁾ =F^L(ρ₂,u₂) +F^Q(u₁,u₁)−A^†(V· ∇)f^L(1). (3.44)

Taking∇ρ₁ = 0 into account and introducing the symmetric part of the velocity derivative

S[u] = (∇u+∇u^T)/2, (3.45)

we find

A^†(V· ∇)f⁽¹⁾ =c⁻²_s S[u₁] :A^†(V⊗Vf^∗). (3.46) Since the trace of S[u1] equals 2∇·u₁ which vanishes due to the incom-pressibility condition in (3.36), V⊗V can be replaced by its traceless part Λ =V⊗V− |V|²/dI. Applying the property (iv) ofA, we have

A^†(V· ∇)f⁽¹⁾=µc⁻⁴_s S[u1] : Λf^∗, (3.47) which gives the last term on the right hand side of (3.44) a more explicit form.

In view of (2.28), (2.32) and (2.35), we can now computep⁽²⁾ from (3.44) which gives rise to

p⁽²⁾ =c²_sρ₂I+u₁⊗u₁−2µS[u₁]. (3.48) With substitution of (3.43) and (3.48) into (3.42), the second equation in (3.40) is also reformulated to an explicit form.

In summary, the sufficient and necessary condition (3.40) is equivalent to the equations

∇·u₂= 0,

∂tu₁+u₁· ∇u₁+c²_s∇ρ₂ =ν∆u₁+G, (3.49) where we have used the definition

µ=ν+1 2c²_s. The function f⁽³⁾ has to be a solution of

f⁽³⁾ =f^L(3)+f^Q(3)−A^†ξ⁽³⁾(f⁽⁰⁾,f⁽¹⁾,f⁽²⁾). (3.50) To obtain a specific form of the possible solutions to (3.50), we assume thatρ₁=

¯

ρ₁is constant, and thatu₁andc²_sρ₂solve our target problem, the incompressible Navier-Stokes equation (2.1). Further, we choose a divergence free field u₂ in the definition of f⁽²⁾. Thenξ⁽³⁾ and f^Q(3) can be computed and we set

f⁽³⁾=F^L(ρ₃,u₃) +F^Q(u₁,u₂)−A^†ξ⁽³⁾(f⁽¹⁾,f⁽²⁾) (3.51) with arbitrary values ofρ₃as well asu₃, which is a solution of (3.50) and hence satisfies (3.39). In this way, also the third order term is removed.

As forthe fourth order residue,f⁽⁴⁾ is searched so that

A(f^eq,(4)(f⁽⁰⁾, . . . ,f⁽⁴⁾)−f⁽⁴⁾)−ξ⁽⁴⁾(f⁽⁰⁾, . . . ,f⁽³⁾) = 0. (3.52) The sufficient and necessary condition to solve (3.52) is consequently

h1,ξ⁽⁴⁾i= 0, h1,Vξ⁽⁴⁾i= 0, (3.53)

3.1.1 Asymptotic expansion for nonlinear LBM 29

which is equivalent to the equation governing the moments ρ₃ andu₃,

∇·u₃=−c²_s∂_tρ₂−1

2∇·G, (3.54)

∂tu₂+ (u₁· ∇)u₂+u₂· ∇u₁+c²_s∇ρ₃ =ν△u₂. (3.55) With arbitrary values of ρ4 and u₄, and ρ3, u₃ as well as u₂ determined by (3.54),(3.55), (3.41) and the other already known lower order moments, the functionf⁽⁴⁾ is defined by

f⁽⁴⁾) =F^L(ρ₄,u₄) +F^Q(u₁,u₃) +F^Q(u₂,u₂)−A^†ξ⁽⁴⁾(f⁽⁰⁾, . . . ,f⁽³⁾). (3.56) It solves the equation

f⁽⁴⁾ =f^eq,(4)(f⁽⁰⁾, . . . ,f⁽⁴⁾)−A^†ξ⁽⁴⁾(f⁽⁰⁾,f⁽²⁾,f⁽³⁾) and removes the 4th order term in (3.12).

Going tothe 5th order termin the residue (3.12), we definef⁽⁵⁾ by

f⁽⁵⁾=F^L(ρ₅,u₅) +F^Q(u₁,u₄) +F^Q(u₂,u₃)−A^†ξ⁽⁵⁾(f⁽⁰⁾, . . . ,f⁽⁴⁾), (3.57) in whichρ₅ andu₅ are taken arbitrarily andρ₄,u₃ as well asu₄ are controlled by the following equations

∂tu₃+u₁· ∇u₃+u₃· ∇u₁+c²_s∇ρ₄ =ν△u₃+R (3.58)

∇·u₄= 0 (3.59)

with

R= (∂tρ₂+ 1

2∇·G)u₁+B₁G+B₂u₁+B₃ρ₂+B₄u₁⊗u₁

+B5ρ3+B6u₂+B7u₂⊗u₂. Here B₁, B₂, B₃, B₄, B₅, B₆ are higher order differential operators: B₁ and B₅ are linear combinations of∂tand∇²,B2 is a linear combination of∂_t^a(∇)^b with 2a+b= 4, andB₃,B₄ andB₆ are linear combinations of∂_t^a(∇)^b with 2a+b= 3; B₇ is a first order spatial derivative operator. The related combination parameters depend on the constant equilibrium f^∗ and the linear mappingA.

Because they act on the solution (u, p) = (u₁, c²_sρ₂) and the data G of the target Navier-Stokes problem, R plays the role of an explicitly given source term in the equation foru₃.

Successively we could define higher order coefficientsf^(m) to remove the higher order residues. As a result, similar equations for the higher order momentsρm

and u_m of f^(m) appear. The emerging source terms in their equations depend on the temporal and spatial derivatives of the lower order moments and the force like (3.54), (3.58) foru₃.

It is remarked that the missing information on initial and boundary conditions for u₁, ρ₂, u₂, ρ₃, u₃ and ρ₄ as well as the other moments will follow from corresponding parts of the lattice Boltzmann algorithm. The smoothness of

these moments is determined by the regularity of the initial and boundary values and source terms in the corresponding differential equations. In particular,f^(m) are smooth enough if these data are sufficiently smooth.

It is noted that in view of u₀= 0, ρ₀ = ¯ρ₀,ρ₁= ¯ρ₁, the moments of the lattice Boltzmann solution have the expansions

ˆ

u=hu₁+h²u₂+h³u₃+. . . , (3.60) ˆ

ρ= ¯ρ₀+h¯ρ₁+h²ρ₂+h³ρ₃+. . . . (3.61) In case u₂ 6= 0, only a first order accurate Navier-Stokes velocity can be ex-tracted from the numerical solution

u= 1

huˆ−hu2−h²u₃+. . . . (3.62) In case thatu₂ has irregular initial or boundary values,f⁽²⁾ cannot be regular, so that ρ₂ is expected to be not smooth either. Consequently, we expect the numerical result to be inconsistent to the Navier-Stokes pressure.

Further, we note that the first order error term u₂ in velocity is controlled by a homogeneous generalized Oseen problem ((3.41), (3.55)). Therefore u₂ = 0 is the solution if the initial and boundary values of u₂ are zero, so that the velocity field of the Navier-Stokes equation is approximated with a second order accuracy.

As for the pressure, the value of ρ₃ behaves as the first order error term. In case that irregular initial or boundary values of ρ₃ or u₃ can not be avoided, f⁽³⁾ has no regularity in return. Therefore we expect that the pressure can be obtained from the numerical solution by

p= c⁻²_s

h² (ˆρ−ρ¯₀−h¯ρ₁)−hc⁻²_s ρ₃+. . . , (3.63) with at most first order accuracy. On the other hand, when a smooth f⁽³⁾ is obtained and if ρ₃ is also accompanied with zero initial values, we find ρ₃ = 0 so that a second order accurate pressure field can be extracted from the lattice Boltzmann solution too.

The next higher order errors u₃ and ρ4 are determined again by equations of Oseen type ((3.54) and (3.58)). Introducing a field ¯u₃ = ∇ϕ−G/2 where ϕ satisfies the equation ∆ϕ=−∂tρ₂, we can rewrite theu₃-equation as equation for the incompressible field w = u₃−u¯₃ (for details see [46]). This equation can then be solved with zero initial value for w and eventually yields u₃ and ρ₄. But sinceR in generally not zero, the solutions u₃ and ρ₄ can not be zero and thus lead to non-trivial error terms. Moreover, the velocity erroru₃ is not incompressible (see equation (3.54)) which actually reflects the compressibility error of the lattice Boltzmann scheme. As a result, the accuracy of the Navier-Stokes solutions extracted from the moments of the lattice Boltzmann solution is restricted by 2nd order for both velocity and pressure.

3.1.1 Asymptotic expansion for nonlinear LBM 31

In summary, we have the following list of possible accuracy constellations Table 3.1: Possible accuracy of Navier-Stokes solutions extracted from the mo-ments of lattice Boltzmann solution.

expected accuracy of (u,p) cases

1, 0 u₂ 6= 0 andρ₂ irregular 2, 1 u₂ = 0 andρ₃ irregular 2, 2 u₂ = 0 andρ₃= 0

For a more careful discussion of possible accuracy orders we refer to section 3.3 and section 3.5.

From the above information, the first several f^(k) (k ≤ 5) in (3.1) can be defined under the assumption that the involved partial differential equations possess sufficiently smooth solutions.

In particular, the coefficients with trivial lower order momentsρ₀= ¯ρ₀,u₀ = 0 and ρ₁ = ¯ρ₁ are of our special interest, which render (u₁ =u,c²_sρ₂ =p) to be the solution of the impressible Navier-Stokes problem (2.1), (u₂,ρ₃) to be the solution of equation (3.41), (3.55), and (u₃,c²_sρ₄) to satisfy the equation (3.54), (3.58). We set u₅ = 0, ρ₅ = 0 and u₄ = 0 as simple, specific choice. Using these moments in the formula of

f^(m) =F^L(ρm,um) + X

k+l=m

F^Q(uk,ul)−A^†ξ^(m), m= 1, . . . ,5, we arrive at the definition of the first five coefficients in the expansion (3.1) for the nonlinear lattice Boltzmann solutions:

f⁽⁰⁾= ¯ρ0f^∗; Combining them, we thus obtain a truncated asymptotic expansion of the lattice Boltzmann solutions

fˆ(n,j) =f⁽⁰⁾(tn,xj) +hf⁽¹⁾(tn,xj) +· · ·+h⁵f⁽⁵⁾(tn,xj). (3.65)

Due to the construction, the truncated expansion fˆremoves the residue terms in (3.10) up to orderh⁵ and thus satisfies the lattice Boltzmann update rule up to terms of orderh⁶. We summarize this result in the following theorem.

Theorem 1 Assume that solutions u, p of the Navier-Stokes equation (2.1), solutionsu₂, ρ₃ of (3.41),(3.55) and solutionsu₃, ρ₄ of (3.54),(3.58) have the following regularity

u∈C⁵([0, T],Ω),¯ p∈C⁴([0, T],Ω),¯

u₂ ∈C⁴([0, T],Ω),¯ ρ₃ ∈C³([0, T],Ω),¯ u₃ ∈C³([0, T],Ω),¯ ρ₄ ∈C²([0, T],Ω),¯ then the prediction function fˆdefined by (3.65) and (3.64) together with arbi-trary constants ρ¯₀ and ρ¯₁ satisfies

fˆi(n+ 1,j+ci) = ˆfi(n,j) +Ji( ˆf)(n,j) +gi(n,j) + ˆri(n,j) (3.66) at any node x_j ∈Ω withx_j+hci ∈Ω and there is a h−independent constant κ≥6 such that

|rˆi(n,j)|=O(h^κ), nh²≤T. (3.67) The proof is very straightforward. Inserting the prediction ˆf into the scheme (2.20) and using Taylor’s theorem at (tn,x_j) for each f_i^(k), k= 1, . . . ,5, with the Lagrange form of the truncation error at order 6−k, we find the uniformly bounded expression

ˆ

ri(n,j) =h⁶

5

X

k=1

D_6−kf_i^(k)(tξ_k,i(n),xη_k,i(j)), (3.68) witht_ξk,i(n) =tn+ξ_k,ih², ξ_k,i∈[0,1] andxη_k,i(j) =x_j+η_k,ihci, η_k,i ∈[0,1].

Im Dokument Analysis of Lattice Boltzmann Boundary Conditions (Seite 39-46)