Pore-scale modeling of phase change in porous media

The combination of high-resolution visualization techniques and pore-scale flow modeling is a powerful tool used to understand multiphase flow mechanisms in porous media and their impact on reservoir-scale processes. One of the main open challenges in pore-scale modeling is the direct simulation of flows involving multicomponent mixtures with complex phase behavior. Reservoir fluid mixtures are often described through cubic equations of state, which makes diffuse-interface, or phase-field, theories particularly appealing as a modeling framework. What is still unclear is whether equation-of-state-driven diffuse-interface models can adequately describe processes where surface tension and wetting phenomena play important roles. Here we present a diffuse-interface model of single-component two-phase flow (a van der Waals fluid) in a porous medium under different wetting conditions. We propose a simplified Darcy-Korteweg model that is appropriate to describe flow in a Hele-Shaw cell or a micromodel, with a gap-averaged velocity. We study the ability of the diffuse-interface model to capture capillary pressure and the dynamics of vaporization-condensation fronts and show that the model reproduces pressure fluctuations that emerge from abrupt interface displacements (Haines jumps) and from the breakup of wetting films.

Figure 1

The van der Waals equation of state near the critical point. (a) Pressure isotherms as a function of density at various temperatures. Below the critical temperature the pressure is nonmonotonic, which reveals the presence of an unstable region that leads to phase separation and phase coexistence. In these pressure isotherms we indicate stable (blue line), metastable (black line), and unstable (red line) regions. (b) The Helmholtz free energy is nonconvex below the critical temperature, with a concave region that corresponds to the unstable region promoting spinodal decomposition. We use the same color code for stable, metastable, and unstable regions.

Figure 2

Phase equilibrium and phase diagram for a van der Waals fluid. (a) At a given temperature, states with homogeneous density are only possible at sufficiently large (liquid) or small (vapor) densities (blue lines). Densities lying inside the spinodal region (red line) are unstable, leading to phase separation. The equilibrium vapor and liquid densities in heterogeneous states, where liquid and vapor coexist, are obtained through the common-tangent construction in the specific free energy. The equilibrium states are commonly known as Maxwell states ${\rho }_l^M$ and ${\rho }_v^M$. The convex regions between the unstable region and the Maxwell states are metastable (black lines). (b) Phase diagram for a van der Waals fluid at different temperatures and average density. The shaded areas indicate metastable vapor (MV), metastable liquid (ML), and unstable states.

Figure 3

Model validation using the Young-Laplace equation for a liquid slug: pressure maps. We compute the equilibrium density field for a liquid bridge in a microchannel. We vary the static contact angle (a) $\theta ={30}^{\circ }$, (b) $\theta ={60}^{\circ }$, (c) $\theta ={120}^{\circ }$, and (d) $\theta ={150}^{\circ }$. The initial condition is such that the liquid bridge occupies the center of the domain at equilibrium. When the liquid phase is wetting $\left(\theta <{90}^{\circ }\right)$, the pressure inside the slug is lower than that of the surrounding vapor. Otherwise, the pressure inside the liquid bridge is larger than the pressure of the vapor.

Figure 4

Model validation using the Young-Laplace equation for a liquid slug. Using the same simulation setup as in Fig. 3, we measure the pressures inside the liquid and vapor phases, at points sufficiently far from the interface, and plot the pressure difference between the vapor and liquid phases, the capillary pressure $\mathrm{\Delta }{p}_r=p_{rv}-p_{rl}$, against $cos\theta$. The linear dependence between capillary pressure and $cos\theta$ indicates that the model honors the Young-Laplace equation.

Figure 5

The interface thickness, in a diffuse-interface model driven by a cubic equation of state $\delta$, is controlled by the Korteweg number $\mathrm{\Lambda }$. We compute the equilibrium density field for a liquid bridge in a microchannel with $\theta ={30}^{\circ }$ and different values of $\mathrm{\Lambda }$. We plot the capillary pressure against $\mathrm{\Lambda }$ and observe that our numerical simulations recover the theoretical scaling $\delta \sim \sqrt{\mathrm{\Lambda }}$. In practice, this strong dependence of surface tension on interface thickness limits the use of these direct simulations to the micrometer scale, as the numerical discretization needs to be fine enough to capture the interface thickness that is consistent with the correct surface tension.

Figure 6

Capillary pressure affects the equilibrium states of liquid and vapor densities. We study the liquid and vapor densities at equilibrium in the context of the liquid bridge simulations analyzed in Figs. 3, 4, 5. (a) For a relatively small Korteweg number $\mathrm{\Lambda }={10}^{-4}$, the static contact angle only slightly alters the equilibrium densities, which remain close to the theoretical Maxwell states. (b) The increase in surface tension as $\mathrm{\Lambda }$ increases promotes much larger deviations from the equilibrium values and the liquid densities enter into the metastable region. The van der Waals pressure curve at this temperature $T_r=0.7$ is plotted as a reference (black line).

Figure 7

Schematic of the validation problem simulating evaporation of a van der Waals fluid in a channel.

Figure 8

Evaporation in a channel. The channel height is $H=0.1$ units, and a density ${\rho }_v^{\text{out}}=0.07$ is imposed at the left boundary. The initial density is ${\rho }_0=2.2$. (a) Reduced density profiles at time $t/t^{*}=10$ for three static contact angles $\theta ={45}^{\circ },{90}^{\circ }$, and ${135}^{\circ }$. (b) Comparison between numerical and analytical vapor fluxes (41) through the left boundary.

Figure 9

Phase separation of a van der Waals fluid in porous media: density maps at dimensionless time $t_D=50$, showing the arrest of coarsening due to the pore geometry and wetting conditions. We show results for the strongly wetting case $\theta =0^{\circ }$ and mean densities (a) $\overline{\rho }=0.75095$, (b) $\overline{\rho }=0.90152$, (c) $\overline{\rho }=1.1525$, and (d) $\overline{\rho }=1.5038$. The temperature is $T_r=0.7$ and the Korteweg number is $\mathrm{\Lambda }=2\times {10}^{-4}$. The initial density fields are random perturbations around the mean density.

Figure 10

Phase separation of a van der Waals fluid in porous media: pressure maps at dimensionless time $t_D=50$, showing the arrest of coarsening due to the geometry and wetting conditions. We show results for the strongly wetting case $\theta =0^{\circ }$ and mean densities (a) $\overline{\rho }=0.75095$, (b) $\overline{\rho }=0.90152$, (c) $\overline{\rho }=1.1525$, and (d) $\overline{\rho }=1.5038$. The temperature is $T_r=0.7$ and the Korteweg number is $\mathrm{\Lambda }=2\times {10}^{-4}$. Note that the final liquid pressures (and therefore densities) vary with the final liquid volume fraction.

Figure 11

Phase separation of a van der Waals fluid in porous media: metastable density states after long-time simulations of phase separation. We plot the equilibrium densities in the van der Waals pressure-density curves. The simulations correspond to those of Figs. 9 and 10 plus other additional cases not shown above. Note that the key driver behind the nonequilibrium densities is the capillary pressure, controlled by surface tension and pore size distribution rather than the mean density or initial conditions.

Figure 12

Simulation of a vaporizing front in a micromodel. The internal domain walls are strongly wetting to the liquid phase $\left(\theta =0^{\circ }\right)$. A closed porous domain is initially filled with liquid of density ${\rho }_{r0}=2.2$ at temperature $T_r=0.7$. We force phase transition by decreasing the pressure at the left boundary, while the others remain no-flow walls. The low boundary pressure, enforced by imposing a density ${\rho }_{rLB}=0.07$, triggers a phase change front where liquid vaporizes at the liquid-vapor interface and flows out of the domain.

Figure 13

Simulation of a vaporizing front in a micromodel. The internal domain walls are strongly wetting to the vapor phase $\left(\theta ={180}^{\circ }\right)$. A closed porous domain is initially filled with liquid of density ${\rho }_{r0}=2.2$ at temperature $T_r=0.7$. We force phase transition by decreasing the pressure at the left boundary, while the others remain no-flow walls. The low boundary pressure, enforced by imposing a density ${\rho }_{rLB}=0.07$, triggers a phase change front where liquid vaporizes at the liquid-vapor interface and flows out of the domain.

Figure 14

Time evolution of mass fluxes for different static contact angles during the advance of a vaporizing front in a micromodel (including the simulations in Figs. 12 and 13). These results indicate that the simulated phase change fronts can be understood as simple diffusion-limited processes.

Figure 15

Vaporizing front patterns are controlled by heterogeneity of the pore space. We show snapshots of a simulation with the same conditions as those of Fig. 13 $({\rho }_{r0}=2.2,{\rho }_{rLB}=0.07,T_r=0.7$, and $\theta =0^{\circ })$ but with a distribution of obstacles that slightly promotes imbibition along the top and bottom boundaries. These patterns of vapor invasion resemble the invasion-percolation-like fronts in capillary fingering, as proposed in the context of drying [83, 84].

Figure 16

Vaporizing front patterns under mixed-wet conditions. We show snapshots of a simulation with the same conditions as those of Fig. 13 $({\rho }_{r0}=2.2,{\rho }_{rLB}=0.07$, and $T_r=0.7)$ but with a distribution of circular and rectangular obstacles. In this case the quadrilateral obstacles are strongly wetting to the vapor phase $\theta ={180}^{\circ }$, while the circular ones are strongly wetting to the liquid phase $\theta =0^{\circ }$. See also video 1 in the Supplemental Material [99].

Figure 17

Vaporizing front patterns under mixed-wet conditions. We show snapshots of a simulation with the same conditions as those of Fig. 13 $({\rho }_{r0}=2.2,{\rho }_{rLB}=0.07$, and $T_r=0.7)$ but with a distribution of circular and rectangular obstacles. In this case the quadrilateral obstacles are strongly wetting to the liquid phase $\theta =0^{\circ }$, while the circular ones are strongly wetting to the vapor phase $\theta ={180}^{\circ }$. See also video 2 in the Supplemental Material [99].

Figure 18

Pressure maps at the end of the simulations (dimensionless time $t_D=12)$ for the two mixed-wetting cases shown in Figs. 16 and 17. The pockets of residual liquid remain at metastable conditions at very low pressures.

Figure 19

Time evolution of mass fluxes for the two mixed-wetting cases shown in Figs. 16 and 17 during the advance of a vaporizing front in a micromodel. These results indicate that the simulated phase change fronts can be understood as simple diffusion-limited processes.

Figure 20

Pressure and density fluctuations associated with the advancement of interfaces and the rupture of liquid films. We show the time evolution of density and pressure at a location at the surface of an obstacle (indicated with a red dot), together with snapshots of the density field immediately before and after the front has reached this location. The model parameters are ${\rho }_{r0}=2.2,{\rho }_{rLB}=0.07,T_r=0.7$, and $\theta ={180}^{\circ }$. The Maxwell densities ${\rho }_{M_1}$ and ${\rho }_{M_2}$, which are the thermodynamic equilibrium densities of vapor and liquid at this temperature, are indicated with dashed blue lines. The initial, left-boundary, and Maxwell pressures $p_0,p_{BC}$, and $p_M$, respectively, are also indicated with dashed blue lines. The large pressure perturbations as the interface advances correspond to the propagation of pressure waves generated when interfaces move or break. In this imbibition case, the transition densities correspond almost exactly to those of theoretical equilibrium (indicated with dashed blue lines). The pressure is higher than the theoretical pressure on the liquid side of the interface $p_M$ due to surface tension.

Figure 21

Pressure and density fluctuations associated with the advancement of interfaces and the rupture of liquid films. We show the time evolution of density and pressure at a location at the surface of an obstacle (indicated with a red dot), together with snapshots of the density field immediately before and after the front has reached this location. The model parameters are ${\rho }_{r0}=2.2,{\rho }_{rLB}=0.07,T_r=0.7$, and $\theta ={90}^{\circ }$. The Maxwell densities ${\rho }_{M_1}$ and ${\rho }_{M_2}$, which are the thermodynamic equilibrium densities of vapor and liquid at this temperature, are indicated with dashed blue lines. The initial, left-boundary, and Maxwell pressures $p_0,p_{BC}$, and $p_M$, respectively, are also indicated with dashed blue lines. In this neutrally wetting case, pressures oscillate around the theoretical transition pressure $p_M$ without major differences between the liquid and vapor pressures (zero capillary pressure).

Figure 22

Pressure and density fluctuations associated with the advancement of interfaces and the rupture of liquid films. We show the time evolution of density and pressure at a location at the surface of an obstacle (indicated with a red dot), together with snapshots of the density field immediately before and after the front has reached this location. The model parameters are ${\rho }_{r0}=2.2,{\rho }_{rLB}=0.07,T_r=0.7$, and $\theta =0^{\circ }$. The Maxwell densities ${\rho }_{M_1}$ and ${\rho }_{M_2}$, which are the thermodynamic equilibrium densities of vapor and liquid at this temperature, are indicated with dashed blue lines. The initial, left-boundary, and Maxwell pressures $p_0,p_{BC}$, and $p_M$, respectively, are also indicated with dashed blue lines. In this case of drainage, pressures inside the vapor phase are significantly higher than in the liquid phase, which is consistent with a large capillary pressure.

Figure 23

Pressure and density fluctuations associated with the advancement of interfaces and the rupture of liquid films. A summary of the time evolution results shows the differences in pressure behavior with the wetting properties of the system.

Figure 24

Propagation of pressure fluctuations generated by the advance of the vaporizing front. We show (a)–(d) density maps at different time levels and (e)–(h) the corresponding pressure maps, illustrating the propagation of a pressure pulse generated by the breakup of a liquid film at the front.