A many-body dissipative particle dynamics model, namely, MDPD, is applied for simulation of pore-scale, multi-component, multi-phase fluid flows in fine-grained, nanoporous shales. Since this model is able to simultaneously capture the discrete features of fluid molecules in nanometer size pores and continuum fluid dynamics in larger pores, and is relatively easy to parameterize, it has been recognized as being particularly suitable for simulating complex fluid flow in multi-length-scale nanopore networks of shales. A remarkable feature of this work is the integration of a high-resolution FIB-SEM (focused ion beam scanning electron microscopy) digital imaging technique to the MDPD model for providing 3D voxel data that contain the invaluable geometrical and compositional information of shale samples. This is the first time that FIB-SEM is seamlessly linked to a Lagrangian model like MDPD for fluid flow simulation, which offers a robust approach to bridging gaps between the molecular- and continuum-scales, since the relevant spatial and temporal scales are too big for molecular dynamics, and too small for computational fluid dynamics with known constitutive models. Simulations ranging from a number of benchmark problems to a forced two-fluid flow in a Woodford shale sample are presented. Results indicate that this model can be used to deliver reasonable simulations for multi-component, multi-phase fluid flows in arbitrarily complex pore networks in shales.

Approximately 75% of the sedimentary rocks on earth are clastic nanoporous rocks, including compacted and diagenetically altered sands (tight sands) and a variety of rocks formed from fine grained sediments such as mudstones, siltstones, claystones, and marlstones, which are often referred to as “shales.” These shales contain most of the world’s oil and possibly most of the natural gas resources as well, according to a recent estimation (Kuuskraa et al., 2013). However, only a small fraction can be recovered using the latest technology—typically about 5% for oil and about 20% for gas (Kuuskraa et al., 2013). This low recovery is at least in part due to the lack of the knowledge and understanding of the relevant fundamental physics that ultimately control the dynamics of fluids organic-rich shales with extremely low permeabilities in the micro- to nano-Darcy range and average pore sizes from a few nanometers (10−9 m) to a few micrometers (10−6 m). Filling these knowledge gaps could potentially facilitate the development and deployment of more effective recovery strategies from tight shales.

Most of the theoretical understanding of fluid flow in porous geomaterials (and the predictive models built upon such understanding) has been based on the concepts of classical continuum fluid dynamics and a rigid porous or fractured solid porous matrix, which determines the “ideal” non-slip boundary conditions (BCs) for fluid flow and transport (Bear, 1972). For most of the part, those concepts and models have proven to be adequate for developing a fundamental understanding of a single- and multi-phase fluid flow in permeable porous media such as aquifers, soils, and conventional oil and gas reservoirs (Bertoncello and Honarpour, 2013). Many pore-scale models have been successfully developed in either Eulerian or Lagrangian frame, based on such continuum computational fluid dynamics (CFD). To list a number of recent works, Pan et al. (2001; 2004) applied a lattice Boltzmann (LB) and pore-network modeling approach to simulate a single- and two-phase flow in complex pore networks; Tartakovsky and Meakin (2005; 2006) introduced a smoothed particle hydrodynamics (SPH) model for immiscible and miscible flow simulations in fractured porous media, and Huang et al. (2005a; 2005b) applied a finite volume method based on the volume-of-fluid (VOF) model for multiphase flow simulations in fracture aperture and fracture junction with complex geometries. However, the behavior of fluids in organic-rich, fine-grained, nanoporous shales is very different from that in coarse-grained porous media. For example, in an idealized spherical pore with a diameter of 100 nm, ∼6% of the fluid is within a distance of 1 nm from the solid surface, and in a pore with diameter of 10 nm, ∼49% is within a distance of 1 nm; the properties of fluids in this thin surface layer are expected to be physically and chemically different from those of bulk fluids. Thus in organic-rich, fine-grained, nanoporous shales, the discreteness of atoms and molecules starts to become important for flow and transport processes at smaller scales; solid organic materials play an important role as mechanical components, sorbents and sources of fluids, and large specific surface areas make surface reactions and surface transport more important. Therefore, understanding large-scale fluid flow and transport behaviors within organic-rich, fine-grained, nanoporous shales requires robust and accurate multiscale computational models that can bridge both the spatial and temporal scales between fluid molecular dynamics and pore-scale flow.

Originally introduced by Hoogerbrugge and Koelman (1992) for microscopic hydrodynamics, the dissipative particle dynamics (DPD) models constitute a relatively new class of mesoscale models that can be used to simulate a single- and a multi-phase fluid flow, e.g., see Tiwari and Abraham (2006), Heldele et al. (2006), Visser et al. (2006), and Liu et al. (2006). A theoretical basis for DPD based on statistical mechanics was discussed in Español and Warren (1995) and Marsh (1998). A recent review of DPD theory, methodology, recent developments, and applications can be found in Moeendarbary et al. (2009) and Liu et al. (2015). Using DPD, a complex system can be simulated with a set of interacting particles, where each particle represents a small cluster of atoms or molecules instead of a single molecule (Figure 1). The particle-particle interaction force in a DPD embodiment consists of a “conservative” (non-dissipative) component, a random (fluctuating) component, and a dissipative component that represents the effect of viscosity. The inclusion of random particle-particle interactions in addition to the conservative particle-particle interactions is the most remarkable feature that distinguishes DPD from the SPH method. The random force generates fluctuations that correspond to thermally driven fluctuation in real fluid molecules. Such fluctuations are only important on microscopic physical scales, e.g., pores of sizes in nanometer ranges. Conversely, the DPD fluid can reproduce the continuum Navier-Stokes equations only on large scales (scales much greater than the particle size) with the effect of thermal fluctuations being negligible. This makes DPD essentially a mesoscale method between the molecular and continuum hydrodynamic scales, and facilitates simulations of complex fluid systems scales spanning a wide range of length and time scales. DPD rigorously conserves the number of particles (equivalently the mass) and momentum of the system. In addition, much larger time steps are possible in DPD simulations than in MD simulations. This is because the DPD particle-particle interactions are much softer than those (e.g., Lennard-Jones potential) experienced in MD. For these reasons, DPD is more efficient than MD for mesoscale flow simulations in terms of computational costs. Among the many existing DPD models, the “many-body” DPD model by Warren (2003), namely, MDPD, has been recognized as being particularly suitable for multi-fluid systems and relatively simple to parameterize for specific types of fluids. This MDPD model has been applied in simulations of liquid-vapor interface, surface tension, and multi-component fluid flows in micro-scale channels, e.g., Pan (2010), Chen et al. (2011), Ghoufi and Malfreyt (2011), Chen et al. (2013), and Chen et al. (2014). More recently, Liu et al. (2007a; 2007b) successfully applied a DPD model using an SPH-like conservative interaction force to simulate unsaturated flows in porous media. However in their works, only pair-wise particle-particle interaction potential and single phase were considered.

FIG. 1.

A conceptual schematic of coarse graining from molecules to mesoscale level DPD particles. Each DPD particle represents a collection of molecules.

FIG. 1.

A conceptual schematic of coarse graining from molecules to mesoscale level DPD particles. Each DPD particle represents a collection of molecules.

Close modal

The major effort of the present work is the application of a robust MDPD model for simulations of pore-scale, multi-component, multi-phase fluid flows within fine-grained, nanoporous shales. This is a relatively unexplored area so far. The novelties of the work lie in the following two aspects. First, this model is able to and has been applied to capture simultaneously the discrete features of fluid molecules in nanometer size pores and continuum fluid dynamics in larger pores. An example in Figure 2 demonstrates the multi-length-scale modeling capability of MDPD, where the same particle-particle interaction force and parameters have been used to model fluid-fluid interfaces in two cases: (a) at the micrometer-scale where the interface is relatively stable; and (b) at the nanometer scale where the interface is diffused due to the Brownian motion of fluid particles. Second, an FIB-SEM (focused ion beam scanning electron microscopy) high-resolution digital imaging technique (Goral and Miskovic, 2015) has been utilized to provide the MDPD model with 3D voxel data that contain the invaluable micro-geometrical and compositional information of shale samples. The process of the digital imaging process is briefly discussed in the  Appendix for clarity. The digitized microstructural information is used to construct pore networks with solid boundaries that enclose pores. To the best of the authors’ knowledge, this is the first time that FIB-SEM has been seamlessly integrated to a Lagrangian particle transport model like MDPD for modeling fluid flows in real multi-length-scale nanoporous shale samples, as compared with earlier DPD fluid flow simulations where only manufactured (or simulated) pore networks were used (Liu et al., 2007b). It is worth of mentioning that although the FIB-SEM machinery has become a popular approach for characterizing shale rock samples since the works by Curtis et al.2010; 2011; 2012a; 2012b, most of the earlier flow and transport simulations in shales were based on continuum CFD models. For example, Dewers et al. (2012) reported a finite element simulation of shale gas transport in Hayneville Formation mudstones by assuming sufficient interconnected large pores in the samples and that Darcy’s Law was valid. Unfortunately, those models are not appropriate for heterogeneous nanoporous shale rocks any more. In comparison, the integrated research method to be presented in this work offers a robust approach to bridging gaps between the molecular- and continuum-scales, as the relevant spatial and temporal scales are too big for molecular dynamics, and too small for continuum CFD with known constitutive models.

FIG. 2.

This example shows the multi-length-scale modeling capability of the MDPD model, where fluid-fluid interfaces between a non-wetting fluid (yellow) and a wetting fluid (blue) in pore space can be modeled by MDPD in (a) a 0.1-μm wide pore and (b) a 2-nm wide pore. The same particle-particle interaction force and parameters were used in these two cases.

FIG. 2.

This example shows the multi-length-scale modeling capability of the MDPD model, where fluid-fluid interfaces between a non-wetting fluid (yellow) and a wetting fluid (blue) in pore space can be modeled by MDPD in (a) a 0.1-μm wide pore and (b) a 2-nm wide pore. The same particle-particle interaction force and parameters were used in these two cases.

Close modal

The rest of the paper is organized as follows. First, the MDPD model is presented in Section II. Then a number of validation test cases are reported in Section III. A pioneering flow simulation in a pore network enclosed in an organic-rich, nanoporous shale is presented in Section IV. Finally, the concluding remarks are given in Section V.

In a generic DPD model, particles interact via pairwise central forces, i.e., 𝐅ij=𝐅ijR+𝐅ijD+𝐅ijC, where 𝐅ijR represents a random force, 𝐅ijD a dissipative force, and 𝐅ijC a conservative force between particle i and j, respectively. If ri and vi are used to denote the position and velocity of particle i, respectively, the random force, 𝐅ijR, and the dissipative force, 𝐅ijD, can be expressed as 𝐅ijR=σwR(rij)ξij𝐫^ij, and 𝐅ijD=γwD(rij)(𝐫^ij𝐯ij)𝐫^ij, where rij = rirj, rij=|𝐫ij|, 𝐫^=𝐫ij/rij,vij = vivj. These forces constitute a thermostat if the amplitude, σ, of the random variable, ξij, and the viscous dissipation coefficient, γ, satisfy a fluctuation-dissipation theorem: σ2=2γkBT and wD(r)=(wR(rij))2, where kBT denotes the desired temperature in the unit of Boltzmann’s constant kB. In the original DPD model, the conservative force, 𝐅ijC, is defined as 𝐅ijC=aijwC(rij)𝐫^ij, where aij denotes the magnitude of the force, and the weight function wC(r) vanishes when the inter-particle distance r is larger than a cutoff range rc. 𝐅ijC is usually derived from a soft and unspecific weight function, wC(rij), thus allowing for a fairly large integration time step. Different weight functions describe different material properties. A common choice for wC(rij) is wC(rij) = 1 − rij/rc, and wR=wC. The standard velocity Verlet algorithm (Revenga et al., 1998) is employed to integrate the resulting equations of motion. A quadratic equation of state (EOS) was obtained with respect to the average particle density ρ, as shown in Figure 3(a). However, the original DPD model is not sufficient to model multi-species or multi-phase fluid flow phenomena such as those associated with liquid-vapor interfaces, liquid-liquid interfaces, and free capillary surfaces. A more complex EOS needs to be represented with the DPD method. To achieve this, a long-range attractive, short-range repulsive conservative force, FC, is required. The approach employed in the present work is the so-called many-body DPD method (Warren, 2003), namely, MDPD. The 𝐅ijC in the standard DPD method is augmented by density-dependent contributions, and the resulting method includes the van der Waals loop in the EOS, as shown in Figure 3(b). In the MDPD model, the conservative force is expressed as

𝐅ijC=AijwC(rij)𝐫^ij+Bij(ρ¯i+ρ¯j)wd(rij)𝐫^ij,
(1)

which consists of a long-range attractive part that is density-independent, and a short-range repulsive part that depends on a weighted average of the local particle density. The attractive component, AijwC(rij)𝐫^ij, can be obtained by simply turning the sign of the original force parameter aij, i.e., Aij < 0, with a cutoff range rC = 1, and Bij(ρ¯i+ρ¯j)wd(rij)𝐫^ij is a many-body repulsive component with Bij > 0, and shorter cutoff wd(rij) = 1 − r/rd, where rd < rC. The averaged local density, ρ¯i, at the position of particle i can be computed as ρ¯i=jiwρ(rij), where the normalized weight function, wρ, needs to satisfy 04πr2wρ(r)dr=1. wρ is defined as follows for a three-dimensional computational domain: wd(r)=152πrd3(1rrd)2.

FIG. 3.

Validation of the EOS: (a) p=ρkBT+0.1aijrC4ρ2 for the original DPD model with kBT = 1, γ=4.5, rC = 1, and aij = 25; (b) p=ρkBT+αAijρ2+2αBijrd4(ρ3cρ2+d) for the MDPD model with kBT = 1, γ=4.5, rC = 1, rd = 0.75, Aij = −40, Bij = 25, α=0.101, c = 4.16, and d = 18. Pressure for each particle number density ρ was obtained by averaging over 1000 time steps after equilibrium, in a 10 × 10 × 10 periodic box.

FIG. 3.

Validation of the EOS: (a) p=ρkBT+0.1aijrC4ρ2 for the original DPD model with kBT = 1, γ=4.5, rC = 1, and aij = 25; (b) p=ρkBT+αAijρ2+2αBijrd4(ρ3cρ2+d) for the MDPD model with kBT = 1, γ=4.5, rC = 1, rd = 0.75, Aij = −40, Bij = 25, α=0.101, c = 4.16, and d = 18. Pressure for each particle number density ρ was obtained by averaging over 1000 time steps after equilibrium, in a 10 × 10 × 10 periodic box.

Close modal

Because of the soft interaction between DPD particles, fluid particles may penetrate through a solid matrix when DPD and MDPD models are used for pore-scale fluid flow simulations. Such penetration is not physical and must be avoided. Most of the previous efforts to properly model solid-fluid interfaces in DPD fluid flow simulations were focused on imposing rigorous macroscopic boundary conditions (BCs), e.g., a non-slip BC at sharply defined impenetrable solid surfaces. That motivation originated from a strict mesoscale interpretation of DPD models, where a single DPD fluid particle represents a cluster of fluid molecules on scales well above the atomistic scale (Groot and Warren, 1997). In order to make the non-slip BC possible, additional forces must be exerted on fluid particles at the vicinity of solid-fluid interfaces. Parameterization must be carried out very carefully in order to avoid undesired spurious behavior, e.g., artificial slip (Pivkin and Karniadakis, 2005), oscillations in temperature (Revenga et al., 1999), and particle layering (Pivkin and Karniadakis, 2006). The disadvantages of those specialized boundary models is that they are only applicable to mostly flat, or spherically curved solid surfaces, e.g., Liu et al. (2007b). To overcome the difficulty of the rigorous non-slip BC, Henrich et al. (2007) devised a boundary model, where a weak external repelling force is imposed on those fluid particles whenever they penetrate into a thin layer in the solid boundary. This prevents fluid particles from further diffusing into solid matrix indefinitely.

However, in nanometer-sized pores, surface slip might be an important transport mechanism as well. Thus in order to account for the multi-length-scale nature of pores in shales, a boundary model is desired to model the continuum no-slip boundary in large pores, the non-continuum discrete motions of fluid particles in small pores, and the transition between the non-continuum and continuum fluid motions.

In this work, a simple boundary model is devised for modeling solid-fluid interfaces in the MDPD model. In this model, a dense layer of solid particles is defined with the MDPD particle attributes, and used to enclose pore space. Distribution of the solid particles is configured to be sufficiently dense to prevent penetration of the MDPD fluid particles that flow in pore space. This is equivalent to the application of a weak external repelling force approach (Pivkin and Karniadakis, 2006). The coordinates of solid particles are determined using 3D voxel data that are obtained through a digital rock imaging process, and these coordinates are fixed during simulations. For clarity, an example of data conversion from a 2D pixel map to coordinates of solid particle is illustrated in Figure 4. Later in Section IV, 3D voxel data of a shale sample obtained from the FIB-SEM digital rock imaging process will be used. For MDPD simulations, Eq. (1) is applied for (a) fluid-fluid, (b) solid-solid, and (c) solid-fluid particle interactions. However, the coefficients in Eq. (1) must be parameterized for each of the three interaction modes, so that fluid properties at the pore scale under consideration can be properly represented, e.g., liquid-vapor interface, miscibility between two types of fluids, and fluid wettability against the solid surface. Notice that this is not a completely new approach for modeling solid boundaries in DPD-type methods; for example, in their work, Pivkin and Karniadakis (2006) mentioned a similar idea, with which though they did not demonstrate a test case. Using the present approach, however, there is a potential artifact of undesired minor temperature oscillations in the vicinity of solid-fluid interfaces. Nevertheless, this simple boundary model is probably so far the only viable way for MDPD to model arbitrarily complex networks of enclosed pore space in nanoporous shales. As already shown in Figure 2, this boundary model is able to render a diffusive solid-fluid interface in a nanometer-scale pore, and also a stable solid-fluid interface in a micrometer-scale pore.

FIG. 4.

This example illustrates the data conversion from pixel map to particle coordinates. (a) First, pixels data obtained from digital imaging process are read in. Pixel cells that represent the locations of solids and pore space are labeled “1” (grey cells) and “0” (white cells), respectively. (b) Second, solid cells that are neighbors to pore cells are tagged; see the black dots. User-defined cell refinement is provided, e.g., a level-1 cell refinement means that each of the tagged solid cells is split into four. (c) Finally, central coordinates of tagged solid cells are saved. They will be used to represent solid walls in pores.

FIG. 4.

This example illustrates the data conversion from pixel map to particle coordinates. (a) First, pixels data obtained from digital imaging process are read in. Pixel cells that represent the locations of solids and pore space are labeled “1” (grey cells) and “0” (white cells), respectively. (b) Second, solid cells that are neighbors to pore cells are tagged; see the black dots. User-defined cell refinement is provided, e.g., a level-1 cell refinement means that each of the tagged solid cells is split into four. (c) Finally, central coordinates of tagged solid cells are saved. They will be used to represent solid walls in pores.

Close modal

The MDPD model described above has been programmed in an in-house code based on the LAMMPS framework (Plimpton, 1995). Before the code could be reliably applied for problems of practical interest, code validation was performed in this section to validate the implementation and key capabilities of the model for simulating multi-component, multi-phase fluid flows. Three test cases are presented, with the particle mass, m, temperature, kBT and many-body attraction cutoff range, rC, chosen to be the reduced units of mass, energy, and length, respectively, i.e., m = (kBT) = rC = 1. In the simulations, a constant time step dt = 0.01 in the reduced time unit along with the constant-energy, constant-volume NVE ensembles was used. The gravity force was neglected.

The objective of this test case is to validate if the implemented MDPD code accurately calculates properties of a specific type of fluid. We followed Ghoufi and Malfreyt (2011) to simulate a water liquid-vacuum interface, where the water density and surface tension were calculated. A 3D computational domain was defined, extending from −10rc to 10rc in the x-direction, from 0 to 5rc in the y-direction, and from 0 to 5rc in the z-direction. A periodic BC was used in all the directions. 1000 DPD particles were then randomly placed in x = [−5rc, 5rc]. The parameters in Eq. (1) are set as Aij = −50, Bij = 25, and rd = 0.75rC for modeling the water (Ghoufi et al., 2010). In this case, one DPD particle approximately represents a cluster of three water molecules, i.e., Nm = 3, and the size of one DPD particle is about 90 Å3. More details of the correspondence between the DPD reduced units and physical values can be found in Ghoufi and Malfreyt (2011). To calculate the thermodynamic properties, 1 × 106 time steps were carried out to equilibrate the system. Figure 5 displays an instantaneous snapshot of the system after equilibrium. Finally, 1 × 106 additional time steps were then conducted to calculate the time-averaged properties.

FIG. 5.

An instantaneous liquid-vacuum interface simulated from an MDPD simulation by 1000 particles at equilibrium. The x-axis is normal to the interface.

FIG. 5.

An instantaneous liquid-vacuum interface simulated from an MDPD simulation by 1000 particles at equilibrium. The x-axis is normal to the interface.

Close modal

Figure 6 shows a time-averaged density profile as a function of the x-axis. The density near x = 0 (the center of the water drop) is 6.88, matching the value reported in an MD study (Ghoufi and Malfreyt, 2011). Due to the simple geometry in this problem, the water-vacuum interfacial tension, γWV, can be calculated by subtracting the mean tangential stresses, σyy and σzz, from the normal stress, σxx:γWV=Lxσxx12(σyy+σzz). The calculated γWV is 12.4, matching the value reported in Ghoufi and Malfreyt (2011). The values for water density and water-vacuum interfacial tension can be converted into the physical units with the equations: rc=rc*(ρ*NmV)13Å, ρ=ρ*NmMNarc3kgm3, and γ=γ*kBTrc2Nm1, where the superscript * denote values in the DPD reduced unit, V is the volume of one water molecule (30 Å), M is the molar weight of a water molecule (18 g mol−1), Na is Avogadro’s number, and kB is Boltzmann’s constant, and T is 298 K. In the physical units, ρ=994kgm3 and γ=70.6×106Nm1, which agree well with the MD results in Ghoufi and Malfreyt (2011). This indicates that the implemented MDPD code delivers accurate predictions of thermodynamic properties for fluids of interest.

FIG. 6.

Time-averaged density of water versus x-axis calculated from an MDPD simulation by 1000 particles.

FIG. 6.

Time-averaged density of water versus x-axis calculated from an MDPD simulation by 1000 particles.

Close modal

Interface stability and interfacial tension between two types of fluids are crucial properties for characterizing fluid-fluid miscibility. To model and simulate multi-fluid flows, it is desirable to have a simple approach for independently calculating the interfacial tension between liquid-liquid pairs. In this numerical experiment, we set up the simulation problem that is shown in Figure 7. It contains an 11 × 10 × 10 non-wetting liquid slab with 4000 particles, and an 11 × 10 × 10 wetting liquid slab with 4000 particles. These slabs are initially attached to each other. The total tension of the system, γtotal, consists of three components: the wetting liquid surface tension, γWV, the non-wetting liquid surface tension, γNV, and the interfacial tension between the wetting and non-wetting fluids, γWN. For this nomenclature, the subscripts “W,” “N,” and “V” denote the wetting fluid, non-wetting fluid, and vacuum, respectively. γWN can be calculated by subtracting the mean tangential stresses, σyy and σzz, from the normal stress, σxx, i.e., γtotal=Lxσxx12(σyy+σzz), where Lx is the length of the simulation box in the x direction. Since γWV and γNV remain unchanged, γWN can then be calculated by subtracting the two surface tensions from the total tension, γtotal, i.e., γWN=γtotalγWVγNV, where γWV=γNV=7.67 is obtained when there is only one fluid type in the system.

FIG. 7.

Schematic of the fluid-fluid interfacial tension test case at the initial status.

FIG. 7.

Schematic of the fluid-fluid interfacial tension test case at the initial status.

Close modal

The MDPD parameters used in this simulation are as follows: AWW = −40, BWW = 25, and rd = 0.75rC between the wetting fluid particles; ANN = −40, BNN = 25, and rd = 0.75rC between the non-wetting fluid particles; and BWN = 25 and rd = 0.75rC between the wetting and non-wetting fluid particles. By changing AWN, different γWN values can be obtained, leading to different miscibility between the two liquids; see Figure 8. By increasing the AWN value from −40 to −5, that is, turning the attractive interaction between the wetting and non-wetting liquids from strong to weak, the liquid–liquid interface changes accordingly from being completely miscible to entirely immiscible. Notice that before AWN approaches 0, e.g., when AWN = −5, the two liquids already tend to detach from each other since the inter-liquid attraction becomes too weak between the two liquids. A correlation between γWN and AWN is shown in Figure 9. This matches the result in Chen et al. (2011). In general, successively altering the attraction parameter, A, while fixing the repulsion parameter, B, is an adequate approach for tuning the interaction of concern. Similarly, a solid–liquid interaction can also be tuned by adjusting the attraction parameter, A, from −40 to 0, which will result in an increase of solid–liquid interfacial tension as well as a decline of wettability of the liquid on the solid.

FIG. 8.

Instantaneous snapshots of the MDPD simulations for fluid-fluid interfaces with different inter-fluid attraction parameter values AWN.

FIG. 8.

Instantaneous snapshots of the MDPD simulations for fluid-fluid interfaces with different inter-fluid attraction parameter values AWN.

Close modal
FIG. 9.

Plot of the interfacial tension γWN versus the inter-fluid attraction parameter AWN.

FIG. 9.

Plot of the interfacial tension γWN versus the inter-fluid attraction parameter AWN.

Close modal

In this test case, a static fluid residing in a nanometer-scale slit pore was simulated using MDPD. The objective is to demonstrate the ability of MDPD to model the static contact angle between a single fluid and bounding solid matter. This is an important requirement for characterizing the wetting properties of fluids in organic-rich, nanoporous shales. MDPD can model the interaction force between solids and liquids by choosing a proper set of parameters, ASL, BSL, and rd in Eq. (1). The subscript “S” and “L” denote the solid and liquid, respectively. For simplicity, we fixed BSL = 25 and rd = 0.75rC, and selected six ASL values to study the dependence of the wettability of the liquid on the solid matter. To model liquid-liquid interactions, the MDPD parameters used in this simulation are ALL = −50, BLL = 25, and rd = 0.75rC. Notice that in a similar simulation (Pan, 2010), different choices of values for rC and rd in Eq. (1) were used. In the present test problem, the simulation domain extended from −30 to 30 in the x direction, from −5 to 5 in the y direction, and from −2.5 to 2.5 in the z direction. Periodic BCs were prescribed in the x and z directions.

To form the solid matter of the slit pore, 4500 solid particles were placed randomly in the two regions of y = [−5, −4] and y = [4, 5], respectively. These were run to thermodynamic equilibrium in these regions using the original DPD method, and then frozen for the rest of the simulations. The width of the slit pore was 8rC, corresponding to 8.616 nm in physical units. Next, 4000 liquid particles were randomly placed in a region bounded by x = [−13, 13] and z = [−4, 4], and run for 10 000 time steps to thermodynamic equilibrium using the MDPD method. The bottom and top solid particle layers of the domain were dense enough to prevent the liquid particles from penetrating indefinitely.

The instantaneous particle distributions corresponding to ASL ranging from −40 to −5 are shown in Figure 10. The larger is the absolute value of ASL, the more wetting is the liquid in the slit pore. For example, the case of ASL = −40 corresponds to a fully wetting effect with a contact angle of θ = 0°. This is shown in Figure 10(a). As the absolute value of ASL was successively reduced, the liquid gradually turned from wetting to non-wetting; see Figures 10(b)10(e). Finally, when ASL = −5, a fully non-wetting effect with a contact angle of θ = 180° is observed in Figure 10(f). Overall, the simulation results presented in this test problem have demonstrated that the MDPD method can be used to conveniently model the contact between solid matter and liquids.

FIG. 10.

Instantaneous particle distribution of a single liquid bounded by solid matter in a nanometer-scale slit pore, simulated by MDPD with different attractive force parameter values ASL.

FIG. 10.

Instantaneous particle distribution of a single liquid bounded by solid matter in a nanometer-scale slit pore, simulated by MDPD with different attractive force parameter values ASL.

Close modal

In this section, we assess the developed capability for problems of practical interest. A workflow that integrates high-resolution digital rock imaging to the MDPD model was utilized for simulating a two-component, two-phase fluid flow in heterogeneous, organic-rich, nanoporous shale. The choice for a region of interest in a shale rock sample and the corresponding fluid flow simulation results are discussed as follows.

Figure 11 displays FIB-SEM reconstructed images of a (5000 nm)3 domain obtained in a Woodford shale sample, which shows three of the main compositions: organic matter, inorganic matter, and pores enclosed in organic matter. As seen in Figure 11(c), most of the pores enclosed in organic matter are relatively small, isolated, and not connected to others. This indicates that no fluid flow is possible between those pores. A fluid flow is only possible in those inter-connected larger pores, which are therefore the candidates of region of interest for simulations. However, a domain like that in Figure 11 is prohibitive for MDPD simulations in terms of storage and computing time. Instead, a smaller subdomain is considered to be representative if it contains sufficient details of interconnected pores with both large and small pore throats enclosed in the organic matter. Notice that it is not trivial to locate such a region in Woodford shale samples. In many individual samples, we found no satisfying region at all. The sample shown in Figure 11 was among the many that we inspected, and was the one we eventually chose.

FIG. 11.

A set of three-phase images of a 5000 × 5000 × 5000 nm3 Woodford shale sample reconstructed by the FIB-SEM process. The inorganic matter is colored green, the organic matter is colored black, and the pores enclosed in the organic matter is colored grey.

FIG. 11.

A set of three-phase images of a 5000 × 5000 × 5000 nm3 Woodford shale sample reconstructed by the FIB-SEM process. The inorganic matter is colored green, the organic matter is colored black, and the pores enclosed in the organic matter is colored grey.

Close modal

We identified a 440 nm × 950 nm × 1600 nm subdomain as a region of interest from the domain in Figure 11, and then used it to prepare the simulations. This subdomain is shown in Figure 12(a), where surface geometries of an organic-matter-hosted pore network are represented by a finite element mesh generated based on the 3D voxel data. The largest pore was denoted as “primary pore” in the rest of the text. Rich details of interconnected large and small pore throats can be seen in the primary pore. It is perhaps only part of a larger pore network, but geometric details of this pore are complex enough for demonstrating fluid flows in nanometer-sized pores. In addition, there is also a large isolated pore in this region, through which we assumed no fluid flow. Figure 12(b) shows 590 608 MDPD solid particles that are used to enclose pore space in the region of interest. The particle coordinates were determined using the voxel map and converted into LAMMPS data format using the procedures described earlier. As one can see, the geometric details of pores are different at multiple locations in this region between Figures 12(a) and 12(b). For example, one can see much more small isolated pores in Figure 12(b), but not in Figure 12(a). This is because the finite elements may not preserve all of the geometric information from the voxel map, depending on the minimum element size specified during the mesh generation. This also indicates the possibility that some small pore throats connecting bigger pores could have been ignored in a finite element mesh, resulting in a false representation of the original geometry of pore networks. However, in the MDPD model, no geometric information was lost by conversion from the voxel map to MDPD particles.

FIG. 12.

Plot of a 440 × 950 × 1600 nm3 block containing a complex pore network enclosed in a Woodford shale sample: (a) the pore surfaces (colored red) represented by finite element mesh, (b) the pore surfaces (colored red) represented by MDPD solid particles.

FIG. 12.

Plot of a 440 × 950 × 1600 nm3 block containing a complex pore network enclosed in a Woodford shale sample: (a) the pore surfaces (colored red) represented by finite element mesh, (b) the pore surfaces (colored red) represented by MDPD solid particles.

Close modal

The values of the attractive interaction parameter, A, for fluid-fluid and fluid-solid interactions are listed in Table I. The Type-1 fluid was defined as a wetting liquid that initially resided in the primary pore, and the Type-2 fluid was defined as a non-wetting liquid to be injected in the primary pore in order to extract the Type-1 fluid. The rest of the MDPD parameters used in this simulation problem are Bij = 25 and rd=0.75rC for all particle-particle interactions between the same and different types. As shown in Figure 12(a), the flow inlet boundary was defined at the open pore throats in the lower left corner of the primary pore, and the flow outlet boundary was defined at the open pore throats in the top right corner. In order to get an initial condition such that the primary pore was fully saturated with the Type-1 fluid, a total of 250 000 Type-1 fluid particles were injected into the empty primary pore with a reflection BC imposed at the flow inlet and flow outlet boundaries to prevent fluid particles from exiting the pore. Additionally 5000 time steps were then run to equilibrate the system. This resulted in a highly pressurized pore space with an averaged fluid pressure of 2.05 in the reduced pressure unit. The coordinates of Type-1 fluid and solid particles were then saved and used as initial conditions for the simulations of Type-2 fluid injection that follow. The injection of Type-2 fluid particles into the primary pore was designed to last for a total time of 30 000 in the reduced time unit, which is roughly equivalent to 20 μs. Simulations were conducted with three Type-2 particle injection rates:

  • Test 1, one particle every five time steps;

  • Test 2, one particle every two time steps;

  • Test 3, one particle every time step.

TABLE I.

Fluid flow simulation in a Woodward shale organic-matter-hosted pore network: specification of the attractive interaction parameters, A.

ASolidType-1 fluidType-2 fluid
Solid  −40 −10 
Type-1 fluid −40 −40 −20 
Type-2 fluid −10 −20 −40 
ASolidType-1 fluidType-2 fluid
Solid  −40 −10 
Type-1 fluid −40 −40 −20 
Type-2 fluid −10 −20 −40 

Upon an approximate conversion of units to the real physical units, these injection rates correspond to ∼2 × 10−12 g/s, ∼5 × 10−12 g/s, and ∼10 × 10−12 g/s, respectively. In this problem, an accurate calibration to match specific types of fluids like in Section III A is out of our scope. To ensure that fluid particles did not arbitrarily exit the pore, the reflection BC was maintained at the flow inlet and flow outlet boundaries. Meanwhile, an ejection BC was used to randomly delete any fluid particle that was very near the outlet boundary at a removal rate equal to the corresponding injection rate. This treatment guaranteed a constant total mass of fluids in the primary pore, and also a roughly constant average fluid pressure during the injection. Notice that this resembled a possible scenario of producing from nanoporous, kerogen-rich shale at reasonable flow rates.

A total of 240 CPU processors (2.66 GHz for each one) were used in parallel computing for each injection test, and the total wall time of each test was about 270 h. Unavoidably, the parallel scalability of the simulations was largely impaired by operations like particle injection at the inlet and particle removal at the outlet, as well as the highly irregular spatial distribution of the geometry. Regular load rebalancing algorithms by re-partitioning subdomains according to the updated particle distribution did not help much in this case. This MDPD system appears containing not a huge number of particles as compared to large-scale MD simulations. If we take the Type-II kerogen (C1512H1764N36O144S18) to be the main composition in the organic matter as an example, where six of such kerogen structures occupy a volume of (3.25 nm)3 (Pawar and Huang, 2016), then an equivalent MD system may contain 3 543 648 kerogen structures and 750 000 fluid molecules. However, the possible physical time evolution by such an MD simulation is severely constrained by the allowable time step size that is at the scale of femtosecond (10−15 s), and by the tremendous MD potential calculations in each step. Therefore, so far it is still practically impossible for MD simulations to model an equivalent system of such size.

Figures 13–15 show a series of instantaneous fluid particle distribution simulated in the three injection tests. Type-1 and Type-2 fluid particles were colored yellow and blue, respectively, and the size of Type-2 fluid particles was set larger to highlight the flow path of the injected fluid. Four scenes are displayed in each figure, with each scene representing a typical moment when the same amount of Type-2 fluid particles had been injected between the three tests. For example, a moment when 1000 Type-2 fluid particles had been injected is depicted in Figures 13(a), 14(a), and 15(a) for each test, respectively. Distribution of the Type-2 fluid particles looks similar in those scenes. Likewise, the scenes in (b)–(d) of each figure look like their counterparts in the other two figures. The equivalence in those flow paths of the injected fluid indicates that the evolution of preferential flow path in the primary pore is relatively independent of the injection rate.

FIG. 13.

MDPD fluid flow simulation in a Woodward shale organic-matter-hosted pore network: instantaneous fluid particle distributions during Test 1.

FIG. 13.

MDPD fluid flow simulation in a Woodward shale organic-matter-hosted pore network: instantaneous fluid particle distributions during Test 1.

Close modal
FIG. 14.

MDPD fluid flow simulation in a Woodward shale organic-matter-hosted pore network: instantaneous fluid particle distributions during Test 2.

FIG. 14.

MDPD fluid flow simulation in a Woodward shale organic-matter-hosted pore network: instantaneous fluid particle distributions during Test 2.

Close modal
FIG. 15.

MDPD fluid flow simulation in a Woodward shale organic-matter-hosted pore network: instantaneous fluid particle distributions during Test 3.

FIG. 15.

MDPD fluid flow simulation in a Woodward shale organic-matter-hosted pore network: instantaneous fluid particle distributions during Test 3.

Close modal

Figure 16 displays the percentage of the Type-1 fluid extracted from the primary pore (ratio of recovery) versus the elapsed injection time. At the beginning, the ratio of recovery for the Type-1 fluid grew linearly with respect to time before the preferential flow path of Type-2 fluid was not fully formed. Once the preferential flow path of Type-2 fluid was developed with the ratio of recovery of the Type-1 fluid at about 60%, the recovery of the Type-1 fluid went to stagnation quickly regardless of the injection rate. In addition, Figure 17 shows the extraction speed of the Type-1 fluid (rate of recovery) versus the elapsed injection time. Notice that the rate of recovery for the Type-1 fluid dropped dramatically soon after the injection of the Type-2 fluid was started, despite the injection rates. The simulation results indicate that it may not be very challenging to reach a relatively low ratio of recovery for the fluid stored in the pores (60% or less). However, efficiency of the recovery is extremely low if 80% or more is pursued. This indicates that smarter geo-techniques have to be developed for further improvement of production. Moreover, the ratio of recovery estimated from the simulations is noticeably higher that those reported in oil and gas industry, e.g., 5% for oil, and about 20% for gas (Kuuskraa et al., 2013). One also has to distinguish the difference that those numbers reported in literature are based on the reservoir-scale production, rather than a tiny representative pore network used in our simulations. Thus from the simulations as well as the data in literature, we could infer that the majority of oil or gas stored in the pores enclosed by organic matters was not recovered from shale reservoirs with existing geo-engineering processes because most of the large pore networks are likely to be stay disconnected and isolated from others in heterogeneous, nanoporous shales (as evidence is already shown in Figure 11(c)). Therefore, it would be no surprise to see those actual low ratios of recovery in industry, if no effective fracturing is performed to create sufficient connections between those large pore networks. This also indicates that more advanced fracturing techniques may be needed.

FIG. 16.

MDPD fluid flow simulation in a Woodward shale organic-matter-hosted pore network: percentage of extracted Type-1 fluid particles versus elapsed injection time.

FIG. 16.

MDPD fluid flow simulation in a Woodward shale organic-matter-hosted pore network: percentage of extracted Type-1 fluid particles versus elapsed injection time.

Close modal
FIG. 17.

MDPD fluid flow simulation in a Woodward shale organic-matter-hosted pore network: extraction rate of Type-1 fluid particles versus elapsed injection time.

FIG. 17.

MDPD fluid flow simulation in a Woodward shale organic-matter-hosted pore network: extraction rate of Type-1 fluid particles versus elapsed injection time.

Close modal

In conclusion, this simulation test problem has completely demonstrated the advantages and robustness of a coherent integration from a high-resolution digital rock imaging process to the MDPD flow model, which enables pore-scale, multi-component, multi-phase fluid flow simulations in real fine-grained, nanoporous shale samples.

A many-body dissipative particle dynamics (MDPD) model has been applied for simulating multi-phase fluid flow in nanopore networks enclosed in fine-grained, nanoporous shales. A distinguished part of this work that sets it apart from the previous ones (Liu et al., 2007a; Liu et al., 2007b and Chen et al., 2011) is the integration of high-resolution FIB-SEM digital imaging technique, which provides 3D voxel data that contain invaluable geometrical and compositional information of shale samples. Although applications of FIB-SEM to studying shales are not new, this is the first time that FIB-SEM has been seamlessly integrated to a Lagrangian particle transport model like MDPD for modeling fluid flows in multi-length-scale nanoporous shale samples, as compared with earlier DPD fluid flow simulations where only manufactured (or simulated) pore networks were used. Since the relevant spatial and temporal scales are too big for molecular dynamics, and too small for computational fluid dynamics with known constitutive models, this integrated research method offers a robust approach to bridging gaps between the molecular- and continuum-scales. Additional remarks worth of noting are as follows.

First, it is not possible for MDPD to properly simulate a fluid dynamics system that contains dissimilar fluids (Ghoufi and Malfreyt, 2011), e.g., water and supercritical CO2, and methane and supercritical CO2. However, those are important scenarios in fine-grained, nanoporous shales. New conservative force formulation must be explored for the MDPD to treat dissimilar fluids. Second, it is computationally prohibitive to use MDPD to simulate a shale sample when the characteristic length of the sample extends from tens of microns to several millimeters. Computer memory required for the digital rock structure reconstruction and data storage could exceed tens of terabytes for a large shale specimen (Goral and Miskovic, 2015). Despite the known limitations, the combination of FIB-SEM and MDPD model offers an unparalleled approach for investigating flow physics in nanopore networks enclosed in organic-rich, heterogeneous, nanoporous shales.

This work is supported by the Idaho National Laboratory (INL) Laboratory Directed Research and Development (LDRD) Program under DOE Idaho Operations Office Contract No. DE-AC07-05ID14517. The authors wish to acknowledge helpful reviews and comments from Ian Walton and John McLennan from Energy and Geoscience Institute. We also thank Jeff Gelb and Jack Kasahara from Carl Zeiss Microscopy for supporting much of this work. Additionally, the authors wish to acknowledge the help and support of Earl Mattson from INL. We are also grateful to U.S. Geological Survey (USGS) Core Research Center for providing access to core materials.

A digital-rock-physics (DRP) technique recently discussed in Goral and Miskovic (2015) is employed to enable a representative and uniquely informative perspective on the structures and chemical compositions of shales. The procedures are presented as follows.

  • Nondestructive micro and nano X-ray microscopy is applied in order to identify a region of interest (ROI) within a shale sample. The selected ROI is then imaged with destructive FIB-SEM serial sectioning.

  • The FIB-SEM data are processed, segmented, and reconstructed into digitized information of organic matter, nonorganic (mineral) matter, and pore space in either of the two phases and between.

  • The ROI of the organic-matter-hosted pore system within the shale rock matrix is identified and separated. The resulting pore network geometry is then gridded with a tetrahedron volume mesh to enable further alternative modeling and simulation studies.

An example of workflow from DRP to numerical modeling is displayed in Figure 18. Notice that image post-processing/segmentation is not a straightforward process, which highly relies on image interpretation. Improper image processing and segmentation can lead to, for example, misidentification of organic matter as pore space, resulting in overestimation of porosity and permeability.

FIG. 18.

A workflow from shale matrix reconstruction and pore network separation to finite element meshes (nonorganic matter—green, organic matter—black, orange—pore network).

FIG. 18.

A workflow from shale matrix reconstruction and pore network separation to finite element meshes (nonorganic matter—green, organic matter—black, orange—pore network).

Close modal
1.
Bear
,
J.
,
Dynamics of Fluids in Porous Media
(
Dover Publications, Inc.
,
New York, USA
,
1972
).
2.
Bertoncello
,
A.
and
Honarpour
,
M. M.
, “
Standards for characterization of rock properties in unconventional reservoirs: Fluid flow mechanism, quality control, and uncertainties
,” in
SPE Annual Technical Conference and Exhibition
(
Society of Petroleum Engineers
,
2013
).
3.
Chen
,
C.
,
Lu
,
K.
,
Zhuang
,
L.
,
Li
,
X.
,
Dong
,
J.
, and
Lu
,
J.
, “
Effective fluid front of the moving meniscus in capillary
,”
Langmuir
29
(
10
),
3269
3273
(
2013
).
4.
Chen
,
C.
,
Lu
,
K.
,
Li
,
X.
,
Dong
,
J.
,
Lu
,
J.
, and
Zhuang
,
L.
, “
A many-body dissipative particle dynamics study of fluid–fluid spontaneous capillary displacement
,”
RSC Adv.
4
(
13
),
6545
6555
(
2014
).
5.
Chen
,
C.
,
Zhuang
,
L.
,
Li
,
X.
,
Dong
,
J.
, and
Lu
,
J.
, “
A many-body dissipative particle dynamics study of forced water–oil displacement in capillary
,”
Langmuir
28
(
2
),
1330
1336
(
2011
).
6.
Curtis
,
M. E.
,
Ambrose
,
R. J.
, and
Sondergeld
,
C. H.
, “
Structural characterization of gas shales on the micro-and nano-scales
,” in
Canadian Unconventional Resources and International Petroleum Conference
(
Society of Petroleum Engineers
,
2010
).
7.
Curtis
,
M. E.
,
Ambrose
,
R. J.
,
Sondergeld
,
C. H.
, and
Rai
,
C. S.
, “
Transmission and scanning electron microscopy investigation of pore connectivity of gas shales on the nanoscale
,” in
North American Unconventional Gas Conference and Exhibition
(
Society of Petroleum Engineers
,
2011
).
8.
Curtis
,
M. E.
,
Cardott
,
B. J.
,
Sondergeld
,
C. H.
, and
Rai
,
C. S.
, “
Development of organic porosity in the Woodford shale with increasing thermal maturity
,”
Int. J. Coal Geol.
103
,
26
31
(
2012a
).
9.
Curtis
,
M. E.
,
Sondergeld
,
C. H.
,
Ambrose
,
R. J.
, and
Rai
,
C. S.
, “
Microstructural investigation of gas shales in two and three dimensions using nanometer-scale resolution imaging
,”
AAPG Bull.
96
(
4
),
665
677
(
2012b
).
10.
Dewers
,
T. A.
,
Heath
,
J.
,
Ewy
,
R.
, and
Duranti
,
L.
, “
Three-dimensional pore networks and transport properties of a shale gas formation determined from focused ion beam serial imaging
,”
Int. J. Oil, Gas Coal Technol.
5
(
2-3
),
229
248
(
2012
).
11.
Espanñol
,
P.
and
Warren
,
P.
, “
Statistical mechanics of dissipative particle dynamics
,”
Europhys. Lett.
30
(
4
),
191
(
1995
).
12.
Ghoufi
,
A.
,
Morineau
,
D.
,
Lefort
,
R.
, and
Malfreyt
,
P.
, “
Toward a coarse graining/all atoms force field (CG/AA) from a multiscale optimization method: An application to the MCM-41 mesoporous silicates
,”
J. Chem. Theory Comput.
6
(
10
),
3212
3222
(
2010
).
13.
Ghoufi
,
A.
and
Malfreyt
,
P.
, “
Mesoscale modeling of the water liquid-vapor interface: A surface tension calculation
,”
Phys. Rev. E
83
(
5
),
051601
(
2011
).
14.
Goral
,
J.
and
Miskovic
,
I.
, “
A workflow for multi-scale modeling and simulation of transport phenomena in Woodford shale rock matrix
,” in
Proceedings of the Unconventional Resources Technology Conference (URTeC)
(
Society of Petroleum Engineers
,
San Antonio, TX, USA
,
2015
).
15.
Groot
,
R. D
and
Warren
,
P. B
, “
Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation
,”
J. Chem. Phys.
107
(
11
),
4423
(
1997
).
16.
Heldele
,
R.
,
Schulz
,
M.
,
Kauzlaric
,
D.
,
Korvink
,
J. G.
, and
Haußelt
,
J.
, “
Micro powder injection molding: Process characterization and modeling
,”
Microsyst. Technol.
12
(
10-11
),
941
946
(
2006
).
17.
Henrich
,
B.
,
Cupelli
,
C.
,
Moseler
,
M.
, and
Santer
,
M.
, “
An adhesive DPD wall model for dynamic wetting
,”
Europhys. Lett.
80
(
6
),
60004
(
2007
).
18.
Hoogerbrugge
,
P. J.
and
Koelman
,
J. M. V. A.
, “
Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics
,”
Europhys. Lett.
19
(
3
),
155
(
1992
).
19.
Huang
,
H.
,
Meakin
,
P.
, and
Liu
,
M.
, “
Computer simulation of two-phase immiscible fluid motion in unsaturated complex fractures using a volume of fluid method
,”
Water Resour. Res.
41
(
12
),
W12413
, doi:10.1029/2005wr004204 (
2005a
).
20.
Huang
,
H.
,
Meakin
,
P.
,
Liu
,
M.
, and
McCreery
,
G. E.
, “
Modeling of multiphase fluid motion in fracture intersections and fracture networks.
,”
Geophys. Res. Lett.
32
(
19
),
L19402
, doi:10.1029/2005gl023899 (
2005b
).
21.
Kuuskraa
,
V.
,
Stevens
,
S. H.
, and
Moodhe
,
K. D.
, “
Technically recoverable shale oil and shale gas resources: An assessment of 137 shale formations in 41 countries outside the United States
,” U.S. Energy Information Administration (
prepared for the U.S. Department of Energy
),
Washington, DC
,
2013
.
22.
Liu
,
M. B.
,
Liu
,
G. R.
,
Zhou
,
L. W.
, and
Chang
,
J. Z.
, “
Dissipative particle dynamics (DPD): An overview and recent developments
,”
Arch. Comput. Methods Eng.
22
(
4
),
529
556
(
2015
).
23.
Liu
,
M.
,
Meakin
,
P.
, and
Huang
,
H.
, “
Dissipative particle dynamics with attractive and repulsive particle-particle interactions
,”
Phys. Fluids
18
(
1
),
017101
(
2006
).
24.
Liu
,
M.
,
Meakin
,
P.
, and
Huang
,
H.
, “
Dissipative particle dynamics simulation of fluid motion through an unsaturated fracture and fracture junction
,”
J. Comput. Phys.
222
(
1
),
110
130
(
2007a
).
25.
Liu
,
M.
,
Meakin
,
P.
, and
Huang
,
H.
, “
Dissipative particle dynamics simulation of pore-scale multiphase fluid flow
,”
Water Resour. Res.
43
(
4
),
W04411
, doi:10.1029/2006wr004856 (
2007b
).
26.
Marsh
,
C.
, “
Theoretical aspects of dissipative particle dynamics
,” Ph.D. thesis,
University of Oxford
,
UK
,
1998
.
27.
Moeendarbary
,
E.
,
Ng
,
T. Y.
, and
Zangeneh
,
M.
, “
Dissipative particle dynamics: Introduction, methodology and complex fluid applications—A review
,”
Int. J. Appl. Mech.
1
(
04
),
737
763
(
2009
).
28.
Pan
,
C.
,
Hilpert
,
M.
, and
Miller
,
C. T.
, “
Pore-scale modeling of saturated permeabilities in random sphere packings
,”
Phys. Rev. E
64
(
6
),
066702
(
2001
).
29.
Pan
,
C.
,
Hilpert
,
M.
, and
Miller
,
C. T.
, “
Lattice-Boltzmann simulation of two-phase flow in porous media
,”
Water Resour. Res.
40
(
1
),
W01501
, doi:10.1029/2003wr002120 (
2004
).
30.
Pan
,
W.
, “
Single particle DPD: Algorithms and applications
,” Ph.D. thesis,
Brown University
,
USA
,
2010
.
31.
Pawar
,
G.
and
Huang
,
H.
, “
A journey into the subnanometer scale realm of organic matter: An application of molecular dynamics simulation approach
,” in
Proceedings of 2016 AIChE Annual Meeting
,
San Francisco, CA, USA
,
13-16 November 2016
.
32.
Pivkin
,
I. V.
and
Karniadakis
,
G. E.
, “
A new method to impose no-slip boundary conditions in dissipative particle dynamics
,”
J. Comput. Phys.
207
(
1
),
114
128
(
2005
).
33.
Pivkin
,
I. V.
and
Karniadakis
,
G. E.
, “
Controlling density fluctuations in wall-bounded dissipative particle dynamics systems
,”
Phys. Rev. Lett.
96
(
20
),
206001
(
2006
).
34.
Plimpton
,
S.
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
(
1
),
1
19
(
1995
).
35.
Revenga
,
M.
,
Zuniga
,
I.
, and
Espanol
,
P.
, “
Boundary conditions in dissipative particle dynamics
,”
Comput. Phys. Commun.
121
,
309
311
(
1999
).
36.
Revenga
,
M.
,
Zuniga
,
I.
,
Espanol
,
P.
, and
Pagonabarraga
,
I.
, “
Boundary models in DPD
,”
Int. J. Mod. Phys. C
9
(
08
),
1319
1328
(
1998
).
37.
Tartakovsky
,
A. M.
and
Meakin
,
P.
, “
A smoothed particle hydrodynamics model for miscible flow in three-dimensional fractures and the two-dimensional Rayleigh-Taylor instability
,”
J. Comput. Phys.
207
(
2
),
610
624
(
2005
).
38.
Tartakovsky
,
A. M.
and
Meakin
,
P.
, “
Pore scale modeling of immiscible and miscible fluid flows using smoothed particle hydrodynamics
,”
Adv. Water Res.
29
(
10
),
1464
1478
(
2006
).
39.
Tiwari
,
A.
and
Abraham
,
J.
, “
Dissipative-particle-dynamics model for two-phase flows
,”
Phys. Rev. E
74
(
5
),
056701
(
2006
).
40.
Visser
,
D. C.
,
Hoefsloot
,
H. C.
, and
Iedema
,
P. D.
, “
Modeling multi-viscosity systems with dissipative particle dynamics
,”
J. Comput. Phys.
214
(
2
),
491
504
(
2006
).
41.
Warren
,
P. B.
, “
Vapor-liquid coexistence in many-body dissipative particle dynamics
,”
Phys. Rev. E
68
(
6
),
066702
(
2003
).