Quantification of Solid-Liquid Mixing Uniformity via Three-Dimensional Ripley’s L Function in DEM-VOF Simulations

1 Introduction

Stirred tanks are core pieces of equipment in the chemical, biochemical, pharmaceutical, polymer, petrochemical, and mineral processing industries, where solid-liquid mixing is crucial for diverse process objectives, including crystallization, polymerization, catalytic reactions, water treatment, dissolution, and leaching [1,2,3]. In these systems, the goal of mixing is not merely particle suspension, but the establishment of a homogeneous solid-liquid distribution to ensure stable mass transfer, heat transfer, and reaction environment [4,5]. However, in practical operations, particle distribution within stirred tanks is often characterized by multiscale non-uniformity, including large-scale stratification or circulation-induced particle accumulation, as well as small-scale clustering and local concentration fluctuations. Such non-uniformity can adversely affect industrial processes. At the macroscopic scale, particle accumulation or segregation may create dead zones, unstable suspension states, and uneven residence time distribution, thereby reducing the overall utilization efficiency of the reactor volume. At the mesoscopic and microscopic scales, localized particle-rich or particle-lean regions can result in heterogeneous mass transfer, uneven heat dissipation, and inconsistent reaction or dissolution environments. For crystallization and polymerization processes, this may cause broad product size distributions or variable product quality; for catalytic and leaching systems, it may lower the effective interphase contact area as well as reaction or extraction efficiency [6,7].

Mixing time is an important parameter for characterizing the hydrodynamic mixing state within stirred reactors and serves as a key indicator of mixing efficiency [8,9,10]. Experimental methods for its measurement are generally classified into two categories, namely intrusive and non-intrusive techniques [11]. Commonly adopted approaches for determining mixing time include colorimetric methods [12], conductivity measurements [13], optical techniques [14], and electrical resistance tomography (ERT) [15]. Conventional criteria for mixing time determination are primarily based on statistical indicators such as concentration standard deviation [16] or relative standard deviation [17]. These methods typically rely on signals obtained from a limited number of measurement points or averaged values over predefined spatial zones, thereby providing only single-scale quantitative information on overall mixing uniformity. Consequently, they are unable to capture detailed features of particle aggregation, stratification, and spatiotemporal evolution, and remain highly sensitive to sampling locations and spatial partitioning schemes. With the development of computational fluid dynamics (CFD) and related numerical simulation techniques, three-dimensional spatial coordinates of large particle populations during stirring can now be directly obtained. The effective exploitation of such spatial data for multiscale and quantitative characterization of solid-liquid mixing states has thus emerged as an important research issue in the field of solid-liquid mixing.

Numerical simulation is a powerful approach for investigating multiphase flow, mixing, and transport phenomena in stirred systems [18]. Existing numerical methods for particulate multiphase flows are generally classified into two categories: Eulerian-Eulerian (E-E) approaches [19] and Eulerian-Lagrangian (E-L) approaches [20]. In the E-E framework, the dispersed phase is treated as a continuum, offering relatively low computational cost but limited access to particle-scale dynamics. In contrast, the E-L framework tracks the motion of individual particles within the carrier fluid, making it more suitable for resolving microscopic behaviors such as particle-fluid interaction, particle-particle collision, dispersion, and spatial redistribution in stirred tanks [21]. For systems involving free surfaces or deformable phase interfaces, accurate interface capturing is also essential, since interfacial evolution strongly influences local flow structure, transport behavior, and particle distribution. The volume of fluid (VOF) method is widely used for this purpose because it effectively reconstructs and tracks transient fluid-fluid interfaces [22,23]. By coupling discrete element method (DEM) for particle motion with VOF-based fluid flow calculation, the DEM-VOF framework enables simultaneous resolution of particle dynamics and interface evolution, thereby providing a more physically representative description of complex multiphase systems. This coupled strategy has been applied for studying bubble-particle interaction, particle entrainment, and interfacial transport mechanisms in gas-liquid-solid or solid-liquid systems [24,25].

To achieve quantitative evaluation of complex multiphase mixing processes, mathematical methods are increasingly being used in the study of mixing behavior. Constructing statistical models and integrating them with algorithmic analyses can enable systematic characterization of mixing efficiency and distribution features. In recent years, such methods have been extensively applied to investigate flow, mass transfer, and dispersion characteristics in multiphase flow systems. The Quality Measurement Method (Q method) has been employed to assess the uniformity of two-dimensional point sets [26], but its applicability is limited by insufficient consideration of spatial correlations. Xu et al. [27] proposed a three-dimensional mixing metric combining the CD and WD indices, while Sun et al. [28] developed a mixing homogeneity evaluation model using the mean bubble distance based on geometrical statistical theory. Zhang et al. [29,30] introduced a PD method to investigate uniformity of two-dimensional and three-dimensional mixing systems. These studies have provided valuable guidance for the quantitative description of mixing processes.

However, existing studies are primarily focused on overall mixing uniformity, while making insufficient use of spatial point-pattern information of the dispersed phase. In particular, they remain limited in revealing the evolutionary features of particle distribution structures and the underlying spatial correlation characteristics. Moreover, in practical stirred systems, it is often difficult to obtain comprehensive three-dimensional spatiotemporal information on particle distributions solely through experimental methods. Owing to limitations in visualization conditions, spatial resolution, and sampling frequency, experiments usually provide only local or two-dimensional particle data, restricting their utility for quantitative analyses based on spatial point patterns.

To address these limitations, DEM-VOF numerical simulations were employed in this study to obtain the time-evolving three-dimensional spatial coordinates of particles in a stirred tank. On this basis, Ripley’s L function from spatial statistics was introduced for quantitative evaluation of solid-liquid mixing homogeneity [31]. Unlike existing evaluation methods that rely on regional averaging or single distance-based statistics, the proposed approach starts from the spatial point pattern of particles, analyses the structural evolution of particle distributions from clustering towards uniformity, and determines an appropriate characterization scale through numerical experiments. This enables quantitative evaluation of the mixing state [32]. The main objective of this study is to establish and validate a quantitative method for evaluating solid-liquid mixing homogeneity based on DEM-VOF and spatial statistical analysis, thus providing a theoretical basis and methodological guidance for the characterization and optimization of mixing processes in stirred tanks.

2 Mathematical Models and Methods

2.1 DEM-VOF Numerical Simulation Methodology

2.1.1 DEM Mathematical Model

In practical stirring processes, the spatiotemporal distribution of particles cannot be comprehensively obtained through experiments. Due to limitations in visualization techniques, spatial resolution, and sampling frequency, most experimental methods provide limited particle information within local regions or on two-dimensional cross-sections. This makes it difficult to construct three-dimensional spatiotemporal point pattern data suitable for Ripley’s L function analysis. Consequently, the direct application of Ripley’s L function to solid-liquid mixing systems is challenging. To overcome these limitations, a DEM-VOF numerical simulation approach is employed in this study to obtain the three-dimensional spatiotemporal distribution of particles within the stirred tank.

The DEM is a numerical technique specifically designed to simulate the dynamic behavior of discontinuous media [33]. To capture the complexity of multiphase flow involving particle swarms in agitated suspensions, a two-way coupled fluid-particle DEM model is adopted, and a comprehensive set of dynamic parameters is specified for the solid phase. This modeling framework enables systematic simulation of particle motion and interaction mechanisms within the flow field. Particle trajectories are governed by Newton’s laws of motion. For each particle, the governing equations for translational and rotational motion can be expressed as follows [34]: midvidt=∑jFc,ij+∑kFlr,ik+FPf.i+Fg,i(1) Iidωidt=∑(Mt,ij+Mr,ij)(2) here, m i , I i , v i , and d w i denote the mass, moment of inertia, linear velocity, and angular velocity of particle i, respectively. F c , i j , F l r , i k , F p f , i , and F g , i represent the contact force, non-contact force, particle-fluid interaction force, and gravitational force, respectively. The tangential and rolling friction torques are denoted by M t , i j and M r , i j , respectively.

A key feature of the DEM model is the calculation of both non-contact and contact forces. The contact force, which consists of normal and tangential components, is described as follows [35]: Fc,ij=fcn,ij+fct,ij(3) Fc,ij=−kn,ijδn,ij−γn,ijδ.n,ij−kt,ijδn,ij−γn,ijδ.t,ij(4) where k n , i j , k t , i j , γ n , i j , δ n , i j , and δ t , i j denote the normal stiffness coefficient, tangential stiffness coefficient, normal damping coefficient, tangential damping coefficient, and the normal and tangential particle overlaps, respectively. δ . n , i j and δ . t , i j represent their corresponding time derivatives. Further details are given in Table 1.

2.1.2 VOF Mathematical Model

To track the interface between the two phases, a volume fraction is assigned to each fluid phase. In the present study, only two phases, namely water and air, are considered in the VOF model. Water is treated as the primary phase and air as the secondary phase, and their volume fractions are denoted by α 1 and α 2 , respectively. Accordingly, each computational cell exists in one of three states: α i = 0 , indicating that phase i is absent from the cell; α i = 1 , indicating that phase i completely fills the cell; and 0 < α i < 1 , indicating that both phases coexist within the cell. The interfacial dynamics are obtained as the solution of the continuity equation for the volume fraction of the secondary phase [36]: ∂α2∂t+(u·∇)α2=0(5) here, α 1 and α 2 denote the volume fractions of the two phases, respectively; u is the velocity vector; and ∇ is the differential operator. And the calculation of the primary-phase fluid is based on the following conservation condition:

In a multiphase system, the density and viscosity of the fluid in each computational cell are determined by combining the physical properties and volume fractions of all continuous phases, as follows:

The mixing process in the reactor is essentially a rotating turbulent flow. Therefore, the Reynolds stress model (RSM) is employed for the turbulence calculation [37]. This model demonstrates good convergence and computational efficiency, and provides accurate predictions of turbulent flow fields at high Reynolds numbers, making it suitable for simulating the multiphase mixing process considered in this study. The governing equation is presented in Eq. (9): ρf[∂(εfu)∂t+∇·(εfuu)]=εf(−∇p+∇·(τ¯−ρfui′uj′¯)+fs+ρfg)+fpf(9) here, ε f denotes the porosity, and the fluid density and velocity are represented by ρ f and u , respectively, while ρ f u i ′ u j ′ ¯ represents the Reynolds stress. The terms on the right-hand side of Eq. (9) represent the pressure gradient, viscous stress, Reynolds stress, surface tension, gravity, and particle-fluid interaction force. The detailed equations are given in Table 2.

2.1.3 DEM-VOF Coupling Calculations

In the DEM-VOF model, only the dynamic interaction forces between the fluid and particles are considered, while the effects of virtual mass forces, Basset force, and lubrication forces are neglected. The governing equation is given as follows: Fpf,i=Fd,i+F∇p,i+F∇τ,i+Fsaff,i+FMag,i(10) where F d , i , F ∇ p , i , F s a f f , i , and F M a g , i denote the drag force, pressure gradient force, Saffman lift force, and Magnus lift force on particle i, respectively, and F ∇ τ , i represents the combined contribution of viscous stress and surface tension forces. The corresponding formulations and parameter definitions are summarized in Table 3.

2.1.4 Coupling and Data Exchange

In this study, the solid particle phase is simulated using EDEM, whereas the fluid phase is solved in ANSYS FLUENT. Fluid-particle coupling is realized through bidirectional data exchange between the two solvers. The coupling procedure can be summarized as follows. First, the stirred-tank geometry is constructed and meshed in FLUENT, and the VOF multiphase model, turbulence model, and other physical parameters are specified. The corresponding computational domain is then created in EDEM, where the particle size, density, contact model, and initial spatial distribution are initialized. After the computation is initiated, FLUENT solves the fluid governing equations at each fluid time step and passes the instantaneous cell-based velocity and pressure fields to EDEM. On this basis, EDEM evaluates the forces acting on the particles, integrates their equations of motion to update particle velocities and positions, and calculates the particle-fluid interaction terms, including drag, lift, and pressure-gradient forces. These interaction forces are subsequently written to a coupling file and transferred back to FLUENT, where they are introduced as source terms in the momentum equations for the subsequent time step to obtain an updated flow field. This FLUENT-EDEM-FLUENT loop, schematically illustrated in Fig. 1, is repeated throughout the simulation, enabling synchronous iterative updates of the fluid field and particle motion and thereby accurately capturing the bidirectional coupling between the complex flow structures and particle distribution in the stirred tank.

Figure 1: Coupling of computational fluid dynamics and discrete element method.

2.1.5 Validation of DEM-VOF Model

In the DEM-VOF framework, particles are treated as finite-volume solids. When particles are introduced into a water-filled vessel, the displacement of liquid leads to a rise in the free surface, making it necessary to verify volume conservation. According to the principle of volume conservation, the theoretical increase in liquid level should be equal to the ratio of the total particle volume to the cross-sectional area of the vessel bottom. To validate volume conservation and assess the ability of the DEM-VOF method to capture free-surface evolution, the entire process of a particle assembly freely settling under gravity from air into water was numerically simulated, from the initial entry into the liquid until all particles were fully submerged and reached a stationary state.

The computational domain consists of a rectangular container with dimensions of 0.05 × 0.05 × 0.20 m3, and the initial liquid level is set to 0.05 m. In the simulation, a total of 8000 particles with a diameter of 2 mm and a density of 2500 kg/m3 are adopted. During the settling process, the particles undergo free fall under gravity. As shown in Fig. 2, at t = 0 s the particle assembly is initially at rest above the free surface. Once settling commences, the particles successively penetrate the free surface and enter the liquid, resulting in a rise in liquid level due to the increasing volume of liquid displaced by the particles.

Figure 2: Snapshots of particles sinking into water.

The snapshots between 0 and 1 s show that, the free surface undergoes pronounced oscillations during the period in which all particles become fully submerged and migrate towards the bottom of the container. The liquid climbs along the vessel wall, while the free surface in the central region is drawn downward by the descending particle cluster. As a result, a parabolic free-surface profile is formed, characterized by a higher liquid level near the vessel wall and a depressed center. The maximum difference in liquid level occurs at t = 1 s. During this stage, only a small fraction of particles is entrained upward by the free-surface motion. Subsequently, as the particles continue to settle, the free surface gradually moves downward. As an increasing number of particles reach the bottom of the container and come to rest, the liquid-level difference diminishes, surface oscillations are attenuated, and the free surface eventually returns to a quiescent state. These findings are highly consistent with the simulation results reported by Sun et al. [39].

Theoretically, based on volume conservation, the rise in the free-surface level should be 13.40 mm. In the present simulation, the computed increase in free-surface height is 63.57 − 50.00 = 13.57 mm, corresponding to a relative error of only 1.26%, which lies within an acceptable range. These results demonstrate that the numerical model developed in this study exhibits good volume-conservation performance and accurately captures the dynamic evolution of system volume during solid-liquid mixing.

To further verify the reliability of the numerical simulation, a water-model experiment was performed using a flat-bottom glass vessel with a height of 240 mm, and the results were compared with the numerical simulation results. POM particles were used in the experiment, with a particle diameter of 2 mm, a density of 1200 kg/m3, and a total number of 20,000 particles. During the experiment, the mixing process was recorded by video, and images at representative time points were extracted to observe the particle distribution and its evolution characteristics. A geometrically scaled numerical model consistent with the experimental setup and operating conditions was established. The comparison between the experimental and simulation results is shown in Fig. 3.

Figure 3: Comparison and verification of entrainment experiment with deposited particles at the bottom of reactor: (a) water modeling experimental results; (b) simulation results.

As seen from Fig. 3, the experimental results are generally consistent with the simulation results in terms of particle distribution patterns and evolution trends. At t = 0 s, the particles are mainly deposited at the bottom of the vessel. When t = 0.2 s, under the entrainment effect generated by the rotating impeller, the bottom particles gradually gather towards the center and form a relatively obvious conical accumulation structure. At t = 0.4 s, more particles are entrained and lifted upward. A small portion of the particles reaches the blade tip region and begins to diffuse radially, while also exhibiting evident axial motion under the action of the main circulation flow. When t = 0.6 s, the particles are further lifted and dispersed into the vessel, and some particles enter the main circulation region. The direction of particle motion near the wall gradually aligns with the main flow direction of the fluid. As the POM particles used in the experiment were limited by machining accuracy, their actual shapes were not perfectly spherical and their diameters also had certain deviations. Therefore, some discrepancy between the experimental and simulation results is reasonable and acceptable. Considering that this experiment was mainly intended to characterize the particle distribution pattern and its evolution process, it can be concluded that the DEM-VOF coupled model reproduces the interaction between particles and fluid in the solid-liquid two-phase system with good accuracy, thereby providing a reliable basis for subsequent numerical simulation studies.

2.1.6 Simulation Setup

In this study, three axial-flow impellers were selected, as illustrated in Fig. 4: a conventional four-blade impeller, a four-blade turbine impeller, and a four-blade wide-blade hydrofoil impeller, which are denoted as impellers A, B, and C, respectively. To investigate the influence of blade installation angle on flow structure and solid-liquid mixing characteristics, each impeller was tested at blade angles of 45° and 90°. This design enabled systematic evaluation of the proposed method’s ability to discriminate between different blade angles and impeller types.

Figure 4: Impeller configurations and dimensions: (a) A-90°; (b) A-45°; (c) B-90°; (d) B-45°; (e) C-90°; (f) D-45°.

As illustrated in Fig. 5, the multiphase mixing system consists of a stirred tank equipped with an impeller. The tank has a height H of 240 mm and a bottom diameter D of 160 mm. The clearance between the impeller center and the tank bottom is denoted by C. The impeller has a blade width of 80 mm, a blade height of 18 mm, and a blade thickness of 2 mm.

Figure 5: Dimensions of the mixing system and mesh discretization: (a) configuration and detailed dimensions of the multiphase mixing system; (b) discretization of the fluid domain.

All six geometric configurations employ the same meshing strategy, dividing the fluid domain into a stationary fluid region and a rotating impeller region. The impeller rotational speed is set to 400 rpm, and the number of particles is fixed at 20,000. In the numerical solution, the pressure-velocity coupling is achieved using the SIMPLE algorithm. The pressure term is discretized using the PRESTO! scheme, the momentum equations are solved using a second-order upwind scheme, and the volume fraction is treated using the Geo-Reconstruct method. The turbulence-related terms are discretized using first-order upwind schemes, while the temporal discretization adopts a first-order implicit scheme. The DEM time step, CFD time step, and coupling time step are selected based on numerical stability considerations and with reference to the study by Wu et al. [40]. The physical properties of the fluid and particles, together with the detailed numerical parameter settings, are summarized in Table 4.

2.1.7 Grid Independence Verification

During domain discretization, significant momentum variations occur near the impeller due to its rotation. Appropriate mesh refinement in this region is therefore essential to ensure computational accuracy. To accurately resolve the flow field, an appropriate type and number of computational cells must be employed. In this study, the trajectories of solid particles are assumed to be largely governed by the local flow field. To ensure that particle motion can be adequately resolved within each computational cell, the minimum cell size is set to at least twice the particle diameter [40].

With respect to mesh sizing, the impeller thickness represents the smallest geometric length scale in the system; therefore, the minimum mesh size is adopted in this region to accurately capture the impeller geometry. In the stationary region, maximum cell sizes of 8 mm, 7 mm, 6 mm, and 5 mm are employed, corresponding to total mesh counts of 244k, 397k, 518k, and 733k, respectively. Mesh accuracy was assessed by comparing the axial velocity profiles extracted from the rotating impeller region, as shown in Fig. 6. Compared with the 733k mesh, the maximum axial velocity error along the selected line is 7.4% for the 244k mesh, 5.3% for the 397k mesh, and 4.7% for the 518k mesh. When the total mesh count reaches 518k, the computational error decreases to below 5%, satisfying the criterion for mesh independence.

Based on the above analysis, the mesh size in the stationary region is set to 6 mm, with a smooth transition applied between adjacent regions. The computational model comprising 518k cells is selected for subsequent simulations, as it achieves a favorable balance between accuracy, reproducibility, and efficiency. After mesh generation, the orthogonal grid quality exceeds 0.4, indicating that the mesh quality satisfies the requirements for reliable numerical simulations.

Figure 6: Mesh independence verification.

2.2 Monte Carlo Method

The Monte Carlo method [41] is a numerical method based on probability and statistical theory, in which solutions are obtained through random sampling. It is particularly suitable for mathematical and engineering problems with complex analytical forms, or problems that are difficult to solve directly using deterministic approaches. The fundamental principle of the Monte Carlo method is to generate a large number of random samples and perform statistical estimation of the target problem, thereby approximating the true solution in a probabilistic sense.

In Monte Carlo simulations, the generation of pseudo-random numbers that conform to prescribed distribution characteristics is critical for obtaining reliable numerical solutions. By constructing random point sets with different distribution forms, physical or geometrical problems involving stochastic effects can be effectively modeled and analyzed. Random number sequences are typically generated using pseudo-random number algorithms with good statistical properties, such as the Mersenne Twister, to ensure adequate randomness and representativeness of the samples. Through repeated sampling and statistical analysis of random point sets with different dimensionalities, geometric domains, and distribution characteristics, a quantitative framework for characterizing spatial distribution features under multiscale geometric configurations can be established. As the number of simulations and sample size increase, the stability and accuracy of the statistical results improve correspondingly.

In this study, random point sets under different distribution conditions are generated using the Monte Carlo method, and numerical estimation is performed with Ripley’s L function. This procedure provides a reference benchmark for quantitative analysis and validation of spatial distribution characteristics in complex systems.

2.3 Ripley’s L Function

In solid-liquid stirred systems, traditional methods for determining mixing time are typically based on concentration standard deviation or relative standard deviation, which rely on limited sampling points or averaged measurements. As a result, they provide only single-scale assessments of bulk uniformity and fail to capture particle aggregation, stratification, and their spatiotemporal evolution. With advances in computational fluid dynamics (CFD) and related numerical techniques, the three-dimensional coordinates of large particle populations can now be directly obtained during stirring. This development has created a demand for methods that can extract multiscale mixing characteristics from spatial point data. Ripley’s K function, a representative multi-distance spatial clustering method, offers such a multiscale statistical framework by quantifying deviations from complete spatial randomness (CSR) based on spatiotemporal information [42,43]. This enables the identification of clustered, random, or regular particle distributions across different characteristic scales, thereby providing a more reliable statistical basis for evaluating solid-liquid dispersion performance and optimizing operating conditions and flow-structure regulation.

Ripley’s K function is a second-order statistic based on pairwise distances, which simultaneously accounts for the number of points and the distances between them. This enables quantitative assessment of deviations from randomness across multiple spatiotemporal scales [44]. The theoretical space-time Ripley’s K function can be defined as the expected number of events within a cylinder of radius r and temporal half-width t, normalized by the point intensity λ r t : Kr,t=E(r,t)λrt(11) where the expected number of neighboring events within the space-time cylinder is denoted by E ( r , t ) , and the point intensity over the studied region and time period is denoted by λ rt . Conceptually, Eq. (11) characterizes the distribution pattern of space-time points by considering a reference cylinder with a radius r, a height of 2t, and a base area of πr2, and by evaluating the expected number of neighboring events located within this cylindrical domain.

Assuming that the point pattern strictly follows complete spatial randomness (CSR), the expected value of the space-time K function is equal to its theoretical value under CSR. If the evaluated value exceeds the CSR expectation, it indicates the presence of clustering within the specified spatial distance r and temporal distance t. Conversely, if the evaluated value is smaller than the CSR expectation, a dispersed space-time distribution pattern is indicated. On the basis of Eq. (12), the space-time K function can be further expressed as: K^(r,t)=A·Dn2∑i=1n∑j≠iIr,t(dij,uij)ωijvij(12) where A denotes the total number of sample points, D represents the area of the study region, dij and uij are the spatial and temporal distances between points i and j, respectively, wij and νij are the edge-correction weights in space and time, respectively, and I r , t ( d i j , u i j ) is the indicator function, which takes the value of 1 if the specified condition is satisfied and 0 otherwise. As shown in Eq. (13) [45]:

When the space-time cylindrical neighborhood centered at a given point extends beyond the boundaries of the study region, the number of neighboring points within the cylinder may be underestimated due to the absence of observations outside the domain. This gives rise to the so-called edge effect, which introduces bias into the statistical results and prevents effective analysis. Consequently, edge correction represents one of the key distinctions between Ripley’s K function and other point pattern analysis methods.

In point pattern analysis, the K function increases approximately as a power function of the spatial scale r, which causes its magnitude to grow rapidly with increasing scale. As a result, direct comparison of clustering intensity across different distance scales becomes unintuitive. To improve interpretability, Ripley [46] proposed the L-function transformation, which rescales the K function into a more readily interpretable form. Ripley’s L function is a variance-stabilized transformation of Ripley’s K function. For three-dimensional space, the transformation can be expressed as: L(r,t)=3K(r,t)4π3−r(14) here, the sign of L ( r , t ) can be used to directly identify the deviation of the point pattern from complete spatial randomness. L ( r , t ) > 0 indicates clustering at scale r, L ( r , t ) < 0 indicates inhibition, and L ( r , t ) ≈ 0 indicates consistency with complete spatial randomness.

A comparison between the scale dependence and interpretability of the K(r) function and its transformed L(r) form is presented in Fig. 7 to clarify the reason for adopting L(r) in this study. As shown in Fig. 7a, K(r) increases nonlinearly with the spatial scale r, and its value is strongly scale-dependent. Although clustering and dispersion can still be distinguished from the upward or downward deviations relative to the CSR reference, the magnitudes of these deviations cannot be compared directly over different scale ranges, which limits the intuitive interpretation of the results. By contrast, after the L(r) transformation, the CSR reference is converted into a zero baseline, as shown in Fig. 7b. In this form, positive values indicate clustering, whereas negative values indicate dispersion. This transformation reduces the effect of scale-dependent growth embedded in K(r), makes the deviations of different spatial point patterns more explicit, and improves comparability of the results across multiple scales.

For these reasons, the L(r) function is adopted in the present work rather than the original K(r) function. Since this study aims to characterize the evolution of particle distribution structures in a stirred tank over different spatial scales, a statistic with clearer physical meaning and better cross-scale interpretability is required. Compared with K(r), the L(r) function provides a more intuitive basis for identifying aggregation, dispersion, and the gradual transition to a spatially uniform state. Therefore, it is more suitable for quantitative analysis of particle distribution characteristics in confined systems such as stirred tanks.

Figure 7: (a) Schematic of the K-function describing clustering and regularity, (b) schematic of the L-function describing clustering and regularity.

3 Results and Discussion

3.1 Numerical Analysis Based on the Monte Carlo Method

A numerical model was developed to investigate particle distribution within a cylindrical container with a diameter of 16 cm and a height of 16 cm. Using the Monte Carlo method, 1000 particles with a diameter of 2 mm were generated, and their spatial coordinates were assigned within the cylindrical domain using pseudo-random sequences following uniform, normal, and binomial distributions. These point sets represent distinct three-dimensional particle distribution patterns. As shown in Fig. 8, the XY and XZ planes were selected to examine the corresponding two-dimensional cross-sectional distributions of the particles.

Figure 8: Particle distribution views: (a) uniform distribution: 3D view, XY view, and XZ view; (b) normal distribution: 3D view, XY view, and XZ view; (c) binomial distribution: 3D view, XY view, and XZ view.

The variation of the L-function deviation, L(r), with spatial scale r for the three spatial point patterns within the cylindrical domain is presented in Fig. 9. For the uniform distribution, the L(r) curve remains close to zero over the range r = 0.4–3.0 cm, indicating no significant deviation from complete spatial randomness (CSR). In contrast, the curves for the normal and binomial distributions lie markedly above the CSR baseline (L(r) = 0) across most spatial scales, implying statistically significant clustering. The binomial distribution exhibits larger positive deviations at small to intermediate scales (approximately 0.4–2.2 cm), indicating stronger local clustering. The two curves become comparable around r ≈ 2.0–2.3 cm, beyond which the normal distribution curve continues to rise, revealing more pronounced clustering at larger scales. At the largest scale considered, the binomial distribution curve shows a slight decline. Overall, the L(r) curves clearly differentiate the uniform, normal, and binomial spatial distribution patterns.

Figure 9: Variation of the L-function with scale (r) for the uniform, normal, and binomial spatial distributions.

To evaluate the dependence of particle spatial distributions in a cylindrical stirred tank on scale r and sample size N, a boundary-corrected statistical framework based on the three-dimensional Ripley’s L function was developed. The study region was a cylindrical domain with diameter D = 16 cm and height H = 16 cm. Particles were uniformly sampled according to the Mersenne Twister random-number generator for N = {500, 1000, 2000, 5000, 10,000, 20,000}, and the corresponding L(r) curves were computed over r ∈ [0.4, 3.0] cm. Based on 99 Monte Carlo simulations, the results are shown in Fig. 10. The red dashed line L(r) = 0 is taken as the CSR reference. As N increases, the curves progressively converge toward the baseline across all scales and their deviations are markedly reduced. The curves for N = 500 and N = 1000 exhibit pronounced fluctuations at small scales; for N = 2000, slight scale-dependent deviations remain. When N = 5000, the curve is already close to the baseline and shows only a small, nearly monotonic offset. For N = 10,000, the curve lies close to the baseline overall and oscillates around it. At N = 20,000, the curve almost coincides with the CSR baseline L(r) = 0 across the entire range, indicating good consistency and stability under large-sample conditions.

Figure 10: L(r) curves for different particle numbers.

To quantitatively characterize the degree of uniformity of particle distributions at different times, six sets of three-dimensional particle coordinates obtained from FLUENT simulations at t = 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 s were selected. The corresponding three-dimensional L(r) curves were computed, as illustrated in Fig. 11. At small spatial scales r, the curves are strongly influenced by local clustering and sampling noise, whereas at large scales they tend to converge due to boundary effects and scale averaging. The separation among the curves reaches the maximum at intermediate scales. Based on the location of maximum separation between the L(r) curves in Fig. 11, a characteristic discrimination scale of r = 4 cm was identified. At this scale, the differences in L(r) values at different times are most pronounced, facilitating the ranking of mixing uniformity. According to the criterion that smaller absolute values of |L(r)| indicate states closer to complete spatial uniformity, the L(r) curves exhibit a progressive downward shift with increasing time; the curve at t = 3.0 s is closest to zero. This trend indicates a continuous improvement in overall mixing uniformity over time. Therefore, in subsequent comparisons and optimization analyses, r = 4 cm is adopted as a single-scale evaluation point, providing an appropriate balance between sensitivity to mixing evolution and statistical robustness.

Figure 11: L(r) as a function of r at different times.

3.2 Mixing Uniformity of Particle Distribution in the Cylindrical Vessel

In Section 3.1, numerical analyses were conducted on three representative spatial point patterns generated via the Monte Carlo method to systematically examine the ability of Ripley’s L function to discriminate among particle distribution modes and to determine its effective discrimination range and multiscale analysis criteria. To further elucidate the spatial evolution of the particle phase during stirring, this statistical framework is applied to a solid-liquid mixing process in a cylindrical stirred tank simulated using the DEM-VOF approach. By analyzing the time-dependent statistical characteristics of particle spatial distributions, mixing uniformity is quantitatively evaluated and its dynamic evolution is characterized.

As illustrated in Fig. 12, a four-blade turbine impeller is selected as a representative case. Three-dimensional particle coordinates are extracted at 0.1 s intervals, yielding a total of 50 datasets. The corresponding L(r) values at the spatial scale r = 4 cm are computed to analyze the temporal evolution of solid-liquid mixing uniformity under this operating condition. The time-resolved maps of particle spatial distributions during the solid-liquid stirring process are shown in Fig. 12, illustrating the dispersion behavior of solid particles over 0–5 s under the action of the four-blade turbine impeller B-90°. The instantaneous particle configurations can be directly correlated with the L(r) trajectory evaluated at r = 4 cm.

The L(r) trajectory for the B-90° impeller over 0–5 s can be divided into three stages (Fig. 12). During the initial stage (0–0.3 s), the L(r) values are markedly high, ranging from 37 to 45. As stirring proceeds from 0.3 to 1.0 s, L(r) decreases from approximately 37 to 30, although the values remain relatively large. In the period from 1.0 to 3.7 s, L(r) continues to decline with time to about 5. Finally, during 3.7–5.0 s, only slight fluctuations around 6 are observed, indicating that the system has approached a quasi-steady mixing regime.

Figure 12: L(r) curve for the B-90° impeller.

The values of L(r) faithfully reflect the instantaneous state of the particle system. As shown in Fig. 13, during 0–0.3 s the particles are mainly accumulated at the bottom of the stirred tank, forming a shallow packed bed and exhibiting pronounced stratification. In this stage, L(r) is large and positive, indicating a strongly clustered spatial arrangement of the particles. As stirring proceeds (0.3–1.0 s), the radial jet generated by the turbine impeller and the associated vertical circulation gradually develop, entraining a large number of particles into the main circulation region. The suspension height and radial spread increase rapidly, and the corresponding L(r) curve drops sharply, revealing a pronounced reduction in spatial clustering and an accelerated mixing process. Between 1 and 3 s, the particles are already distributed throughout most of the tank, although noticeable local concentration gradients persist. During this period, L(r) decreases slowly with time, indicating a transition from macroscopic blending toward microscopic homogenization. At around t = 3.7 s, L(r) reaches a relatively low level and exhibits only small fluctuations. Combined with the observation that, between 3 and 4 s, the particles are stably dispersed without obvious zonation, this suggests that the system has approached an approximately uniform state, marking the onset of a steady mixing regime.

Figure 13: Temporal evolution of particle spatial distribution under the B-90° operating condition.

From the perspective of flow mechanism, the three-stage evolution of L(r) is directly related to the force balance acting on the particles and the development of the flow field in the stirred tank. In the initial stage, the particles are mainly governed by gravity and remain deposited at the tank bottom, while the impeller-induced flow is not yet fully established. As a result, the system exhibits pronounced bottom accumulation and stratification, corresponding to relatively high L(r) values. As stirring proceeds, the radial jet generated by the rotating impeller gradually strengthens and induces axial circulation. Under the combined effects of entrainment, lifting, and redistribution, the bottom particles rapidly enter the main circulation region. The original particle-rich zones are continuously disrupted, leading to a sharp decrease in L(r). In the later stage, particles are distributed throughout most of the tank, and the dominant process shifts from large-scale redistribution to the gradual attenuation of local concentration fluctuations. Accordingly, the rate of decrease in L(r) becomes significantly smaller and gradually approaches a stable level. These results indicate that the variation of L(r) not only reflects the evolution of particle spatial distribution, but also characterizes the dynamic process from bottom accumulation to suspension and finally to quasi-uniform dispersion under the evolving flow field.

Comparison between Fig. 12 and Fig. 13 demonstrates that, at the spatial scale of r = 4 cm, the mixing state and mixing time evaluated using Ripley’s L function are in good agreement with the particle distribution evolution obtained from the FLUENT simulations. This agreement confirms that the proposed statistical method possesses clear physical interpretability and strong practical applicability.

3.3 Evaluation of Mixing Performance of Different Impellers Based on Ripley’s L Function

To quantitatively evaluate the mixing performance of different impellers, Ripley’s L function is employed to describe the spatial uniformity of particle distribution during stirring. By computing and comparing L(r) values for six impeller configurations under identical simulation conditions, differences in the mixing mechanisms and performance are elucidated from the perspectives of spatial statistics and flow structure. The temporal evolution of L(r) for the three pairs of impeller configurations: A-45° vs. A-90°, B-45° vs. B-90°, and C-45° vs. C-90° is presented in Fig. 14a–c. In all three cases, the L(r) values for the 45° pitched-blade impellers decrease more rapidly than those for the corresponding 90° impellers, indicating faster disruption of particle clusters and more efficient mixing.

Specifically, as shown in Fig. 14a, the L(r) values under both operating conditions exhibit an initial rapid decline followed by stabilization, indicating that particle aggregates are progressively broken up and the system becomes increasingly uniform. For the A-45° impeller, L(r) first drops to a low and nearly steady level at t = 3.6 s, suggesting that the system reaches a relatively uniform state at this time. In contrast, the A-90° impeller achieves a comparable level only at t = 4.3 s, and its L(r) curve remains consistently above that of the A-45° impeller throughout the entire process. This indicates that the A-45° impeller provides more effective particle dispersion and superior mixing performance. Similarly, the B-45° impeller reaches a relatively uniform state at t = 2.4 s, whereas the B-90° impeller does so at t = 3.7 s. Likewise, the C-45° impeller reaches a relatively uniform state at t = 3.2 s, while the C-90° impeller reaches it at t = 3.9 s. Taken together, for all three impeller types, the 45° pitched-blade cases not only reduce L(r) more rapidly but also maintain it at a lower level (Fig. 14a–c). This further confirms that pitched-blade geometries offer a clear advantage over straight blades in promoting a uniform particle distribution. Fig. 14d shows the temporal evolution of L(r) for the three impellers under the 45° operating condition. L(r) for B-45° decays the fastest, C-45° exhibits the second-fastest decay with a slightly higher final plateau, and A-45° displays the slowest decrease in L(r). Throughout the entire process, the L(r) curves for B-45° and C-45° remain below that of A-45°.

The above results reveal differences in the decay rate of L(r) as well as distinct kinetic behaviors of the solid-liquid two-phase flow generated by the three impeller pairs. For impeller A, the 45° pitched blades introduce axial pumping in addition to radial discharge, promoting upward transport of particles from the lower accumulation zone and enhancing large-scale circulation throughout the vessel. Consequently, particle clusters formed near the tank bottom are progressively broken up, and the corresponding L(r) decreases earlier than in the A-90° case. In contrast, the A-90° impeller is dominated by a radial jet, so particle entrainment is concentrated near the impeller plane, while the upper and lower regions are refreshed more slowly. This weakens global solid-liquid exchange and delays homogenization. For impeller B, the difference between the 45° and 90° configurations is more pronounced. The B-45° impeller combines strong discharge capacity with effective axial circulation, thereby generating both intense local shear and efficient bulk transport. This dual effect accelerates the breakup of dense particle clusters and rapidly redistributes particles over the entire tank volume, resulting in the fastest decrease in L(r) among all cases. By contrast, although the B-90° impeller produces a relatively strong jet in the impeller region, the absence of sufficient axial pumping limits particle transfer between the lower bed and the upper liquid domain. Therefore, local dispersion is strong, but global homogenization is slower.

For impeller C, the hydrofoil geometry yields a smoother and more stable flow. In the C-45° case, the inclined hydrofoil blades generate uniform axial circulation with lower hydraulic resistance, favoring continuous suspension and redistribution of particles. This produces a clear reduction in L(r), although the decay rate is slightly slower than that of B-45° due to weaker local shear intensity. In the C-90° case, the flow becomes more radial and the circulation loop weakens, delaying particle lifting and vertical exchange. Overall, the kinetic response shown in Fig. 14 demonstrates that the particle distribution evolution is governed by the combined effects of local dispersion intensity and global circulation capacity. The 45° impellers outperform the 90° impellers by simultaneously enhancing particle suspension, axial transport, and cross-regional exchange, whereas the 90° impellers mainly intensify local radial motion.

Figure 14: L(r) curves for the six impeller configurations: (a) comparison of L(r) for the A-45° and A-90° impellers; (b) comparison of L(r) for the B-45° and B-90° impellers; (c) comparison of L(r) for the C-45° and C-90° impellers; (d) comparison of L(r) for the A-45°, B-45°, and C-45° impellers.

These results indicate that, at the same blade angle, the solid-liquid mixing performance is highest for impeller B, followed by impeller C, with impeller A being the least effective. This ranking is consistent with the separate comparisons between the 45° and 90° cases, further confirming the sensitivity of the proposed evaluation method to differences in impeller geometry.

To quantitatively characterize the mixing homogeneity from the perspective of particle concentration fluctuations, the relative standard deviation (RSD) [47] of particle counts in spatial cells was introduced as a supplementary evaluation index and compared with L(r), as shown in Fig. 15. Since RSD reflects the dispersion of particle counts in different spatial cells relative to the mean, a smaller RSD indicates a more uniform particle concentration distribution. Therefore, the temporal evolution of RSD describes the transition of the system from an initially non-uniform state to a uniform one.

Figure 15: RSD curves of particle concentration distribution for the six impeller configurations: (a) comparison of C-45°, B-45°, and A-45°; (b) comparison of C-90°, B-90°, and A-90°.

As shown in Fig. 15, the RSD curves of all six impeller configurations exhibit a common trend: a rapid decline from an initial peak value, followed by a slower decay, and finally small-amplitude fluctuations. This trend is consistent with that of L(r), indicating that, at the initial stage of stirring, both particle agglomerates and large-scale concentration inhomogeneities are rapidly reduced. As the mixing process proceeds, the system gradually evolves from an overall redistribution stage to a fluctuation-dominated stage characterized by local dispersion. Further comparison shows that, for the same impeller type, the 45° configuration generally exhibits a faster overall decay in RSD and a lower minimum value than the corresponding 90° configuration. This indicates that the 45° impeller has a greater advantage in reducing particle concentration fluctuations and improving the spatial uniformity of particle distribution. Moreover, the optimal mixing times determined from the RSD curves for different impeller types are consistent with those obtained from the L(r) curves; the time at which each impeller configuration reaches its minimum RSD value corresponds well to the time at which L(r) reaches its minimum, as summarized in Table 5.

L(r) and RSD differ in statistical basis and descriptive perspective—RSD reflects the overall fluctuation characteristics of particle concentration distribution, whereas L(r) emphasizes the structural evolution of particle spatial distribution. Despite these differences, both indices show good agreement in the identification of mixing states. In particular, both methods are highly consistent in describing the evolution trend of mixing, the relative order in which different impeller configurations approach homogenization, and the determination of the final homogeneous state. These results confirm that the use of L(r) to characterize particle mixing homogeneity is both reliable and reasonable.

The observed quantitative differences can be attributed to the distinct flow structures generated by the different impellers. The flow fields produced by the conventional four-blade impeller, the four-blade turbine impeller, and the four-blade wide hydrofoil impeller at installation angles of 45° and 90° are shown in Fig. 16. All three 45° pitched-blade configurations develop a dominant axial circulation loop spanning the entire tank, indicating favorable global transport. Among these, the B-45° impeller generates the most pronounced axial flow, with streamlines extending over a wide region and forming well-defined circulation paths, resulting in superior mass-transfer and mixing efficiency. In contrast, the flow structures associated with the three 90° straight-blade impellers display typical radial-flow characteristics, in which the high-velocity region is largely confined near the impeller and energy dissipation is relatively localized. This flow pattern results in significantly lower velocities in the upper and lower regions of the tank and insufficient streamline coverage, thereby promoting the formation of particle hold-up zones and hindering uniform mixing.

From the perspective of flow-generation mechanisms, the 45° pitched-blade impeller induces a strong axial circulation loop throughout the tank due to blade inclination, which is a key contributor to its superior mixing efficiency. In contrast, the 90° straight-blade impeller is dominated by radial jets, and the absence of effective axial transport limits overall mixing performance. The observed differences among the three impellers mainly arise from their geometries and the resulting flow patterns. The turbine impeller (impeller B) generates intense radial jets and high turbulence, which rapidly break up and disperse particles in the impeller region. Consequently, its dispersion capability exceeds that of the conventional four-blade impeller (impeller A), which produces weaker axial circulation and limited tip-shear. These characteristics promote particle re-agglomeration and result in slower homogenization. By contrast, the four-blade wide hydrofoil impeller (impeller C) with curved blades induces stable, low-resistance axial flow, suppresses flow separation, and enhances pumping capacity. This accelerates particle dispersion and suspension, allowing L(r) to reach a steady level at an earlier stage. Overall, the turbine and wide hydrofoil impellers, relying on turbulence-intensified dispersion and efficient axial pumping, respectively, outperform the conventional four-blade impeller in promoting uniform particle distribution.

Figure 16: Flow fields and particle trajectories at the time of uniform mixing for the six impeller configurations: (a) A-45°; (b) A-90°; (c) B-45°; (d) B-90°; (e) C-45°; (f) C-90°.

To further quantify the flow fields shown in Fig. 16, the average liquid-phase velocity within the same slice region was calculated. The results show that the average velocities for A-45°, A-90°, B-45°, B-90°, C-45°, and C-90° are 0.53, 0.60, 0.49, 0.56, 0.47, and 0.49 m/s, respectively. Overall, the three 90° straight-blade impellers exhibit higher slice-averaged velocities than their 45° pitched-blade counterparts, indicating that the 90° configurations generate stronger local flow near the impeller. However, when combined with the velocity contours and streamline distributions in Fig. 16, it is evident that the high-velocity regions of the 90° straight-blade impellers are concentrated near the blades, showing a typical radial-jet feature. The energy distribution is therefore localized, while the flow in the upper and lower regions of the tank remains weak, preventing the establishment of an effective circulation loop covering the entire vessel. By contrast, although the 45° pitched-blade impellers exhibit slightly lower slice-averaged velocities, they generate more pronounced axial circulation loops and broader streamline coverage, resulting in stronger overall transport capacity. Consequently, they are more effective for particle redistribution and uniform mixing throughout the tank.

In particular, although the average velocity of the B-45° impeller (0.49 m/s) is lower than those of A-90° and B-90°, its flow structure is the most conducive to forming a tank-spanning circulation loop, thereby achieving the best mixing performance. Thus, mixing performance is determined not only by the magnitude of local velocity, but also by the organization of the flow field and its ability to promote global transport. In addition, although the type C impellers show relatively low average velocities, their velocity distributions are more uniform, and the transition between local high-velocity and low-velocity regions is smoother, indicating better flow stability. In particular, the C-45° impeller benefits from the low-resistance and efficient axial pumping effect generated by its wide-blade hydrofoil structure. This structure produces stable overall circulation flow under lower energy intensity, thereby alleviating local flow separation and particle retention. Although its particle dispersion rate is not as high as that of the B-45° impeller, it shows clear advantages in sustaining particle suspension and maintaining balanced spatial distribution.

These flow-field analysis results are consistent with earlier findings based on Ripley’s L function and further corroborate the qualitative observations of the flow structures. Taken together, the flow-structure and particle-distribution perspectives jointly demonstrate the superior mixing performance of the 45° pitched-blade impellers and provide a sound theoretical basis for impeller selection and optimization. Thus, the present study confirms that Ripley’s L function is an effective quantitative tool for discriminating and comparing the mixing performance of different impeller designs.

4 Conclusions

In this study, the feasibility and effectiveness of using Ripley’s L function to characterize three-dimensional particle mixing states were systematically validated through a Monte Carlo method and DEM-VOF simulations. The main conclusions are summarized as follows:

• (1)A three-dimensional Ripley’s L function from spatial statistics was introduced as a quantitative tool for evaluating solid-liquid mixing homogeneity, and a corresponding three-dimensional spatial statistical analysis framework was established. The Monte Carlo method demonstrated that the L(r) curves can clearly distinguish among three typical spatial distribution patterns, namely uniform, normal, and binomial distributions. This confirms the effectiveness of the proposed method in identifying differences in particle clustering and dispersion.
• (2)Ripley’s L function was found to effectively capture the temporal evolution of particle distribution during stirring. The separation among the L(r) curves at different times was most pronounced over an intermediate range of spatial scales, and a characteristic scale of r = 4 cm was identified as an appropriate evaluation point. At this scale, the decrease in L(r) corresponded well to the physical evolution of particles in the stirred tank from clustering to dispersion.
• (3)Comparison of Ripley’s L function with the concentration-based RSD index confirmed the reliability and rationality of the spatial-statistical method based on Ripley’s L function. The optimal mixing times determined from the L(r) curves for different impeller configurations were consistent with those obtained from the RSD curves. Specifically, the B-45° impeller exhibited the best mixing performance, reaching an approximately uniform state at about 2.4 s; the C-45° impeller ranked second and reached a similar state at about 3.2 s; while the A-45° impeller required 3.6 s. For the 90° configurations, the corresponding times were about 3.7 s for B-90°, 3.9 s for C-90°, and 4.3 s for A-90°. These results indicate that, at the same blade angle, the overall performance of different impeller types in promoting particle dispersion and reducing spatial non-uniformity follows the order B > C > A, and the 45° impellers consistently outperform the 90° impellers.