Project Topics
Get Complete Project Material File(s) Now! »
Multiphase flow
Multiphase flow occurs when two or more distinct fluids coexisting are separated by a sharp interface. Each fluid is considered, then, as a phase with its own physical and chemical properties. Taking into account the contact angle between the phases in the presence of a solid phase (θ), the phases can be divided into wetting (θ < 90°) and non-wetting (θ > 90°). Due to the surface tension, there is a pressure change between phases defined with the capillary pressure ( ):
In the case where two fluids are miscible and no interface occurs between them, they are considered as part of the same phase and are called components. One phase can contain one or more components (e.g. a gas phase composed by CO2 and CH4) and each component can be present in one or more phases (i.e. CO2 can be present in a system as a gas or dissolved in liquid).
Mathematical description
The number of phases and components to be described in the governing equations depends on the system. This thesis is focused on applications where the natural systems are composed by two phases and two components. Three-phase flow models are applied mainly for water, oil and gas systems in the hydrocarbon industry. A complete description of the three-phase flow can be found at Chen et al. (2006) and Peaceman (1977) among others.
The mathematical description of two-phase systems can be found in the literature (Bastian, 1999; Bear, 1972; Chen et al., 2006; Garcia, 2003; Hassanizadeh and Gray, 1979; Helmig, 1997); and it is comprehensively described in articles related to multiphase flow simulation codes: MUFTE (Assteerawatt et al., 2005), DUMUx (Flemisch et al., 2011), OpenGeoSys (Goerke et al., 2011; O. Kolditz et al., 2012b), PFLOTRAN (Lu and Lichtner, 2005), TOUGH2 (Pruess et al., 1999; Pruess and Spycher, 2007), CodeBright (Vilarrasa, 2012), HYTEC (Sin et al., 2017), ECLIPSE (ECLIPSE, 2004), ELSA (Nordbotten et al., 2009, 2004b), FEHM (Pawar et al., 2005), GPRS (Cao, 2002; Jiang, 2008), IPARS (Wheeler and Wheeler, 2001), COORES (Le Gallo et al., 2006; Trenty et al., 2006), MoReS (Farajzadeh et al., 2012; Wei, 2012) and RTAFF (Sbai, 2007).
The mathematical model can vary widely in each application depending on the specific questions to be answered, the relevant processes, and the time scale and length associated to those processes. This thesis does not consider the geomechanical effects and the heat conservation in the multiphase flow. Simplifications of the multiphase flow such as vertical equilibrium models (Nordbotten and Celia, 2006) or macroscopic percolation models (Cavanagh and Haszeldine, 2014) are found to greatly reduce the computational complexity of the system. These simplified models are out of the scope of this thesis. A review of these methods can be found in Bandilla et al. (2015), Celia et al. (2015) and Celia and Nordbotten (2009).
This section starts with the description of the more general compositional model, where mass transfer between phases is taken into account. The component mass conservation system of equations and the constitutive and state relations needed to solve the system are explained. Then, some immiscible formulations are presented. In systems where the component behavior is not relevant, the multiphase flow can be represented by immiscible formulations with a much lower computational cost.
Compositional model
The mass conservation of each component ( ) in a system with phases can be written as:
where denotes the fluid phase, is the porosity of the medium ( 3 ⋅ −3), is the phase saturation ( 3 ⋅ −3), is the phase density ( ⋅ −3), is the molality of the component in the phase ( ⋅ −1), is the component molar weight ( ⋅ −1), corresponds to the phase Darcy’s flow in its multiphase form ( 3 ⋅ −2 ⋅ −1), is the phase diffusion/dispersion tensor ( 2 ⋅ −1), is the component mass source term in the phase ( ⋅ −3 ⋅ −1) and is the component interphase mass transfer term ( ⋅ −3 ⋅ −1).
The last four terms can be extended as:
where is the permeability tensor ( 2), is the relative phase permeability (−), is the phase dynamic viscosity ( ⋅ −1), is the pressure ( ) and is the gravity vector ( ⋅ −2), is the species-independent molecular diffusion coefficient ( 2 ⋅ −1), and are, respectively, the longitudinal and transverse dispersivity ( ), | | is the Euclidean norm of Darcy’s flow, ( ) is the orthogonal projection along velocity and ⊥ = − ( ), is the phase mass source term ( ⋅ −3 ⋅ −1), ( −1) is the interphase mass transfer of component and ̅ ( ⋅ −1) is the solubility limit of component in the phase.
The system of equations in (2-2) has 2 + dependent variables: 2 from the phase pressures and saturations and components molalities. The system is constrained by the phase saturation equation, the phase mass fraction and ( − 1) phase equilibrium relations:
Constitutive relations
To close the system ( − 1) capillary pressure relations needs to be taken into account. Generally, the capillary pressure is described as function of the saturation of the higher wettability (Niessner and Helmig, 2007). The most common approaches of the capillary pressure are Brooks and Corey (1964) and van Genuchten (1980). Recently, a more complex model where capillary pressure is function of the specific interfacial area apart from the saturation has been developed (Niessner and Hassanizadeh, 2008); however, it has been rarely implemented in multiphase flow models (Tatomir et al., 2015). In the multiphase flow equations implemented in this thesis the user can define whether Brooks and Corey or Van Genuchten capillary approach is used.
Additional equations relate the phase density and viscosity with its pressure, temperature and composition. The influence of pressure and temperature in the liquid density and viscosity is limited, while gas phases present much more variability. Usually, experimental relationships are used to determine the density and viscosity of the components. In this thesis, apart from the relations already implemented in COMSOL (Comsol, 2016), the properties of the CO2-brine system implemented by Pruess and Spycher (2007) in TOUGH2 are incorporated. The gas density can also be determined following the ideal gas law and the Peng-Robinson equation of state (Robinson and Peng, 1978). Other relationship could be easily added by the user to the system.
The partitioning of components between the phases is usually considered in equilibrium, which is a reasonable assumption when the time scale of mass transfer is large compared to the time scale of flow. Compositional multiphase formulations based on local equilibrium are extensively described in the literature (Firoozabadi, 1999; Lu and Lichtner, 2005; Michelsen et al., 2008; Pruess et al., 1999). In gas-liquid systems the Henry’s law and the saturation vapor pressure are used to define the mole fraction in each phase.
If component partitioning is not described by equilibrium, a kinetic expression is used for the interphase mass transfer (equation 2-6). Several approaches are described in the literature. Van Antwerp et al. (2008) extended the dual domain approach for kinetic mass transfer during air sparing from the laboratory scale (Falta, 2003, 2000) to field scale. Niessner and Hassanizadeh (2009) developed a method where the interphase mass transfer was function of the mass balance of the specific interfacial area. Other authors use a lumped mass transfer rate calculated through the non-dimensional modified Sherwood number ( ℎ) to solve the equation (Agaoglu et al., 2015; Imhoff et al., 1994; Miller et al., 1990; Powers et al., 1994, 1992; Zhang and Schwartz, 2000):
where is the kinetic mass transfer rate ( −1), 50 is the mean size of the reservoir grains ( ), is the Reynold’s number (-), is the saturation of the phase and , , are dimensionless fitting parameters. This approach is followed in the compositional multiphase formulation implemented in this thesis.
Choice of primary variables and system of equations Different linear combinations of equations and state variables can be used to diminish the stiffness of the compositional multiphase flow system (Binning and Celia, 1999; Chavent, 1981; Chen et al., 2000; Chen and Ewing, 1997a; Coats, 1980). The most classic approach is to sum over components and choose the so-called natural variables (pressure and saturation of the phases and molality of components) as dependent variables and apply a local variable switching depending on the phases present in each cell (Class et al., 2002; Coats, 1980; Lu and Lichtner, 2005; Pruess et al., 1999).
Phase disappearance is one of the major problems in multiphase flow simulation. It leads to the degeneration of the equation satisfied by the saturation, making natural variables inappropriate. Variable switching is able to cope with this process varying the dependent variables depending on the phases present in each cell. Nevertheless, this approach has some drawbacks as generates perturbations on the Newton iteration when phase changes occur and affect the structure of the Jacobian matrix (Lu et al., 2010).
Several alternative approaches are proposed in the literature where the dependent variables remain constant throughout the domain. Jaffré and Sboui (2010) use complementarity conditions. Abadpour and Panfilov (2009) allow the extension of saturation to negative values. Moreover, various combinations of dependent variables are described for two phase flow: liquid phase pressure and water mass concentration (Bourgeat et al., 2013); gas pressure and liquid pressure (Angelini et al., 2011); pressure, saturation and fugacities (Lauser et al., 2011); capillary pressure and gas pressure (Neumann et al., 2013);… Comparison between some of these methods can be found in literature (Gharbia et al., 2015; Lu et al., 2010; Masson et al., 2014; Neumann et al., 2013; Voskov and Tchelepi, 2012).
In this thesis, the compositional multiphase flow is implemented and applied to carbon capture and storage (section 3.2) using the natural variables as dependent variables. Water miscibility in the CO2 rich phase is low, on the order of one per cent (Spycher and Pruess, 2005). Therefore, water dissolution in supercritical CO2 is neglected, reducing the system of equations (2-2) to three. Thus, only one interphase mass transfer ( 2 ) has to be defined through the Sherwood number approach. No variable switching is applied and it is considered that the liquid phase remains throughout the system, even in small amounts. A new linear combination of the system of equations (2-2) is used to describe the three equations of the system:
where is the external CO2 concentration. These equations are derived from:
1. the total mass conservation equation (2-13) obtained by summing over all equations of system of equations (2-2);
2. an equation for the CO2 mass conservation (2-14) obtained by adding the equations for the CO2 component ( = 2) in both phases ( = , );
3. equation (2-15), derived by subtracting equation (2-13) multiplied by
2 2 from the aqueous CO2 mass conservation (equation (2-2) for = and = 2).
The dependent variables solved are liquid pressure ( ) in equation (2-13), gas saturation ( ) in equation (2-14) and dissolved CO2 molality ( 2 ) in equation (2-15).
The approaches followed to model the fluids density and viscosity, the relative permeability and capillary pressure are used-defined.
Immiscible formulations
In cases where no mass transfer between phases or when the behavior of the components is not relevant for the purpose of the model; the mass balance equations can be written for the bulk phase rather than for individual components. The phase mass conservation is achieved by summing equation (2-2) over components and assuming no internal component gradient (no diffusion-dispersion within the phases):
where is the phase mass source term ( ⋅ −3 ⋅ −1). The Darcy’s flow is described in equation (2-3). The system of equations (2-16) has 2 unknowns that can be solved taking into account the saturation constrain (equation2-7) and the ( − 1) capillary pressure relations (equation 2-1).
As for the compositional formulation, various dependent variables and linear combinations of equations can be selected (Aziz and Settari, 1979; Binning and Celia, 1999; Chavent, 1981; Chen et al., 2000; Hassanizadeh and Gray, 1979). Some of the formulations described in the literature and implemented during the thesis are summarized below for the two-phase flow systems.
Coupled Pressure-Saturation formulation
In this formulation, the general system of two equations (2-16) is solved without any linear combination. One equation takes the phase pressure as the dependent variable, while the other solves the phase saturation. Integrating the Darcy’s velocity and taking wetting saturation ( ) and non-wetting phase pressure ( ) as unknowns the system is written as:
where subscript and denote the wetting phase and non-wetting phase respectively and is the phase mobility ( = ⁄ ).
Similarly, the wetting pressure/non-wetting saturation ( − ) can be developed. This system of equations is strongly coupled and non-linear.
Decoupled Pressure-Saturation formulations
Dividing the equations in (2-16) by the density, the mass conservation equation is transformed into a volumetric equation. Summing, then, these equations over the two phases, the temporal derivative over saturation is removed, diminishing the coupling of the system:
Equation (2-19) is called the pressure equation and it can be solved using different pressures as dependent variable (Chen et al., 2006; Hoteit and Firoozabadi, 2008; O. Kolditz et al., 2012b). If is explicitly evaluated, the pressure can be first computed and, then, the saturation can be solved with one of the equations of system in (2-16). This is the implicit pressure-explicit saturation (IMPES) scheme that has been extensively used in two-phase flow simulations and was first introduced by Sheldon and Cardwell (1959) and Stone and Garder (1961).
This formulation with the wetting pressure/non-wetting saturation ( − ) as dependent variables is applied for the model in section 4. The multiphase flow for these unknowns has the following form:
Similar systems of equations can be developed for the − unknowns or using the wetting phase equation instead of (2-21).
Global pressure formulation
The global pressure ( ) concept was first introduced by Antontsev (1972) and Chavent and Jaffré (1986) and it has been frequently used in the modelling of multiphase flow since (Bastian, 1999; Chen and Ewing, 1997b; Sin, 2015). Global pressure is defined as:
Instead of using one of the phase pressures as dependent variable as in the formulations presented before, this formulation solves the equation (2-19) using the global pressure as unknown. The system of equations is completed with the wetting phase equation (2-17) solved with for the wetting saturation written in terms of total velocity:
The global pressure and total mobility are always continuous, even in the absence of one phase, favoring the convergence of the models. Furthermore, the total velocity is smoother than the phase velocities. However, the global pressure variable is not an actual physical property and initial and boundary conditions do not have a direct physical meaning with this formulation.
Capillary pressure formulation
The formulation that uses capillary pressure and the non-wetting pressure as dependent variables is also relevant in the literature (O. Kolditz et al., 2012b; Neumann et al., 2013). The system of equations is formed by the wetting phase equation (2-17) and the non-wetting phase equation (2-18) solved with the non-wetting pressure ( ) and the capillary pressure ( ) as unknowns respectively:
Similar to the global pressure formulation, the use of the capillary pressure as dependent variable has the advantage that is continuous independently of the number of phases present in the cell. As a counterbalance, it is difficult to assign meaningful initial and boundary conditions with this formulation.
Numerical methods
The multiphase flow mathematical models are implemented in COMSOL Multiphysics (Comsol, 2016), a general-purpose finite element (FE) platform for the modelling and simulation of all kind of partial differential equations (PDE), ordinary differential equations (ODE) and algebraic differential equations (DAE) systems. It allows a flexible configuration of the shape functions, system properties, constitutive relations, solution schemes and solvers of the equations systems. COMSOL has a very intuitive and robust graphical user interface (GUI) with many options for geometry delineation, mesh generation and post-processing and visualization of the results. The multiphase flow PDE systems with their associated DAEs can be further combined with a set of physical interfaces for common physics applications areas available in COMSOL, accounting for coupled and multiphysics phenomena.
The main numerical characteristics (spatial discretization, solution scheme and linear equations system solvers) of the implemented multiphase flow equations are described in this section.
Spatial discretization
Many numerical methods have been developed for solving the convection-diffusion like equations through finite element methods (Chen et al., 2006; Helmig, 1997). A standard-Galerkin or Bubnov-Galerkin method based on linear continuous Lagrange shape or space functions is implemented in this thesis. The Galerkin method is based on the mean weighted residuals method where shape functions are equal to the weighting (also known as test) functions.
Alternative methods use two different finite element spaces, one for dependent scalar variables (e.g., pressure) and other for vector variables (e.g., fluid velocity). These methods are called mixed finite element (MFE) and should satisfy the Babuska–Brezzi (Brezzi et al., 1985) condition, which states that the shape functions for pressure must be of lower order than the shape functions for velocity. They were first introduced by (Nedelec, 1980; Raviart and Thomas, 1977) and many MFE were further developed in the literature. Two of the more common methods are the streamline upwind/Petrov-Galerkin (SUPG) method (Brooks and Hughes, 1982; Westerink and Shea, 1989; Zienkiewicz et al., 2005) and the Galerkin least-squares (GLS) method (Codina, 1998). The advantage of these methods is that they approximate both the scalar and vector variables simultaneously, giving a high order approximation. However, their implementation is more complex and rarely used (Hinkelmann, 2005). The discretization of PDEs with discontinuous shape functions across inter-element boundaries, i.e. discontinuous Galerkin (DG) method (Reed and Hill, 1973), are also an option for solving advective dominated problems from which continuous FEM lack robustness.
Sharp saturation fronts of the convection-diffusion equation can become unstable in Galerkin finite element methods (Carrayrou et al., 2003; De Windt et al., 2003; Helmig, 1997). The instabilities can be overcome with the help of stabilization methods. There are various stabilization techniques. Streamline diffusion and crosswind diffusion (based on the SUPG and GLS methods) are among the most common (Hauke, 2001; Hauke and Hughes, 1994; Hughes and Mallet, 1986). In the miscible model presented in this thesis, an isotropic diffusion stabilization method is used. The streamline diffusion and crosswind diffusion stabilizations are also implemented.
Solution scheme
The nonlinear equations system of the multiphase flow is solved through the Newton iterative method. Newton’s iterations in time can be solved with explicit and implicit methods.
Implicit schemes are unconditionally stable but computationally expensive as a nonlinear system of equations must be solved in each time step. Explicit schemes demand less computational resources; however, they come with stability time-stepping restrictions. In nonlinear systems as the multiphase flow equations, the time steps required may be too small to solve a problem explicitly within a reasonable computational time. The implicit schemes allow larger times but they can have severe damping effects (i.e., numerical dispersion) (Allen, 1984; Settari and Aziz, 1975).
An intermediate scheme, the implicit pressure-explicit saturation (IMPES), has been used in the literature to achieve better stability without having as much computational cost (Sheldon and Cardwell, 1959; Stone and Garder, 1961). They have been used widely for the solution of incompressible decoupled Pressure-Saturation formulations. However, this scheme is not efficient with systems with strong nonlinearities.
The solution of all the coupled nonlinear equations simultaneously and implicitly is known as the fully implicit method (FIM) or the simultaneous solution (Douglas et al., 1959). Variables are updated every Newton iteration, leading to a stable scheme that allows relatively large time steps. However, in systems with many chemical components, the Jacobian matrix may be large, requiring high computational resources.
Other relevant implicit methods are sequential (SM) (MacDonald and Coats, 1970), where the saturation dependent functions (e.g., relative permeability) have the values of the previous Newton iteration. These schemes are more stable than IMPES and more computationally efficient than FIM. When using these methods, it should be taken into account that the accuracy on the evaluation of the saturation dependent functions influences the accuracy and stability of the numerical methods (Aziz and Settari, 1979; Peaceman, 1977; Settari and Aziz, 1975).
Finally, an alternative solution scheme is the adaptive implicit method (AIM), which seeks an efficient middle point between the FIM and the IMPES scheme (Russell, 1989; Thomas and Thurnau, 1983). This time-stepping algorithm switches between the implicit and explicit depending on the problem and gradients of each grid block. Various comparative studies of these solution schemes are available in the literature (Bastian, 1999; Chen et al., 2006; Marcondes et al., 2009).
The solution scheme for the implemented multiphase flow equations is user-defined. The user can apply the solution schemes provided by COMSOL (all the previously explained except the IMPES). The models presented in this thesis are solved with the backward differentiation formulas (BDF) fully implicit scheme. The generalized-method has been also successfully applied in FIM and implicit sequential schemes in applications not shown in this thesis.
Table of contents :
1 Introduction
2 Numerical tools
2.1.1 Mathematical description
2.1.2 Numerical methods
2.1.3 Verification
2.2.1 Mathematical model
2.2.2 Numerical method
2.2.1 Verification:
3 Carbon capture and storage Interactions of CO2 gravity currents, capillarity, dissolution and convective mixing in a syncline anticline system.
3.2.1 Introduction
3.2.2 Model Description
3.2.3 Results and Discussion
3.2.4 Concluding remarks
3.3.1 Introduction
3.3.2 Model description
3.3.3 Results
3.3.4 Conclusions
4 Underground hydrogen storage Quantitative assessment of seasonal underground hydrogen storage from surplus energy in Castilla-León (north Spain)
4.2.1 Introduction
4.2.2 Problem definition
4.2.3 Mathematical Model
4.2.4 Results and discussion
4.2.5 Conclusions
5 General conclusions and perspectives
References
