Thermal simulation of flow in porous media is used in many applications where the coupling between temperature and the other variables is important, from heavy oil recovery to CO sequestration. To access huge reserves of heavy oil (Briggs et al., 1988) there has been a strong interest in thermal recovery methods such as steam injection, Steam Assisted Gravity Drainage (SAGD) and In-Situ Combustion (ISC) (Prats, 1982; Burger and Sahuquet, 1977; Lake, 1989). These methods provide heat to the reservoir in the form of hot fluids (steam) or generate heat in-situ. ISC processes inject air into a heated reservoir, reaching spontaneous ignition and using the exothermic oxidation reactions to lower the oil viscosity, allowing it to be displaced by the injected fluids. Although these recovery processes have been successfully applied at many fields, the complexity of the physical phenomena makes predictive numerical simulation extremely challenging. This work tackles one critical need in particular, the necessity of having a robust and scalable linear solver strategy.
The first thermal-compositional-reactive simulation models were developed in the late 1970s and early 1980s (Burger and others, 1976; Crookston et al., 1979; Coats and others, 1980; Youngren and others, 1980; Young et al., 1983; Rubin et al., 1985). Due to limitations on the available computing power at the time, those models typically used a small number of components and simplified physical models for phase behavior, as well as basic solvers. The development of fast and accurate algorithms for flash calculations (Michelsen, 1982a, b; Whitson and Michelsen, 1989)
paved the way for integrated simulators capable of simulating the tightly coupled system of non-linear, time dependent Partial Differential Equations (PDEs)(CMG, 2016; Schlumberger, 2015; Cao, 2002; Lapene et al., 2010). Discretizing these PDEs typically leads to a series of large, non-symmetric, ill-conditioned linear systems. Direct linear solvers (Li, 2005; Kourounis et al., 2018) can be used to obtain exact solutions, but their memory requirements severely limit the number of components and grid size that can be handled. The industry-standard approach for reservoir simulation is now to use an iterative Krylov subspace solver (Saad, 2003), such as Generalized Minimum RESidual (GMRES) (Saad and Schultz, 1986), paired with a suitable preconditioning scheme.
Recently, the reservoir simulation community has focused heavily on novel preconditioners for multiphysics problems, particularly coupled flow and geomechanics (White and Borja, 2011; Haga et al., 2012; Gries, 2015; White et al., 2016; Gaspar and Rodrigo, 2017; White et al., 2019). Although the physical processes are different, the idea of using approximate Schur-complements to extract easy-to-address subsystems is at the core of all multi-stage preconditioners. In the field of reservoir simulation, most isothermal models use the two-stage Constrained Pressure Residual (CPR) preconditioner introduced in Wallis and others (1983); Wallis et al. (1985). Initially designed for Black-Oil models, CPR decouples the pressure unknown, which exhibits essentially elliptic behavior, from saturation unknowns, which exhibit essentially hyperbolic behavior. The most common implementation, as described in Lacroix et al. (2003) or Cao et al. (2005), uses Algebraic Multigrid (AMG) (Ruge and Stüben, 1987; Stüben, 2001) for the first stage and an Incomplete LU (ILU) factorization for the second stage. AMG is well-suited to the elliptic pressure system, while ILU readily corrects local errors associated with the saturation field. CPR in isothermal compositional simulations is also known to perform well (Cusini et al., 2015; Moncorgé et al., 2018; Møyner et al., 2018).
For thermal compositional simulations, heat conduction through the rock can be a dominant mechanism, leading to CPR convergence degradation. Treating the temperature variable in the ILU second stage poorly approximates long-range temperature coupling. Li et al. (2015, 2017)
illustrated severe convergence issues using CPR on steam injection and SAGD cases. Their solution, called Enhanced CPR (ECPR), is to retain more variables in the first stage according to a given heuristic derived from the coupling strength between variables. While an improvement over traditional CPR, this approach has two drawbacks: (1) the expanded first stage subsystem loses any physical meaning, and (2) there is no guarantee the resulting matrix will be well suited for AMG. In another relevant work,Roy et al. (2019a, b) recently presented a detailed study of preconditioning strategies for two-phase dead-oil hot-water injection, including a new Schur-complement expression for the temperature subsystem arising from the specific physics they were studying.
In this paper we propose a family of multi-stage preconditioners suitable for general thermal reservoir simulation, showing reduced sensitivity to the thermal regime and good convergence behavior. We test these preconditioners on two dimensional, thermal-compositional cases with no reactions, as well as on reactive (combustion) cases. We compare these results to traditional CPR as a baseline. Another key contribution is a specific ordering strategy for the unknowns that guarantees well-posed primary and secondary systems of equations.
The paper is structured as follows. Section 2 presents the governing formulation. Section 3 presents an unknown ordering and static-condensation strategy. Section 4 provides a description of traditional CPR and proposes a new family of CPR-like preconditioners. Section 5 presents numerical results for non-reactive and reactive test cases across different thermal regimes for both homogeneous and heterogeneous media. We conclude in Section 6 with some avenues for future work.
2 Problem Statement
For the physical model considered here, the coupled mass and energy conservation equations can be cast as a coupled system of Advection-Diffusion-Reaction (ADR) equations,
with the conserved variable, the diffusivity, v the velocity field, and the source term representing both chemical reactions and well sources/sinks.
2.1 Mass Conservation Equations
Mass conservation for each component across phases , with components and phases, reads
Here, is the porosity; is the mole fraction of component in phase ; , , , and are the molar density, saturation, velocity, and volumetric flow rate of phase ; and is the source term from reactions. The equation is dimensional and has units of [mol/day]. We do not consider mass diffusion. In this work, we can also have solid, immobile components in the solid phase (with ).
2.2 Energy Conservation Equation
For thermal simulations, we also consider the energy conservation,
Here, is the temperature; is the rock volumetric internal energy; is the thermal conductivity; and are the internal energy and enthalpy of phase ; and is the source term from reactions. The equation is dimensional and has units of [J/day]. Thermal diffusion (also called conduction) usually cannot be neglected due to the large heat conductivity of the rock matrix, leading to one more term in the energy equation than in the mass equations.
2.3 Local Constraints
where and are phases indices, and is the fugacity of component in phase . Fugacity has units of pressure [bars]. These fugacity constraints give equations. In each phase, the sum of all molar fractions must be equal to 1, giving equations,
indicating the mole fraction of component in different fluid phases (oil, water, gas). We will also separate the solid components from the others and denote by the number of mobile components.
2.4 Treatment of Water
with , and the oil, vapor and water phase molar fractions and the overall mole fraction of component . These equations state that the hydrocarbon components cannot be dissolved in the water phase (Eq. 6) and that the water component cannot be dissolved in the oil phase (Eq. 7). This removes fugacity equality constraints by construction, as well as the water phase constraint equation . We remove the unknowns that are identically zero (all except for water, and ), as well as since it is identically one. The final size of the global system is .
We use the Automatic-Differentiation General Purpose Research Simulator (AD-GPRS (Voskov et al., 2012)) to solve the system of equations using a finite volume, fully implicit discretization and the natural variables formulation. More details about the reservoir simulation implementation can be found in Aziz and Settari (1979); Coats and others (1980); Cao (2002); Voskov et al. (2009). We omit details here as they are not the central focus of this work. At a given time step, we simply note that this discretization leads to a nonlinear set of residual equations,
for the algebraic vector of grid-cell unknowns. This system is solved using Newton’s method with an Appleyard chopping algorithm (Schlumberger, 2015)
to improve convergence robustness. We use limits on the relative variable changes, as well as a check to ensure all the bounded variables remain in the physical range. Given a solution estimate, an improved estimate is determined by
Here, is the Jacobian system evaluated at and the vector is the Newton update for timestep .
3 Static-Condensation Strategy
To lower the computational cost of solving the linear system, we first note that it can be partitioned into primary and secondary unknowns, as
For a suitable choice of partitioning (described below) the secondary unknowns are only coupled within their own grid cell, but are not coupled across cell boundaries. The block is then block-diagonal, and the secondary unknowns can be readily eliminated through static condensation. This leads to a reduced system
There are several subtleties that must be addressed to create a successful condensation strategy. In particular, phases and components can disappear or reappear in a grid cell as the simulation evolves. If a phase appears or disappears, we expand or shrink the set of variables according to the natural-variable formulation (see Cao (2002)). In all cases, the number of primary variables should be , corresponding to the number of globally-coupled equations. There is no unique guidance, however, on the specific primary/secondary partitioning adopted.
For a general three-phase flow simulation (ignoring the immobile solid phase) we can have seven different phase states in a cell. Due to the specific nature of the problem studied here, however, we only encounter three cases:
Gas (G) cell. The water has been vaporized and the oil was either vaporized or burnt. There are then unknowns: pressure, temperature, solid concentrations and vapor mole fractions.
Oil-Gas (OG) cell. The water component has been vaporized and is only present in the gas phase, but the oil phase is present. There are unknowns, adding the oil mole fractions for all non-water components and the gas saturation.
Oil-Water-Gas (OWG) cell. All mobile phases are present. There are unknowns, adding the oil saturation.
In all cases, we need to pick an appropriate subset of primary unknowns. A mandatory requirement is that the resulting block is invertible. This is achieved by aligning equations with properly chosen unknowns within each grid cell. We also note that it is convenient to avoid pivoting to the extent possible when applying Gaussian-Elimination to these small diagonal blocks. We therefore try to avoid zero diagonal entries appearing in the alignment strategy.
A few choices are immediately clear. The solid problem has equations and unknowns. The simplified model adopted here has only one solid concentration. We always align the solid conservation equation with the solid concentration, guaranteeing a non-zero diagonal value. Due to the different nature (and scaling) of the energy conservation equation, we always align it with temperature.
We now work on the flow unknowns to satisfy the invertibility requirement. When a component is not present, its conservation equation will only show non-zero values in the columns corresponding to that component’s mole fractions ( and potentially ) and those values are 1. The same is true for the fugacity constraints. These considerations lead us to align equations with their corresponding unknowns as much as possible. However, recall that we need to pick primary unknowns in the flow problem. We always want to retain pressure as a primary unknown. The saturation columns in the Jacobian only have non-zero values in the conservation equations, since fugacities and phase constraints do not depend on saturations. If we consider saturations as secondary variables, we can introduce a full row of zeros in , making it singular. Therefore, if saturations are part of the global set of unknowns (i.e. we have two or three phases present), they need to be selected in the primary unknowns. Conversely, when water is present as a phase, we have a fugacity constraint for water. Since we do not retain the water mole fractions in the free-water context, the water K-value is only a function of , (both primary variables) and . We need to have as a secondary variable to avoid another row of zeros in for the 3-phase cells.
With these general considerations in mind, we now describe a recommended partitioning strategy for each of the three possible phase states.
3.1 Gas Cell
The pure gas case happens upstream of the combustion front, where the water has been vaporized and moved downstream, and the oil has been burnt. In our implementation, we then only have one secondary equation, the vapor phase constraint . There are no saturations in the set of flow unknowns, but we retain pressure. We can pick any of the vapor mole fraction as the sole secondary unknown, since it will be a non-zero value in the phase constraint equation. We do, however, need to be careful that this component does not disappear, since it would lead to a zero value on the diagonal for the first row (given the lack of dependence on pressure). We align pressure with the nitrogen conservation equation since it is the least likely component to disappear in a pure gas cell. Figure 1 shows the corresponding Jacobian block sparsity pattern for a gas cell.
3.2 Oil-Gas Cell
Oil-Gas (OG) cells are located between the combustion front and the water front. The temperature is high enough that water can only be present in the vapor phase, or can be absent altogether. If is it absent, we need to align it with its own mole fraction , so we identify water as the third component (). As all of the oil has not yet burned or displaced, the flow unknowns also include the gas saturation and the oil mole fractions, . Note that we do not have since the water component cannot be present in the oil phase. We pick pressure, gas saturation, the vapor mole fraction of water and oil mole fractions as primary unknowns, and all other vapor mole fractions plus the remaining two oil mole fraction as secondary unknowns. We choose for the primary block, so that we can align them with their corresponding components. In an OG cell, the least likely component to disappear is the heaviest hydrocarbon component. We align it with the gas saturation to make sure that we do not encounter a zero-diagonal value. We align the secondary equations with their respective fugacity equations, and the two remaining mole fractions (one vapor and one oil) to their respective phase constraints. Figure 2 shows the resulting Jacobian entries.
3.3 Oil-Water-Gas Cell
Finally, Oil-Water-Gas cells are downstream of the water condensation front, in the cold zone. All three phases are present, and we now add the oil saturation as well as the fugacity equation for water since it can now be present in two phases. The ordering is very similar to the OG case; we keep the oil saturation as the third primary variable and move to the set of secondary variables. We previously mentioned that we need in the set of secondary equations, and we now can be assured that the water conservation equation will never be an issue since water is present by construction in OWG cells. We order the secondary unknowns similarly to the OG case, making sure the fugacity constraints are properly aligned to their respective components. Figure 3 shows the final ordering in the Jacobian.
In summary, we designed an ordering scheme for unknowns in the natural formulation to ensure that the Jacobian matrix is well-suited to the numerical methods we use to solve the linear system. Table 1 summarizes the set of unknowns in all three cases (G, OG, OWG).
|Case||Flow Primary Unknowns||Flow Secondary Unknowns|
We note that although this ordering scheme has been designed for combustion cases, it will perform well for thermal-compositional cases (no reactions) as well. All of the previously mentioned considerations still apply since both phases saturations and components mole fractions can be driven to zero by displacement or vaporization. With this particular partitioning, the reduced system may be readily assembled and solved with an appropriate preconditioning strategy.
Reservoir problems are typically driven by well boundary conditions, which have an important impact on the linear system and resulting solver strategy. In this work, we eliminate all of the well unknowns using an exact Schur complement decomposition as suggested in Zhou et al. (2013). All of our cases use two wells with single perforations, allowing us to exactly invert the 2x2 block while preserving the sparsity pattern of the reservoir unknowns. More sophisticated well treatments, while important, lie outside the scope of the current work.
4 Multi-Stage Preconditioning
We now explore several preconditioning schemes for the primary system . These schemes differ in the specifics, but all of them may be cast as a split preconditioning,
with left-preconditioning operator and a right-preconditioning operator . The system matrix has a lower-upper block factorization . As a guiding principle, all of the preconditioners below attempt to approximate and . The details of how these approximations are made, however, will differ from one scheme to the next.
For all schemes, the left “preconditioner” is just a scaling operation that may be explicitly applied before entering the Krylov solver. We denote the scaled matrix and right hand side as
This cheap scaling pushes the system matrix closer to upper-block-triangular form. After scaling, a multi-stage right-preconditioner is then applied within the Krylov iterations to solve the preconditioned system
Multi-stage preconditioners are frequently used to efficiently tackle coupling in multi-physics problems. The global form of a multi-stage preconditioner can be formally written as,
with the number of stages, the preconditioner for , and
the identity matrix. In this work, we consider both two- and three-stage variants.
The system matrix contains three important groups of variables. It is therefore convenient to partition into a block-system as
where the subscript denotes saturations and mole-fractions, denotes pressure, and denotes temperature. We then consider three preconditioner variants:
A two-stage Constrained Pressure Residual (CPR) method
A two-stage Constrained Pressure Temperature Residual (CPTR) method
A three-stage Constrained Pressure Temperature Residual (CPTR3) method
The three-stage preconditioner works directly on the partitioning above. In the two-stage variants, some variables are re-grouped to produce a partitioning instead.
To unify the presentation, we adopt the following notational conventions: indicates an individual block of the original system; identifies a block of a first-level Schur complement approximation; identifies a block of a second-level Schur complement approximation; is a block that may be well approximated with a zero block; and denotes a (block-)diagonal matrix.
4.1 Constrained Pressure Residual (CPR)
CPR leverages the different nature of the strongly elliptic (pressure) and strongly hyperbolic (saturations) variables in reservoir simulation. Initially introduced in Wallis and others (1983); Wallis et al. (1985) and extended in Lacroix et al. (2003) and Cao et al. (2005), it was designed for isothermal, Black-Oil cases. It is also used in compositional (Cusini et al., 2015; Moncorgé et al., 2018; Møyner et al., 2018) and even thermal simulations (Li et al., 2015, 2017), where it tends to show convergence issues. For thermal cases, it treats temperature as a strongly-hyperbolic variable.
To be precise, let the subscript to denote the union of saturation and temperature variables, with the re-partitioned system
Note that in the matrix saturation and temperature variables are actually stored in an interleaved order, with small, dense blocks appearing for each grid cell. The first step (left-scaling) used for CPR is:
Here, and are block-diagonal approximations of and , respectively. Also,
If the block diagonal approximations were replaced with their exact equivalents, the effect of left-scaling would be to form an upper-block-triangular factor—i.e. . The first stage of the right-preconditioner is then given by
with —i.e. the approximation of via an AMG preconditioner. To correct errors associated with the as-yet-untouched saturation variables, the second stage is an ILU(0) sweep based on the full matrix. It is given as,
In the second-stage, it is often convenient to permute the matrix to interleave the pressure and hyperbolic variables, so that all unknowns are ordered cell-wise. This interleaving is particularly useful if block-ILU (BILU) variants are chosen.
In this work, we extract the diagonal of block matrices, such as , to use as approximations in the Schur-complements. The objective of the diagonal approximation is to preserve the sparsity pattern of the initial matrix in Schur-complement operations.
In thermal simulations, the temperature variable can exhibit strongly elliptic behavior in the presence of high thermal diffusivity. In that case, treating it as a hyperbolic variable in the second stage of CPR is unlikely to be a good approximation. In the next two preconditioner variants, we explore alternative treatments of the energy equation and temperature unknowns.
4.2 Two-stage Constrained Pressure-Temperature Residual (CPTR)
We first consider a two-stage approach, where we now group temperature and pressure as elliptic variables, denoted with subscript . This leads to the partitioning,
Similar to the CPR case, we use a lower-block triangular scaling to perform an approximate Schur reduction,
To precondition this system, we need to approximate the inverse of the subsystem , which is a 2x2 block,
To tackle this block monolithically, we would need a well-suited algebraic preconditioner designed for systems of PDEs. Example options include System-AMG (SAMG) (Gries et al., 2014) or BoomerAMG (Baker et al., 2010) from the Hypre library (Falgout and Yang, 2002). The latter has shown promising results for coupled PDEs arising from multiphase flow (Bui et al., 2017, 2018). The first-stage preconditioner is then simply
The second-stage preconditioner is identical to CPR,
This two-stage variant will be denoted CPTR in the remainder of this work.
4.3 Three-stage Constrained Pressure-Temperature Residual (CPTR3)
Finally, we consider a three-stage variant. In the previous approach, the first left scaling step leads to
To get closer to an upper block-triangular system, we perform a second scaling:
The first stage preconditioner is then
with . The second stage is
with . The last stage is the local ILU(0) sweep of the scaled matrix,
This three-stage variant will be denoted CPTR3 in the remainder of this work. Compared to CPR, we have one more subsystem solve for the temperature matrix. It is important to note that this algorithm can be readily added on top of any existing CPR implementation. It only requires block matrices and scalar AMG preconditioners, both prerequisites for CPR.
4.4 Summary Algorithms
5 Numerical Results
5.1 Thermal-Compositional Homogeneous Cases
We start by studying thermal-compositional cases, with no reactions and homogeneous properties (porosity and permeability). For all cases, we run the simulator (AD-GPRS) on a laboratory scale, 2D case using parameters for an extra-heavy oil from Venezuela (Lapene, 2010). Our compositional description has seven hydrocarbon components: nitrogen, oxygen, carbon dioxide, water and one coke solid species. This gives and , for a total of 14 unknowns per cell. Table 2 shows the main parameters of our simulations. We inject hot air into a mixture of oil, water and nitrogen for 100 minutes. The water phase gets vaporized by the hot temperature front, and light hydrocarbon components will be stripped from the oil phase. Although no chemical reactions are taking place, this case shows a strong coupling between the mass and energy transport through viscosity, density and relative permeability calculations.
|Initial Oil Saturation||0.4791||–|
|Initial Water Saturation||0.2048||–|
As representative linear test problems, we have output the Jacobian matrix and residual vector from the simulator at multiple timesteps at the first Newton iteration. All of the convergence results are then obtained with a right-preconditioned GMRES (Saad and Schultz, 1986) algorithm applied to the left-scaled system given in Equation (13), using a relative tolerance of 1e and no restart. Using homogeneous properties allows us to conduct a mesh refinement study and quantify the impact of the grid size on the linear solver convergence.
The absence of an enthalpy source (from reactions) in the domain allows us to quantify the relative magnitude of thermal advection versus diffusion in a straightforward way. The effects can be compared through the dimensionless thermal Péclet number (Incropera et al., 2007),
with a characteristic length, a characteristic velocity, a characteristic flow rate, and the density and the specific heat capacity of the fluid carrier and the global thermal conductivity. With no conduction (), Pe is infinite, and with no convection () Pe is zero. The only heat source in the absence of reactions is the injection well, where we know the composition, rate, density and heat capacity of the fluid. The characteristic length is the length of the domain, and the thermal conductivity an input parameter depending on the rock properties. Table 3 summarizes the values we use to compute our Péclet numbers. Note that to vary Pe, we will use the thermal conductivity of the rock () while keeping all fluid properties constant.
Figure 4 plots the temperature profile for three different Péclet numbers along a cross-section through the domain. We can clearly observe different regimes, transitioning from advection dominated with sharp temperature fronts (blue curve) to a diffusion dominated case with a smooth profile and clear long-range interactions (purple curve). The orange curve is a unit Péclet number case, in the transition regime. The nature of the energy equation will change according to the Péclet number, from hyperbolic for pure advection cases to elliptic for pure diffusion cases, with intermediate cases being parabolic.
We start by studying the quality of the various Schur complement approximations themselves. To do so, we use a direct solver (rather than AMG) to apply the exact operators,
Using a direct solver also allows for a straightforward implementation of CPTR, since we do not have to worry about a multi-PDE AMG routine. Figure 5 summarizes results for the min timestep, and a complete performance profile is given in Table 4.
We compare the multi-stage preconditioners with a single-stage ILU(0) preconditioner. The latter approach does not take any physics into account nor does it scale well for this problem. Unless the grid size is very small, it will not converge in a reasonable number of iterations for any of our cases. We remark that the finer () grid size exhibits strong grid-orientation effects when the Péclet number is far from unity, and therefore linear solver performance trends are likely corrupted for these cases. The reader is referred to Kozdon et al. (2009) for more details about grid-orientation problems for transport in porous media. While important, these challenges are outside the scope of the current work.
|Péclet||Grid Size||Matrix Size||Non-zeros||ILU0||CPR||CPTR||CPTR3||CPR-AMG||CPTR3-AMG|
With respect to convergence behavior for increasing grid resolution, we observe that both CPTR and CPTR3 are virtually unaffected. As expected, CPR struggles to converge when the Péclet number is low, since the conduction part of the temperature dominates and the ILU(0) second stage cannot efficiently reduce the low frequency modes in the energy equation. Those convergence issues worsen with the grid size.
CPTR here solves the pressure/temperature subsystem (of size 2) exactly using a direct solver. One then observes excellent iteration performance—though not timing performance due to the cost of the direct solves—regardless of the Péclet number and the grid size. The ILU(0) second stage is able to efficiently reduce the remaining high frequency error modes in the saturation/mole fraction variables. For the grid, CPTR outperforms CPR in terms of iteration counts by 76%, 48% and 83% for high, unit and low Péclet numbers respectively. CPTR3 introduces one more level of approximation than CPTR, and as a result it shows a slightly higher number of iterations. However, it is still able to perform well and for an overall lower computational cost. For the grid, CPTR3 outperforms CPR in terms of iteration counts by 57%, 46% and 80% for high, unit and low Péclet numbers respectively. These results demonstrate that a specific treatment of the energy equation and temperature unknown coupling significantly increases the performance and robustness of the preconditioners. The performance of the proposed methods is always better than CPR across all grid sizes and thermal regimes.
We now consider a more scalable implementation of the preconditioners, using AMG for the key subsystems. We denote the AMG versions of the preconditioners by adding the suffix “-AMG”. Given the good results using CPTR3 in the previous section, and the fact it only requires a scalar AMG implementation, we focus on this approach in the remainder of the work. CPTR remains an appealing strategy—and could perhaps outperform CPTR3—but it requires a robust system-AMG implementation. For our comparison, we use a classic AMG method (Ruge and Stüben, 1987) as implemented in the HSL_MI20 package (HSL, 2002)
with default parameters: symmetric Gauss-Seidel smoother, single V-cycle, and a direct coarse solver. We set the maximum coarse problem size to 1000 degrees of freedom. Figure6 shows the results for min, and a complete performance summary is also given in Table 4.
CPTR3-AMG is able to keep the number of GMRES iterations below 15 in all cases up to the grid size. Both the pressure and temperature subsystems are well approximated and the scaling with respect to grid size is good. Using the grid, CPTR3-AMG outperforms CPR-AMG by 64%, 36% and 74% for high, unit and low Péclet numbers respectively. CPR-AMG shows an increased number of iterations compared to CPR for unit and low Péclet numbers of about 50%. However, for the high Péclet number case, we see a much larger increase, leading to 150% more iterations compared to the direct sub-solver version.
To investigate further, we compute the eigenvalues of the pressure matrix for the grid and the CPR preconditioner. Figure 7 (top) shows the spectra for the high (left), unit (middle) and low (right) Péclet numbers using the CPR preconditioner. For the high Péclet number we find both positive and negative eigenvalues. In the CPR case, the first pressure decoupling step includes the energy equation, since the temperature is part of the variables (Eq. 18). It seems doing so corrupts the ellipticity of the matrix. We observe better behavior with the matrix we get with the CPTR3 preconditioner, for which the pressure decoupling is done purely with mole conservation equations. That matrix is negative definite as desired, as shown in Figure 7 (bottom).
A very important feature for a preconditioner is its ability to deal with different physical regimes that may arise over the course of the simulation. In thermal reservoir simulations, different parts of a reservoir can show different property values and local Péclet numbers can greatly vary. If the reservoir is composed of different rock types, the thermal conductivity can also vary and impact the local Péclet number. Figure 8
illustrates that the sensitivity to thermal regime is greatly reduced for the proposed preconditioners, both using direct solvers (top), and using AMG preconditioners (bottom). To provide a quantitative metric, we compute the coefficient of variation (CV), defined as the standard deviation divided by the mean. For CPR-AMG, we get 41.6%, but that number reduces to 12.8% for CPTR3-AMG.
5.2 Thermal–Compositional Heterogeneous Cases
So far we have considered lab scale cases with homogeneous properties (leading to a fairly uniform velocity profile), but in a typical reservoir the fluid velocities can exhibit differences of several orders of magnitude.
We now test our method on heterogeneous cases to confirm the trends and performance we observed in the previous section. We generate two different problems based on the SPE10 Model 2 dataset (Christie et al., 2001). These cases are identical to the homogeneous cases for all parameters and properties except for the porosity and permeability fields. We crop the initial SPE10 layers ( cells) to 60 cells in the x-direction to create two square patterns. We retain a quarter-five spot pattern, injecting at the bottom left corner and producing at the top right corner. Both squares are 0.36 meters wide and show a similar length scale and pore volume compared the homogeneous cases. Note that we do not vary the grid refinement here, in order to avoid any issues with property downscaling for a consistent comparison.
We start with the top layer, modeling a prograding near-shore environment. The property variations are smooth but the heterogeneity is quite strong, with the permeability ranging across four orders of magnitude. The problem has 50,400 primary unknowns with 1.6 millions non-zeros in the system matrix. We conduct the same numerical study as for the homogeneous case. Figure 10 shows the performance of both new preconditioners using direct solvers for the subsystems (top) and using AMG preconditioners (bottom). These results confirm the trend we observed in the homogeneous case, and we actually see even more improvement. The addition of heterogeneity strengthens the coupling between transport and thermal effects, since now the viscosity reduction is spatially dependent. For all Péclet numbers, the number of iterations remains under 31 with an average of 24 using CPTR3-AMG, whereas CPR-AMG requires at least 40 iterations and an average of 57.
We also test the bottom layer, representing a fluvial, channelized environment. The property variations are sharp between channels and levees, and the permeability spans more than six orders of magnitude. Figure 11 shows the performance of both new preconditioners using direct solvers for the subsystems (top) and using AMG preconditioners (bottom). SPE10’s bottom layer is one of the more challenging benchmark cases available regarding permeability variations. The coupling is even stronger than for the top layer, so the direct solver cases perform admirably. The subsystems become significantly more challenging, however, for AMG. Nevertheless, for the high Péclet case, CPR-AMG does not converge in under 200 iterations (it takes 227 iterations) but CPTR3-AMG makes that case tractable again (62 iterations). The average number of iterations across Péclet numbers for CPR-AMG is 123, while CPTR3-AMG achieves an average of only 49 iterations with a maximum of 62.
5.3 In-Situ Combustion Homogeneous Case
One recovery method for heavy and extra-heavy oil is In-Situ Combustion (ISC), which involves chemical reactions (see Crookston et al. (1979) and Coats and others (1980) for details). Hot air injection is only efficient if one can leverage the chemical properties of the oil and oxidize part of it. Oxidation reactions are highly exothermic and allow the remainder of the oil to be displaced by lowering the viscosity. For reference, the oil we used throughout this study is 6.5API and the viscosity at reservoir conditions (48C) is close to 10,000 cP, rendering it virtually immobile for secondary processes like waterflooding. We consider the reaction scheme,
Although simple, this scheme allows us to study the ideal combustion case of laying out the solid fuel (Coke) from the heaviest oil fraction (Oil) and burning the fuel in the presence of oxygen. We do not crack the oil and generate lighter products here, but simply use the enthalpy of reaction to lower the viscosity. In this case, one can completely consume both the oxygen and the fuel at the combustion front; the heat generated by combustion is located at the oxygen concentration front. We use a standard Arrhenius (Arrhenius, 1889) model and oxygen partial pressure to compute the reaction rates (see Nissen et al. (2015) for more details).
There are three main differences between the previous cases and combustion. First, we generate heat inside the domain. The Péclet number is no longer sufficient to fully describe the thermal regime, but should be used in conjunction with the Damköhler numbers. These compare reaction and advection () and reaction and diffusion (). We will not perform a full dimensionless analysis in this case, but will present results for a typical ISC case using properties of unconsolidated sand for the thermal conductivity and parameters adapted from Dechelette et al. (2006) for the reactions. If we still use the flow rate to compute the Péclet number, its value is , corresponding to a slightly diffusion dominated case. The second difference is that the reaction parameters (pre-exponential factor and activation energy) tend to be grid-sensitive and should be different across grid sizes to give similar temperature plateaus and front speeds (Kovscek et al., 2013; Nissen et al., 2015). Lastly, the beginning of the simulation is the ignition step, which is an unusual regime that we will not study in detail. To ensure the comparability of results, we ensure all cases for all grid sizes fall in the developed combustion regime and that the temperature profiles are similar.
Figure 12 shows the results using direct solvers (top) and the AMG preconditioners (bottom). The CPR solver struggles, partly because of significant thermal diffusion. Another important effect is the much stronger coupling between the temperature and the saturations/mole fractions due to the reaction terms. The combustion front consumes oil and oxygen to release carbon dioxide, along with a large amount of energy. Failing to capture that coupling leads to significantly higher number of iterations for CPR. CPTR and CPTR3 outperform CPR by 81% and 74%, respectively. We observe a significant increase in the number of iterations for all grid sizes when using AMG preconditioners. Both the temperature, and especially the pressure subsystems, are much more challenging in the ISC case. Even with those more challenging subproblems, both CPTR-AMG and CPTR3-AMG still outperform CPR-AMG by around 55% for the and grid sizes. In the context of thermal-compositional-reactive simulations, the non-linear solver is likely to struggle with the very tight coupling and the multi-scale nature of the physical phenomena. With many time-step cuts and small time steps, the number of linear solves per non-linear iterations will grow significantly, making it all the more important and valuable to be able to rely on a robust and fast preconditioning technique. CPTR3-AMG achieves that goal on our ISC test case, keeping the number of iterations below 22.
6 Conclusion and Discussion
In this paper we have presented two new multi-stage preconditioners for thermal-compositional(-reactive) flow in porous media. Using a dedicated treatment of temperature, either in a pseudo-two-stage (CPTR) or three-stage (CPTR3) fashion, the convergence of GMRES was significantly improved for all of our tests cases. The reduction in the number of iterations is 40-85% compared to the industry standard CPR two-stage method, and at least one order of magnitude compared to a single-stage ILU(0) preconditioner on practical grid sizes for both homogeneous and heterogeneous cases.
The sensitivity of the proposed methods to the thermal regime, described by the thermal Péclet number in the absence of reactions, is greatly reduced and yields a much more robust preconditioner for a variety of reservoir simulation conditions. The use of AMG for the pressure and temperature subsystems and an ILU(0) final stage in CPTR3 shows good performance across all cases.
An interesting direction for future work is to apply a system AMG strategy for the elliptic subsystem of CPTR. Promising results have been obtained using SAMG for flow and transport problems (Gries et al., 2014), and BoomerAMG for linear elasticity problems Baker et al. (2010). This work has also tackled just one element of improving the computational speed and robustness of thermal-compositional-reactive simulations, which exhibit many well-known computational bottlenecks.
The authors want to thank Prof. Hamdi A. Tchelepi and Prof. Margot G. Gerritsen for interesting discussions and suggestions. Funding was provided by TOTAL S.A. through the FC-MAELSTROM project. Portions of this work were performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07-NA27344.
- Über die reaktionsgeschwindigkeit bei der inversion von rohrzucker durch säuren. Zeitschrift für physikalische Chemie 4 (1), pp. 226–248. Cited by: §5.3.
- Petroleum reservoir simulation. Chapman & Hall. Cited by: §2.5.
Improving algebraic multigrid interpolation operators for linear elasticity problems. Numerical Linear Algebra with Applications 17 (2-3), pp. 495–517. Cited by: §4.2, §6.
- Development of heavy-oil reservoirs. Journal of petroleum technology 40 (02), pp. 206–214. Cited by: §1.
- Algebraic multigrid preconditioners for multiphase flow in porous media. SIAM Journal on Scientific Computing 39 (5), pp. S662–S680. Cited by: §4.2.
Algebraic multigrid preconditioners for two-phase flow in porous media with phase transitions. Advances in water resources 114, pp. 19–28. Cited by: §4.2.
- Spontaneous ignition in oil reservoirs. Society of Petroleum Engineers Journal 16 (02), pp. 73–81. Cited by: §1.
- Les methodes thermiques de production des hydrocarbures-combustion in-situ, principes et etudes de laboratoire. Revue de L’Inst. Franrais du Petrole (2), pp. 141–88. Cited by: §1.
- Parallel scalable unstructured cpr-type linear solver for reservoir simulation. In SPE Annual Technical Conference and Exhibition, Cited by: §1, §4.1.
- Development of techniques for general purpose simulators. Ph.D. Thesis, Stanford University Stanford, CA. Cited by: §1, §2.5, §3.
- Tenth spe comparative solution project: a comparison of upscaling techniques. In SPE Reservoir Simulation Symposium, Cited by: §5.2.
- STARS. Technical description. Cited by: §1.
- In-situ combustion model. Society of Petroleum Engineers Journal 20 (06), pp. 533–554. Cited by: §1, §2.5, §5.3.
- A numerical simulation model for thermal recovery processes. Society of Petroleum Engineers Journal 19 (01), pp. 37–58. Cited by: §1, §5.3.
- Constrained pressure residual multiscale (cpr-ms) method for fully implicit simulation of multiphase flow in porous media. Journal of Computational Physics 299, pp. 472–486. Cited by: §1, §4.1.
- Air injection-improved determination of the reaction scheme with ramped temperature experiment and numerical simulation. Journal of Canadian Petroleum Technology 45 (01). Cited by: §5.3.
- Hypre: a library of high performance preconditioners. In International Conference on Computational Science, pp. 632–641. Cited by: §4.2.
- On the fixed-stress split scheme as smoother in multigrid methods for coupling flow and geomechanics. Computer Methods in Applied Mechanics and Engineering 326, pp. 526–540. Cited by: §1.
- Preconditioning for efficiently applying algebraic multigrid in fully implicit reservoir simulations. SPE Journal 19 (04), pp. 726–736. Cited by: §4.2, §6.
- System-amg approaches for industrial fully and adaptive implicit oil reservoir simulations. Ph.D. Thesis, Universität zu Köln. Cited by: §1.
- A parallel block preconditioner for large-scale poroelasticity with highly heterogeneous material parameters. Computational Geosciences 16 (3), pp. 723–734. Cited by: §1.
- A collection of fortran codes for large scale scientific computation, 2002. Cited by: §5.1.
- Fundamentals of heat and mass transfer. Wiley. Cited by: §5.1.
- Phase equilibrium computations are no longer the bottleneck in thermal compositional eos based simulation. In SPE Reservoir Simulation Symposium, Cited by: §2.4.
- Toward the next generation of multiperiod optimal power flow solvers. IEEE Transactions on Power Systems 33 (4), pp. 4005–4014. Cited by: §1.
- Improved predictability of in-situ-combustion enhanced oil recovery. SPE Reservoir Evaluation & Engineering 16 (02), pp. 172–182. Cited by: §5.3.
- Robust multi-d transport schemes with reduced grid orientation effects. Transport in porous media 78 (1), pp. 47–75. Cited by: §5.1.
- Iterative solution methods for modeling multiphase flow in porous media fully implicitly. SIAM Journal on Scientific Computing 25 (3), pp. 905–926. Cited by: §1, §4.1.
- Enhanced oil recovery. Cited by: §1.
- Three-phase free-water flash calculations using a new modified rachford–rice equation. Fluid Phase Equilibria 297 (1), pp. 121–128. Cited by: §1, §2.4.
- Etude expérimentale et numérique de la combustion in-situ d’huiles lourdes. Ph.D. Thesis. Cited by: §5.1.
- Enhanced constrained pressure residual ecpr preconditioning for solving difficult large scale thermal models. In SPE Reservoir Simulation Conference, Cited by: §1, §4.1.
- A parallel linear solver algorithm for solving difficult large scale thermal models. In SPE Reservoir Simulation Symposium: Society of Petroleum Engineers, doi, Vol. 10. Cited by: §1, §4.1.
- An overview of SuperLU: algorithms, implementation, and user interface. ACM Transactions on Mathematical Software 31 (3), pp. 302–325. Cited by: §1.
- The isothermal flash problem. part i. stability. Fluid phase equilibria 9 (1), pp. 1–19. Cited by: §1.
- The isothermal flash problem. part ii. phase-split calculation. Fluid Phase Equilibria 9 (1), pp. 21–40. Cited by: §1.
- Sequential fully implicit formulation for compositional simulation using natural variables. Journal of Computational Physics 371, pp. 690–711. Cited by: §1, §4.1.
- A mass-conservative sequential implicit multiscale method for isothermal equation-of-state compositional problems. SPE Journal. Cited by: §1, §4.1.
- Upscaling kinetics for field-scale in-situ-combustion simulation. SPE Reservoir Evaluation & Engineering 18 (02), pp. 158–170. Cited by: §5.3, §5.3.
- Thermal recovery. Cited by: §1.
- A block preconditioner for non-isothermal flow in porous media. Journal of Computational Physics. Cited by: §1.
- A constrained pressure-temperature residual (CPTR) method for non-isothermal multiphase flow in porous media. arXiv preprint arXiv:1907.04229. Cited by: §1.
- A general purpose thermal model. Society of Petroleum Engineers Journal 25 (02), pp. 202–214. Cited by: §1.
- Algebraic multigrid. In Multigrid methods, pp. 73–130. Cited by: §1, §5.1.
- GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical computing 7 (3), pp. 856–869. Cited by: §1, §5.1.
- Iterative methods for sparse linear systems. Vol. 82, siam. Cited by: §1.
- ECLIPSE reservoir simulator. Technical description. Cited by: §1, §2.5.
- A review of algebraic multigrid. In Numerical Analysis: Historical Developments in the 20th Century, pp. 331–359. Cited by: §1.
- Technical description of ad-gprs. Energy Resources Engineering, Stanford University. Cited by: §2.5.
- General nonlinear solution strategies for multiphase multicomponent eos based simulation. In SPE Reservoir Simulation Symposium, Cited by: §2.5.
- Constrained residual acceleration of conjugate residual methods. In SPE Reservoir Simulation Symposium, Cited by: §1, §4.1.
- Incomplete gaussian elimination as a preconditioning for generalized conjugate gradient acceleration. In SPE Reservoir Simulation Symposium, Cited by: §1, §4.1.
- Block-preconditioned newton–krylov solvers for fully coupled flow and geomechanics. Computational Geosciences 15 (4), pp. 647. Cited by: §1.
- A two-stage preconditioner for multiphase poromechanics in reservoir simulation. Computer Methods in Applied Mechanics and Engineering 357, pp. 112575. Cited by: §1.
- Block-partitioned solvers for coupled poromechanics: a unified framework. Computer Methods in Applied Mechanics and Engineering 303, pp. 55–74. Cited by: §1.
- The negative flash. Fluid Phase Equilibria 53, pp. 51–71. Cited by: §1.
- A generalized compositional approach for reservoir simulation. Society of Petroleum Engineers Journal 23 (05), pp. 727–742. Cited by: §1.
- Development and application of an in-situ combustion reservoir simulator. Society of Petroleum Engineers Journal 20 (01), pp. 39–51. Cited by: §1.
- A scalable multistage linear solver for reservoir models with multisegment wells. Computational Geosciences 17 (2), pp. 197–216. Cited by: Remark 1.