Batch Solvent Extraction of Bitumen from Oil Sand. Part 1: Theoretical

Apr 3, 2017 - In the following results, effects of the solvent fraction on viscosity and density of diluted bitumen were taken into account. These pro...
1 downloads 6 Views 3MB Size
Download PDF
Article pubs.acs.org/EF

Batch Solvent Extraction of Bitumen from Oil Sand. Part 1: Theoretical Basis Merouane Khammar* and Yuming Xu CanmetENERGY, Natural Resources Canada, 1 Oil Patch Drive, Devon, Alberta T9G 1A8, Canada ABSTRACT: A model for flow through a porous medium with variable saturation is presented in this work to predict the volume and flow rate of extracted solvent-diluted bitumen from mineral solids during centrifugal filtration. The numerical method used for solving the flow equation is based on the mass-conservative-modified Picard iteration scheme and the Anderson acceleration algorithm. It is found that bitumen extraction involves filtration and desaturation phases, and the transition between them can be identified from the evolution with time of extracted liquid volume and flow rate. An analytical description of the filtration phase is developed to predict the evolution of the liquid layer thickness above the cake during centrifugation. Numerical and analytical models presented in this work can be used in a fitting procedure to determine hydraulic properties of oil sand cake. Effects of centrifugal force and solvent mass fraction on the predicted fraction of extracted diluted bitumen are investigated.

1. INTRODUCTION Bitumen extraction from mined oil sands is currently based on the Clark hot-water extraction process. This technology can have significant impacts on the environment as a result of high water consumption, generation of large tailing ponds, and high energy requirement. The solvent extraction process is currently being developed as an alternative non-aqueous technology to reduce the environmental impact of bitumen extraction from oil sands. In this process, solvent is mixed with oil sand to dissolve bitumen. Solvent-diluted bitumen is then separated from sand particles before being disposed to the environment. An important challenge associated with solvent bitumen extraction is disposing of the separated sand to the environment with the lowest possible bitumen and solvent contents.1 Solvent can easily be recovered from diluted bitumen. However, the extraction of diluted bitumen from the wet sand was found to be challenging because capillary and adsorptive forces bind liquid to the porous structure formed by the mineral particles. After complete bitumen dissolution in solvent, the bitumenrecovered fraction corresponds to the solvent-diluted bitumenextracted fraction. Thus, the problem of maximizing bitumen extraction in a solvent extraction process is analogous to the development of an efficient continuous solid−liquid separation process of solvent-diluted bitumen from oil sand. Meadus et al.2 highlighted the importance of solid−liquid separation in the solvent extraction process. It was demonstrated that the efficiency of solid−liquid separation governs the overall recovered fraction of bitumen and solvent from oil sand. Two separation processes were compared: screening operation and spherical agglomeration process. Screening operation is a solid−liquid separation process based on gravity drainage, where a mixture of diluted bitumen + sand is drained on a stationary screen mesh under the effect of gravitational force. Spherical agglomeration is based on adding water to a mixture of solvent and diluted bitumen in a rotating mill. Mineral particles are bound by water to form granules, and diluted bitumen in the pores of the granules is expelled both by water displacement and granule compaction. Published XXXX by the American Chemical Society

In the screening operation, liquid is retained in the pores by capillary and adsorptive forces in the porous structure and gravity alone is unable to completely desaturate liquid retained in the pores. Centrifugation enhances the magnitude of gravitational force by orders of magnitude. Wu and Dabros3 investigated batch centrifugal extraction of solvent-diluted bitumen from mineral solids. A vessel was designed and built; it was divided into two chambers by a filter paper separating them. In the upper chamber, an oil sand sample was placed and solvent was added. Oil sand was mixed with solvent for 1 h in a wrist action shaker before being placed in a centrifuge and spun for 10 min. High bitumen recovery was achieved after 10 min of batch centrifugation. Effects of the solvent/bitumen ratio, centrifugal force, and number of extraction stages on bitumen extraction were investigated, and it was found that extraction efficiency increased with centrifugal force and solvent/bitumen ratio. Faradonbeh et al.1 addressed the problem of diluted bitumen recovery trapped in mineral particles. They investigated experimentally and theoretically liquid extraction from a packed bed of sand. A mathematical model was developed to predict transient and equilibrium liquid saturation profiles along the sand pack for both gravity drainage and pressure desaturation. Results for pressure desaturation experiments concluded that the use of a semi-impermeable membrane at the bottom of the sand bed decreased significantly average liquid saturation and improved liquid recovery. A continuous process of solvent-diluted bitumen extraction from oil sands based on centrifugal filtration should maximize diluted bitumen recovery before directing the wet mineral particles to a solvent recovery process and disposing of the dry particles to the environment. The design of a continuous process of solvent-diluted bitumen extraction from oil sands requires the development of the theoretical basis governing the Received: December 7, 2016 Revised: March 6, 2017 Published: April 3, 2017 A

DOI: 10.1021/acs.energyfuels.6b03252 Energy Fuels XXXX, XXX, XXX−XXX

Energy & Fuels

Article

Figure 1. Schematic representation of the two stages of diluted bitumen extraction: (a) filtration and (b) desaturation.

solved in a gravitational field, the solution produced significant mass balance errors. They proposed a modified Picard’s massconservative discretization method to solve the mixed form of the Richards equation expressed in terms of both the liquid content and pressure head. Simunek and Nimmo8 used a numerical discretization of the mixed form of the Richards equation in a centrifugal field based on Celia et al. massconservative method.11 The numerical solver simulated multirotational transient flow to estimate the hydraulic parameters of a packed sand sample. Walker et al.12,13 and Lott et al.14 showed that the combination of the Anderson acceleration algorithm and the modified Picard iteration significantly improves the convergence speed and robustness of the iterative scheme when applied to solving the Richards equation in a gravitational field. This algorithm was, however, not applied in a centrifugal field. The determination of hydraulic properties using batch centrifugal filtration requires both modeling predictions and online measurements. In this work, a model for solvent-diluted bitumen extraction in the centrifugal field is developed. The sample modeled corresponds to introducing an oil sand cake sample into a filtration cell and adding a layer of diluted bitumen above it. The pores of the oil sand cake are initially saturated with solvent-diluted bitumen. The extracted fraction of solvent-diluted bitumen is predicted for different rotation speeds and solvent fractions. A method is developed to identify the transition time between filtration and desaturation and obtain hydraulic properties of the diluted bitumen−oil sand cake system.

dependence of the extraction efficiency upon cake properties, centrifugal force, and solvent/bitumen ratio. The numerical solution used to solve the conservation equation should be mass conservative to accurately predict the fraction of extracted liquid and pore saturation distribution. Model results of batch extraction combined with experimental data can be used to extract flow properties of solvent-diluted bitumen through oil sand. These properties are required for the prediction of bitumen extraction efficiency in a continuous extraction process. Properties of liquid flow through porous media with variable saturation are characterized by two nonlinear relationships: the liquid retention curve and hydraulic conductivity function. The liquid retention curve describes the relationship between the volumetric liquid content, θ, and the matric potential, ψ.4 The volumetric liquid content, θ, characterizes pore saturation, (S = θ/θs), whereas the matric potential, ψ, is the potential energy per unit volume that results from capillary and adsorptive forces between liquid and the solid matrix.4 These forces bind liquid to the porous matrix and lower its potential energy below that of bulk liquid.4 The hydraulic conductivity function K defines the dependence of the hydraulic conductivity upon matric potential, ψ, or liquid content, θ. The hydraulic conductivity function decreases with the liquid content and is valid only when the liquid phase is continuous in the porous layer. Batch centrifugal filtration has been used as an accelerated experimental technique to determine the hydraulic properties of liquid flow through porous media.5−8 In a centrifugal field, pressure and pore saturation can vary significantly between the top and bottom of a porous sample. Several researchers have worked on solving the equations of liquid flow through porous media in a centrifugal field. Bear et al.9 solved the equations of fluid flow through a deformable porous medium in the presence of a centrifugal field using the finite difference method. They calculated liquid content and pressure head profiles for a variety of simulations. Porosity changes during centrifugation were found to be negligible in most applications. Nimmo10 numerically solved the Richards equation for flow through incompressible porous medium with variable saturation. Good agreement was reported between experimental water content data and results of the numerical solution of the flow equation during centrifugal filtration. Celia et al.11 showed that when the Richards equation, expressed in terms of pressure head, was

2. THEORY The model described in this work is developed and solved for the prediction of the flow rate and volume of solvent-diluted bitumen extracted by batch centrifugal filtration. The following assumptions are used in the model: (1) Before centrifugation, the oil sand−solvent mixture is assumed to consist of a compact incompressible sand cake saturated and immersed in diluted bitumen. (2) Diluted bitumen is assumed to be an incompressible fluid. (3) Two stages of bitumen extraction are considered (Figure 1): (a) Filtration: a layer of diluted bitumen exists above the saturated sand cake. The thickness of this layer decreases with time during centrifugation. (b) Desaturation: there is no layer of diluted bitumen above the sand cake. The cake of sand desaturates during centrifugation until an equilibrium distribution of pore saturation is obtained in the cake. (4) Unsaturated flow B

DOI: 10.1021/acs.energyfuels.6b03252 Energy Fuels XXXX, XXX, XXX−XXX

Energy & Fuels

Article

properties of diluted bitumen through the sand cake are assumed to be described by the van Genuchten equations. (5) van Genuchten parameters used in this work are independent of the solvent mass fraction in bitumen. However, the model can be applied to cases where the parameters are known for different solvent mass fractions in bitumen. For one-dimensional flow in a centrifugal field, the flow rate q through a porous cake sample is given by the following Darcy− Buckingham equation:8

q=−

⎞ K (θ ) ⎛ dψ ⎜ − ρω2r ⎟ ⎠ ρg ⎝ dr

approximation in time and the modified Picard iteration scheme. Superscripts n and k indicate time and iteration level, respectively. The discretization of the left- and right-hand sides of eq 3 is given by the following equations:

(1)

(5)

q n + 1, k + 1 − qin + 1, k + 1 ∂q = i+1 ∂z Δz

(6)

Subscript “i” indicates spatial discretization in the z direction. The term θni + 1,k + 1 is expanded in a truncated Taylor series with respect to h around the point hni + 1,k as follows:

where ψ is the matric potential given by ψ = P − P0, in which P0 and P are the atmospheric pressure and liquid pressure, respectively, ρ is the density, ω is the rotational speed, K is the hydraulic conductivity, and r is the radial coordinate. The Darcy−Buckingham equation can be expressed in terms of the pressure head in the z coordinate, originating from the filter radial position, as follows (Figure 2):

⎛ dh ω2(r0 − z) ⎞ ⎟ q = − K (θ)⎜ + g ⎠ ⎝ dz

⎛ θ n + 1, k + 1 − θ n ⎞ ∂θ i ⎟⎟ = ⎜⎜ i t ∂t Δ ⎠ ⎝

θin + 1, k + 1 = θin + 1, k +

⎛ dθ ⎞n + 1, k n + 1, k + 1 ⎜ ⎟ (hi − hin + 1, k ) + o(δ 2) ⎝ dh ⎠i

(7) If all nonlinear terms are neglected, the discretization of the time derivative of the liquid content term becomes

(2)

⎛ θ n + 1, k + 1 − θ n ⎞ ⎛ θ n + 1, k − θ n ⎞ ⎛ C n + 1, k ⎞ i i ⎟⎟ + ⎜⎜ i ⎟⎟ = ⎜⎜ i ⎜⎜ i ⎟⎟ Δt Δt ⎠ ⎝ Δt ⎠ ⎠ ⎝ ⎝

where h = ψ/ρg is the pressure head.

(hin + 1, k + 1 − hin + 1, k )

(8)

where

C=

dθ dh

(9)

A staggered grid is used for the flow rates as shown in Figure 3. The flow rates around node i are discretized as follows:

Figure 2. Control volume in a rotating cake with variable saturation.

Figure 3. Discretization of the numerical domain across cake thickness.

The continuity equation of liquid flow through the cake is given by Figure 2

∂q ∂θ =− ∂t ∂z

⎧⎛ n + 1, k + 1 − h n + 1, k + 1 ⎞ i 1, k ⎪ hi + 1 ⎟⎟ ⎨⎜⎜ qin++11, k + 1 = − K in++1/2 ⎪ Δz ⎠ ⎩⎝ ⎫ 2 ⎪ ω ⎬ + (r0 − zi + 1/2)⎪ g ⎭

(3)

where θ is the cake volumetric liquid content. After substitution of eq 2 into eq 3, the following mixed form of the Richards equation is obtained:

⎛ dh ω (r0 − z) ⎞⎞ ∂θ ∂⎛ ⎜⎜K (θ )⎜ ⎟⎟⎟ = + g ∂t ∂z ⎝ ⎝ dz ⎠⎠

(10)

2

⎧⎛ n + 1, k + 1 − h n + 1, k + 1 ⎞ 1, k ⎪ hi i−1 ⎜⎜ ⎟⎟ ⎨ qin + 1, k + 1 = − K in−+1/2 ⎪ Δz ⎠ ⎩⎝

(4)

Equation 4 is the general equation that governs flow through a cake during batch centrifugal filtration and desaturation. It is highly nonlinear as a result of the dependence of the saturation and the hydraulic conductivity upon the pressure head. The discretization of the mixed form of the Richards equation leads to a system of coupled nonlinear equations that must be solved at each time step, with hydraulic conductivity changing by orders of magnitude with time and across cake thickness. Therefore, a robust numerical algorithm should be adopted for the solution of eq 4. 2.1. Numerical Solution of the Flow Equation. The numerical scheme adopted in this work to solve the Richards equation for flow in a centrifugal field follows the modified Picard mass-conservative method proposed by Celia et al.11 This algorithm linearizes the saturation to a first-order Taylor expansion with respect to the pressure head h. The discretization of the mixed form of the Richards equation in a centrifugal field is developed using the backward Euler

+

⎫ ⎪ ω2 (r0 − zi − 1/2)⎬ ⎪ g ⎭

(11)

where ⎧ z + zi , z = i+1 ⎪ ⎪ i + 1/2 2 ⎨ zi + zi − 1 ⎪ , ⎪ zi − 1/2 = ⎩ 2

Ki+1 + Ki 2 K + Ki−1 = i 2

K i + 1/2 = K i − 1/2

(12)

Using eqs 5−11, eq 3 gives the following discretization: C

DOI: 10.1021/acs.energyfuels.6b03252 Energy Fuels XXXX, XXX, XXX−XXX

Energy & Fuels

Article

For the case of cake desaturation (Figure 1b): The flow rate qB at the top of the cake is zero. The boundary condition is given by

1, k 1, k ⎫ n + 1, k ⎛ K n + 1, k ⎞ ⎧ K in−+1/2 + K in++1/2 ⎪C ⎪ ⎜⎜− i − 1/2 ⎟⎟(hin−+11, k + 1 − hin−+11, k ) + ⎨ i ⎬ + 2 2 ⎪ ⎪ t Δ z z Δ Δ ⎝ ⎠ ⎩ ⎭

qz = z = qB = 0

⎛ K n + 1, k ⎞ i + 1/2 ⎟(h n + 1, k + 1 − hin++11, k ) (hin + 1, k + 1 − hin + 1, k ) − ⎜⎜ 2 ⎟ i+1 Δ ⎝ z ⎠ ⎛ θ n − θ n + 1, k ⎞ ⎛ K n + 1, kω2 ⎞ i ⎟⎟ + ⎜⎜ i + 1/2 ⎟⎟(r0 − zi + 1/2) = ⎜⎜ i t Δ ⎠ ⎝ g Δz ⎠ ⎝ ⎛ K n + 1, kω2 ⎞ ⎛ K n + 1, k ⎞ i − 1/2 ⎟⎟(r0 − zi − 1/2) + ⎜⎜ i − 1/2 ⎟⎟hin−+11, k − ⎜⎜ 2 ⎝ Δz ⎠ ⎝ g Δz ⎠ ⎛ K n + 1, k + K n + 1, k ⎞ ⎛ K n + 1, k ⎞ i − 1/2 i + 1/2 ⎟⎟hin + 1, k + ⎜⎜ i + 1/2 ⎟⎟hin++11, k − ⎜⎜ 2 Δz 2 ⎝ ⎠ ⎝ Δz ⎠

After discretization, the boundary condition is given by the following discretized equation: k ⎫ ⎛ K n + 1, k ⎞ ⎧ C n + 1, k KNn+−1,1/2 ⎜⎜ − N − 1/2 ⎟⎟(hNn+−1,1 k + 1 − hNn+−1,1 k ) + ⎨ N ⎬ + 2 2 ⎪ ⎪ Δz ⎠ Δz ⎭ ⎝ ⎩ Δt ⎪

⎛ 2 ⎞⎛ K n + 1, k ⎞ ⎛ θ n − θ n + 1, k ⎞ q ω N ⎟⎟ − B − ⎜⎜ ⎟⎟⎜⎜ N − 1/2 ⎟⎟(r0 − zN − 1/2) = ⎜⎜ N Δt ⎝ ⎠ Δz ⎝ g ⎠⎝ Δz ⎠ ⎛ K n + 1, k ⎞ N − 1/2 ⎟(hn + 1, k − hNn+ 1, k ) + ⎜⎜ 2 ⎟ N−1 ⎝ Δz ⎠

(13)

Aδ k = b

⎛ K n + 1, k ⎞ i + 1/2 ⎟ Ai , i + 1 = − ⎜⎜ 2 ⎟ ⎝ Δz ⎠

(17b)

where r0 is the radial position at the filter (z = 0), Hc is the height of the porous sand cake, and Hliq is the height of the liquid layer above the cake. The radial position at the liquid surface increases as the liquid layer thickness above the cake decreases. The liquid layer height is updated at each time step using the flow rate at the interface between the liquid layer and the cake as follows:

A1,1 = 1,

⎡⎛ hn + 1 − hn + 1 ⎞ ⎛ θ n − θ n+1 ⎞ N N−1 ⎟ ⎟ − KNn+−11/2⎢⎜ N qB = Δz⎜ N ⎢⎣⎝ Δt Δz ⎠ ⎝ ⎠ +

⎤ ω2 (r0 − zN − 1/2)⎥ ⎥⎦ g

AN , N

Hliq(t ) =

∫0

(19)

Vliquid layer(t ) πD 2 4

( )

⎧ K n + 1, m ⎪− N − 1/2 , if Hliq = 0 =⎨ Δz 2 ⎪ ⎩ 0, if Hliq > 0

⎧⎛ n ⎛ 2 ⎞⎛ K n + 1, k ⎞ n + 1, k ⎞ q ω ⎪ ⎜ θN − θN ⎟⎟ − B − ⎜⎜ ⎟⎟⎜⎜ N − 1/2 ⎟⎟ ⎜ ⎪⎝ z Δt Δ ⎠ ⎝ g ⎠⎝ Δz ⎠ ⎪ n 1, k ⎞ + ⎛ ⎪ KN − 1/2 ⎟(hn + 1, k − hNn+ 1, k ), bN = ⎨ (r0 − zN − 1/2) + ⎜⎜ 2 ⎟ N−1 ⎪ ⎝ Δz ⎠ ⎪ ⎪ if Hliq = 0 ⎪ ⎩ 0, if Hliq > 0

t

qB(τ ) dτ

b1 = 0

m ⎧ C n + 1, m KNn+−1,1/2 ⎪ N + , if Hliq = 0 = ⎨ Δt Δz 2 ⎪ ⎩1, if Hliq > 0

AN , N − 1

(18)

The volume and height of the liquid layer are obtained from the following equations:

⎛ πD 2 ⎞ Vliquid layer(t ) = Vliquid layer(t = 0) − ⎜ ⎟ ⎝ 4 ⎠

(24)

⎛ θ n − θ n + 1, k ⎞ ⎛ K n + 1, kω2 ⎞ i ⎟⎟ + ⎜⎜ i + 1/2 ⎟⎟(r0 − zi + 1/2) bi = ⎜⎜ i Δt ⎠ ⎝ g Δz ⎠ ⎝ ⎛ K n + 1, kω2 ⎞ ⎛ K n + 1, k ⎞ i − 1/2 ⎟⎟(r − zi − 1/2) + ⎜⎜ i − 1/2 ⎟⎟hin−+11, k − ⎜⎜ 2 0 ⎝ g Δz ⎠ ⎝ Δz ⎠ ⎛ K n + 1, k ⎞ ⎛ K n + 1, k + K n + 1, k ⎞ i − 1/2 i + 1/2 ⎟⎟hin + 1, k + ⎜ i + 1/2 ⎟hin++11, k − ⎜⎜ ⎜ Δz 2 ⎟ Δz 2 ⎝ ⎠ ⎝ ⎠

Rcake and Rliq are radial positions at the cake and liquid layer surfaces, respectively. They are given by the following equations:

R liq = r0 − Hc − Hliq

⎪

⎛ K n + 1, k ⎞ i − 1/2 ⎟⎟ Ai , i − 1 = ⎜⎜ − z2 ⎠ Δ ⎝

(16)

(17a)

k

⎪

and the pressure head at node N is given by the following equation:

R cake = r0 − Hc

k+1

1, k 1, k ⎫ ⎧ C n + 1, k K in−+1/2 + K in++1/2 i ⎬ Ai , i = ⎨ + ⎪ ⎪ Δz 2 ⎩ Δt ⎭

(15)

P(z = zN ) ρg

(23)

−h. where δ = h The matrix A and the vector b are defined as follows: For I = 2, N − 1 k

At the top surface of the cake: z = zN. The boundary condition applied depends upon the extraction phase. For the case of liquid layer filtration through the cake (Figure 1a): The pressure head at the surface of the cake is caused by the liquid layer pressure applied above the cake. The pressure applied by the liquid layer on the surface of the cake is approximated to be equal to the pressure at node N

hN (z = zN ) =

(22)

In matrix format, the discretization equations can be summarized in the following equation:

(14)

1 P(z = zN ) = ρω2(R cake 2 − R liq 2) 2

⎪

(hNn+ 1, k + 1 − hNn+ 1, k )

The following boundary conditions are applied: At the bottom surface of the cake: z = z1 = 0. The hydraulic resistance of the filter paper is assumed negligible.8 Thus, liquid pressure at the bottom of the cake is atmospheric, and the pressure head is zero: h(z = 0) = h1 = 0

(21)

N

(20)

where D is the diameter of the cylindrical cake. D

DOI: 10.1021/acs.energyfuels.6b03252 Energy Fuels XXXX, XXX, XXX−XXX

Energy & Fuels

Article

The analytical solution to eq 31 is given as follows:

It is found through initial numerical tests that the mass-conservative numerical scheme developed by Celia et al. for the solution of the variably saturated flow equation lacks robustness when used to solve the Richards equation in a centrifugal field. In this work, the combination of the modified Picard iteration and the Anderson acceleration algorithm12−14 are used to solve the Richards equation in a centrifugal field. Both stages of diluted bitumen extraction, filtration and desaturation, are modeled. 2.2. Analytical Solution of Liquid Filtration through Saturated Cake. In the section above, the numerical discretization used for the solution of liquid flow through porous cake for both filtration and desaturation phases was presented. In this section, an analytical solution is presented for the case of centrifugal liquid filtration through an incompressible saturated cake. This analytical solution can be used in fitting experimental data to obtain the saturated hydraulic conductivity, Ks. The problem of centrifugal filtration through a porous cake saturated with incompressible liquid is given by the following equation and boundary conditions: ω2(r0 − z) ⎞⎞ ∂ ⎛ ⎛ dh ⎜⎜K s⎜ ⎟⎟⎟ = 0 + ∂z ⎝ ⎝ dz g ⎠⎠

Hliq(t ) =

σ1 =

⎛ ⎛ ω2 ⎞⎛ z 2 ⎞ ω2 ⎜⎛ 1 ⎞⎡⎢ Hliq(t ) h(z , t ) = ⎜ ⎟⎜ ⎟ + (2R cake − Hliq(t )) ⎜ ⎟ g ⎜⎝⎝ Hc ⎠⎢⎣ 2 ⎝ g ⎠⎝ 2 ⎠

(28)

The flow rate per unit surface area at the cake surface is given by the following equation:

(29)

The flow rate per unit surface area at the cake surface also corresponds to the rate of liquid level decrease with time.

dHliq dt

(30)

After eqs 28, 29 and 30 are combined, the following differential equation is obtained:

dHliq dt

+ AHliq(t )(B − Hliq(t )) = C

θ(h) − θr 1 = (1 + |αh|n )m θs − θr

(35) (36)

variable

value

modeling time (min) rotation speed (rpm) cake thickness, Hc (mm) liquid layer above the cake, Hliq (mm) radial position at the filter paper, r0 (mm) cake diameter, Dcake (mm) saturated liquid content, θs, or porosity residual liquid content, θr van Genuchten parameter α (m−1) van Genuchten parameter n van Genuchten parameter l sand particle diameter (μm)

60 2000 10 5 61.7 48.5 0.28 0.047 3.28 2.03 0.5 100

Genuchten soil hydraulic parameters presented in Table 1 were reported by Simunek and Nimmo8 for a soil sample made up of a packed sand−water system. In the present paper, sand−water hydraulic parameters are used to describe a system of sand and toluene-diluted bitumen for the purpose of using the numerical model to predict the transition time from filtration to desaturation and describe qualitatively the effects of centrifugal force and solvent fraction on the extraction process. In a companion paper (10.1021/ acs.energyfuels.6b03253),19 the hydraulic properties of solvent-diluted bitumen−sand system are obtained from experimental data. The saturated hydraulic conductivity of the cake was estimated using the Kozeny−Carman equation

(31)

where ⎧ K ω2 ⎪A = s 2gHc ⎪ ⎪ ⎨ B = 2R cake ⎪ K sω2 ⎛ H 2⎞ ⎪ ⎜r0Hc − c ⎟ ⎪C = − gHc ⎝ 2 ⎠ ⎩

(34)

Table 1. Parameters Used in the Batch Filtration and Desaturation Model

2 ⎞⎤

q(z = Hc) =

A2 B2 4

where Se(h) is the effective liquid content as a function of the pressure head h. K(θ) is the hydraulic conductivity function of the volumetric liquid content θ, which is the volume fraction of liquid in the porous cake and can be estimated by θ = Sϕ, where ϕ is the cake porosity and S is the pore saturation that varies from 0 to 1. θr and θs are residual and saturated liquid contents, respectively. l is a pore connectivity parameter assumed to be 0.5.16 α, n, and m = 1 − 1/n are empirical parameters, where α is related to the mean pore size of the porous matrix and n is related to the non-uniformity of the pore size distribution.17 The numerical model is used to predict filtration and desaturation of an oil sand cake. Density and viscosity values for Athabasca bitumen diluted with toluene are obtained from Guan.18 Operational and geometric parameters used in the model are presented in Table 1. van

where H0 is the initial height of the liquid layer on the cake. The pressure at the cake surface, P(t, z = Hc), is given by eq 15. The solution of eq 25 gives the head distribution in the porous medium as follows:

c

(33)

2A

K (θ) = K sSel[1 − (1 − Se1/ m)m ]2

(27)

⎡ ∂h ⎤ ω2 q(z = Hc) = − K s⎢ (r0 − z)⎥ + g ⎣ ∂z ⎦z = H

AC −

Se(h) =

with the initial condition

⎞ ⎛ H + ⎜r0Hc − c ⎟⎥ − r0⎟⎟z 2 ⎠⎥⎦ ⎝ ⎠

σ1

Equation 33 represents time variations of liquid layer thickness above the cake during centrifugal filtration. 2.3. Hydraulic Properties and Model Parameters. The hydraulic properties of the sand cake with variable diluted bitumen saturation are assumed to be described by the following van Genuchten−Mualem equations:15

(26)

Hliq(t = 0) = H0

⎞ σ1⎟⎟σ1 + AB ⎟ ⎠ ⎠

⎛ AB − 2AH0 ⎞ ⎞ ⎟ atan⎜ ⎝ 2σ1 ⎠ ⎟

where

(25)

⎧ h(t , z = 0) = 0 ⎪ ⎨ P(t , z = Hc) ⎪ h(t , z = L) = ρg ⎩

⎛⎛ 2 tan⎜⎜⎜t − ⎜ ⎝⎝

(32) E

DOI: 10.1021/acs.energyfuels.6b03252 Energy Fuels XXXX, XXX, XXX−XXX

Energy & Fuels ⎛ ρg ⎞ d p2θs 3 K=⎜ ⎟ ⎝ μ ⎠ 180(1 − θs)2

Article

(37)

where dp is the sand particle diameter and θs is the porosity of the cake.

3. RESULTS AND DISCUSSION 3.1. Centrifugal Filtration and Desaturation. Results are first presented for the case of toluene-diluted bitumen with toluene mass fraction of 0.6 extracted at a rotation speed of 2000 rpm. Effects of rotation speed and solvent fraction on extraction efficiency are presented in the subsequent sections. An example of the evolution of the residual error with iteration number at different modeling times is presented in Figure 4.

Figure 5. (Left y axis) Evolution with time of liquid layer thickness above the cake for toluene mass fraction of 0.6 and 2000 rpm: (dots) numerical solution and (line) analytical solution. (Right y axis) Evolution with time of flow rate for toluene mass fraction of 0.6 and 2000 rpm. The vertical dashed line corresponds to the time at the transition from filtration to desaturation.

filtration to desaturation corresponds simultaneously to the time at which the liquid layer thickness above the cake recedes to zero and to a significant slope change in the logarithmic plot of flow rate versus time. The extracted liquid fraction at any time t is calculated using the following equation: t ⎧ 1 qfiltration(t ) dt t < t0 ⎪ ⎪ Vliq 0 x(t ) = ⎨ V θ (1 − S ̅ (t )) ⎪ Vliq layer t ≥ t0 + cake s ⎪ V Vliq ⎩ liq

∫

Figure 4. Convergence behavior of the numerical scheme for toluene mass fraction of 0.6 and rotation speed of 2000 rpm at different times indicated in the legend.

The residual at the kth iteration is defined as the 2-norm of the vector δk. residual = || δ k ||2

(39)

where Vliq and Vliq layer are the initial sample total liquid volume and liquid volume layer above the cake, respectively. t0, is the transition time from filtration to desaturation. S̅(t) is the average cake saturation at time t; it is calculated by averaging the saturation distribution through the cake thickness as follows:

(38)

Figure 4 indicates the robustness of the numerical scheme used for solving the Richards equation in a centrifugal field. Figure 5, with y axis on the left, shows the evolution with time of the liquid layer above the cake during centrifugation. Dots in this figure represent results of the numerical model, where the line represents the analytical solution. It is clear that good agreement is obtained between the analytical (eq 33) and numerical solutions during filtration. The vertical dashed line represents the transition time from filtration to desaturation, which is 2.4 s in this case. It corresponds to the time where the liquid level above the cake drops to zero. Figure 5, with y axis on the right, shows that liquid flow rate evolution with time can be described in two stages: filtration and desaturation. In the filtration stage, the flow rate decreases slowly with time as the liquid layer above the saturated cake recedes. On a logarithmic scale, the filtration stage is presented as a nearly horizontal flow rate line. In the desaturation stage that follows, the flow rate decreases by orders of magnitude because of the significant decrease in hydraulic conductivity with pore saturation. The transition between filtration and desaturation clearly coincides with the significant change in the slope of the flow rate variations plotted with time. Thus, the transition time from

S ̅ (t ) =

1 Hc

∫0

Hc

θ (z , t ) dz θs

(40)

Figure 6, with y axis on the right, shows the evolution with time of the extracted liquid fraction on a logarithmic scale. This figure shows that the transition from filtration to desaturation corresponds to the inflection point in the plot of extracted liquid volume fraction versus log time. The slope of this plot corresponds to the product (flow rate × time). This product can be interpreted as the volume of liquid extracted if the flow rate at time t was kept constant from time t = 0 to time t. The evolution with time of the product (flow rate × time) is presented in Figure 6, with y axis on the left, where the transition time from filtration to desaturation coincides with the maximum in the plot of (flow rate × time) versus time. The evolution of the pressure head distribution with time is presented in panels a and b of Figure 7. Figure 7a presents the pressure head contours in the coordinates distance z and log(time), whereas Figure 7b presents the time evolution of the F

DOI: 10.1021/acs.energyfuels.6b03252 Energy Fuels XXXX, XXX, XXX−XXX

Energy & Fuels

Article

Figure 6. (Left axis) Product (flow rate × time) for toluene mass fraction of 0.6 and 2000 rpm. The vertical dashed line corresponds to the time at the transition from filtration to desaturation. (Right axis) Fraction of extracted liquid on a logarithmic time scale.

pressure head profiles at different times. It is clear that, with increasing time, the pressure head at the surface decreases from positive to negative values to reflect the transition from filtration to desaturation. The pressure head at the cake surface is positive when a liquid layer applies centrifugal pressure above the surface of the cake; it becomes negative when the liquid layer thickness vanishes to zero and cake desaturation is initiated. The pressure head profile converges with time to an equilibrium profile represented by a parabolic decrease from the bottom to the top of the cake. Panels a and b of Figure 8 show the evolution with time of the pore saturation profile across the cake thickness. Saturation contours are presented in Figure 8a in the coordinates distance z and log(time) and the evolution of saturation profiles with time is presented in Figure 8b. Pore saturation remains at 1.00 until the transition from filtration to desaturation when the porous cake is desaturated progressively from the top surface. Pore saturation profiles converge with time to a final equilibrium profile. Figure 9 shows the evolution with time of free surface and average pore saturation in the cake. The average saturation is higher than the saturation at the free surface because saturation increases across the cake thickness and reaches 1 at the bottom of the cake. The pore saturation profile stabilizes with time at equilibrium. 3.2. Effect of Rotation Speed and Cenrifugal Force. For the case of oil sand mixed with toluene with solvent fraction of 0.6 in diluted bitumen, the effect of centrifugal force is investigated by increasing rotation speed from 1000 to 6000 rpm. This rotation speed range corresponds to a centrifugal force range from 63.4g to 2284.2g. For industrial conical filtering centrifuges, used for continuous solid−liquid separation, the maximum centrifugal force at the basket lip varies in the range from 2250g to 2890g.20 Figure 10 shows the effect of rotation speed on the extracted liquid flow rate. Increasing rotation speed increases flow rates during the filtration phase and accelerates the transition from filtration to desaturation. During desaturation, the flow rate depends upon both rotation speed and saturation distribution. At any given time, local cake saturation decreases with increasing rotation speed. This causes,

Figure 7. Pressure head distribution in the cake, h(m), for toluene mass fraction of 0.6 and rotation speed of 2000 rpm: (a) pressure head contour versus time and cake thickness and (b) variations with time of the head profile across cake thickness.

over a time period, desaturation flow rates to be lower at higher rotation speed. Figure 11 shows the extracted liquid fractions as a function of the extraction time at different centrifugal speeds. This figure shows that the transition time from filtration to desaturation, corresponding to the inflection point, is shortened with increasing centrifugal force. The fraction of extracted liquid increases with increasing centrifugal force at any given time, and a higher level of cake desaturation is obtained at the final time. 3.3. Effect of the Solvent Mass Fraction. The effect of the toluene mass fraction in diluted bitumen is investigated in this section for a rotation speed of 2000 rpm. The initial height of the diluted bitumen layer above the porous cake is selected to be 5 mm for all cases of the solvent fraction investigated. The purpose of this calculation is to show the effect of diluted bitumen properties on filtration and desaturation. In the following results, effects of the solvent fraction on viscosity and density of diluted bitumen were taken into account. These G

DOI: 10.1021/acs.energyfuels.6b03252 Energy Fuels XXXX, XXX, XXX−XXX

Energy & Fuels

Article

Figure 9. Evolution of pore saturation with time for toluene mass fraction of 0.6 and rotation speed of 2000 rpm: (continuous line) average saturation and (dashed line) surface saturation.

Figure 8. Pore saturation in the cake for toluene mass fraction of 0.6 and rotation speed of 2000 rpm: (a) pore saturation contour versus time and cake thickness and (b) variations with time of the pore saturation profile across cake thickness.

Figure 10. Flow rate of extracted liquid by centrifugation for different rotation speeds and toluene mass fraction of 0.6: (continuous line) 1000 rpm (63.4g), (dashed line) 2000 rpm (253.8g), (dash-dotted line) 3000 rpm (571.0g), (thick dotted line) 4000 rpm (1015.2g), (thick continuous line) 5000 rpm (1586.2g), and (thick dashed line) 6000 rpm (2284.2g).

properties affect both filtration and desaturation. van Genuchten parameters were kept constant as a result of the lack of prior data of the diluted bitumen−oil sand system. These hydraulic parameters control the liquid retention curve and hydraulic conductivity function (eqs 35 and 36), and they are expected to have a dependence upon the solvent fraction. Figure 12 shows the effect of the toluene fraction on the time evolution of the extracted liquid flow rate. Flow rates increase during the filtration phase, and the transition from filtration to desaturation is accelerated with an increasing toluene fraction. This can be explained by the fact that the flow rate is significantly influenced by viscosity. The latter decreases by 3 orders of magnitude as the toluene mass fraction is increased from 0.1 to 0.6. When desaturation is initiated, flow rates decrease by orders of magnitude. In this case, desaturation flow rates also depend upon pore saturation. At any given time during desaturation, flow rates are lower for a higher toluene fraction as a result of lower average cake pore saturation.

Figure 13 shows the evolution with time of the extracted liquid volume fraction for an increasing toluene mass fraction from 0.1 to 0.6. The higher the mass fraction of toluene, the faster the equilibrium value of the extracted liquid fraction is reached. However, final equilibrium values of the extracted liquid fraction appear to converge with an increasing toluene mass fraction. This trend was also observed in experimental data.3 In the calculation performed in this work, this result is expected because van Genuchten parameters are assumed to be independent of the solvent mass fraction. These parameters describe the effect of capillary and adsorptive forces through the liquid retention curve and the hydraulic conductivity function. Equilibrium pore saturation profiles depend upon balance between capillary and adsorptive forces and centrifugal forces and not on viscous forces. H

DOI: 10.1021/acs.energyfuels.6b03252 Energy Fuels XXXX, XXX, XXX−XXX

Energy & Fuels

Article

Figure 11. Total extracted liquid fraction for different rotation speeds: (continuous line) 1000 rpm (63.4g), (dashed line) 2000 rpm (253.8g), (dash-dotted line) 3000 rpm (571.0g), (thick dotted line) 4000 rpm (1015.2g), (thick continuous line) 5000 rpm (1586.2g), and (thick dashed line) 6000 rpm (2284.2g).

Figure 13. Predicted evolution of (a) extracted liquid volume and (b) fraction of extracted liquid volume for rotation speed of 2000 rpm (253.8g) and different toluene mass fractions.

hydraulic conductivity, Ks, can be determined by fitting experimental data of the extracted liquid volume during the filtration phase. van Genuchten parameters, θr, θs, α, and n, of diluted bitumen flow in the porous cake can be determined from experimental data of the extracted liquid volume during the desaturation phase. The hydraulic parameters can be determined by minimization of the objective function defined as the sum of the squared deviations between measured and calculated liquid volumes in the filtrate reservoir. Using experimental volume data of the filtration phase, the saturated hydraulic conductivity can be calculated by minimizing the following objective function: obj(K s) = m

njfilt

experimental predicted (ti , ωj) − Vfiltrate (ti , ωj , K s)]2 ∑ ∑ [V filtrate j=1 i=1

(41)

Experimental volume data of the desaturation phase can be used to obtain van Genuchten hydraulic parameters by minimizing the following objective function:

Figure 12. Predicted evolution of the diluted bitumen flow rate from the cake to the reservoir for rotation speed of 2000 rpm (253.8g) and different toluene mass fractions.

obj(param) =

Oil sand ore grade has an impact on both the liquid retention curve and hydraulic conductivity function. These functions have a significant effect on the fraction of liquid recovered by centrifugation. For the case of water, the relationship between the particle size distribution and the liquid retention curve has been investigated by several authors.21,22 In this work, van Genuchten parameters, which characterize these functions, were kept constant. Experimental investigations for different grades of oil sand should be performed to characterize the effect of ore grade on the hydraulic properties of the diluted bitumen−oil sand cake system. 3.4. Parameter Optimization. The model described in this work can be used to predict the flow rate and filtrate volume during batch centrifugal filtration of solvent-diluted bitumen from mineral particles. Given experimental filtrate volume data, a criterion for the determination of the transition time from filtration to desaturation was presented. Saturated

desat m nj

experimental predicted (ti , ωj) − Vfiltrate (ti , ωj , param)]2 ∑ ∑ [V filtrate j=1 i=1

(42)

where m is the number of constant rotational speed desat experiments. nfilt are the number of data measured j and nj at a particular rotational speed during filtration and desaturation phases, respectively. param is the parameter vector (θr, θs, α, and n).

4. CONCLUSION The numerical solution of the Richards equation for flow through a porous medium with variable saturation was developed. The numerical method was based on the massconservative-modified Picard iteration scheme and Anderson acceleration algorithm. This numerical scheme proved to be I

DOI: 10.1021/acs.energyfuels.6b03252 Energy Fuels XXXX, XXX, XXX−XXX

Energy & Fuels

Article

International Symposium of the Society of Core Analysts; Calgary, Alberta, Canada, Sept 10−13, 2007; SCA 2007-22. (7) Ruth, D. Analysis of Centrifuge Relative Permeability Data. Proceedings of the International Symposium of the Society of Core Analysts; Calgary, Alberta, Canada, Sept 7−10, 1997; SCA 9711. (8) Šimůnek, J.; Nimmo, J.R. Estimating soil hydraulic parameters from transient flow experiments in a centrifuge using parameter optimization technique. Water Resour. Res. 2005, 41 (4), W04015. (9) Bear, J.; Yavuz Corapcioglu, M.; Balakrishna, J. Modeling of centrifugal filtration in unsaturated deformable porous media. Adv. Water Resour. 1984, 7, 150−167. (10) Nimmo, J. R. Experimental Testing of Transient Unsaturated Flow Theory at Low Water Content in a Centrifugal Field. Water Resour. Res. 1990, 26 (9), 1951−1960. (11) Celia, M. A.; Bouloutas, E. T.; Zarba, R. L. A General Mass Conservative Numerical Solution for the Unsaturated Flow Equation. Water Resour. Res. 1990, 26 (7), 1483−1496. (12) Walker, H. F. Anderson Acceleration: Algorithms and Implementations; Mathematical Sciences Department, Worcester Polytechnic Institute: Worcester, MA, 2011; Research Report MS-615-50. (13) Walker, H. F. W.; Carol, S.; Yang, U. M. An Accelerated FixedPoint Iteration for Solution of Variably Saturated Flow. Proceedings of the XVIII International Conference on Water Resources; Barcelona, Spain, June 21−24, 2010. (14) Lott, P. A.; Walker, H. F.; Woodward, C. S.; Yang, U. M. An accelerated Picard method for nonlinear systems related to variably saturated flow. Adv. Water Resour. 2012, 38, 92−101. (15) van Genuchten, M. Th. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 1980, 44, 892−898. (16) Mualem, Y. A New Model for Predicting the Hydraulic Conductivity of Unsaturated Porous Media. Water Resour. Res. 1976, 12 (3), 513−522. (17) Miller, C. T.; Williams, G. A.; Kelley, C. T.; Tocci, M. D. Robust solution of Richards’ equation for nonuniform porous media. Water Resour. Res. 1998, 34 (10), 2599−2610. (18) Guan, J. Physical Properties of Athabasca Bitumen and Liquid Solvent Mixtures. M.S. Thesis, Chemical and Petroleum Engineering, University of Calgary, Calgary, Alberta, Canada, Jan 2013. (19) Khammar, M.; Xu, Y. Batch Solvent Extraction of Bitumen from Oil Sand. Part 2: Experimental Development and Modeling. Energy Fuels 2017, DOI: 10.1021/acs.energyfuels.6b03253. (20) Wakeman, R.; Tarleton, S. Solid−Liquid Separation: Scale-up of Industrial Equipment; Elsevier: Amsterdam, Netherlands, 2005. (21) Fredlund, M. D.; Wilson, G. W.; Fredlund, D. G. Use of the grain-size distribution for estimation of the soil−water characteristic curve. Can. Geotech. J. 2002, 39 (5), 1103−1117. (22) Simms, P. H.; Yanful, E. K. Predicting soil−water characteristic curves of compacted plastic soils from measured pore-size distributions. Geotechnique 2002, 52 (4), 269−278.

robust in all cases investigated. The numerical method can be used to solve for flow through oil sand cake during filtration and desaturation phases of the centrifugal solvent extraction process. The predicted filtrate liquid flow rate plotted on a log− log scale shows clear delineation between the filtration and desaturation phases of the extraction process. It was found that the transition time from filtration to desaturation corresponds to the inflection point in the plot of extracted liquid volume versus log(time). Thus, the transition time corresponds to the time at the maximum of the plot of (flow rate × time) versus time. The effects of the rotation speed and toluene mass fraction in diluted bitumen on the extracted liquid fraction were investigated. It was found that diluted bitumen extraction is accelerated when the toluene mass fraction and/or centrifugal force are increased. Increasing the centrifugal force also increases the final fraction of diluted bitumen extracted. Evolutions of the extracted liquid fraction with the centrifugal force and solvent fraction are qualitatively similar to those obtained experimentally for batch centrifugal extraction of solvent-diluted bitumen from oil sands. A fitting method is proposed for the determination of the porous cake hydraulic properties based on experimental data of extracted liquid volumes. Saturated hydraulic conductivity of the cake can be obtained by fitting experimental liquid volume data during the filtration phase, and van Genuchten parameters for the unsaturated cake can be obtained by fitting experimental liquid volume data during the desaturation phase.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Merouane Khammar: 0000-0003-0077-931X Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge the late distinguished scientist Dr. Tadek Dabros, who made significant contributions to oil sands research and development. The authors also acknowledge support from the Government of Canada’s Interdepartmental Program of Energy Research and Development (PERD).

■

REFERENCES

(1) Faradonbeh, M. R.; Dong, M.; Harding, T. G.; Abedi, J. Experimental and modeling study of residual liquid recovery from spent sand in bitumen extraction processes from oil sands. Environ. Sci. Technol. 2013, 47 (4), 2109−2116. (2) Meadus, F. W.; Bassaw, B. P.; Sparks, B. D. Solvent extraction of athabasca oil-sand in a rotating mill Part 2. Solidsliquid separation and bitumen quality. Fuel Process. Technol. 1982, 6, 289−300. (3) Wu, J.; Dabros, T. Process for Solvent Extraction of Bitumen from Oil Sand. Energy Fuels 2012, 26 (2), 1002−1008. (4) Hillel, D. Introduction to Environmental Soil Physics; Elsevier: Amsterdam, Netherlands, 2004. (5) Forbes, P. Centrifuge data analysis techniques: An SCA survey on the calculation of drainage capillary pressure curves from centrifuge measurements. Proceedings of the International Symposium of the Society of Core Analysts; Calgary, Alberta, Canada, Sept 7−10, 1997; SCA 9714. (6) Fernø, M. A.; Treinen, R.; Graue, A. Experimental measurements of capillary pressure with the centrifuge techniqueEmphasis on equilibrium time and accuracy in production. Proceedings of the J

DOI: 10.1021/acs.energyfuels.6b03252 Energy Fuels XXXX, XXX, XXX−XXX