What is DARTS?
Delft: we belong to Civil Engineering and Geoscience (CEG) Department at Civil Engineering Faculty of TU Delft. The development team directly linked to the new GeoEnergy program which connects Geology, Geophysics and Petroleum Engineering sections of the department.
Advanced: the simulation framework is based on the recently proposed Operator-Based Linearization (OBL) approach, which helps to decouple the complex nonlinear physics and advanced unstructured discretization from the core simulation engine. The framework is targeting the solution of forward and inverse problems.
Research: the development team includes five PhD students from the CEG department and multiple MSc students working on their thesis project on DARTS platform. The simulation platform is developed within Delft Advanced Reservoir Simulation (DARSim) program and linked to multiple research in the area of reservoir simulation, inverse modeling and uncertainty quantification.
Terra: the developed framework is utilized for forward and inverse problems in petroleum engineering, low- and high-enthalpy geothermal applications, subsurface storage and subsurface integrity. The primary focus and developed capabilities are currently cover low-enthalpy geothermal operations which include multicomponent multiphase flow of mass and heat with complex chemical interactions. Another focus is thermal-compositional processes for Enhanced Oil Recovery.
Simulator: the main simulation kernel implemented on multi-core CPU and many-core GPU architectures. Advance multiscale nonlinear formulation improves the performance of forward and inverse models. Additional afford invested in representative proxy models for complex subsurface processes including multiphase multicomponent flow.
Introduction
DARTS is constructed within the Operator-based Linearization (OBL) framework for general-purpose reservoir simulations. The ‘super-engine’ in DARTS contains the governing equations characterizing the general thermal-compositional-reactive system should be expressed in form of operators, which are the functions of the thermodynamic state. These operators are then utilized to construct the residual and Jacobian for the solution of a nonlinear system using interpolation in thermodynamic parameter space. The values and derivatives of these operators at the supporting points will be evaluated adaptively during the simulation or pre-calculated in form of tables through the physics implemented in either C++ or Python.
The assembly of resulting Jacobian matrix in C++ is generalized and improved by the OBL approach, which greatly facilitates code development. The constructed linear system is passed to the dedicated linear solver for a solution. More details about the DARTS architecture are described in [1].
To enable the convenient usage of DARTS, the functionalities in C++ are exposed to users via a Python interface. At the same time, the Python interface provides the possibility to integrate different external modules (e.g., packages for physical property calculation) into DARTS.
Table of Contents
Conservation Equations
Mass and heat transfer involves a thermal multiphase flow system, which requires a set of equations to depict the flow dynamics. In this section, the governing equations and detailed spatial and temporal discretization and linearization procedures are introduced.
Governing Equations
For the investigated domain with volume \(\Omega\), bounded by surface \(\Gamma\), the mass and energy conservation can be expressed in a uniformly integral way, as
Here, \(M^c\) denotes the accumulation term for the \(c^{\mathrm{th}}\) component (\(c = 1, \ldots, n_c\), indexing for the mass components, [e.g., water, \(\mathrm{CO_2}\)] and \(c = n_c + 1\) for the energy quantity); \(\bf{F}_c\) refers to the flux term of the \(c^{\mathrm{th}}\) component; \({\bf{n}}\) refers to the unit normal pointing outward to the domain boundary; \(Q_c\) denotes the source/sink term of the \(c^{\mathrm{th}}\) component.
The mass accumulation term collects each component distribution over \(n_p\) fluid phases in a summation form,
where \(\phi\) is porosity, \(s_j\) is phase saturation, \(\rho_j\) is phase density \([\mathrm{kmol/m^3}]\) and \(x_{cj}\) is molar fraction of \(c\) component in \(j\) phase.
The energy accumulation term contains the internal energy of fluid and rock,
\begin{equation} \begin{aligned} M^{n_c+1} = \phi\sum\limits^{n_p}_{j=1}\rho_js_jU_j + (1 - \phi)U_r, \end{aligned} \end{equation}
where \(U_j\) is phase internal energy \([\mathrm{kJ}]\) and \(U_r\) is rock internal energy \([\mathrm{kJ}]\). The rock is assumed compressible and represented by the change of porosity through:
\begin{equation} \phi = \phi_0 \big(1 + c_r (p - p_\mathrm{ref}) \big), \end{equation}
where \(\phi_0\) is the initial porosity, \(c_r\) is the rock compressibility [1/bar] and \(p_\mathrm{ref}\) is the reference pressure [bars].
The mass flux of each component is represented by the summation over \(n_p\) fluid phases,
\begin{equation} \begin{aligned} \bf{F^c} = \sum\limits_{j=1}^{n_p}x_{cj}\rho_j \bf{u_j} + s_{j}\rho_{j} \textbf{J}_{cj}, \quad c = 1, \ldots, n_c. \end{aligned} \end{equation}
Here the velocity \(\bf{u_j}\) follows the extension of Darcy’s law to multiphase flow,
where \(\mathbf{K}\) is the permeability tensor \([\mathrm{mD}]\), \(k_{rj}\) is the relative permeability of phase \(j\), \(\mu_j\) is the viscosity of phase \(j\) \([\mathrm{mPa\cdot s}]\), \(p_j\) is the pressure of phase \(j\) [bars], \(\bf{\gamma_j}=\rho_j\bf{g}\) is the specific weight \([\mathrm{N/m^3}]\) and \(z\) is the depth vector [m].
The \(\textbf{J}_{cj}\) is the diffusion flux of component \(c\) in phase \(j\), which is described by Fick’s law as
where \(\textbf{D}_{cj}\) is the diffusion coefficient [m\(^2\)/day].
The energy flux includes the thermal convection and conduction terms,
\begin{equation} \begin{aligned} \bf{F}^{n_c+1} = \sum\limits^{n_p}_{j=1}h_j\rho_j \bf{u_j} + \kappa\nabla{T}, \end{aligned} \end{equation}
where \(h_j\) is phase enthalpy \([\mathrm{kJ/kg}]\) and \(\kappa\) is effective thermal conductivity \([\mathrm{kJ/m/day/K}]\).
Finally, the source term in mass conservation equations can be present in the following form
\begin{equation} \begin{aligned} {Q}^{c} = \sum\limits_{k=1}^{n_k}v_{ck}r_k, \quad c = 1, \ldots, n_c, \end{aligned} \end{equation}
where \(q_j\) is the phase source/sink term from the well, \(v_{ck}\) is the stoichiometric coefficient associated with chemical reaction \(k\) for the component \(c\) and \(r_{k}\) is the rate for the reaction. Similarly, the source term in the energy balance equation can be written as
\begin{equation} \begin{aligned} {Q}^{n_c+1} = \sum\limits_{k=1}^{n_k}v_{ek}r_{ek}. \end{aligned} \end{equation}
Here \(v_{ek}\) is the stoichiometric coefficient associated with kinetic reaction \(k\) for the energy and \(r_{ek}\) is the energy rate for kinetic reaction.
The nonlinear equations are discretized with the finite volume method using the multi-point flux approximation on general unstructured mesh in space and with the backward Euler approximation in time. For the \(i^{\mathrm{th}}\) reservoir block, the governing equation in discretized residual form reads:
Here \(V_i\) is the volume of the \(i^{th}\) grid block, \(\omega_{i}\) refers to state variables at the current time step, \(\omega^{n}_i\) refers to state variables at previous time step, \(A_l\) is the contact area between neighboring grids.
Operator Form of Governing Equations
DARTS includes capabilities for the solution of forward and inverse problems for subsurface fluid and heat transport. The OBL approach is employed in DARTS for the solution of highly nonlinear problems. It was proposed recently for generalized complex multi-phase flow and transport applications and aims to improve the simulation performance [2, 3]. For spatial discretization, a finite volume fully implicit method in combination with two-point flux approximation on unstructured grids is implemented in DARTS. Besides conventional discretization in temporal and spatial space, DARTS also utilizes discretization in physical space using the OBL approach.
With the OBL approach, the governing equations are written in form of state-dependent operators. The state-dependent operators can be parameterized with respect to nonlinear unknowns in multi-dimension tables under different resolutions. The values and derivatives of the operators in the parameter space can be interpolated and evaluated based on supporting points. For the adaptive parameterization technique [4], the supporting points are calculated `on the fly’ and stored for later re-usage, which can largely save time for parameterization in high-dimension parameter space (i.e. in multi-component compositional simulations). At the same time, the Jacobian assembly becomes flexible with the OBL, even for very complex physical problems.
Parameterization of the operators in 2D (pressure & enthalpy) space with a predefined OBL resolution [adapted from][5]. Here, the size of the quadrilateral represents the resolution for an operator interpolation.
Conservation of mass and energy
Pressure, temperature and overall composition are taken as the unified state variables in a given control volume in general-purpose thermal-compositional simulation. Upstream weighting of the physical state is used to determine the flux-related fluid properties determined at the interface \(l\). The discretized mass conservation equation in operator form for girdblock (here we omit \(i\)) reads:
\begin{equation} V\phi_0[ \alpha_c (\omega) -\alpha_c( \omega_n)]-\Delta t\sum_{l\in L(i)}\sum_{j=1}^{n_p}[\Gamma^l\beta_{cj}^l(\omega^u)\Delta\psi_j^l + \Gamma_d^l\gamma_{j}^l(\omega)\Delta \chi_{cj}]+\Delta t V \delta_c(\omega)=0 . \end{equation}
where \(V\) is the control volume, \(\omega_n\) is the physical state of block \(i\) at the previous timestep, \(\omega\) is the physical state of block \(i\) at the new timestep, \(\omega^{u}\) is the physical state of upstream block, \(\Gamma^l\) and \(\Gamma_d^l\) are the fluid and diffusive transmissibilities respectively and \(L(i)\) is a set of interfaces for gridblock \(i\).
Here we defined the following state-dependent operators,
The phase-potential-upwinding (PPU) strategy for OBL parametrization is applied in DARTS to model the gravity and capillary effect [4, 6]. The potential difference of phase \(j\) on the interface \(l\) between block 1 and 2 can be written as:
\begin{equation} \Delta \psi^l_{j} = p_1-p^c_{j}(\omega_1) - (p_2-p^c_{j}(\omega_2))-\frac{\rho_j(\omega_1)+\rho_j(\omega_2)}{2}g(z_2-z_1), \end{equation}
where \(p^c_{j}\) is the capillary pressure.
The discretized energy conservation equation in operator form can be written as:
where:
In addition, for accounting the energy of rock, three additional operators should be defined:
\(\alpha_{eri}\) and \(\alpha_{erc}\) represent the rock internal energy and rock conduction, respectively. \(U_r\) is a state-dependent parameter, thus these two rock energy terms are treated separately.
This agglomeration of different physical terms into a single nonlinear operator simplifies the implementation of nonlinear formulations. Instead of performing complex evaluations of each property and its derivatives with respect to nonlinear unknowns, operators can be parameterized in physical space either at the pre-processing stage or adaptively with a limited number of supporting points. The evaluation of operators during the simulation is based on multi-linear interpolation, which improves the performance of the linearization stage. Besides, due to the piece-wise representation of operators, the nonlinearity of the system is reduced, which improves the nonlinear behavior [3, 4]. However, to delineate the nonlinear behavior in the system, especially strong nonlinearity (e.g., at high-enthalpy conditions), it is necessary to select a reasonable OBL resolution to characterize the physical space. Too coarse OBL resolution may lead to a large error in the solutions [2].
Well treatment
A connection-based multi-segment well is used to simulate the flow in the wellbore [7]. The communication between well blocks and reservoir blocks is treated in the same way as between reservoir blocks. Besides, the top well block is connected with a ghost control volume, which is selected as a placeholder for the well control equations. The bottom hole pressure (BHP), volumetric and mass rate controls are available in DARTS to model various well conditions.
As for the BHP well control, the injector and/or producer will operate under fixed bottom-hole pressure. A pressure constraint is defined at the ghost well block:
\begin{equation} p-p^{target} = 0. \end{equation}
The volumetric rate control in DARTS is implemented through the volumetric rate operator \(\zeta^{vol}_p(\omega)\):
\begin{equation} \Gamma^l\zeta_j^{vol}(\omega)\Delta{p} - Q^{target} = 0, \end{equation}
where
where \(Q^{target}\) is the target volumetric flow rate at separator conditions \([\mathrm{m^3/day}]\), \(\beta_{cj}(\omega)\) is the mass flux operator as shown in (1), \(\hat{s}_j\) and \(\hat{\rho}_t(\omega)\) are the saturation and total fluid density respectively at separator conditions.
Any of the described well controls can be coupled with energy boundary conditions, defined by the temperature or enthalpy of the injected fluid at the injection well. Since temperature is the function of the thermodynamic state, it is expressed in operator form and the temperature well control reads:
\begin{equation} \chi(\omega) - T^{target} = 0, \end{equation}
where \(\chi(\omega)\) is defined by (2) and \(T^{target}\) is the target temperature of injected fluid. Alternatively, the enthalpy of the injected fluid can be defined:
\begin{equation} h(\omega) - h^{target} = 0, \end{equation}
where \(h\) is the enthalpy of the well control block, \(h^{target}\) is the target enthalpy of injected fluid. For the production well control, enthalpy is taken equal to that of the upstream well block.
Treatment of porosity
The porosity \(\phi\) depends on the concentrations of the minerals according to the relationship:
where \(M\) is the number of reactive minerals, \({\cal M}_m\) is the molar mass of mineral \(m\), \(\rho_m\) is the mass density of mineral \(m\) and \(c_{ms}\) represents and the molar concentration of mineral \(m\).
Let’s represent the total bulk volume of control element in simulation
\begin{equation} V = V_f + V_r + V_{nr}, \end{equation}
and \(V_r\) denotes reactive volume and \(V_{nr}\) represents the non-reactive volume (not altered by any chemical reaction). Dividing this by the total (bulk) volume gives
\begin{equation} 1 = \phi + \phi_r + \phi_{nr} = \phi^T + \phi_{nr}, \end{equation}
where \(\phi_r\) represents the reactive volume fraction, \(\phi_{nr}\) is the non-reactive volume fraction, and \(\phi^T\) is the total porosity defined as the sum of the fluid porosity and reactive volume fraction. Since only the reactive volume and fluid porosity can change due to chemical reactions, it follows directly that the total porosity remains constant throughout simulation (when neglecting compressibility). This and the changes in volume fractions due to precipitation and dissolution is illustrated in the following figure.
Schematic of the different volumes in the domain. The domain consists of three distinct regions, particularly the fluid volume which is occupies by all the mobile phases (liquid and gaseous in the case of two phase flow), the reactive volume which consist of solid phases that can react or precipitate, and finally the nonreactive volume (the part of the control volume which doesn’t participate in any chemical reaction)
Note that the fluid porosity can always be obtained with the following constitutive equation
where \(M\) is the number of solid phases (occupying the reactive volume fraction) and \(\hat{s}_m\) is the saturation of solid phase. Please note that the \(s_{\alpha}\) is the fluid saturation (defined over the pore volume) while \(\hat{s}_{m}\) is the {solid saturation of mineral phase \(m\)} (defined over the pore and reactive rock volume).
References
M. Khait. Delft Advanced Research Terra Simulator: General Purpose Reservoir Simulator with Operator-Based Linearization. PhD thesis, TU Delft, 2019. doi:10.4233/uuid:5f0f9b80-a7d6-488d-9bd2-d68b9d7b4b87.
D. V. Voskov. Operator-based linearization approach for modeling of multiphase multi-component flow in porous media. Journal of Computational Physics, 337:275–288, 2017. doi:10.1016/j.jcp.2017.02.041.
M. Khait and D. V. Voskov. Operator-based linearization for efficient modeling of geothermal processes. Geothermics, 74:7–18, 2018. doi:10.1016/j.geothermics.2018.01.012.
M. Khait and D. V. Voskov. Adaptive parameterization for solving of thermal/compositional nonlinear flow and transport with buoyancy. SPE Journal, 23:522–534, 2018. doi:10.2118/182685-PA.
M. Khait and D. V. Voskov. Operator-based linearization for general purpose reservoir simulation. Journal of Petroleum Science and Engineering, 157:990–998, 2017. doi:10.1016/j.petrol.2017.08.009.
X. Lyu, M. Khait, and D. Voskov. Operator-based linearization approach for modelling of multiphase flow with buoyancy and capillarity. SPE Journal, pages 1–18, 2021. doi:https://doi.org/10.2118/205378-PA.
M. Khait and D. V. Voskov. Integrated framework for modelling of thermal-compositional multiphase flow in porous media. In SPE Reservoir Simulation Conference. 2019.