**2. Methodology**

#### *2.1. Computational Domain and Meshing*

In order to verify the accuracy of the numerical method compared with the experimental data [14], a single-channel pump consisting of the inlet pipe, impeller, and volute was selected in this work to study the inner flow field. This single-channel pump model was kept the same as that used in experiments, and the basic parameters of the centrifugal pump are as follows: design discharge Q = 210 m<sup>3</sup>/h, head H = 20 m, rotational speed n = 1800 rpm, inlet diameter D1 = 50 mm, impeller diameter D2 = 138 mm, pump outlet diameter D3 = 110 mm. The polyhedron meshes are generated in the entire computational domain, as shown in Figure 1a. In addition, for the turbulent flow simulation, an appropriate resolution of the near-wall region is needed, and 5 prism layers are created next to all the wall surfaces (see Figure 1b) to improve the accuracy of the flow solution. Table 1 shows the results of the mesh dependency test for the head of the pump. It shows that the pump head remains steady when the grid number exceeds 535,665. Given the grea<sup>t</sup> amount of CFD-DEM coupling calculations, relatively few meshes are significant for simulation efficiency in the subsequent optimization process. Hence, the optimal number of mesh cells was determined as 535,665. As a kind of real-mass particle, the acceleration of gravity (*g* = 9.81 m/s2) is also taken into consideration, with its direction opposite to the *y*-axis.

 **Figure 1.** (**a**) Computational domain; (**b**) mesh cells.


## *2.2. Governing Equations*

The equations of continuity and momentum (Navier–Stokes equation) that govern the liquid phase are as follows in the tensor form:

$$\frac{\partial}{\partial t}(\alpha\_f \rho\_f) + \frac{\partial}{\partial x\_j}(\alpha\_f \rho\_f u\_j) = 0,\tag{1}$$

$$\frac{\partial}{\partial t}(a\_f \rho\_f u\_i) + \frac{\partial}{\partial \mathbf{x}\_j}(a\_f \rho\_f u\_i u\_j) = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \left[ a\_f \mu\_{eff} \left( \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right) \right] + a\_f \rho\_f \mathbf{g} + \mathbf{F}\_{\mathbf{s}\_i} \tag{2}$$

where ρ*f* is the fluid density, *u* is the fluid velocity, *p* is the pressure of the fluid, μ*e*ff is the e ffective viscosity, *x* is the coordinates, **g** is the acceleration of gravity, and **F***s* is the drag force between the particles and the liquid. <sup>α</sup>*f* represents the porosity around the particle, which can be calculated as:

$$\alpha\_f = 1 - \sum\_{i=1}^{n} V\_{p,i} / V\_{\text{cell}\prime} \tag{3}$$

where *Vp,i* represents the volume of particle *i* in a CFD cell, *n* represents the number of particles inside the cell, *Vcell* represents the volume of the cell.

Particle trajectory and particle behavior are controlled by the combined influence of hydrodynamic interactions and external force, which are described by Newton's laws of motion. According to DEM, the translational and rotational movements of a particle can be calculated with Newton's kinetic equation:

$$m\frac{d\mathbf{v}}{dt} = m\mathbf{g} + \sum \mathbf{F}\_{\mathbf{c}} + \mathbf{F}\_{\text{drag}} \tag{4}$$

$$\mathbf{I}\frac{d\mathbf{w}}{dt} = \sum \mathbf{T}\_c + \mathbf{T}\_f \tag{5}$$

where **F***c* represents the contact force, **<sup>F</sup>***drag* represents fluid drag force, m is the particle mass, and **I** is the moment of inertia of the particle. *d***v**/*dt* is the translational acceleration of the particle, *d*ω/*dt* is the angular acceleration of the particle, **T***c* is the contact torque, and **T***f* is the torque caused by the fluid.

#### *2.3. CFD-DEM Coupling*

During the two-way coupling adopted in this study, the modeling of the particle by DEM code is at the individual particle level, while the liquid flow by the CFD solver is at the computational cell level [15]. The coupling is initiated by analyzing the liquid flow based on CFD simulation. Once the iterative calculation in the CFD simulation is converged within a given time step, the DEM simulation then starts to calculate the instantaneous fluid–particle interaction force exerted on the individual particles from the CFD results. Thereafter, the new position and velocities of all the particles in the next fluid time step are determined. Inputting the updated particle information into the CFD solver, the forces on the fluid from particles are subsequently introduced into the liquid for the next loop through a series of momentum source terms [16].

#### *2.4. Particle Model*

In order to factor the size and shapes of particles, the simulation was divided into two groups. One group (mixed-size particles group) represents spherical particles with a certain range of diameters, and another (di fferent-shaped particles group) represents three types of particles with di fferent shapes. For the mixed-size particles group, the size of the particles varied from 1.0 to 5.0 mm in diameter in a random manner. For the di fferent-shaped particles group, the shapes of particles in the simulation were set as a cylinder, pie, and sphere (see Figure 2a–c, respectively), which are the most three representative shapes of foreign matter in sewage. In addition, the cylindrical and pie-shaped particles were formed from rigid sphere clusters, where the adhesion force between spheres was set to infinity, and they would not separate during the simulation. The density of the particles was taken as 2000 kg/m3, and the particles were selected in a diameter of 3.0 mm, while the particle flow rate was set to 4000 particles/s at the pump inlet. When the particles contact each other or solid boundaries, their momentum and energy are exchanged. This necessitates a DEM phase interaction model to describe the particle–particle and particle–wall interactions. In this model, the Hertz–Mindlin contact model [17] was adopted to solve the contact forces of particles, which were described by a soft-sphere model [18]. The Hertz–Mindlin contact model is the standard model that was used for describing the particle–particle and particle–wall interactions. The collision parameters used in the models are summarized in Table 2.

**Figure 2.** Particle shape models: (**a**) cylinder particle, (**b**) pie-shaped particle, (**c**) spherical particle.


**Table 2.** Collision parameters in the Hertz–Mindlin contact model.

#### *2.5. Fluid Phase Setup*

In this work, the transient fluid phase flow was analyzed through CFD simulations, based on solutions of the transient Reynolds averaged Navier–Stokes (RANS) equation. The realized two-layer k-ε model and "high-y+ wall treatment" in the STAR-CCM+ platform were employed for turbulence modeling. This wall treatment assumes that the near-wall cell lies within the region of the boundary layer. In the present simulations, the values of y<sup>+</sup> ranged from 0.31 to 259.74, and the average y<sup>+</sup> values for the impeller and volute surfaces are 80.38 and 109.64, respectively. The distribution of y<sup>+</sup> for the impeller and volute is shown in Figure 3. The flow was assumed to be isothermal and incompressible with the properties of water (μ = 8.887 × 10−<sup>3</sup> Pa·s, ρ = 998 kg/m3). The pump walls as well as particle surfaces were defined as no-slip walls, and the pump outlet was defined as a pressure outlet with *p* = 1.0bar. At the velocity inlet, a constant profile was specified. The absolute convergence criterion for the calculated residuals was set as 10−<sup>4</sup> by default, together with monitoring the variations in average pressure at the pump inlet during the computations, as shown in Figure 4.

**Figure 3.** Distribution of y<sup>+</sup> for the impeller and volute.

**Figure 4.** Variations in the average pressure at the pump inlet.

#### **3. Results and Discussion**

## *3.1. Validation*

For verifying the reliability of the CFD-DEM coupling method, the numerical results of the total head were compared with the existing experimental data [14] for the single-channel pump, as shown in Figure 5. The flow rate, size, and rotational speed of the computational pump model are the same as those of the pump model 22 used in the experiment. The total computational head of the pump is 18.41m, showing a good agreemen<sup>t</sup> with that of the experimental one—19.85 m. Additionally, this indicates the feasibility of the CFD-DEM coupling method. Meanwhile, the head predicted by means of CFD-DEM was a little bit lower than the experimental one. This discrepancy could result from the clean water used in the experiment compared to the solid-liquid two-phase flow in the simulation.

**Figure 5.** Comparison of the total head between the present simulation and the experimental data [12].

#### *3.2. Mixed-Sized Particle*

#### 3.2.1. Particle Trajectory and Distribution

Figure 6 shows the trajectory and size distribution of spherical particles whose diameter varied from 1.0 to 5.0 mm in a random manner (see Figure 6a). Remarkably, to reduce the impeller weight and unbalanced mechanical force, the blade is normally hollow (see Figure 6b). In this case, the hollow area is large enough that the flow status and particle motion in it cannot be ignored. In the impeller passage, the particle number density is different in different regions, with the characteristic of larger density in the vicinity of the inlet section and the blade pressure side.

(**a**)

(**b**)

**Figure 6.** Particle size distribution (**a**) and trajectory (**b**) mixed-size particles.

Besides this, to distinguish the trajectories of different-sized particles, the larger particles are presented in warmer colors, while smaller particles are shown in cooler colors. In general, the particles tended to maintain a steady trajectory towards the volute that corresponds with the shape of the impeller blade. After entering the volute, the particles tended to cluster along the volute outside wall and move downstream towards the outlet. However, the smaller particles had a much more uniform distribution in the passage of the impeller and volute than the larger ones. The most likely cause of this was the lower inertia of the smaller particles.

#### 3.2.2. Particle Velocity Distribution

Figure 7 shows the scatter plot of the particle velocity and diameter. It can be seen that the smaller size the particles possess, the greater the velocity distribution range and velocity peak they get. Generally, the majority of the particle velocity would be between 1 and 7.5 m/s. In this case, when the particle size is tiny enough, the maximum speed reaches 14.7 m/s. With the increase in particle size, the velocity distribution range becomes smaller. This can be manifested by the decrease in the maximum speed—that is, the smaller particles are more likely to acquire a higher velocity.

**Figure 7.** Scatter plot of particle velocity and diameter.

Figure 8 presents the variations in the volume-average velocities of both the liquid phase and solid phase in the pump domains. The results show that the average velocity of the liquid phase approached a nearly steady value of 4.7 m/s, while that of the solid phase reached a considerably lower value of 4.1 m/s. There are apparently slip velocities between the liquid and solid phase flows inside the pump. Furthermore, the average velocity of the liquid phase approached a steady value after t = 0.15 s, while the average velocity of particles spent 0.24 s attaining a stable value. Preliminary analysis suggests that the dominant influence on particles could be exerted by the liquid phase. Coupled with the inertia of the particles, the average velocity of the particles would take longer to reach a plateau.

**Figure 8.** Variations in the volume-average velocities of the liquid and particles in flow domains.

#### 3.2.3. Particle Contact Force Distribution

Figure 9 shows the histogram of the particle size to contact force at t = 0.5 s. It is obvious that the larger the particle size is, the greater the contact force particles would be. The small particles (diameter less than 2.4 mm) had a fairly weak but even contact force compared to the large particles. Although the contact force of the large particles increased significantly due to the greater surface area and mass, the force was rather uneven, especially for those particles whose diameters ranged from 2.4 to 4.0 mm.

**Figure 9.** Histogram of the particle size to contact force.

#### *3.3. Di*ff*erent-Shaped Particles*

#### 3.3.1. Particle Trajectory and Distribution

Figure 10 shows the trajectories and distribution of three distinct-shaped particles at the same particle flow rate. In the comparison of these three figures, it can be seen that the trajectory of cylindrical particles in the impeller was close to the pressure side (see Figure 10a). After acquiring a significant speed from the impeller, cylindrical particles tended to scatter in the volute passage. With regard to the trajectory of the pie-shaped particles, it was similar compared to that of the cylindrical particles but closer to the pressure side of the impeller passage (see Figure 10b). Meanwhile, the pie-shaped particles were distributed more uniformly than the cylindrical particles in the volute passage. In contrast, the distribution of spherical particles became most uniform in the impeller and volute compared to the other two cases (see Figure 10c). This may indicate that the motion of spherical particles would be mainly under the influence of the liquid phase so that they can disperse uniformly in the passage.

**Figure 10.** Trajectory and distribution of three distinct-shaped particles (**a**) cylinder practice; (**b**) pie-shaped particle; (**c**) spherical particle.

#### 3.3.2. Velocity Field of Liquid Phase

Figure 11 shows the relative velocity distribution of the liquid phase at the mid-span of the impeller for the three shape cases. It can be seen that there is little difference in the relative velocity fields of the three cases, which may reveal the subtle effect on liquid exerted by particles with different shapes under this particle flow rate. Additionally, it can be clearly seen that a wide low-velocity region is generated in the vicinity of the pump inlet. This can be explained by the fact that the inlet angle of the blade is pretty small, causing stagnation points both on the pressure side and suction side near the inlet edge. Thus, a certain range of flow separation occurs along the downstream direction of the inlet edge. Simultaneously, there is a vortex field in the hollow area of the blade, causing a backflow at the interface of the impeller volute. These second flows would significantly impact the pump performance.

**Figure 11.** Relative velocity field of the liquid phase at the mid-span of the impeller for the three shape cases: (**a**) cylinder particle, (**b**) pie-shaped particle, (**c**) spherical particle.

On the other hand, compared with the results of the particle distribution, it can be reasoned that the liquid phase might possess the dominant influence on particles. Especially for those particles in the hollow area of the blade, the backflow would carry the particles back to the hollow area.

#### 3.3.3. Collision Between Particle and Wall

The following figure illustrates the particle–wall contact force for three shape cases (refer to Figure 12). It is clear that the pie-shaped particles have the most serious collision, and spherical particles have the least in total. When it comes to the collision of particles with hub and shroud, it can be obtained that both the hub and shroud suffer a minor contact force, whose value was under around 1.0 N, by all the three types of particles. On the contrary, the blade and volute wall both sustained considerable contact force from particles, especially from the cylindrical and pie-shaped particles. This indicates that these solid particles may cause severe wear problems on the blade and volute areas. Hence, in the further optimization of the single-channel pump, the blades and volute can be designed to have more wear-resistance than the other parts. Moreover, before entering the pump, the particles or foreign bodies could be made as round as possible to reduce wear.

**Figure 12.** Particle–wall contact force for the three shape cases.
