A viscosity adaption method for Lattice Boltzmann simulations

Accepted Manuscript A viscosity adaption method for Lattice Boltzmann simulations Daniel Conrad, Andreas Schneider, Martin Böhle

PII: DOI: Reference:

S0021-9991(14)00554-3 10.1016/j.jcp.2014.08.008 YJCPH 5418

To appear in:

Journal of Computational Physics

Received date: 9 January 2014 Revised date: 19 June 2014 Accepted date: 1 August 2014

Please cite this article in press as: D. Conrad et al., A viscosity adaption method for Lattice Boltzmann simulations, J. Comput. Phys. (2014), http://dx.doi.org/10.1016/j.jcp.2014.08.008

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.

A Viscosity Adaption Method for Lattice Boltzmann Simulations Daniel Conrad∗, Andreas Schneider, Martin B¨ohle Institute of Fluid Mechanics and Fluid Machinery, Technische Universit¨ at Kaiserslautern, Gottlieb-Daimler-Str. 67663 Kaiserslautern, Germany.

Abstract In this work, we consider the limited fitness for practical use of the Lattice Boltzmann Method for non-Newtonian fluid flows. Several authors have shown that the LBM is capable of correctly simulating those fluids. However, due to stability reasons the modeled viscosity range has to be truncated. The resulting viscosity boundaries are chosen arbitrarily, because the correct simulation Mach number for the physical problem is unknown a priori. This easily leads to corrupt simulation results. A viscosity adaption method (VAM) is derived which drastically improves the applicability of LBM for non-Newtonian fluid flows by adaption of the modeled viscosity range to the actual physical problem. This is done through tuning of the global Mach number to the solution-dependent shear rate. We demonstrate that the VAM can be used to accelerate LBM simulations and improve their accuracy, for both steady state and transient cases. Keywords: Lattice Boltzmann Method, simulation acceleration, accuracy, adaptive time step, stability, non-Newtonian, pulsatile

1. Introduction By now, Lattice Boltzmann Methods (LBM) are established as an alternative approach to classical computational fluid mechanics methods. The LBM has been shown to be an accurate and efficient tool for the numerical simulation of weakly compressible or incompressible fluids. Fields of application reach from turbulent simulations through thermal problems to acoustic calculations among others. The regular grid based scheme makes the LBM ideally suitable for simulations involving complex solids. Such geometries are common in the food processing industry, where fluids are mixed by static mixers or agitators. Those fluid flows are often laminar and non-Newtonian. Several authors have shown that the LBM is capable of simulating those fluids (e.g.[1]). Special care has to be taken for non-Newtonian Lattice Boltzmann simulations in order to keep them stable. A straight forward way is to truncate the modeled viscosity range by numerical stability criteria [2]. This is an effective approach, but from the physical point of view the viscosity bounds are chosen arbitrarily. ∗

Corresponding author. Email address: [email protected] (Daniel Conrad )

Preprint submitted to Journal of Computational Physics

August 6, 2014

Moreover, these bounds depend on and vary with the grid and time step size and therefore with the simulation Mach number. Thus, the modeled viscosity range may not fit to the range of the physical problem (see Fig.1), because the correct simulation Mach number is unknown a priori. A way around is, to perform precursor simulations on a fixed grid to determine a possible time step size and simulation Mach number, respectively [3]. These precursor simulations can be time consuming especially for complex cases. This makes the LBM unattractive for use in practical simulations of non-Newtonian fluids. To our best knowledge there currently does not exist a consistent technique to adapt the numerical viscosity range to the range of the physical problem. Therefore, in this paper, a viscosity adaption method (VAM) is proposed. Chapter 2 describes the LBM model. In Chapter 3 the adaption method is presented. We verify the proposed viscosity adaptive Lattice Boltzmann Method (VALBM) in Chapter 4 for steady state channel flow and transient pulsatile pipe flow of non-Newtonian fluids. 2. Lattice Boltzmann Method 2.1. Model description The Lattice Boltzmann Method is a special technique to solve the Boltzmann equation numerically in the continuum limit. It is a mesoscopic method due to simplified gas dynamics in which the gas molecules are mathematically captured through a distribution function f = f (x, t). Discretization in space is done with an equidistant lattice and the molecules are restricted to travel to the nearest neighbor lattice nodes, which results in a discretization of the velocity space. In this work we rely on the widely used D3Q19 discretization model [4]. Throughout the paper we use dimensional form of the variables. The distribution function f at a lattice node evolves through advection in direction i of the discretized velocity space and subsequent collision locally on the lattice nodes fi (x + ei Δt, t + Δt) = fi (x, t) + Ci .

(1)

The Boltzmann collision operator can be replaced due to Bhatnagar, Gross, and Kroog (BGK) with a relaxation of f towards a local equilibrium [5]. The BGK approximation of the collision operator reads: Ci = Ω(fieq (x, t) − fi (x, t)) Ω ∈ [0, 2],

(2)

with the dimensionless collision frequency Ω = Δt/τ and the local equilibrium fieq , which is a function of the macroscopic variables density ρ and velocity v. The components e of the molecular velocity vector ei are defined as Δx/Δt. This equilibrium distribution can be viewed as the discretized and truncated Maxwell distribution with the weights wi accounting for the different lattice directions: fieq = wi ρ · (1 +

ei · v (ei · v)2 |v|2 + − 2 ). c2s 2 · c4s cs 2

(3)

The speed of sound cs is deduced from the simulation Mach √ number built with a characteristic flow velocity M a = U/cs at which the relationship 3 · cs = e holds. Through the Chapman-Enskog expansion, it can be shown that the LBM recovers the Navier-Stokes equations for small Mach and Knudsen numbers [4]. An other result of this analysis is the relationship between the relaxation parameter Ω and the kinematic viscosity of the fluid, which eliminates discrete lattice effects to represent the correct transport coefficient in a simulated gas c2s Δt Ω= . (4) ν + 0.5 · c2s Δt The link between the microscopic state and the macroscopic variables is established through the calculation of the moments m of the distribution function. The physically meaningful moments in this model are density and momentum, which are conserved quantities of the system. They are given as the zeroth and first order moments:  ρ= fi (5) i

ρv =



ei f i

(6)

i

The pressure is coupled with the density through the isothermal equation of state p = c2s · ρ. From the non-equilibrium part of the second order moment tensor Παβ the strain rate tensor αβ can be determined αβ = −

 Ω Ω (1) Π ≈ − eiα eiβ (fi − fieq ), Δt · ρ · c2s αβ Δt · ρ · c2s i

(7)

where the distribution function may be split in an equilibrium part and a non-equilibrium part fi = fieq + fineq [6]. In a Newtonian fluid the viscous stress is related to the strain rate through the dynamic viscosity η = ν · ρ, which yields the stress tensor σαβ : σαβ ≈ −

Ων  eiα eiβ (fi − fieq ). 2 Δt · cs i

(8)

Because Ω is the only free parameter in (2), this model is known as the single relaxation time (SRT) model. An improvement in terms of numerical stability is the multiple relaxation time (MRT) model. In this model, the collision operation (2) is reformulated to take place in the moment space [7]: Ci = M −1 S · (M fieq − M fi ). (9) The definition of the Matrix M for the linear transformation is given in Appendix A. In (9) the relaxation parameter is extended to the relaxation Matrix S with the non-zero entries representing independent collision rates sa−e ∈ [0, 2]: S = diag(0, sa , sb , 0, sc , 0, sc , 0, sc , Ω, sd , Ω, sd , Ω, Ω, Ω, se , se , se ) 3

(10)

2.2. Generalized Newtonian model In contrast to Newtonian fluids, where the molecular viscosity is assumed to be constant for isothermal flows, the viscosity of a non-Newtonian fluid varies with the local strain rate. A common approach is the power-law model by Ostwald de Waele [8] τw = k ||n .

(11)

Combined with the linear relationship for wall shear stresses of Newtonian fluids the effective kinematic viscosity of a generalized Newtonian fluid becomes ν=

k n−1 || , ρ

(12)

where k is the flow consistency index [P asn ] and n is the dimensionless flow behavior index. Fluids with n < 1 are called pseudo plastics (shear thinning) and n > 1 are dilatants (shear thickening). Note, that k can be identified as the dynamic viscosity for Newtonian fluids where n = 1. These parameters are fluid dependent and must be determined through curve fitting of experimental data. Incorporation of non-Newtonian behavior can be done by substitution of the kinematic viscosity in (4) with the effective viscosity related to the magnitude of the strain rate tensor according to the constitutive law used. One of the major advantages of the lattice Boltzmann method is the ability to evaluate the stress and strain tensors locally with second order accuracy [9]. The magnitude of the strain rate tensor (7) is calculated through the second invariant:   1 Ω (1) (1) || = · Παβ Παβ ∀α = β. (13) 2 Δt · ρ · cs αβ 2 The kinematic viscosity is then computed from (12). Technically speaking, this expression can not be evaluated directly, because the kinematic viscosity itself appears in the collision frequency (4). For an exact solution the root of this implicit equation has to be computed. However, to keep the algorithm simple the viscosity from the previous time step is used to compute the collision frequency in (13), which has proofed to be a good approximation. It is well known that the LBM suffers from numerical instabilities when low viscosities are present and the relaxation parameter approaches 2. Hence, for the case of shear thinning fluids, the viscosity obtained through (12) must be limited with a lower bound in regions with vanishing magnitude of the strain rate tensor. Similarly, to ensure hydrodynamic behavior, it is necessary to introduce an upper bound on ν [10]. In the case of the power law model, the effective viscosity becomes ⎧ ⎪  > max → Ωmin ⎨νmax n−1 k ν = ρ || (14)  ⎪ ⎩  < min → Ωmax . νmin This effectively results in a lower bound Ωmin and an upper bound Ωmax of the dimensionless collision frequency. In our simulations, we use an upper bound of Ωmax = 1.965. The value 4

for Ωmin is set to 1.0. In fact, the shear rate bounds min and max are computed from Ωmin and Ωmax , respectively. Actual physical boundaries, if known, can easily be implemented in the VALBM by setting the appropriate values in eq.(14). For example, the real upper bound on the viscosity for a shear thinning material corresponds to the horizontal line of the shaded area (physical problem) in Figure 1. In practical simulations, the collision frequency built with the real upper viscosity bound of the material (if known) may lead to values Ωmin  1, which in turn leads to instability of the LBM. Thus, a smaller viscosity bound must be chosen, such that, e.g., Ωmin = 1.0. If the real upper bound results in a collision frequency of Ωmin ≥ 1.0, that value must be chosen for the limits in eq.(14). Hence, the implementation of the VALBM is straight-forward for other non-Newtonian models,e.g., the Carreau-Yasuda or Cross model. 3. Viscosity adaption Method 3.1. Theoretical considerations As already stated above, the correct simulation Mach number is not known beforehand. The modeled viscosity range of a non-Newtonian fluid is restricted through the stability criteria (14). Whereas the level of this range in terms of shear rate is determined by the time step size and therefore by the simulation Mach number. Strictly speaking, the considerations made in this paper are only valid for flows were the Mach number is irrelevant and the simulation Mach number is a free parameter. From the computational view, the Mach number is used to deduce a time step for the computation, which ensures, that the simulation is conducted in the ”quasi compressible” regime. In the LBM, an appropriate time dependent change of the density, and therefore the hydrodynamic pressure, is used to render the continuity equation divergence free. Based on the choice of the simulation Mach number, the modeled range may not lie in the physical range of the problem. This is depicted in Fig. 1. The objective of the VAM is to shift one of the boundaries of Ω to coincide with one of the boundary values of the physical problem. For this purpose, a rescaled time step size Δt∗ is introduced which leads to rescaled lattice values: e∗ =

Δx Δt∗

e∗ c∗s = √ 3

M a∗ =

U c∗s

Ω∗ =

c∗s 2 Δt∗ , ν ∗ + 0.5 · c∗s 2 Δt∗

(15)

where Ω∗ must be rescaled at each lattice site. In order to yield the same macroscopic moments in the new system, the distributions fi must be rescaled. This is done with the linear transformation matrices M in the old system and M ∗ in the new system according to fi∗ = M ∗ −1 M fi .

(16)

One can conclude, that the moments and especially the second order moment tensor Παβ is unchanged by this operation. However, the macroscopic shear rate (7) changes because of the prefactor. Consequently, the kinematic viscosity changes through the constitutive law used in the simulation. It is desired that the rescale operation (16) leaves the shear rate and 5

Kinematic viscosity ν

Shear thinning power-law fluid Physical problem

Ωmin

Ωmax modeled

Shear rate Figure 1: Physical problem and numerical approximation

therefore the viscosity untouched. We achieve this, through rescaling of fi∗ neq in (7) to give the same shear rate, hence ν ∗ = ν, fi∗∗

=

fi∗ eq

+

Δx2 6·Δt∗ Δx2 + 6·Δt

ν+ ν

· (fi∗ − fi∗ eq ),

(17)

where fi∗ eq is computed from (3) in the new system. 3.2. Choice of time step size In order to adapt to the physical problem, the time step size Δt∗ is determined from the minimal or maximal shear rate in the modeled system, which is a solution-dependent quantity. Following the example situation in Fig.1, the numerical approximation of the physical problem is improved, if we shift the maximal collision frequency towards the maximal shear rate of the system. Therefore, the target collision frequency Ωt should take the value of Ωmax built with the kinematic viscosity νm obtained from the maximum shear rate. Thus, we choose the new time step size Δt∗ according to ∗

Δt =

( Ω1t − 12 ) · Δx2 3 · νm

.

(18)

Since the value of νm is solution-dependent, it can be expected to converge to the physical value of the system, if the time step size is recursively adapted. From (15) it follows that an increase in Ω causes a decreased time step size and therefore a decrease in Mach because 6

Δt M a ∝ Δx . Hence, we can expect a prolonged simulation time with an enhanced time resolution. An other possibility is to set the target relaxation parameter by Ωmin . If Δt∗ is determined by (18) with Ωt = Ωmin and νm from the minimal shear rate the time step size as well as the Mach number is increased. Therefore the simulation is accelerated. Note, that in a general flow situation it may then be possible to cut off viscosities from high shear rates due to the Ωmax limit with this approach. The optimal choice for Ωmin will be part of future work. For the simulations conducted in this paper, we use a value of Ωmin = 1.0.

3.3. Mach number limitation As described in 3.2, the viscosity adaptive LBM makes use of a variable time step size on a fixed grid. Given that the Mach number is coupled with the time step size, we have to make sure the lattice values remain valid throughout the simulation. To this end, we limit the possible choice of the time step size through physically meaningful Mach number limits ⎧ max ⎪ M a∗ > M amax ⎨M a ∗∗ M a = M a + λ · (M a∗ − M a) λ ∈ [0, 1] (19) ⎪ ⎩ min ∗ min Ma Ma < Ma . It is common practice to consider flows with Mach numbers M amax below 0.3 to be weakly compressible [11]. Additionally we set M amin = 10−4 to stay in the weakly compressible regime because the hydrodynamic limit does not exist for uncorrected BGK dynamics [12, 13]. In practice, this also prevents the time step size to get uneconomically small. For enhanced stability properties of the viscosity adaptive LBM we compute the simulation Mach number for the next time step M a∗∗ through under-relaxation of the Mach number M a∗ obtained from (15) and (18). The final lattice values, including the time step size, have to be recomputed from M a∗∗ . For stationary flows, the simulation Mach number can be expected to converge as the target shear rate converges to the physical problem. Hence, we are able to introduce a convergence criterion for the simulation Mach number, |M a∗∗ − M a| < δ.

(20)

If the expected Mach number change is small the adaption process can be abandoned. 3.4. Algorithm To precis the procedure, we consider a stationary flow of a shear thinning fluid (Fig.1) with the aim to accelerate the simulation. Starting the algorithm for the viscosity adaption method at an arbitrary time step t0 with an arbitrary simulation Mach number, it can be summarized as follows: - Save the minimal shear rate m present in the system. - Compute νm with m according to the constitutive law e.g (12). - Compute new time step size Δt∗ with the target collision frequency Ωt and νm from (18). 7

- Compute new lattice values from (15). - Relax Ma with under-relaxation factor λ or truncate it (19). - Recompute final adapted lattice values (15) with the time step size obtained from Mach number M a∗∗ and the new viscosity bounds for (14). - Finally rescale distributions through (16) and (17). The algorithm must be repeated in the subsequent time steps, until the Mach number M a∗∗ converges and the desired target collision frequency is reached. This can also be realized in adaption intervals allowing the solution to develop in between VAM steps. 4. Verification 4.1. Steady state case For the verification of the proposed algorithm, we investigate the flow through a 3D duct. We impose periodic boundary conditions in streamwise direction and the flow is driven through a body force [14]. The lattice is constructed as to exactly retain the channel width and height (h), while using the simple bounce back boundary condition. For the sake of simplicity, we model the side walls of the duct to be frictionless. Therefore we can use the plane Hagen-Poiseuille solution for power-law fluids. The velocity distribution for this case is given by [15] n+1 n+1 dp 1 1 n ua (y) = ·( ) n · (h n − y n ), (21) n + 1 dx k where y ∈ [−h, h] is perpendicular to the flow direction. The pressure gradient is substituted dp for the body force dx = ρ · a. To quantify the quality of the simulation, we measure its accuracy with respect to the analytic solution over all lattice sites N,  1  E= · (u(x) − ua (x))2 . (22) N x The simulations in the steady state case are run until the following condition is satisfied, |E(t) − E(t − Δt)| < 10−12 .

(23)

Mach number investigation In a first approach, we want to investigate the convergence properties of the viscosity adaption method and its robustness to the starting Mach number which is unknown user input. A series of simulations with starting Mach numbers from 0.1 to 10−4 built with U=1.0 m/s are conducted for a shear thickening fluid with k = 250, n = 1.25, ρ = 1000 and an acceleration of unity. The target collision frequency for the VAM is set to Ωmin = 1.0. We set the relaxation parameter λ = 0.05 to gain an appropriate resolution of the Mach number 8

10−4

Mach 0.1 Mach 0.01 Mach 0.001 Mach 0.0001

Mach number

10−3

10−2

10−1

100

0

100

200 300 Time steps

400

500

Figure 2: Mach number development over time for different staring Mach numbers.

development in this case. The VAM is started after 50 time steps and repeated every 10 time steps thereafter until the Mach number M ∗∗ converges according to (20) with a residuum of δ = 10−8 . Simulation continues until (23) is satisfied. Figure 2 illustrates the Mach number development for the VALBM. The Mach number is plotted against the number of time steps. It is visible, that the Mach number is adapted to the actual physical shear rate as the solution develops. All simulations yield a converged simulation Mach number M a∗∗ = M aC = 0.0439772 independent of the Mach number at the beginning of the simulation. To verify the theoretical considerations in 3.2, we rerun the same series of simulations without viscosity adaption. Additionally a standard LBM simulation is conducted for the converged Mach number M aC . Figure 3 shows the error in velocity for the different starting Mach numbers. It can be seen, that the simulation Mach number has a major impact on the error of the standard LBM approach. In fact the simulation is performed more or less outside of the discrete physical problem (Fig.1). The quality of the simulation is improved as the simulation Mach number of the standard LBM approaches the regime of the converged Mach number M a∗∗ which is unknown a priori. All viscosity adaptive LBM simulations yield the same error independent of the starting Mach number, since the Mach number is changed to the converged Mach number by VAM (Fig.2). For the converged Mach number, the standard LBM and the VALBM yield the same error. This confirms numerically that the viscosity is not influenced by the rescale operation as proposed in 3.1.

9

101 VALBM Standard LBM

Global absolute error

100 10−1 10−2 10−3 10−4

10−4

10−3 10−2 Mach number

10−1

Figure 3: Error for Mach numbers corresponding to Fig.2 for Standard LBM and the viscosity adaptive LBM.

Accuracy Next, the accuracy of the viscosity adaptive LBM is investigated. We perform simulations for a shear thinning n = 0.75, shear thickening n = 1.25, and a Newtonian fluid n = 1.0 with the boundary conditions and convergence creteria corresponding to the setup in 4.1. The lattice is successively refined to test the spatial convergence. Figure 4 shows a log-log plot of the global absolute error over the grid resolution. Again, the respective result on each lattice is independent of the starting Mach number. Therefore, only already Mach-converged results are shown. It is visible that the viscosity adaptive LBM ensures diffusive scaling of the lattice values [16]. This results in a global second order convergence as indicated by a line corresponing to a −2 slope. 4.2. Transient case For the investigation of the viscosity adaptive LBM in transient cases, we consider the pulsatile flow through a 3D pipe with radius R. Again, the fluid is driven through a spatially constant body force and the simulation domain is periodic in streamwise direction. We impose a second order accurate no-slip boundary condition to model the pipe wall [17]. In this case, the Hagen-Poiseuille solution for the velocity of power-law fluids is given in radial coordinates [15] n+1 n+1 dp 1 1 n ua (r) = ·( ) n · (R n − r n ). (24) n + 1 dx 2k 10

Global absolute error

10

n=0.75 n=1.25 n=1 Slope -2

−2

10−3

10−4

10−5

101

102 Grid resolution N

Figure 4: Accuracy of the viscosity adaptive LBM with Ωt = Ωmin = 1.0 for shear thickening, shear thinning and Newtonian fluids on a successively refined lattice.

The pressure gradient is substituted for the body force is made pulsatile with frequency ω through a(t) = cos(ω · t).

dp dx

= ρ · a at which the acceleration (25)

An important number in this context is the dimensionless dynamic similarity parameter α. It is known as the Womersley number and describes the relation of a characteristic time periodic flow frequency to viscous effects. It can be used to judge transient flows,  ω α=R· . (26) νmin In particular, a pulsatile flow can be considered quasi-steady if α < 1 and the instantaneous velocity distribution is given by (24) if the instantaneous pressure gradient is used [18]. Because the viscosity is not constant in non-Newtonian simulations we compute the maximal Womersley number through the minimal kinematic viscosity in the system to characterize the transient flow. The quality of the transient simulations is measured over one period T of the oscillation and computed in accordance with (22):  1 1  Et = · (u(x, t) − ua (x, t))2 (27) T t N x 11

We are interested in the properties of the VALBM in transient flow situations. For this purpose, we choose a shear thinning material with n = 0.75, k = 250 and ρ = 1000 to simulate the Womersley flow. The pulsation frequency is set to 3 · 10−3 . In contrast to the standard LBM, the simulation Mach number is expected to vary because the VAM adjusts to the varying shear rates of the pulsatile flow. The target collision frequency is set to Ωt = 1.0. The adaption process is repeated every 10 time steps with an under-relaxation factor of λ = 0.8. Like in the steady state case, we conduct standard LBM simulations for constant Mach numbers 0.1, 0.05, 0.01 and 0.005. Figure 5 shows the results for the viscosity adaptive LBM simulation over one pressure oscillation depicted in the upper plot. In the middle plot the development of the Mach and Womersley number is shown. The maximal Womersley number as defined in (26) is α = 0.123625. With this value we can justify the quasi-steady flow approximation. The maximal Mach number yields a value of 0.023859. In order to check if the choice of the adaption interval and under-relaxation factor is reasonable, a viscosity adaptive steady state simulation with the full pressure gradient is attached. It yields a converged Mach number of M aC = 0.0238591. An additional standard LBM simulation for M aC is conducted for the transient case. The plot at the bottom of Fig. 5 shows the development of the global absolute error as defined in (22). The error tends to zero as the pressure gradient vanishes and rises again with increasing values of the pressure gradient. However, the error decreases again when Δp approaches the turning points at 1/2 T and T . This can be attributed to the quasi-steady approximation, because the derivative of the pressure gradient becomes zero. Hence the VAM is able to converge comparable to the steady state case which results in a more accurate result. Simulation quality To judge the performance of the viscosity adaptive LBM, we compute the error Et from (27) and compare the results with standard LBM solutions. Figure 6 shows the time accumulated global absolute error plotted against the number of time steps needed to complete one pressure oscillation. The simulation with the Mach number of 0.1 shows the highest error. This error is composed of a high viscosity error. Also there could be a heightened time discretization error present [13]. The simulations with the Mach numbers 0.01, 0.005 and M aC show comparable results, whereas the latter needs the lowest computational effort. The viscosity adaptive LBM yields the most accurate result with a slightly increased number of time steps. In this case, the additional computational effort in terms of the number of time steps as opposed to the standard LBM simulation with the converged Mach number amounts to 33.8%. Therefore, the accuracy is enhanced by 44.5%. Again, it should be pointed out that the Mach number M aC is unknown a priori. If the viscosity adaptive simulation is compared to the standard LBM simulations with M a = 0.01 and M a = 0.005, the VALBM simulation is more accurate as shown above and accelerated in terms of time steps by 178.2% and 336.4%, respectively.

12

Δp

1000

Mach number

0.025 0.02 0.015 0.01 0.005

Error

-1000

10−5 10−7 10−9

0.16 0.12 0.08 Ma α

0.04 0

0

1/4T

1/2T Period

3/4T

Womersley Number

0

T

Figure 5: Top: Pressure oscillation, Middle: Development of the Mach and Womersley number, Bottom: Development of the global absolute error over one period of oscillation.

Global absolute error Et

100 Standard LBM Ma=0.1 Ma=0.05 Ma=0.01 Ma=0.005 M aC = 0.0238591 VALBM

10−1

10−2

10−3

10−4

106

107 Number of time steps

Figure 6: Plot of the error measured over one period against the number of time steps for the standard LBM with different simulation Mach numbers and VALBM.

13

5. Conclusion For practical simulations of non-Newtonian fluids the correct simulation Mach number for which the modeled viscosity range describes the physical system is unknown a priori. The modeled range must be truncated because of stability criteria and the resulting viscosity bounds are chosen arbitrarily from a physical point of view. Therefore, longsome precursor simulations are necessary to determine a suitable simulation Mach number. In this paper we have derived a viscosity adaption method for Lattice Boltzmann simulations which makes it possible to consistently fit the modeled viscosity range to the actual physical problem. The viscosity adaption method ensures validity of the modeled viscosity range even if the physical problem changes in a transient case. We verified the proposed viscosity adaptive Lattice Boltzmann Method for a steady state channel flow and a transient pulsatile pipe flow of non-Newtonian fluids. As we have demonstrated in the verification cases that an improper choice of the simulation Mach number leads to corrupt results. For the steady state case the convergence properties and robustness of the VALBM regarding the starting simulation Mach number was investigated. Moreover, it was shown that the viscosity adaptive LBM ensures diffusive scaling for the target collision frequency Ωt = 1.0, which results in second order spatial convergence. In the transient flow situation it could be numerically verified that the simulation Mach number is adapted to the varying shear rates in the pulsatile flow. The results show that the viscosity adaptive LBM possesses superior accuracy to the standard LBM simulations. Moreover, VALBM could remain a larger mean time step size which results in a strong acceleration of the transient simulation. Although the viscosity adaptive LBM was employed here in a SRT manner, it is naturally incorporated in MRT approaches with virtually no extra cost, because the rescale process makes use of the same transformation matrix. Acknowledgements The simulations were partly conducted on the high performance cluster ”Elwetritsch” at TU Kaiserslautern which is part of the ”Alliance of High Performance Computing RheinlandPfalz” (AHRP). We kindly acknowledge for the support.

14

Appendix A. Transformation Matrix ⎡

(1 1 1 1 1 1 ⎢(−1 0 0 0 0 0 ⎢ ⎢ (1 −2 −2 −2 −2 −2 ⎢ ⎢ (0 1 −1 0 0 0 ⎢ ⎢ (0 −2 2 0 0 0 ⎢ ⎢ (0 0 0 1 −1 0 ⎢ ⎢ (0 0 0 −2 2 0 ⎢ ⎢ (0 0 0 0 0 1 ⎢ ⎢ (0 0 0 0 0 −2 ⎢ ⎢ (0 2 2 −1 −1 −1 ⎢ ⎢ (0 −2 −2 1 1 1 ⎢ ⎢ (0 0 0 1 1 −1 ⎢ ⎢ (0 0 0 −1 −1 1 ⎢ ⎢ (0 0 0 0 0 0 ⎢ ⎢ (0 0 0 0 0 0 ⎢ ⎢ (0 0 0 0 0 0 ⎢ ⎢ (0 0 0 0 0 0 ⎢ ⎣ (0 0 0 0 0 0 (0 0 0 0 0 0

1 1 0 1 −2 1 0 1 0 1 0 1 0 1 −1 0 2 0 −1 1 1 1 −1 1 1 1 0 1 0 0 0 0 0 1 0 −1 0 0

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 −1 1 −1 1 −1 −1 1 −1 1 −1 −1 −1 1 0 0 −1 −1 1 0 0 0 0 0 1 −1 0 0 0 1 −1 1 1 1 1 1 1 1 1 1 1 1 1 1 −1 −1 1 1 1 −1 −1 1 −1 −1 0 0 0 0 0 0 0 0 0 0 1 1 −1 1 −1 −1 1 1 1 −1 0 0 0 0 0 1 −1

1 1 1 1 1 0 0 −1 −1 1 1 −1 −1 0 0 −1 −1 0 −1

⎤ 1 1 1 1 1) 1 1 1 1 1 1) e2 ⎥ ⎥ 1 1 1 1 1) e4 ⎥ ⎥ −1 0 0 0 0) e ⎥ ⎥ −1 0 0 0 0) e3 ⎥ ⎥ 0 1 −1 1 −1) e ⎥ ⎥ 0 1 −1 1 −1) e3 ⎥ ⎥ 1 1 −1 −1 1) e ⎥ ⎥ 1 1 −1 −1 1) e3 ⎥ ⎥ 1 −2 −2 −2 −2) e2 ⎥ ⎥ 1 −2 −2 −2 −2) e4 ⎥ ⎥ −1 0 0 0 0) e2 ⎥ ⎥ −1 0 0 0 0) e4 ⎥ ⎥ 0 0 0 0 0) e2 ⎥ ⎥ 0 1 1 −1 −1) e2 ⎥ ⎥ −1 0 0 0 0) e2 ⎥ ⎥ 1 0 0 0 0) e3 ⎥ ⎥ 0 1 −1 1 −1) e3 ⎦ 1 −1 1 1 −1) e3

Matrix M for the linear transformation from distribution space to moment space.

References [1] J. Boyd, J. Buick, S. Green, A second-order accurate lattice Boltzmann non-Newtonian flow model, Journal of physics A: Mathematical and General 39 (46) (2006) 14241. [2] S. Gabbanelli, G. Drazer, J. Koplik, Lattice Boltzmann method for non-Newtonian (power-law) fluids, Phys. Rev. E 72 (2005) 046312. [3] R. Ouared, B. Chopard, Lattice Boltzmann simulations of blood flow: Non-Newtonian rheology and clotting processes, Journal of Statistical Physics 121 (2005) 209–221. [4] S. Succi, The lattice Boltzmann equation: For fluid dynamics and beyond, Oxford university press, 2001. [5] P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collison processes in gases, Phys. Rev. 94 (1954) 511. [6] T. Kr¨ uger, F. Varnik, D. Raabe, Shear stress in lattice Boltzmann simulations, Physical Review E 79 (4) (2009) 046704. [7] D. d’Humi`eres, I. Ginzburg, M. Krafczyk, P. Lallemand, L.-S. Luo, Multiple–relaxation–time lattice Boltzmann models in three dimensions, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 360 (1792) (2002) 437–451. [8] J. Harris, Rheology and non-Newtonian flow, Longman London, 1977. [9] T. Kr¨ uger, F. Varnik, D. Raabe, Second-order convergence of the deviatoric stress tensor in the standard Bhatnagar-Gross-Krook lattice Boltzmann method, Physical Review E 82 (2) (2010) 025701. [10] D. J. Holdych, D. R. Noble, J. G. Georgiadis, R. O. Buckius, Truncation error analysis of lattice Boltzmann methods, Journal of Computational Physics 193 (2) (2004) 595–619.

15

[11] J. H. Ferziger, M. Peri´c, Computational methods for fluid dynamics, Vol. 3, Springer Berlin, 1996. [12] P. J. Dellar, Incompressible limits of lattice Boltzmann equations using multiple relaxation times, Journal of Computational Physics 190 (2) (2003) 351–370. [13] J. Latt, Hydrodynamic limit of lattice Boltzmann equations, Ph.D. thesis, Universit´e de Gen´eve (03/09 2007). [14] Z. Guo, C. Zheng, S. Baochang, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys Rev E 65 (2002) 046308. [15] R. B. Bird, W. E. Stewart, E. N. Lightfoot, Transport phenomena, John Wiley & Sons, 2007. [16] M. Junk, A. Klar, L.-S. Luo, Asymptotic analysis of the lattice Boltzmann equation, Journal of Computational Physics 210 (2) (2005) 676–704. [17] R. Mei, W. Shyy, D. Yu, L.-S. Luo, Lattice Boltzmann method for 3-D flows with curved boundary, Journal of Computational Physics 161 (2) (2000) 680–699. [18] C. Loudon, A. Tordesillas, The use of the dimensionless Womersley number to characterize the unsteady nature of internal flow, Journal of theoretical biology 191 (1) (1998) 63–78.

16