Sensitivity analysis of Immersed Boundary Method simulations of ﬂuid ﬂow in dense polydisperse random grain packings

. Polydisperse granular materials are ubiquitous in nature and industry. Despite this, knowledge of the momentum coupling between the ﬂuid and solid phases in dense saturated grain packings comes almost exclusively from empirical correlations [2–4, 8] with monosized media. The Immersed Boundary Method (IBM) is a Computational Fluid Dynamics (CFD) modelling technique capable of resolving pore scale ﬂuid ﬂow and ﬂuid-particle interaction forces in polydisperse media at the grain scale. Validation of the IBM in the low Reynolds number, high concentration limit was performed by comparing simulations of ﬂow through ordered arrays of spheres with the boundary integral results of Zick and Homsy [10]. Random grain packings were studied with linearly graded particle size distributions with a range of coe ﬃ cient of uniformity values ( C u = 1.01, 1.50, and 2.00) at a range of concentrations ( φ ∈ [0 . 396 , 0 . 681]) in order to investigate the inﬂuence of polydispersity on drag and permeability. The sensitivity of the IBM results to the choice of radius retraction parameter [1] was investigated and a comparison was made between the predicted forces and the widely used Ergun correlation [3].


Introduction
Multi-phase flows involving granular materials occur in many industrial processes for example in the oil and gas, construction and pharmaceutical industries as well as in natural processes such as landslides, scour and erosion, and sedimentation.Knowledge of the momentum coupling between the fluid and solid phases is therefore of great commercial and societal value yet a general treatment over a wide range of solid volume fractions, φ, and Reynolds number, Re, is at present not possible due to the complex nature of these systems.
For flows involving saturated soil and sand packings such as those occurring in internal erosion and liquefaction processes, two independent terms, the pressure gradient or buoyancy force, F b , and the drag force, F d , dominate the fluid-particle interaction.Stokes law gives the drag force acting on a sphere of diameter D p moving through a fluid with dynamic viscosity μ in the limit φ → 0, Re → 0, In the context of grain packings, U is the superficial or apparent flow velocity which can be expressed in terms of u, the interstitial flow velocity, as U = (1 − φ)u.As the concentration of solids and the Reynolds number are increased the drag force varies with a complex dependence on both φ and Re.
e-mail: christopher.knight13@imperial.ac.uk Previous experimental and numerical efforts to characterise the dependence of F d on φ and Re have focussed on monodisperse systems.Broadly speaking the experimental approaches are either based on fluidisation experiments [3] or sedimentation experiments [2].Hill et al. used the Lattice Boltzmann (LB) method to study the drag force in ordered and random arrays at low to moderate Re [4].More recently Immersed Boundary Methods (IBM) [8] have been used to provide more accurate results, improving on these monodisperse drag correlations.
In reality granular media typically exhibit not just a single grain size but a range of grain sizes.This size polydispersity affects the packing density and pore space configuration of dense granular materials significantly with implications for the fluid-particle interaction forces.Due to the prevalence of polydisperse media in nature and industry there is a clear need for improved understanding of the relationship between grain size distributions and fluidparticle momentum coupling.
An outline of the IBM is given in section 2, followed by a description of validation work conducted for the IBM in section 3.In section 4 results of IBM simulations with random packings are presented and finally, conclusions are given in section 5. tries, such as the pore space in soils.Fluid is represented by a fixed structured grid of Eulerian cells in which the Navier-Stokes equations for mass and momentum conservation are solved locally.Solid surfaces are represented by rigid meshes of points inserted into the fluid domain at which no-slip, no-penetration (ns-np) boundary conditions are imposed through additional forces applied locally to the fluid equation system.These points can move relative to the fixed Eulerian cells adopting arbitrary positions and are therefore referred to as Lagrangian points.An illustration of the two meshes is shown in figure 1.In the following, lower and upper case symbols are used to refer to Eulerian and Lagrangian variables respectively.In the present work a multi-direct forcing type IBM was used to solve the Navier-Stokes equations for the flow of fluid in the pore spaces between particles and determine the fluidparticle interaction forces.
Uhlmann [9] developed a direct forcing type IBM applicable to structured, uniform, cubical fluid meshes in which the additional Lagrangian forces are smoothed out across the Eulerian domain by interpolation and spreading processes using a discrete Dirac delta function (DDF).Forces to be applied at the Lagrangian points on the surface of an immersed object are determined explicitly from the current velocity of the fluid in the local neighbourhood of the point by requiring that a no-slip condition is enforced.The required velocity, U d (X l i ), at position X l i of Lagrangian point l on the surface of object i is determined by the rigid body motion of the object, where x c i and v i are the position and translational velocity vectors of the centre of mass of object i respectively and ω i is the rotational velocity vector about x c i .Expressing the fluid momentum conservation equation in terms of the fluid-particle interaction force and inserting the desired velocity from Equation (2) into the transient term gives the force required to impose the no-slip condition at each Lagrangian point, Here U N (X l i ) is the current value of the fluid velocity interpolated at the Lagrangian point at time level N, Δt is the fluid time step and RHS N (X l i ) represents the interpolated values of the derivative fields corresponding to the pressure, viscous and convective terms in the momentum conservation equations.The Lagrangian forces are then spread across the fluid cells in the support region of the DDF of each point and summed up to give the additional forcing in each Eulerian fluid cell.The force f applied to Eulerian cell m with centre at x m is then given by where δ h is the three point DDF of Roma [6] in the present work, ΔV l i is the Lagrangian point volume and the sums are taken over all N P particles and all N L Lagrangian points of each particle.As a consequence of the overlapping support regions, there is some interference of one Lagrangian point with another.This means that in practice the ns-np condition is not well enforced after a single round of forcing.Luo et al. [5] proposed the multi-direct forcing scheme in which forces are applied iteratively until a convergence criteria or maximum number of iterations has been reached.
Another side effect of the interpolation of fluid velocities and spreading of forces is the enhancement of particle size which occurs as fluid in the immediate vicinity of an immersed object is frozen in the diffuse cloud of forces surrounding the object.In the case of the Roma DDF [6], with a support region three cells wide, this means that the smeared out particle interface extends a distance 3Δx/2 into the fluid where Δx is the Eulerian grid spacing.As F d is proportional to the projected area of the particle this enhancement of size has the undesirable consequence of increasing F d .Breugem [1] suggested the radius retraction scheme in which the Lagrangian points are placed a small distance r d = bΔx inside the particle interface in order that the actual particle size is better approximated after the interpolation and spreading procedures by suitable choice of the parameter b.Breugem found empirically that the choice b = 0.3 was optimal, giving second order convergence of the predicted drag forces with increasing grid resolution.Tang et al. [7] showed that b = 0.3 is only optimal in the low concentration limit previously considered by Breugem, with the optimal value becoming more sensitive to φ and Re at higher concentrations.

Validation of the IBM with ordered packings
Flow through ordered arrays of monodisperse spheres in the low Re limit has been used for validation purposes by several researchers for LB [4] and IBM simulations [1,8].Zick and Homsy [10] considered flow through periodic simple cubic (SC), body centred cubic (BCC) and face centred cubic arrays (FCC) in an infinite space in the limit Re → 0 from dilute to maximum packing fraction in each case (SC: φ max S C = 0.5236; BCC: φ max BCC = 0.6802; FCC: φ max FCC = 0.7405).The Navier-Stokes equations were expressed in terms of a set of two-dimensional integral equations of the hydrodynamic stress at the sphere surfaces.For each configuration and concentration considered, Zick and Homsy determined a drag coefficient K = F d /F s d numerically by evaluating these integral equations using increasing numbers of terms and extrapolating the results to infinity.The results obtained by Zick and Homsy were shown to compare well with other values available in the literature.
The steep dependence of K on φ in the densely packed limit means that any small variation in the effective particle diameter as φ → φ max leads to large degrees of error in the predicted values of K. To establish confidence in the predicted drag forces the influence of r d on K was investigated using the results of Zick and Homsy as benchmark values [10].SC, BCC and FCC cases were considered.In each case an inlet boundary condition was applied along the x direction with periodic boundary conditions in the ŷ and ẑ directions.Several layers of particles (4 for SC and 3 for BCC and FCC cases) were arranged along the x direction with only the forces acting on internal layers considered in order to minimise boundary effects.All simulations were conducted at Re ≤ 0.1 with the fluid time step chosen to yield a Courant number of C ≈ 0.1.For each configuration several values of φ were considered by varying the particle diameter whilst the centre positions and domain dimensions remained fixed.Δx was varied to satisfy the target resolutions D p /Δx = 16, 32, 48, and 64 for each configuration and concentration considered.
Representative results of the variation of K with φ for the BCC case for b = 0.0 and 0. mal, as suggested by Breugem [1], giving agreement to within 2% of the benchmark values at the low resolution of D p /Δx = 16.The results at φ = φ max BCC show that b = 0.3 leads to an underprediction of K with the optimal value appearing to lie somewhere in between b = 0.0 and b = 0.3.Figure 3 shows representative results for the variation of K with φ for the FCC case for b = 0.0 and 0.3 at D p /Δx = 16 and 48.Results at φ ≈ φ max BCC confirm that the optimal value of b lies in the range 0.0 < b < 0.3 whilst results at higher concentration suggest that b → 0.0 as φ → φ max FCC .

Drag in polydisperse packings
Random grain packings with linearly graded size distributions with coefficient of uniformity values C u = 1.01, 1.50 and 2.00 were prepared using the Discrete Element Method (DEM).Keeping the minimum diameter fixed at D min p = 0.50 mm, ensembles of particle diameters were generated using a pseudo-random number generator with standard deviations of the errors between the generated and target size distributions for each C u of less than 2%.A disperse cloud of particles was produced by choosing random initial positions with no particles overlapping.Isotropic compression to 1.0 × 10 2 kPa was then performed with periodic boundaries in all directions.The resulting particle positions and domain dimensions were then used to generate the Lagrangian points and fluid domain dimensions respectively for IBM simulations with an inlet/outlet condition along one coordinate direction and periodic boundary conditions on the Lagrangian points and fluid equations in the two orthogonal directions.A grid resolution of D min p /Δx = 24 was used and the inlet velocity, U, was chosen to satisfy Re ≈ 0.1 in all cases, with the fluid time step, Δt, chosen to yield a Courant number of C ≈ 0.1.In order to not over constrain the velocities and pressures in the pore spaces at the inlet and outlet, additional fluid cells were added in these regions such that all particles lay entirely inside the fluid domain along the inlet/outlet direction.Results for the total drag on each packing normalised by domain volume V against concentration for all C u values with b = 0.00 and 0.30 are shown in figure 4. For comparison the predictions of the widely used Ergun correlation [3] are shown where the characteristic diameter is taken as the average diameter by mass for each packing.It can be seen that in the range φ ∈ [0.39, 0.52] the Ergun correlation and the IB simulation results are in good agreement.At the highest concentrations considered there is a large discrepancy between the Ergun correlation and the IB results.As C u increases the surface area to volume ratio decreases and hence the total drag is seen to decrease for approximately equal concentrations.The results become less sensitive to changes in b as C u increases which can be explained by the decrease in dF d /dD p as D p increases.
Drag forces on individual particles are shown against diameter in figure 5.The Ergun correlation provides good predictions for the smallest particles, however the predic-tion diverges from the IB results as D p increases.In a polydisperse packing the larger particles experience a locally higher concentration than they would in a monodisperse packing at the same φ as small particles can fill the gaps between larger particles.A first order indication of this increase in local φ is given by the number of interparticle contacts, or coordination number, N z , which is seen to increase with increasing particle size.Monodisperse correlations do not account for this local increase in φ which explains the discrepancy.Particles in the inlet and outlet boundary regions where the flow is developing experience systematically less drag and are not shown in figure 5.

Conclusions
Results with ordered packings shown in section 3 support the findings of Breugem [1] that b = 0.30 is optimal for disperse systems whilst results at higher concentrations indicate that b → 0 as φ → 1 confirming the findings of Tang et al. [7] that the optimal value of b is dependent on φ. Results for random packings in section 4 show that the Ergun correlation [3] underpredicts the drag acting on individual particles in polydisperse packings with an error which increases with particle diameter which is expected as the effect of the local increase in φ with D p is not taken into account.To remove boundary effects for cases presented in section 4, alternate boundary conditions along the inlet/outlet direction will be explored.

Figure 1 .
Figure1.Illustration of Lagrangian points (red circles) immersed within an Eulerian fluid grid (blue lines).The Lagrangian points are placed a small distance r d from the true surface of the particle (solid black line).It can be seen that some overlap of the DDF support regions occurs for points in close proximity.

3 D 3 Figure 2 .
Figure 2. Drag coefficient K against φ for the BCC case.

3 Figure 3 .
Figure 3. Drag coefficient K against φ for the FCC case.

Figure 4 . 6 FFigure 5 .
Figure 4. Drag force per unit volume F d /V against φ for all values of C u considered with b = 0.00 and 0.30.Solid lines show the predictions of the Ergun correlation.