Universit¨at Stuttgart - Institut f¨ ur Wasserbau Lehrstuhl f¨ ur Hydromechanik und Hydrosystemmodellierung Prof. Dr.-Ing. Rainer Helmig

Master Thesis

Modelling of Multiphase Multicomponent Flow and Transport in Heap Leaching of Copper Ores

Submitted by

Jing Li Matrikelnummer 2144375

Stuttgart, 8th November, 2005

Examiners: Prof. Dr.-Ing. Rainer Helmig Supervisor: Dr.-Ing. Holger Class, Dipl.-Ing. Andreas Kopp, M.Sc.

Acknowledgement This study is supported by many people to whom I would like to express my deep gratidute. I would like to say a lot of thanks to Prof. Rainer Helmig who gave me the chance to study this interesting topic through which I obtained precious knowledge and experiences on modelling. I also benefited a lot from his encouragement and advice. I particularly thank my supervisors: Dr. Class and Andreas Kopp who were always patient to answer my questions and gave me a lot of warmhearted help to solve all the problems. I would like to thank also all my friends in Wasserbau (Anozie, Leo, Harmut, Ulrich, Iryna . . . ). They are always willing to help me when I have resort to them. My very special thanks goes to my beloved husband, Min Liu, who supported me morally and helped me from all the aspects.

Task of Jing Li’s Master Thesis Modelling of Multiphase Multicomponent Flow and Transport in Heap Leaching of Copper Ores The configuration considered in this master thesis consists of a ore heap that is irrigated on top with a dilute solution of sulfuric acid. As the solution permeates through the ore material it dissolves the minerals (the copper ions). In this way the copper is leached out of the ore and can be separated later on from the leaching solution and be further processed. Two phenomenas of interest in this work are the multiphase flow problem through the ore heap and the physico-chemical reactions. The reactions of the acid leaching solution with the copper will be simulated as a desorption problem from the ore into the solution. The mass transfer of the copper into the solution is in our case not dependend on the acid concentration. Furthermore the acid concentration does not deplete due to the leaching process. In detail, the following tasks have to be fullfilled: • Simulate flow through the heap by a 2 Phase flow model in 2 dimensions. Necessary steps include grid construction, definition of boundary conditions, define parameter settings like irrigation rate etc. • Simulate 2 phase 3 component flow through the heap by an extendet model setup. Components in the liquid phase include copper, sulfuric acid and air. • A sensitivity study has to be conducted. The parameters under investigation can be the Brooks-Corey or vanGenuchten constants in the capillary pressure saturation relationship, irrigation rate, absolute permeability, relative permeability, boundary conditions, initial settings, desorption coverning parameters or any other related process or parameter with the agreemen of the supervisor. • Alternatively or additionally the simulation can be transfered to the 3 dimensional space.

I

Author’s Statement I herewith certify that I have prepared this master thesis independently, and that only those sources, aids, and advisors that are duly noted herein have been used and / or consulted.

Signature

Date

Contents 1 Introduction 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Concepts and Steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 1 2 2

2 Model Concepts 2.1 Two-phase (2p) Model . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Mathematical Formulations . . . . . . . . . . . . . . . 2.1.2 Supplementary constraints - Constitutive relationships 2.2 Two-phase Three-component (2p3c) Model . . . . . . . . . . 2.2.1 Advective/Convective Transport (Viscous Flow) . . . . 2.2.2 Diffusive Transport . . . . . . . . . . . . . . . . . . . . 2.2.3 Mathematical Formulations . . . . . . . . . . . . . . . 2.2.4 Supplementary Constraints . . . . . . . . . . . . . . . 2.3 Numerical Scheme . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Spatial Discretization - Fully Upwind Box Method . . . 2.3.2 Time Discretization - Implicit Euler Scheme . . . . . . 2.3.3 Mass Lumping . . . . . . . . . . . . . . . . . . . . . . 2.3.4 Linearization - Newton-Raphson Method . . . . . . . . 2.3.5 MUFTE-UG . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . .

4 4 5 7 9 11 11 12 14 19 19 21 22 23 24

3 Simulations 3.1 Definition of the Domain . . . . . . . . . . 3.1.1 Initial and Boundary Conditions . . 3.1.2 Initial and Boundary Conditions for 3.1.3 Initial and Boundary Conditions for 3.2 Simulation Results . . . . . . . . . . . . . 3.2.1 Evolvement of water saturation . . 3.2.2 Evolvement of copper concentration

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

26 26 26 28 30 32 32 34

4 Sensitivity Test 4.1 Influence of the Irrigation Ratio . . . . . . . . . . . . . . . . . . . . . . 4.2 Influence of the Intrinsic Permeability . . . . . . . . . . . . . . . . . . .

37 37 41

III

. . . . . . . . . . . . . . . . . . the 2p Model . the 2p3c Model . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . .

. . . . . . .

CONTENTS

4.3 4.4

4.5 4.6

Influence of the Porosity . . . . . . . . . . . . . Influence of the Pc − Sw Curve . . . . . . . . . . 4.4.1 Influence of the residual water saturation 4.4.2 Influence of the Van Genuchten α . . . . 4.4.3 Influence of the Van Genuchten n . . . . 4.4.4 Analysis and Interpretation . . . . . . . Influence of the Adsorption Coefficient . . . . . Conclusions . . . . . . . . . . . . . . . . . . . .

IV

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

45 48 48 51 53 54 55 57

5 Simulations under Heterogeneous Conditions

58

6 Summary and Outlook

62

List of Figures 1.1

Configuration of Heap Leaching . . . . . . . . . . . . . . . . . . . . . .

2.1 2.2 2.3 2.4 2.5 2.6 2.7

Concept of REV. . . . . . . . . . VG pc -Sw and krα -Sw relationship Phase transitions . . . . . . . . . Henry’s law and Raoult’s law . . Box method . . . . . . . . . . . . Newton-Raphson method . . . . . MUFTE-UG . . . . . . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

5 8 10 17 21 24 25

3.1 3.2 3.3 3.4 3.5 3.6 3.7

Boundary conditions of the 2p model. . . . . . . . . Boundary conditions of the 2p3c model. . . . . . . Color contour of water saturation evolvement. . . . Approximate location of the sample points . . . . . Sw − t curve. . . . . . . . . . . . . . . . . . . . . . Colour contour of copper concentration evolvement. Xwc − t curve. . . . . . . . . . . . . . . . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

29 31 32 33 34 35 36

4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15

Approximate location of the sample points . Influence of irrigation ratio on Sw . . . . . . Influence of irrigation ratio on Xwc . . . . . . Influence of absolute permeability on Sw . . Influence of absolute permeability on Xwc . . Influence of porosity on Sw . . . . . . . . . Influence of porosity on Xwc . . . . . . . . . Influence of Swr on Sw . . . . . . . . . . . Influence of Swr on Xwc . . . . . . . . . . . Influence of Van Genuchten α on Sw . . . . Influence of Van Genuchten α on Xwc . . . . Influence of Van Genuchten n on Sw . . . . Influence of Van Genuchten n on Xwc . . . . pc − Sw curves for sensitivity test. . . . . . . Influence of adsorption coefficient on Sw and

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

38 39 41 43 44 46 47 49 50 51 52 53 54 55 56

V

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Xwc

. . . . . . . . . . . . . . .

3

LIST OF FIGURES

5.1 5.2 5.3 5.4 5.5

Set-up of the heterogeneous domain. . . . . . . Evolvement of Sw in the heterogeneous domain Sw − t curves in the heterogeneous domain . . . Evolvement of Xwc in the heterogeneous domain Xwc − t curves in the heterogeneous domain . . .

VI

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

59 59 60 61 61

List of Tables 2.1 2.2

Parameters of pc curve implemented in the simulation . . . . . . . . . . List of secondary parameters . . . . . . . . . . . . . . . . . . . . . . . .

7 14

3.1 3.2 3.3

Values of common parameters used in the models . . . . . . . . . . . . Initial and boundary conditions for the 2p model . . . . . . . . . . . . Initial and boundary conditions for 2p3c model . . . . . . . . . . . . .

27 29 31

5.1

Parameters of the heterogeneous domain . . . . . . . . . . . . . . . . .

58

VII

Nomenclature Symbol A c ccw Dκα Dκpm,α Fsc g Hwa Ja Jd K Kef f Kf krα M Mc Mw n Nj (x) pα pc pg pg0 pκg pw sat qα qc qg

Unit [-] [kg/m3 ] [kg/m3 ]

Definition coefficient matrix volumetic concentration volumetic concentration of copper in water phase 2 [m /s] molecular diffusion coefficient of component κ in phase α 2 [m /s] macroscopic diffusion constant of component κ in phase α in a porous media [kg/kg] copper content in the solid matrix [m/s2 ] gravity vector (0, −g)T [Pa ] Henry coeffient of component air in water phase 2 [kg/(s · m )] advective flux 2 [kg/(s · m )] diffusive flux [m2 ] absolute permeability [m/s] effective permeability [m/s] hydraulic conductivity [-] relative permeabiltiy of phase α [-] mass matrix [kg/mol] molecular mass of copper [kg/mol] molecular mass of water [-] Van Genuchten parameter [-] linear ansatz function at node j [P a] pressure of phase α [P a] capillary pressure [P a] gas phase pressure [P a] initial gas phase pressure [P a] partial pressure of component κ in gas phase [P a] saturation vapor pressure [s−1 ] sink and source term for phase α 3 [kg or mol/m · s] Neumann flux of copper [kg or mol/m3 · s] Neumann flux of gas

Nomenclature

Symbol qw ra rc rw Sg Sw u vα W Xακ x, y α α %α %mol,α µα λα Ω Γ τ

IX

Unit [kg or mol/m3 · s] [mol/m3 · s] [mol/m3 · s] [mol/m3 · s] [-] [-] [-] [m/s] [-] [mol/mol] [m] [-] [P a−1 ] [kg/m3 ] [mol/m3 ] [kg/(m · s)] [m · s/kg] [-] [-] [N/m2 ]

Definition Neumann flux of water sink and source term for component air sink and source term for component copper sink and source term for component water saturation of gaseous phase saturation of water phase primary variable Darcy velocity of phase α weighting function mole fraction of component κ in phase α coordinates phase indication Van Genuchten parameter mass density of phase α molar denstiy of phase α dynamic viscosity of phase α mobility of phase α indication of the whole domain boundaries of the domain sheer stress

Chapter 1 Introduction 1.1

Background

Copper plays an important role in modern society. It is quite close to our life, and appears in coins, kitchen utensils, furniture, paints and musical instruments. It is also present in cars, planes, ships and even space vehicles, either in the form of pure element or as part of alloys. Many manufacturers of electronic equipment use copper because it is more efficient in conducting electricity and it lasts longer than other materials. [14] Nowadays, most of the available copper is sparsely distributed over large areas, mixed with mineralized materials and rock. These are the porphyric deposits, which could only be exploited with metallurgical skills separating and recovering the metal. [14] “Copper is mainly found associated to sulfide minerals, but also to oxide minerals. These two types of minerals require different productive processes, but in both cases the starting point is the same: the extraction of the material from open pit or underground mines, which requires the fragmentation and transportation of the material that has been previously identified by geological surveys. The extracted mineral goes first through a milling process. In the case of oxide minerals, the processing involves submitting the material to a leaching solution, which will produce solutions of copper sulphate, which enter a process of solvent extraction followed by electrowinning, which final result is a copper cathode 99.99% pure. The sulfide minerals go first through crushing and milling, followed by a classification process to obtain the copper concentrate, with 30% copper. Later purification stages are carried out in furnaces that generate blister or anode copper 99Finally, electro-refining transforms the anodes into 99.99% pure cathodes.” [14]

1.2 Motivation

1.2

2

Motivation

Currently, heap leaching is considered as a superior alternative to other conventional processing of the ore in copper mining industry. It is a cheaper and more effective alternative to e.g. flotation, agitation and vat leaching [7]. The effectiveness could be further improved by numerical simulation of the processes. In the framework of this thesis, simulations with numerical simulator MUFTE-UG (Multiphase Flow, Transport and Energy, Unstructured Grids) are carried out to estimate the copper recovery rate and provide a preliminary step for industrial application of the model, e.g. to optimize the operation conditions through sensitivity studies.

1.3

Concepts and Steps

The configuration of the heap is a trapezoid consists of fragmented ore material. The heap can be considered as a homogeneous domain which is irrigated on top with a diluted sulfuric acid solution which desolves and also desorps the copper when it permeates through. In this way the copper is leached out of the heap and collected at the bottom for a later separation. During the irrigation process, the flow of the two fluids (liquid and gas) and transport of copper are triggered (see Figure 1.1). These are the two phenomena of interest of this thesis and the multiphase flow is simulated by both the two-phase and two-phase three-component models, while the transport process is simulated by the two-phase three-component model. As mentioned, the extraction of copper from the solid matrix includes two processes: reaction and adsorption. I must state now that only the adsorption-desorption process is considered in our model. That is to say we don’t distinguish the difference between water and the acid solution. We assume that water has the same desorption ability, therefore the component acid is simply replaced by the component water in the two-phase three-component model. This is the most important simplification made in the thesis. As a first step, in order to have a clear description of the flow behavior of the heap leaching process, a two-phase flow model is employed. The two phases refer to the water phase (leach solution is simplified by water) and gaseous phase. At this stage, the change of water saturation over time within the domain is studied. Since the concentration change of copper is of our greatest interest, a two-phase three-component model is built up to simulate both the transport processes of copper and flow behavior of the phases. As mentioned before, since it seems too ambitious to model both the chemical reaction and adsorption process within the frame work of one thesis, we only consider the adsorption-desorption of copper between the solid matrix and the liquid phase and drop out the dissolution process of copper from ore bed into the liquid phase by sulfuric acid. Therefore the component of sulfuric acid can be

1.3 Concepts and Steps

3

Irrigation Process

Air

5m

H2SO4

Air H2O

Cu

45o

45o 25 m Figure 1.1: Configuration of Heap Leaching

neglected, and a two-phase (water and gas) three-component (water, copper and air) system can be established. In this step, both the change of copper concentration and the change of water saturation are investigated. The latter one is compared with the result of the two-phase model. Then, in order to test the performance of the models under different conditions, a sensitivity study is carried out by varying the parameters of interest, for instance, the irrigation rate, absolute permeability, adsorption-desorption coefficient and so on. And initially, it is intended to compare the results of our models with those of the models which are developed by a Chilean researcher Mr. E. Cariaga as a verification step. But his models still have problem with convergence, therefore this step can not be done. Subsequently, since a layered configuration is sometime used in the mining industries in order to enhance the efficiency of leaching process, our models are further extended to simulate this process with a heterogeneous domain which is divided into three horizontal layers. Finally, an outlook is given about the future research work to improve the model and adapt it to industrial applications.

Chapter 2 Model Concepts In this chapter the mechanisms of Multiphase multicomponent flow and transport considered in this thesis are presented. For both the two phase model and the two-phase three-component model, a basic description along with the simplifications involved and the mathematical and numerical methods are discussed. The descriptions of the underlying mathematical formulation and numerical methods of simulation are taken from Helmig (1997) [11].

2.1

Two-phase (2p) Model

Maybe at first it is necessary to mention here that the terminology phase “is used to differentiate between two or more fluid continua, separated by a sharp interface, across which discontinuities in fluid properties exist [11].” Multiphase model reflects the reality with the concept of REV (Representative Elementary Volume), within which the exact distribution of inidividual phases seperated by sharp interfaces is simplified by phase saturation (see Figure 2.1). Then we get macroscopic parameters such as saturation, relative permeability, capillary pressure and so on which contribute together to the description of flow behavior the phases. All the balance equations which will be introduced afterwords take the REV as the control volume. The heap leaching process can be abstracted to a water-gas two-phase flow model. The flow type of the two phases in the porous medium (ore bed) is unsaturated flow, in which the flow permeability depends on the degree of saturation and capillary pressure [16]. During the irrigation process, water phase displaces the gaseous phase and fills up the pores gradually. The two fluids have different mobilities and therefore distint velocities, but they interact on each other through capillary pressure-saturation relationship and relative permeability-saturation relationship [16]. Some simplifications are made in the 2p model:

2.1 Two-phase (2p) Model

5

microscale

REV soil matrix liquid phase

averaging gas phase

Figure 2.1: Concept of REV.

- Since the temperature change is negligible, isothermal condition is assumed. - Only the gaseous phase is considered compressible and the ideal gas law is assumed to be valid. - The compressibility of the porous medium and thus a pressure dependent change of the porosity is supposed to be negligible. - Effects of hysteresis are not considered, which means gravity is always in equilibrium with capillarity.

2.1.1

Mathematical Formulations

After the conceptual model is set up, we need the mathematical formulations to describe the physical processes. The basic mathematical formulations in the 2p model include the conservation of mass, momentum and energy (energy balance is neglected here, due to the assumption of isothermal condition), as well as the system dependent equations of state (constitutive relationships). The derivation of the formulations are based on the continuum theory which corresponds to an averaging process of substance properties. The properties at the micro scale are averaged within the REV (see Figure 2.1). In our case, a REV consists of the solid matrix and two fluid phases: water and gaseous. The composition of the REV is characterized by the porosity φ [−] of the solid matrix and saturation Sα [−] of the phases. They are defined as follows: φ=

pore volume total volume of the porous medium

2.1 Two-phase (2p) Model

Sα =

6

volume of fluid phaseα pore volume

α ∈ {w, g}

in which w indicates the water phase and g indicates the gaseous phase. The transition from micro scale to macro scale produces new equations, e.g. Darcy equation which is derived from the basic Navier-Stokes-equation at micro scale, and macroscopic parameters such as saturations and relative permeabilities of different fluids. Moreover, the integration of the boundary and initial conditions for the primary variables is also very important for the mathematical model (for details, see Chapter 3). 2.1.1.1

Mass Balance

First a mass balance for each phase has to be set up: ∂(φSα %α ) | ∂t {z }

+

∇ · (%α vα ) | {z }

−

advection term

storage term

% α qα |{z}

=0

(2.1)

sink/source term

where α ∈ {w, g}, vα stands for the velocity of each phase and qα for any sink and source terms. In our problem we don’t have any sink or source. 2.1.1.2

Momentum Balance

The conservation of momentum is based on Darcy’s law for single phase flow. However a large number of experiments related to multiphase processes has shown that the Darcy velocity for each phase in a porous medium can be described by the generalized form of Darcy’s law [11] [18]: vα = −

krα K · (grad pα − %α g) µα

(2.2)

where krα is the relative permeability ranging from 0 to 1 for the respective phase. 2.1.1.3

pg -Sw Formulation

Since the available initial and boundary conditions are given for gas phase pressure pg and water phase saturation Sw , they are chosen as the primary variables and the pg Sw formulation is selected. At first, the following modifications (based on the closure equations) are required: ∇pw = ∇(pn − pc )

(2.3)

2.1 Two-phase (2p) Model

7

∂Sg ∂ ∂Sw = (1 − Sw ) = − (2.4) ∂t ∂t ∂t Inserting equation (2.2) and modification of equation (2.3) and (2.4) into the mass balance equation (2.1), one obtains for the water phase: ∂(φSw %w ) − ∇ · (λw %w K(∇pg − ∇pc − %w g)) − %w qw = 0 ∂t

(2.5)

the gaseous phase: −

∂(φSw %g ) − ∇ · (λg %g K(∇pg − %g g)) − %g qg = 0 ∂t

where λw = kµrw and λg = w phase respectively.

2.1.2

krg . µg

(2.6)

They are the mobilities of liquid phase and gaseous

Supplementary constraints - Constitutive relationships

In order to solve the basic mathematical equations, supplementary constraints or socalled constitutive relationships have to be set up to describe the dependence of secondary variables on the primary variables. 2.1.2.1

Capillary Pressure

On the macro scale, the most usual correlations of pc − Sw for a two-phase system can be described by Brooks and Corey (see Brook and Corey (1964) [3]) or by Van Genuchten [9]. The latter one is used in this thesis: Se (pc ) =

Sw − Swr = [1 + (α · pc )n ]m 1 − Swr

for pc > 0

where Se [-] is the effective water saturation, Swr [-] is the residual water saturation, n [-], m [-] and α [1/Pa] are VG paramters. All these paramters are based on parameters characterizing the pore geometry and determined by fitting to experimental data [11]. The following table shows the value of each parameter used in our simulations: parameter value

Swr [-] 0.0

Sgr [-] α [1/Pa] 0.0 1.35 · 10−4

n [-] 1.411

Table 2.1: Parameters of pc curve implemented in the simulation The left picture of Figure 2.2 shows the capillary pressure against water saturation using the values in Table 2.1.

2.1 Two-phase (2p) Model

8

Figure 2.2: VG pc -Sw and krα -Sw relationship 2.1.2.2

Relative Permeability

In case of two or more fluids flowing simultaneously through a porous medium, a relative permeability krα [-] for each of the fluids can be defined as the ratio of effective permeability Kef f [m/s] to intrinsic permeability K [m2 ]. The intrinsic permeability of the porous medium only depends on the properties of the solid matrix, whereas the effective permeability depends on both the structure of the porous medium and the respective saturations. Furthermore, the hydraulic conductivity Kf accounts for the properties of both solid matrix and corresponding fluid, as well as the saturation of this fluid. It is defined as: Kf = Kkrα

%α g %α g = Kef f µα µα

Then the Darcy velocity of a certain fluid phase is calculated as: vα = −Kf ∇( vα = −K

pα + z) %α g

krα (∇pα − %α g) µα

Together with the capillary pressure, the relative permeability characterizes the interaction between liquid phase and solid phase. The relative permeability and saturation relationship can also be approximated by functions after Van Genuchten and Brooks and Corey. The former one is chosen since the Van Genuchten relation is applied for the description of capillary pressure-saturation relationship. krw

 m i 2 p h 1 m = Se 1 − 1 − S e

(2.7)

2.2 Two-phase Three-component (2p3c) Model

krn

9



i p  13 h 1 2m m = 1 − Se 1 − Se

(2.8)

The right picture of Figure 2.2 shows the relative permeability-saturation relationship according to the values in Table 2.1. 2.1.2.3

Density

As stated in Section 2.1, only the gaseous phase is considered compressible. Therefore the density of the water phase is a constant in the 2p model whose value is given in chapter 3. Another assumption made in Section 2.1 is the validation of the ideal gas law, from which the density of the gaseous phase can be derived: pg = with

n R Ma Ra

nMa Ra T nRT = , V V

(2.9)

number of moles [mol] universal gas constant [J/(mol K)] molar mass of air [kg/mol] individual gas constant of air [J/(kg K)]

from equation 2.9, one can obtain the gas density %g [kg/m3 ]: %g = 2.1.2.4

pg Ra T

(2.10)

Viscosity

The dynamic viscosity µα [Pa s] describes the ratio of the shear stress τ [N/m2 ] and velocity gradient in a fluid phase α [13]: τ =µ

du dy

The pressure dependence of the viscosities of both phases is neglected in the 2p model. Although, normally, their temperature dependence should be considered, they have constant values in the 2p model due to the assumption of isothermal condition. The values of the viscosities of both phases are give in Chapter 3.

2.2

Two-phase Three-component (2p3c) Model

If we are interested not only in the flow behavior of the fluids, but also in the transport processes of a certain component, we need the multiphase and multicomponent model. In our case, the concentration change of copper in the water is of the greatest interest, therefore a two-phase three-component model (2p3c model) is set

2.2 Two-phase Three-component (2p3c) Model

10

up. The term component “stands for constituents of the phases which can be associated with a unique chemical species (e.g., hydrogen, oxygen) or a molecular substance (e.g., pure water), as well as a mixture of different substances (e.g., air) [11] [6].” The two phases in the 2p3c model still refer to the water and gaseous phases. The gaseous phase consists of two components: air (a) and water (w), whereas the water phase consists of three components: air, water, and copper (c). As mentioned in Chapter 1, the component sulfuric acid is replaced by component water, because the chemical reaction between sulfuric acid and copper or the dissolution of copper by sulfuric acid is not considered in our model. Only the adsorption-desorption process is modeled, and no difference is distinguished between the sulfuric acid and water. In addition, the component air is actually a pseudo-component with averaged fluid properties without accounting for the fractions of different constituents (N2 , O2 , etc.) [11]. A component can be transferred from one phase into another by the change of the thermodynamic state (e.g., by vaporization, condensation, dissolution, etc.). This phenomenon is called phase transition, Figure 2.3 illustrates the occurring phase transitions in the 2p3c model. evaporation condensation

Liquid Phase

Gaseous Phase water (w)

water (w) air (a)

degassing

air (a)

copper (c) dissolution so

de n

n

tio

rp

tio

so

rp

ad

Solid Phase

Figure 2.3: Phase transitions

Besides the assumptions made in the 2p model, three additional assumptions are added. Due to the low flow velocity, local equilibrium of phase transitions are assumed. In other

2.2 Two-phase Three-component (2p3c) Model

11

words, the equilibrium of phase transitions (e.g. adsorption-desorption, degassingdissolution) is built up much faster than the flow velocities of the fluids, then phase transitions are at equilibrium at every point within the domain. Due to the same reason, dispersion is neglected and only diffusion is considered. Moreover, as aforementioned, chemical reaction is not considered, only adsorption-desorption is modelled by the linear sorption isotherm.

2.2.1

Advective/Convective Transport (Viscous Flow)

The advective or convective transport of a component is driven by the viscous flow which is caused by gravity or a pressure gradient. If a fluid phase consists of different components, the moving fluid transports all of them at the same velocity, which is given by the generalized Darcy law. According to Cirpka (1997) [5] The advective mass flux density Ja of a component is: Ja = vc in which v [m/s] is the Darcy velocity (see equation (2.2)) of the fluid, and c [g/l3 ] is the volumetric concentration of the respective component.

2.2.2

Diffusive Transport

“Diffusion is a mass-transfer process caused by the random Brownian motion of solute particles in fluids [5].” The concentration gradient of a component leads to a diffusive flux Jd of the respective component from locations of higher to lower concentrations [11]. It can be described by Fick’s first law : Jd = −Dκpm,α ∇c

(2.11)

in which Dκpm,α [m2 /s] is the macroscopic diffusion constant of component κ within phase α. It is derived from the molecular diffusion coefficient Dκα [m2 /s]. In the porous media, the molecular diffusion coefficient is smaller than in a free fluid because the diffusive path lines are tortuous [5]. In a multiphase system within a porous media, the diffusion process is further hindered by the unsaturated condition of the fluid. Thus Dκpm,α is dependent on the porosity, tortuosity [-] of the porous medium and saturation Sα of the respective fluid: Dκpm,α = φτ Sα Dκα (2.12) The tortuosity itself is a function of the saturation and porosity, it increases when the saturation and porosity increase. In this two-phase three-component model, the model of Millington and Quirk (1961) is applied [11] [15]: 1

7

τ = φ 3 Sα3

(2.13)

2.2 Two-phase Three-component (2p3c) Model

12

insert equations (2.13) and (2.12) into equation (2.11), the diffusive flux is: 4

10

Jd = −φ 3 Sα3 Dκα ∇c

(2.14)

The molecular diffusion coefficient of each component in different phases are given by empirical value which is specified in Chapter 3.

2.2.3

Mathematical Formulations

The basis of the mathematical description of the two-phase three-component model is also the mass and momentum balance. Instead of considering the flow and transport equations separately, we will write one balance equation for each component. Therefore the mass transfer of the phase transitions of the components from one phase into another by e.g. dissolution, vaporization and condensation (see Figure 2.3) is automatically included in the formulation. That is to say the sink and source terms of components in each phase are embedded in the formulation.[11] In the two-phase three-component model, the mole fraction is used to describe the compositions of the phases. xκα [-] denotes the mole fraction of component κ in phase α. It is the ratio between the number of moles of a component and the total number of moles of all the components in the corresponding phase. for component air (a): ∂ [φ(%mol,w xaw Sw + %mol,g xag Sg )] ∂t +div(%mol,w xaw vw + %mol,g xag vg ) −div[Dapm,w %mol,w gradxaw

+

Dapm,g %mol,g gradxag ] −ra

storage advective transport diffusive transport sink and sources

=0

(2.15)

for component water (w): ∂ w [φ(%mol,w xw w Sw + %mol,g xg Sg )] ∂t w +div(%mol,w xw w vw + %mol,g xg vg ) w −div[Dw pm,w %mol,w gradxw

+

w Dw pm,g %mol,g gradxg ]

−rw =0

storage advective transport diffusive transport sink and sources (2.16)

2.2 Two-phase Three-component (2p3c) Model

13

for component copper (c): ∂ [φ%mol,w xcw Sw + (1 − φ)%s %mol,w kd xcw ] ∂t +div(%mol,w xcw vw ) −div[Dcpm,w %mol,w gradxcw ] −rc

storage advective transport diffusive transport sink and sources

=0

(2.17)

in which: %mol,w /%mol,g : a c xw w /xw /xw : a xw g /xg : w Dpm,w /Dapm,w /Dcpm,w : a Dw pm,g /Dpm,g :

%s kd

molar density of water/gaseous phase [mol/m3 ] mole fraction of water/air/copper in the water phase [-] mole fraction of water/air in the gaseous phase [-] macroscopic diffusion constant of water/air/copper in the liquid phase [m2 /s] macroscopic diffusion constant of water/air in the gaseous phase [m2 /s] bulk density of the porous medium [kg/m3 ] adsorption-desorption coefficient [m3 /kg]

vw and vg are the Darcy velocities of water and gaseous phases: krw K(∇pw − %w g) µw krg K(∇pg − %g g) = − µg

vw = − vg

(2.18)

in which %w and %g are the mass densities of water and gaseous phases in [kg/m3 ]. In equation (2.17), the part of storage of copper in the solid phase is derived as follows: (1 − φ)%s Fsc Mc

is the copper stored in the porous medium in mol/m3

(2.19)

of copper in which, Mc is the molecular weight of copper, and Fsc = kg kg which is the of solid matrix mass fraction of copper. Its value is given in our case, and the copper concentration in the water phase is calculated from it via the linear partitioning law or so called linear sorption isotherm: Fsc = kd ccw (2.20)

in which kd [m3 /kg] stands for the patition coefficient between the solid and the aqueous phase, namely the adsorption coefficient, and ccw [kg/m3 ] is the volumetric concentration of copper in the water phase whose relationship with the mole fraction of copper

2.2 Two-phase Three-component (2p3c) Model

14

in water phase is: xcw =

ccw /Mc ccw = %w /Mw Mc %mol,w c ⇒ cw = xcw Mc %mol,w

(2.21)

in which, %w [kg/m3 ] and Mw [kg/mol] is the mass density and molecular weight of the water phase. Inserting equation (2.21) and (2.20) into (2.19), the term in formulation (2.17) is obtained: (1 − φ)%s %mol,w kd xcw

2.2.4

Supplementary Constraints

According to the given initial and boundary conditions, we choose the gas pressure pg , water saturation Sw and mole fraction of copper in the water phase xcw as our primary variables. Again we need to set up all the constitutive relationships to solve the basic equations. Besides those used in the 2p model, additional constraints are added for the 2p3c model. First, the closure equations read as follows:

Sw + S g = 1

a c xw w + xw + xw = 1

a xw g + xg = 1

p w = pg − pc

Second, we need to formulate equations of state which describe the coefficients of the transport equations as functions of the primary variables. The secondary parameters which depend on the primary variables are listed in Table 2.2. parameter saturation relative permeability viscosity mass density capillary pressure water mass fraction air mass fraction water molecular diffusivity air molecular diffusivity copper molecular diffusivity Henry constant

liquid phase gaseous phase Sg (Sw ) kr w(Sw ) kr g(Sw ) µw (T ) µg (xw c ,T) c %g (pg , xw , T ) %w (T ) pc (Sw ) c xw (p , x , T ) xw w g w g (pg , T ) a c xw (pg , xw , T ) xag (pg , T ) Dw pm,g (Sw ) a Dpm,g (Sw ) c Dpm,w (Sw ) a Hgw (T )

Table 2.2: List of secondary parameters

2.2 Two-phase Three-component (2p3c) Model

15

We should notice that the temperature T is assumed to be constant in the model. In addition, switching of the primary variables due to phase change may influence the dependence of the secondary variables on specific primary variables. But no phase change occurs in our model, which means we always have two phases. Additional constitutive relationships are discussed in details in the following pages. 2.2.4.1

Capillary Pressure and Relative Permeability

The capillary pressure-saturation and relative permeability-saturation relationships for the 2p3c model remain the same as for the 2p model. 2.2.4.2

Density

As in the 2p model, the density of the gaseous phase is variable and that of the water phase is assumed to be a constant. The former one is influenced by the constitution of the gaseous phase. Its mass density is calculated from the molar density: a %g = %mol,g (xw g Mw + x g Ma )

for the gaseous phase

(2.22)

in which, Mw , Ma are the molecular weight of water and air respectively. And the calculations of molar density again conforms to the ideal gas law: pg RT in which R is the universal gas constant with the value of 8.314 [J/(molK)]. %mol,g =

(2.23)

From the above equations, it is obvious that the mass density of gaseous phase changes with the mole fractions of components which constitute the gas phase. But its influence on the flow behavior of the fluids should be insignificant, since the percentage of water vapor in the gaseous phase is very small. Then since the mass density of the water is given, the molar density can be obtained easily from: %w %mol,w = (2.24) Mw 2.2.4.3

Viscosity

As aforementioned, viscosity of fluid is normally considered as a function of temperature. With increasing temperature, the dynamic viscosity of liquid phase decreases whereas the dynamic viscosity of gaseous phase increases. But this effect has no significance in our model because of the assumption of isothermal condition.

2.2 Two-phase Three-component (2p3c) Model

16

Normally, in a multiphase multicomponent model, the dynamic viscosity of the gaseous phase µg [P a · s] is further dependent on the phase constitution. It is a combination of viscosity of the component air µa and that of the component water vapor µwvap . According to Haefner et al. (1985) [8], µg of gas mixtures is: µg = xag µa + xw g µwvap

(2.25)

One needs to find out a way to get the dynamic viscosity of air and water vapor. For example, the rule after Reid, Prausnitz and Poling (1987) [17] can be used to calculate the dynamic viscosity of water vapor. The influence of phase composition on viscosity is much smaller for the water phase and is often neglected. However, the values of dynamic viscosities of both phases in our models are set to fixed values, which are given in Chapter 3, according to E.Cariaga. 2.2.4.4

Phase Transitions between Water and Gas

As illustrated in Figure 2.3, the contents of a component in different phases are correlated via the corresponding phase transition processes. The relationship obeys a certain partitioning law of compounds between different phases. Because the component copper does not participate in the phase transition between water and gaseous phases and its concentration is very low, here, we first consider the liquid phase as a binary mixture of air and water system. For this kind of liquid mixture of two compounds that differ strongly in their chemical structure, at the same time one component is in very small concentration (xaw → 0) and the other one is in very high concentration (xw x → 1), we call it an ideally diluted solution. Under this condition, Henry’s law is valid for the solute, in our case the component air, and Raoult’s law is valid for the solvent which in our case is the component water. Figure 2.4 shows that Henry’s law and Raoult’s law do good approximations at the two ends of extreme mole fractions: If the water saturation is less than the residual water saturation, the partitioning of water between liquid and gaseous phase is calculated by Raoult’s law : w w pw g = xw psat

(2.26)

w in which pw g is the partial pressure of water vapor in gas phase and psat is the saturation vapor pressure of water which is a function of temperature.

2.2 Two-phase Three-component (2p3c) Model

17

Figure 2.4: Partial pressures of the two components in equilibrium within a real binary mixture. [6]

To calculate the mole fraction of water and air in the gaseous phase, Dalton’s law is applied: X X pκg κ κ κ pg = pg xg = 1 xg = (2.27) pg κ κ Since water alway exists as a liquid phase in our model, the saturation vapor pressure is in equilibrium state equal to the partial pressure of water vapor in the gas phase: w pw g = psat

therefore xw g =

pw g pg

(2.28)

The partitioning of air between liquid and gaseous phase conforms to Henry’s law : xaw

pag = a Hw

pag

= pg −

pw g

xag

pag = pg

(2.29)

in which pag is the partial pressure of air in the gas phase and Hwa [Pa ] is the Henry coefficient for the solution of the component air in the liquid phase. The mole fraction of copper in the liquid phase xcw is considered as one of the primary c a variables, from equation (2.28), (2.29) and the restriction of xw w + xw + xw = 1, one can set up the whole interrelations of the mole fractions of each component in each phase.

2.2.4.5

Phase Transition between Solid and Water - Adsorption

The copper dissolved in the water phase is originated from the solid matrix or the heap. As stated before, only adsorption is considered as the extraction process

2.2 Two-phase Three-component (2p3c) Model

18

of copper in our model. The flowing water carries the desorbed copper away and dilutes its concentration in the original place, which breaks the sorption equilibrium instantaneously. In order to build up the new equilibrium again, more copper is desorbed from the solid phase immediately. In this way the copper in the solid can be depleted gradually. The adsorption-desorption process is described by the linear adsorption isotherm (see equation (2.20)). From the balance equation for component copper (Eq. (2.17)), the storage term can be reformulated to:   (1 − φ)%s kd ∂(%mol,w Sw xcw ) φ+ (2.30) Sw ∂t in which the term

(1 − φ)%s kd (2.31) Sw is the so called retardation factor R. This retardation factor, R, is always larger than one and shows that sorption at equilibrium leads to an enhancing factor in the time derivative of the aqueous mole fraction compared to the case without sorption. This means that the advective-diffusion terms lead to a slower copper concentration change over time in the sorbing system than in the non-sorbing one, because more mass is stored than just the mass in the liquid phase. Therefore the term ”retardation” means ”slowing down”. [5] φ+

This effect is examined in the sensitivity test in Chapter 4. 2.2.4.6

Diffusion

From section 2.2.2, one can conclude that the diffusivity of a component in a certain phase is dependent on the porosity, tortuosity of the porous medium and is a function of the saturation of the respective phase. The effective diffusion coefficient of a component in a certain phase can be calculated from the molecular diffusion coefficient: 4

10

Dκpm,α = −φ 3 Sα3 Dκα

(2.32)

Furthermore, the diffusion fluxes of the components in the two phases have the following relationship: Jww = −Jwa − Jwc

Jga = −Jgw

(2.33)

which means that the amount of the water (Jww ) diffused out is equal to the amount of the air (Jwa ) and copper (Jwc ) diffused into the water phase, and the amount of air (Jga ) diffused out equals to the amount of the water vapor (Jgw ) diffused into the gaseous phase.

2.3 Numerical Scheme

2.3

19

Numerical Scheme

The balance equations described in Chapter 2 represent a system of coupled and nonlinear differential equations and therefore can not be solved without the use of numerical methods. Numerical methods include the discretization schemes which replace the differential terms by algebraic terms, the linearization of the algebraic system of equations, and the solution of the obtained linear system. In this thesis, the finite difference method is used for the time discretization, the fully upwind box method for the spatial discretization and the Newton-Raphson method for the linearization of the equation system. As aforementioned, the software program MUFTE-UG is used to implement the numerical schemes. They are introduced one by one in the following sections.

2.3.1

Spatial Discretization - Fully Upwind Box Method

The Fully Upwind Box Method is a Finite Volume formulation with piecewise linear shape functions including fully upwinding of the upstream mobilities.The box scheme is based on the Weighted Residual Method. In the spatial discretization, the exact value of the primary variables in the PDE (partial differential equation) is approximated by interpolation functions (in our case, the linear ansatz function N (x)) between the discrete nodal values:  1 if j = k Nj x = δjk = 0 if j 6= k Then:

u b nno

u e=

nno X j=1

u bj Nj

(2.34)

discrete nodal value of the primary variable number of nodes within the domain Ω

Inserting the approximation of the primary variables into the balance equation produces a non-zero value ε (which is called residual or defect) on the right hand side of the PDE. This residual ε is minimized within the domain by using a weighting function (test function), which is done by the weighted residual method: Z ! Wi εdΩ = 0 with i = 1, 2, 3, ..., nno (2.35) Ω

For example, applying this to the weak form (which is obtained by integrating the PDE

2.3 Numerical Scheme

20

over a domain Ω) of the pg -Sw formulation for the liquid phase: Z X ∂ Wi φ (ρw Sbw )j Nj dΩ ∂t Ω j∈ηi   Z X − Wi ∇ · (λw ρw )ij K (b pgj − pcj + ρw ij gzj )∇Nj dΩ − ηi (·)ij

Z

Ω

j∈ηi

!

Wi ρw qw dΩ = 0

(2.36)

Ω

set of neighboring nodes of node i denotes that one needs to decide where the influence coefficient (·) is evaluated

Applying the Green-Gaussian Theorem to the flux term of the balance equation: Z = =

Z Z



Wi ∇ · (λw ρw )ij K Ω

j∈ηi



∇ · Wi (λw ρw )ij K Ω

Wi (λw ρw )ij K Γ

X X j∈ηi

X



Ψwj ∇Nj dΩ 

Ψwj ∇Nj dΩ −

Ψwj ∇Nj ndΓ −

j∈ηi

Z

Z



∇Wi (λw ρw )ij K Ω



∇Wi (λw ρw )ij K Ω

X j∈ηi

X j∈ηi



Ψwj ∇Nj dΩ 

Ψwj ∇Nj dΩ (2.37)

where Ψw is the total potential of the liquid phase: Ψw = pw + ρw gz In the box method, a control volume Bi with corners at the midpoints of element edges and at the element’s centroids - as illustrated in Figure 2.5 is assigned to each node i and a piecewise constant weighting functions Wi is chosen [11]:  1 in the box Bi Wi = 0 outside of the box Bi R Since Wi is a constant,R ∇Wi = 0. Thus, the volume integral Ω disappears and only the boundary integral Γ remains. equation (2.36) can then be transformed to: Z X ∂ Wi φ (ρw Sbw )j Nj dΩ ∂t Ω j∈ηi Z X − Wi (λw ρw )ij K Ψwj ∇Nj ndΓ −

Z

Γ

j∈ηi

!

Wi ρw qw dΩ = 0 Ω

(2.38)

2.3 Numerical Scheme

21

Due to the hyperbolic characteristic of the advective flux-term/advective transport term, a fully upwind weighting is used to calculate the mobility term:  λαi if (Ψαj − Ψαi ) ≥ 0 F UB λαij = λαj if (Ψαj − Ψαi ) ≤ 0 This assures a correct description of the physical problem. control volume subdomain, Bi














































































































































 









































































































































i









  









































































 








center of gravity of element, e

finite element mesh j∈η

element e

Figure 2.5: Box method [1]

2.3.2

Time Discretization - Implicit Euler Scheme

Because the instationary balance equations are parabolic in time, i.e. the solution at a specific time step only depends on the state of the previous time step, a solution method which steps forward in time can be applied [11]. Considering the ordinary differential equation of the form: ∂u = f (u) ∂t

(2.39)

Applying a finite difference method yields: un+1 − un = f (um ) := f m ∆t

(2.40)

when m = n first order explicit method when m = n + 1 first order implicit method (implicit Euler) introducing a weighting parameter Θ: un+1 = un + ∆t(Θf n+1 + (1 − Θ)f n )

(2.41)

2.3 Numerical Scheme

22

Θ=0 first order explicit Euler Θ=1 first order implicit Euler
Θ = 0.5 second order Crank-Nicholson O (∇t)2

For the sake of stability aspect, the implicit Euler method is used in this model. We first rewrite equation (2.38) to:   db u + A[b uj ] = R, M dt db u = M−1 (R − Ab u) (2.42) dt in which M=

Z

Wi Nj dΩ

A=−

Ω

Z

v · gradNj dΓ Γ

The right hand side of the second equation corresponds to f in equation (2.40). Then the implicit Euler formulation can be derived: un+1 − un n+1 = [M−1 (R − Ab u)] ∆t

(2.43)

And R is the vector including the Neumann boundary flux terms and source/sink terms inside the domain Ω [11].

2.3.3

Mass Lumping

”Studies by Neumann et al. (1974) [19] for nonlinear problems have shown that the application of a lumped mass matrix has considerable computational advantages [11].” From equation (2.36), all coefficients of a row in the mass matrix Z Mij := Wi Nj dΩ Ω

are added up (lumped) to the diagonal entry of the mass matrix. Then the diagonalized lumped mass matrix coefficients can be expressed by: R  R Wi dΩ = Ω Ni dΩ =: Vi if i = j Ω Mij = (2.44) 0 if i 6= j In addition, the source/sink terms can be lumped in the same way: Wi ρw qw dΩ =: (ρw qw )i Vi

(2.45)

Finally, after all the above transformations (spacial and time discretization, GreenGaussian Theorem, upwinding method and mass lumping), one obtains the discretized

2.3 Numerical Scheme

23

weak form of the balance equation for the water phase:   φVi (n+1) n (ρw Sw )i − (ρw Sw )i ∆t Z X n+1 − (λw ρw )n+1 K ∇Nj ndΓ(Ψn+1 wj − Ψ wi ) ups(i,j) Γi

j∈ηi

−(ρw qw )n+1 Vi i

=0

(2.46)

The discretized balance equation of the gaseous phase can be derived in the same way.

2.3.4

Linearization - Newton-Raphson Method

The nonlinearity of the equation system is mainly caused by the constitutive relationships that depend on the saturation, e.g. the capillary pressure-saturation relation pc (Sw ) and the relative permeability-saturation relation krα (Sα ). Therefore a linearization of the algebraic system of equations is necessary and is done by the Newton-Raphson method [10] in this model. The coupled system of nonlinear PDEs resulted from the mathematical description of the multiphase/multiphasemulticomponent problems in the previous sections can be written in simplified functional form: F(x) = 0 (2.47) where x represents the the vector that holds the primary variables, for which equation (2.47) has to be solved. If we neglect the higher-order terms of the Taylor series, the Newton algorithm yields for the iteration step m + 1 at time level k + 1:   ∂F k+1,m+1 k+1,m F(x ) ≈ F(x )+ · (xk+1,m+1 − xk+1,m ) (2.48) ∂x k+1,m As F(xk+1,m+1 ) must become zero, we transform the above equation into: K(xk+1,m ) · u = −F(xk+1,m )

(2.49)

where K = ∂F/∂x stands for the Jacobian matrix, u = xk+1,m+1 − xk+1,m is the vector holding the corrections of the primary variables and F(xk+1,m ) represents the defect term at time level k + 1 and iteration step m. Analytical computation of the derivatives for the coefficient Kij of the Jacobian matrix is rather complex and therefore often computed numerically by: Kij =

∂Fik+1,m

∂xk+1,m j Fi (..., xj−1 , xj + ∆xj , xj+1 , ...) − Fi (..., xj−1 , xj − ∆xj , xj+1 , ...) (2.50) ≈ 2∆xj

2.3 Numerical Scheme

24

F(u)

∂F ∂u

F(un )

F(un+1)

exact solution, u un+1

un

u

Figure 2.6: Newton-Raphson method [1]

where ∆xj = δ · xj and δ is an arbitrary small increment, e.g. δ = 10−8 . Figure 2.6 illustrates the iteration concept of the Newton-Raphson linearization method.

2.3.5

MUFTE-UG

MUFTE-UG (MUltiphase Flow Transport and Energy model on Unstructured Grids) is a numerical simulator which is used for the simulation of the multiphase multicomponent flow and transport in porous medium. The part of MUFTE comprises the discretization schemes used for the multiphase-flow partial differential equations and a set of constitutive relationships. Currently, the Fully Upwind Box Method and the Control-Volume Finite-Element Method are the implemented discretization techniques. The former one is a Finite Volume formulation with a piecewise constant weighting function including fully upwinding of the upstream mobilities, whereas the latter one is a mass-conservative formulation on a discrete patch including a first order upwinding scheme (see Helmig:1997) [11]. The model concepts include the extended Darcy law for the multiphase flow in

2.3 Numerical Scheme

25

isotropic or anisotropic, homogeneous or heterogeneous media. Furthermore, a large set of constitutive relationships for the capillary pressure-saturation relationship and relative permeability-saturation relationship as well as state equations for densities, viscosities, etc. are available. (see Helmig et al (1994) [12]) UG is a toolbox for the solution of partial differential equations on unstructured grids. It provides a collection of ”state-of-the-art” numerical methods for a fast and efficient solution of system of equations. It comprises robust multi-grid solvers for both linear and non-linear problems as well as adaptive local grid refinement and coarsening techniques for a desired accuracy. It is able to handle complex geometries in two and three dimensions because of the data structures designed for unstructured grids. In addition, its functional parallelization enables possibility to solve complex problems efficiently. Finally, UG possesses several pre- and post-processing tools, which are applicable for parallel computers, e.g. a powerful graphical user interface. (see Bastian et al (1997) [2]). Both MUFTE and UG parts have a strictly modular structure which allows a simple combination of different modules special for a practical problem. Figure 2.7 shows the concept of MUFTE-UG.

(Helmig et. al 1997, 1998)

(S. Lang, K. Birken, K. Johannsen et. al 1997)

(Bastian et. al 1997, 1998)

Institute for Hydraulic Engineering (IWS)

- problem description - constitutive relationships - physical-mathematical models - discretization methods - numerical schemes - refinement criteria - physical interpretation MUFTE (Helmig)

Interdisciplinary Center for Scientific Computing (IWR)

- multigrid data structures - local grid refinement - solvers (multigrid, etc) - r,h,p-adaptive methods - parallelization - user interface - graphic representation UG (Wittum, Bastian)

Figure 2.7: MUFTE-UG

Chapter 3 Simulations The simulation of the copper leaching process is under the two dimensional condition and divided into two steps. The first step is to simulate the leaching process with the 2p model in order to give an insight on the flow behavior of both phases by investigation of water saturation evolvement. The second step is to simulate both the flow behavior of the fluids and the transport processes of copper by implementing the concept of the 2p3c model. In this step both the saturation evolvement and concentration evolvement are investigated, and the former one is compared with that obtained from the 2p model. The references of all the parameters used in the following simulations are taken from the Chilean researchers C. Ortiz A. and E. Cariaga. The common parameters used in the 2p model and the 2p3c model are listed in Table 3.1.

3.1 3.1.1

Definition of the Domain Initial and Boundary Conditions

After setting up the mathematical model and choosing the appropriate numerical methods, one needs to specify the initial and boundary conditions of the primary variables in order to find a unique solution of the equation system of PDEs and run the simulation. Initial conditions must be given for the whole domain. For the balance equations of the 2p model, it means: pg0 = f (x, y, t0 )

Sw0 = g(x, y, t0 )

x, y ∈ Ω

(3.1)

For the balance equations of the 2p3c model, aside from the above conditions, one more constraint must be added due to one more existing degree of freedom which requires

3.1 Definition of the Domain

name porosity mass density mass density absolute permeability liquid viscosity gas viscosity initial water saturation residual water saturation temperature heap slope heap width heap height

27

symbol φ %w %g K µw µg Sw0 Swr T θ W H

value 0.33 1011 pg /RT 1.78 · 10−12 10−3 1.85 · 10−5 0.4343 0.0 298.15 π/4 25 5

unit [-] [kg/m3 ] [kg/m3 ] [m2 ] [kg/(m · s)] [kg/(m · s)] [-] [-] [K o ] radians [m] [m]

Table 3.1: Values of common parameters used in both two-phase and two-phase twocomponent model

one more primary variable: xcw0 = h(x, y, t0 )

x, y ∈ Ω

(3.2)

The boundary conditions must be specified for the whole simulation time period, because the regarded process is transient. There are three types of boundary conditions which are implemented in the our models. 3.1.1.1

Dirichlet Boundary Condition

For this boundary condition, the fixed values of the primary variables must be given at the boundary. However, the Dirichlet conditions can also be time dependent [10]. In addition, Dirichlet conditions overwrite the initial conditions at the respective boundaries, if the latter one is different as the former one [10]. For the 2p model, the following conditions must be specified at a Dirichlet boundary: pgΓ1 = f (x, y, t)

SwΓ1 = g(x, y, t)

x, y ∈ Γ1

(3.3)

Again, one additional constraint must be added for the 2p3c model: xcwΓ1 = h(x, y, t) 3.1.1.2

x, y ∈ Γ1

(3.4)

Neumann Boundary Condition

For this type of boundary condition, the values of the normal derivative of the solution function should be given. It can be time dependent. Assignment of Neumann boundary

3.1 Definition of the Domain

28

conditions means the specifications of the mass fluxes across the respective boundary for the 2p model, and the molar fluxes for the 2p3c model.  kg/(m3 · s) for the 2p model qg = f (x, y, t) in (3.5) mol/(m3 · s) for the 2p3c model qw = g(x, y, t) in



kg/(m3 · s) for the 2p model mol/(m3 · s) for the 2p3c model

qc = h(x, y, t) in mol/(m3 · s) for the 2p3c model

x, y ∈ Ω x, y ∈ Ω

(3.6) (3.7)

in which, qg , qw and qc stand for the gas flux, water flux and copper flux, respectively. And all the fluxes correspond to the source/sink terms in the balance equations [10]. 3.1.1.3

Free Flow Boundary Condition

At this type of boundary, the fluids break through into the ambient environment at atmospheric pressure. The conditions (in our case, the targeted variables such as water saturation and copper concentration) at this boundary change over time during the flow process. This boundary cannot be described by either a Dirichlet or a Neumann boundary condition. In order to solve this problem, we extended the domain from this boundary and made a quasi-simulation of the ambient conditions. Virtual properties must be assigned to the extended part to ensure that the real boundary is influenced minimally by the arbitrary Dirichlet boundary conditions at the dummy boundary. [6] The details about the assignment of initial and boundary conditions are discussed in the following for the 2p model and the 2p3c model separately.

3.1.2

Initial and Boundary Conditions for the 2p Model

The whole domain of the heap is under atmospheric pressure and has the water saturation of 0.4343 at the very beginning. The western and eastern boundary are assumed as Neumann no-flux boundaries for both the water and gaseous phases. At the bottom of the heap, the fluid phases are in contact with the atmosphere, therefore a free flow boundary is assigned to the bottom. The domain is extended in the vertical direction by a factor of two. In the extended domain, the capillary pressure is set to zero in order to avoid the influence of capillary rise on the flow behavior, and the absolute permeability is enlarged one hundred times to allow the fluids to flow freely as if in the ambient environment. The assignment of initial and boundary conditions are listed in Table 3.2.

3.1 Definition of the Domain

x, y ∈ Ω

initial conditions

x, y ∈ Γi

Neumann boundary Dirichlet boundary

x, y ∈ Γl

Neumann no flux

x, y ∈ Γr

Neumann no flux

x, y ∈ Γd

Dirichlet boundary Neumann boundary

29

Sw0 (x, y, t0 ) = 0.4343 pg0 (x, y, t0 ) = 101325 qw (x, y, t) = 0.003 pg (x, y, t) = 101325 qw (x, y, t) = 0.0 qg (x, y, t) = 0.0 qw (x, y, t) = 0.0 qg (x, y, t) = 0.0 Sw (x, y, t) = 0.9 qg (x, y, t) = 0.0

[−] [Pa ] [kg/(m2 · s)] [pa ] [kg/(m2 · s)] [kg/(m2 · s)] [−] [kg/(m2 · s)]

Table 3.2: Initial and boundary conditions for the two-phase model

Figure 3.1 illustrates the boundary conditions and settings of properties in the extended domain.

Figure 3.1: Boundary conditions of the 2p model.

3.1 Definition of the Domain

3.1.3

30

Initial and Boundary Conditions for the 2p3c Model

All the initial and boundary conditions for the gas phase pressure and water phase saturation in the 2p3c model are the same as in the 2p model. The initial value of copper mole fraction in the water phase xcw is calculated from the given initial copper content F cs in the solid phase under the assumption of sorption equilibrium. It is given as follows: F cs = 2.53 · 10−6

kg mass of copper kg mass of solid phase

(3.8)

with the adsorption coefficient kd = 3.5 · 10−5 [m3 /kg]. According to the linear sorption isotherm: F cs =

kd ccw

⇒

ccw

2.53 · 10−6 F cs = = 0.0723 [kg/m3 ] = −5 kd 3.5 · 10

(3.9)

Subsequently, from the relation between Xwc and ccw derived in equation (2.21): xcw

0.0723 ccw = = 2.0113 · 10−5 [mol/mol] = Mc · %mol,w 0.064 · 56166.7

(3.10)

The top and both western and eastern sides of the heap are all Neumann no flux boundaries for the variable xcw , whereas at the bottom the free flow boundary condition is valid. Likewise, the domain is extended twice as high as the real heap to simulate the free flow boundary condition. The same settings for the absolute permeability and capillary pressure are used for the 2p3c model as for the 2p model. In addition, much higher diffusion coefficients are specified in the upper part of the extended domain, then the diffusion for each component is switched off at a lower level which is 4cm away from the bottom dummy boundary. By doing this, we ensure that the arbitrary Dirichlet value of xcw has tiny influence on the actual condition of copper concentration at the real bottom of the heap and allow the components to diffuse freely as if they are in contact with the atmosphere. The adsorption coefficient is set to zero in the extended domain. The assignment of the initial and boundary conditions for the 2p3c model are listed in Table 3.3 and illustrated in Figure 3.3.

3.1 Definition of the Domain

x, y ∈ Ω

initial conditions

x, y ∈ Γi

Neumann boundary Dirichlet boundary Neumann no flux

x, y ∈ Γl

Neumann no flux

x, y ∈ Γr

Neumann no flux

x, y ∈ Γ

d

Dirichlet boundary Neumann boundary Dirichlet boundary

31

Sw0 (x, y, t0 ) = 0.4343 pg0 (x, y, t0 ) = 101325 xcw0 (x, y, t0 ) = 2.0113 · 10−5 qw (x, y, t) = 0.003 pg (x, y, t) = 101325 qc (x, y, t) = 0.0 qw (x, y, t) = 0.0 qg (x, y, t) = 0.0 qc (x, y, t) = 0.0 qw (x, y, t) = 0.0 qg (x, y, t) = 0.0 qc (x, y, t) = 0.0 Sw (x, y, t) = 0.9 qg (x, y, t) = 0.0 xcw = 2.0 · 10−5

[−] [Pa ] [mol/mol] [kg/(m2 · s)] [pa ] [kg/(m2 · s)] [kg/(m2 · s)]

[kg/(m2 · s)] [−] [kg/(m2 · s)] [mol/mol]

Table 3.3: Initial and boundary conditions for the two-phase three-component model

Figure 3.2: Boundary conditions of the 2p3c model.

3.2 Simulation Results

3.2 3.2.1

32

Simulation Results Evolvement of water saturation

With all the aforementioned boundary and initial conditions, the simulation is executed The evolvement of water saturation is shown in Figure 3.3. The subfigures on the left hand side show the result from the 2p model, whereas those on the right hand side show the result from the 2p3c model. Due different simulation runs, the program is not able to produce the plot files with exactly the same time step for the two models, we can only compare the results of the two models with slightly different time steps. This will not be explained again. t = 8.95 hr

t = 8.58 hr

t = 23.86 hr

t = 23.49 hr

t = 47.72 hr

t = 47.35 hr

t = 77.55 hr

t = 77.18 hr

t = 110.36 hr

t = 109.98 hr

t = 170.01 hr

t = 171.13 hr

Sw: 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85

Figure 3.3: Color contour of water saturation evolvement. As shown in Figure 3.3, due to the irrigation process, the heap is wetted gradually and reaches its steady state with the water saturation of about 0.95 in the main part of the heap after 170 hours. Because of the capillary force, at the bottom of the heap,

3.2 Simulation Results

33

the front smears into the two corners. Comparing the results from the two models, no significant differences are displayed for the early time steps, only when the front propagates to the bottom of the heap, some differences between the shape of the fronts can be detected. Except the calculation of gas density, exactly the same values of other parameters are used for both models. As stated in section 2.2.4, the mole fraction of the water vapor in the gas phase is so small that its influence on the fluid properties can be neglected. In addition, diffusion can not account for this difference either, since the diffusion coefficient of each component in the water phase is very small. Therefore those aspects can not explain the difference. But it might be explained from the point of numerical accuracies. Because in the 2p model, we have two primary variables, while in the 2p3c model, we have three primary variable, that might lead to different numerical accuracies. Investigation of the change of water saturation over time are also done in three sample points. Their approximate locations are illustrated in Figure 3.4. The coordinates of these points from top to bottom are (12.5, 4.7), (12.5, 2.5) and (12.3, 0.3), respectively.

Figure 3.4: Approximate location of the sample points

The corresponding curves of water saturation against time are shown in Figure 3.5. The dash-dot lines show the result of the 2p model, while the solid lines show the result of the 2p3c model. The little difference between the two models can be observed from the curves.

3.2 Simulation Results

34

0.95 0.9 0.85 0.8

from the 2p model from the 2p3c model

Sw [-]

0.75 0.7

0.65 0.6 0.55 0.5 0.45 0.4

0

20

40

60

80

100

t [hr]

(a) top point

0.95

0.95

0.9

0.9

0.85

0.85

Sw [-]

0.8 0.75

Sw [-]

0.8 0.75 0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

80

100

0.4

0

20

40

t [hr]

(b) middle point

60

80

100

t [hr]

(c) bottom point

Figure 3.5: Change of saturation over time at the three sample points.

3.2.2

Evolvement of copper concentration

Figure 3.6 shows the color contour of the copper concentration evolvement obtained from the 2p3c model. Figure 3.7 shows the change of copper mole fraction in the water phase over time at the three sample points. We can see that the copper concentration in the water phase is always decreasing. That is because we assume the adsorption equilibrium from the initial time step on. The solid phase contains most copper at the very beginning, and thus copper concentration in the water reaches its maximum value at the initial time step. Afterwards it is diluted and washed out by the irrigated water contininuously.

3.2 Simulation Results

Xcw:

35

t = 0.0 hr

t = 109.98 hr

t = 8.58 hr

t = 141.30 hr

t = 47.35 hr

t = 171.13 hr

t = 77.18 hr

t = 224.8

1E-06 3E-06 5E-06 7E-06 9E-06 1.1E-05 1.3E-05 1.5E-05 1.7E-05 1.9E-05

Figure 3.6: Colour contour of copper concentration evolvement. In addition, comparing Figure 3.5 and Figure 3.7, we find out that the change of copper concentration is much slower than that of the water saturation. It is because copper is first stored in the solid phase and extracted through desorption process step by step after the already desorbed copper is washed away.

3.2 Simulation Results

36

2E-05

Xcw [mol/mol]

1.5E-05

1E-05

5E-06

0

0

50

100

150

200

250

300

t [hr]

2E-05

2E-05

1.5E-05

1.5E-05

Xcw [mol/mol]

Xcw [mol/mol]

(a) top point

1E-05

1E-05

5E-06

0

5E-06

0

50

100

150

200

t [hr]

(b) middle point

250

300

0

0

50

100

150

200

250

300

t [hr]

(c) bottom point

Figure 3.7: Change of copper mole fraction in the water phase over time at the three sample points.

Chapter 4 Sensitivity Test In order to test the ability and performace of the model to simulate the leaching process under different conditions, and also to provide a basis for optimization of operation conditions, some sensitivity tests are carried out by changing the parameters of importance one by one. In this chapter the results of different values of each parameter are compared by the curves of water saturation versus time and mole fraction of copper in the water phase versus time at the three sample points (see Figure 3.4). Simultaneously, the effects of the same sensitivity test on the 2p model and the 2p3c model are compared. In the following sections, all the figures about water saturation evolvement have two columns of pictures. The left column of three pictures shows the results of the 2p model, while the right one shows the results of the 2p3c model. This will not be explained again. In addition, the parameters under investigation take the original values for the solid lines we used in the simulations in Chapter 3. They are increased and decrease by a certain factor for the dashed and dash-dot lines respectively.

4.1

Influence of the Irrigation Ratio

Figure 4.2 shows the curves of water saturation over time under different irrigation ratios. The irrigation ratio for the solid lines is 10.68 l/(hr · m2 ). It is doubled for the dashed lines and halved for the dash-dot lines. It is obvious that the influence of reducing the irrigation ratio is larger than that of increasing it. With the middle irrigation ratio, the heap is nearly fully saturated. When the irrigation ratio is doubled, the maximum saturation is only slightly increased. If the irrigation ratio is further enlarged, the water saturation at the top of the heap always increases with time instead of reaching a horizontal plateau, and

4.1 Influence of the Irrigation Ratio

38

after a certain period of simulation run, the value of saturation exceeds 1. That is because under high irrigation rate, certain amount of water will pond on the top before it can entirely permeate into the heap. In the reality, the ponding water will flow down along the two sides of the heap. But our model has no function to simulate this phenomenon, therefore we get unreasonable values for the primary variables. Increase of the irrigation ratio enlarges the maximum saturation, because more water has to be developed within the same time, when the permeability remains unchanged, the additional water cannot permeate through in time and thus increases the saturation. A larger irrigation ratio also leads to a shorter lower plateau of the curve which reflects the time required for the front to propagate to this point. That means the flow velocity is increased. In addition, the slope of the curve is also steeper, which means a larger time derivative of saturation. From the mathematical point of view, irrigation corresponds to the Nuemann boundary flux. It goes into the flux term of the balance equation, which influences both the maximum saturation and flow velocity. A larger irrigation ratio not only increases the flux or flow velocity, but also the gradient of the flux (see Figure 4.1), thus the storage is increased as well. This can also be w inferred from equation 2.5. In the storage term, the part ∂S is reflected by the ∂t slope of the curve. That is why we see a steeper slope under the a higher irrigation ratio.

Fin storage

dx

Fin −Fout = storage dx

Fout

Figure 4.1: Approximate location of the sample points

In Figure 4.1, Fin is increased by the Neumann boundary flux, and Fout is initially produced by the drainage process. In addition, the effects of varying irrigation ratio on the results of the 2p model and of the 2p3c model show no big differences.

4.1 Influence of the Irrigation Ratio

39

0.95

0.95

0.9

0.9

0.85

0.85

0.8

0.8

0.75

0.75

Sw [-]

Sw [-]

(a). top point

0.7 0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

r=10.68 l/(hr. m2) r=5.34 l/(hr. m2) r=21.36 l/(hr. m2)

0

25

50

75

100

125

150

0.4

175

0

20

40

60

80

t [hr]

100

120

140

160

180

t [hr]

(b). middle point 0.95

0.95

0.9

0.9

0.85

0.85

0.75

Sw [-]

0.8

0.75

Sw [-]

0.8

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

25

50

75

100

125

150

0.4

175

0

25

50

75

t [hr]

100

125

150

175

125

150

175

t [hr]

(c). bottom point 0.95

0.95

0.8

0.8

0.75

0.75

Sw [-]

0.9 0.85

Sw [-]

0.9 0.85

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

25

50

75

100

t [hr]

125

150

175

0.4

0

25

50

75

100

t [hr]

Figure 4.2: Influence of irrigation ratio on water saturation

4.1 Influence of the Irrigation Ratio

40

Figure 4.3 shows the curves of copper mole fraction in the water phase versus time under different irrigation ratios at the three points. We can see that the copper mole fraction in the water phase decreases more rapidly under higher irrigation ratios than under lower ones. Moreover, just like the change of saturation, the change of copper mole fraction is more sensitive to the decrease of the irrigation ratio. This can be explained by the dependence of concentration change on saturation change. The wash out and desorption processes of copper are influenced by the wetting process of the heap. The faster the heap is wetted or the fresh water is renewed, the faster the copper in the water phase is washed out and the sorption equilibrium is broken and rebuilt. Moreover, from the transport equation, it can be derived that a larger irrigation ratio enhances the advective transport and at the same time enlarges time derivative of mole fraction or concentration. Both effects contribute to a faster decrease of concentration.

4.2 Influence of the Intrinsic Permeability

4.2

41

Influence of the Intrinsic Permeability

Figure 4.4 and Figure 4.5 show the influences of the intrinsic permeability on water saturation and copper concentration in the water phase. The intrinsic permeability is 1.78 · 10−12 m2 for the solid lines, and is raised three times (5.34 · 10−12 m2 ) for the dashed lines, and reduced to one third (5.933 · 10−13 m2 ) for the dash-dot lines. For the long dashed lines, the value is ten times higher (1.78 · 10−11 m2 ) than the solid ones. But it cannot be reduced to less than one third of the solid lines, because otherwise the same problem will occur as when the irrigation ratio is increased too much. That is to say we will also get water saturation larger than one at the top of the heap.

2E-05

1.5E-05

Xcw [mol/mol]

r=10.68 l/(hr. m2) r=5.34 l/(hr. m2) r=21.36 l/(hr. m2)

1E-05

5E-06

0

0

100

200

300

400

500

600

t [hr]

2E-05

2E-05

1.5E-05

1.5E-05

Xcw [mol/mol]

Xcw [mol/mol]

(a) top point

1E-05

1E-05

5E-06

0

5E-06

0

100

200

300

400

t [hr]

(b) middle point

500

600

0

0

100

200

300

400

500

600

t [hr]

(c) bottom point

Figure 4.3: Influence of irrigation ratio on the mole fraction of copper in the water phase

4.2 Influence of the Intrinsic Permeability

42

Contrast to the influence of the irrigation ratio, the models are more sensitive to the increase of the intrinsic permeability. In addtion, a larger intrinsic permeability leads to a smaller maximum saturation and slope. But just like increase of irrigation ratio, increase of the intrinsic permeability also leads to a shorter lower plateau, which means a faster propagation of the front. This is understandable, because if the instinsic permeabiliy remains constant and the irrigation ratio is decreased, not enough pressure is exerted to push the water to fill up the heap to the original maximum saturation, on the other hand, if the irrigation ratio is unchanged, but the intrinsic permeability is increased, more water will permeate through before the heap can be saturated to the original extent. In addition, because the flow velocity is directly proportional to the permeability, it is clear that the front propagates more rapidly under a larger intrinsic permeability. But on the other hand, the storage is reduced when the permeability is increased, because the drainage process is enhanced, while the inflow is not increased by the same factor, which means the gradient of the flux (see Figure 4.1) is reduced. Therefore a flatter slope can be observed. One can also think about this physically. A higher irrigation ratio just exerts a higher pressure which pushes the front stronger into the heap. When the permeability of the porous medium is relatively small compared with the driving force exerted by irrigation, the front behaves more like a shock wave. But when the permeability is relatively high, water can permeate through more easily, therefore on a certain level, it needs more time for the water to accumulate up to the maximum saturation which it can reach and the front behaves more like an expansion wave. The results from the 2p model and the 2p3c model are similar as shown in Figure 4.4.

4.2 Influence of the Intrinsic Permeability

43

(a). top point 0.9

0.85

0.85

0.8

0.8

0.75

0.75

Sw [-]

0.95

0.9

Sw [-]

0.95

0.7

0.7

0.65

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

0.4

80

K=1.78 .10 -12 m2 K=5.933 .10 -13m2 K=1.78 .10 -11m2 K=5.34 .10-12m2

0

20

40

t [hr]

60

.

80

t [hr]

(b). middle point 0.95

0.95

0.9

0.9

0.85

0.85

0.75

Sw [-]

0.8

0.75

Sw [-]

0.8

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

0.4

80

0

20

40

t [hr]

60

80

60

80

t [hr]

(c). bottom point 0.95

0.95

0.8

0.8

0.75

0.75

Sw [-]

0.9 0.85

Sw [-]

0.9 0.85

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

t [hr]

80

0.4

0

20

40

t [hr]

Figure 4.4: Influence of absolute permeability on water saturation

4.2 Influence of the Intrinsic Permeability

44

From Figure 4.5, we can see that, similar to the influence of irrigation ratio, the copper concentration decreases more rapidly under a higher intrinsic permeability. But it is evident that the influence of permeability on the concentration change is not as large as that of the irrigation ratio. It is reasonable, because a larger irrigation ratio not only enhances the advective transport but also the time derivative of concentration, whereas a larger intrinsic permeability only enhances the former one but reduces the latter one. Therefore the total effect of the intrinsic permeability on concentration is weakened. 2E-05

Xcw [mol/mol]

1.5E-05

1E-05

K=1.78 .10 -12 m2 K=5.933 .10 -13m2 K=1.78 .10-11m2 K=5.34 .10-12m2

5E-06

0

0

50

100

150

200

250

.

300

t [hr]

2E-05

2E-05

1.5E-05

1.5E-05

Xcw [mol/mol]

Xcw [mol/mol]

(a) top point

1E-05

1E-05

5E-06

0

5E-06

0

50

100

150

200

t [hr]

(b) middle point

250

300

0

0

50

100

150

200

250

300

t [hr]

(c) bottom point

Figure 4.5: Influence of absolute permeability on the mole fraction of copper in the water phase

4.3 Influence of the Porosity

4.3

45

Influence of the Porosity

In Figure 4.6 and Figure 4.7, the porosity is set to be 0.33 for the solid lines, and is increased and reduced by 36% for the dashed lines and the dash-dot lines respectively. It can be inferred from Figure 4.6 that porosity has no influence on the maximum saturation. It only has an influence on the speed of wetting process. It is becuase porosity does not play a role in the flux term which determines the maximum saturation. However, it is in the storage term, a larger porosity leads to a smaller time derivative of saturation or a flatter slope. Moveover, a porous medium with a higher porosity has larger storage capacity for the fluid, therefore the front propagates more slowly. Again, the results of the 2p and 2p3c models are similar.

4.3 Influence of the Porosity

46

(a). top point 0.95

0.95

0.8

0.8

0.75

0.75

Sw [-]

0.9 0.85

Sw [-]

0.9 0.85

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

porosity=0.33 porosity=0.45 porosity=0.21

0

20

40

60

80

100

0.4

120

0

20

40

t [hr]

60

80

100

120

80

100

120

80

100

120

t [hr]

(b). middle point 0.95

0.95

0.9

0.9

0.85

0.85

Sw [-]

0.8 0.75

Sw [-]

0.8 0.75 0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

80

100

0.4

120

0

20

40

t [hr]

60

t [hr]

(c). bottom point 0.95

0.95

0.9

0.9

0.85

0.85

Sw [-]

0.8 0.75

Sw [-]

0.8 0.75 0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

t [hr]

80

100

120

0.4

0

20

40

60

t [hr]

Figure 4.6: Influence of porosity on water saturation

4.3 Influence of the Porosity

47

From Figure 4.6, we can see that the smaller the porosity, the faster the decrease of copper concentration. Because the porosity is involved in the retardation factor (see formulation (2.31)), and a smaller porosity leads to a smaller retardation faster, thus a larger time derivative of concentration. Moreover, since the front propagation is accelerated by a smaller porosity, the advective transport of copper is enhanced as well. Since the value of porosity is much larger than that of the concentration of copper in the water phase, any small change of the porosity can lead to a considerable change on the time derivative of the concentration, and at the same time the value of porosity also influence the advective transport of copper. Those two factors influence the concentration change in the same direction. Therefore the change of copper concentration is more sensitive to the porosity than to other parameters. 2E-05

1.5E-05

Xcw [mol/mol]

porosity=0.33 porosity=0.45 porosity=0.21

1E-05

5E-06

0

0

50

100

150

200

250

300

350

400

t [hr]

2E-05

2E-05

1.5E-05

1.5E-05

Xcw [mol/mol]

Xcw [mol/mol]

(a) top point

1E-05

1E-05

5E-06

0

5E-06

0

50

100

150

200

250

t [hr]

(b) middle point

300

350

400

0

0

50

100

150

200

250

300

350

400

t [hr]

(c) bottom point

Figure 4.7: Influence of porosity on the mole fraction of copper in the water phase

4.4 Influence of the Pc − Sw Curve

4.4

48

Influence of the Pc − Sw Curve

There are three parameters which determine the capillary pressure-saturation curve. They are residual water saturation, Van Genuchten parameter α [-] and Van Genuchten parameter n [-]. For the sensitivity test, these three parameters are changed one by one and the corresponding results are pesented separately and an overall interpretation is given at the end of this section.

4.4.1

Influence of the residual water saturation

Since the original value of residual water saturation Swr is zero, we can only increase it. In Figure 4.8 and 4.9, the value is 0.0 for the solid lines, 0.1 for the dashed lines and 0.2 for the dash-dot lines.

4.4 Influence of the Pc − Sw Curve

49

(a). top point 0.95

0.95

0.8

0.8

0.75

0.75

Sw [-]

0.9 0.85

Sw [-]

0.9 0.85

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

Swr=0.0 Swr=0.1 Swr=0.2

0

20

40

60

0.4

80

0

20

40

t [hr]

60

80

60

80

60

80

t [hr]

(b). middle point 0.9

0.85

0.85

0.8

0.8

0.75

0.75

Sw [-]

0.95

0.9

Sw [-]

0.95

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

0.4

80

0

20

40

t [hr]

t [hr]

(c). bottom point 0.9

0.85

0.85

0.8

0.8

0.75

0.75

Sw [-]

0.95

0.9

Sw [-]

0.95

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

t [hr]

80

0.4

0

20

40

t [hr]

Figure 4.8: Influence of Swr on water saturation

4.4 Influence of the Pc − Sw Curve

50

2E-05

Swr=0.0 Swr=0.1 Swr=0.2

Xcw [mol/mol]

1.5E-05

1E-05

5E-06

0

0

50

100

150

200

250

300

t [hr]

2E-05

2E-05

1.5E-05

1.5E-05

Xcw [mol/mol]

Xcw [mol/mol]

(a) top point

1E-05

1E-05

5E-06

0

5E-06

0

50

100

150

200

t [hr]

(b) middle point

250

300

0

0

50

100

150

200

250

300

t [hr]

(c) bottom point

Figure 4.9: Influence of Swr on the mole fraction of copper in the water phase

4.4 Influence of the Pc − Sw Curve

4.4.2

51

Influence of the Van Genuchten α

In Figure 4.10 and Figure 4.11, the value of α for the solid lines are 1.35 · 10−4 [1/P a], the values for the dashed lines and dash-dot lines are increased and decreased by 22%. (a). top point 0.95

0.95

0.8

0.8

0.75

0.75

Sw [-]

0.9 0.85

Sw [-]

0.9 0.85

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

VG_alpha=1.35 .10-4 Pa-1 VG_alpha=1.05 .10 -4 Pa-1 VG_alpha=1.65 .10 -4 Pa-1

0

20

40

60

0.4

80

0

20

40

t [hr]

60

80

60

80

60

80

t [hr]

(b). middle point 0.95

0.95

0.85

0.8

0.8

0.75

0.75

Sw [-]

0.9

0.85

Sw [-]

0.9

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

0.4

80

0

20

40

t [hr]

t [hr]

(c). bottom point 0.95

0.95

0.9

0.9

0.85

0.85

Sw [-]

0.8 0.75

Sw [-]

0.8 0.75 0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

t [hr]

80

0.4

0

20

40

t [hr]

Figure 4.10: Influence of Van Genuchten α on water saturation

4.4 Influence of the Pc − Sw Curve

52

2E-05

1.5E-05

Xcw [mol/mol]

VG_alpha=1.35 .10 -4 Pa-1 VG_alpha=1.05 .10 -4 Pa-1 VG_alpha=1.65 .10 -4 Pa-1

1E-05

5E-06

0

0

50

100

150

200

250

300

t [hr]

2E-05

2E-05

1.5E-05

1.5E-05

Xcw [mol/mol]

Xcw [mol/mol]

(a) top point

1E-05

1E-05

5E-06

0

5E-06

0

50

100

150

200

t [hr]

(b) middle point

250

300

0

0

50

100

150

200

250

300

t [hr]

(c) bottom point

Figure 4.11: Influence of Van Genuchten α on copper mole fraction in the water phase

4.4 Influence of the Pc − Sw Curve

4.4.3

53

Influence of the Van Genuchten n

In Figure 4.12 and 4.13, n is set to 1.411 [−] for the solid lines, and increased and decreased by 7% for the dashed lines and dash-dot lines. (a). top point 0.95

0.95

0.8

0.8

0.75

0.75

Sw [-]

0.9 0.85

Sw [-]

0.9 0.85

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

VG_n=1.411 VG_n=1.311 VG_n=1.511

0

20

40

60

80

100

0.4

120

0

20

40

t [hr]

60

80

100

120

80

100

120

80

100

120

t [hr]

(b). middle point 0.95

0.95

0.85

0.8

0.8

0.75

0.75

Sw [-]

0.9

0.85

Sw [-]

0.9

0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

80

100

0.4

120

0

20

40

t [hr]

60

t [hr]

(c). bottom point 0.95

0.95

0.9

0.9

0.85

0.85

Sw [-]

0.8 0.75

Sw [-]

0.8 0.75 0.7

0.65

0.7

0.65

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

20

40

60

t [hr]

80

100

120

0.4

0

20

40

60

t [hr]

Figure 4.12: Influence of Van Genuchten n on water saturation

4.4 Influence of the Pc − Sw Curve

54

2E-05

1.5E-05

Xcw [mol/mol]

VG_n=1.411 VG_n=1.311 VG_n=1.511

1E-05

5E-06

0

0

50

100

150

200

250

300

t [hr]

2E-05

2E-05

1.5E-05

1.5E-05

Xcw [mol/mol]

Xcw [mol/mol]

(a) top point

1E-05

1E-05

5E-06

0

5E-06

0

50

100

150

200

250

300

0

0

t [hr]

(b) middle point

50

100

150

200

250

300

t [hr]

(c) bottom point

Figure 4.13: Influence of Van Genuchten n on copper mole fraction in the water phase

4.4.4

Analysis and Interpretation

From the above figures, we can conlude that the Van Genuchten n has the biggest influence on the evolvement of water saturation, while the residual water saturation has the smallest influence. In the balance equation, capillary pressure-saturation relation also influences the flux term via relative permeability and pressure gradient, therefore it has impact on both the maximum saturation and flow velocity. It can also be inferred from the transport equation that capillary pressure-saturation relation influences the advective transport and the time derivative of concentration as well. The lower the capillary pressure, the faster the flow velocity and the quicker the concentration decrease. The deviations of the Sw − t curves and the Xwc − t curves under different

4.5 Influence of the Adsorption Coefficient

55

pc − Sw relations are not distinct, because the change of the pc − Sw parameters are very small. We do this because with the original parameter set, the capillary pressure almost goes into infinity when water saturation is about 0.38. Therefore we can not change the parameters too much to avoid a too high capillary pressure at a relatively large saturation, which is not realistic. And particularly, one must notice that the time scale for the Xwc −t curves is quite large. That is why almost no differences between the Xwc − t curves can be detected. Figure 4.14 shows the pc − S − w curves with different values which are used for the pc − Sw parameters. with original values: Swr=0.0 [-] alpha=1.35E-4 [Pa-1]; n=1.411 [-] Swr changed to 0.2 Swr changed to 0.1 alpha changed to 1.05E-4 alpha changed to 1.65E-4

50000

40000

pc [Pa]

n changed to 1.311 n changed to 1.511

30000

20000

10000

0

0

0.2

0.4

Sw [-]

0.6

0.8

Figure 4.14: pc − Sw curves for sensitivity test.

4.5

Influence of the Adsorption Coefficient

In Figure 4.15, The left column of three pictures shows the results of water saturation evolvement with different adsorption coefficients, and the right one shows those of the copper mole fraction in the water phase. The value of the adsorption coefficient kd is set to 3.5 · 10−5 m3 /kg for the solid lines, 8.67 · 10−5 m3 /kg for the dashed lines and zero for the dash-dot lines. For the non-adsorption case (represented by dash-dot lines), we assumed the same initial copper mole fraction in the water phase as the original value (the value for the solid lines). Then the comparison between the solid lines and the dash-dot lines gives us an insight about the retardation effect of the adsorption equilibrium as discussed in section 2.2.4. Whereas we assume that the solid lines and the dashed lines have the same initial copper content in the solid matrix. Then according to adsorption equilibrium and linear sorption isotherm (see equation (2.20)), the initial copper mole fraction in the water phase for the dashed lines can be

4.5 Influence of the Adsorption Coefficient

56

recalculated. (a). top point 2E-05

0.95 0.9 0.85

1.5E-05 0.8

-5

kd=3.5 .10 m3/kg kd=0.0 m3/kg kd=8.67 .10 -5m3/kg

Xcw [mol/mol]

Sw [-]

0.75 0.7

1E-05

0.65 0.6 0.55

5E-06

0.5 0.45 0.4

0

20

40

60

80

0

100

0

100

200

t [hr]

300

400

300

400

t [hr]

(b). middle point 2E-05

0.95 0.9 0.85

1.5E-05 0.8

Xcw [mol/mol]

Sw [-]

0.75 0.7

1E-05

0.65 0.6 0.55

5E-06

0.5 0.45 0.4

0

20

40

60

80

0

100

0

100

200

t [hr]

t [hr]

(c). bottom point 2E-05

0.95 0.9 0.85

1.5E-05 0.8

Xcw [mol/mol]

Sw [-]

0.75 0.7

1E-05

0.65 0.6 0.55

5E-06

0.5 0.45 0.4

0

20

40

t [hr]

60

80

100

0

0

50

100

150

200

250

300

350

400

t [hr]

Figure 4.15: Influence of adsorption coefficient on water saturation and copper mole fraction in the water phase

4.6 Conclusions

57

It can be inferred from Figure 4.15 that the change of the adsorption coefficient has almost no influence on the behavior of fluid flows, but a considerable influence on the change of mole fraction. It is because that although the adsorption coefficient is involved in the storage term (or retardation factor, see equation (2.30)), its value is several magnitude smaller than that of the saturation but in the same range as the value of mole fraction. In addition, the retardation factor is directly proportional to the adsorption coefficient (see formulation (2.31)). Therefore the higher the adsorption coefficient the slower the decrease of copper concentration in the water phase. This is reflected by the slope of the Xwc − t curves.

4.6

Conclusions

From the above results and discussions of the sensitivity test, we can draw the following conclusions: • Water saturation is most sensitive to the irrigation ratio. • Copper concentration is most sensitive to the porosity. • There are some correlations between the effects of the parameters on water saturation and copper concentration. • The value of irrigation ratio, permeability and capillary pressure have influences on both the maximum saturation and flow velocity, while porosity can only influence the flow velocity. • The value of adsorption coefficient has almost no influence on the water saturation but a considerable influence on the copper concentration. • The effects of each parameter on the 2p model are similar as on the 2p3c model.

Chapter 5 Simulations under Heterogeneous Conditions Because sometimes a layered configuration of heap is used in some mining companies to enhance the efficiency of leaching process, our model is further extended to simulate the process with a heterogeneous domain which is divided into three horizontal layers with different material properties. Since a complete set of measured parameters for a layered heap are not available, we take the reference of parameters from literatures. The values of the parameters which are different for the three layers are listed in Table 5.1. Other parameters, such as fluid density, viscosity and so on, take the values from Table 3.1 in Chapter 3. The set-up of the domain is shown in Figure 5.1. parameter φ [-] K [m2 ] kd [m3 /kg] Swr [-] n [−] α [1/P a]

zone-1 0.13 1.78 · 10−12 8.2 · 10−5 0.01 1.85 8.0 · 10−4

zone-2 0.14 2.15 · 10−12 5.8 · 10−5 0.01 1.90 8.58 · 10−4

zone-3 0.15 2.50 · 10−12 3.5 · 10−5 0.01 1.95 9.02 · 10−4

Table 5.1: Values of parameters describing the properties of the heterogeneous domain. [4] The same boundary and initial conditions (except that the initial Xxc has to be recalculated for different layers according to the sorption equilibrium) as defined in Chapter 3 are used in the following simulations. Figure 5.2 shows the color contour of evolvement of water saturation. The left column of six pictures shows the results from the 2p model, whereas the right one shows that of the 2p3c model.

59

1.5 m

zone−1 zone−2

2.5 m

zone−3

1m 5m 45

o

45

o

25 m

Figure 5.1: Set-up of the heterogeneous domain.

t = 5.22 hr

t = 4.29 hr

t = 9.69 hr

t = 8.76 hr

t = 15.66 hr

t = 14.73 hr

t = 24.60 hr

t = 23.67 hr

t = 41.01 hr

t = 40.08 hr

Sw: 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85

Figure 5.2: Evolvement of water saturation.

60

0.85

0.85

0.8

0.8

0.75

0.75

0.7

0.7

0.65

0.65

Sw [-]

Sw [-]

In Figure 5.3, the left picture shows the change of the water saturation over time at the three sample points (the coordinates of which are the same as in Chapter 3) of the 2p model, and the right one shows that of the 2p3c model.

0.6

0.6

0.55

0.55

0.5

0.5

0.45

0.45

0.4

0

10

20

t [hr]

30

0.4

0

10

20

30

t [hr]

Figure 5.3: Change of saturation over time at the three sample points. From the color contour of water saturation evolvement, we can distinguish the different layers by the shape of the fronts propagating downwards. When the fluid flows into the lower layers with higher permeabilities and lower capillary pressure (see Table 5.1), the width of the front is getting narrower. It is because with a higher permeability, the vertical velocity becomes larger, thus the fluid can flow downwards more quickly and at the same time the horizontal flow induced by capillary suction becomes less significant. Moreover, as mentioned in the sensitivity test, a higher permeability leads to a smaller extent of saturation. That is why from the top to bottom, the saturation of the layers becomes smaller and smaller. This can also be proved by Figure 5.3. Both Figure 5.2 and Figure 5.3 show that the results from the 2p and 2p3c models are very consistent. Figure 5.4 shows the color contour of the evolvement of copper mole fraction in the water phase of the 2p3c model. Figure 5.5 shows the change of copper mole fraction in the water phase at the three sample points. One cannot find a fixed rule to tell the differences of the shapes between those curves by the influences of the parameters. Because more than one parameter takes different values in the three layers, and some of them influence the curves in opposite directions. For example, a higher adsorption coefficient leads to a slower decrease of copper concentrations, while smaller porosity leads to a faster decrease.

61

t = 0.00 hr

t = 77.73 hr

t = 8.76 hr

t = 109.05 hr

t = 47.54 hr

t = 225.38 hr

Sw: 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85

Figure 5.4: Evolvement of copper mole fraction in the water phase.

2E-05

Xcw [mol/mol]

1.5E-05

1E-05

5E-06

0

0

50

100

150

t [hr]

Figure 5.5: Change of copper mole fraction in the water phase at the three sample points.

Chapter 6 Summary and Outlook Copper is a very important element to human life. Therefore, the extraction efficiency of copper receives much attention and improvement of the process is expected. There are different alternatives to acchieve this target and modelling is a relatively cheap and efficient way to evaluate the recovery rate of copper and help the optimization of operation conditions. For this purpose, a two-phase model and a two-phase three-component model are set up to simulate the flow behavior and transport of copper with adsorption in the heap leaching process which is considered as a superior way to other copper processing methods. However, the existing chemical reaction between the sulfuric acid and copper ore during the leaching process is not considered in our model, only adsorption is modelled by the linear sorption law. Thus the acid solution is simpified by the water phase and the component acid is simplified by the component water. Chapter 2 gives the basics of the model concepts and mathematical formulations as well as selected numerical schemes for both the two-phase model and the two-phase three-component model. In Chapter 3, proper initial and boundary conditions are assigned to the two dimenstional and homogeneous domain. Especially, at the bottom of the domain, since in the reality the fluids are in contact with the atmosphere, a free flow boundary condition is assigned and a quasi-simulation of the ambient environment is made by extention of the domain with virtual properties in the extended part. The simulation results of water saturation evolvement from both the two-phase model and two-phase three-component model are very consistent with each other. The little difference can be explained by different numerical accuracies due to different number of primary variables taken in the two models. In addition, the two-phase three-component model shows reasonable result about the evolvement of copper concentration in the water phase. Initially, the results of our models were supposed to be compared with the results of the model developed by the Chilean reseacher, Mr. E. Cariaga. But

63

unfortunately, his model does not work due to convergency problem. Chapter 4 investigates the effects of some parameters on the fluid flow behavior and the concentration change of copper by a sensitivity study. This step can also be considered as the basis for the improvement of operation conditions of heap leaching process. The results show that water saturation is most sensitive to the irrigation ratio. If a mining company wants to enhance the wetting speed and saturation extent of the heap which in turn accelerate the concentration change of copper in the water phase, the most efficient way is to enlarge the irrigation ratio. And copper concentration is most sensitive to the porosity. However, we can not yet give any advice about the most efficient way to enhance the recovery rate of copper at this stage, since the chemical reaction plays an important role in the real extraction process, which is not included in our model. The effects of the paramters under investigation on the change of water saturation showed by the two-phase model and the two-phase three-component model are very similar. Finally, in order to test the simulation ability of our models under a layered configuration which is sometimes used in the mining industry, our model is further extended to model the leaching process with a heterogeneous domain which is divided into three horizontal layers. Since a complete measured parameter set is not available, we referred to the literatures to set up an arbitrary layered domain. Again both the two-phase and two-phase three-component models give consistent and reasonable results. All of the above simulations indicate that the flow behavior and transport of copper with adsorption in the heap leaching process can be well simulated by our models under different conditions. But this is only a preliminary trail of the industrial application of the model, still a lot of improvement can be done in the future development: - Our model cannot yet produce the real copper concentration change over time, because of neglection of the reaction process. So the inclusion of the reaction between sulfuric acid and copper ore is of the most significance to evaluate the extraction of copper from solid matrix to the liquid phase. - More realistic parameter values and boundary conditions should be investigated. - Inclusion of heterogeneities e.g. fractures and shear zones is necessary to investigate the preferential flow or funnel flow phenomena in the real heap. - The influence of different dimensions of the heap on the flow pattern and extraction efficiency can be studied to optimize the shape of the heap which can be of industrial interest.

64

- Aside from the chemical reaction process, other processes like bio-reaction which has recently received much attention for heap leaching process of mines and can be included. - Finally, We need the validation of the model to test the performance of all the above mentioned functions. The validation should be based on experiments.

Bibliography [1] E. Anozie. Thermal Effects of Carbon Dioxide Sequestration in the Subsurface. Diplomarbeit, 2005. [2] Peter Bastian, K. Birken, K. Johannsen, S. Lang, K. Eckstein, N. Neuss, H. RentzReichert, and C. Wieners. UG - A Flexible Software Toolbox for Solving Partial Differential Equations. Computing and Visualization in Science, 1(1), S. 27-40, 1997. [3] R.H. Brooks and A.T. Corey. Hydraulic Properties of Porous Media. In Fort Collins, editor, Hydrol. Pap., volume 3, Colorado State University, 1964. [4] Olaf A. Cirpka. Environmental Fluid Mechanics I — Flow in Natural Hydrosystems. Universit¨at Stuttgart · Institut f¨ ur Wasserbau, 2003. [5] Olaf A. Cirpka. Environmental Fluid Mechanics II — Solute and Heat Transport in Natural Hydrosystems. Universit¨at Stuttgart · Institut f¨ ur Wasserbau, 2004. [6] H. Class. Lecture Notes of ”Application of Single-phase and Multiphase Models to Environmental and Technical Problems”. 2005. [7] J.E.Kiel D J.A.van Zyl, I.P.G.Hutchison. Introduction to Evaluation, Design and Operation of Precious Metal Heap Leaching Protects. Society of Mining Engineering, Inc, Colorado, 1988. [8] H¨afner et al. Geohydrodynamische Erkundung von Erd¨ol-, Erdgas- und Grundwasserlagerst¨atten. Technical report 1, WTI des ZGI, Berlin, 1985. [9] M.T. Van Genuchten. A Closed-Form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Sci. Soc. Am. J., 44:892–898, 1980. [10] H. Class R. Helmig and P. Bastian. Numerical Simulation of Non-isothermal Multiphase Multicomponent Processes in Porous Media. Advances in Water Resources, 25:533–550, 2002. [11] R. Helmig. Multiphase Flow and Transport Processes in the Subsurface. SpringerVerlag, 1997. 65

BIBLIOGRAPHY

66

[12] R. Helmig, C. Braun, and M. Emmert. MUFTE - A numerical modelf or simulation of multiphase flow processes in porous and fractured-porous media. Programmdokumentation (HG 208). Technical Report 94/3. Technical report, Institut f¨ ur Wasserbau, Universit¨at Stuttgart, 1994. [13] R. Helmig and H.Class. lecture notes of Hydromechanics. 2004. [14] http://www.codelco.cl. [15] R.J. Millington and J.P. Quirk. Permeability of Porous Solids. Department of Agronomy and Agricultural Chemisty, Trans. Faraday, Soc., 57:1200–1207, 1961. [16] S. Orr. Enhanced heap leaching–Part 1: Insights. Mining Engineering, 55(9):49– 55, 2002. [17] J.M. Prausnitz R.C. Reid and B.E. Poling. Properties of Gases and Liquids. McGraw-Hill, 1987. [18] A.E. Scheidegger. The Physics of Flow Through Porous Media. University of Toronto Press, 1974. [19] R.A. Feddes S.P. Neumann and E. Bresler. An Approach for Analyzing Transport in Strongly Saturated-Unsaturated Soils Considering Water Uptake by Plants. Developments of methods, tools and solutions for unsaturated flow, third annumal report (part 1), 1974.