# Investigation of Movement and Deposition Behaviors of Solid Particles in Hydraulic Water Reservoir via the CFD-DEM Coupling Method – Chinese Journal of Mechanical Engineering

Sep 7, 2022

### Solving Process

To analyze the movement and deposition of solid particles in the reservoir, the CFD-DEM coupling method is implemented in this study. The flow field is simulated by ANSYS Fluent©, and the particle motion is simulated by EDEM©. The CFD-DEM coupling simulation schematic diagram is shown in Figure 3.

In the coupling process, fluid in the reservoir is regarded as the continuous phase, and particles are regarded as the discrete phase. Fluent simulates the transient flow field by solving the time-averaged Navier-Stokes equation with the SST kω model until the steady state. Then, the coupling calculation of particle motion and flow field is started. The forces calculated in the flow field on particles are transmitted to EDEM. It will update the speed and position information of particles and then transfer the reverse force from particles to the flow field. This process will be repeated on the next time step until the final time step is completed.

### CFD Model

In this case, the fluid is assumed to be isothermal and incompressible. The continuity and momentum equations for the fluid phase are given as follows:

$$frac{{partial rho_{{text{f}}} }}{partial t} + nabla cdot (rho_{{text{f}}} {varvec{u}}_{{text{f}}} ) = 0,$$

(1)

$$frac{{partial left( {rho_{{text{f}}} {varvec{u}}_{{text{f}}} } right)}}{partial t} + nabla cdot left( {rho_{{text{f}}} {varvec{u}}_{{text{f}}} {varvec{u}}_{{text{f}}} } right) = – nabla p + nabla cdot {varvec{tau}} + rho_{{text{f}}} {varvec{g}} + {varvec{F}}_{{text{s}}} ,$$

(2)

where ρf, uf, p, τ, and g indicate fluid density, fluid velocity, pressure, stress component, and gravitational acceleration. In addition, Fs represents the interaction term that involves the effect of the forces on particles.

### DEM Model

In the DEM model, the translational motion and rotational motion of a particle can be described according to Newton’s laws of motion [19, 20], which are:

$$m_{mathrm{p}} frac{text{d} {varvec{u}}_{mathrm{p}}}{text{d} t}={varvec{F}}_{mathrm{f}}+sum {varvec{F}}_{mathrm{c}} ,$$

(3)

$$I_{mathrm{p}} frac{text{d} {varvec{{omega}}}_{mathrm{p}}}{text{d} t}=sum {varvec{T}}_{mathrm{c}} ,$$

(4)

where mp, up, Ip, and ωp indicate particle mass, particle velocity, inertia tensor, and rotational velocity, respectively.

The fluid force on the particle Ff is given by:

$${varvec{F}}_{mathrm{f}}={varvec{F}}_{mathrm{G}}+{varvec{F}}_{mathrm{B}}+{varvec{F}}_{mathrm{P}}+{varvec{F}}_{text {Drag }}+{varvec{F}}_{mathrm{VR}}+{varvec{F}}_{mathrm{Saff}} ,$$

(5)

where FG, FB, Fp, FDrag, FVR, and FSaff indicate gravity, buoyancy force, pressure gradient force, drag force, virtual mass force, and Saffman force.

The contact force Fc includes the normal and tangential component, which is given by:

$${varvec{F}}_{mathrm{c}}={varvec{F}}_{mathrm{c}, mathrm{n}}+{varvec{F}}_{mathrm{c}, mathrm{t}} ,$$

(6)

where Fc,n is normal contact force, and Fc,t is a tangential contact force.

The contact torque Tc is generated by tangential contact force and rolling friction, which is given by:

$${varvec{T}}_{mathrm{c}}={varvec{T}}_{mathrm{t}}+{varvec{T}}_{mathrm{r}} ,$$

(7)

where Tt is tangential contact torque, and Tr is generated by rolling friction.

### Fluid-solid Interaction Force

There are several types of interaction forces between fluid and particles such as buoyancy force, pressure gradient force, drag force, virtual mass force, and Saffman force.

These interaction forces are given by:

$${varvec{F}}_{{text{p}}} = – frac{1}{6}uppi d_{{text{p}}}^{3} frac{{partial {varvec{p}}}}{partial x},$$

(8)

$${varvec{F}}_{text {Drag }}=frac{1}{8} C_{mathrm{D}} pi d_{mathrm{p}}^{2} rho_{mathrm{f}}left({varvec{u}}_{mathrm{f}}-{varvec{u}}_{mathrm{p}}right)^{2},$$

(9)

$$C_{{text{D}}} = frac{24}{{Re}} cdot left( {1 + 0.15 cdot Re^{0.687} } right),$$

(10)

$$Re,=,frac{{rho_{{text{p}}} d_{{text{p}}} left| {{varvec{u}}_{{text{f}}} – {varvec{u}}_{{text{p}}} } right|}}{mu },$$

(11)

$${varvec{F}}_{{{text{VR}}}},=,- frac{1}{2}V_{{text{p}}} rho_{{text{f}}} frac{{text{d}{varvec{u}}_{{text{p}}} left( t right)}}{{text{d}t}},$$

(12)

$${varvec{F}}_{{{text{Saff}}}} = 1.61left( {mu rho_{{text{f}}} } right)^{1/2} left( {frac{{d_{{text{p}}} }}{2}} right)^{2} left( {{varvec{u}}_{{text{f}}} – {varvec{u}}_{{text{p}}} } right)left| {frac{{{text {d}}{varvec{u}}_{{text{f}}} }}{{{text {d}}y}}} right|,$$

(13)

where CD, dp, Vp, ρp, μ, and Re indicate drag factor, particle diameter, particle volume, particle density, fluid kinetic viscosity, and Reynolds number.

### Particle Size and Shape Distribution Description

As fundamental physical parameters of particles, particle size and shape could affect the motion of particles and the interaction between particles and fluid. Therefore, it is essential to consider the actual shape and size of particles in the simulation [21].

Two hundred iron particles and sand particles samples in the experiment are randomly selected to represent the shape and size of the particles. Then, multiple SEM photographs are taken to observe the details of particles.

The typical SEM photographs for iron particles and sand particles are shown in Figure 4(a). The shape and size of iron particles are regular and uniform; their particle diameter is 500 μm. However, sand particles have various sizes and shapes. Hence, an average statistical method combined with projection is used to describe sand particles since the size and shape distribution of particles follow the normal distribution [22].

Figure 4(b) shows that the shape of sand particles is mainly cuboid. The aspect ratio can describe the two-dimensional shape of those particles. Its statistical mean value is 1.3 for quartz sand [22]. Due to the shape diversity of sand particles and the coarse size in Section 2.1., it is not facile to directly correspond with the DEM model. Therefore, photograph measurement results and projection area equality are used to represent the size of sand particles.

As the projection area equality is shown in Figure 5, the fitted rectangular area is equivalent to the actual projected area of sand particles. l and w indicate the best-fit length and width of sand particles. A matching circle is constructed to maintain the area equality with the fitted rectangular. The diameter of the circle is evaluated, which matches the normal distribution.

Among the SEM photographs shown in Figure 4(b), the largest and smallest sand particles in this image are taken as the up and down limit of the particle size to evaluate the matching circle equivalent diameter de. It is subject to the normal distribution, which is:

$$d_{{text{e}}} sim N(830,150^{2} ){upmu text{m}}.$$

(14)

Model particle diameter dm distribution in EDEM for sand particles could be conducted from the aspect ratio and circle diameter de, which is:

$$d_{{text{m}}} sim N(645,117^{2} ){upmu text{m}}.$$

(15)

The iron and sand particles models are consequently established in EDEM, shown in Figure 4(c).

### Stokes Number

Particle properties for following the fluid could be measured by Stokes number. It is the ratio of particle relaxation time and flow characteristic time [23]. When Stokes number is more significant, the time for particles to respond to flow is longer, which means the particles have a weaker ability to follow the fluid flow. The expression of Stokes number is:

$$S t=frac{tau_{mathrm{p}}}{tau_{mathrm{f}}}=frac{rho_{mathrm{p}} d_{mathrm{p}}^{2} {varvec{u}}_{infty}}{18 mu L},$$

(16)

where τp, τf, u, and L indicate particle relaxation time, flow characteristic time, flow characteristic velocity, and flow characteristic length.

### Simulation Setup

The structured CFD mesh of the reservoir fluid domain is shown in Figure 6. The liquid-particle mixing fluid flows into the inlet pipe and out of the outlet pipe to the water tank. Therefore, the inlet boundary condition is defined as a velocity-inlet condition of 1.24 m/s, and the outlet boundary condition is set as an outflow condition. The upper surface of the liquid in the reservoir is not restrained. As a result, the boundary condition for the upper surface of the model is set as symmetry. According to the experimental situation, the rest parts of the model are selected as the wall boundary condition.

Besides, it is essential to verify the grid independence of the reservoir flow field so that the grid number does not affect simulation accuracy [24]. By monitoring and recording the outlet flow data under different grid numbers, grid independence is guaranteed with the 4 mm grid size and the total of 1482233 grid numbers.

According to the experiment, the material of the reservoir model is polymethyl methacrylate. As a result, the wall material properties in EDEM include density (2500 kg∙m−3), Poisson’s ratio (0.25), and shear modulus (2.2×108 Pa). Because the particles flowing into the reservoir in experiments have a random direction and random velocity, the dynamic particle factory is defined on the inlet plane of the flow field. The particle’s release speed from the dynamic particle factory is set from 1 m/s to 1.5 m/s in random directions. Excepting parameters of particles and model, other the DEM setting parameters for iron and sand particles are listed in Table 1 [25,26,27,28].

In the CFD-DEM coupling simulation, the particle phase mesh does not affect the liquid phase mesh. EDEM mesh is used for particle tracking and retrieval, which affects the simulation speed. EDEM grid size is generally set to 3–5 times particle radius to maintain particle simulation accuracy [29]. In this work, the grid size is set to be three times particle radius to ensure simulation accuracy. Furthermore, the time step in Fluent should be greater than the EDEM time step. Typically, the time step in Fluent should be an integral multiple of time step in EDEM to capture the collision in particle motion [30]. In this work, the time step in Fluent is set to be 100 times the time step of EDEM. The total number of particles is 2000 until the particles are settled. The details of simulation parameter settings are shown in Table 2.