The chaotic rotation of elongated particles in anoscillating Couette ow

(1)

Master Thesis

The chaotic rotation of elongated particles in an

oscillating Couette flow

Kungliga Tekniska H¨ogskolan

Author:

Niclas Berg

niber@kth.se Supervisors: Fredrik Lundell fredrik@mech.kth.se Tomas Ros´en rosen@mech.kth.se October 1, 2013

(2)

Abstract

Flows with suspended particles occur in a large range of engineering appli-cations and biological systems. A fundamental understanding of the fluid-particle coupling can be obtained by analyzing single fluid-particles in simple flow situations. Nilsen and Anderson (2013) studied the dynamics of a pro-late spheroid rotating in the flow-gradient plane of an oscillating creeping shear flow. They found that given a high enough particle inertia the ro-tation could become chaotic. This thesis extends on this work by taking fluid inertia into account and studying the rotation of an elliptical cylinder. Numerical simulations with the Lattice Boltzmann EBF method reveal a large variety of rotational states. Further, chaotic solutions are identified with a new Lyapunov exponent based algorithm. These chaotic solutions are found to be present for heavy particles (compared to the fluid) and small oscillation frequencies (compared to the shear rate in the channel).

(3)

Acknowledgements

First and foremost, I would like to thank my supervisors Tomas Ros´en and Fredrik Lundell for their never-ending stream of enthusiasm and for the knowledge they have passed on. I would further like to thank everyone at the fluid mechanics lab, in particular those in the Lundell group, for making my time as a thesis student very pleasant. Last, but not least, I would like to thank my family for their love and support.

(4)

Introduction

1.1 Background

Flows with suspended particles occur naturally in a large array of engineer-ing applications and biological systems. Consider for instance the process of paper making. In this process a mixture of wood fibers and water is emitted through a contraction onto a long screen where the water is dried off, leaving a web of interweaved fibers [1]. The quality of the finished paper is greatly influenced by how of the fibers are oriented with respect to each other.

Examples from the biological world include the blood circulation and the respiratory system where blood cells or dust represents the particles in the flow.

In all of the above examples, the particles and the flow will interact, which in turn will have an effect on e.g. the rheology of the flow. A funda-mental understanding of the coupling can be gained by analyzing the case of a single particle in simple flow situations.

The nature of the fluid-particle interaction depends on the inertia of the fluid and the particle, whose relative influences are quantified by the particle Reynolds number and the Stokes number, respectively. These numbers are given by

Rep =

Ga2

ν , St = αRep, (1.1) where G is the magnitude of the shear rate close to the particle, a the particle length scale, ν the kinematic viscosity of the fluid and α the solid-to-fluid density ratio.

The governing equations simplify greatly when fluid inertia is negligible (Rep 1). Jeffery [2] considered this case for an ellipsoid in a linear shear

flow, and derived an analytical expression for the torque on the particle. He further found that if the effect of both fluid and particle inertia are small (St 1, Rep 1), a spheroid (axisymmetric ellipsoid) will rotate

(7)

CHAPTER 1. INTRODUCTION

initial orientation and aspect ratio. These orbits are known as Jeffery orbits and resemble the motion of a kayak paddle.

Lundell and Carlsson [3] coupled Jeffery’s expression for the torque with the angular momentum equation, and evaluated numerically the effect of particle inertia on the orbits for a prolate spheroid (”a dragged out sphere”) in a linear shear flow. For low but finite values of St, the orbits resembled the ”kayaking” state but with a constant drift toward the plane in which the shear acts (the ”flow-gradient plane”). In these parameter regions, the rotation will thus eventually end up in this plane. For low values of Rep, it as

been seen that particle inertial effects dominate even for neutrally bouyant particles (α = 1) leading to a similar drift [4].

This drift to a planar rotation motivates the study of a 2D particle since the dynamics are equivalent, a fact already explored by Ding and Aidun [5] and Ros´en et al. [4]. Ding and Aidun found also that the influence of fluid inertia increases the rotation period which diverges to infinity past a critical Reynolds number, Rec. The particle thus enters a steady (motionless) state.

Particle inertia can also lead to a transition in rotation period in the fluid inertia-free case. Lundell and Carlsson, and later also Nilsen and Andersson [6], considered the case of a prolate spheroid with a rotation constrained to the flow-gradient plane. They both found that the period remained rather constant with St for low and high values of St, and with a rapid transition in the neighborhood of a critical value depending on the aspect ratio.

Nilsen and Andersson further considered the rotation in an oscillating linear shear flow. They found that the rotation could become chaotic given a high enough Stokes number, which was always greater than the critical value in the constant shear case.

For the chaotic rotation, small changes in parameters can yield a large change in the long-term statistics such as the average rotational energy. It is therefore possible that the transition to chaos can have an influence on the rheology. This effect can be further investigated but will not be part of this thesis.

It has already been seen that particle inertia affects the dynamics both in constant and oscillatory shear. However, although fluid inertial effects are well explored in constant shear, the effects in oscillating shear is still unknown.

In this thesis the effect of fluid inertia on the transition to chaos will be examined. The rotation will be restricted to the flow-gradient plane, and the particle is represented by an infinitely long elliptical cylinder. This reduces the spatial dimension of the problem to 2D which leads to a substantial reduction in computational cost, but the dynamics should still be equivalent to a prolate spheroid.

(8)

CHAPTER 1. INTRODUCTION

1.2 Structure of the thesis

In chapter 2, the physical model will be introduced and analyzed. The flow-particle solver is then introduced in chapter 3. In order to classify chaotic solutions some concepts from non-linear dynamics, in particular the Lyapunov spectrum, is introduced in chapter 4. Chapter 5 is devoted to the calculation of the Lyapunov spectrum. Some of the methods used in this thesis are evaluated in chapter 6. In chapter 7 and 8 results are presented and conclusions are then drawn in chapter 9.

(9)

Chapter 2

Model and preliminary

analysis

2.1 Physical model

The setup considered in this thesis is shown in figure 2.1 and consists of an elliptical cylinder placed in the center of an infinite channel containing a Newtonian fluid. The two channel walls are separated by a distance 2d and are moving in opposing directions with the time dependent velocity Uwg(t), where g(t) is non-dimensional. The main focus in this thesis will

be the case when the walls are oscillating with a constant frequency, f , and g(t) = sin(f t).

The wall-induced flow is assumed to be incompressible and can thus be described by the incompressible Navier-Stokes equation with accompanying continuity equation, ∂u ∂t + (u · ∇)u = − 1 ρf ∇p + ν∇2u, x ∈ Ωf, (2.1) ∇ · u = 0, x ∈ Ωf, (2.2)

where ρf is the fluid density, ν the kinematic viscosity, Ωf the region

oc-cupied by the fluid, and u = (u, v) and p the velocity and pressure fields, respectively.

No slip conditions are imposed on all boundaries, whence

u = Uwg(t)ex, y = d (2.3)

u = −Uwg(t)ex, y = −d (2.4)

u = ˙θez× x, x ∈ Γp (2.5)

where ˙θ is the angular velocity of the particle. Note that it has been assumed that the motion of the particle is purely rotational, since it is placed in the middle of the domain and will thus not experience any net force.

(10)

CHAPTER 2. MODEL AND PRELIMINARY ANALYSIS U_wg(t) Uwg(t) 2a 2b x x’ y y’ θ 2d

Figure 2.1: Problem setup

The particle is further assumed to be rigid and its surface, Γp, is thus

given by

x02 a2 +

y02

b2 = 1, (2.6)

where the prime indicates that the variable is expressed in the particle fixed coordinate system.

The surface of the particle will experience a shear stress from the fluid,

σ = −pI + µ(∇u + ∇uT) (2.7)

where I is the identity tensor and ∇u is to be interpreted as the Jacobian of u. The resulting force per unit length on a small line segment dΓ with outward pointing normal n, will be

df0 = σ0· ndΓ0, (2.8)

which, integrated over the circumference of the particle, yields the torque per unit length,

T = Z Γp x0× df0 = Z Γp x0× (σ0· n)dΓ0 = µ Z Γp x0× ((∇0u0+ ∇0u0T) · n)dΓ0. (2.9) The pressure only acts in the normal direction and does not contribute to the torque.

The resulting rotation of the particle is governed by Euler’s equation of motion

(11)

CHAPTER 2. MODEL AND PRELIMINARY ANALYSIS

where Izz is the moment of inertia per unit length along the z direction,

which for an elliptical cylinder equals to

Izz =

1

4πρsab(a

2_{+ b}2_), _(2.11)

where ρs is the density of the particle.

2.2 Scale analysis and non-dimensionalization

The governing equations can be cast into a more revealing form by rescaling all quantities so that all terms in the equations are of order 1. There does exist two length scales imposed by the geometry, i.e. the particle and the channel size, and the dominating physical effects may be different on these two scales. Let κ denote the scale separation, or the confinement of the particle, such that

κ = a

d. (2.12)

where the particle major axis half-length and the channel half-height has been used as scales.

The physics on the channel scale will be considered in the next section, and the discussion will now be focused on the particle scale.

The torque on the particle is, as indicated by eqn. (2.9), determined by the shear rate in the neighboring flow. Let G ∼ k∇uk be the scale of the shear rate near the particle. G has the dimension of time−1, which indicates that a suitable timescale is 1/G.

The scale of the velocity near the particle, Up, is revealed by eqn. (2.2),

which shows that

Up ∼ kuk = k ˙θez× xk ∼ aG. (2.13)

The pressure forces close to the particle are assumed to be of the same order of magnitude as the viscous, hence

k 1 ρf

∇pk ∼ kν∇2_uk _(2.14)

and p ∼ Gµ, where µ = ρfν is the dynamic viscosity.

The equations can now be non-dimensionalized through the transforma-tion

x ← ax, t ← Gt, u ← Gau, p ← Gµp, (2.15) which yields the set of equations

Rep ∂u ∂t + (u · ∇)u = −∇p + ∇2u, (2.16) ∇ · u = 0, (2.17) Stq ¨θ = Z Γp x0× ((∇0u0+ ∇0u0T) · n)dΓ, (2.18)

(12)

where Rep = Ga2/ν is the particle Reynolds number, q = Izz/ρsa4 the

non-dimensional inertia and St = αRep the Stokes number. Here α = ρs/ρf is

the solid-to-fluid density ratio.

The non-dimensional inertia for the elliptical cylinder is given by

q = 1

4πkb(1 + k

2

b), (2.19)

where kb = b/a is the aspect ratio. Throughout this thesis, this number as

well as the elongation,

= 1 − k

2 b

1 + k_b2, (2.20) will both be used to describe the particle shape.

The system (2.16)-(2.18) will be referred to as the FI (fluid inertial) system.

2.3 Flow without particle

Some understanding of the flow field can be gained by analyzing the case when no particle is present. Assuming solutions of the form

u = u(t, y)ex, p = p(y) (2.21)

gives the (dimensional) governing equation

∂u ∂t = ν ∂2u ∂y2, (2.22) u(d) = Uwg(t), (2.23) u(−d) = −Uwg(t), . (2.24)

The physical parameters on this scale can be combined to form a channel Reynolds number

Red=

Uwd

ν . (2.25)

2.3.1 Steadily moving walls

When the walls are moving with constant velocity, i.e. g(t) = 1, the solution is given by

u(t, y) = Uw

y

d. (2.26)

which corresponds to a linear velocity profile, with a constant shear rate, γsteady = Uw/d. Setting G = γsteady gives the following particle Reynolds

number: Rep = Ga2 ν = a2 d2 Uwd ν = κ 2_Re d (2.27)

(13)

CHAPTER 2. MODEL AND PRELIMINARY ANALYSIS −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 u/Uw y/d (a) β = 0.1 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 u/U w y/d (b) β = 2 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 u/U w y/d (c) β = 10

Figure 2.2: Velocity field in an oscillating Couette flow for different values of β and at times f t = 0, π/8, ..., π/2. As β is increased, a shear layer is formed close to the walls.

2.3.2 Oscillating walls

There exists an analytic solution [7] also for the case when the walls are oscillating. Let g(t) = sin(f t) for some frequency f . The system (2.22)-(2.24) then has the time-periodic solution

u(t, y) = Uw· Im sinh((1 + i)β(y/d)) sinh((1 + i)β) e if t , (2.28) where β =pfd2_/(2ν).

The effect of β can be seen in figure 2.2. For low values the flow has an almost linear profile. As β increases, a shear layer is formed near the walls in which a substantial momentum dissipation takes place. The result is a reduction in magnitude of the velocity and the shear rate near the middle of the channel.

The shear rate can be found explicitly by differentiating the velocity field

γ(t, y) = ∂u ∂y =

Uw

d · Im

cosh((1 + i)β(y/d))

sinh((1 + i)β) β(1 + i)e

if t

. (2.29)

A typical shear rate can be chosen as the maximum value during an oscillation period in the center of the channel, giving the particle Reynolds number Reosc._p = κ2Red max f t∈[0,2π]Im β(1 + i)eif t sinh((1 + i)β) . (2.30)

In the rest of this thesis, the particle Reynolds number based on the linear Couette flow will be used since it allows the oscillatory and the constant shear flow to be compared in the same non-dimensional framework. This Rep can be related to the oscillating case by

Reosc._p = Rep max f t∈[0,2π]Im β(1 + i)eif t sinh((1 + i)β) . (2.31)

(14)

The scaling function can be found plotted in figure 2.3. For large values of β, it is clear that the particle Reynolds number overestimates the ”actual” Reynolds number near the particle.

The parameter β can be rewritten in non-dimensional quantities as β =p(f /G)Rep/(2κ2). In order to be able to simulate high values of Reosc.p ,

the frequency needs to be kept low or the confinement high.

0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 β Re p osc. / Re p

Figure 2.3: The quotient between the particle Reynolds number in the os-cillating and the steady case.

2.4 The creeping flow limit

A substantial simplification of the governing equations can be made when the particle Reynolds number is small (creeping flow). In this case the non-linear term in the Navier-Stokes equation drops out, and the flow is governed by

∇p = ∇2_u. _(2.32)

Jeffery [2] analyzed this case for a linear shear flow around an ellipsoid, and derived an expression for the torque. In the notation of Lundell & Carlsson [3], but with a small correction, the non-dimensional form reads

T = 16π 3   (K2k−1_b + K3k−1c )−1[(kb2− kc2)D02,3+ (k2b + k2c)(V2,30 − ω01)] (K3k−1c + K1)−1[(k2c− 1)D03,1+ (k2c+ 1)(V3,10 − ω02)] (K1+ K2k−1_b )−1[(1 − k_b2)D01,2+ (1 + k2b)(V 0 1,2− ω03)]  , (2.33) where kb and kcare the aspect ratios of the ellipsoid, (ω10, ω02, ω30) the angular

velocities and D0 and −V0 the symmetric and anti-symmetric part of the velocity gradient tensor ∇0u0 of the base flow expressed in the body-fixed coordinate system. The Ki-coefficients can be evaluated by K1 = K(k2_b, kc2),

(15)

CHAPTER 2. MODEL AND PRELIMINARY ANALYSIS K2= K(k−2_b , k2ck−2b ) and K3 = K(k −2 c , kb2kc−2), with K(C1, C2) = Z ∞ 0 dλ (1 + λ){(1 + λ)(C1+ λ)(C2+ λ)}1/2 . (2.34)

The expression for the torque will now be used to derive the equation of motion for a particle rotating in the flow gradient plane (the x-y-plane) under the creeping flow assumption. The derivation will first be done for an ellipsoid, which will then translate to the elliptical cylinder case.

2.4.1 Rotation of an ellipsoid in the flow gradient plane

When the particle is rotating in the x-y-plane, the ω1 and ω2 components

are zero. Further, the transformation between the body and the world fixed coordinate systems is given by the rotation matrix

R =   cos θ − sin θ 0 sin θ cos θ 0 0 0 1   (2.35) such that x0 = Rx.

The velocity gradient tensor in the world fixed coordinate system

∇u =   0 g(t) 0 0 0 0 0 0 0   (2.36) hence transforms to ∇0u0= R(∇u)RT = g(t)  

− sin(θ) cos(θ) cos2(θ) 0 − sin2_(θ) _{cos(θ) sin(θ)} ₀

0 0 0



. (2.37)

Extracting the symmetric and anti-symmetric part of ∇0u0, and inserting the result into (2.33) yields

T = 16π 3   0 0 (K1+ K2k_b−1)−1[(k_b2− 1) sin2(θ)g(t) − k_b2g(t) − (1 + k2_b) ˙θ]   (2.38) The z-component can now be inserted into the Euler equation

Stq ¨θ = 16π 3 (K1+ K2k −1 b ) −1_[(k2 b− 1) sin2(θ)g(t) − kb2g(t) − (1 + kb2) ˙θ] (2.39)

and with the use of the non-dimensional moment of inertia of an ellipsoid,

q = 4

15ρsπkbkc(1 + k

2

(16)

one finally obtains

Stkc(kbK1+ K2) 20 ¨ θ = 1 − k 2 b 1 + k2_b 1 2− sin 2_θ g(t) −1 2g(t) − ˙θ. (2.41)

This can be restated with the use of some algebra to

f

St¨θ = −1

2[1 − cos(2θ)]g(t) − ˙θ. (2.42) where the elongation = (1 − k2_b)/(1 + k2_b) and the modified Stokes number f

St = Stkc(kbKα+ Kβ)/20 have been introduced as in ref. [6]. This system

will be referred to as the CF (creeping flow) system.

There are symmetries in (2.42) that should be noted. If θ(t) is a solution so will θ(t) + π. Further, if g(t) = sin(f t) or g(t) = cos(f t) then −θ(t + π/f ) will also be a solution. Hence, if the system has a periodic solution for which −θ(t + π/f ) 6= θ(t) then there will exist at least two different possible periodic solutions.

2.4.2 Rotation of an elliptic cylinder in the flow gradient

plane

By letting kc go to infinity, the equation of motion corresponds to that of

an elliptical cylinder. The integrals in the modified Stokes number were evaluate numerically using the quadgk method in Matlab. In the order to evaluate the limit of kc → ∞ the integrals were evaluated repeatedly

with the value of kc successively doubled. In this way, and given that the

integrands decrease with increasing kc, a converged value could be obtained.

The convergence criteria used was that the difference between two successive values should be less than 10−6.

The numerical evaluation showed that

lim kc→∞ kc(kbK1+ K2) 20 = kb 10, (2.43)

(17)

Chapter 3

The Lattice Boltzmann EBF

method

In this chapter, different techniques of simulating flows will first be discussed. Then a description of the Lattice Boltzmann method follows, which is the method of choice for this project. The implementation details for the current case are then provided in the last section.

3.1 From micro to macroscale methods

There exists a large array of methods for simulating flows. When classifying different methods, one starting point can be the length and times scales they consider. On the lowest level one could consider the motion of each molecule in the flow, which is known as molecular dynamics simulations (MD). In these simulations a potential function expressing the inter-molecular forces is assigned and Newton’s laws of motion are solved for each molecule. New physical properties can be introduced quite easily by adding extra terms to the potential.

The downside of MD is the shear amount of information that needs to be tracked. In an ordinary flow the number of molecules is at least of the order of the Avogadro number (NA∼ 1023). The computational complexity

that this induces is of course a large challenge, as well as the question of how one should set initial conditions. In addition, does one always need to know the motion of each molecule in a flow?

In most scientific and engineering applications the usual quantities of interrest are, for instance, shear stresses on bodies in the flow, the mean motion of the fluid and the pressure field. These effects are due to the collective motion of molecules and are given by time and space averages, taken over suffiently small domains. Up until this point in this thesis, all of the physical modelling have indeed been discussed in terms of macroscopic quantities.

(18)

CHAPTER 3. THE LATTICE BOLTZMANN EBF METHOD

In conventional computational fluid dynamics (CFD) the starting point is the macroscopic conservation laws for mass, momentum and energy (e.g. the incompressible Navier-Stokes equation). These equations are then dis-cretized in space and time using resulting in (large) algebraic equation sys-tems to be solved.

In recent years a new approach has been developed to simulating flows, called the Lattice Boltzmann method. The method has its origins in kinetic theory of gases but attempts to discard as much information as possible while making sure that the correct macroscopic description of the system is retained.

3.2 Kinetic theory

A fluid is composed of a large amount of molecules. Therefore it appears reasonable to consider the statistics of the molecular motion. Let f (x, c, t) denote the probability of finding a particle at position x moving with velocity c at time t. Knowing f , macroscopic quantities such as density, momentum and internal energy can be obtained through the velocity moments

ρ(x, t) = Z mf (x, c, t)dc, (3.1) ρ(x, t)u(x, t) = Z mcf (x, c, t)dc, (3.2) ρ(x, t)e(x, t) = = Z mc2f (x, c, t)dc, (3.3)

where m is the mass of one molecule. These expressions constitute the link between the macroscopic and mesoscopic description of a flow.

3.2.1 The Boltzmann transport equation

Ludvig Boltzmann derived in 1872 an equation that bears the name of the Boltzmann equation and describes the time evolution of the distribution function of a gas. In absence of external forces, the statistics will follow

∂

∂t+ c · ∇x

f = Ω[f ], (3.4)

where Ω[f ] is a functional representing the local redistibution in momentum due to the collision of molecules.

3.2.2 The Bhatnagar-Gross-Krook collision model

The general expression for the collision operator, Ω[f ], is quite complex and is not suitable for computational purposes. In order to circumvent

(19)

this, several models have been developed to describe the collision process. One commonly used model is the one derived by Bhatnagar, Gross and Krook [8] which states that collisions cause the local state to relax towards an equilibrium distribution at a single timescale. The equilibrium distribution is the state for which Ω[f ] = 0, and is given by

feq = ρ (2πRT )D/2 exp −(c − u) 2 2RT , (3.5)

where T is the temperature, R the ideal gas constant and D the dimension of the space.

The BGK collision operator is defined as

Ω[f (x, c, t)] = −1

τ[f (x, c, t) − f

eq_{(x, c)],} _(3.6)

where τ is the relaxation time.

3.3 The Lattice Boltzmann method

The Lattice Boltzmann method (LBM) was developed during the late 1980s and is a rather new method for simulating hydrodynamic problems. It has its origins in the Lattice Gas Automata [9] and was originally proposed as a way to solve statistical noise issues therein. It has later been proven that the method can be obtained as a second order discretization of the Boltzmann equation [10].

In LBM the distribution function is represented in a discrete configura-tion space. Space is divided into a regular lattice with posiconfigura-tions represented by integers, so that x = (i, j, k) with i ∈ [0, Nx − 1], j ∈ [0, Ny − 1] and

k ∈ [0, Nz− 1]. Each lattice point (or node) is linked to its neighbors, and

these linkages define a set of discrete velocities eα (i = 0, 1, ..., q − 1). A

node at x thus has a neighbor at x + eα.

The discrete distribution function (the ”populations”) evolves over time in the configuration space in two steps. First, the collisions between molecules cause a local redistribution in momentum (fig. 3.1b) according to

f_k∗(x, t + 1) = Ωk[fk(x, t)] (3.7)

where Ωk is a discrete collision operator. For the purpose of this thesis, a

discrete BGK-like model will be used.

The post-collision populations then propagate to neighboring nodes (fig. 3.1c),

fk(x + ek, t + 1) = fk∗(x, t). (3.8)

This step is known as the streaming step and represents the movement of the molecules in the fluid.

(20)

CHAPTER 3. THE LATTICE BOLTZMANN EBF METHOD f1 f5 f2 f6 f3 f7 f4 f8

(a) Before (b) Collision (c) Streaming

Figure 3.1: A Lattice Boltzmann iteration. The length of the arrow indicate the size of the populations at the node and the directions the associated velocities.

A complete Lattice Boltzmann BGK iteration can be written as

fk(x + ek, t + 1) = fk(x, t) −

1

τ(fk− f

eq

k ), (3.9)

where τ is the relaxation time which is, as we will see later, related to the viscosity of the fluid and f_keq is a discrete equilibrium distribution.

The macroscopic variables can, as in the continuous case, be obtained from the velocity moments of the distribution function,

ρ =X k fk, (3.10) ρu =X k fkek, (3.11) Π = ∇(ρu) + ∇(ρu)T =X α eαeα(fα− fαeq), (3.12)

where eαeα is to be interpreted as the outer product between the vectors.

The hydrostatic pressure can be found by the equation of state

p = c2_sρ (3.13)

where cs is the speed of sound. For incompressible flows this is a constant

that depends on the lattice.

It has been shown [11,12] that the Lattice Boltzmann BGK scheme yields the weakly compressible Navier-Stokes equation in the low Knudsen and Mach number limit, given that an appropriate lattice and equilibrium distribution function is provided. The lattice used in this thesis will be described in the next section.

Before concluding this section, it should be stated that the dimension of the quantities yielded by (3.10)-(3.12) will not be laboratory units. LBM

(21)

CHAPTER 3. THE LATTICE BOLTZMANN EBF METHOD e1 e2 e3 e4 e5 e6 e7 e8

Figure 3.2: The D2Q9 lattice.

employs its own unit system, known as lattice units, in which the distance between two neighboring nodes as well as the duration of a timestep are both unitary. Dimensionless parameters are however the same. Prior to using LBM, all quantities should thus be converted in a way so that the relevant non-dimensional numbers are obtained.

3.4 The D2Q9 lattice

The first step when using LBM is to choose a suitable lattice. This choice may at first seem arbitrary, but there are actually very few arrangements that can produce reliable results.

One of the main concerns is that the continuous equilibrium distribution must be able to be represented accurately in order for the method to be stable [13]. Furthermore, certain symmetries of the lattice is required in order for the model to approximate the Navier-Stokes equation [14]. A systematic way of constructing proper lattices is given in ref. [10].

Lattices are commonly classified according to the DnQm notation, where n is the dimension of the space and m is the number of discrete velocities. For incompressible hydrodynamic problems in 2d , the most commonly used lattice is the D2Q9 lattice which is shown in figure 3.2. The discrete veloci-ties are given by

ek =      0, k = 0 cos(π₂(k − 1)), sin(π₂(k − 1) , k = 1, 2, 3, 4 √ 2 cos(π₂(k − 1)), sin(π₂(k − 1) , k = 5, 6, 7, 8 . (3.14)

with a lattice speed of sound of cs = 1/

√ 3.

(22)

The equilibrium distribution is evaluated by

f_keq= wkρ(1 + 3ek· u + 9 2(ek· u) 2₋3 2u 2_). _(3.15) with weights wk=      4/9, k = 0 1/9, k = 1, 2, 3, 4 1/36, k = 5, 6, 7, 8 . (3.16)

3.5 Boundary conditions

In order for a hydrodynamic problem to be well posed, boundary condi-tions must be supplied. These are usually defined in terms of macroscopic quantities and must be related to the distribution functions in LBM.

At each boundary node, there will be populations whose associated ve-locity direction is pointing out of the domain. After each streaming step, these populations will be unknown. The aim of the boundary conditions in LBM is thus to find appropriate values for these populations, and possibly also replace the other ones, so that the wanted macroscopic state is obtained.

3.5.1 Periodic boundary condition

Periodic boundary conditions are rather straightforward to implement, and can be included into the streaming step. Assume for now that a periodic boundary condition is to be implemented in the x-direction. The popula-tions at the left and right boundaries are now neighbors. The left-going populations at the left wall (3, 6, 7) should thus propagate to the corre-sponding ones at the right wall, and vice versa.

This can be stated as:

f1(0, j) = f1(Nx− 1, j), (3.17) f5(0, j + 1) = f5(Nx− 1, j), (3.18) f8(0, j − 1) = f8(Nx− 1, j), (3.19) f3(Nx− 1, j) = f3(0, j), (3.20) f6(Nx− 1, j + 1) = f6(0, j), (3.21) f7(Nx− 1, j − 1) = f7(0, j). (3.22)

3.5.2 Velocity boundary condition

When dealing with boundaries in the incompressible Navier-Stokes equation only the velocity needs to prescribed in order for the problem to be well-posed. The pressure is a mere functional whose only purpose is to keep the velocity field divergence free. In contrast, the LBGK scheme models

(23)

the weakly compressible case and the density (and therefore the pressure) is a part of the system’s state. The boundary condition implementation must thus both find appropriate values for the density and the unknown distributions.

The key equations are (3.10) and (3.11) which, in 2d, form a set of 3 equations. A wall node in the D2Q9 lattice has 3 outgoing populations, and assigning values to these as well as the density is thus an underdetermined problem. In this section three implementations which addresses this issue will be discussed. For the sake of simplicity, it will be assumed that the boundary node is on the northern wall.

All the boundary conditions introduced in this section are supposed to be followed by a collision and streaming step in order to yield accurate results.

Evaluating the density at the wall

The way in which the density at the wall is evaluated is almost the same for all boundary conditions in LBM. Consider the density eqn. (3.10) and the velocity eqn. (3.11) in the wall normal direction

ρ =X k fk= f0+ f1+ f2+ f3+ f5+ f6+ f4+ f7+ f8, (3.23) ρv =X k fkek,2= f2+ f5+ f6− (f4+ f7+ f8). (3.24)

The unknown populations f4, f7 and f8 can be eliminated by taking the

sum of the above equations, and the density is thus given by

ρ = 1

1 + v[f0+ f1+ f3+ 2(f2+ f5+ f6)] . (3.25)

The bounce back condition

The bounce back condition is a way to model the interaction between a stationary wall and the flow (no slip). This is achieved by letting the in-coming populations, as the name of the condition implies, bounce back in the direction they came from.

Let opp(k) denote the opposite direction of the population k, such that eopp(k)= −ek. Then for all inwards going population, one sets

f_opp(k) = fk. (3.26)

It can be easily verified that this yields a zero velocity in (3.11).

The Zou and He condition

This boundary condition developed by Zou and He [15] and allows a non-zero velocity to be specified at the boundary. The bounce back rule is first

(24)

applied to the non-equilibrium part of the population in the wall normal direction,

f4− f₄eq= f2− f₂eq. (3.27)

Using the expression for the equilibrium distribution, one finds

f4 = f2−

2

3ρv (3.28)

The two remaining unknown populations can then be found with the aid of (3.10) and (3.11), f7 = f5+ 1 2(f1− f3) − 1 6ρv − 1 2ρu, (3.29) f8 = f6+ 1 2(f3− f1) + 1 2ρu − 1 6ρv. (3.30) This boundary condition yields a state at the boundary that has a con-sistent density and velocity.

The Regularized condition

The regularized condition introduced by Latt [12,14] replaces all the pop-ulations at the boundary. Non-equilibrium bounce back is first applied to assign values to all the unknown populations. The strain tensor, Π, at the node is then evaluated using eqn. (3.12). Using the strain tensor, the off-equilibrium components of the populations can be evaluated using the expression

f_kneq = wi 2c4 s

(ekek− c2sI) : Π, (3.31)

and the final replacement value at the node is thus

fk = f_keq+

wi

2c4 s

(ekek− c2sI) : Π. (3.32)

This boundary condition yields not only appropriate values for ρ and u but also for the stress tensor at the wall [12].

3.6 Body forces

The Lattice Boltzmann method can be extended to include body forces in the way proposed by Latt [12]. Let g denote the force density, the complete scheme then reads

fα(x + eα, t + 1) = fα(x, t) −

1

τ[fα(x, t) − f

eq

(25)

3.7 Implementation

3.7.1 Problem setup and initialization

The geometry used for simulating the FI system consists of a quadratic domain with Nx = Ny = 128 nodes, if not otherwise stated. A periodic

boundary condition is imposed in the x-direction (i = 0 and Nx− 1) and a

velocity condition at the lower and upper walls, with velocities U_wLBMg(t) and −U_wLBMg(t), respectively. A value of U_wLBM = 0.05 has been used.

The non-dimensional quantities of this problem can be stated in lattice units Rep = GLBMa2_LBM νLBM , (3.34) α = ρ LBM s ρLBM_f , (3.35) κ = 2aLBM Ny− 1 , (3.36) f = fLBM GLBM , (3.37) kb = bLBM aLBM , (3.38)

where GLBM = 2UwLBM/(Ny − 1) and ρf = 1. By providing the

non-dimensional numbers to the solver, all the parameters in lattice units can be obtained from the above relations.

Prior to the simulation, all the populations are initialized to the equi-librium distribution for the given velocity field u and with a density of ρLBM_f = 1.

3.7.2 Fluid-solid interaction

The interaction between the solid and the fluid will in this thesis be resolved using the External Boundary Force (EBF) method developed by Wu and Aidun [16]. In this method, the fluid is represented on an Eulerian grid and the solid on a Lagrangian and the no slip condition on the surface is modeled with a force acting on the fluid and solid.

Let Ωf and Ωs denote the solid and fluid domains, respectively, and

ΩM_f = {xf_i} and ΩM

s = {xsi} the discretizations thereof. The boundary

between the domains is denoted Γ, with discrete counterparts Γf and Γs.

Now let g(x, t) and Ff si(x, t) = −g(x, t) represent the force densities acting on the fluid and solid, respectively. These are to be chosen so that the velocity at Γ is kept the same for both the solid and fluid nodes. If Uf(xs_i, t) and Us(xs_i, t) are the fluid and solid velocities at a solid boundary

(26)

node, then prior to an iteration of the flow and solid solver it will, due to the no slip condition, hold that

Uf(xs_i, t − ∆t) = Us(xs, t − ∆t), (3.39)

where ∆t = 1 is the LBM timestep.

After the iteration there will be a difference between the velocities which is counteracted by imposing a force on the solid nodes given by

Ff si(xs_i, t) = ρf[U

f_(xs

i, t) − Us(xsi, t − ∆t)]

∆t (3.40)

where ρf is the density of the fluid.

The nodes on the two grids will in general not coincide and an interpo-lation must be performed to find the fluid velocity at the node. This is done with the use of a discrete delta function, D(x), so that

Uf(xs_i, t) = Z Ωf Uf(xf, t)δ(xf− xs_i)dxf ≈ X xf j∈ΩMf Uf(xf_j, t)D(xf_j − xs_i)∆Vf (3.41)

where ∆Vf = 1 is the volume (in 3D) or area (in 2D) occupied by a fluid node. For 2d problems, the discrete Dirac function is approximated by

D(x =    1 16h2 1 + cosπx 2h 1 + cosπy 2h if kxk ≤ 2h, 0 otherwise. (3.42)

where h = 1 is the spacing between the fluid nodes.

Once the force density has been obtained, the total force, F , and torque, T , on the solid can be found through

F (t) = X xi∈Γs Ff si(xs_i, t)∆V_is, (3.43) T (t) = X xi∈Γs (xs_i − xc) × Ff si(xsi, t)∆Vis, (3.44)

where xc the center of mass of the solid and ∆Vis is the volume (or area)

occupied by the ith solid node.

The torque can now be inserted into the equations of motion. In our case this corresponds to solving

(27)

The integration was performed with the explicit Euler method. If one lets ω = ˙θ, this can be expressed as

θ(t + ∆t) = θ(t) + ∆tω(t), (3.46)

ω(t + ∆t) = ω(t) + ∆t Jzz

Tz(t). (3.47)

The angle θ(T + ∆t) can now used to update the position of the solid nodes, and these new node positions are used to find the force on the fluid,

g(xf_i, t) = X xs

j∈ΩMs

D(xs_j − xf_i)Ff si(xs_j, t)∆V_js, (3.48)

which was then included in the LBM via eqn. (3.33).

The above algorithm was implemented in the C programming language and parallelized with OpenMP.

(28)

Chapter 4

Concepts from nonlinear

dynamics

4.1 Introduction

The theory of non-linear dynamics is concerned with systems on the form

˙

x = f (x; µ1, µ2, ...), (4.1)

where f depends (non-linearly) on x and on a set of parameters {µi}i. The

entity x is usually refered to as the state of the system and the set of all possible states is called the phase space.

If the state is known at a time t, an infitesimal time later the state will have evolved to x(t + dt) = x(t) + f (x(t))dt. The solutions to a dynamical system can thus be interpreted as streamlines of a flow field f (x). These lines are known as trajectories or orbits of the system.

A complete description of a system would, in theory, require knowledge of all the trajectories. This is very infeasible since it is only in rare cases that the system can be integrated analytically, and one must in general rely on numerical methods. Qualitative features of the long term behavior can however usually be extracted with the aid of geometric and physical reasoning.

A useful fact is that many real world systems are dissipative. These systems most commonly contain internal friction and in order for the system to not end up in rest, energy must be constantly provided through a driving force. The combination of the driving force and the dissipation tends to eliminate transients. If the system is allowed to evolve over a long time, all trajectories will end up in confined regions of the phase space, the limit set. This also implies that volumes in phase space contract over time.

The classification of a system usually boils down to tracing out the limit set of the system and analyzing the stability of its’ subsets. This greatly simplifies matters since the limit set often has lower dimension than the

(29)

CHAPTER 4. CONCEPTS FROM NONLINEAR DYNAMICS

phase space itself. The difficulty of finding different subsets of the limit set depends on the neighboring flow. Attractors are stable sets to which nearby trajectories are approaching whereas repellors are unstable and orbits in their neighborhood are repelled. There may also exist sets that act as attractors in some directions and repellors in others.

Before concluding this section it should be pointed out that the flow and the topology of the limit set will vary with the parameters. If a qualitative change takes places due to the change of a parameter, the system is said to undergo a bifurcation.

4.2 The Lyapunov spectrum

While classifying the limit sets of a system, techniques such as linear stability theory [17] and Floquet theory are commonly used for fixed points and periodic orbits, respectively. These methods provide growth rates of small disturbances applied to the set. In linear stability theory the growth rate is provided by eigenvalues of the Jacobian matrix while Floquet theory uses Floquet exponents. If all of the growth rates are negative, the set is stable. The Lyapunov spectrum provides a generalization of the above methods to arbitrarily time-varying systems. An introduction is provided by Wolf et al. [18], and only the key points will be stated here.

The spectrum can be understood as follows: For a continuous n-dimensional system, the main idea is to consider the action of the system on an infinites-imal n-sphere of initial conditions. Over time the trajectories will approach each other in some directions and possibly separate in others, thus deform-ing the volume into an n-ellipsoid. If kδxi(t)k is the length of one of the

principal axes of the hyperellipsoid at time t, the corresponding Lyapunov exponent (LE) is defined as

λi= lim t→∞ 1 t ln kδxi(t)k kδx_i(t0)k , i = 1, ..., n, (4.2)

where the Lyapunov exponents are ordered by magnitude, so that λ1 ≥ λ2 ≥

... ≥ λn.

The signs of the Lyapunov exponents provide a qualitative description of the attractors of the system. A positive (negative) exponent is associated with an axis along which orbits are on average diverging (converging) expo-nentially fast. For instance, if all exponents are negative, trajectories in all directions are converging exponentially from all directions toward a point. The system then has an asymptotically stable fixed point [19].

For a time-dependent dissipative three-dimensional systems the only pos-sible combinations of signs and associated attractors are [18]:

(0, -, -) Limit cycle (a closed loop) (0, 0, -) Two-torus, quasi-periodic orbit (+, 0, -) Strange attractor

(30)

The strange attractor is particularly interesting. On one hand, nearby orbits are diverging at an exponential rate in one direction near the attractor. At the same time, the system is volume contracting and the attractor must thus be bounded. In order for this to be possible the phase space near the attractor must have been repeatedly stretched and folded [17,18]. This decorrelates nearby orbits, making the attractor very sensitive to initial conditions. This is one of the properties of chaos, which will be discussed further in the next section.

As a final remark, it should be noted that the Lyapunov spectrum is a global property of the system and is independent of the choice of reference trajectory. This is a consequence of the Osedelets ergodicity theorem [20] and holds only in the limit of time tending toward infinity.

4.3 Chaos

There does not as of now exist a universally agreed on definition of chaos. In this thesis the one due to Strogatz [17] will be used, which tells that a chaotic system should have the following properties:

1. Aperiodic long-term behavior - There exists orbits that does not end up in fixed points, periodic orbits or quasiperiodic orbits as time tends to infinity.

2. Deterministic - The system does not contain any random elements.

3. Sensitive dependence on initial conditions - The system has at least one positive Lyapunov exponent.

The last property is one of the key features of chaos. To see the implica-tions, consider two orbits separated by a distance δx0 at time t0. Over time

this distance will grow as

δx(t) ≈ δx0eλ(t−t0). (4.3)

The exponential separation rate of nearby orbits introduces an unpre-dictability in the system, especially from a practical point of view. This is due to the fact that measuremens of real world systems are never exact. If δx0 represents the accuracy of the instruments used, the uncertainity of

the state at a later time will grow exponentially in time, in accordance with equation (4.3). The state of the system can thus only be predicted for a finite time horizon.

As an example, Nilsen et al. [6] found a maximal Lyapunov exponent for the creeping flow system, (2.42), of λmax = 0.053 for a certain parameter

combination. They concluded that an initial uncertainty of 10−6 could grow to as much as 0.1 after only 7.7 wall oscillation periods.

(31)

As one can see, in order to classify a chaotic system, only the knowl-edge of the largest Lyapunov exponent (LLE) is required. This will greatly simplify matters in the next chapter.

(32)

Chapter 5

Computation of the largest

Lyapunov exponent

In the last chapter, the importance of the Lyapunov spectrum and in par-ticular the largest exponent (the ”LLE”) was established. In this chapter several methods for its computation will be introduced and compared. Some useful properties will first be derived in the mathematical preliminary. This leads naturally to an algorithm for the creeping flow system, which is then extended to the FI system. Finally, an a posteriori method is introduced.

5.1 Mathematical preliminary

5.1.1 A formal definition

In the last chapter the Lyapunov spectrum was introduced as the long-term evolution of the principal axes of an infinitesimal sphere of initial conditions. A more formal definition can be given in terms of a generalized Jacobian matrix X, defined so that

δx(t) = X(t; t0)δx(t0), (5.1)

for any arbitrary, infinitesimal disturbance δx(t0) (i.e. a point on the surface

of the sphere).

If Λ(X(t; t0)TX(t; t0)) = {αi(t)}i, is the set of (time-dependent)

eigen-values of XTX, the ith Lyapunov can be defined as

λi= lim t→∞

1

2tln αi(t), i = 1, 2, ..., n. (5.2) To illustrate the link to the geometric introduction in the last chapter, let vi(t) be the normalized eigenvector corresponding to αi. It then holds

that

(33)

CHAPTER 5. COMPUTATION OF THE LARGEST LYAPUNOV EXPONENT

Taking the eigenvector as direction of the ith principal axis of the ellipsoid, vi(t) = δxi(t0)/kδxi(t0)k, gives αi = δxi(t0)T kδxi(t0)k (XTX) δxi(t0) kδxi(t0)k = kXδxi(t0)k 2 kδxi(t0)k2 (5.4) and thus λi = lim t→∞ 1 2tln kXδxi(t0)k2 kδxi(t0)k2 = lim t→∞ 1 t ln kδxi(t)k kδxi(t0)k . (5.5) 5.1.2 Disturbance convergence

Assume now that one wants to evaluate the Lyapunov spectrum for a system. This would require the knowledge of how orbits near a reference trajectory are behaving, which can be obtained, for instance, by solving the system for a number of slightly different initial conditions. Another way would be to linearize the equation around a reference trajectory and solve this system in parallel with the non-linear equation.

Now consider what would happen to an initial, infinitesimal perturbation in a unit direction n0. Expanding this in the eigenbasis {vi} of XTX gives

n0 =

X

i

(n0· vi)vi. (5.6)

The length of the perturbation will, at a later time t, have scaled to

kn(t)k2= X i (n0· vi)vTi ! (XTX)   X j (n0· vj)vj   = X i (n0· vi)vTi !  X j (n0· vj)αivj   =X i (δx0· vi)2αi = (n0· v1)2α1 1 +(n0· v2) 2 (n0· v1)2 α2 α1 + O α3 α1 ,

since α1 > α2≥ α3 ≥ ... ≥ αn for large times.

From the definition of the Lyapunov exponents it follows that αi(t) ≈

e2tλi _{and thus,} kn(t)k2 _{≈ (n} 0· v1)2e2tλ1 1 +(n0· v2) 2 (n0· v1)2 e−2t(λ1−λ2)_{+ O}_e−2t(λ1−λ3) . (5.7) All the terms in the braces goes to zero at an exponential rate and any infinitesimal perturbation that is not perpendicular to the fastest growing

(34)

eigendirection will thus align itself with it. This fact helps out greatly if only the largest exponent is sought, as will be shown in the next section. In the case when the complete spectrum is wanted, on the other hand, this presents a great challenge. Even if the directions of the initial perturbations have been chosen with care, numerical errors will most likely produce small components in the fastest growing direction, which will then dominate. This is usually handled by repeatedly reorthogonalizing the disturbances during the calculation [18], or by related methods based on matrix factorizations. A review of such methods is given by Geist et al. [21].

5.2 Finding the largest Lyapunov exponent

As was shown in the last paragraph, basically any small initial perturbation will evolve at the rate corresponding to the largest Lyapunov exponent. One can thus solve the non-linear equations for two initial conditions x0 and x00

separated by a small (but finite) distance δx0 and monitor how this distance

evolves over time.

The distance between the orbits can however not grow/shrink infinitely. An upper bound is imposed by the diameter of the attractor, and the lower by the limited machine accuracy. This can be solved by normalizing the distance regularly during the integration.

In order to present this idea, split the time interval [t0, t] into smaller

in-tervals, so that t0 < t1 < ... < tn= t. Given the definition of the generalized

Jacobian, it follows that

δx(tn) = X(tn; tn−1)X(tn−1; tn−2) . . . X(t1; t0)δx0, (5.8)

and hence δx(tk) = X(tk; tk−1)δx(tk−1).

Now let δx(tk) = δxk and introduce the following iteration scheme:

δx∗_k= X(tk; tk−1)δ ˜xk−1, (5.9) dk= kδx∗_kk kδx0k (5.10) δ ˜xk= δx∗_k dk (5.11)

where δx∗₀= δ ˜x0 = δx0. Note that kδ ˜x∗kk = kδx0k for all k.

Combining these expressions gives dkδ ˜xk = X(tk; tk−1)δ ˜xk−1. Insertion

into (5.8) and recursive evaluation then yields

δxn= X(tn; tn−1)X(tn−1; tn−2) . . . X(t2; t1)X(t1; t0)δx0

= X(tn; tn−1)X(tn−1; tn−2) . . . X(t2; t1)d1δ ˜x1

= ...

= dndn−1...d1δ ˜xn.

(35)

A finite time Lyapunov exponent can thus be found by

λmax = 1 tn− t0 lnkδx(tn)k kδx(t0k = 1 tn− t0 ln n Y k=1 dk ! = 1 tn− t0 n X k=1 ln dk, (5.13) which will converge to the largest Lyapunov exponent when time goes to infinity, by the arguments of the last section.

5.3 Implementation for the CF system

The CF system is given by ˙

x = f (x), (5.14)

where x = (θ, ω, t)T with initial condition x(t0) = x0 = (θ0, ω0, t0). The

flow function equals

f (x) =     ω 1 f St −1 2[1 − cos(2θ)]g(t) − ω 1     . (5.15)

An infinitesimal perturbation δx0 = (δθ, δω, δt) applied at t = t0, will

evolve in accordance with the linearized system ˙

δx = J (x)δx, (5.16) where J is the Jacobian of f , which is given by

J (x) =   0 1 0 − sin(2θ)g(t)/fSt −1/fSt −[1 − sin(2θ)]g0(t)/(2fSt) 0 0 0   (5.17)

The last row consists of zeros, and the magnitude of an initial perturba-tion in the t direcperturba-tion will stay constant over time. This corresponds to a zero-exponent, and it is easy to see that the associated direction is the one tangential with the flow.

The initial perturbation was chosen as δx0 = n0δx0, where the unit

direction n0 was chosen as

n0 = 1 pf1(x0)2+ f2(x0)2   f2(x0) −f1(x0) 0   (5.18)

i.e. in a direction orthogonal to the flow - the zero-exponent direction. The introduction of the unit direction allows rewriting (5.16) to

˙

n = J (x)n, (5.19)

(36)

which removes the requirement of specifying a small distance δx0- the

small-ness is granted by the fact that the evolution is governed by the linearized equations.

The eqns. (5.14) and (5.19) were integrated simultaneously with the fourth order Runge-Kutta-Fehlberg method implemented in the GSL library [22]. The initial condition for the reference orbit was generated from a zero state integrated until t = 2000fSt.

The perturbation direction was normalized after each integration. The normalization procedure used was equivalent to the one in the last section but with the growth factors evaluated as dk = kn(tk)k, and the length of

the vector normalized to 1.

At each time that was an integer multiple of the period, 2π/f , the Lya-punov exponent estimate was saved. The exponent was considered to have converged if the standard deviation of the last 20 values was less than a thousandth of the mean of the values.

The final value was averaged over 5 different initial perturbations evenly distributed in one period.

5.4 Extension to the FI system

For the FI system, the state is given by not only the state of the particle (θ, ω) but also by the velocity and pressure fields of the flow. These fields are continuous functions in space, rather than discrete values, and a suitable norm must be chosen to measure the divergence rate in phase space. One option would be to use the L2-norm for the fields and the Euclidean for the

scalar quantities, so that the measure associated with a state x = (θ, ω, u, p) would be given by

kxk2 = θ2+ ω2+ Z

Ωf

(kuk2₂+ |p|2)dΩ. (5.21)

where the integral is to be taken over the fluid domain.

Having defined a measure, the distance between two states x and x0 is given by kx − x0k. In order to evaluate the Lyapunov exponent, two initially slightly perturbed states would have to be considered. Over time, the fluid domains for the two solutions would however not in general coincide and there is no obvious way to evaluate the integral.

To circumvent this, an algorithm was developed in which the divergence rate was measured only in the particle variables (θ, ω). The flow was consid-ered as a mere response function which would yield a torque on the particle depending on the time and the particle state.

The LBM EBF solver was first ran for a long time in order to eliminate the transient part of the solution. Two separate solver instances were then

(37)

executed in parallel, both starting from the post-transient state but with one of them with a slightly perturbed particle state, so that

kx − x0k₂ =p(θ − θ0₎2_{+ (ω − ω}0₎2 _{= δ}

0 (5.22)

where x = (θ, ω) and x0 = (θ0, ω0) are the two particle states and δ0 the

initial distance.

At regular time intervals, the particle state of the perturbed solution was normalized in accordance with the iteration scheme introduced in section 5.2:

δ ← p(θ − θ 0₎2_{+ (ω − ω}0₎2 δ0 , (5.23) θ0 ← θ + θ 0_{− θ} δ , (5.24) ω0← ω + ω 0_{− ω} δ (5.25)

The flow state was left unchanged. It was assumed that if the normalization was performed often enough, the particle state would change only slightly and the associated flow state would be consistent.

The δ-values were saved during the simulation and the Lyapunov expo-nent was evaluated with eqn. (5.13).

5.5 An adaption of the Rosenstein method

The method due to Rosenstein et al. [23] is commonly used to evaluate the LLE from experimental data. When one has performed an experiment, the result is most likely to be a long timeseries of some measurable quantity. This quantity will depend on the state of the system. In order to find the LLE, it is however not the the time evolution of the state itself that is of importance but how nearby trajectories in phase space are behaving.

The main idea of the Rosenstein method is that a chaotic trajectory will over time travel through a large part of the phase space. Time data can thus be regarded as space data. When the complete state of the system is not known, the attractor is usually reconstructed using a technique known as delay coordinates. In the case of this thesis, the complete state of the particle is known from the simulation data.

The first step is to collect statistical data about how nearby trajectories in the timeseries diverge over time. For this purpose, a nearest neighbor search is first performed for each point, where the distance is given by

kx_j− x_ik =q(θj− θi)2+ (ωj − ωi)2, (5.26)

where xj = (θ(tj), ω(tj)). It is further required that the time of the two

(38)

In order to avoid numerical errors, another restriction is imposed by

rmin< kxj− xik < rmax, (5.27)

where rmin = 10−8 and rmax = 0.1 has been used.

Due to the symmetry of the problem, the phase space is invariant to translations in θ by multiples of π. All angles have therefore been mapped to the interval [0, π) prior to the neighbor search.

The separation at a later time between the point xi and its neighbor xj

is found by

δi(tk) = kxi+k− xj+kk, (5.28)

The LLE can now be found by fitting a curve to the average logarithmic growth rate:

hln δi(tk)i = λtk+ const. (5.29)

where the mean hln δi(tk)ii is taken over all δi at time tk.

This fit was made for a number of timesteps corresponding to one wall oscillation period. The goodness of the fit, or the R2-value, will play a crucial role in the classification as will be shown in the next chapter.

(39)

Chapter 6

Method evaluation and

validation

6.1 Velocity boundary condition comparison

In chapter 3, two different velocity boundary conditions were introduced -the Zou & He and -the Regularized condition. A recent study [12] compared these two for a steady 2D Couette flow, and found a difference in numerical stability. Furthermore, a difference in accuracy was found for the case of a channel flow driven by a pulsating pressure gradient.

In order to asses the accuracy of the two conditions, the flow problem with oscillating walls was solved and compared with the analytical solution from chapter 2. The frequency was set to f = 0.1 and two different channel Reynolds numbers were used, respectively yielding β = 2 and β = 5. The relaxation time was set to τ = 0.55 for β = 2 and τ = 0.51 for β = 5, and the wall velocity was adapted to fit the Reynolds numbers.

All simulations were made on a quadratic domain with Nx= Ny = N for

various values of N . The simulations were started from a zero velocity profile and was each ran for five and nine wall oscillation periods, respectively for each β. The discrete L2 norm of the local errors was averaged over the last

period, and the total error was thus evaluated as

e = 1 Nt Nt X i=0 v u u t 1 NxNy Nx X j=0 Ny X k=0 uLBM(ti, xj, yk) ULBM w − u_analytic(ti, yk) 2 (6.1)

where Nt is the number of LBM iteration steps in the last period.

Figure 6.1 shows the result of the simulations. Visual inspection shows that both conditions yields a second order convergence in space. Further, the Zou and He condition is more accurate than the Regularized. The former appeared however to be very unstable when coupled with the EBF. The Regularized condition will therefore be used in the remaining part of

(40)

CHAPTER 6. METHOD EVALUATION AND VALIDATION 101 102 103 10−5 10−4 10−3 10−2 N average error

Zou and He condition Regularized condition (a) β = 2 101 102 103 10−4 10−3 10−2 10−1 N average error

Zou and He condition Regularized condition

(b) β = 5

Figure 6.1: Accuracy of the Zou & He and the Regularized boundary con-dition for the oscillating Couette flow case.

this thesis, even though it should be noted that it is less accurate than its counterpart.

6.2 Largest Lyapunov Exponent calculation

6.2.1 Calculation using two parallel simulations

The new LLE method was applied to a large range of parameter values. The method showed promising results for heavy particles. Figure 6.2 shows a bifurcation diagram and the LLEs found for a range of -values and for a case with a density ratio of α = 100. The bifurcation diagram shows all the values of the angular velocity at times that are an integer multiple of the wall oscillation period. If the system is periodic with the wall oscillation period, only one point will appear. A chaotic orbit will be aperiodic and the angular velocity will therefore be different at each period crossing, yielding a dense line in the diagram. Positive Lyapunov exponents are found for all the dense regions, which indicates that the algorithm yields correct results. The method did however struggle for lighter particles. An example is provided in figure 6.3 which was evaluated for α = 10. The phase plot shows that the system has reached a stable limit cycle which would not exist if the system was chaotic, but the LLE estimate is converging to a positive value. The main issue with the method is that the flow field is not updated when the particle state is normalized. When the particle is heavy, the timescale of the flow is much faster than that of the particle and the flow will have time to adjust to the updated particle state before influencing the particle rotation too much. For lighter particles this is no longer the case.

(41)

CHAPTER 6. METHOD EVALUATION AND VALIDATION 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 ω 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.02 0 0.02 0.04 ε λ

Figure 6.2: Bifurcation diagram (top) and the largest Lyapunov exponent found with the new method (bottom). Parameters used are κ = 0.5, Rep =

1, f = 0.1, St = 10. −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 cos(2θ) ω 200 300 400 500 600 −0.2 0 0.2 0.4 0.6 0.8 t λ

Figure 6.3: Phase plot (left) and convergence history of the Lyapunov ex-ponent (right). κ = 0.5, Rep = 1, St = 100

(42)

CHAPTER 6. METHOD EVALUATION AND VALIDATION 100 samples/period 500 samples/period f λref λ R2 λ R2 0.0500 -0.0130 0.0043 0.5637 0.0071 0.6009 0.1000 0.0121 0.0183 0.6499 0.0159 0.8582 0.1500 0.0356 0.0390 0.9614 0.0391 0.9968 0.2000 0.0335 0.0338 0.9507 0.0340 0.9958 0.4000 -0.1663 0.0083 0.1836 0.0271 0.4428

Table 6.1: Largest Lyapunov Exponents estimated with the Rosenstein method and the goodness of fit parameter R2. The reference values, λref

are those obtained with the algorithm from section XX for the FI system.

6.2.2 The Rosenstein method

The implementation of the Rosenstein method will now be evaluated. The LLEs for the CF system with fSt = 3.0, = 0.8 and for f = 0.05, 0.1, 0.15, 0.2 and 0.4 will be used as reference. The reference values are the ones found with the algorithm from section 5.3.

For the Rosenstein implementation, a timeseries with equidistant sam-ples was generated using the ode45 function in Matlab. The timeseries consisted of 100 periods and the number of samples per period was varied. The estimated LLEs with the two methods as well as the R2-values for the curve fit are shown in table 6.1.

It can be seen that the method systematically overestimates the value of the Lyapunov exponent. The positive LLE values found are however of correct order. For the negative values, positive Lyapunov exponents are found, but the goodness of the curve fits are quite low. In order for a solution to be accepted as chaotic, a R2-value of at least 0.95 will be required.

It should however be noted that the R2-requirement will make the algo-ritm not accept the f = 0.1 case as chaotic. The difference between this case and the other two chaotic cases is that the Lyapunov exponent is smaller. One can therefore conclude that the algorithm should be able to classify chaotic solutions that have a large enough Lyapunov exponent.

(43)

Chapter 7

Effect of aspect ratio, fluid

inertia and particle inertia in

a steady Couette flow

The case when the walls of the channel are moving with a constant velocity has been previously studied by, among others, Ding and Aidun [5]. It has been found that the rotation period of the particle diverges when the particle Reynolds number exceeds a critical value Rec. By regarding the particle

dynamics from a non-linear system point of view, Ding and Aidun identified this transition as a saddle node bifurcation. The nature of this bifurcation will be explained below.

In this chapter the effect of aspect ratio and particle inertia on this state transition will be investigated.

7.1 The transition from tumbling to steady state

for a neutrally buoyant particle

In order to show the nature of the state transition, phase plots for a cylinder with aspect ration kb = 0.25 ( = 0.9), density ratio α = 1 and for two

different values of Rep are provided in figure 7.1. Due to the symmetry

of the problem, all angles have been mapped to the [0, π)-interval and the orbits leaving the plot on one side are to be thought of as continuing on the other.

The lines in the phase plot show how different initial conditions evolve over time. For Rep = 1, all initial conditions converge to a limit cycle,

and the cylinder is rotating with a tumbling motion. For Rep = 10, the

trajectories are attracted to a single point - an attracting fixed point. This fixed point is on the ω-axis, and the cylinder will thus always eventually end up in a state of no rotation.

(44)

CHAPTER 7. EFFECT OF ASPECT RATIO, FLUID INERTIA AND PARTICLE INERTIA IN A STEADY COUETTE FLOW

0 0.5 1 1.5 2 2.5 3 −1 −0.8 −0.6 −0.4 −0.2 0 θ mod π ω (a) Rep= 1 0 0.5 1 1.5 2 2.5 3 −1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 θ mod π ω (b) Rep= 10

Figure 7.1: Phase plots showing how initial conditions evolve over time for Rep = 1 and 10. The parts of the lines with ω < 0 (> 0) are going left

(right). For Rep = 1 all intial conditions on the ω axis are converging

to a limit cycle. For Rep = 10 two fixed points are present of which one

is attracting (filled o) and one is semi-repelling (o). All initial conditions converge to the stable fixed point in this case.

Another fixed point on the ω-axis can be identified by considering the orbits near θ ≈ 0. These orbits are repelled from a point close to θ ≈ π. Through physical reasoning, Ding and Aidun identified this point as a saddle node (repelling in one direction and attracting in the other).

When the Reynolds number decreases, the fixed points approach each other. There will thus be a critical value, Rec, at which the fixed points

merge and then disappear.

Let θsand θusdenote the angle of the stable and the unstable fixed point,

respectively. The values of these angles can be evaluated using an interval halving algorithm developed by Ros´en et al. [24]. By evaluating these angles for a range of Rep-values, the critical particle Reynolds number where the

bifurcation takes place can be found.

7.2 Effect of particle inertia

When the solid-to-density ratio is increased, the system undergoes a second bifurcation. Figure 7.2 shows phase plots for different values of α and for Rep = 10. The position of the fixed points are not affected by α.

For α = 10, one can clearly see that some initial conditions converge to a limit cycle and some to the stable fixed point. As α decreases, the limit cycle approaches the unstable fixed point and somewhere between α = 7.4 and 7.6 the two will be connected. Below this critical value of the density

The chaotic rotation of elongated particles in anoscillating Couette ow

Master Thesis