# A Coupled Model for Solution Flow and Bioleaching Reaction Based on the Evolution of Heap Pore Structure

Authors: Shenghua Yin; Xun Chen; Chi Ma; Leiming Wang; Aixiang Wu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Based on the basic seepage law, equations have been derived to descript the solution flow within the copper ore heap which is treated as anisotropy porous media. The relationship between heap permeability and pore ratio has been revealed. Given the consideration of cover pressure and particle dissolution, pore evolution model has been set up. The pore evolution mechanism, due to the process of dissolution, precipitation, blockage, collapse, and caking, has been investigated. The comprehensive model for pore evolution and solution flow under the effect of solute transport and leaching reaction has been established. A trapezoidal heap was calculated, and the estimated results show that permeability decreases with the decreasing of pore ratio. Therefore, the permeability of the heap with small particles is relatively low because of its low pore ratio. Furthermore, permeability and height are found to be the two main factors influencing the solution flow.


---

## Body

## 1. Introduction

The classic percolation mechanism assumes that the porous media within which the solution flows are completely rigid, that is, the solid matrix develops no elastic or plastic deformation as the evolution of the fluid pore pressure; thus, the seepage can be treated as uncoupling issues. For heap leaching of low-grade copper ores, however, some unexpected variations are likely to occur throughout, such as mineral dissolution, solution precipitation, and microparticle migration. Unfortunately, on the one hand, these variations result in structure evolution, which leads to heap deterioration such as pore ratio and permeability. On the other hand, the mentioned variations in turn affect the pore fluid and pressure distribution [1]. Therefore, we must consider the seepage mechanism within the porous media and its impact on the inner structure of the heap, for example, the coupling effect between the seepage and the mass transfer as well as multiphase reaction.Normally, the interaction of the porous media and the seepage in leaching is complicated due to physical, chemical reactions and temperature effect [2, 3]. Therefore, the evolution of the structure mainly reflects on the pore ratio development within the heap, which in turn results in the state changes of the porous flow channel and flow velocity [4, 5]. As a result, instead of a constant, only variables of the permeability and osmotic coefficient can be obtained macroscopically [6].


## 2. Basic Seepage Law for Solution Flow within the Heap

### 2.1. Darcy’s Law

When kinetic inertia force is ignorable and the hydraulic gradient is greater than the initial gradient, the flow within the anisotropic porous media agrees well with linear resistance according to Darcy’s law(1)
v
x
=
-
k
x
μ
∂
p
∂
x
;
v
y
=
-
k
y
μ
∂
p
∂
y
;
v
z
=
-
k
z
μ
∂
p
∂
z
.Given the consideration of the gravity, the differential form of Darcy’s law becomes(2)
v
x
=
-
k
x
μ
∂
p
∂
x
;
v
y
=
-
k
y
μ
∂
p
∂
y
;
v
z
=
-
k
z
μ
∂
p
∂
z
+
ρ
w
g
,
where v
x,  v
y, and v
z are the flow velocity parallel to x,  y, and z axes, respectively, m/s; k
x,  k
y, and k
z are effective permeability parallel to x,  y, and z axes, respectively, m2; μ is the coefficient of flow viscosity, kg/(m·s); p is the transfusion pressure, Pa; ρ
w is the fluid density, kg/m3; g is the gravity acceleration, m/s, while (2) can be implied by a gradient form as
(3)
v
→
=
-
k
μ
∇
p
-
ρ
w
g
∇
D
.


### 2.2. Continuous Equation for Saturated Solution Flow

In a certain porous media region within the heap, assume that an infinitesimal with the length ofd
x, d
y, and d
z centers on point P
(
x
,
y
,
z
), and the infinitesimal volume is V
=
d
x
d
y
d
z, and then the molar difference of the inlet flow and outlet flow in unit time can be demonstrated as
(4)
Δ
Q
=
-
∂
ρ
w
v
x
∂
x
+
∂
ρ
w
v
y
∂
y
+
∂
ρ
w
v
z
∂
z
d
x
d
y
d
z
.If we substitute (1) into (4), then we draw
(5)
Δ
Q
=
ρ
w
∂
∂
x
k
x
μ
∂
p
∂
x
+
∂
∂
y
k
y
μ
∂
p
∂
y
+
∂
∂
z
k
z
μ
∂
p
∂
z
k
y
μ
∂
p
∂
y
d
x
d
y
d
z
.Moreover, if the pore ratio isϕ, the solution occupied volume in the infinitesimal therefore is n
d
x
d
y
d
z, and the corresponding mass M is n
ρ
d
x
d
y
d
z
=
n
ρ
V. Thus, variation rate of mass M according to time t is
(6)
∂
M
∂
t
=
∂
ϕ
ρ
w
V
∂
t
.On the basis of law of mass conservation, the accumulated rate of the solution mass is the same as the rate according to time. As stated by(7)
ρ
w
∂
∂
x
k
x
μ
∂
p
∂
x
+
∂
∂
y
k
y
μ
∂
p
∂
y
+
∂
∂
z
k
z
μ
∂
p
∂
z
d
x
d
y
d
z
=
ρ
w
V
∂
ϕ
∂
t
+
ϕ
ρ
∂
V
∂
t
+
ϕ
V
∂
ρ
w
∂
t
,
the left terms are indicative of the accumulated rate of the solution mass in the infinitesimal, therein the first term on the right is compression rate of the ore particle, the second term is the changing rate of the porous volume, and the third term is the changing rate of the fluid intensity. In particular, if we do not take account of the ore particle deformation and solution density variation, we obtain the basic differential equation of the heap flow from (7):
(8)
∂
∂
x
k
x
μ
∂
p
∂
x
+
∂
∂
y
k
y
μ
∂
p
∂
y
+
∂
∂
z
k
z
μ
∂
p
∂
z
=
∂
ϕ
∂
t
.


### 2.3. Solution Flow Model for Unsaturated Heap

Gas phase and liquid phase are hardly miscible due to an obvious gas-liquid interface and interfacial tension therein that separate the two phases from each other. Since the gas-liquid two-phase flow models involve comprehensive parameters (e.g., solution saturation, gas saturation, and relative permeability), we usually postulate that saturations of each phase are constants, and pore ratio as well as osmotic coefficients of each phase stays stable during leaching. Considering the gravity effect, we characterize the gas-liquid two-phase flow in the heap with gas-liquid dynamic partial differential equation:(9)
∂
∂
x
k
r
w
k
a
x
B
w
μ
w
∂
p
w
∂
x
+
∂
∂
y
k
r
w
k
a
y
B
w
μ
w
∂
p
w
∂
y
+
∂
∂
z
k
r
w
k
a
z
B
w
μ
w
∂
p
w
∂
z
+
ρ
w
g
=
∂
∂
t
ϕ
S
w
B
w
,
where k
r
w is relative permeability of the solution, B
w is volumetric coefficient of the solution, p
w is the solution pressure, Pa, μ
w is viscosity coefficient of the solution, kg/m·s, and k
a
x,  k
a
y, and k
a
z are absolute permeability along x,  y, and z axes, respectively, m2.


## 2.1. Darcy’s Law

When kinetic inertia force is ignorable and the hydraulic gradient is greater than the initial gradient, the flow within the anisotropic porous media agrees well with linear resistance according to Darcy’s law(1)
v
x
=
-
k
x
μ
∂
p
∂
x
;
v
y
=
-
k
y
μ
∂
p
∂
y
;
v
z
=
-
k
z
μ
∂
p
∂
z
.Given the consideration of the gravity, the differential form of Darcy’s law becomes(2)
v
x
=
-
k
x
μ
∂
p
∂
x
;
v
y
=
-
k
y
μ
∂
p
∂
y
;
v
z
=
-
k
z
μ
∂
p
∂
z
+
ρ
w
g
,
where v
x,  v
y, and v
z are the flow velocity parallel to x,  y, and z axes, respectively, m/s; k
x,  k
y, and k
z are effective permeability parallel to x,  y, and z axes, respectively, m2; μ is the coefficient of flow viscosity, kg/(m·s); p is the transfusion pressure, Pa; ρ
w is the fluid density, kg/m3; g is the gravity acceleration, m/s, while (2) can be implied by a gradient form as
(3)
v
→
=
-
k
μ
∇
p
-
ρ
w
g
∇
D
.


## 2.2. Continuous Equation for Saturated Solution Flow

In a certain porous media region within the heap, assume that an infinitesimal with the length ofd
x, d
y, and d
z centers on point P
(
x
,
y
,
z
), and the infinitesimal volume is V
=
d
x
d
y
d
z, and then the molar difference of the inlet flow and outlet flow in unit time can be demonstrated as
(4)
Δ
Q
=
-
∂
ρ
w
v
x
∂
x
+
∂
ρ
w
v
y
∂
y
+
∂
ρ
w
v
z
∂
z
d
x
d
y
d
z
.If we substitute (1) into (4), then we draw
(5)
Δ
Q
=
ρ
w
∂
∂
x
k
x
μ
∂
p
∂
x
+
∂
∂
y
k
y
μ
∂
p
∂
y
+
∂
∂
z
k
z
μ
∂
p
∂
z
k
y
μ
∂
p
∂
y
d
x
d
y
d
z
.Moreover, if the pore ratio isϕ, the solution occupied volume in the infinitesimal therefore is n
d
x
d
y
d
z, and the corresponding mass M is n
ρ
d
x
d
y
d
z
=
n
ρ
V. Thus, variation rate of mass M according to time t is
(6)
∂
M
∂
t
=
∂
ϕ
ρ
w
V
∂
t
.On the basis of law of mass conservation, the accumulated rate of the solution mass is the same as the rate according to time. As stated by(7)
ρ
w
∂
∂
x
k
x
μ
∂
p
∂
x
+
∂
∂
y
k
y
μ
∂
p
∂
y
+
∂
∂
z
k
z
μ
∂
p
∂
z
d
x
d
y
d
z
=
ρ
w
V
∂
ϕ
∂
t
+
ϕ
ρ
∂
V
∂
t
+
ϕ
V
∂
ρ
w
∂
t
,
the left terms are indicative of the accumulated rate of the solution mass in the infinitesimal, therein the first term on the right is compression rate of the ore particle, the second term is the changing rate of the porous volume, and the third term is the changing rate of the fluid intensity. In particular, if we do not take account of the ore particle deformation and solution density variation, we obtain the basic differential equation of the heap flow from (7):
(8)
∂
∂
x
k
x
μ
∂
p
∂
x
+
∂
∂
y
k
y
μ
∂
p
∂
y
+
∂
∂
z
k
z
μ
∂
p
∂
z
=
∂
ϕ
∂
t
.


## 2.3. Solution Flow Model for Unsaturated Heap

Gas phase and liquid phase are hardly miscible due to an obvious gas-liquid interface and interfacial tension therein that separate the two phases from each other. Since the gas-liquid two-phase flow models involve comprehensive parameters (e.g., solution saturation, gas saturation, and relative permeability), we usually postulate that saturations of each phase are constants, and pore ratio as well as osmotic coefficients of each phase stays stable during leaching. Considering the gravity effect, we characterize the gas-liquid two-phase flow in the heap with gas-liquid dynamic partial differential equation:(9)
∂
∂
x
k
r
w
k
a
x
B
w
μ
w
∂
p
w
∂
x
+
∂
∂
y
k
r
w
k
a
y
B
w
μ
w
∂
p
w
∂
y
+
∂
∂
z
k
r
w
k
a
z
B
w
μ
w
∂
p
w
∂
z
+
ρ
w
g
=
∂
∂
t
ϕ
S
w
B
w
,
where k
r
w is relative permeability of the solution, B
w is volumetric coefficient of the solution, p
w is the solution pressure, Pa, μ
w is viscosity coefficient of the solution, kg/m·s, and k
a
x,  k
a
y, and k
a
z are absolute permeability along x,  y, and z axes, respectively, m2.


## 3. Relationship between Heap Permeability and Pore Ratio

With the assumption that neither elastic/plastic deformation nor the solid structural distorts within the solid matrix occur, the classic percolation mechanism commonly overlooks the developments of the permeability and pore ratio. However, microstructure may comprehensively evolve due to the cover pressure, mineral dissolution, precipitation, and particle migration and then disturbs the seepage low [7, 8]. In order to set up the mathematical simulation of seepage process under the coupling effect of mass transfer and leaching reaction, it is important to build the relational model between permeability and pore ratio and thus establish the synthetical expression of the pore ratio evolution with mineral dissolution, precipitation, and cover pressure [9].It is hypothesized that heap deformation occurs due to cover pressure, mineral dissolution, and so forth. According to the hydromechanical continuous equation, we obtain the continuous equation of the solid matrix as follows:(10)
ρ
s
1
-
ϕ
∇
·
V
s
+
1
-
ϕ
∂
ρ
s
∂
t
-
ρ
s
∂
ϕ
∂
t
=
0
,
where ρ
s is the particle density, kg/m3, and V
s is the migration rate of the particle, m/s.However, from (10), we have
(11)
1
-
ϕ
∂
ε
V
∂
t
+
1
K
s
∂
p
w
∂
t
=
∂
ϕ
∂
t
,
where ε
V is the volumetric strain of the ore particle and K
s is the bulk modulus of the particle, MPa.Noting that the integral equation of (11), we have
(12)
1
-
ϕ
=
C
·
e
-
ε
V
+
Δ
p
w
/
K
s
.Ifε
V
=
0,  K
s
≥
1, then ϕ
=
ϕ
0, and we substitute this term into (12) and assume that C
=
1
-
ϕ
0, so we obtain
(13)
ϕ
=
ϕ
ε
V
,
p
=
1
-
1
-
ϕ
0
·
e
-
ε
V
-
Δ
p
/
K
s
,
where ϕ
0 is the initial pore ratio.Moreover, when we approximately expend the exponential term in (13), we get
(14)
ϕ
=
1
-
1
-
ϕ
0
1
-
Δ
p
w
/
K
s
1
+
ε
V
.Considering the temperature effect, (14) can also be expressed as
(15)
ϕ
=
1
-
1
-
ϕ
0
1
-
Δ
p
w
/
K
s
+
3
β
s
Δ
T
1
+
ε
V
.However, the porous permeability is not a constant in the same vein but evolves continuously due to many influential factors in leaching. According to Kozeny equation in percolation mechanism, we deduce the permeability equation as(16)
k
=
k
0
1
+
ε
V
1
+
ε
V
ϕ
0
-
3
β
s
Δ
T
+
Δ
p
w
/
K
s
1
-
ϕ
0
ϕ
0
3
.Therefore, we establish the relational model between permeability and pore ratio within the heap, as shown by(17)
k
=
k
0
1
+
ε
V
×
1
+
ε
V
-
1
-
ϕ
1
+
ε
V
-
1
-
ϕ
0
2
Δ
p
w
/
K
s
-
1
ϕ
0
3
.Equation (17) elucidates that permeability k mainly is influenced by pore ratio ϕ and covers pressure Δ
p; it increases as the ascending of ϕ while it decreases as the ascending of Δ
p. On the one hand, cover pressure Δ
p gradually increases with heap heightening and repeatedly mechanical lamination, which lowers the heap permeability. On the other hand, precipitation, micrograin blockage, and caking occur as leaching proceeds, which aggravates the permeability within the heap.


## 4. Pore Evolution under the Coupled Effect of Solute Transport and Leaching Reaction

As aforementioned, the cover pressure and temperature gradient influence the pore ratio in some extent. Additionally, the mineral dissolution, chemical precipitation, particle blockage, and caking also reflect on the gradual evolution of the heap [10]. Firstly, when we consider the migration equation during the seepage process according to law of mass conservation, we have
(18)
∂
∂
t
ϕ
C
1
=
-
∂
∂
x
i
v
C
1
+
∂
∂
x
i
D
∂
C
1
∂
x
j
-
q
1
ρ
w
-
R
1
,
∂
∂
t
ϕ
C
2
+
1
-
ϕ
ρ
s
ρ
w
G
=
-
∂
∂
x
i
v
C
2
+
∂
∂
x
i
D
∂
C
2
∂
x
j
-
q
2
ρ
w
+
R
2
,
where C
1 is the concentration of the leaching reagent, kg/L, q
1 is the inlet or outlet velocity of the leaching reagent, kg/m2·s, R
1 is the consumption rate of the leaching reagent, kg/s, C
2 is the metal ion concentration, kg/L, G is the grade of the ore, %, q
2 is the inlet or outlet velocity of the ore, kg/m2·s, and R
2 is the consumption rate of the ore, kg/s.Generally, leaching reactions occur between soluble reagentA and solid ore particle B and then generate soluble M and insoluble solid S, as shown by
(19)
n
A
A
solution
+
n
B
B
solid
=
n
M
M
solution
+
n
S
S
solid
.In this case, the dissolution rate of the solid ore is(20)
d
C
2
d
t
=
4
π
r
2
δ
ρ
s
G
λ
n
B
C
1
4
π
r
2
δ
ρ
w
ϕ
n
A
=
ρ
s
G
λ
n
B
ρ
w
ϕ
n
A
C
1
,
where r is the radial of the particle, m, δ is the diffusion thickness of the particle, m, λ is reaction rate constant, and n
A and n
B are stoichiometric coefficient.Assume that (19) is reversible; namely, the reaction proceeds toward both left and right. Thus, the dissolution rate of the solid ore can also be described as
(21)
d
C
2
d
t
=
λ
f
C
1
-
λ
b
C
2
,
where λ
f and λ
b are forward reaction rate and backward reaction rate, respectively.Additionally, when the reaction reaches a relative equilibrium state, (21) can be simplified as
(22)
λ
f
C
1
=
λ
b
C
2
,
in which equilibrium constant of the reaction is
(23)
λ
eq
=
λ
f
λ
b
=
C
2
C
1
.Similarly, if we postulate thatλ
f
=
(
n
B
/
n
A
)
λ
1 and R
2
=
d
C
2
/
d
t, we can draw (24) from (21) and (23):
(24)
R
2
=
n
B
n
A
λ
1
1
-
C
2
C
1
λ
eq
C
1
=
n
B
n
A
λ
1
1
-
Q
m
λ
eq
C
1
=
γ
C
1
,
where Q
m is the activity coefficient of the ions; simultaneously, the precipitation coefficient is defined as
(25)
γ
=
ρ
s
G
λ
ϕ
ρ
w
n
B
n
A
1
-
Q
m
λ
eq
.Therefore, the pore ratio of the heap is presented by(26)
ϕ
=
ρ
s
G
λ
γ
ρ
w
n
B
n
A
1
-
Q
m
λ
eq
.In summary, the pore ratio of the heap is mainly influenced by solution precipitation, leaching rate, and the equilibrium coefficient.


## 5. Solution Flow Model Given Consideration to Solute Transport and Leaching Reaction

In the wake of seepage process, concomitant phenomena such as heap compaction, mineral dissolution, precipitation, and micrograin migration always take place [11]. On the one hand, the unexpected process causes the variation of the microstructure and the pore ratio and thus exacerbates the permeability macroscopically. On the other hand, the variation results in uneven distribution of the seepage in the heap (e.g., overland flow, preferential flow, and channeling) and then hinders the target mineral from sufficient dissolution [12, 13]. In the view of this phenomenon, we established the seepage model within the porous media under the coupling effect of mass transfer and leaching reaction ((27), according to (9), (17), and (26)):
(27)
∂
∂
x
k
r
w
k
a
x
B
w
μ
w
∂
p
w
∂
x
+
∂
∂
y
k
r
w
k
a
y
B
w
μ
w
∂
p
w
∂
y
+
∂
∂
z
k
r
w
k
a
z
B
w
μ
w
∂
p
w
∂
z
+
ρ
w
g
=
∂
∂
t
ϕ
S
w
B
w
,
k
=
k
0
1
+
ε
V
×
1
+
ε
V
-
1
-
ϕ
1
+
ε
V
-
1
-
ϕ
0
2
Δ
p
w
/
K
s
-
1
ϕ
0
3
,
ϕ
=
ρ
s
G
λ
γ
ρ
w
n
B
n
A
1
-
Q
m
λ
eq
.


## 6. Calculated Example

The model was applied to calculate the certain parameters of a two-dimensional heap of copper ores. For the onsite heap, it is commonly long and narrowly designed and characterized of its section along the short axis. For simplicity, the special section was assumed to be isosceles trapezoid, the upper side length of which is 40 m, lower side length is 80 m, and height is 40 m (Figure1). Additionally, the two-dimensional model is consisted of porous ore particles, where the gas-liquid two-phase flow fluids intraparticles.Figure 1
A two-dimensional heap model for calculation.The model was solved by using the simulation tool COMSOL Multiphysics directly. An adaptive finite element mesh was used for the computations with a comparatively fine meshing along the boundaries and a comparatively coarse meshing in center of the heap. In this range of elements the computations were insensitive to the grid size. Following the established model, when the initial pore ratio is 0.3, the corresponding permeability is 5.36 m2, while when the cover pressure is 3.5 MPa, the corresponding volumetric stain of the ore is 8.275 × 10−7. The bulk modulus of the ore particle herein is 90 MPa, and the relative curve between the permeability and the pore ratio is shown in Figure 2. It implies that the permeability increases gradually as the pore ratio increases. The ascending process, however, can be divided into two stages. The first stage exerts a relative slow ascending and, especially, the permeability increases positively but slightly when the pore ratio is less than 0.5. On the contrary, the second stage witnesses a severe rapid ascending and, especially, the permeability increases positively and severely when the pore ratio is greater than 0.5. However, thein situ heap permeability is generally lower than 0.5, let alone the heaps with micrograins or/and high silt content; thus, the low permeability therein leads to a relatively worse permeability.Figure 2
The relative curve between the permeability and the porosity based on the two-dimensional model.In heap leaching, the physical and chemical reactions generally change the heap permeability by means of pore ratio variation and hence change the flow velocity within the heap. Given this fact, the sectional boundary of the model includes inlet boundary (the upper side of the heap), outlet boundary (the lower side of the heap), and left and right boundaries (along the heap slope).The mathematical simulation was carried out according to the following hypothesis.(1)
Since there is hardly any solution flows along the heap slope, the boundary of the ore particles is(28)
v
left
and
right
boundary
=
0
.
(2)
The solution is bearing the uniform pressure distribution on the inlet and outlet boundaries; moreover, we have(29)
p
upper
boundary
=
p
0
,
p
lower
boundary
=
0
.According to the model, the flow velocity in the center of the heap is demonstrated by Figure3. As shown, the flow velocity at the height of 40 m from the bottom up (namely, heap surface) is the same as the spray intensity (1.0 mm/s). The solution flow downwards and around, with its velocity, linearly decreases with heap heightening due to various resistance and, not unexpectedly, reaches the minimum speed when it arrives at the bottom. Meanwhile, the permeability affects the flow velocity comprehensively, where in the flow velocity at the same height increases as the ascending of the permeability on the condition of constant spray intensity at the heap surface. In addition, the curve slope increases as the permeability decreases, which indicates that the lower the permeability is, the faster the flow velocity lowers down. Particularly, the flow velocity at the heap bottom is as low as 0.25 mm/s when the permeability is 2 m2. Furthermore, the velocity profile of the heap with permeability k
=
6 m2 is shown in Figure 4.Figure 3
The relation curve of the flow velocity in the center of the heap and the heap height.Figure 4
The velocity profile of the dump with permeabilityk
=
6 m2.


## 7. Conclusions

(1)
The pore structure within the heap affects comprehensively the seepage during the leaching. The conventional research, however, normally treats the porous heap as rigid media, wherein the pore structure remains stable. The basic equations for the seepage mainly include Darcy’s law, continuous equation of the saturated solution, and the unsaturated seepage model.
(2)
A relational expression between the pore ration and permeability within the heap is established. Our present study indicates that the static permeability is mainly influenced by pore ratio and cover pressure, and it increases as the ascending of the pore ratio, whereas it decreases as the ascending of the cover pressure. Additionally, the pore structure within the heap witnesses a large extent evolution due to mineral dissolution, solution precipitation, and microparticle migration; thus, it disturbs the seepage law.
(3)
A calculation example was conducted following a two-dimensional model, the results of which indicate that the permeability decreases gradually as the descending of the pore ratio. To the heaps with micrograins or/and high silt content, its low permeability leads to a relatively worse permeability state. Besides, the flow velocity is greatly affected by the permeability and heap height.



---

*Source: 342537-2014-11-03.xml*