The Density-Enthalpy Method Applied to Model Two–phase Darcy Flow

The Density-Enthalpy Method Applied to Model Two–phase Darcy Flow

Ibrahim, D.

Vuik, C. (promotor)

Electrical Engineering, Mathematics and Computer Science

Delft Institute of Applied Mathematics

2012-05-30

In this thesis, we use a more recent method to numerically solve two-phase fluid flow problems. The method is developed at TNO and it is presented by Arendsen et al. in [1] for spatially homogeneous systems. We will refer to this method as the densityenthalpy method (DEM) because the density-enthalpy phase diagrams play an important role in this approach. Multiphase flow occurs in numerous natural and industrial processes. These processes (or flow systems) are typically modeled by one or more sets of PDEs. In the literature, a huge variety of mathematical models for flow and transport in porous media are presented and used to simulate these processes. Many authors classify them into moving grid/free boundary methods and fixed grid methods. The method we use falls in the latter category. Although the DEM is developed for multiphase flow problems but this thesis is limited only to two–phase fluid flow of one substance (Propane). As the name indicates, density and enthalpy are our primary variables. The mathematical model for our approach consists of a mass balance, an energy balance, Darcy’s law and other thermodynamic relations. We solve the mass and energy balances for two state variables, the density (p) and the enthalpy (h). Other solution variables (such as pressure, temperature, and gas mass fraction, etc) are obtained from given p-h phase diagrams. These diagrams are obtained from thermodynamic properties of a substance (in our case the substance is Propane). For the spatial discretization, we use the finite-element method. An Euler Backward method is used for the time integration. The finite-element method is used for the spatial discretization of the system over 1D and 2D grids. This method is selected because of its ability to handle complex domain geometries. In particular, SUPG (Streamline Upwind Petrove-Galerkin) is used in the initial chapters. The use of SUPG is related to the numerical wiggles as discussed in the coming chapters. Later on, a standard Galerkin algorithm is applied where no spurious oscillations are observed. We use piecewise (bi-)linear basis functions to approximate solution variables and test functions appearing in the weak forms of the PDEs. The backward Euler method is used for the time integration of the PDEs. The use of piecewise (bi-)linear basis functions and Euler backward time integration scheme implies a numerical error of order O(?x2 + ?t), where ?x and ?t are the spatial and time steps, respectively. This is verified numerically. The 0D model works fine as shown by Arendsen [1, 2, 3]. The first attempt to solve a multidimensional flow system by using density-enthalpy phase diagrams was made by Abouhafc¸ in [12]. He used the IMEX (IMplicit-EXplicit) scheme for system linearization. For his solver, it was necessary to use very small time-steps. This also means that a large computation time is required for a relatively small process time. At the start of our thesis, the real cause behind a small time-step necessity was unknown. We started with IMEX to verify the results obtained by Abouhafc¸. Next, we applied Euler backward with the Picard iteration method but the challenge of a small timestep remained there. Later on, we show that our criterion for choosing the initial guess, as required at the start of each Picard iteration, makes it equivalent to IMEX linearization scheme, which imposes an upper bound on the time-step for convergence. This motivated us to use the Newton-Raphson method for the system linearization. In general, the Newton method is more sensitive to the initial guess. In our case, we use the system variables from the previous time-step as initial guess at the start of a fresh Newton-Raphson iteration. We have shown that this scheme successfully allows a reasonably large time-step (independent of the spatial grid size). Another challenge at the start of our thesis was to deal with solution variables exhibiting nonphysical steep gradients and sudden variations near the boundary. This problem also appeared in [12] but it was left untreated. We have shown that the use of nonhomogeneous boundary conditions makes the problem ill-posed. This is responsible for nonphysical behavior of the solution variables. We propose some guidelines to use nonzero boundary fluxes to keep the problem well-posed. From the simulation results of this fluid system, we conclude that the new method can successfully be applied to numerically solve multi-phase fluid systems. However we need to consider certain aspects regarding this approach. One issue is that density-enthalpy phase diagrams are not widely available for many multi-phase systems. There are also certain issues with non-homogenous boundary conditions (such as well-posedness and selection criterion of certain parameters). These issues are treated partially but are also included in the recommendation for future work.

http://resolver.tudelft.nl/uuid:be9b5fc5-fcd0-4d12-b0d8-a8e028b5d64b

978-94-6203-049-7

Institutional Repository

doctoral thesis

(c) 2012 Ibrahim, D.