Dynamics of settling disks in a cross-flow

Physical properties of rocks depend strongly on its grain packing which is determined by the natural process involving fluid-structure interaction during deposition. Here, we address the dynamics of settling solids with different sizes in a cross-flow numerically. This is an important aspect in rock formation process. Computational fluid-structure interaction usually involves the use of body-conformed grid. For the case of solid moving in a viscous fluid, the body-conformed grid has to be reconstructed in every time iteration which can be time consuming. The so-called Immersed Boundary method can be used to eliminate this grid reconstruction process. In this study, we use the couple of Immersed Boundary and Lattice Boltzmann Method (usually termed as IBLBM). In lattice Boltzmann, one does not require solving Poisson equation for pressure and a computation of intermediate field as in the usual fractional step procedure since velocity and pressure field can be obtained from the locally solved distribution function. These advantages make IBLBM can be easily employed. In this study, we perform simulations on the dynamics of a cluster of disks with two different values of diameter in a cross-flow. We compute the velocity distribution of the disks to study the sedimentation dynamics and we compare the results with a sedimentation process in a quiescent fluid. For the solid-fluid density ratio of 5 and incoming flow Reynolds number of 12, the velocity distribution is close to the normal distribution as also observed in the case of sedimentation in quiescent fluid.


Introduction
Motion of solids in viscous fluid are very important aspects in many natural phenomena and have many applications in engineering.An important mechanism where they become very challenging is in sedimentation process which involves in rock formation.This process governs the final structure or packing of formed rocks and in consequence their pore structure.Here we would like to study numerically how a number of solids moves in viscous fluid under the influence of a body force, gravity for instance.Since viscous effect should be taken into account, particle motion in inviscid fluid is not used here.Difference between particle motion in inviscid fluid and in viscous Newtonian fluid is discussed thoroughly by [1].Also, the relative velocity between the solid and the fluid may take values such that inertial effect cannot be neglected as in the Stokes approximation.
Transportation of many solids in fluids can be viewed as a convection-dispersion process which can be described by considering the evolution of particle concentration and mass of deposited particles in time and space, as was done by [2] and [3] for particle transportation in porous medium with uniform flow under saturated conditions.With the increasing of computer speed and memory, numerous studies have been performed on the deposition of finite size particles to capture both the particles dynamics and the hydrodynamics around the particles, for instance by [4], [5], [6].Recently, [7] performed a numerical study on the dynamics of square particle using Immersed Boundary Adaptive Mesh Refinement technique to identify different settling dynamics due to initial inclination angle and solid-fluid density ratio.One may solve the governing equations for fluid dynamics, namely the Navier-Stokes equations, numerically by imposing moving boundary conditions to capture both the dynamics of the solids and the fluid flow around them.As an alternative, the flow can be viewed in mesoscopic scale by solving the discretised Boltzmann transport equation which is known as the Lattice Boltzmann method.It is an approach to the incompressible Navier-Stokes equations under the assumption of low Mach number [8].In this study, we use the Lattice Boltzmann method to obtain the flow fields.
In Lattice Boltzmann, fluid is regarded as distributions of particles which are evolved via two main processes which are streaming and collision.Velocity and pressure field are macroscopic variables which can be computed from the particle distribution function once it is obtained.Since transportation of fluid is regarded as the evolution of particle distribution function via streaming and collision, the distribution function can be computed locally which eliminates the need of solving linear system of equations and thus gives less time computing.Description of collision in Boltzmann transport equation can be very complicated.Thus, a collision model is often used such as the Bathnagar-Gross-Krook (BGK) model [9] which greatly simplifies the collision description.Also, the computation can be further simplified since one does not need to solve a Poisson equation to compute pressure field.These features make Lattice Boltzmann method an efficient alternative to conventional Navier-Stokes solver.More on theoretical background and implementation of Lattice Boltzmann method can be found in literature, for instance in [9].
For more general geometry, body-conformed grid can be used to improve smoothness of the solid-fluid interface at which the no-slip condition has to be imposed.A regeneration of the body-conformed grid must then be performed in each time step if the solid moves.It is then natural to use a non-body-conformed method to eliminate the time-consuming grid regeneration and to simplify numerical implementation.As an alternative, one may use the idea of immersed boundary method proposed by [10] in order to eliminate the grid generation.This method involves the use of a set of Lagrangian points which represent discretised solid surface.One then assigned forces at each Lagrangian point.The fluid flow is solved in an Eulerian grid.The assigned forces are then spread into the fluid by adding forcing term in the governing equation of the fluid flow so the fluid and the solid can interact.This idea has been extensively used and developed to study solid-fluid interaction (see [11] for vast implementations of immersed boundary method).Recent study on the variants of immersed-boundary was presented by [12] in the study of fixed packed bed reactor.
In Lattice Boltzmann method however, one imposes the no-slip condition by performing a bounce-back procedure (see [9] for details).Therefore, it is then natural to use the idea of bounce-back procedure to include solid-fluid interaction.Nevertheless since the solid-fluid interface is moving, the bounce-back procedure has to be adapted, for instance as was proposed by [13].The use of bounce-back procedure, however, produce a sudden change from non-fluid node to fluid nodes (vice-versa) as the solid moves.The hydrodynamic forces computed at the interface is then fluctuating during the movement of the solid.The use of Immersed Boundary method to account the presence of the solid thus eliminates the sudden change of solid-fluid nodes since the fluid flow is solved in an Eulerian coordinate, with an additional forcing term to mimic the presence of the solid.
Earlier attempt to couple Immersed Boundary to the Lattice Boltzmann method was performed by [14].In their work, the force assigned on each immersed boundary point is a kind of restoration force which is a linear spring model.The assigned forces are then spread into the fluid.The rigidity of the solid is then controlled by the choice of proportionality coefficient in the spring model.A fiber model can also be used to assign forces at each Lagrangian point as was done by [15].The Lagrangian forces are computed by modelling tension, bending, and fastening.The spreading of Lagrangian forces are done by the use of a discretised Dirac delta function (see [16] on how to construct numerical approach of the Dirac delta function).Once the spreading procedure is performed, one then may use the spread force in the forcing term for the Lattice Boltzmann equation.There are several formulation of the Lattice Boltzmann forcing term, for instance one may use the formulation proposed by [17] or [18].
Here, we use a coupling procedure described in [19].In this coupling procedure, the forces assigned at each Lagrangian point are computed by means of interpolation of the macroscopic fluid velocity.The no-slip condition at the interface is imposed by taking the fluid velocity at the location of the Lagrangian points equals to the velocity of the solid boundary.The spreading procedure is however modified.To better ensure the no-slip condition, [19] used a modified spreading operator proposed by [20].Direct computation of forces from the velocity thus eliminates the need of calculating predictive velocity as commonly used in Immersed Boundary-Lattice Boltzmann coupling.Also, since fluid velocity and pressure are computed directly from the fluid particle distribution function, one requires no pressure-correction step as usually done when Immersed boundary is coupled with conventional Navier-Stokes solver based on finite differences or finite volumes which simplifies the computation even further as stated in [21].
In this article, we first describe the general idea of coupled Immerse Boundary-Lattice Boltzmann method proposed by [19].We then show the application of the coupling to simulate moving solids in fluid from single to many disks in the context of sedimentation process.The last part would be the discussions on sedimentation dynamics of multiple disks in the presence of cross flow.This study is performed in two-dimensional configuration under the assumption that the solids take circular shape (disks).

Immersed Boundary-Lattice Boltzmann coupling
In this section, we introduce briefly the idea of Immersed Boundary and Lattice Boltzmann method as well as their coupling we use in this study.The coupling gives us the advantages from both methods.

Lattice Boltzmann
The Lattice Boltzmann method involves streaming and collision of fluid particle density function n which depends on function of position (x), momentum (p), and time (t).Suppose that the density function is transported from one state at time t to another state at t + dt.The evolution of the distribution function with Bathnagar-Gross-Krook (BGK) collision model is given by In relation (1), the right hand side represents collision that might happen during transportation of the distribution function on the interval dt called the Bathnagar-Gross-Krook (BGK) collision model.This model regards the collision process as a relaxation of the distribution function towards a local equilibrium n eq (x, p, t) where τ gives the relaxation time.This equation is solved in its discrete form, namely the Lattice BGK model which is known as the Lattice Boltzmann method.As stated in previous section, the Lattice Boltzmann method can be regarded as an approximation to the Navier-Stokes equations for low Mach number [8].The space discretisation for the Boltzmann transport equation can be taken as a regular grid; we take equal spacing of ∆x = ∆y = 1.Momentum can be discretised as some possible moving directions a distribution function can take.For two-dimensional configuration, a popular discretisation model is the D2Q9 model in which nine possible directions are considered for momentum.The D2Q9 model is used throughout this study.The Lattice Boltzmann equation for D2Q9 model is given by Note that e k gives the k-th momentum directions.In general those directions can be written as e = (e x , e y ) where e x = {0, 1, 0, −1, 0, 1, −1, −1, 1} and e y = {0, 0, 1, 0, −1, 1, 1, −1, −1}.Here ∆t is taken as unity which gives lattice speed c = ∆x/∆t = 1.For the D2Q9 model, the equilibrium distribution function n eq k can be computed by the following relation where c s is the lattice speed of sound and is given by c s = c/ √ 3. The weighting factor ω k is given by Viscosity of the fluid is related to the relaxation time parameter via the following relation It is clear from relation ( 5) that τ has to be larger than 0.5 to get a physically meaningful kinematic viscosity.The last term in the r.h.s. of relation ( 2) is the forcing term.In this study, we use the forcing formulation proposed by [18] which is given by where f is the force exerted on the fluid.The macroscopic velocity can be computed from the distribution function by using the following relation with ρ(x, t) = 8 k=0 n k (x, t).Since the goal in Lattice Boltzmann method is to compute the distribution function, one does not need to solve a Poisson equation to compute pressure field as in the usual Navier-Stokes solver.Together with the use of regular grid, this makes Lattice Boltzmann can easily be implemented.

Immersed Boundary
Instead of using bounce-back procedure on moving boundary, we make use the idea of Immersed Boundary method proposed by [10] (see also [16]).We follow the calculation of forces F IB (x s , t) at each Lagrangian point where x s = x (X s , t)(see figure 1) by using relation 7 as proposed by [19].Supposed that we have obtained the value of n k and ρ at t + ∆t, F IB (x s , t) can be calculated as follows where U(x s , t+∆t) is the velocity of the s-th immersed boundary points.Next, F IB (x s , t + ∆t) will then be spread into the surrounding fluid.An influence domain of the force is required.It is defined by a Dirac delta function δ (x − x s ). Figure 1 shows the support domain of a Lagrangian point X s .We use an approximation of the Dirac delta function proposed by [22].As stated before, the computed Lagrangian forces are spread into the fluid by using a spreading operator proposed by [20].The idea is to define characteristic strip-width ϵ(x s , t + ∆t).The forces acting on the fluid f (x, t + ∆t) is given by the following relation A requirement ∆s = ∆x has to be satisfied to obtain a good value for ϵ (see [20] for further details).Note that for rigid bodies ϵ(x s , t + ∆t) can be computed once at the beginning of computation.The dependency of ϵ with respect to the movement of the immersed boundary for rigid and flexible bodies can be found in [21].In this study, ϵ is computed only at initial time since there is no significant change on the computed ϵ for rigid body which saves significant computing time.Now that f (x, t+∆t), u(x, t+∆t) are obtained via relation (7), we can proceed to the next time iteration.

Moving rigid body under gravity.
The equation of motion derived from Newton's second law has to be solve for translational and rotational movements which are given by (see [23] for details): Here, a simple Euler explicit algorithm is used to advance equations ( 10) and (11) in time.
In the case of sedimentation, the solids may have varying velocity.Therefore, it is natural to introduce a dimensionless number other than the Reynolds number.Here we use the Archimedes number which is defined as Throughout the study, ρ rb ρ f = 5 and Ar = 92.16are used.

Sedimentation of N -disks.
Here we show the application of the coupled Immersed Boundary-Lattice Boltzmann method for sedimentation of N -disks in viscous fluid.In this study, we use 149 disks with two values of diameter.There are 105 disks having diameter D 1 = D and 44 disks having diameter D 2 = 2D where D = 40∆x.The width L x of the domain is 75D and its height L y is 45D.The problem configuration is given in figure 2. During the sedimentation process, the disks may have interaction amongst them or with walls.Therefore we use the repulsive force formulation by [24] to include disk-disk collisions.We also put limits such that the disks cannot go towards the open boundaries of the domain by applying the same interaction at x/D = 18.75 and x/D = 56.25.
Two cases are considered in this study.First, the disks are settling in an initially quiescent fluid.The second one is sedimentation in a cross-flow originating from the left open boundary in figure 2 with velocity u in = (u 0 , 0) where u 0 is taken as 0.05 lu/ltu ('lu' is lattice units and 'ltu' is lattice time units).Note that this value of u 0 , along with D, will be used as normalisation factor for both cases.
Let us consider now the first case where the fluid is initially at rest, i.e., no crossflow imposed.In this case, extrapolation conditions are imposed at all open boundaries.The extrapolation is used to compute the unknown distribution function on the boundaries by means of second order polynomial [25].For instance, on the outlet (right side boundary) we have The no-slip condition for fluid in contact with walls is imposed by using the bounce-back procedure.At early stages of the sedimentation process, the fluid is forced to move around the pack.Disks on the left and right ends of the pack are then forced to move upwards by the fluid.The interaction also creates two counter-rotating flows behind the pack as its center part moves downwards.The pack is then settling down as the disks arrange themselves at the bottom of the domain.Some of the disks, which moved upwards at the early stages, are then moving downwards.In order to find a pattern on the movement of the disks, we construct the velocity distribution of the solids.The distribution of velocity component perpendicular and parallel to gravity (normalised by the maximum value of each component) are given in figure 3 and figure 4 for different time.Note that here we only take four snapshots between the early stages of the sedimentation and right before the disks are settling at the bottom of the domain.The velocity distributions are all compared to the normal distribution.We also extract the information of mean value at each snapshot.
The velocity component perpendicular to gravity shows a good fit to the normal distribution for the four snapshots as shown in figure 3. The mean value changes during the sedimentation which indicates different overall horizontal motion even if there could be a small number of disks which moves in a direction opposite to the overall motion.At tu 0 /D = 31.25, it has mean value of 0.11 which means that most of the disks are moving to the right.However, there are a small number of disks (initially at the right part of the packing) moving to the left due to the fluid portions moving upwards at the right end of the packing and recirculating back in a form of a counterclockwise vortex.The mean value is very small at tu 0 /D = 62.5 although it is negative.This does not mean that there is no horizontal dynamics, it is merely because the number of disks moving to the left is almost balanced by the number disks moving to the right.It is depicted by an almost symmetric distribution.For tu 0 /D = 93.75 and tu 0 /D = 125 most of the disks have small horizontal velocities.
The mean value for the velocity parallel to gravity are negative for tu 0 /D = 31.25,tu 0 /D = 62.5 and tu 0 /D = 93.75 which show that the population move downwards.However, at tu 0 /D = 125 the mean value is positive since most of the disks bounce at the bottom boundary.It can also be said from figure 4, that during the downward movement that small amount of disks move upward due to the interaction with the portion of fluid which moves upward.An example of flow visualisation is shown in figure 5.
In the second case, a constant velocity boundary condition u in = (u 0 , 0) is imposed at the left

Conclusions
Simulations of settling circular solids with 2 groups of diameters in quiescent fluid and in a crossflow Reynolds number of 12 (defined as Du 0 /ν) have been done based on Immersed Boundary-Lattice Boltzmann method.In the two conditions, the solid velocity distributions are closed to normal distribution.Interestingly, we do not observe very different tendencies that the velocity distribution would obey other kind of distribution when the sedimentation process is happening in a cross-flow.This may be due to the cross-flow is still laminar and the solid-fluid density ratio of 5 is sufficiently large for the flow to easily affect the solid movement.
For future work, imposing turbulent flow may cause different dynamics due to the presence of multi-scale eddies.Also, the population has only two groups of diameters.A wide variety of diameter may cause richer dynamics during sedimentation process.Another issue for further development is the spatial resolution as any two disks approach each other or when any disk approach a solid wall.The dynamics of the fluid may be under-resolved at the space between the contact.One may naturally think to do a mesh refinement in that particular region.This requires remeshing at many near-collision regions which could be anywhere inside the

10thFigure 1 :
Figure 1: Discretized immersed boundary points (red dots) rectangular grid.Dashed rectangle represents the support domain of the s-th Lagrangian point X s .Grey dots are the points where flow fields are computed.

Figure 2 :
Figure 2: Configuration of the problem: Boundary conditions used at open boundaries and initial positions of the disks as well as the initial arrangement of the disks.

Figure 5 :
Figure 5: Vorticity visualisation during the sedimentation process at tu 0 /D = 62.5 for sedimentation in initially quiescent fluid.The colour scale on the bottom right corner is the normalised vorticity ωD/u 0 .

Figure 8 :
Figure 8: Vorticity visualisation during the sedimentation process at tu 0 /D = 18.75 in crossflow.The colour scale on the bottom right corner is the normalised vorticity ωD/u 0 .