# Numerical Study of Percolation and Seepage Behaviors in Ion-Adsorption-Type Rare Earth Ore Leaching Process

Authors: E. Dianyu; Zhongfang Su; Jia Zeng; Liuyimei Yang; Jing Li; Qiang Xu; Lin Wang; Jiaxin Cui
Journal: Mathematical Problems in Engineering
 (2023)

---

## Abstract

Ionic rare earth ore is a type of featured rare earth ore in China. Its mining process suffers from a long leaching cycle and considerable consumption of leaching agents. Improving mining efficiency requires a sound physical understanding of the leaching process. In this study, the CFD-based numerical model is used to analyze the physical process of leaching through porous media formed by particles. The simulation results indicate that a lower packing porosity and smaller particles packed granular porous medium result in much larger energy dissipation during seepage, and the energy dissipation increases with seepage velocity. It is found that when the seepage velocity increases to a certain high value, the energy dissipation exceeds the value predicted by Darcy’s law, which is mainly caused by liquid turbulence. Additionally, the effect of particle shape is examined. The results show that the granular medium composed of prolate particles causes larger energy dissipation than oblate particles, and spherical particles play the least role. This phenomenon may result from the particle shape affecting the area of the frontal contact surface between particles and liquid. The results provide new insights into the fundamental understanding of percolation and seepage behaviors in the ion-adsorption-type rare earth ore leaching process.


---

## Body

## 1. Introduction

Ionic rare earth ore, also known as weathering crust leaching rare earth ore, is a type of featured rare earth ore in China. It plays a significant role in high-tech and new material industries due to its various advantages, e.g., diverse sources, abundant reserves, and including all the light, medium and heavy rare earth, particularly rich in terbium, dysprosium, and europium [1]. Ionic rare earth ore is formed by the adsorption of rare earth elements in the form of ions on clay minerals [2]. There are considerable reserves of rare earth elements in the earth’s crust. However, due to its unique properties and the limitation of the technical level, although lots of rare earth elements have been explored in many areas, they cannot be well exploited and utilized. In the mining process, several key issues are a concern, such as the long leaching cycle and the large consumption of leaching agents. Therefore, it is of great significance to improve the understanding and optimize the extraction process of rare earth minerals.The in-situ leaching process incorporates a typical fluid seepage in a porous medium accompanied by an ion exchange reaction, which is influenced by the hydrodynamics in the porous medium and the leaching reaction kinetics. Many researchers have conducted extensive experiments to improve the leaching efficiency [3, 4]. To date, previous experimental work mainly generates input and output data but lacks analyses of fluid flow, mass transfer, pressure, and velocity distributions. In terms of numerical studies, the lattice Boltzmann method is mostly adopted in previous numerical research, which is regarded as a more suitable method for a small number of lattices but not for engineering calculation [4]. Zhao et al. [5] used Fluent software to simulate the migration of the water phase in three different filling structures under different initial water velocities and dynamic viscosities and indicated between the fixed particles (soil skeleton) in the filling structure, the looser structure, and the greater influence of gravity on the water phase seepage. Wang et al. [6] study the grouting seepage of microfractures under different grouting pressures and fracture opening conditions, and the grouting pressure and fracture opening changes are nonlinearly attenuated along the grout seepage direction. Li et al. [7] examined the effects of particle size on flow regimes and the validity of Darcy law in the seepage process. Luo et al. [8] carried out column leaching experiments to consider the impact of permeability in the rare earth leaching process and analyzed the mechanism of the change of permeability coefficient. Dixson and Hendrix [9] proposed a mathematical model for the leaching of various solid reactants from ore pellets heap assuming that the size, density, and packing porosity of the spherical particles were the same and the leaching process was isothermal and irreversible. Ma et al. [10, 11] established a one-dimensional radial three-phase flow model of water-rock-silt and obtained the temporal-spatial distribution of hydraulic properties by numerical simulation. Li et al. [12] established a flow model for the seepage process under gravity in the in-situ leaching process and simulated the vertical infiltration process of the leachate, which provided knowledge for further exploration of the flow mechanism of the leachate. Li et al. [13] quantitatively investigated the pore structure and hydraulic properties of the broken rock mass by the nuclear magnetic resonance technology and the non-Darcy models. Additionally, a logistic regression model was proposed to characterize the effects of compressive stress and GSD on the permeability and the characteristic parameter of nonlinear flow in the broken rock mass. Huang et al. [14] investigated the effect of surface roughness on laminar flow in macro tubes and conducted head drop experiments in six organic glass tubes with different relative roughness. In general, although the applicable critical value of Darcy’s law has been measured in previous studies, the composition of experimental materials is relatively simple, and the research is limited to simplified situations, which are difficult to represent the most practical scenarios. The related work is mainly based on experiments, which can qualitatively analyze the process, while the relevant mathematical models mostly focus on the effect of packing particle size, ignoring some factors such as particle profile.This study would apply a single-phase flow model and the Euler–Euler two-fluid model to simulate the seepage process of the fluid through a porous medium formed by particle packing. The effects of most related factors, such as seepage velocity and particle properties, on seepage efficiency will be analyzed, aiming to provide a fundamental understanding of this leaching process.


## 2. Model Description

In this work, simulations were conducted using the commercialized Fluent software package. The effects of porosity, particle size, and particle shape on energy loss during seepage were investigated. More specifically, the Euler-Euler two-fluid model, which is widely used to describe the fluid flow in porous media, was adopted in the seepage process simulations, while the single-phase flow model was employed to examine the effect of particle size and shape for better capturing pore scale information. In the two-fluid model, the solid and liquid phases are considered fully interpenetrating continua coupled with an interaction term. Since the volume occupied by one phase can no longer be occupied by other phases, the concept of phase volume fraction is introduced. The volume fraction is a continuous function of time and space, and the sum of the volume fractions of all phases equals 1.0. The two phases of transport phenomena are resolved with similar mass and momentum conservation equations. For the liquid phasel, the mass conservation equation is given as follows:(1)∂∂tαlρl+∇⋅αlρlvl=∑p=1nm˙ls−m˙sl+Sl,where αl, ρl, and vl represent the volume fraction, density, and velocity of the liquid phase, respectively; m˙ls stands for the mass transfer from the liquid phase to the solid phase; m˙sl denotes the mass transfer from the solid phase to the liquid phase; and Sl is the source term. The momentum conservation equation of the liquid phase can be expressed as follows:(2)∂∂tαlρlvl+∇⋅αlρlvlvl=−αl∇p−∇pl+∇⋅τl+αlρlg−m˙lsvls,where p is pressure, τl is the stress-strain tensor of the liquid phase, and vls is the relative velocity of the liquid and solid phases.Similarly, the mass and momentum conservation equations of the solid phase are given as follows:(3)∂∂tαsρs+∇⋅αsρsvs=∑p=1nm˙sl−m˙ls+Ss,∂∂tαsρsvs+∇⋅αsρsvsvs=−αl∇p−∇ps+∇⋅τs+αsρsg−m˙lsvls.As for the interaction force, namely, the drag force, the Gidaspow model is adopted. It is based on the Wen-Yu model and combines it with the Ergun equation [15]. When αl > 0.8, the liquid-solid exchange coefficient equation is described as follows:(4)Ksl=18αlRes1+0.15αlRes0.687αsαlρlvs−vldsαl−2.65.Whenαl < 0.8, the liquid-solid exchange coefficient equation is given as follows:(5)Ksl=150αs1−αlμlαlds2+1.75ρlαsvs−vlds.The standardk-epsilon model, which is suitable for fluids with fully developed turbulence [16], is adopted for turbulence description. The expression of the turbulence equations, including the turbulent kinetic energy k and the turbulent dissipation rate ε are given as follows:(6)∂∂tρk+∂∂xiρkui=∂∂xjμ+μfσk∂k∂xj+Gk+Gb−ρε−YM,∂∂tρε+∂∂xiρεui=∂∂xjμ+μtσε∂ε∂xj+C1εεkGk+C3εGb−C2ερε2k,where μt is the turbulent viscosity; Gk is the turbulent kinetic energy; Gb is the turbulent kinetic energy caused by buoyancy; YM is the change in the total dissipation rate caused by pulsating expansion; μ is the fluid viscosity; σk and σε are the turbulent Prandtl number of k and ε, respectively; and ui is the average fluid velocity. In general, C1ε = 1.44, C2ε = 1.92, C3ε = 0.99, σk = 1.0, and σε = 1.3 as reported elsewhere [17, 18].


## 3. Simulation Conditions

### 3.1. Geometric Structure and Meshing

In this study, the simulation simplifies the physical model (a tube, referring to a previous experimental study [19]) with a 2D cross-section to improve the computing efficiency in considering the effect of pack bed porosity. The calculation domain is in the dimensions of 350 mm (length) × 100 mm (width, tube diameter) as shown in Figure 1. To speed up the calculation convergence and eliminate the inlet and outlet effects, two buffer sections with a length of 50 mm were set at the inlet and outlet, respectively. Note that the size of the established numerical model is slightly shorter than that of the experimental model. Therefore, the sufficient and stable flow of the fluid phase should not be affected. Additionally, the corresponding 3D geometry is given in Figure 2. In particular, the tube diameter is 47 mm and the length of the middle section of the porous medium is 100 mm. The inlet and outlet buffer sections were set up with a length of 150 mm.Figure 1
(a) 2D geometry of the simulation domain and (b) grid arrangement for the study on the effect of porosity.
(a)(b)Figure 2
(a) 3D geometry of the simulation domain and (b) grid arrangement for the study on the effect of the particle size and shape.
(a)(b)The CFD mesh is critical to the reliability and accuracy of the simulation results. Figure3 shows the effect of grid number on the simulated outlet velocity. According to this grid independence test, the grid numbers with the values of 8,750 and 371,110 were selected in the 2D and 3D simulations, respectively, to study the effect of porosity, particle size, and shape.Figure 3
The effects of different grid numbers on the simulated outlet velocities: (a) 2D simulation to study the effect of porosity; (b) 3D simulation to study the effect of particle size and shape.
(a)(b)


### 3.2. Boundary Condition

The “Velocity Inlet” boundary condition was adopted, with fluid speed values of 0.001-0.011 m/s. The outlet adopted the widely used “Pressure Outlet” boundary condition. Additionally, the operating temperature was ambient temperature, set to 298 K. The pressure reference algorithm was selected, and the implicit solver was used for the transient solution. The time step size was 10−3 s. Pressure-speed coupling chose the SIMPLE algorithm, and the difference scheme adopted second-order upwind; the wall was processed according to the “no-slip” boundary condition. The Gidaspow model was adopted as drag force calculation, and the residuals for simulation calculation were less than 10−3.


## 3.1. Geometric Structure and Meshing

In this study, the simulation simplifies the physical model (a tube, referring to a previous experimental study [19]) with a 2D cross-section to improve the computing efficiency in considering the effect of pack bed porosity. The calculation domain is in the dimensions of 350 mm (length) × 100 mm (width, tube diameter) as shown in Figure 1. To speed up the calculation convergence and eliminate the inlet and outlet effects, two buffer sections with a length of 50 mm were set at the inlet and outlet, respectively. Note that the size of the established numerical model is slightly shorter than that of the experimental model. Therefore, the sufficient and stable flow of the fluid phase should not be affected. Additionally, the corresponding 3D geometry is given in Figure 2. In particular, the tube diameter is 47 mm and the length of the middle section of the porous medium is 100 mm. The inlet and outlet buffer sections were set up with a length of 150 mm.Figure 1
(a) 2D geometry of the simulation domain and (b) grid arrangement for the study on the effect of porosity.
(a)(b)Figure 2
(a) 3D geometry of the simulation domain and (b) grid arrangement for the study on the effect of the particle size and shape.
(a)(b)The CFD mesh is critical to the reliability and accuracy of the simulation results. Figure3 shows the effect of grid number on the simulated outlet velocity. According to this grid independence test, the grid numbers with the values of 8,750 and 371,110 were selected in the 2D and 3D simulations, respectively, to study the effect of porosity, particle size, and shape.Figure 3
The effects of different grid numbers on the simulated outlet velocities: (a) 2D simulation to study the effect of porosity; (b) 3D simulation to study the effect of particle size and shape.
(a)(b)


## 3.2. Boundary Condition

The “Velocity Inlet” boundary condition was adopted, with fluid speed values of 0.001-0.011 m/s. The outlet adopted the widely used “Pressure Outlet” boundary condition. Additionally, the operating temperature was ambient temperature, set to 298 K. The pressure reference algorithm was selected, and the implicit solver was used for the transient solution. The time step size was 10−3 s. Pressure-speed coupling chose the SIMPLE algorithm, and the difference scheme adopted second-order upwind; the wall was processed according to the “no-slip” boundary condition. The Gidaspow model was adopted as drag force calculation, and the residuals for simulation calculation were less than 10−3.


## 4. Results and Discussion

### 4.1. Model Validation

To validate the numerical model, the variation of the hydraulic gradient with seepage velocity was depicted as shown in Figure4. The hydraulic gradient indicates the ratio of the head drops along the seepage path to its length, given as follows:(7)J=∆H∆L,where ∆H is the head drop of the fluid flowing through the porous media region, and ∆L is the length of the porous media region. It can be seen that the simulated results are in good agreement with the experimental results [14], demonstrating the effectiveness of the present CFD model.Figure 4
The developed (a) 2D and (b) 3D model predictions against experimental measurements.
(a)(b)


### 4.2. Effect of Porous Medium Porosity

Ionic rare earth ore is normally packed as a porous medium. Normally, its porosity varies according to its different categories. To investigate its influence, the porosity of 0.37, 0.40, 0.43, 0.46, and 0.49 is adopted in the simulation. The pressure contours under different porosity are depicted in Figure5. The results show that the pressure distribution varies with porosity under a constant seepage velocity of 0.001 m/s. Figure 6 shows the hydraulic gradient increases with seepage velocity, which results from the small vortices formed during the fluid flowing through the porous medium. As the seepage velocity increases, the resulting vortex increases correspondingly, which is a key factor for energy dissipation. Additionally, under the same seepage velocity, a larger porosity results in a smaller hydraulic gradient. This is because porosity affects the resistance of fluid flow. More specifically, the flow resistance is larger with a smaller porosity, which requires a larger hydraulic gradient to achieve a certain flow rate.Figure 5
Contours of the pressure under different porosity. The black lines close to the inlet and outlet are the boundaries between the buffer section and the porous medium.Figure 6
Variation of energy dissipation as a function of seepage velocity under different porosity.


### 4.3. Effect of Packed Bed Particle Size

For the size of most ionic rare earth ores, around two-thirds of the ore particles’ diameters are smaller than 7.75 mm [1]. In this study, five different values of spherical particle diameter are adopted for comparison analyses, namely, 1.0 mm, 1.5 mm, 2.0 mm, 2.5 mm, and 3.0 mm. The pressure contours of different particle-size packed beds are shown in Figure 7. The smaller the particle size, the larger the pressure drop. The variation of the hydraulic gradient with seepage velocity is shown in Figure 8. The results indicate that under the same seepage velocity, the smaller the particle size, the larger the corresponding hydraulic gradient, which indicates that the fluid flow in a medium packed with fine particles is more difficult than that packed with coarse particles.Figure 7
Contour of the pressure drop under different particle sizes.Figure 8
Variation of energy dissipation with seepage velocity under different particle sizes.According to Darcy’s law, the cross-sectional areaA, the seepage flow Q through A, the head drop (h1 − h2), and the seepage path length L are in the following relationship [19]:(8)Q=K∙A∙h1−h2L,which is simplified as follows:(9)vs=K∙J,where K is the permeability coefficient, determined by the solid and fluid materials; vs is the seepage velocity; and J is the hydraulic gradient. The corresponding curves based on Darcy’s law are shown in Figure 8. With the increase in seepage velocity, the seepage phenomenon for each of the five types of particles gradually deviates from Darcy’s law. This is because the seepage process can be regarded as laminar fluid flow under a small seepage velocity. With a larger seepage velocity, the seepage velocity and hydraulic gradient are no longer linear, and the seepage process changes into a turbulent fluid flow. A similar deviation from Darcy’s law was also observed at high velocity [20, 21].


### 4.4. Effect of Packed Bed Particle Shape

Since the ore particles are generally nonspherical ones in practical industries, this work also examines the effect of particle shape on the seepage process, such as prolate particles and oblate particles. Note that each particle with different shapes (ellipsoid aspect ratio) has the same volume in this study. The aspect ratio of a spherical particle is 1.0, while it is 1.25 for a prolate particle and 0.82 for an oblate particle as shown in Figure9. Figure 9(a) represents spherical particles with a radius of 0.57 mm. The other ellipsoidal particles’ projections are perpendicular to the inlet direction presenting a circle with a radius of 0.61 mm and an ellipsoid with a long/short semiaxis of 0.75 and 0.5 mm, respectively.Figure 9
Packing of (a) spherical particles; (b) oblate particles, and (c) prolate particles.
(a)(b)(c)Figure10 shows the variation of energy dissipation with seepage velocity for the systems packed with differently shaped particles. It is shown that the energy dissipation is not significantly influenced by the particle shape under a small seepage velocity. However, the seepage resistance for a spherical particle system is smaller than that of a nonspherical particle system at a larger seepage velocity. For the nonspherical systems, the prolate particle system has larger seepage resistance and energy dissipation than that of the oblate particle system. The difference can also be observed in the contours of fluid pressure on the particle surfaces shown in Figure 11 and the corresponding reasons are discussed later in the article.Figure 10
Variation of energy dissipation with seepage velocity under different particle shapes.Figure 11
Contours of the pressure on varied particle surfaces (seepage velocity is 0.011 m/s).The mean curvature, referring to the ratio of the angle between the tangents at the ends of the arc to the length of the arc, and the frontal contact surface area of the particle are calculated for different particle shapes as shown in Table1. Since the selected arc corresponds to the long axis, the angle between the ends of the arc is 180°. The frontal contact surface area in the situation of this work equals half of the particle surface area. As listed in Table 1, the oblate and spherical particles have almost the same curvature. In contrast, the energy dissipation under larger seepage velocities for oblate particles is larger than that for spherical particles. This may reveal that the curvature does not determine the energy dissipation during leaching. The difference in frontal contact surface area corresponds better with the difference in energy dissipation, indicating that the frontal contact surface area may be the main factor determining the energy dissipation.Table 1
Curvature and frontal contact surface area for each particle shape.
Particle shapeMean curvatureFrontal contact surface area (mm2)Spherical particle0.89534.0828Oblate particle0.89554.1138Prolate particle1.03504.1888


## 4.1. Model Validation

To validate the numerical model, the variation of the hydraulic gradient with seepage velocity was depicted as shown in Figure4. The hydraulic gradient indicates the ratio of the head drops along the seepage path to its length, given as follows:(7)J=∆H∆L,where ∆H is the head drop of the fluid flowing through the porous media region, and ∆L is the length of the porous media region. It can be seen that the simulated results are in good agreement with the experimental results [14], demonstrating the effectiveness of the present CFD model.Figure 4
The developed (a) 2D and (b) 3D model predictions against experimental measurements.
(a)(b)


## 4.2. Effect of Porous Medium Porosity

Ionic rare earth ore is normally packed as a porous medium. Normally, its porosity varies according to its different categories. To investigate its influence, the porosity of 0.37, 0.40, 0.43, 0.46, and 0.49 is adopted in the simulation. The pressure contours under different porosity are depicted in Figure5. The results show that the pressure distribution varies with porosity under a constant seepage velocity of 0.001 m/s. Figure 6 shows the hydraulic gradient increases with seepage velocity, which results from the small vortices formed during the fluid flowing through the porous medium. As the seepage velocity increases, the resulting vortex increases correspondingly, which is a key factor for energy dissipation. Additionally, under the same seepage velocity, a larger porosity results in a smaller hydraulic gradient. This is because porosity affects the resistance of fluid flow. More specifically, the flow resistance is larger with a smaller porosity, which requires a larger hydraulic gradient to achieve a certain flow rate.Figure 5
Contours of the pressure under different porosity. The black lines close to the inlet and outlet are the boundaries between the buffer section and the porous medium.Figure 6
Variation of energy dissipation as a function of seepage velocity under different porosity.


## 4.3. Effect of Packed Bed Particle Size

For the size of most ionic rare earth ores, around two-thirds of the ore particles’ diameters are smaller than 7.75 mm [1]. In this study, five different values of spherical particle diameter are adopted for comparison analyses, namely, 1.0 mm, 1.5 mm, 2.0 mm, 2.5 mm, and 3.0 mm. The pressure contours of different particle-size packed beds are shown in Figure 7. The smaller the particle size, the larger the pressure drop. The variation of the hydraulic gradient with seepage velocity is shown in Figure 8. The results indicate that under the same seepage velocity, the smaller the particle size, the larger the corresponding hydraulic gradient, which indicates that the fluid flow in a medium packed with fine particles is more difficult than that packed with coarse particles.Figure 7
Contour of the pressure drop under different particle sizes.Figure 8
Variation of energy dissipation with seepage velocity under different particle sizes.According to Darcy’s law, the cross-sectional areaA, the seepage flow Q through A, the head drop (h1 − h2), and the seepage path length L are in the following relationship [19]:(8)Q=K∙A∙h1−h2L,which is simplified as follows:(9)vs=K∙J,where K is the permeability coefficient, determined by the solid and fluid materials; vs is the seepage velocity; and J is the hydraulic gradient. The corresponding curves based on Darcy’s law are shown in Figure 8. With the increase in seepage velocity, the seepage phenomenon for each of the five types of particles gradually deviates from Darcy’s law. This is because the seepage process can be regarded as laminar fluid flow under a small seepage velocity. With a larger seepage velocity, the seepage velocity and hydraulic gradient are no longer linear, and the seepage process changes into a turbulent fluid flow. A similar deviation from Darcy’s law was also observed at high velocity [20, 21].


## 4.4. Effect of Packed Bed Particle Shape

Since the ore particles are generally nonspherical ones in practical industries, this work also examines the effect of particle shape on the seepage process, such as prolate particles and oblate particles. Note that each particle with different shapes (ellipsoid aspect ratio) has the same volume in this study. The aspect ratio of a spherical particle is 1.0, while it is 1.25 for a prolate particle and 0.82 for an oblate particle as shown in Figure9. Figure 9(a) represents spherical particles with a radius of 0.57 mm. The other ellipsoidal particles’ projections are perpendicular to the inlet direction presenting a circle with a radius of 0.61 mm and an ellipsoid with a long/short semiaxis of 0.75 and 0.5 mm, respectively.Figure 9
Packing of (a) spherical particles; (b) oblate particles, and (c) prolate particles.
(a)(b)(c)Figure10 shows the variation of energy dissipation with seepage velocity for the systems packed with differently shaped particles. It is shown that the energy dissipation is not significantly influenced by the particle shape under a small seepage velocity. However, the seepage resistance for a spherical particle system is smaller than that of a nonspherical particle system at a larger seepage velocity. For the nonspherical systems, the prolate particle system has larger seepage resistance and energy dissipation than that of the oblate particle system. The difference can also be observed in the contours of fluid pressure on the particle surfaces shown in Figure 11 and the corresponding reasons are discussed later in the article.Figure 10
Variation of energy dissipation with seepage velocity under different particle shapes.Figure 11
Contours of the pressure on varied particle surfaces (seepage velocity is 0.011 m/s).The mean curvature, referring to the ratio of the angle between the tangents at the ends of the arc to the length of the arc, and the frontal contact surface area of the particle are calculated for different particle shapes as shown in Table1. Since the selected arc corresponds to the long axis, the angle between the ends of the arc is 180°. The frontal contact surface area in the situation of this work equals half of the particle surface area. As listed in Table 1, the oblate and spherical particles have almost the same curvature. In contrast, the energy dissipation under larger seepage velocities for oblate particles is larger than that for spherical particles. This may reveal that the curvature does not determine the energy dissipation during leaching. The difference in frontal contact surface area corresponds better with the difference in energy dissipation, indicating that the frontal contact surface area may be the main factor determining the energy dissipation.Table 1
Curvature and frontal contact surface area for each particle shape.
Particle shapeMean curvatureFrontal contact surface area (mm2)Spherical particle0.89534.0828Oblate particle0.89554.1138Prolate particle1.03504.1888


## 5. Conclusions

This work applies a validated CFD model to fundamentally study the seepage features in the ionic rare earth ore leaching process. The results indicate that a smaller packing porosity and a smaller particle size of the granular porous medium lead to larger energy dissipation during seepage. The energy dissipation increases with seepage velocity. Due to turbulence, the energy dissipation becomes larger than that predicted by Darcy’s law when the seepage velocity increases to a large value. Additionally, the granular porous medium packed with prolate particles causes more energy dissipation, followed by the oblate and spherical particles, which may be mainly caused by the varied particle frontal contact surface areas.



---

*Source: 5916294-2023-02-11.xml*