Modeling and simulation of multiphase flow in porous media have been a major effort in reservoir engineering and in environmental study. Petroleum engineers use reservoir simulation models to manage existing petroleum fields and to develop new oil and gas reservoirs, while environmental scientists use subsurface flow and transport models to investigate and compare for example various schemes to inject and store CO2 in subsurface geological formations, such as depleted reservoirs and deep saline aquifers. One basic requirement for accurate modeling and simulation of multiphase flow is to have the predicted physical quantities sit within a physically meaningful range. For example, the predicated saturation should sit between 0 and 1 while the predicated molar concentration should sit between 0 and the maximum value allowed by the equation of state. Unfortunately, popular simulation methods used in petroleum industries do not preserve physical bounds. A commonly used fix to this problem is to simply apply a cut-off operator (say, to the computed saturation) at each time step, i.e., to set the saturation to be zero whenever it becomes negative, and to set it to one whenever it becomes larger than one. However, this cut-off practice does not only destroy the local mass conservation but it also damages the global mass conservation, which seriously ruins the numerical accuracy and physical interpretability of the simulation results. In the talk, we will present our recent work on bound-preserving discretization and solvers for subsurface flow models based on a fully implicit framework. We reformulated a few subsurface flow models using variational inequalities that naturally ensure the physical feasibility of the physical quantities including saturations and concentrations. We applied a mixed finite element method to discretize the model equations for the spatial terms, and the implicit backward Euler scheme with adaptive time stepping for the temporal integration. The resultant nonlinear system arising at each time step was then solved in a monolithic way by using a Newton–Krylov type method, where the resultant nonlinear system was solved by a generalized Newton method, i.e., active-set reduced-space method, and then the ill-conditioned linear Jacobian systems were solved with an effective preconditioned Krylov subspace method. The used nonlinear preconditioner was built by applying overlapping additive Schwarz type domain decomposition and nonlinear elimination. Numerical results will be presented to examine the performance of the newly developed algorithm on parallel computers. It was observed from numerical tests that our nonlinear solver overcomes the severe limits on the time step associated with conventional methods, and it results in superior convergence performance, often reducing the total computing time by more than one order of magnitude.
This presentation is based on the joint work with Haijian Yang (Hunan University), Chao Yang (Beijing University), and Yiteng Li (KAUST).