Three-Dimensional Printing of Synthetic Core Plugs as an Alternative to Natural Core Plugs: Experimental and Numerical Study

Abstract

This paper proposes three-dimensional (3D) additive fabrication of synthetic core plugs for core flooding experiments from spheres and grains of Berea Sandstone using a digital particle packing approach. Samples were generated by systematically combining the main textural parameters of the rock reservoir to design synthetic core plugs Numerical flow simulation was per-formed using the lattice Boltzmann method (LBM) to verify the flow distribution and permeability for comparison with the experimentally measured permeability and to that obtained from correlations in the literature. The digital porosity of the sample was compared to the porosity measured using an HEP-P helium porosimeter. The numerical and experimental results for permeability and porosity differed by a maximum of 18%.

1. Introduction

Many areas of applied science and engineering involve the study of fluid flow in porous media. Fields as diverse as geology, hydrogeology, and petroleum engineering rely heavily on it to solve problems. In the oil field, understanding multiphase flow in porous media is the basis for developing recovery scenarios and oil/gas production strategies. The oil and gas industry has traditionally used core flooding tests to experimentally measure permeability, relative permeability, saturation change, porosity, and fluid–rock interactions in reservoir and outcrop rock. These tests involve forcing fluids at high confinement pressure through a core plug within a pressure vessel (core holder) [1]. The fluids used may include nitrogen, reservoir brine, crude oil, drilling fluids, and drilling mud filtrate. Other fluids specifically designed to improve or enhance oil recovery (EOR) are also used in core flooding tests to evaluate the effect of treatment and characterize fluid mobility for injection test design.

EOR is a technique for recovering crude oil from a reservoir after primary production has been exhausted due to natural reservoir pressure. EOR processes consist of injecting fluids into the porous medium of the reservoir to create a physical or chemical stimulation effect that increases oil recovery. Various fluids, such as gas, steam, foams, polymers, surfactants, liquid CO₂, alkalis, microbial products, or combinations of these fluids, are injected into the reservoir [2]. These fluids are selected according to the desired recovery effect to be induced in the reservoir. The effectiveness of a potential EOR process can be tested experimentally via core flooding tests prior to expensive field implementation.

However, there are several problems associated with core flooding tests. One of the main problems is obtaining intact cuttings from the reservoir rock to obtain reliable experimental results. Reservoir rock is cut into cylinders to produce core plugs for experimental testing. This is difficult because of the destruction of rock samples and contamination from drilling mud particles that occur during the oil well drilling process, which make it difficult to obtain sufficient material in sufficiently good condition to produce reliable core plug samples. Each core has a specific mineralogy and a unique pore network configuration, even for cores from the same well, formation, and orientation. As a result, researchers must use cores of different morphologies and geometries to extend or repeat their experiments, and the inherent uncertainties involved make it difficult to interpret experimental results. In addition, many of the core flooding tests are destructive (e.g., reactive flow). Therefore, it is impossible to repeat and design experiments because rock samples can only be used once for a single test. In addition, most tests require sample decontamination to remove residual oil, accumulated salts, drilling mud, and other contaminants that can alter the original morphology of the rock. Because of the difficulty of obtaining natural rock cuttings, the oil industry has used outcrop rock or sandstone cuttings as substitutes for natural rock cuttings [3]. However, even in sandstones, which have a very consistent morphology, it is impossible to obtain samples with identical pore networks for repeated testing. Efforts have been made for decades to represent the geometry of porous rock reliably. Starting from very simplified arrangements using clusters of spheres or bundles of capillary tubes of the same diameter [4].

In recent years, a branch of study has focused on reproducing the real microstructure of porous rocks to generate digital rock samples, using techniques such as microtomography and spatial microscopy and disciplines such as statistical physics [5,6]. The combination of 3D microtomographic image reconstruction of natural rocks and numerical simulation of pore-level flow complements, and in some cases replaces, traditional laboratory core flood testing [7]. Another approach to the reconstruction technique involves the development of porous media packages with different geometric, morphological, and anisotropic configurations using physics-based algorithms to simulate the sedimentation, compaction, and cementation processes that lead to rock formation [8,9,10]. This approach is very useful in developing of morphologically manipulated porous samples for research purposes, where it is necessary to control specific morphological parameters to tune or validate analytical and numerical models. This approach is enhanced by combining it with additive manufacturing, which offers the possibility of 3D printing digital representations of porous structures obtained from CT images and physics-based algorithms [9].

The purpose of study was to investigate whether 3D additive manufacturing of porous samples could reproduce the morphological characteristics of real rocks and serve as viable samples for core flooding experiments. The design of the digital samples was based on the manipulation of typical morphological parameters of the rock, such as grain size and shape, grain size variation (sorting), and cementation. The morphological manipulation of the porous rock provides the opportunity to produce many digital and 3D-printed specimens on an ad hoc basis for research and academic studies, and for tuning numerical and analytical models of flow in porous media.

We used lattice Boltzmann simulations as a complementary method to calculate the flow and predict the permeability of the digitally designed sample prior to 3D printing. The LBM provides an accurate, high-throughput method for solving fluid flow problems in porous media with complex geometries, such as those generated via digital rock physics. LBM is widely used to model pore-scale flow in porous structures [11,12,13,14,15,16,17,18,19,20]. For a general introduction to the application of lattice Boltzmann theory in porous media, see [11].

This paper describes (i) the digital construction of core plugs for core flooding experiments by packing samples with spheres of different sizes and realistic grains with different shapes, obtained via micro-CT of Berea sandstone; (ii) numerical simulation of flow in a digital core plug by means of the LBM method to determine the permeability and verify the connectivity of the pore network; and experimental determination of the permeability and porosity of 3D-printed core plugs via core flooding experiments and a helium porosimeter (HEP-P), respectively.

2. Methods

2.1. Digital Core Plug Design

The digital construction of the samples was based on the main morphological parameters of the rock. The systematic combination of these parameters made it possible to control the two main petrophysical properties of the porous rock medium: permeability and porosity. These are the main macroscopic properties used to characterize and classify porous media and describe fluid flow. Permeability and porosity were determined numerically and experimentally for the digital and 3D-printed samples, respectively, to permit a comparison of the values obtained. The results were used to evaluate the potential of 3D-printed core plugs for use in core flooding tests. Some advantages of producing 3D-printed replicas of digitally generated samples while controlling the morphological parameters of the core samples are as follows.

-

The 3D-printed replicas can be designed ad hoc according to the experimental requirements by manipulating the morphological parameters to produce specific samples to experimentally calibrate numerical and analytical models of fluid flow in a porous medium.

-

The availability of 3D-printed samples with the same pore network allows the design of experiments (DoE) where multiple input factors are manipulated to determine their effect on a desired outcome (response). The DoE can also include intentional changes in the morphological structure of the sample to allow for multilevel factorial experiments.

-

The 3D replication of the samples allows for the repeatability of experiments. In contrast, the pore network of natural rock is regularly destroyed during experimental tests [6,7,21,22].

-

3D-printed replicas offer a low-cost alternative to expensive original samples [23].

-

It is possible to directly compare numerical simulation results and experimental results for samples with the same digital and 3D-printed pore networks, respectively.

-

Morphological manipulation at the digital level allows the development of samples with predicted future morphologies derived from natural changes in the rock structure of reservoirs.

-

During core flooding experiments, unconsolidated rock samples can be 3D-printed to create replicas that can withstand extreme pressures and temperatures [23].

-

The development of digitals specimens, combined with 3D printing, allows for the resizing of specimens for enhanced oil recovery testing, where fluid mobility characterization requires a minimum characteristic length.

Sample packing is performed using a digital particle packing approach based on the physical processes of sedimentation, compaction, and cementation that rocks undergo during formation [10]. In this approach, the particles are a collection of voxels that move in a grid with six orthogonal and 20 diagonal movements. For the digital construction of samples with the morphology of real rocks, digital grains of natural sandstone from Berea were used as particles for the packing of the samples. Digital grains were extracted from X-ray computed tomography of a cube of Berea sandstone with a resolution of 4.87 µm/voxel [9].

In nature, a porous rock medium of the rock is affected by a process of diagenesis, characterized by the consolidation of material between the grains or particles of the porous medium, caused by a dynamic process of maturation of a sediment as it passes through the rock (lithification). This process is characterized by the formation of mineral sediments in the spaces between the rock grains, which causes, among other things, the cementation of the porous medium.

Cementation causes a reduction in pore space, which reduces the porosity and permeability of the porous medium. A voxel-based geometric cementation process was used to simulate this phenomenon during the construction of the specimens. This process consisted of digitally reconstructing bridges between grains using particles with different geometries to simulate the geometric configuration of diagenesis in real rocks [24]. Digital voxel-based geometric cementation allows control of the porosity of the samples by adding or reducing bridges between the interstitial spaces of the grains. The steps for developing core plug samples are as follows:

-

The center of mass and voxel coordinates for each particle were stored in a database. This information was used to create a set of connected 3D voxels for each particle. Particles were extracted from this database to form the packed samples. The geometries ranged from spheres to Berea sandstone grains. The texture of the grains could be modified by changing the voxelization to a convex hull (rounded grains) [9].

-

The sizes and shapes of the particles used to build the sample were selected according to the morphological characteristics desired for the sample. We started with the digital development of typical porous media samples, such as those formed by homogeneous spheres, to represent the porous medium in a simplified way and establish a reference according to Darcy’s law. The samples evolve in their digital construction through the systematic combination of morphological parameters until samples with realistic rock morphology are obtained.

-

Sample packing was performed using the DigiPac digital particle packing tool [10], which uses stochastic and deterministic approaches to particle packing. Digital voxel-based geometric cementation was then added using the digital tool. During this process, the porosity of the sample was determined from the voxelization of the sample.

-

Prior to 3D printing, the permeability of the sample was determined by numerical simulation of fluid flow using LBM.

-

Once the morphological characteristics of the sample and the required petrophysical properties of permeability and porosity were determined, 3D printing was performed. The samples can be rescaled depending on the resolution of the printer and the needs of the experiment. Figure 1 shows the most complex digital packages for the study cases presented in this work, corresponding to study cases 4, 6, 7, and 8 described in Section 3.3. Digital packaging e (Figure 1) was used to calibrate the parameters of the sedimentation models used to develop the digital packages.

2.2. Lattice Boltzmann Method

The lattice Boltzmann method is used to calculate the flow and predict the permeability of a digital design before 3D printing. The fluid is represented by a set of fictitious particles moving independently (flow) at different velocities (impulse discretization) in each grid domain (space discretization) along time steps (time discretization) [13]. LBM consists of two main steps, propagation, and collision of particles. The propagation step represents the movement of particles by which their densities are shifted from one node to another according to the values of the distribution function for each of the possible directions. The collision step represents the interaction between the particles, which change their directional velocities when they reach the next node. The lattice Boltzmann equation (Equation (1)) describes the evolution of a discretized particle distribution function f_i (x, $\overrightarrow{{\mathit{e}}_{\mathit{i}}}$, t), which represents the probability of finding a particle moving and colliding in a lattice, with a given discrete velocity $\overrightarrow{{\mathit{e}}_{\mathit{i}}}$ at time t at a given location in the domain x.

$${\mathit{f}}_{\mathit{i}}\left(\overrightarrow{\mathit{x}}+\overrightarrow{{\mathit{e}}_{\mathit{i}}}\mathsf{\Delta }\mathit{t},\mathit{t}+\mathsf{\Delta }\mathit{t}\right){-\mathit{f}}_{\mathit{i}}\left(\overrightarrow{\mathit{x}},\mathit{t}\right)=-\frac{\mathsf{\Delta }\mathit{t}}{\mathit{\tau }}\left[{\mathit{f}}_{\mathit{i}}\left(\overrightarrow{\mathit{x}},\mathit{t}\right)-{\mathit{f}}_{\mathit{i}}^{\mathit{e}\mathit{q}}\left(\overrightarrow{\mathit{x},}\mathit{t}\right)\right]$$

(1)

where f_i ($\overrightarrow{\mathit{x}}$ t) is the ith direction (i = 0, 1, 2… 18 in our case) density distribution function at the lattice site x and time t with a discrete velocity $\overrightarrow{{\mathit{e}}_{\mathit{i}}}$, ${\mathit{f}}_{\mathit{i}}^{\mathit{e}\mathit{q}}$ ($\overrightarrow{\mathit{x}}$ t) is the equilibrium distribution function, τ is the dimensionless relaxation time related to the viscosity, and Δt is the lattice time step, which is in the lattice units Δt = 1.

The LBM algorithm is implemented in two steps: first, particle collision controls relaxation to equilibrium, and second, particle propagation moves the distribution functions to adjacent lattice cells.

$$\mathrm{collision}:{\mathit{f}}_{\mathit{i}}^{\mathbf{\prime }}\left(\overrightarrow{\mathit{x}},\mathit{t}\right)={\mathit{f}}_{\mathit{i}}\left(\overrightarrow{\mathit{x}},\mathit{t}\right)-\frac{1}{\mathit{\tau }}\left[{\mathit{f}}_{\mathit{i}}\left(\overrightarrow{\mathit{x}},\mathit{t}\right)-{\mathit{f}}_{\mathit{i}}^{\mathit{eq}}\left(\overrightarrow{\mathit{x}},\mathit{t}\right)\right]+{\mathit{F}}_{\mathit{i}}\left(\overrightarrow{\mathit{x}},\mathit{t}\right){\mathbf{\Delta }}_{\mathit{t}}$$

(2)

$$\mathrm{propagation}:{\mathit{f}}_{\mathit{i}}\left(\overrightarrow{\mathit{x}}+\overrightarrow{{\mathit{e}}_{\mathit{i}}}\mathbf{\Delta }\mathit{t},\mathit{t}+\mathbf{\Delta }\mathit{t}\right)={\mathit{f}}_{\mathit{i}}^{\prime }\left(\overrightarrow{\mathit{x}},\mathit{t}\right)$$

(3)

constant external force F_i, defined as a constant pressure gradient to drive the fluid flow, is added to the RHS of Equation (2), which is defined as ${\mathit{F}}_{\mathit{i}}=-\mathbf{\nabla }\mathit{p}/{\mathit{\rho }}_{\mathit{o}}$ and the single relaxation time (SRT) τ is related to the kinematic lattice viscosity, ν.

$$\mathbf{\nu }=\left(\mathsf{\tau }-\mathbf{0.5}\right){\mathit{c}}_{\mathit{s}}^{\mathbf{2}}$$

(4)

where c_s is the nondimensional sound speed (c_s = 1/√3). In this study, we consider the Bhatnagar–Gross–Krook (BGK) approximation for the collision implementation (right side of Equation (1) and defined in Equation (2)), which models a collision as a linear relaxation of the distribution function toward equilibrium in the Lattice Boltzmann equation. The local equilibrium distribution function ${\mathit{f}}_{\mathit{i}}^{\mathit{e}\mathit{q}}$ is given as follows:

$${\mathit{f}}_{\mathit{i}}^{\mathit{e}\mathit{q}}={\mathit{w}}_{\mathit{i}}\mathit{\rho }\left(\mathbf{1}+\frac{\overrightarrow{{\mathit{e}}_{\mathit{i}}}·\mathit{u}}{{\mathit{c}}_{\mathit{s}}^{\mathbf{2}}}+\frac{{\left(\overrightarrow{{\mathit{e}}_{\mathit{i}}}·\mathit{u}\right)}^{\mathbf{2}}}{\mathbf{2}{\mathit{c}}_{\mathit{s}}^{\mathbf{4}}}-\frac{{\left|\mathit{u}\right|}^{\mathbf{2}}}{\mathbf{2}{\mathit{c}}_{\mathit{s}}^{\mathbf{2}}}\right)$$

(5)

where ${\mathit{w}}_{\mathit{i}}$ is the weight associated with the velocity, $\overrightarrow{{\mathit{e}}_{\mathit{i}}}$; and $\mathit{\rho }$ and u are the density and macroscopic velocity, respectively, that must satisfy the requirement for a low Mach number, i.e., u/c_s = Ma << 1. These macroscopic quantities (density and velocity) can be computed in terms of the moments of the velocity distribution functions as follows:

$$\mathit{\rho }={\sum }_{\mathit{i}}{\mathit{f}}_{\mathit{i}}$$

(6)

and

$$\mathit{u}=\sum _{\mathit{i}}\overrightarrow{{\mathit{e}}_{\mathit{i}}}{\mathit{f}}_{\mathit{i}}+\mathit{\tau }\mathit{\rho }{\mathit{F}}_{\mathit{i}}^{\prime }$$

(7)

where, ${\mathit{F}}_{\mathit{i}}^{\prime }={\mathit{F}}_{\mathit{i}}/{\mathit{\rho }}_{\mathbf{0}}$, the body force is restricted to be along the X axis. D3Q19 model (D is for dimension and Q for the number of discrete directions) is used in the present work, which has the velocity vectors ${\overrightarrow{\mathit{e}}}_{\mathit{i}}=\left[{\overrightarrow{\mathit{e}}}_{\mathit{i}\mathit{x}},{\overrightarrow{\mathit{e}}}_{\mathit{i}\mathit{y}},{\overrightarrow{\mathit{e}}}_{\mathit{i}\mathit{z}}\right]$, with ${\mathit{w}}_{\mathit{i}}$, a lattice set of weighting coefficients, which are: ${\mathit{w}}_{\mathbf{0}}$ = 1/3; ${\mathit{w}}_{\mathbf{1}–\mathbf{6}}$ = 1/18 and ${\mathit{w}}_{\mathbf{7}–\mathbf{18}}$ = 1/36.

Packed samples are incorporated directly into the LBM code. The particles (grains or spheres) are fixed during the simulation procedure. The boundary conditions applied were bounce back (for the internal solid no-slip boundary) and periodical (for the external boundary). Numerical simulation of fluid flow in digital rocks using the LBM is described in detail in the references [11,12,13,14,15,16,17,18,19,20].

The input parameter (body force) was varied by trial and error until surface velocities below 10⁻⁴ were obtained to calibrate the simulation. The numerical simulation of the flow for each case study presented in this work is performed in a computational domain larger than the Representative Element Volume (REV) of the specimen, which is the smallest computational domain, to ensure the reliability of the numerical results [15,16]. One of the advantages of numerical simulation is the visualization of the flow distribution in the pore network of the specimen, which allows the identification of errors in the digital design through pore connectivity and flow channeling. In addition, the LBM method can be directly integrated with X-ray micro-tomography or digitally generated microstructures.

2.3. Numerical Permeability

The permeability of the sample can be determined using the Darcy equation, expressed in lattice units (LU) and time step (TS), as shown in Equation (8) [9].

$$\mathit{k}=\frac{\mathit{U}\mathbf{\nu }}{\mathit{f}}$$

(8)

where k is the permeability of the core plug (LU²); ν = (2τ − 1)/6 is the kinematic viscosity (LU² TS⁻¹); f is a body force (LU TS⁻²); and U is the calculated velocity averaged over the weighted pore area obtained from the numerical simulation over the entire flow domain, including the solids (LU TS⁻¹). The input parameters for the LBM numerical simulation are the relaxation parameter (τ) and the body force (${\mathit{F}}_{\mathit{i}}$), which must be greater than 0.5 and less than 0.015, respectively, for numerical stability. Various external forces, including gravity, Lorentz, and Coriolis, can be applied to the fluid. The external pressure gradient applied in the x-direction can also be considered an external force field. In this study, we used the body force ${\mathit{F}}_{\mathit{i}}$, as $\left(\mathit{d}\mathit{p}/\mathit{d}\mathit{x}\right)/\mathit{\rho }$ (where $\mathit{\rho }$ is the mass density of the fluid). Figure 2 shows the velocity profiles obtained for case studies 5, 6, and 7 in LU.

To interpret the simulation results independent of the physical units, it is common to write the permeability equation in dimensionless form by scaling with the characteristic length. Therefore, the permeability $\mathit{k}$ obtained from Equation (8) is normalized using the squared mean diameter D², as shown in Equation (9).

$${\mathit{k}}^{\mathbf{\prime }}=\frac{\mathit{k}}{{\mathit{D}}^{\mathbf{2}}}$$

(9)

where the mean particle (grain or sphere) diameter D serves as the characteristic length in LU or physical units, which is calculated as the harmonic average of the individual particle sizes [8], for each grain packing configuration defined in the case studies in Section 3.3.

$$\mathit{D}=\frac{{\sum }_{\mathit{i}}^{\mathit{n}}{\mathit{f}}_{\mathit{i}}{\mathit{V}}_{\mathit{i}}}{{\sum }_{\mathit{i}}^{\mathit{n}}\frac{{\mathit{f}}_{\mathit{i}}{\mathit{V}}_{\mathit{i}}}{{\mathit{D}}_{\mathit{i}}}}$$

(10)

where f_i is the number fraction of the ith particle with diameter D_i and volume V_i.

The volume-equivalent diameter, defined as the diameter of the sphere with the same volume as the grain (Equation (10)), was used to determine the individual diameter of each grain (D_i = ∛((6V_i)/π)), with LU or physical units. For each case study, the numerical permeability was compared with the dimensionless permeability obtained from correlations reported in the literature. The correlations used have the structure of Equation (11), where the coefficients a and b depend on the characteristic morphological parameters of the granular medium, such as the shape and distribution of the grains and the cementation in the pore spaces [8].

$$\frac{\mathit{k}}{{\mathit{D}}^{\mathbf{2}}}=\mathit{a}{\mathbf{\varnothing }}^{\mathit{b}}$$

(11)

Torskaya et al. [8] used sedimentation algorithms to construct grain packings of various geometries. They constructed a base case consisting of sandstone grain packing with a realistic shape and distribution. The packing grains were replaced by spheres, spherical grains, and ellipsoidal grains. They used fluid flow modeling via finite difference approximation and tomographic imaging to determine the permeability of the packings and obtain the values of the parameters a and b in Equation (11). Torskaya’s work was used as a reference to compare the dimensionless permeability results obtained in this study.

2.4. Three-Dimensional Additive Manufacturing of Core Plug

MultiJet Printing (MJP) technology, which uses an inkjet printing process with piezoelectric print heads to deposit photo-curable plastic resin and molten wax layer by layer, was used in this study to additively manufacture artificial core plugs. The MJP 3600 Max (3D Systems, Rock Hill, SC, USA) industrial printer was used with a resolution of 16 µm/layer. This wax-assisted technology keeps the pores free of resin material to improve the representation of voids as rock pore networks. The core plug was digitally embedded in a housing (core holder) that confines the core plug particle packing, preventing particle deformation and disintegration during core flooding tests. The core holder is a pressurized chamber that directs fluids through the sample. The array core plug and core holder were digitally assembled and saved as an STL file for 3D printing. This array can be scaled to other dimensions while maintaining the morphological structure of the porous sample. This is useful for adapting to the resolution of the printer used or to the needs of the experimental design.

The procedure for printing the core plug consists of four steps. The first step is to design a target 3DP pattern using particle stacking and densification tools such as Digipac [10], AutoCAD^® V 2023, or proprietary digital tools. The second step is to convert the digital structure of the 3D rock model into an STL file that most 3D printers can recognize. In this step, the core plug and core housing assembly can be rescaled to the requirements of the experimental design (DoE) or the resolution of the printer. The third step is to print a target solid rock model. The fourth step is to remove support material from the physical rock models (core plugs). Figure 3 shows the selection of grains for the core plug with different levels of voxelization, from the digitized grains of the Berea Sandstone microtomography to smoothed and rounded grains (convex hull). The grains are packed to form the digital core. The figure al-so shows a natural Berea sandstone core for comparison to the digitally constructed core.

One of the greatest challenges in 3D printing core plugs using the FDM printing technique is the removal of the support material (wax) from the pore cavities. First, the sample was pre-treated with an ultrasonic cleaning system at a temperature of 90.0 °C (363.15 K) to remove as much wax as possible. The core block was then placed in the core flooding system to circulate a commercial cleaning solution (RESINAWAY^®; Monocure 3D, Chicago, IL, USA) at temperatures ranging from 70.0 °C (343.15 K) to 90.0 °C (363.15 K) using a computerized positive displacement pump (PDP). The core flooding system is thermoregulated using an air bath oven (maximum operating temperature of 363.15 K). The cleaning times are long and depend on sample porosity, temperature, and materials used. The characterization of the removal of support material is beyond the scope of this work. Each sample was allowed to drain for days prior to experimental testing to evaporate any remaining cleaning solution prior to the core flooding test. Figure 4 shows a 3D-printed core plug, Figure 4a shows only the core plug, and Figure 4c show the core plug is embedded in a core holder.

2.5. Core Flooding Experimental Setup

The experimental setup of the core flooding system used in this study for the determination of the experimental permeability of the manufactured core plugs is shown in Figure 5. The setup consisted of the following components: (i) a BFSP-1000-15 Vinci syringe computerized positive displacement pump (PDP) (Vinci Technologies, Nanterre, France), (ii) a stainless-steel high-pressure cylinder (500 mL capacity) equipped with a floating piston, (iii) a 3D-printed core block, (iv) a differential pressure transmitter (DPT- EJX115A; Yokogawa, Mexico City, Mexico), and (v) a back-pressure regulator (BPR). The DPT was connected to the inlet and outlet ports of the core block, which was calibrated against a deadweight balance with an uncertainty of ±0.01% at full scale. A graphical user interface (GUI) and data acquisition system was used to control and monitor the pressure, volume, and temperature in the core flooding system.

The key element in this approach is the 3D-printed core block, which consists of the core holder and the core plug (porous medium). The advantage of 3D printing the core block is that the core plug particle assembly is contained within the core holder, preventing the deformation or disaggregation of the particles during experimental testing. The core block caps are removable to facilitate the removal of the support material (wax) from the porous media sample during the cleaning process (Figure 6). The caps have a built-in plug and distribution plate that can be designed in various configurations. The core block is a “plug and play” element, designed ad hoc for the needs of the experimental design, which can be inserted and replaced directly into the experimental setup and replaced as needed, and the same core block can be printed as many times as necessary to validate and repeat experiments. A constant overburden pressure may be required on the core holder to withstand the confined pressures required during the execution of experimental tests, and various commercially available 3D printing materials can be used to achieve the required pressures and temperatures. Today, 3D printing material technology has advanced rapidly, and the market offers a wide range of available materials, such as ULTEM 9085 for pressures up to 69 MPa and temperatures up to 153.0 °C (426.15 K).

3. Results

3.1. Experimental Permeability

Eight core flood experiments were conducted using the morphological configurations described in Section 3.3. The samples used in the experimental tests were 20 mm (0.788 in) in diameter and up to 200 mm (7.874 in) long. Figure 6 shows the 3D-printed array core plug array and core holder for the core flooding tests.

The core block was installed horizontally in the experimental setup (Figure 5) and the PDP was operated in constant-flow (steady-state) mode for all cases studied. Pressure ports are located at both ends of the core block to measure the differential pressure at the inlet-outlet sample. The data was transferred to the computer and plotted over time. When the differential pressure has stabilized and showed gradual variation on the graph, it was used with the volumetric flow rate data to calculate the average gas permeability (k) according to Darcy’s modified equation (Equation (12)) for compressible gas flow [25,26]. This equation is valid for laminar flow when there is a linear relationship between the pressure gradient and gas volume flow.

$$\frac{{\mathit{P}}_{\mathbf{2}}^{\mathbf{2}}-{\mathit{P}}_{\mathbf{1}}^{\mathbf{2}}}{\mathbf{2}{\mathit{P}}_{\mathbf{1}}\mathit{L}}=\frac{\mathit{\mu }}{\mathit{k}}\frac{\mathit{Q}}{\mathit{A}}$$

(12)

where k is the permeability of the sample, μ is the gas viscosity (Pa-s), L is the length of the core plug (m), Q is the volumetric flow rate (m³/s), and P₁ and P₂ are the pressures at the low- and high-pressure sides of the 3D-printed core holder (Pa), respectively. In cases of turbulent flow (Reynolds numbers greater than 10) and very low sample permeabilities (less than 10 md), inertial effects and gas slip effects (the so-called Klinkenberg effect), respectively, must be considered in Equation (12) [27,28]. The experimental permeability obtained using Equation (12) is divided by the square of the harmonic mean diameter of each sample’s individual particle sizes (Equation (10)) to obtain the dimensionless experimental permeability (Equation (9)). The core flooding experiments were performed over a volume flow range of 600 to 1800 mL/h to obtain Reynolds numbers less than 1, and the permeabilities of the samples were greater than 10 mD, so that the Klinkenberg effects were not included in the Equation (12). The dimensionless experimental permeability obtained for each of the study cases is presented in Section 3.3.

3.2. Experimental Porosity

In natural rocks, the original (primary) porosity is modified via post-depositional processes to produce secondary porosity. These processes result from rock compaction and chemical and biological processes that lead to the formation of mineral deposits or cementation in the pore space. Cementation causes a reduction in pore space, which reduces the porosity and permeability of the rock. To simulate cementation in the constructed samples, the original porosity of the sample was modified by developing voxel-based geometric cementation by inserting particles of different geometries between the grains of the sample (hyperboloids, ellipsoids, and cylinders). This was achieved by building digital bridges in the pore space of the sample using a proprietary voxel-based geometric algorithm developed in FreeCAD software (Version 0.18.3-2019) [29]. Cementation allows porosity to be created in the sample according to the needs of the experimental design. Digital voxel-based geometric cementation allows control of the porosity of the samples by adding or reducing the bridges between the interstitial spaces of the grains.

In this study, a Vinci HEP-P helium porosimeter was used to determine the effective porosity (${\mathbf{\varnothing }}_{\mathit{e}\mathit{f}\mathit{f}})$ of the samples, which is a measure of the fluid storage capacity of a rock and is defined as the fraction of the rock’s total volume corresponding to spaces that can store fluids. The limits of its values for any porous medium range from zero to one. Most porosimeters measure the effective porosity (volume of interconnected pores), which is of interest for estimating oil and gas in place and is mathematically defined as shown in Equation (13)

$${\mathit{\varnothing }}_{\mathit{e}\mathit{f}\mathit{f}}=\frac{{\mathit{V}}_{\mathit{i}\mathit{v}}}{{\mathit{V}}_{\mathit{b}}}$$

(13)

where ${\mathit{\varnothing }}_{\mathit{e}\mathit{f}\mathit{f}}$ is the effective porosity and ${\mathit{V}}_{\mathit{i}\mathit{v}}$ and ${\mathit{V}}_{\mathit{b}}$ are the volume of the interconnected voids and the bulk volume, respectively.

The result of the experimental process is the measurement of the volume of the solids placed in the expansion chamber of the porosimeter. For 3D-printed samples, the total volume (bulk volume) of the solid consists of the volume of the particles, spheres or grains, and the volume of the cementation, plus the volume of the sample holder. The expansion porosimeter is based on Boyle’s ideal gas law, which describes isothermal expansion for gases with ideal behavior. The pressure range typically used is 90 to 100 psi, which allows the gas to be modeled as ideal. The system was calibrated using steel discs to obtain a reference volume. Each measurement was repeated three times and averaged to an accuracy of 0.1%. Helium gas has a low molecular weight, which allows the molecules to penetrate the pore space due to its high diffusivity. In addition, helium does not adhere to the walls of the sample and core holder. The samples used for porosity measurements had a core diameter of 1.5 in and a length of up to 2 in. Figure 7 shows the digital porosity profile obtained via voxelization of the core plug and the 3D-printed sample used to find the porosity of case study 4. Details of the helium porosimeter porosity measurement procedure can be found in References [30,31]. The digital porosity results and the experimental helium porosity results are shown in the next section.

3.3. Study Cases

The first case study corresponds to a sample formed f-homogeneous sphere, which has been widely used to represent porous media and was taken as a reference case. In case 2, the spheres are increased in size to embed each other and achieve a reduction in porosity. Case 3 is a mixture of spheres of different sizes to manipulate the permeability of the sample. The morphological complexity of the subsequent case studies was progressively increased through mixing particles of different sizes and shapes. In addition, smaller particles were incorporated between the pore spaces to control the porosity of the sample, simulating the natural cementation of the rock. For cases 5 and 6, grains obtained from microtomography of Berea sandstone were used as packing particles. In case 5, the core plugs were formed with the same grain (size and shape), while in case 6, the porosity was reduced by using ellipsoidal bridges between the grains of case 5. Finally, for cases 7 and 8, samples were formed with grains of different sizes and shapes, with geometric diagenesis simulated by placing ellipsoidal bridges between grains to generate realistic samples.

The numerical and experimental permeabilities were made dimensionless by dividing them by the square of the harmonic mean diameter of each sample (Equation (9)) for the purpose of comparing the numerical and experimental results. The numerical dimensionless permeability was also compared to the dimensionless permeability obtained via the correlation shown in Equation (11). The results are presented in

Table 1. Comparison of results obtained for dimensionless permeability and porosity of 3D-printed cores.

Case	Sample Morphology	Mean Part. Diameter μm	Porosity			Dimensionless Permeability (k/D2) × 104			Differ. Numer–Analyt. (%)	Differ. Numer–Exptl. (%)
			Digit.	HeP	Differ. (%)	Num.	Analyt.	Exptl.
1	Mono-sized spheres	30	0.41	0.38	7	9.6120	9.0959 †	8.8564	5.36	7.86
2	Mono-sized embedded spheres	12.63	0.277	0.26	6	2.4980	2.2924 †	2.1473	8.23	14.03
3	Mixed spheres	26.58	0.407	0.29	7	9.7700	8.8713 †	8.5691	9.19	12.29
4	Mixed spheres with cementation	31.4	0.331	0.30	10	4.7846	4.3115 †	3.9369	9.88	17.71
5	Mono-sized grains	22.85	0.419	0.40	5	10.674	9.6544 ‡	9.500	9.55	10.99
6	Mono-sized grains with ellipsoids	54.61	0.350	0.33	6	5.6915	5.1265 ‡	4.710	9.92	17.24
7	Multi-sized grains	38.95	0.412	0.40	3	9.8362	9.1038 ‡	8.718	7.44	11.35
8	Multi-sized grains with cementation	48.98	0.300	0.27	10	2.8944	2.9935 ‡	2.378	3.42	17.84

Coefficients a and b in Equation (11) for cases 1–4 ^† are 0.021 and 3.52, respectively. For cases 5–8 ^‡, coefficients a and b are 0.020 and 3.49, respectively.

. The coefficients a and b in Equation (11) for cases 1–4 are 0.021 and 3.52, respectively. These coefficients correspond to packed spheres with the same grain size distribution of sandstone. For cases 5–8, the coefficients a and b are 0.020 and 3.49, respectively, for packed samples with spherical grains. Table 1 also shows the porosity values obtained digitally and experimentally with the Vinci HEP-P helium porosimeter (HeP) (Vinci Technologies, Nanterre, France).

The experimental dimensionless permeability values are lower than the dimensionless permeability values obtained through numerical simulation and those obtained through the correlation (Equation (11)). For the non-cemented samples (cases 1, 3, 5, and 7), the percentage differences between the experimental and numerical dimensionless permeabilities are close to 2%. For the cemented samples (cases 2, 4, 6, and 8), the percentage differences reach values close to 18%. The experimental results for permeability show the same trend as those for porosity, with values lower than the numerical results. Although porosity does not appear in Equation (12), permeability is dependent on porosity. A factor to be considered in reducing the porosity of the samples and, consequently, the permeability is the presence of residues of the support material in the pore space of the sample. In particular, the percentage reductions in permeability and porosity are higher for the cases with cementation between the grains (cases 2, 4, 6, and 8) than for those without cementation (cases 1, 3, 5, and 7) and even higher than for the results obtained from the numerical simulation, in which, of course, support material in the pore space is not considered.

4. Discussion

The core flooding experiments were conducted with the overall objective of determining whether 3D-printed core plugs have potential as a substitute for real rock core plugs. To perform this evaluation, we proposed the digital design and 3D additive manufacturing of eight samples (core plugs), varying the main morphological parameters for each case study, such as the grain size and shape, sorting, and cementation. These parameters dominate the sample porosity and permeability. These petrophysical properties are the most representative of the porous medium. Therefore, we determined the porosity and permeability value of each of the 3D-printed samples by three means—numerical flow simulation, an analytical model reported in the literature, and experimental core flooding tests—to permit direct comparisons between the results obtained and assess whether porosity and permeability are reproducible in 3D-printed cores designed ad hoc for experimental research purposes.

We believe that 3D printing is a good way to produce cores with morphologies that can be adapted to research needs. In EOR processes, in particular, adequate characterization of the mobility of injection fluids in the reservoir is crucial for the design and control of the injection test. Having samples with the same pore network and morphological configuration susceptibility for which the porosity and permeability can be manipulated represents an opportunity to extend the experimental design (DoE). The numerical and experimental results for dimensionless permeability and porosity show differences of 18% and 10%, respectively. We consider these acceptable degrees of difference that can be reduced with technological advances in 3D printing and the development of new materials.

5. Conclusions

The strength of this approach is the digital construction of the sample by manipulating morphological parameters such as the shape and size of the particles (grains or spheres) and the degree of cementation between them to achieve the desired permeability and porosity. Digital rock physics, combined with particle packing algorithms and 3D printing technology, allows the physical construction of porous samples on an ad hoc basis for experimental designs required for purely experimental research as well as for the validation of analytical/numerical formulations of fluid flow in hydrocarbon production systems. In particular, the samples presented in the present work are oriented towards fluid mobility studies in laboratory-scale EOR tests (Figure 6). The disadvantage of the MJP 3D printing technique used in this work is the difficulty in removing the support material, especially for samples with interparticle bridges (cementation cases). However, the development of new 3D printing and removal materials is increasing the reliability of 3D printing as an alternative replacement for reservoir rock samples.

In nature, it is virtually impossible to find porous media with the same morphological characteristics, even in sandstones such as Berea, which has a slightly heterogeneous morphology. Historically, the inability to find porous media with the same morphological characteristics has severely limited experimental work. The generation of synthetic and 3D-printed rocks offers the possibility of having replicas to perform a series of experiments under different process conditions in the same pore network. This fact will change how flow is studied in porous media created ad hoc according to the morphology required by the target research, since it will be possible to compare numerical and experimental results in the same pore network. Having 3D-printed replicas will allow the design of factorial experiments (DoE). In most experiments in porous media related to hydrocarbon recovery, the porous sample is permanently damaged. The approach presented in this work allows the development of samples with more complex morphologies, with fractures and vugular structures, for the digital and physical reproduction of naturally fractured rock reservoirs.