### Recent Submissions

• #### Base level changes based on Basin Filling Modelling: a case study from the Paleocene Lishui Sag, East China Sea Basin

(Petroleum Science, Springer Science and Business Media LLC, 2020-07-17) [Article]
Estimation of base level changes in geological records is an important topic for petroleum geologists. Taking the Paleocene Upper Lingfeng Member of Lishui Sag as an example, this paper conducted a base level reconstruction based on Basin Filling Modelling (BFM). The reconstruction was processed on the ground of a previously interpreted seismic stratigraphic framework with several assumptions and simplification. The BFM is implemented with a nonlinear diffusion equation solver written in R coding that excels in shallow marine stratigraphic simulation. The modeled results fit the original stratigraphy very well. The BFM is a powerful tool for reconstructing the base level, and is an effective way to check the reasonableness of previous interpretations. Although simulation solutions may not be unique, the BFM still provides us a chance to gain some insights into the mechanism and dynamic details of the stratigraphy of sedimentary basins.
• #### Fully discrete energy stable scheme for a phase-field moving contact line model with variable densities and viscosities

(Applied Mathematical Modelling, Elsevier BV, 2020-03-04) [Article]
In this study, we propose a fully discrete energy stable scheme for the phase-field moving contact line model with variable densities and viscosities. The mathematical model comprises a Cahn–Hilliard equation, Navier–Stokes equation, and the generalized Navier boundary condition for the moving contact line. A scalar auxiliary variable is employed to transform the governing system into an equivalent form, thereby allowing the double well potential to be treated semi-explicitly. A stabilization term is added to balance the explicit nonlinear term originating from the surface energy at the fluid–solid interface. A pressure stabilization method is used to decouple the velocity and pressure computations. Some subtle implicit–explicit treatments are employed to deal with convention and stress terms. We establish a rigorous proof of the energy stability for the proposed time-marching scheme. A finite difference method based on staggered grids is then used to spatially discretize the constructed time-marching scheme. We also prove that the fully discrete scheme satisfies the discrete energy dissipation law. Our numerical results demonstrate the accuracy and energy stability of the proposed scheme. Using our numerical scheme, we analyze the contact line dynamics based on a shear flow-driven droplet sliding case. Three-dimensional droplet spreading is also investigated based on a chemically patterned surface. Our numerical simulation accurately predicts the expected energy evolution and it successfully reproduces the expected phenomena where an oil droplet contracts inward on a hydrophobic zone and then spreads outward rapidly on a hydrophilic zone.
• #### Sorption and Diffusion of Methane and Carbon Dioxide in Amorphous Poly(alkyl acrylates): A Molecular Simulation Study

(The Journal of Physical Chemistry B, American Chemical Society (ACS), 2020-01-29) [Article]
Molecular simulations were carried out to understand the structural features and the sorption and diffusion behavior of methane and carbon dioxide in amorphous poly(alkyl acrylates) in the temperature range of 300-600 K. The hybrid Monte Carlo/molecular dynamics approach was employed to address the effects of polymer swelling and framework flexibility on the gas sorption. Simulations show that the glass-transition temperature decreases with the side-chain length of poly(alkyl acrylate), consistent with experiments. This is due to the fact that the shielding of the polar ester groups increases with the side-chain length. The simulated sorption isotherms for methane and carbon dioxide were in agreement with the experimental data. The polymer swelling becomes more pronounced, especially in the case of sorption of carbon dioxide. A significant swelling occurs, possibly because of the greater interaction between carbon dioxide and the polar ester groups in the polymers. The uptake of methane and carbon dioxide by poly(alkyl acrylates) generally increases with the side-chain length. Our simulations confirm the experimental findings that the diffusion coefficients of methane and carbon dioxide in poly(alkyl acrylates) increase with the side-chain length. Interestingly, the activation energies of gas diffusion decrease with the side-chain length. The diffusion coefficients of the penetrants have an exponential relationship with the void fraction, which is in agreement with the free volume theory.
• #### A Novel Energy Factorization Approach for the Diffuse-Interface Model with Peng--Robinson Equation of State

(SIAM Journal on Scientific Computing, Society for Industrial & Applied Mathematics (SIAM), 2020-01-07) [Article]
The Peng--Robinson equation of state (PR-EoS) has become one of the most extensively applied equations of state in chemical engineering and the petroleum industry due to its excellent accuracy in predicting the thermodynamic properties of a wide variety of materials, especially hydrocarbons. Although great effort has been made to construct efficient numerical methods for the diffuse interface models with PR-EoS, there is still not a linear numerical scheme that can be proved to preserve the original energy dissipation law. In order to pursue such a numerical scheme, we propose a novel energy factorization (EF) approach, which first factorizes an energy function into a product of several factors and then treats the factors using their properties to obtain the semi-implicit linear schemes. We apply the EF approach to deal with the Helmholtz free energy density determined by PR-EoS, and then propose a linear semi-implicit numerical scheme that inherits the original energy dissipation law. Moreover, the proposed scheme is proved to satisfy the maximum principle in both time semidiscrete form and cell-centered finite difference fully discrete form under certain conditions. Numerical results are presented to demonstrate the stability and efficiency of the proposed scheme.
• #### Special Issue: Advanced numerical modeling and algorithms for multiphase flow and transport

(Journal of Computational Physics, Elsevier BV, 2020-01-06) [Article]
• #### Bulk and interfacial properties of decane in the presence of carbon dioxide, methane, and their mixture.

(Scientific reports, Springer Science and Business Media LLC, 2019-12-24) [Article]
Molecular dynamics simulations were performed to study the bulk and interfacial properties of methane + n-decane, carbon dioxide + n-decane, and methane + carbon dioxide + n-decane systems under geological conditions. In addition, theoretical calculations using the predictive Peng-Robinson equation of state and density gradient theory are carried out to compare with the simulation data. A key finding is the preferential dissolution in the decane-rich phase and adsorption at the interface for carbon dioxide from the methane/carbon dioxide mixture. In general, both the gas solubility and the swelling factor increase with increasing pressure and decreasing temperature. Interestingly, the methane solubility and the swelling of the methane + n-decane system are not strongly influenced by temperature. Our results also show that the presence of methane increases the interfacial tension (IFT) of the carbon dioxide + n-decane system. Typically, the IFT of the studied systems decreases with increasing pressure and temperature. The relatively higher surface excess of the carbon dioxide + n-decane system results in a steeper decrease in its IFT as a function of pressure. Such systematic investigations may help to understand the behavior of the carbon dioxide-oil system in the presence of impurities such as methane for the design and operation of carbon capture and storage and enhanced oil recovery processes.
• #### A locally and globally phase-wise mass conservative numerical algorithm for the two-phase immiscible flow problems in porous media

(Computers and Geotechnics, Elsevier BV, 2019-12-14) [Article]
In this work, we introduce a novel numerical method to solve the problem of two-phase immiscible flow in porous media that is conservative to both phases. In the widely used implicit pressure, explicit saturation (IMPES) scheme, the conservation of mass of both the two phases are summed to form an equation involving the total Darcy's velocity. In the discretization of such an equation it becomes difficult to enforce the conservation of mass of each phase. To guarantee the conservation of mass of both phases locally and hence globally, we introduce a scheme in which the time discretization of the mass conservation equations is considered separately. Cell-centered finite difference (CCFD) methods are adopted for spatial discretization, where the variables of fluid properties (i.e. relative permeability and mobility) are upwinded separately according to the velocity of each phase and not according to the total velocity. Furthermore, this new scheme updates all phase velocities and uses them to update the corresponding phase saturation. In addition, a two-scale of time-splitting methods are adopted for pressure equation and saturation equations to improve the computational efficiency. For the sake of simplicity, we show a number of examples of two-phase system in two-dimensional geometry solved using the new scheme. It is shown that the new scheme is more embracing the physics and it can be more accurate than traditional IMPES scheme, particularly for the cases in which the phase velocities are in opposite direction, and conventional IMPES schemes fail.
• #### A phase-field moving contact line model with soluble surfactants

(Journal of Computational Physics, Elsevier BV, 2019-12-04) [Article]
A phase-field moving contact line model is presented for a two-phase system with soluble surfactants. With the introduction of some scalar auxiliary variables, the original free energy functional is transformed into an equivalent form, and then a new governing system is obtained. The resulting model consists of two Cahn-Hilliard-type equations and incompressible Navier-Stokes equation with variable densities, together with the generalized Navier boundary condition for the moving contact line. We prove that the proposed model satisfies the total energy dissipation with time. To numerically solve such a complex system, we develop a nonlinearly coupled scheme with unconditional energy stability. A splitting method based on pressure stabilization is used to solve the Navier-Stokes equation. Some subtle implicit-explicit treatments are adopted to discretize convection and stress terms. A stabilization term is artificially added to balance the explicit nonlinear term associated with the surface energy at the fluid-solid interface. We rigorously prove that the proposed scheme can preserve the discrete energy dissipation. An efficient finite difference method on staggered grids is used for the spatial discretization. Numerical results in both two and three dimensions demonstrate the accuracy and energy stability of the proposed scheme. Using our model and numerical scheme, we investigate the wetting behavior of droplets on a solid wall. Numerical results indicate that surfactants can affect the wetting properties of droplet by altering the value of contact angles.
• #### Generalized multiscale approximation of mixed finite elements with velocity elimination for subsurface flow

(Journal of Computational Physics, Elsevier BV, 2019-11-25) [Article]
A frame work of the mixed generalized multiscale finite element method (GMsFEM) for solving Darcy's law in heterogeneous media is studied in this paper. Our approach approximates pressure in multiscale function space that is between fine-grid space and coarse-grid space and solves velocity directly in the fine-grid space. To construct multiscale basis functions for each coarse-grid element, three types of snapshot space are raised. The first one is taken as the fine-grid space for pressure and the other two cases need to solve a local problem on each coarse-grid element. We describe a spectral decomposition in the snapshot space motivated by the analysis to further reduce the dimension of the space that is used to approximate the pressure. Since the velocity is directly solved in the fine-grid space, in the linear system for the mixed finite elements, the velocity matrix can be approximated by a diagonal matrix without losing any accuracy. Thus it can be inverted easily. This reduces computational cost greatly and makes our scheme simple and easy for application. Comparing to our previous work of mixed generalized multiscale finite element method [E. T. Chung, Y. Efendiev, and C. S. Lee. Mixed generalized multiscale finite element methods and applications. Multiscale Modeling and Simulation, 13(1):338-366, 2015.], both the pressure and velocity space in this approach are bigger. As a consequence, this method enjoys better accuracy. While the computational cost does not increase because of the good property of velocity matrix. Moreover, the proposed method preserves the local mass conservation property that is important for subsurface problems. Numerical examples are presented to illustrate the good properties of the proposed approach. If offline spaces are appropriately selected, one can achieve good accuracy with only a few basis functions per coarse element according to the numerical results.
• #### Recent Progress on Phase Equilibrium Calculation in Subsurface Reservoirs Using Diffuse Interface Models

(Springer International Publishing, 2019-11-17) [Book Chapter]
Compositional multiphase flow in subsurface porous media is becoming increasingly attractive due to issues related with enhanced oil recovery, greenhouse effect and global warming, and the urgent need for development in unconventional oil/gas reservoirs. One key effort prior to construct the mathematical model governing the compositional multiphase flow is to determine the phase compositions of the fluid mixture, and then calculate other related physical properties. In this paper, recent progress on phase equilibrium calculations in subsurface reservoirs have been reviewed and concluded with authors’ own analysis. Phase equilibrium calculation is the main approach to perform such calculation, which could be conducted using two different types of flash calculation algorithms: The NPT flash and NVT flash. NPT flash calculations are proposed early, well developed within the last few decades and now become the most commonly used method. However, it fails to remain the physical meanings in the solution as a cubic equation, derived from equation of state, is often needed to solve. Alternatively, NVT flash can handle the phase equilibrium calculations as well, without the pressure known a priori. Recently, Diffuse Interface Models, which were proved to keep a high consistency with thermodynamic laws, have been introduced in the phase calculation, incorporating the realistic equation of state (EOS), e.g. Peng-Robinson EOS. In NVT flash, Helmholtz free energy is minimized instead of Gibbs free energy used in NPT flash, and this energy density is treated with convex-concave splitting technique. A semi-implicit numerical scheme is designed to process the dynamic model, which ensures the thermodynamic stability and then preserve the fast convergence property. A positive definite coefficient matrix is designed to meet the Onsager Reciprocal Principle so as to keep the entropy increasing property in the presence of capillary pressure, which is required by the thermodynamic laws. The robustness of the proposed algorithm is verified via two numerical examples, one of which has up to seven components. In the complex fluid mixture, special phenomena could be capture from the global minimum of TPD functions as well as the phase envelope resulted from the phase equilibrium calculations. It can be found that the boundary between the single-phase and vapor–liquid phase regions will move in the presence of capillary pressure, and then the area of each region will change accordingly. Some remarks have been concluded at the end, as well as suggestions on potential topics for future studies.
• #### Layer Charge Effects on Adsorption and Diffusion of Water and Ions in Interlayers and on External Surfaces of Montmorillonite

(ACS Earth and Space Chemistry, American Chemical Society (ACS), 2019-10-28) [Article]
Molecular simulations using classical force fields were performed to study the adsorption and diffusion properties of water and ions in the high-charge (Arizona-type) montmorillonite clays with varying relative humidity (RH) at 298.15 K. The simulation results of basal distances derived from swelling free energy curves and of water uptake are in good agreement with experiments. Overall, the simulated self-diffusion coefficients of the interlayer species are in reasonable agreement with experiments and lower than those estimated for the external surfaces. Influence of the magnitude of the layer charge was studied by comparing these simulation results with those obtained for the low-charge (Wyoming-type) montmorillonite clays. Most importantly, these comparisons confirm the experimental finding that the high-charge clay generally shifts swelling transitions toward lower RH values. Therefore, the adsorption and dynamics of water and ions are significantly different in the low- A nd high-charge clays near the transition RH values. We find that the amount of water in an interlayer hydration state is mostly independent of the layer charge, probably due to steric reasons. In contrast, the externally adsorbed water content increases with increasing layer charge. Furthermore, the mobility of water and ions is generally lower in the high-charge montmorillonite than in the corresponding low-charge system. However, mobility of cations in the mesopores of high-charge montmorillonite is typically higher than that of the corresponding low-charge system which might be attributed to the presence of the tetrahedral substitutions in the latter case.
• #### Low salinity effect on the recovery of oil trapped by nanopores: A molecular dynamics study

(Fuel, Elsevier BV, 2019-10-25) [Article]
Low salinity waterflooding (LSW) is an effective method for enhancing the oil recovery from many reservoirs, and its success has been traced to a host of low salinity effects. In this work, we perform molecular dynamics simulations to study the feasibility of recovering oil trapped by nanopores by lowering the reservoir salinity. The oil is initially trapped by a slit nanopore, with a portion of the oil protruding from the pore entrance. After the reservoir salinity is lowered, the thin brine films that separate the oil and pore walls become thicker to drive some of the trapped oil out of the pore. We quantify the free energy profile of this process and clarify the underlying molecular mechanisms. Interestingly, the brine film growth is dominated by the water transport from the brine reservoir into the pore rather than by the depletion of ions from the brine film. These results provide molecular evidence that low salinity brines benefit the recovery of the oil trapped by nanopores. They highlight that when ion depletion from thin brine films is suppressed, the osmosis of water can play a fundamental role in the expansion of the brine films; thus, the enhanced oil recovery. The slow osmosis of water through thin brine films and thus the slow displacement of oil from the pore may help explain the anomalously slow oil recovery reported in micro-modeling experiments of LSW.
• #### A full multigrid multilevel Monte Carlo method for the single phase subsurface flow with random coefficients

(arXiv, 2019-10-10) [Preprint]
The subsurface flow is usually subject to uncertain porous media structures. In most cases, however, we only have partial knowledge about the porous media properties. A common approach is to model the uncertain parameters of porous media as random fields, then the statistical moments (e.g. expectation) of the Quantity of Interest(QoI) can be evaluated by the Monte Carlo method. In this study, we develop a full multigrid-multilevel Monte Carlo (FMG-MLMC) method to speed up the evaluation of random parameters effects on single-phase porous flows. In general, MLMC method applies a series of discretization with increasing resolution and computes the QoI on each of them. The effective variance reduction is the success of the method. We exploit the similar hierarchies of MLMC and multigrid methods and obtain the solution on coarse mesh $Q^c_l$ as a byproduct of the full multigrid solution on fine mesh $Q^f_l$ on each level $l$. In the cases considered in this work, the computational saving due to the coarse mesh samples saving is $20\%$ asymptotically. Besides, a comparison of Monte Carlo and Quasi-Monte Carlo (QMC) methods reveals a smaller estimator variance and a faster convergence rate of the latter approach in this study.
• #### Thermodynamically consistent modelling of two-phase flows with moving contact line and soluble surfactants

(Journal of Fluid Mechanics, Cambridge University Press (CUP), 2019-09-27) [Article]
Droplet dynamics on a solid substrate is significantly influenced by surfactants. It remains a challenging task to model and simulate the moving contact line dynamics with soluble surfactants. In this work, we present a derivation of the phase-field moving contact line model with soluble surfactants through the first law of thermodynamics, associated thermodynamic relations and the Onsager variational principle. The derived thermodynamically consistent model consists of two Cahn-Hilliard type of equations governing the evolution of interface and surfactant concentration, the incompressible Navier-Stokes equations and the generalized Navier boundary condition for the moving contact line. With chemical potentials derived from the free energy functional, we analytically obtain certain equilibrium properties of surfactant adsorption, including equilibrium profiles for phase-field variables, the Langmuir isotherm and the equilibrium equation of state. A classical droplet spread case is used to numerically validate the moving contact line model and equilibrium properties of surfactant adsorption. The influence of surfactants on the contact line dynamics observed in our simulations is consistent with the results obtained using sharp interface models. Using the proposed model, we investigate the droplet dynamics with soluble surfactants on a chemically patterned surface. It is observed that droplets will form three typical flow states as a result of different surfactant bulk concentrations and defect strengths, specifically the coalescence mode, the non-coalescence mode and the detachment mode. In addition, a phase diagram for the three flow states is presented. Finally, we study the unbalanced Young stress acting on triple-phase contact points. The unbalanced Young stress could be a driving or resistance force, which is determined by the critical defect strength.
• #### Darcy-scale phase equilibrium modeling with gravity and capillarity

(Journal of Computational Physics, Elsevier BV, 2019-09-05) [Article]
The modeling of multiphase fluid mixture and its flow in porous media is of great interest in the field of reservoir simulation. In this paper, we formulate a novel energy-based framework to model multi-component two-phase fluid systems at equilibrium. Peng-Robinson equation of state (EOS) is used to model the bulk properties of each phase, though our framework works well also with other equations of state. Our model reduces to the conventional compositional grading if restricted to one spatial vertical dimension together with the assumption of monodisperse pore-size distribution (all pores being one size). However, our model can be combined with a general distribution of pore size, which can generate interesting behaviors of capillarity in porous media. In particular, the model can be used to predict the capillary pressure of two-phase fluid as a function of saturation, with a given pore-size distribution. This model is the quantitative study of the first time in the literature for the capillarity of a two-phase fluid with partial miscibility. We proposed an unconditional-stable energy-decay numerical algorithm based on convex-concave splitting, which has been demonstrated to be both robust and efficient using numerical examples. To verify our model, we simulate the compositional grading of a binary fluid mixture consisting of carbon dioxide and normal decane. To demonstrate powerful features of our model, we provide an interesting example of fluid mixture in a porous medium with wide pore size distribution, where the competition of capillarity and gravity is observed. This work represents the first effort in the literature that rigorously incorporates capillarity and gravity effects into EOS-based phase equilibrium modeling.
• #### Heat and mass transfer in a viscous nanofluid containing a gyrotactic micro-organism over a stretching cylinder

(Symmetry, MDPI AG, 2019-09-05) [Article]
This work consists of a theoretical boundary layer analysis of heat and mass transport in a viscous fluid-embracing gyrotactic micro-organism over a cylinder. The flow governing equations are modeled through boundary layer approximations. The governing non-linear partial differential equations are lessened to a set of nonlinear ordinary differential equations using similitude transformation. The boundary layer equations are elucidated numerically, applying the spectral relaxation method with the aid of the computational software MATLAB. The impact of several pertinent parameters on flow convective characteristic phenomena are explored through the use of graphs and tables and are discussed with in-depth physical descriptions. In addition, the friction factor, the rate of heat transfer, rate of mass transfer, and the density number of the motile microorganism are also presented with respect to the above controlled parameters. It is noticed that for the increasing values of the magnetic parameter with reductions and enhancements, the density of the motile microorganism is a declining function of, and the concentration field enhances with the strengthening of, whereas it reduces with the rise of. Furthermore, the streamline patters are emphasized for the impact of controlled flow variables. Current outcomes are compared with the available results from previous cases and are observed to be in agreement.
• #### Controlling Factors of Degassing in Crosslinked Polyethylene Insulated Cables

(Polymers, MDPI AG, 2019-09-02) [Article]
Here, we analyze the degassing process of a byproduct (methane) formed during the peroxide-induced crosslinking of polyethylene. A diffusion model based on Fick’s law is used to obtain the controlling factors of degassing in a crosslinked polyethylene (XLPE) insulated power cable (132 kV with 18 mm of insulation). We quantitatively analyze different scenarios of the diffusion of methane through the XLPE insulation and two semiconductor layers under various in situ degassing conditions. The analyzed degassing conditions include heat transfer and its effect on the diffusion properties, the different transport and boundary conditions due to the free spaces within the cable conductor, and the nonuniform distribution of methane concentrations within the insulation layers. Our simulation results clearly demonstrate that the free spaces between the copper strands in the cable conductor significantly affect the degassing efficiency. However, the temperature-diffusion coupling has a relatively minor effect on the overall degassing efficiency due to the rapid temperature increase of the polymer layers during the initial stages of degassing. Moreover, we find that the nonuniform distribution of methane in the initial stages also plays an important role in degassing in the cable, but this effect varies significantly during the degassing process.
• #### Advances in Gaussian random field generation: a review

(Computational Geosciences, Springer Science and Business Media LLC, 2019-08-05) [Article]
Gaussian (normal) distribution is a basic continuous probability distribution in statistics, it plays a substantial role in scientific and engineering problems that related to stochastic phenomena. This paper aims to review state-of-the-art of Gaussian random field generation methods, their applications in scientific and engineering issues of interest, and open-source software/packages for Gaussian random field generation. To this end, first, we briefly introduce basic mathematical concepts and theories in the Gaussian random field, then seven commonly used Gaussian random field generation methods are systematically presented. The basic idea, mathematical framework of each generation method are introduced in detail and comparisons of these methods are summarized. Then, representative applications of the Gaussian random field in various areas, especially of engineering interest in recent two decades, are reviewed. For readers’ convenience, four representative example codes are provided, and several relevant up-to-date open-source software and packages that freely available from the Internet are introduced.
• #### Structure, Thermodynamics, and Dynamics of Thin Brine Films in Oil-Brine-Rock Systems.

(Langmuir : the ACS journal of surfaces and colloids, American Chemical Society (ACS), 2019-08-02) [Article]
Thin brine films are ubiquitous in oil-brine-rock systems such as oil reservoirs and play a crucial role in applications such as enhanced oil recovery. We report the results of molecular simulations of brine films that are confined between model oil (n-decane) and rock (neutral or negatively charged quartz slabs), with a focus on their structure, electrical double layers (EDLs), disjoining pressure, and dynamics. As brine films are squeezed to ∼0.7 nm (∼3 water molecule layers), the structures of the water-rock and water-oil interfaces change only marginally, except that the oil surface above the brine film becomes less diffuse. As the film is thinned from ∼1.0 to ∼0.7 nm, ions are enriched (depleted) near the rock (oil) surface, especially at a bath ion concentration of 0.1 M. These changes are caused primarily by the reduced dielectric screening of water and the weakened ion hydration near water-oil interfaces and, to a smaller extent, by the increased confinement. When the brine film is ∼1.0 nm thick, hydration and EDL forces contribute to the disjoining pressure between the charged rock and the oil. The EDL forces are reduced substantially as the ion concentration increases from 0.1 to 1.0 M, and the magnitude of the reduction is close to that predicted by the Poisson-Boltzmann equation. When the brine film is thinned from ∼1.0 to ∼0.7 nm, the disjoining pressure increases by ∼10 MPa, which is mostly due to an increase in the hydration forces. The first layer of water on the rock surface is nearly stagnant, even in 0.74 nm-thick brine films, whereas the viscosity of water beyond the first layer is bulk-like, and the slip coefficient of oil-water interfaces is close to that under unconfined conditions. The insights that are obtained here help lay a foundation for the rational application of technologies such as low-salinity waterflooding.
• #### Parallel reservoir simulators for fully implicit complementarity formulation of multicomponent compressible flows

(Computer Physics Communications, Elsevier B.V., 2019-07-30) [Article]
The numerical simulation of multicomponent compressible flow in porous media is an important research topic in reservoir modeling. Traditional semi-implicit methods for such problems are usually conditionally stable, suffer from large splitting errors, and may accompany with violations of the boundedness requirement of the numerical solution. In this study we reformulate the original multicomponent equations into a nonlinear complementarity problem and discretize it using a fully implicit finite element method. We solve the resultant nonsmooth nonlinear system of equations arising at each time step by a parallel, scalable, and nonlinearly preconditioned semismooth Newton algorithm, which is able to preserve the boundedness of the solution and meanwhile treats the possibly imbalanced nonlinearity of the system. Some numerical results are presented to demonstrate the robustness and efficiency of the proposed algorithm on the Tianhe-2 supercomputer for both standard benchmarks as well as realistic problems in highly heterogeneous media.