The flow of fluids adheres to fundamental physical conservation laws, primarily including the conservation of mass, energy, and momentum. The governing equations are the mathematical descriptions of these conservation laws.

The mass conservation equation, referred to as continuity equation (Mahmut, 2023), is given as Equation 1:

Where, ρ is the fluid density; u is the velocity vector.

The law of momentum conservation, which is based on Newton’s second law, is governed by the following conservation Equations 2–4 in the three Cartesian directions (x, y, z):

Where, p represents pressure, μ denotes dynamic viscosity; Su, Sv, Sw are the generalized source terms of the momentum conservation equations.

The set of these three governing conservation equations collectively forms the Navier-Stokes equations.

The Navier-Stokes equations for incompressible flows are given as Equation 5:

The Navier-Stokes equations for inviscid flows are given as Equation 6:

This study employs the Reynolds-Averaged Navier-Stokes (RANS) approach. The standard k−ϵ model is selected as the turbulence closure, which characterizes turbulent motion via two transport equations. One equation governs the turbulent kinetic energy (k), while the other governs its dissipation rate (Tamburini et al., 2011).

The transport equation for the turbulent kinetic energy (k) is given as Equation 7:

The transport equation for the turbulent dissipation rate (ϵ) is given as Equation 8:

For addressing multiphase flow problems, the Volume of Fluid (VOF) model is employed in this work. It is commonly used to simulate multiphase flows with distinct interfaces, and its governing equations are given as Equation 9:

Where, ρq; uq represent the density and velocity of phase q, respectively; Sαq is the mass source term; mpq denotes the mass transfer from phase p to phase q, and mqp signifies the mass transfer from phase q to phase p.

In the present study, the fluid system comprises a primary phase and a secondary phase. Their respective volume fractions within each control volume are denoted by α1 and α2, governed by the Equation 10:

Where, α1 represents the volume fraction of the primary phase; α2 represents the volume fraction of the secondary phase.

Within each control volume, the effective fluid density and viscosity are determined by the sum of the products of each fluid phase’s physical properties and its respective volume fraction (Mark et al., 2015). The expressions are given as Equations 11, 12:

The momentum conservation equation is expressed as Equation 13:

Where, u represents the velocity vector; ρ denotes the fluid density; p is the pressure.

In the continuity equation, the density and volume fraction of each distinct phase in the gas-liquid two-phase flow are obtained by solving for the phase fractions. However, in the momentum equation, the two fluid phases share a single, unified equation. Consequently, when using the VOF model, only the velocity of the mixture phase is resolved, and the individual phase velocities cannot be directly obtained.

In this study, an Eulerian-Lagrangian framework is adopted within the coupled CFD-DEM approach, where the Discrete Element Method (DEM) numerically solves Newton’s second law of motion to compute individual particle motion, with particle-particle interactions calculated via contact models (force analysis illustrated in Figure 2).

Particles are subjected to various forces (fluid-induced forces, gravitational force, and interparticle contact forces) as well as torques arising from tangential forces and rolling friction, leading to the governing equation for particle motion are given as Equations 14, 15:

Where, mi denotes the mass of the i−th particle; vi represents the velocity vector of the i−th particle; Fc,ij signifies the contact force between particles i and j; Flr,ik indicates the non-contact force between particles i and j; Fpf,i is the fluid-induced force acting on particle i. Furthermore, wi is the angular velocity of particle i; Ii is its moment of inertia, while Mt,ij and Mr,ij represent the tangential friction torque and rolling friction torque, respectively, between particles i and j.

This investigation focuses on a dual-impeller mixing system, whose tank assembly comprises a vessel and two impellers—each fitted with two 45-degree pitched blades. The particle inlet and water inlet were set as velocity-inlet boundaries, while the outlet was defined as a pressure-outlet boundary. Its structural configuration is illustrated in Figure 3, with corresponding geometric parameters detailed in Table 1.

The simulations in this study were conducted with ANSYS FLUENT 2023 R1. The sliding mesh technique was used to simulate the rotational motion of the impeller. The computational mesh is shown in Figure 4, with other specified parameters listed in Table 2.

This study employs a one-factor-at-a-time strategy, adhering to the principle of controlled variables, to systematically investigate the influence mechanisms of key parameters on mixing characteristics. The time step was set to 0.001 s. The parameter settings for all investigated cases are detailed in Table 3.

This study adopts the homogeneity evaluation method for particle distribution proposed by Zhang (2020), defining the Relative Standard Deviation (RSD) as Equation 16. A lower RSD value indicates a superior degree of particle mixing homogeneity.

Where, M is the number of sampling cells; Ni is the number of particles within the i−th sample cell; Nav is the average number of particles per sample cell.

Mesh size and quantity are critical factors influencing computational accuracy and precision. Therefore, a grid independence study is essential prior to simulation to ensure the results are independent of the mesh resolution. Simulations were conducted using three distinct mesh configurations with approximately 0.3 million, 0.4 million, and 0.5 million elements, respectively, at an impeller speed of 150 rpm. The velocity profiles along the central axis of the tank were compared for these meshes, as shown in Figure 5.

As shown in Figure 5, the three mesh configurations exhibit minimal influence on the computed velocity field. Nevertheless, to maintain a balance between computational accuracy and cost efficiency, the mesh with approximately 400,000 elements was selected for all subsequent simulations.

The impact of particle density on mixing effectiveness was investigated with the agitation speed, particle size, and operational duration fixed at 100 r/min, 5 mm, and 20 s, respectively. The corresponding results are presented in Figure 6.

Figure 6A reveals similar velocity variation trends across the three density levels within the axial plane. Axial velocity remains around 2 m/s between x = -0.15 m and 0.15 m, while tangential velocity exhibits a pronounced negative peak at x = -0.15 m and a symmetric positive peak at x = 0.15 m, indicating strong rotational flow near the impellers. Radial velocity, in contrast, changes more smoothly along the plane. As shown in Figure 6B, velocity magnitudes fluctuate most noticeably near x = ± 0.15 m, coinciding directly with the impeller shafts. The distribution at ρ = 1000 kg/m³ shows consistently more pronounced variations than those at higher densities, with all velocity components gradually stabilizing in the upper tank region as the flow becomes more uniform.

For turbulent kinetic energy (TKE), two heights levels (y = 0.5 m and 0.9 m) were compared. In the core agitation zone (y = 0.5 m, Figure 7A), TKE displays a negative correlation with particle density: the peak TKE is 0.007 for 1000 kg/m³ particles, drops to 0.004 for 1100 kg/m³ (a 42.9% reduction), and declines for 1200 kg/m³. This attenuation results from the greater inertia of denser particles, which consume more turbulent energy for suspension and exert stronger damping on fluid fluctuations. In the weaker agitation region (y = 0.9 m, Figure 7B), the TKE baseline is lower and the influence of density is less pronounced.

Mixing homogeneity, evaluated through the Relative Standard Deviation (RSD), further confirms the clear negative correlation with particle density (Figure 8). This trend is evident across the tested range: At 1000 kg/m³, RSD reaches its minimum (2.7742), as the near-neutral buoyancy allows particles to closely follow turbulent eddies and thereby achieve a highly uniform dispersion. At 1100 kg/m³, RSD rises notably to 5.7708, marking a 108% increase, which indicates a developing tendency for local particle accumulation. At the highest tested density of 1200 kg/m³, RSD reaches 7.8463—a further 36% rise and 183% above the optimum. Here, particle gravity markedly exceeds fluid drag, leading to pronounced settling and significant aggregation in low-flow zones, which severely compromises mixing uniformity.

Thus, lower particle density promotes turbulent suspension and enhances mixing uniformity, as particles with density close to that of the fluid can better follow turbulent eddies and maintain stable dispersion. Conversely, higher densities impair homogeneity due to increased inertia and gravitational settling, which lead to particle accumulation, especially in low-flow zones of the tank.

The influence of particle diameter on mixing performance was investigated with other parameters held constant (agitation speed: 100 r/min, particle density: 1000 kg/m³, operational duration: 20 s), and the results are shown in Figure 9.

As observed in Figure 9, velocity profiles exhibit similar trends for all particle sizes in both radial and axial distributions. In the axial plane (Figure 9A), radial and tangential components show more pronounced variations, with velocity peaks occurring at Y = ± 0.15 m. Smaller particles (3 mm) display marginally higher peak velocities, though differences are limited. On the radial plane (Figure 9B), velocity is highest near the impeller shafts and decreases toward the tank wall. Again, 3 mm particles maintain slightly elevated velocities compared to their larger counterparts, reinforcing the observation of enhanced momentum and energy transfer to the upper tank region.

Turbulent kinetic energy (TKE) analysis further clarifies the role of particle size. In the core agitation zone (y = 0.5 m, Figure 10A), a clear negative correlation between particle size and TKE is observed: 3 mm particles achieve a peak TKE of 0.012, which is about 42% higher than that for 4 mm and 5 mm particles (TKE ≈ 0.007). This enhancement stems from the larger number of smaller particles (at constant total mass), which increases inter-particle collision frequency and promotes turbulent fluctuations, thereby enhancing the overall mixing dynamics in the system. This higher turbulence level directly improves the suspension and distribution of solid materials. In the weaker agitation region (y = 0.9 m, Figure 10B), the influence of particle size diminishes, and TKE distributions become more uniform across sizes.

Mixing homogeneity, quantified by the Relative Standard Deviation (RSD), confirms the advantage of smaller particles (Figure 11). The RSD decreases from 2.7742 at 5 mm to 2.6931 at 4 mm (a 2.9% reduction), and reaches its minimum of 2.5795 at 3 mm (a further 4.2% decrease, totaling 7.0% improvement over 5 mm). At the smallest size, particles are fully dispersed throughout the flow field without noticeable concentration zones, indicating an optimal state of suspension. This uniform distribution is crucial for ensuring effective contact in subsequent application stages. This superior mixing performance arises because smaller particles experience lower gravitational settling and offer a larger total interfacial area for fluid-particle interaction, thereby achieving more uniform suspension.

To conclude, reducing particle size leads to a consistent improvement in mixing homogeneity. This enhancement is primarily attributed to the increased total particle surface area and number count at a given mass, which effectively intensifies inter-particle collisions and significantly strengthens coupling with the turbulent flow field. Consequently, smaller particles are more readily suspended and distributed throughout the vessel, minimizing dead zones and localized accumulation. Under the conditions studied, the optimum performance was achieved with a particle size of 3 mm, which yielded the most uniform dispersion and the lowest RSD value.

The influence of agitation speed on mixing performance was investigated with particle density, particle size, and operational duration fixed at 1000 kg/m³, 5 mm, and 20 s, respectively (Figure 12).

Velocity distributions on the axial plane (Figure 12A) show consistent variation trends across the three agitation speeds, but with distinct differences in magnitude. Within the region from -0.15 to 0.15 m, axial and radial velocity peaks increase substantially with speed: at 100 r/min they reach 2 m/s and 1.5 m/s, respectively; at 150 r/min they rise to about 3 m/s and 2 m/s; and at 200 r/min they attain 4 m/s and 2.5 m/s. The steeper velocity gradients at higher speeds reflect stronger impeller-induced shear and more concentrated kinetic energy. On the radial plane (Figure 12B), the axial velocity maximum similarly doubles from 2 m/s at 100 r/min to 4 m/s at 200 r/min, and the recirculation zone expands by about 60%, indicating enhanced overall flow activity.

As shown in Figure 13a consistent trend of significantly enhanced turbulent kinetic energy (TKE) with increasing agitation speed is observed in both the core agitation region at y = 0.5 m and the weak agitation region farther from the impeller at y = 0.9 m. The TKE elevation is more pronounced with larger increments in rotational speed. In the core agitation region (y = 0.5 m), the influence of agitation speed on TKE is more dominant. At 200 r/min, the TKE profile lies entirely above those for the 150 and 100 r/min cases across the entire x-axis range. The difference is most prominent at the peak, where the highest speed achieves a 62.5% increase compared to the lowest speed. While the overall TKE baseline in the weak agitation region (y = 0.9 m) is lower than in the core region, the positive correlation between speed and TKE persists. However, the TKE differences between the speed cases are smaller than those in the core region, and the radial distribution is more uniform. This reflects the attenuation of the speed’s influence due to energy dissipation during transport.

Mixing homogeneity improves markedly with agitation speed (Figure 14). The RSD decreases from 2.7742 at 100 r/min to 2.2070 at 150 r/min, and reaches its minimum of 1.6231 at 200 r/min. At the lowest speed, particles tend to accumulate in the lower tank and flow-dead zones, whereas at 200 r/min they become uniformly dispersed without noticeable concentration. This improvement is attributed to the enhanced fluid drag force at higher speeds, which more effectively counteracts particle gravity and inertia, thereby preventing sedimentation and promoting full suspension.

Overall, higher agitation speeds significantly improve mixing uniformity by strengthening turbulent agitation and fluid-particle interaction, with 200 r/min providing the optimum performance under the studied conditions.

The influence of operational duration on mixing performance was investigated with particle density, particle size, and agitation speed kept constant at 1000 kg/m³, 5 mm, and 100 r/min, respectively. The corresponding results are presented in Figure 15.

As shown in Figure 15A, radial velocity distributions on the axial plane exhibit the most intense fluctuations at 25 s: axial velocity forms a distinct high-velocity zone between -0.15 m and 0.15 m, dropping to near zero at the edges (peak-to-peak difference ≈ 3m/s), while tangential velocity extremes approach ±2 m/s. By 30 s, all velocity components become smoother and more continuous; the axial velocity peak-to-peak magnitude shrinks to 2 m/s, tangential fluctuations decrease by over 30%, and radial fluctuations nearly vanish, indicating that the flow field has reached a dynamic equilibrium. This equilibrium state signifies a stable and energy-efficient mixing regime.

On the radial plane (Figure 15B), the axial velocity at 25 s shows two fluctuations near the impeller shafts, with a maximum of 3 m/s. At 30 s, the peak axial velocity decreases by 33%, transitions become smoother, and both radial and tangential velocities stabilize near zero, reflecting uniform axial particle distribution. Such uniformity is essential for achieving consistent homogenization results.

Turbulent kinetic energy (TKE) decreases with increasing operational time in both the core agitation region (y = 0.5 m) and the weaker region (y = 0.9 m), stabilizing notably after 25 s (Figure 16). In the core region, peak TKE drops from 0.0065 at 20 s to 0.0045 at 25 s, and the curve nearly overlaps completely with that at 30 s. Similarly, in the weak agitation region, peak TKE decreases from 0.006 to 0.003 at 25 s and remains stable thereafter. These consistent trends confirm that the turbulent fluctuations attenuate and reach a steady state after 25 s.

Mixing homogeneity improves significantly with time up to 25 s, then stabilizes (Figure 17). At 20 s, the RSD is 2.7742, indicating uneven particle distribution with accumulation in lower tank and dead zones.

At 25 s, the RSD reaches its minimum (2.1089), corresponding to uniform particle dispersion without concentration zones and a 31.2% reduction in grid-particle-count deviation. At 30 s, the RSD shows a slight increase to 2.1581 (a 2.3% rise relative to 25 s), attributed to weak particle-particle collisions causing minor re-agglomeration, but the overall distribution remains uniform.

Therefore, mixing homogeneity improves with operational time until a dynamic equilibrium is reached at 25 s, at which point the particles achieve a fully uniform suspended state within the flow field. Further extension of mixing time yields no substantial improvement and may even induce slight re-agglomeration, which is likely caused by weak inter-particle collisions, indicating that 25 s represents the optimal operational duration under the given conditions.