These formulations treat an arbitrary number of components. Reeves and Abriola [1994] proposed the set-iterative two-phase model described above. The mass-lumped finite-element method used variably-weighted time stepping, with linear kinetic mass exchange. Successive-substitution iteration solved the phase-equilibrium equations. Adenekan et al. [1993] (building on Falta et al. [1992]) and Sleep and Sykes [1993a] detailed three-dimensional three-phase n-component models that allowed every component to appear in every phase, with equilibrium mass transfer.
Adenekan et al. [1993] treated arbitrary phase appearance and disappearance by variable substitution without fictitious saturations. Gas properties were based on a Soave-Redlich-Kwong equation of state, ideal mixing in an incompressible oil phase was assumed, and water properties depended on pressure and temperature only. Advection by capillary, gravity, and viscous forces, diffusion, linear adsorption, and heat transfer were included. A fully implicit upstream-weighted integrated finite-difference discretization that allowed for implicit pressure-controlled wells yielded nonlinear equations solved by Newton-Raphson iteration, with a sparse direct or conjugate-gradient-based iterative linear-equation solver.
Sleep and Sykes [1993a,b] also used variable substitution in a complicated set of cases. Dispersion as well as diffusion was included. Products of phase saturations and dispersion coefficients were averaged harmonically, setting dispersive flux of a phase to zero if the phase did not exist on one side of a cell face. The nonlinear finite-difference equations were solved by Newton-Raphson iteration, with linear equations treated by Orthomin preconditioned by a D4-ordered ILU decomposition. Two numerical features of this model are notable. First, it allowed for fully implicit, IMPES, or adaptive implicit time stepping. The implicit/explicit switching criteria were based only on changes in state variables, so that implicit cells could not be allowed to become explicit. Stability-based criteria do not suffer from this drawback. Second, in explicit cells, second- and third-order upstream weighting schemes were offered as alternatives to the standard first-order approach. The goal was to reduce numerical dispersion, especially critical in NAPL applications because of the importance of the arrival time of low concentrations. It should be noted that the finite propagation speed of an IMPES formulation, advancing the leading edge of a front at most one cell per time step, can be helpful here. With 2 to 3 Newton-Raphson iterations for IMPES and 5 to 6 for fully implicit time steps, and with more-accurate weighting available in explicit cells, IMPES is likely to be more efficient in terms of overall accuracy per unit of processing time, an observation similar to that of Young and Russell [1993], where adaptive implicit in many cases was still more efficient.
Representative of the state of the art in compositional petroleum reservoir simulation is the report of Young [1992]. The three-dimensional generalized compositional formulation [ Young and Stephenson, 1983] separates flow and transport computations from fluid-property calculations, simulating different types of systems by substituting different fluid-property modules that can allow any number of components to flow in any phase. The mass-conservative finite-difference equations, allowing for bordered matrices resulting from implicit wells, are solved by an IMPES or adaptive implicit procedure. With vectorization of 98 to 99 percent of the computations, overall speeds of over 500 MFLOPS (500 million floating-point operations per second) have been achieved on supercomputers. Particular attention has been paid to vectorization of equation-of-state equilibrium calculations, which typically occupy over half of the computing time [ Young, 1991]. An IMPES run with 20,000 cells, 1,000 time steps, and 9 components was completed in 23 minutes on a Cray X-MP [ Young, 1991], and black-oil simulations with over 1,000,000 cells have since been performed in minutes [ Young, 1992]. As is typical for petroleum models, dispersion, sorption, reactions, and biodegradation are not included; the paper indicated how some of this could be done. Hinkley and Killough [1992] applied the Young-Stephenson formulation to three-phase groundwater flow and transport in the vadose zone, finding that gas-phase dynamics was important in their examples.
Recent innovations in compositional reservoir modeling primarily centered around more efficient or robust alternatives to traditional successive-substitution algorithms [e.g., Mehra et al., 1982] for phase equilibria. Litvak [1993] proposed direct minimization of Gibbs free energy, as opposed to methods that find a state where its first partial derivatives vanish. Sams et al. [1993] used equations based on scaling theory from condensed-matter physics to obtain a more robust algorithm in near-critical regions.