Assessing the failure of continuum formula for solid-solid drag force using discrete element method in large size ratios

In loose or moderately-dense particle mixtures, the contact forces between particles due to successive collisions create average volumetric solid-solid drag force between different granular phases (of different particle sizes). The derivation of the mathematical formula for this drag force is based on the homogeneity of mixture within the calculational control volume. This assumption especially fails when the size ratio of particles grows to a large value of 10 or greater. The size-driven inhomogeneity is responsible to the deviation of intergranular force from the continuum formula. In this paper, we have implemented discrete element method (DEM) simulations to obtain the volumetric mean force exchanged between the granular phases with the size ratios greater than 10. First, the force is calculated directly from DEM averaged over a proper time window. Second, the continuum formula is applied to calculate the drag forces using the DEM quantities. We have shown the two volumetric forces are in good agreement as long as the homogeneity condition is maintained. However, the relative motion of larger particles in a cloud of finer particles imposes the inhomogeneous distribution of finer particles around the larger ones. We have presented correction factors to the volumetric force from continuum formula.


Introduction
The solid-solid drag force plays a key role in the hydrodynamics of gas-solid systems such as fluidized beds [1].This force along with the gas-solid drag force determine how the particle clusters form, merge or break at any instant.The clusters are responsible in hydrodynamics as well as heat and mass transfer processes.In computer simulations of gas-solid systems using Eulerian approach, the solid-solid drag force is represented by a continuum formula known as the Syamlal equation [2].This formula is derived based on the assumption of homogeneous mixtures of granular phases.In other words, it is assumed that the distribution of particles within any local volume remains uniform.This assumption could be precise enough when there is monodisperse packing, however, any large size ratio between the granular phases can violate the homogeneity assumption used in the Syamlal equation.Obviously, the large size ratio is associated with the presence of big particles in a cloud of small particles, which could normally result in the inhomogeneous distribution of cloud particles around any big particle.Then the distribution of solid volume fraction (of small particles) and the collision forces on the big particle would be dependent of the relative motion of the big particle with respect to the particle cloud.Such anisotropy in the distribution of collisional forces causes the deviation of mean collisional forces from the predicted values by the Syamlal equation.The focus of this work is to present a correction method of the relevant factors such as the solid volume fraction in terms of the direction of the relative motion of the two granular phases, which consequently lead to a corrected value of the solid-solid drag force.

Methods
In this study, the Lagrangian approach of discrete element method (DEM) is employed.The details of this method can be found in Refs.[3][4].In DEM, the instantaneous positions, velocities and contact forces are calculated for each particle at successive time steps.This helps to calculate the averaged values including the interaction forces.The particle contact is modelled by the Hertzian contact theory where the contact forces are derived based on the spring, dashpot and frictional slider assembly between two overlapping particles.The normal and tangential stiffness coefficients are expressed as functions of particles diameters, their Young's moduli (or shear moduli) and Poisson's ratios [3].The forces at any contact point include the forces in normal and tangential directions.A particle can have multiple contacts at the same time, therefore, the net contact force and its torque are the summation of all contact forces and torques.The net force and torque on any particle determine the translational and rotational accelerations within the corresponding time step.Then the velocities (translational and rotational) are obtained by integrating the accelerations.Similarly, the positions are calculated by the integration of the velocities within the time step.These calculations are performed in each time step for all particles in the system and the trajectory information of particles are stored for further calculations and averaging.The averaged contact force within a volume will result in volumetric drag force between the two granular phases, which can be calculated from DEM simulations as, , .
) ' ( where the mass of particle involved in collision i is mi with the post-collisional and pre-collisional velocities of v´i and vi, where the total collisions of Nml occurs within the volume 'V and time Gt.On the other hand, Syamlal [2] derived the following equation to calculate the volumetric force by knowing the solid volume fraction, velocity and particle diameter of each granular phase along with the material properties including the coefficient of restitution and friction coefficient,

H H H H
Here, Hg is the volume fraction of the void space, which is 1-Hb-Hs.The simulation setup includes a rectangular box within which certain number of particles are randomly positioned to provide the desired initial volume fraction of the small particle size.The second granular phase includes particles with a bigger size.There is either a single big particle or a cluster of them up to 20 particles.Since the size ratio of particles is over 10, the number of particles in each size considerably differs.Figure 1a demonstrates a typical configuration initially in which the volume fraction of small particles (grey) is 0.08.The red particle is the big particle that is initially located at the top of the system.Its initial velocity is given downward along z-direction (here, it is 0.1m/s), as Fig. 1b shows its position after a short time.Only a slice of small particles is represented with the thickness of 0.014m, where the entire width is 0.04m.Here, the diameter of big particle (red) is 12mm, where it is 1mm for the small particles (grey).

Force from DEM
The mean force exerted on the big particle is calculated via Eq.( 1).The instantaneous force is displayed in Fig. 2a.As revealed in Fig. 2b, the mean force exhibits reduced spikes as the averaging given by Eq. ( 1) smooths out the ones observed in Fig. 2a.

a) b)
Figure 2. a) Instantaneous force on big particle vs. time reported every 0.5Ps.b) Time-averaged force on big particle vs. time using Eq. ( 1) with Gt=2.5ms.
Note that the force obtained from Eq. ( 1) is a volumetric measure, which raises an important question on how it can be linked to any continuum formulation such as Eq. ( 2).In Fig. 2b, the volume for averaging is taken out, so it measures the time-averaged force only on a single particle.

Continuum force
In case of single big particle, the interpretation of DEM force (Eq. 1) in a control volume is fuzzy.So, we implement Eqs.(2)(3)(4) to estimate the force on a single big particle as, where Vb is the volume of the big particle as S db 3 /6.Also, Hb is set to a very small value in Eq. ( 4). Figure 3 depicts the variation of the force on big particle Fb versus the relative velocity between the granular phases.The thin solid line stands for the time-averaged force from DEM (FDEM.'V)with the averaging time window of 2.5 ms, which fluctuates and the fitted line is shown by the dash line represented by the exponential function described in the caption of the figure.The thick solid line represents the continuum formula as described by Eqs.(2-4) with the parameters shown in Table 1.Note that 'usb=|us-ub|., where a=8.6x10 -5 , b=0.1465, c=1.3x10 -8 , d=4.8.
There is a considerable discrepancy between the trends of the deterioration of drag force Fb with the relative velocity of granular phases corresponding to the DEM and continuum approaches.The continuum formula demonstrates a smooth decay of Fb with the relative velocity, while the DEM shows a sharp decay upon the start of the particle movement inside the cloud of small particles.Then the attenuation of Fb takes a gentle slope.The two approaches show the same starting points, which indicates the initial homogeneity in DEM setup yields the same force predicted by the continuum formula.However, the inhomogeneity driven by the size difference shows its effect in steepening the slope of force-velocity curve immediately.

Correction to continuum prediction
In order to create any correction scheme to the continuum prediction of solid-solid drag force, depicted in Fig. 3, one must first quantify the inhomogeneity due to the size ratio.Figure 4 shows the vertical profile of the solid volume fraction of the granular phase of small size particles at the moment the big particle has penetrated to a middle position.Note that the force variation is identical for both predictions by the DEM and the continuum formula, however, the initial sharp drop of the force is a transitional event in which the granular pattern develops around the big particle as started from a homogeneous state.

Summary and conclusions
Continuum formulation for solid-solid drag force is vital in Eulerian simulations of particulate flows, which can be found in gas-solid multiphase systems.The key assumption in the derivation of the continuum formula is the homogeneity of the mixture, however, this assumption fails when the size ratio between granular phases grows beyond 10.In this study, we employed DEM simulations to computationally measure the solidsolid drag force and test the performance of the continuum formula.It was observed that, 1) if the system is homogeneous, the predictions of drag are fairly close, 2) the size difference imposes a certain region within which the solid volume fraction obtains a profile indicating the inhomogeneous surrounding of the particle.The region has advection within the medium while the volume fraction profile remains fixed around the intruder.By obtaining the effective value of solid volume fraction, we found out that it can be used in the continuum formula to result the same trend of drag force variation with the penetration velocity.However, the correcting solid volume fraction requires a thorough investigation on various size ratios in a variety of solid volume fraction of the bulk.This is currently pursued as future study in continuation of the current study.Ultimately, a correction table will be created in which the solid volume fraction (e.g., the homogeneous value in a computational cell) is regulated depending its bulk value and the size ratio of particles.Single big particle is studied in the first step though clusters of big particles will be investigated in next steps.
which u is the Eulerian velocity within the given volume with the subscripts s and b standing for the small and big particles.The drag coefficient is expressed by, e and Pf represent the coefficient of restitution and the coefficient of friction for the particles.Other parameters include solid volume fraction H, material density U, and particle diameter d with the same subscripts for granular phases described above.The radial distribution function at contact g0 is expressed by,

Figure 3 .
Figure 3. Force on big particle versus relative velocity between granular phases.The thick solid line represents the continuum formula whereas the thin solid line displays the time-averaged force from DEM with Gt=2.5ms.The dash line is the exponential fit as: sb sb u d u b b e c e a F ' ' . .

Figure 4 .
Figure 4. Profile of solid volume fraction of small grains when big particle is at z=0.070m.The area corresponding to each slice is 0.02x0.02m.The big particle is seen at the center of this slice.The affected region is marked by dash lines that gives an average solid volume fraction of 0.045.The affected region marked by the dashed ellipse in Fig. 4 is persistent which has advection in the direction of the arrow.Therefore, the simplest way to correct the performance of the continuum formula is to modify the solid volume fraction within the affected region.The modified value of Hs can be used in Eqs.(2-4) to obtain new prediction of Fb. Figure 5 illustrates the forcerelative velocity diagram after correcting Hs to the effective value of 0.045.The dash line represents the continuum formula prediction after correction.

Figure 5 .
Figure 5. Modified continuum prediction of force on big particle shown by the dash line.The dots stand for the timeaveraged values from DEM results and solid line is the exponential fit also shown in Fig.3.