Modeling particle-fluid interaction in a coupled CFD-DEM framework

In this work, we present an alternative methodology to solve the particle-fluid interaction in the resolved CFDEM coupling framework. This numerical approach consists of coupling a Discrete Element Method (DEM) with a Computational Fluid Dynamics (CFD) scheme, solving the motion of immersed particles in a fluid phase. As a novelty, our approach explicitly accounts for the body force acting on the fluid phase when computing the local momentum balance equations. Accordingly, we implement a fluid-particle interaction computing the buoyant and drag forces as a function of local shear strain and pressure gradient. As a benchmark, we study the Stokesian limit of a single particle. The validation is performed comparing our outcomes with the ones provided by a previous resolved methodology and the analytical prediction. In general, we find that the new implementation reproduces with very good accuracy the Stokesian dynamics. Complementarily, we study the settling terminal velocity of a sphere under confined conditions.


Introduction
The recent advances in particle-fluid simulations allowed the description of many industrial and natural processes [1][2][3][4]. For instance, the sedimentation of particles in dilute and dense conditions and the flow of particles through submerged hoppers have been recently examined [1,2,5]. In general, experimental and simulated data agreed very well when exploring intermediate Reynolds number regimes. It suggests that the CFD-DEM coupling can resolve the hydrodynamics interactions in those scales with enough accuracy.
In the field of particle suspension modeling, researchers have been able to reproduce classical results, applying a variety of methods combining Eulerian and Lagrangian approaches [1,2,6]. When examining simulation domains that are much larger than the dimensions of the particle, unresolved methods are typically applied. In these methods, the local interactions are determined by auxiliary fields, which represent the particles in the CFD domain. In contrast, in resolved methods, the individual particles are fully discretized, and their dynamics are addressed with better accuracy. They are often used for smaller system sizes.
The resolved coupled CFD-DEM approach described in this paper is an extension of a previous Fictitious Domain Method (cfdemSolverIB) [2]. In this implementation, there is no body force acting on the fluid, which implies that the local pressure profile is solely determined by the local variations of the velocity field. Accordingly, to account for the buoyancy acting on the particle due to the presence of the fluid, an exact term (ρ f gV i ) is added, explicitly. * e-mail: jfonceca@unav.es In our approach, the gravitational field is directly included in the incompressible Navier-Stokes equation as a body force term. Consequently, the pressure profile also contains the gravitational contribution. That is why we no longer need to impose an explicit buoyant force acting on the particles. Thus, our approach does rely on the assumption that the use of absolute pressure to calculate the particle dynamical response is relevant when numerical corrections must be implemented at this scale. Although not relevant differences will be expected when a single particle is analyzed, such differences might be relevant for very dense suspensions where the local pressure varies due to the dynamical coupling between particles.
The present work is structured as follows: in Sections 2.1 and 2.2, we briefly explain the used CFD-DEM scheme. Moreover, in Section 2.5, we briefly comment on a theoretical framework that describes the settling of a sphere in confined conditions. Finally, in Section 3, we present some numerical results highlighting the validity of our methodology.

Fluid phase
The governing equation for the fluid phase is the incompressible Navier-Stokes equation: Our approach explicitly considers the body force term, which was not considered in the previous implementations [1,2,7]. When discretized for a CFD time step ∆t, equation 1 becomes: wherep is the estimated pressure,Û the interim solution for U and U t−∆t is the solution found in the previous time step. Likewise, the continuity equation holds: For each particle, we determine its occupation region T i in the CFD grid by applying a correction to the velocity field in all the corresponding cells. In the cell c ∈ T i , we retrieve the velocityÛ c and pressure p c fields.
Once we apply the velocity correction in region T i , the divergence-free condition for the velocity is violated in the rest of the simulation domain and equation 3 is valid only if we define U =Û − ∇Φ, where Φ is a scalar field. From equation 3, ∇ 2 Φ = ∇ ·Û, which is equivalent to apply a force-term to equation 1. This approach and its potential to resolve the motion of particles in liquids has been extensively discussed by Kloss et al. [8].

DEM Method
The Discrete Element Method was first presented in 1979 by Cundall [9]. One of the main characteristics of DEM simulations is its efficacy when resolving the granular medium at the particle scale. The DEM is a Lagrangian method, meaning that all particles in the computational domain have their trajectories solved explicitly in every time step. Thus, the DEM is capable of simulating a wide range of phenomena, such as dense and dilute particulate systems, as well as rapid and slow granular flow. Here, we briefly discuss the main features of the model, and a more detailed description of the method was published by Goniva et al. [1].
For a particle i sedimenting under gravity in a liquid, its trajectory is calculated using the following force and torque balances: and where i represents the index of the particle, m i is its mass, x i its position vector and I i its moment of inertia. N p represents the number of neighboring particles, and N w the number of near walls.
In Table 1 we describe the forces and torques from equations 4 and 5. The particle-particle F i, j and the particle-wall F i,w contact forces are written in terms of their ortoghonal and parallel components with respect to the plane of contact. Table 1. Summary of forces and torques acting on the particle.
For a full description of the contact forces, refer to [1].
In the method originally presented in [2], the buoyant force was explicitly accounted for in the total fluid force acting on the particles. The current implementation explicitly considers the gravitational field acting on the fluid phase. Thus, we represent the fluid-particle interaction force by the last equation in Table 1. It generalizes the original approach presented by Shirgaonkar et al. [10]. Importantly, the new approach includes a multiplicative term (1 − φ c ) in the pressure gradient accounting for the shape of the particle. Note, the void fraction field φ c represents the fluid occupation in the volume element c.

CFD-DEM coupling
The coupling procedure can then be described by the following steps: -Determine the position and velocity of each particle in the domain from the DEM state.
-Assign the particles to the cells that contain their centers.
-Solve the incompressible Navier-Stokes equations for the fluid phase using the PISO algorithm.
-Correct the velocity found in the previous step in the grid cells where particles are located.
-Calculate F i,D for all particles, this force is applied in the next set of DEM time steps.
-Lastly, the boundary conditions are applied in the fluid phase, and we return to the first step.

Simulation Details
Here, we simulate the motion of a single particle in a 3D viscous fluid. The radius and density of the particle are R = d/2 = 1 mm and ρ p = 7850 kg/m 3 , respectively. The fluid has density ρ f = 970 kg/m 3 and several kinematic viscosities are explored ν = [1000, 2000, 5000, 12500] in cSt. Specifically, we examine a sphere left to settle from rest in a quiescent fluid. In this scenario, we expect the particle to settle with terminal velocity v S = 2 9µ (ρ p − ρ f )gR 2 and characteristic time τ = 2 9µ ρ p R 2 . We set a grid resolution equal to ∆x = d/8 and CFD and DEM time steps ∆t = τ/35. In all cases, the total simulation time was 20τ seconds and the Reynolds numbers were inferior to 0.1.

Sphere settling between plane parallel walls
Complementarily, we study the settling terminal velocity of a sphere under confined conditions. Back in 1923, Faxén studied the effect of plane parallel walls on reducing the settling velocity of spheres [11]. In a low Reynolds number regime, for a sphere centered between two plane parallel walls separated by a distance H, Faxén proposed that the drag force applied to the particle by: R the sphere radius and v the settling velocity. This expression approximates the actual force acting on the sphere when it is located far from the walls (s −1 > 2). A more general solution for the drag force was later introduced by Ganatos et al., covering configurations that were out of the validity of Faxén's solution [12].

Validation
As a first step, we explore a system with dimensions 10d × 10d × 25d. Initially, the particle is centered in the xy-plane at the location of z = 17.5d. Then, it is left to settle from rest, it accelerates and reaches its terminal velocity. Figure 1 illustrates the system state obtained after 20τ seconds, using both approaches. It compares the kinematic pressure (Figs. 1a and 1b) and velocity (Figs. 1c and 1d) fields. Note that the cfdemSolverIB approach (Fig. 1a) does not account for the pressure gradient imposed by the external field in the fluid. It is also noticeable that the kinematic pressure perturbation of the particle movement is significantly smaller than the hydrostatic pressure gradient (Fig. 1b). Remark that when using cfdemSolverIB, the local pressure field is solely defined by the corresponding variations in the velocity field, as a result of the momentum balance equations. In our approach, apart from that local pressure contribution, we also have the gravitational field acting on the fluid. However, no relevant impact is detected in the velocity fields (Figs. 1c and 1d).
Next, we investigate the stability of our method when computing very low Reynolds number conditions. Thus, we perform a systematic study, running a set of distinct scenarios increasing the fluid viscosity ν. Fig. 2a and Fig. 2b illustrate the outcomes obtained using both approaches, in comparison with the analytical prediction of Stokes' law. Our results indicate that both methods yield a reasonable description of the settling process and the particle reaches a viscous-dependent terminal velocity v ′ S at the long-time limit. As expected, the results obtained using both methods collapse in a single curve using v S = 2 9µ (ρ p − ρ f )gR 2 and τ = 2 9µ ρ p R 2 . It is noticeable that both approaches slightly fails to reproduce the analytical solution during the initial acceleration process (see the insets in Figs. 2c and 2d).
However, assuming a negligible impact of the finite size boundary conditions implemented (in line with the results introduced in [12]), the present approach compares  Note that a more accurate asymptotic limit velocity is reached in our method. All solid lines represent the expected Stokesian dynamics.
better with the asymptotic limit behavior. We speculate that the difference arises from the impact of dynamic pressure gradients acting on the particle surface for the original particle discretization. Hence, when a volumetric correction is applied to compute the pressure gradients to calculate the force F i, f , we account for the local solid fraction, resolving the corresponding pressure gradient more accurately.

Neighboring walls effect
Complementarily, we examine the impact of the confinement on the dynamics of a settling particle. Specifically, Figure 3. a) Settling velocity of a sphere between two planes scaled by the terminal velocity in a large system when varying the distance between walls H for a liquid with ν = 1000 cS t. The solid red line represents the inverse of the denominator of Equation 6. b) When analyzing the same velocities ratios of a) varying the viscosity, we observe a similar pattern for the confinement effect on the terminal velocity.
we explore the behavior of a particle located between two parallel walls. The particle is centered in the xy-plane of a box with dimensions 10d × Hd × 25d and is left to settle from rest. We systematically vary the distance between the walls H/d = [1.2, 1.4, 1.8, 2.6], in terms of particle diameter. Fig 3a shows the terminal velocity of the sphere moving between two walls and obtained for different confinement conditions, fixing the fluid viscosity ν = 1000 cS t. For simplicity sake, the data values are rescaled using the terminal velocity v ′ S that corresponds to the infinite system size limit. We observe that the simulations can reproduce a decrease in the terminal velocity when reducing the distance between the confining walls. Although, the sphere settles with a velocity larger than the analytical solutions provided by Faxén [11] and Ganatos et al. [12]. Meanwhile, Fig 3b shows that our implementation has a consistent behavior for the range of viscosities investigated here. A similar result was previously obtained for a sphere in a creeping flow condition using an incompressible method based on OpenFOAM software [13].

Conclusions
We introduce an alternative methodology to solve the particle-fluid interaction in the resolved CF-DEM®coupling framework. The new approach explicitly accounts for the body force acting on the fluid phase, when computing the local momentum balance equations. The proper representation of the pressure profile in the fluid phase led to a single force model that no longer requires an explicit buoyant force to account for the gravitational interaction. Our method was able to outperform the cfdemSolverIB method for the settling velocity of a sphere, exhibiting more accurate results.
Complementarily, we study the dependence of the settling velocity of a sphere under confined conditions. The results indicated a consistent velocity decrease when reducing the distance between confining walls. Nevertheless, our numerical outcomes did not reproduce the analytical solutions found in the literature.