A finite volume method for density driven flows in porous media

ESAIM: PROCEEDINGS, December 2012, Vol. 38, p. 376-386 F. Coquel, M. Gutnic, P. Helluy, F. Lagoutière, C. Rohde, N. Seguin, Editors A FINITE VOLUME METHOD FOR DENSITY DRIVEN FLOWS IN POROUS MEDIA Danielle Hilhorst 1 , Huy Cuong Vu Do 2 and Yushan Wang 3 Abstract. In this paper, we apply a semi-implicit finite volume method for the numerical simulation of density driven flows in porous media; this amounts to solving a nonlinear convection-diffusion parabolic equation for the concentration coupled with an elliptic equation for the pressure. We compute the solutions for two specific problems: a problem involving a rotating interface between salt and fresh water and the classical but difficult Henry’s problem. All solutions are compared to results obtained by running FEflow, a commercial software package for the simulation of groundwater flow, mass and heat transfer in porous media. Introduction We describe here results which have been obtained in the context of an exploratory project of CNRS (PEPS ECODEVA) on the numerical simulation of flows with variable density for the production of lithium batteries. More precisely, the purpose of this project is related to the exploitation of lithium deposits in salt lakes, also known as ”Salars”. In recent years, lithium has become a strategic element for industrial countries because it is the basic element of lithium-ion batteries used for hybrid and electric vehicles. Therefore its production has become of high interest for all major groups involved in the car industry as well as suppliers of these groups. Currently the largest deposit in the world is the Salar Uyuni, in the department of Potosı́, in South-West Bolivia. This deposit represents a third of the world resources. In March 2008, Bolivia has authorized the exploitation, however reserving this right to its nationals. Chile has the second largest deposit with the Salar Atacama and it has become the world’s largest exporter since 1997, with the German company Chemettal as main operator. Argentina also has a lithium deposit, the Salar Hombre Muerto, which is located in the North-West of the country. Other salar areas of the Altiplano of Argentina provide mining exploitation concessions to foreign companies, among whom European groups. Other deposits are exploited, including salt lakes in Tibet as well as mines in Australia, Russia and the United States. They are not accessible to European operators. The largest deposits are either clusters of crystallized salt (solid) or lenses of supersaturated salt water created by evaporation under endorheic conditions (which are not led to a superficial network reaching the sea). The latter type of deposit is that of salars of the Andean altiplano. Rational exploitation implies mastering these special aqueous flows whose density depends on the concentration of salts (lithium included). An operating technique consists in sweeping the reservoir with fresh 1 Laboratoire de Mathématiques, CNRS et Université de Paris-Sud (Bâtiment 425) 2 Early Stage Researcher of the European Marie Curie ITN Project FIRST, Laboratoire de Mathématiques, Université de Paris-Sud (Bâtiment 425) 3 Laboratoire de Recherche en Informatique, Université de Paris-Sud (Bâtiment 650) EDP Sciences, SMAI 2012 © Article published online by EDP Sciences and available at http://www.esaim-proc.org or http://dx.doi.org/10.1051/proc/201238021 377 ESAIM: PROCEEDINGS water in order to obtain a maximal recovery without earthworks and with a minimal impact on fluid levels, and thus a minimal impact on the environment. This explains the need of implementing research methodologies and techniques from hydrogeology. The purpose here is to extract salt water which contains lithium. In a later stage, the lithium will be separated from the salt water. From a mathematical viewpoint, it amounts to study a coupled system describing the interaction between flow and transport in a porous medium, whose density ρ is a strictly increasing function of the concentration c of the salt water; typically we have that ρ(c) = ρ0 (1 + αc). More specifically, we solve a nonlinear parabolic convection-diffusion equation for the concentration coupled with an elliptic equation for the pressure, which is derived from Darcy’s law    ∂(ρ(c)c)  + ∇ · qρ(c)c − ρ(c)D∇c = 0, Φ    ∂t  ∂ρ(c) + ∇ · (qρ(c)) = 0, Φ  ∂t   k  q = − (∇p − ρ(c)g), µ in Ω × (0, T ], in Ω × (0, T ], (1) in Ω × (0, T ], together with suitable boundary conditions. The porosity Φ is the fraction of the voids (empty spaces) over the total volume. In the third equation of (1), namely, Darcy’s law, q is the flux (discharge per unit area), p is the pressure, k is the permeability, µ is the dynamic viscosity, and g is the gravity. The organization of this paper is as follows. In Section 1, we present the numerical algorithm which is applied to solve System (1). It is based upon the standard finite volume method for the spatial discretization and a semi-implicit scheme for the time discretization [9], [6], [7], [8]. Because of the form of the velocity vector q, we simultaneously use two upwind directions for the discretization of the convection terms. In Section 2, we describe Henry’s problem, together with suitable boundary conditions. The main difficulty is that we have to compute pressure and velocity, in a case that the velocity changes direction through one of the boundaries. In Section 3, we study the time evolution of the interface between fresh and salt water in an aquifer in a case that the interface slowly rotates before reaching a stable equilibrium where the fresh water lies above the salt water. A technical difficulty is that the pressure is not uniquely defined since Neumann boundary conditions are prescribed on the boundaries. This leads us to slightly modify System (1) by adding a small parabolic term to the pressure equation. All the solutions which we compute throughout this paper are compared to results obtained by means of FEflow, a commercial software for the simulation of groundwater flow and mass and heat transfer in porous media, which is based upon the finite element method. Finally, we present some conclusions in Section 4. We refer to [5], [11], [2], [4] for modelling aspects. System (1) has already been solved in [10] by means of a vertex centered finite volume method as well as by the mixed finite element method. They directly applied their numerical method to a real case, while we perform a rather detailed study of the application of our numerical scheme to two rather intricate test cases. 1. The finite volume discretization In this section, we discretize System (1)    ∂(ρ(c)c)  Φ + ∇ · qρ(c)c − ρ(c)D∇c = 0,    ∂t  ∂ρ(c) Φ + ∇ · (qρ(c)) = 0,  ∂t    q = − k (∇p − ρ(c)g), µ in Ω × (0, T ], in Ω × (0, T ], in Ω × (0, T ], 378 ESAIM: PROCEEDINGS together with the boundary conditions  c = cD (x, t),     p = pD (x, t), ∂c  = c̄N (x, t),     ∂n q · n = q̄N (x, t), (...truncated)