The influence of the void fraction on the particle migration: A coupled computational fluid dynamics–discrete element method study about drag force correlations

Granular soils subjected to flow through their soil skeleton can show a behaviour in which fine particles migrate through the pore space between coarser particles. This process is called internal instability or suffusion. This contribution deals with the numerical analysis of the migration of fine particles in a soil column subjected to fluid flow with unresolved coupled computational fluid dynamics–discrete element method (CFD–DEM) with special regards to the used drag force correlation. The contribution investigates the influence of the Schiller–Naumann model and its extension with a voidage term on the migration behaviour of fine particles. The voidage term is further varied with a parameter, which controls the impact of the change of the void fraction on the drag force. It could be observed that the Schiller–Naumann model does not yield in a suffusive behaviour while the extended models show significant particle migration. Thereby, increasing the impact of the void fraction on the drag force results in stronger particle migration. These results reveal the need for good validation techniques. They indicate how the drag force correlation can be adapted to depict the correct particle migration behaviour.


| INTRODUCTION
Suffusion is a process in a cohesionless gap-graded soil where a fraction of the fine particles does not receive the effective stresses and can move freely between the coarser fraction. The channels between the coarse particles have to be broad enough for suffusion to occur. Suffusion in soil structures, such as earth dams and dykes, is a serious problem for long-term stability. Suffusion increases the permeability and porosity of the soil while the bulk density decreases, which leads to a minor overall resistance of the soil. A variety of filter criteria is proposed to assess suffusion. Semar et al. 1 Abbreviations: CFD, computational fluid dynamics; DEM, discrete element method; DNS, direct numerical simulation; IB, immersed boundary; LBM, Lattice-Boltzmann method. evaluate these criteria proposed in the literature. They conclude that these approaches are limited in their usability as they are mainly valid for soils, which are comparable with those analysed. Furthermore, local effects and structural changes, which can lead to settlements or influence the hydrodynamic conditions, are entirely neglected. Hence, the need for further research is very high. In general, suffusion has been widely studied at a macroscopic scale. These efforts improved our understanding of the macroscopic effects of suffusion, but the mechanisms at the microscopic scale associated with the initiation and evolution of suffusion remain unclear.
Numerical studies can, therefore, give valuable insight into the leading mechanisms of suffusion and its consequences on the soil matrix. With increasing computational power, the numerical methods enabling the computation at a microscale level gain more and more importance, and the interest in both research and industry increases continuously. These numerical methods allow the investigation of effects, which cannot be depicted by pure continuum approaches, for example, the interaction between the particles and fluid-particle interactions. For the latter, numerical methods are necessary that resolve fluid flow and particle motion. One of these methods is the coupled computational fluid dynamics (CFD)-discrete element method (DEM) (CFD-DEM), where the fluid motion is computed in a Eulerian way with the CFD and the particle motion is tracked in a Lagrangian way with the DEM. In the CFD, a mesh is required, which is usually composed of finite volumes. Within these, the flow field is computed using the volumeaveraged Navier-Stokes equations. Likewise, the velocity of the fluxes is computed at the cell faces and interpolated to the cell centre. The DEM solves the particle motion explicitly at the individual particle level. The particle velocity and position are determined using Newton's second law of motion. Combining these two methods, one needs to define the scale level between the computational cells of the CFD and the particles. In general, there exist two methods to meet this problem: the resolved and the unresolved coupled CFD-DEM (cf. Figure 1). In the framework of the resolved CFD-DEM, the particles are covered by several CFD cells. This enables an accurate computation of the flow field around the particles, but a large number of CFD cells is required. As a result, only small problems with some thousands of particles are feasible concerning the computational costs. To handle larger problems, the unresolved CFD-DEM is the appropriate method. Thereby, the CFD cells are notably larger than the particles in a way that several particles are located in one cell. Thus, problems with some millions of particles can be computed, but one has to deal with a loss of information concerning the fluid flow between the particles. For that reason, fluid-particle interactions are taken into account via transfer coefficients described by closure models. The accuracy of these closure models determines the reliability of the results. Fluid-particle interaction forces consist of lift forces (e.g., Saffmann and Magnus lift force), unsteady forces (Basset history force and virtual mass force) and steady forces like the pressure gradient force, the buoyancy force and the drag force. This contribution will focus on the drag force as it is one of the driving forces in particle-fluid systems.
An analytical formulation of the drag force in a particle assembly does not exist. The drag force for particle-fluid systems has a strong dependence on the pore structure and the pore level physics, which generally requires it to be estimated experimentally or through the use of empirical correlations. In general, two methods have been used so far to find correlations for the drag force in particle-fluid systems: i. Based on empirical correlations for either pressure bed drop 3,4 or fluidisation 5 ii. Based on numerical simulations at microscale with direct numerical simulation (DNS) or with the Lattice-Boltzmann method (LBM). 6 (A) (B) F I G U R E 1 Resolved (A) and unresolved (B) computational fluid dynamics-discrete element method (CFD-DEM) 2 Li and Kuipers 7 investigated the effect of the voidage function exponent χ on gas-fluidised bed reactors. They found that with increasing χ the bed height of the solid-fluid mixture increased dramatically. The increase of χ provides a stronger drag force acting on the particles. Applying this result onto the suffusion problematic, it is interesting to analyse the increase of χ as the void ratio changes are small yet crucial.
Recently, Tang et al. 8 performed DNS simulations with an immersed boundary (IB) method modelling flow through a static monodisperse spherical bed. They established an accurate correlation for this static problem based on correlations they found combined with correlations for low Re-flows and single sphere flows. However, they found 20% of deviation when performing tests with freely moving particles. Hence, accurate drag correlations for static monodisperse arrays of spheres are not a priori applicable to problems with freely moving particles.
To sum up, there exist a variety of models available in the literature. The question is how good they can be applied to problems like suffusion. Empirical correlations (like Ergun 3 ) are obtained for idealised conditions and do not apply for all situations. Also, the majority of the research regarding drag correlations is focussed on fluidised bed reactors with high velocities and groups of particles moving together. [9][10][11][12] The void fraction changes are significant, which is not the case in suffusion. Additionally, the migration of the individual particles is of interest and not the movement of a group of particles. Recently, numerical studies on suffusion with unresolved coupled CFD-DEM were carried out. Hu et al. 13 investigated seepage flow in gap-graded and well-graded soil samples with unresolved coupled CFD-DEM employing the drag force model of Di Felice. 14 They performed a parametric study on the effect of initial fines content and the hydraulic gradients on the erosion kinetics. Their results showed that the ultimate fines content is independent of the imposed hydraulic gradient once erosion is initiated. However, Hu et al. 13 did not perform a validation of the model, and gravitation on the particles was neglected during seepage flow. Hu et al. 15 investigated then the differences in the mechanical response between suffusion numerically modelled by particle removal, unresolved coupled CFD-DEM and reconstituted samples. Reconstituted samples have a particle size distribution (PSD) of already eroded soil. They first applied a seepage flow on the soil samples and then carried out a triaxial test on the eroded soil. They found that compared with the two other methods, the coupled CFD-DEM simulation resulted in more fines that are likely to erode. The seepage flow through the sample leads to uneven erosion processes in the soil specimen, which cannot be captured by the removal of particles or reconstituted soil samples. This result emphasises the need for accurate simulations with coupled CFD-DEM as the seepage flow clearly influences the overall mechanical response of the investigated soil sample. Hence, there exists an essential need for research regarding the applicability of the existing drag force correlations. Kanitz and Grabe 16 investigated the effect of different drag force correlations on the particle migration in suffusion. Significant differences in both the particle paths and the amount of washed out fines could be found. Often, the existing drag force correlations are a priori validated against the drop of a single particle or fluidised bed and afterwards applied to a very different kind of problem. It is uncertain if the validation against cases for which the drag correlations were actually established lead to an accurate result when applied to a different type of problem.
The goal of this contribution is to investigate further the effects of the drag force correlations on the particle migration in suffusion problems. Therefore, existing correlations are used and varied to scrutinise, how the drag correlations act on the particles. This paper is organised as follows: we will first give a brief introduction into the theoretical background of the drag correlations and introduce the numerical model to investigate suffusion. The first simulations are carried out with the conventional drag force models to identify possible differences. Then, based on the analysis of the drag force model of Di Felice, 14 an extension of the drag force model of Schiller and Naumann 17 is presented with varying values for the voidage function to identify its effect on the erosion of fines in a suffusive particle assembly. This analysis will help to gain a deeper understanding of the influence of the drag force and especially of the voidage function in suffusion problems. We will present the results for the varied models and finish with a conclusion that briefly discusses the major findings. The numerical simulations are carried out with the open-source software package CFDEM ® coupling, which combines the discrete element code LIGGGHTS ® with CFD solvers based on OpenFOAM ® . 18 2 | THEORETICAL BACKGROUND

| Coupled CFD-DEM
The DEM was originally developed by Cundall and Strack. 19 In this method, every particle of the solid phase is modelled explicitly based on Newton's second law of motion. The translational and rotational motion of a spherical particle i with a mass m is calculated by where v i and ω i describe the translational and angular velocity of the particle i. F c ij and M ij represent the contact force and the torque acting on particle i by particle j or the walls. F f i and F g i are the particle-fluid interaction force and the gravitational force. I i is the momentum of inertia of the particle i. The motion of the fluid phase in this coupled approach is solved using the CFD. It is based on the conservation of mass and momentum. Hence, the continuity equation and the volume-averaged Navier-Stokes Equations (NSE) 20 are used to solve the pressure p and the velocity u f of the fluid in the computational cells τ f and ρ f are the fluid viscous stress tensor and the density of the fluid. In Equations (3) and (4) two additional terms can be found that are related to the coupling with the DEM: the void fraction ε and the fluid-particle momentum exchange term R pf . Including the void fraction into the computation of the fluid velocity, the presence of the particles and their effect on the fluid flow is accounted for. With the term R pf , the interaction forces between the fluid and the particulate phase are included in the computation of the fluid velocity and pressure. It is determined by the following equation: with the interaction force coefficient K pf . This coefficient is defined as the superposition of the particle-fluid interaction forces. The interaction forces are composed by the fluid drag exerted on the particles, the buoyancy force and the pressure gradient force as well as unsteady forces. The unsteady forces like the virtual mass force, the Basset history force and lift forces, for example, Saffman and Magnus lift force, are usually very small compared with the drag force and therefore not included in the simulations presented in this contribution. The next section will give an introduction to some widely used drag force correlations. In a coupled CFD-DEM simulation, the void fraction in the computational cells determines the drag force acting on the particles. Thus, it is crucial to accurately determine the void fraction in the individual cells for a realistic calculation of the drag force. Thereby, a void fraction scheme should provide the following criteria: 21 1. conserve the total mass of the particulate phase, 2. predict a grid-independent velocity field, 3. produce a smooth void fraction field within a densely packed bed.
Several void fraction schemes exist, which will be briefly summarised in the following. The simplest void fraction model is the so-called particle-centred method (PCM). Thereby, the entire particle volume is located at the centroid of the particle. 22 However, this simplification results in large fluctuations when the particle centroid crosses the cell boundaries. 21 The PCM should only be used in cases where the fluid cell is large compared with the particle size. Another void fraction scheme is the divided particle volume method (DPVM). In this scheme, the volume of the particle is subdivided into smaller volumes. The centroid of each subvolume is then used to assign the volume of the particle to the corresponding computational cell. In the framework of CFDEM ® coupling, the particle is divided in 29 nonoverlapping regions of equal volume. 23 Furthermore, statistical methods can be used in void fraction schemes. In the statistical kernel method, the solid volume fraction is determined by distributing the particle volume by a weighting function that decays away from the particle centre. 24 In the presented numerical simulations, the DPVM is used.

| Drag force correlations
In general, when fluid flows through an assembly of particles, each particle experiences two types of forces: i. hydrostatic force: buoyancy-type force due to the average pressure gradient around the particle, ii. hydrodynamic force: force resulting from local friction losses, which can be referred to as the drag force.
The reaction force from particles onto fluid results in energy dissipation. As a consequence, the pressure in the fluid reduces while it passes through the bed. The most widely known theory to calculate this pressure drop is Darcy's law: with the macroscopic parameter of the permeability K, the dynamic viscosity μ f and the fluid velocity u. It describes a linear relationship between the pressure gradient Δp and the flow field at the macroscopic level. However, it has been shown that Darcy's law is restricted to creeping flows and only valid for homogeneous, isotropic, unbounded and nondeformable porous media. 25 The drag force for an isolated particle can be determined by where C d is the drag coefficient, d p the particle diameter and v the particle velocity. C d is thereby closely related to the Reynolds number and can be determined through Schiller and Naumann: 17 Re p 1 + 0:15Re 0:687 p � � , for Re p < 1000 0:44, for Re p ≥ 1000: Re p is the particle Reynolds number and defined by Re p = ρ f d p ε|u − v|/μ f . ε is the void fraction of the fluid. Wen and Yu 4 extended the formulation for an isolated particle for the pressure drop in a dilute particle assembly by including a momentum exchange coefficient β between the fluid and the solid phase. They derived a correlation for β by correlating their experimental data for bed pressure drop experiments. This limits its applicability to equilibrium or near-equilibrium states. Gidaspow 26 combined the Ergun equation, which is suitable for dense systems, 3 with the correlation of Wenand Yu, 4 only suitable for dilute systems. At the boundary between the two regimes at ε = 0.8, Gidaspow 26 adopted a step change, which was not elegant from a numerical point of view. 27 To overcome this problem, Di Felice 14 proposed a continuous single-function correlation. 9 He investigated how the drag force acting on a particle assembly needs to be modified from the drag force on a single particle using the drag coefficient of Dallavalle. 28 Di Felice 14 introduced a voidage function f(ε). Evaluating f (ε) for fluidised bed and fixed bed experiments, Di Felice 14 derived a function for the parameter χ, which depends solely on the Reynolds number. An overview of these mostly used drag correlations based on empirical correlations can be found in Table 1.

| NUMERICAL SETUP
The influence of the different drag force correlations on suffusion is investigated in a small numerical setup. The details of the setup can be depicted in Figure 2

| Geometry and initial conditions
Two kinds of materials are used for the investigation: glass beads with a diameter of 5 mm and glass beads with a diameter of 1 mm. In this way, it is ensured that the packing is suffusive according to the filter criterion of Burenkova. 29 The T A B L E 1 Some frequently used drag force correlations Re p 1 + 0:15Re 0:687 Re p 1 + 0:15Re 0:687 when ε < 0.8: Note: It is important to clarify that u − v is the relative and not the superficial velocity. ϕ denotes the solids fraction.

F I G U R E 2
Numerical setup for the investigation of suffusion (left) and mesh-to-particle ratio (right) [Colour figure can be viewed at wileyonlinelibrary.com] height of the column is 0.1 m with a diameter of 0.08 m. The diameter is taken from a suffusion test setup used in the laboratory. This diameter ensures that clogging of the coarse particles between the outer walls can be omitted. Thus, an oedometric stress state can be assumed at the beginning of the simulation. The different glass beads are initialised in a layered way (c.f. Figure 2). In this way, it can be monitored how many particles mitigate into the upper and consequently coarser layer of 5-mm particles. Additionally, the void fraction changes are more pronounced in this way, which allows investigating the influence of significant changes in the void fraction on the drag force models. The particles are initialised in a pure DEM simulation, without the coupling to CFD. This lowers the computational costs and does not influence the subsequently coupled simulation. After the creation of the particles, the particles are allowed to settle until the kinetic energy in the simulation is lower than 5 Á 10 −7 Nm. The material parameters were calibrated against material tests and are listed in Table 2. The epsd2 rolling friction model implemented in LIGGGHTS ® proposed by Ai et al. 30 is employed to account for the rolling resistance of the particles. The simulated fluid has the properties of water at a temperature of 20 C. Its density is ρ = 1000 kg/m 3 and the kinematic viscosity is ν f = 1.0 Á 10 −6 m 2 /s.

| Boundary conditions
A velocity boundary condition is initialised at the bottom with a vertical velocity of U z = 0.035 m/s. The pressure boundary condition at the bottom is a zero-gradient boundary. At the top of the column, the velocity boundary condition is zero-gradient, and the pressure boundary condition is set to 0. A slip boundary condition is initialised at the walls. At the bottom of the column, an additional layer of CFD cells without particles is initialised. There, the void fraction is equal to 1. In this way, it is ensured that the applied velocity boundary is the superficial velocity.

| Coupled simulation with CFD-DEM
The coupling interval between the DEM and the CFD is 100, for example, after 100 DEM time steps, the particle information are passed to the CFD solver and the NSE are solved. The DEM time step is t DEM = 1 × 10 −5 s, and consequently, the CFD time step is t CFD = 1 × 10 −3 s. To ensure the stability of the numerical calculation, it is checked in each time step whether the Courant-Friedrich-Lewy (CFL) number is below the selected tolerance of 1. The divided void fraction scheme is used for the determination of the void fraction in the computational cells as it has been widely tested in applications where the particle diameter is in the range of the CFD cell size. 23 The simulated time is 10 s.

| Mesh sensitivity analysis
The mesh size is 1.6 × 2.0 × 1.0 cm. In coupled CFD-DEM simulations, the volume of the particles in the individual CFD cells must be sufficiently large that the average value does not depend on the size of the cell. The limiting volume has to contain a sufficient amount of particles such that variations are maintained below a minimum value, so the flow properties within the volume can be regarded as point values. 31 When it comes to the NSE, the question is how large the volume of the CFD cells has to be to ensure that the variations of the value for K pf are minimal. Pirker et al. 32 and Marshall and Sala 33 suggested 3 Á d p as the lower boundary for the mesh grid size. The mesh-to-particle-ratio is depicted in Figure 2 (right). A mesh sensitivity analysis was carried out to verify if the chosen discretisation leads to accurate results without computational errors. Figure 3 shows the numerical results after 1 s with a coarse (a), a medium (b), and a fine (c) discretisation. It can be observed that the positions of the fine particles vary using a different degree of discretisation. The amount of eroded fines into the upper coarse particle matrix is presented in Figure 4 (right). While the coarse and the fine mesh lead to the same number of eroded fines after 0.8 s, more fines are eroded using a medium discretisation. The pressure drop across the bed is investigated and plotted against the Ergun equation to verify whether the analysed discretisations lead to erroneous results. Therefore, the simulation is carried out with fixed particles via a 'freeze'-command 34 to prevent the particles from moving. Mobile particles lead to deviations from the Ergun pressure drop. 8 Figure 4 (left) presents the pressure drop Δp across the bed with increasing fluid velocity u. It shows that both the medium and the fine mesh result in significant fluctuations. Additionally, the medium discretisation leads to a strong deviation from the Ergun equation. Hence, the amount of eroded fines is overestimated using a medium discretisation. The coarse discretisation results in a smooth pressure drop curve representing the Ergun equation. Therefore, the presented simulations of suffusion are carried out with the coarse discretisation.

| RESULTS AND DISCUSSION
First, the erosion of the fine particles into the upper coarse particle matrix using the conventional models is analysed in the following. To further investigate the influence of the drag force on the erosion of the fines, an extended version of the Schiller and Naumann model 17 with a voidage function is introduced. This enables the analysis of the influence of the void fraction term and accordingly the drag force on suffusion.

| Migration behaviour in suffusion using the conventional drag force models
The migration of the fines using the drag force correlation of Di Felice, 14 Gidaspow, 26 Koch and Hill 6 and Schiller and Naumann 17 is presented in Figure 5.
The particle distribution after 10 s of seepage flow using the drag force model of Di Felice, 14 Gidaspow and 26 Koch and Hill 6 shows slight differences. Using the models of Di Felice, Gidaspow and Koch and Hill, we can observe groups of particles that are washed into the upper coarse particle matrix layer and single particles. The particle groups begin to clog and are not migrating further into the column. The single particles, however, are eroded significantly further upwards in the upper coarse layer. The drag force model of Schiller and Naumann results in significantly reduced erosion of fines. No single eroded particles migrate into the upper coarse particle matrix. Only by investigating the amount of fines, which is washed into the coarse upper layer during seepage flow (see Figure 6), we can identify deviations in the numerical results using the different drag force models. While the drag force model of Schiller and Naumann leads to an erosion of ≈1000 fine particles, the amount of fines using the drag force model of Di Felice, Gidaspow and Koch and Hill is significantly higher. Thereby, the model of Koch and Hill results in the lowest erosion rate among those three. The erosion is smooth, whereby the most substantial migration is finished after 0.3 s. Afterwards, clogging of the particles occurs, and the amount of fines remains constant. The drag force model of Gidaspow results in a higher erosion in comparison with the model of Koch and Hill, but the migration behaviour is comparable. Only after 6.5 s, a blow-out of particles is observable by the sudden increase of eroded fines. The drag force model of Di Felice results in more changes between clogging and blow-outs of the fine particles visible by the step-wise increase of eroded fines after 1.2 s. After 4.5 s, the model of Di Felice yields a constant value of eroded fines, which is the highest among the investigated drag models. In conclusion, these conventional drag models lead to slightly different erosion dynamics in the particle column. The differences are not as significant as observed by Kanitz and Grabe, 16 but visible.

| Extended drag force correlation for suffusion
The problem of suffusion could be classified as a near-equilibrium problem as the particle motion remains quite limited compared with, for example, fluidised bed reactors. Nevertheless, the flow in the narrow channels between the grains plays a crucial role in particle migration. Hence, the direct application of the existing drag correlation without quantifying the influence of the small void fraction changes could be uncertain. Comparing the drag correlations listed in Table 1, it is peculiar that there are two correlations for the drag coefficient C d , whereby the most used is the one from Schiller and Naumann. 17 The only differences then are the various extensions with either a voidage function f(ε) or a different correlation for the momentum exchange coefficient β. It has been found that β mainly influences the pressure drop in the bed while for the particle motion, the volume fraction has a subtle but complex influence on the drag force. 10 In the drag force model of Di Felice, 14 the voidage function is dependent on the Reynolds particle number Re p . Figure 7 depicts the evolution of Re p over the column height and the according value for (− χ + 2) in the formulation of The investigations of the conventional drag force models show that the drag model of Schiller and Naumann 17 does not lead to a significant suffusive behaviour of the particles. Accordingly, this model can serve as a starting point for further investigations. To analyse the influence of the void fraction and the drag force on the particle migration behaviour in a suffusive column, the model of Schiller and Naumann 17 is extended with a voidage function. As the investigations of Re p and the exponent χ in the voidage function of Di Felice 14 reveal, the influence of Re p on the exponent of the voidage function is minor in the investigated range of Re p . This allows to simplify the voidage function with a constant exponent instead of an exponent dependent on Re p . Hence, to analyse the influence of the value of χ on the erosion behaviour of the fine particles, the voidage function is formulated as f ðεÞ = ε − χ with χ 2 f0:75,1:0, 1:25, 1:5g: The influence of the void fraction on the drag force is consequently varied with a constant factor χ. The extended drag force correlation for suffusion concludes to the following expression: With increasing χ, the void fraction has an increasing influence on the drag force and the drag force is enhanced. It is important to mention that if χ < 0, the presence of the particles reduces the particle-fluid interaction forces, which would lead to wrong physical behaviour. Hence, χ needs to be positive. To verify if the presented extension of the drag force model of Schiller and Naumann lead to a correct physical behaviour of the particles, a submerged onedimensional (1D) consolidation test is carried out according to the consolidation test setup presented in Zhao and Shan. 35 They showed that the numerical results using the drag force model of Di Felice 14 correspond to the 1D consolidation theory by Terzaghi. 36 Thereby, particles of 1-mm diameter are first settling under gravity and buoyancy forces and then consolidated with a surcharge load of 100 Pa. Figure 8 depicts the settling and consolidation process of the particles using the drag force model of Di Felice, Schiller and Naumann and the presented extension of the Schiller and Naumann drag model with χ = 0.75,1.0,1.25 and 1.5. The different drag force models result in a comparable response of the particle velocity both in the settlement and in the consolidation phase. Figure 9 depicts the overall erosion of fines in the upper coarse particle matrix using the model of Schiller and Naumann, the model of Di Felice and the extended drag force model of Schiller and Naumann with varying χ. It is shown that the model of Schiller and Naumann results in the lowest erosion with around 1000 particles. After 4 s, the erosion process is finished, and clogging occurs. This behaviour changes significantly by using the extended drag force models, where the major erosion finishes earlier. Table 3 depicts the moment when clogging occurs and the major erosion finishes. In general, increasing the drag on the particles leads to a faster clogging of the particles. The same applies to the number of eroded fines, which increases with the increase of the drag force. For the models with χ = 0.75 and χ For an in-depth analysis of the erosion processes into the upper coarse particle layer, Figure 10 presents the amount of eroded fines per 0.5-cm height in the upper coarse particle matrix over time for the extended drag models with varying values for χ. This enables the analysis of the clogging process and the erosion distance of the fine particles into the particle column.

| Influence on the particle migration and particle loss
In the first layer above the fine particle layer (z = 0.055 m to 0.06 m), a small peak followed by a plateau can be observed for all values of χ. The amount of fines in this layer is similar for χ = 0.75 and χ = 1.0. The drag force models with χ = 1.25 and χ = 1.5 result in a significantly higher amount of eroded fines. The peak at the beginning of the seepage flow signifies further migration of the particles into the overlying layers. The models with χ = 0.75 and χ = 1.0 show a clogging behaviour as the amount of particles stays constant while the amount of fines increases slightly using χ = 1.25 and χ = 1.5.
In the overlying layer from z = 0.06 to 0.065 m, the first peak is more elaborate for all investigated models. Using χ = 1.5, it can be observed that this peak is followed by an increase after 2 s of seepage flow. This increase shows a blowout of particles from the fine layer as this increase occurs at the same time in the underlying layer, but less pronounced. Hence, the increased drag force by using χ = 1.5 results in a fast erosion of the fines from the fine layer into this layer In the next layer from z = 0.065 to 0.07 m, the amount of fines is similar for all investigated values of χ. Clogging occurs in this layer, and no significant increases are observable.
In the layer of z = 0.07 to 0.075 m, it is striking that despite the decreased drag force, the amount of eroded fines is higher for χ = 1.0 compared with χ = 1.5. The drag force models with χ = 0.75 and χ = 1.25 yield a similar amount of fines in this layer. Because of the increased travel distance for the fines to reach this layer, the peak widens and increases to 1.5 s.
The overlying layer from a height of z = 0.075 to 0.08 m shows a different erosion behaviour compared with the underlying layers. While clogging occurs for the drag force models with χ = 1.0 and χ = 1.25, visible by the constant amount of fines after 2 s, the amount of fines increases throughout the simulation using χ = 0.75. The steady increase is a result of the erosion of single particles into this layer. As a result, χ = 0.75 leads to a higher amount of fines in this layer compared with the other values for χ. A blow-out of fines from the underlying layers can be observed for the drag force model using χ = 1.5 after 2 s of seepage flow indicated by the sudden increase of fines.
In the next layer from z = 0.08 to 0.085 m, the drag force models with χ = 0.75, χ = 1.25 and χ = 1.5 reach a constant value while the drag force model with χ = 1.0 results in a steady increase of eroded fines. This drag force model thus leads to the highest amount of fines in this layer compared with the other values for χ. The increase rate is strongest until 2 s of seepage flow, which corresponds with the decline of the peak in the underlying layer. The steady increase after 2 s signifies that from this moment on, only single particles are eroded into this layer. A sudden increase in the amount of fines corresponds hence to groups of particles migrating into the layer. The amount of fines using the drag force model with χ = 0.75 reduces significantly in comparison with the underlying layer. The same observation applies for the drag model with χ = 1.5. However, the analysis of the amount of fines in the next layer reveals different causes for this decreased number of fines. While no erosion of fines into the layer of z = 0.085 to 0.09 m occurs using the drag force model with χ = 0.75, the amount of fines using χ = 1.5 shows a significant increase. Hence, the reduced amount of fines in the underlying layer for χ = 1.5 is a result of partially clogging particles accompanied by a steady migration of fines into the overlying layer. The drag force using χ = 0.75 is not strong enough for the fines for travelling further upwards, so they assemble in the underlying layers. Accordingly, also, the amount of fines using χ = 1.0 is reduced significantly. The drag force model with χ = 1.25 results in a similar amount of fines compared with the drag force model with χ = 1.5. However, the particles require more time to travel into this layer in comparison with the results with χ = 1.5. It is interesting to notice that the amount of fines using χ = 1.5 reduces after 8 s, which corresponds to an increase of fines in the underlying layer. Hence, the particles are moving at the boundary of these layers. No clogging of the particles occurs, but they move freely in the coarse pore space.
To sum up, the following conclusions can be drawn from the investigation of the particle migration depending on the use of different drag force models: i. The default Schiller-Naumann drag model does not yield a significant suffusive behaviour of the fine particles. ii. Increasing χ leads to a stronger migration into the upper layer. iii. Overall clogging occurs earlier with increasing drag. iv. Groups of particles assemble in lower layers; single particles are washed further into the column. v. Increasing χ and accordingly the drag force leads to more blow-outs of fines. vi. Increasing χ leads to an increased erosion distance of the fines. However, the results of the drag force model with χ = 1.25 and χ = 1.5 show similar results. There is hence an upper boundary in regard to the maximum erosion distance. Figure 11 shows the pressure, the fluid velocity, the void fraction and the interaction force coefficient K pf evolution over the soil column after 1 s of simulated time for all used drag force models. By investigating the particle migration over time in Section 4.3, it could be observed that after 1 s, the majority of particle movement has finished. For this reason, the flow field is investigated at this time step. While at the inlet the fluid velocity, the void fraction and K pf seem to be independent of the used drag force correlation, the magnitude of the fluid pressure at the inlet is higher with increasing χ. Likewise, the usage of the extended drag force model with χ = 1.5 leads to an inlet pressure of 1.55 kPa while the default Schiller and Naumann model has an inlet pressure of 0.35 kPa. However, this increase has no significant influence on the fluid flow or the void fraction at the inlet. With increasing height, we can observe a continuous decrease of the pressure and a slight increase of K pf . Thereby, the pressure using the model of Di Felice decreases in the same order as with the extended model with χ = 1.25. The fluid velocity and the void fraction fluctuate moderately caused by the fine particles, which have fallen into the coarse bottom layer. Just before entering the fine particle layer (highlighted in light green), the void fraction decreases in general. When entering the fine particle layer, the interaction between the particles and the fluid increases strongly, which can be observed by the sudden increase of K pf . Thereby, K pf increases with increasing χ. This interaction and, consequently, the transformation of kinetic energy of the fluid onto the particles result in a substantial decrease of the fluid velocity despite the decrease of void fraction in this layer. It is interesting to notice that the decrease of the fluid velocity is less significant using the extended model with χ = 1 and χ = 1.25 in comparison with the other drag force models. Also, the pressure decreases significantly. The subsequent layer, highlighted in light blue, is the so-called 'migration layer', into where the fine particles travel. In this layer, the velocity decreases again accompanied by a decrease of the pressure. The influence of the particles which are washed into this layer results in a substantial decrease of void fraction. The interaction between fluid and particles increases in this layer in comparison with the fine layer indicated by the increased value for K pf . The difference in the magnitude of K pf using different values for χ is striking. This high magnitude of K pf with increasing χ explains the increasing travel distance of the particles. Due to augmenting χ, more energy of the fluid is transferred to the particles. The concurrent higher reaction force of the particles also adds up to the loss of pressure, which is observed to be higher with increasing χ. As a consequence, the void fraction above the 'migration' layer still continues to decrease for χ = 1.5. The amount of fines in the 'migration' layer is generally higher compared with the other drag force models. The void fraction using the other drag force models remains constant except for the extended model with χ = 0.75, where the void fraction begins to increase compared with the migration layer as a result of the reduced migration distance. At the top of the soil column, the influence on the pressure and K pf of the used drag force models on the fluid flow vanishes. The differences in the void fraction of the different drag force models results in slight differences of the fluid velocity at the top of the column.

| Influence on the flow field
In summary, the following conclusions can be made by investigating the flow field over the soil column height: i. Increasing χ leads to higher initial fluid pressure.
ii. Increasing χ results in an increased magnitude of K pf and consequently the increase of the reaction force of the particles and their movement. iii. The pressure drop due to the movement of the particles increases with increasing χ. iv. The freely moving particles lead to a high amount of kinetic energy loss of the fluid. This results in a decrease in the fluid velocity, although the void fraction decreases.

| Drag force acting on the particles
To investigate the course of the drag force on the fine fraction of the particles, the drag force on several fine particles was evaluated and the median value calculated. The evolution of the drag force can be seen in Figure 12.
The extension of the drag force model by a void fraction term increases the magnitude of the drag force. The Schiller and Naumann drag force remains constant at 0.5 mN during the whole simulation. The magnitude of the drag force when using χ = 0.75 is increased by a factor of 2 compared with Schiller and Naumann. This factor reflects the increase caused by the void fraction term with χ = 0.75. In contrast to the model of Schiller and Naumann, a slight increase of the drag force at the beginning of the seepage flow is observable. The drag force remains constant during the rest of the simulation. The magnitude of the drag force with χ = 1.0 fluctuates slightly, and the increase at the beginning of the simulation is greater. The drag force reacts to changes in the void fraction. Because its data were evaluated directly at the particles, these fluctuations indicate particle migrations. Accordingly, either the sampled particle itself has migrated into another computational cell or other particles have migrated into its cell, which is why the void fraction and thus the amount of drag force changes. The same applies for χ = 1.25. The fluctuations are more elaborate as the sensitivity for void fraction changes increases with increasing χ. In contrast to χ = 1 and χ = 1.25, the amount of drag force does not fluctuate as much with χ = 1.5. In this model, however, a sudden increase of drag force up to 0.5 s can be detected. After a slight decline, the drag remains constant, with a small increase at 4.5 s. The drag force model of Di Felice results in a smooth evolution of the drag force during the simulation. Like the extended drag force model, the model of Di Felice increases at the beginning of the simulation. However, fluctuations in the drag force cannot be identified. In summary, a higher tendency towards fluctuations of the drag force can be observed with increasing χ. As the particle movements are stronger, they also change the void fraction in the computational cells. Additionally, with increasing χ, the drag force is more sensitive to changes in the void fraction. Overall, the increase of drag force at the beginning of the simulation increases with a greater drag force. Figure 13 depicts the normal contact forces between the particles at the beginning and at the end of the simulation for all used drag force models. The resulting forces are illustrated as tubes, and the force scales their radius. Hence, greater forces result in a greater radius of the tube, and when a loss of contact occurs, the tubes are no longer visible. The contact force F n is computed using the following equation:

| Influence on the normal contact forces
with K n = 2E c r i r j where K n represents the stiffness of the particles, δ n is the overlap of the contacting particles, E c is the stiffness modulus and r i and r j are the radii of particle i and j. Likewise, the magnitude of the contact force is dependent on the radius of the particle. For this reason, the contact between the coarse particles has a larger magnitude than the contact of the fine particles in Figure 13. At the beginning of the seepage flow, the normal contact force decreases with increasing height caused by the decreasing weight of the overlying particles. The different layers can be well recognised by the different radii of the tubes. After 10 s of seepage flow, Figure 13 clearly shows the influence of the usage of the different drag force models on the particle contact forces.
While the normal contact forces in the bottom part using the Schiller and Naumann drag model and the extended drag models with χ = 0.75 and χ = 1 are similar, with a further increase of χ the contact forces have reduced significantly. The results with χ = 1.5 even show a significant reduction of the contact between the bottom part of the column and the fine particle layer.
In general, the usage of the extended drag models shows a trend towards the loss of internal stability in the bottom part of the column accompanied by a gain of internal stability in the top part of the column due to the F I G U R E 1 3 Normal contact forces in soil column at the beginning of the seepage flow and after 10 s for the Schiller and Naumann and Di Felice drag model and the extended drag models with varying χ [Colour figure can be viewed at wileyonlinelibrary.com] fine particle migration into that layer. This can be identified by the denser contact network especially at the boundary between the fine particle layer and the upper coarse particle layer. This trend increases with increasing χ as the fine particles migrate further into the upper layer and the contact network intensifies. The consolidation of the particles explains the high contact forces in the bottom part using the Schiller and Naumann drag force model and the extended drag models with χ = 0.75 and χ = 1. The fluid force is lower than the weight of the particles. Hence, the fluid force is not majorly reducing the normal particle-particle contact forces. In Section 4.4, it could be concluded that the fluid pressure increases with increasing χ. Accordingly, the effective stresses between the particles are reduced significantly, resulting in a loss of contact between the particles in the bottom part for χ = 1.5. In the top part, the pressure field does not vary significantly from each other. Therefore, the leading mechanism for the gain of stability is the increase in contact partners due to the higher amount of particles in this area caused by migration. The drag force model of Di Felice results in a similar normal contact force network than the extended drag force model with χ = 1.25.
Likewise, the following results can be concluded by investigating the contact force network resulting from the use of the different drag force models: i. With a higher value of χ, the bottom part loses stability due to an increase of the pressure field. ii. The weight of the particles and the resulting particle-particle contact forces dominate the effect of the pressure increase in the bottom part using the Schiller and Naumann drag model and lower values of χ, iii. The increase of χ leads to a greater migration of particles into the upper layer, which is accompanied by a gain of stability in the top part of the soil column caused by a densification of the contact particle network.

| CONCLUSIONS
The influence of the different drag force models on the suffusion in a soil column was investigated with the unresolved coupled CFD-DEM. A voidage term in the form of ε −χ was introduced to the Schiller and Naumann drag force model to analyse the effect of the change of the void fraction on the migration of fine particles through the soil column. An analysis of the exponent χ in the voidage function introduced by Di Felice 14 showed a saddle point for the relevant values of Re p leading to almost constant values for the exponent χ. This allowed to simplify the exponent of the voidage function as a constant value. Four constant values for χ were hence investigated. The results were compared to the Schiller and Naumann drag force model, which does not include such a voidage term. Additionally, the results were compared with the drag force model of Di Felice, which contains a voidage term with an exponent χ dependent on the Reynolds particle number Re p . The impact of using an extension in the form of a voidage term was found to be significant. The variations of the parameter χ influence the migration of the particles and the fluid pressure field as the change of the void fraction have an increasing impact on the computation of the drag force with increasing χ. In conclusion, the following findings can be generalised: I. Without the voidage term, no significant suffusive behaviour could be observed. By extending the Schiller and Naumann drag force with a voidage term, the migration of the fine particles increased significantly. Thereby, it turned out that the increase of the parameter χ leads to a greater distance covered by the particles. The amount of the particles, which are washed out, is dependent on χ. It increases accordingly. It was found that an increase in χ leads to more blow-outs of fine particles during the seepage flow. II. The pressure drop in the particle bed increases with greater χ. At the same time, also the pressure itself in the particle bed increases, while the velocity at the inlet does not alter noticeably using the different drag force correlations. This pressure increase is comparable with the results by Li and Kuipers, 7 who experienced a similar increase in fluidised bed reactors. III. The fluid-particle interaction term increases significantly with increasing values for χ. This leads to a higher amount of kinetic energy loss in the fluid and a subsequent reduction of the fluid velocity. At the same time, it leads to stronger particle migration. IV. The drag force increases according to the increase of the drag force models by the extension with the voidage function ε −χ . The stronger particle movement with increasing value for χ results in fluctuations of the drag force acting on the particles. Additionally, the increase of drag force at the beginning of the simulation increases with increasing drag force.
V. The increase of the parameter χ and, accordingly, the increase of pressure in the bottom part of the soil column leads to a loss of stability in this area. Concurrently, it occurs with a gain of stability in the top part of the column due to the number of fine particles migrating into that layer. The analysis of the contact force network revealed a gain of stability in the top part of the column with increasing the parameter χ.
This summary shows that the choice of the drag force model used in numerical simulations of suffusion with the unresolved coupled CFD-DEM has a crucial impact on the results. The migration behaviour of the particles changes significantly. Hence, already slight differences in the drag force correlations cause different results with the identical numerical setup.
This contribution furthermore reveals the need for a thoughtful validation of the numerical simulations carried out with unresolved CFD-DEM. To verify if the extended drag force models depict the general, correct physical behaviour of particles in a fluid, a consolidation test was carried out. The results of the particle settlement of the different drag force models were similar. However, the different drag force models showed very different behaviour in the numerical simulation of suffusion. This emphasises the fact that the validation of a numerical model cannot be completed with standard validation tests alone (such as the sedimentation of a particle) as the responsible physical actions vary. These standard validation tests are often carried out without further comparison with the actual problem. The results of this contribution indicate that this is problematic. The analysis of the drag force model of Di Felice shows an almost constant value for the exponent χ in the voidage function in the investigated regime of Re p . The presented investigations on this parameter reveal its significant influence on the suffusion processes in a particle column. In future work, the drag force model of Di Felice needs to be verified carefully against suffusion problems. It is necessary to verify if calibrations of the exponent χ are required. This contribution reveals that χ can be adapted to reflect experimental results. By comparing the particle migration or measuring the flow field, the parameter χ can be calibrated so that the particle behaviour in the numerical model depicts the reality.