### Recent Submissions

• #### A decoupled scheme to solve the mass and momentum conservation equations of the improved Darcy–Brinkman–Forchheimer framework in matrix acidization

(AIP Advances, AIP Publishing, 2021-12-01) [Article]
Matrix acidization simulation is a challenging task in the study of flows in porous media due to the changing porosity in the procedure. The improved Darcy–Brinkman–Forchheimer framework is one model to do this simulation. In this framework, the mass and momentum conservation equations are discretized to form a pressure–velocity linear system. However, the coefficient matrix of the linear system has a large condition number, and solving the linear system belongs to the saddle point problem. As a result of that, convergence is hard to achieve when solving it with iterative solvers. It is well known that the scale of the linear systems in matrix acidization simulation is large, and therefore, the usage of iterative solvers is required. Thus, a decoupled scheme is proposed in this work to decouple the pressure–velocity linear system into two independent linear systems: one is to solve for pressure, and the other one is to solve for velocity. It is emphasized that both of the linear systems are discretized from the elliptical partial differential equations, which guarantees fast convergence can be achieved by iterative solvers. A numerical experiment is carried out to demonstrate the correctness of the decoupled scheme and its higher computing efficiency. After that, the decoupled scheme is applied in investigating the factors that cannot change the optimal injected velocity and the dissolution pattern in matrix acidization.
• #### Thermodynamics-Informed Neural Network (TINN) for Phase Equilibrium Calculations Considering Capillary Pressure

(Energies, MDPI AG, 2021-11-18) [Article]
The thermodynamic properties of fluid mixtures play a crucial role in designing physically meaningful models and robust algorithms for simulating multi-component multi-phase flow in subsurface, which is needed for many subsurface applications. In this context, the equation-of-state-based flash calculation used to predict the equilibrium properties of each phase for a given fluid mixture going through phase splitting is a crucial component, and often a bottleneck, of multi-phase flow simulations. In this paper, a capillarity-wise Thermodynamics-Informed Neural Network is developed for the first time to propose a fast, accurate and robust approach calculating phase equilibrium properties for unconventional reservoirs. The trained model performs well in both phase stability tests and phase splitting calculations in a large range of reservoir conditions, which enables further multi-component multi-phase flow simulations with a strong thermodynamic basis.
• #### An energy stable linear numerical method for thermodynamically consistent modeling of two-phase incompressible flow in porous media

(Journal of Computational Physics, Elsevier BV, 2021-11-17) [Article]
In this paper, we consider numerical approximation of a thermodynamically consistent model of two-phase flow in porous media, which obeys an intrinsic energy dissipation law. The model under consideration is newly-developed, so there is no energy stable numerical scheme proposed for it at present. This model consists of two nonlinear degenerate parabolic equations and a saturation constraint, but lacking an independent equation for the pressure, so for the purpose of designing efficient numerical scheme, we reformulate the model forms as well as the free energy function, and further prove the corresponding energy dissipation inequality. Based on the alternative reformulations, using the invariant energy quadratization approach and subtle semi-implicit treatments for the pressure and saturation, we for the first time propose a linear and energy stable numerical method for this model. The fully discrete scheme is devised combining the upwind approach for the phase mobilities and the cell-centered finite difference method. The unique solvability of numerical solutions and unconditional energy stability are rigorously proved for both the semi-discrete time scheme and the fully discrete scheme. Moreover, the scheme can guarantee the local mass conservation for both phases. We also show that the upwind mobility approach plays an essential role in preserving energy stability of the fully discrete scheme. Numerical results are presented to demonstrate the performance of the proposed scheme.
• #### Interfacial properties of the alkane+water system in the presence of carbon dioxide and hydrophobic silica

(Fuel, Elsevier BV, 2021-11-10) [Article]
Molecular dynamics simulations were carried out to understand the interfacial properties of the alkane+water system in the presence of CO and hydrophobic silica at temperatures from 323 to 443 K and pressures up to about 200 MPa. The simulation data were compared to predictions from density gradient theory. Our results of the interfacial tension (IFT) of the alkane+water and alkane+CO+water systems were in reasonable agreement with the experimental data. At a given temperature and pressure, the IFT of the alkane+water system almost linearly increases with the number of carbon atoms in the alkane molecule . The IFTs of the alkane+CO+water system are relatively similar to those reported for the corresponding alkane+water system. The addition of CO decreased the IFT of the alkane+water system. For a given , the IFT is approximately equal for linear, branched, and cyclic alkanes in the presence of water and CO. The water contact angle obtained from simulations of the alkane+water+silica system is in the range of about 117–139. This contact angle decreases with pressure, and in general, the higher the temperature, the more pronounced is this pressure effect. Overall, the contact angle is higher for lower and cyclic alkanes, but branching has no noticeable effect on the contact angle. The contact angles of the CO+water+silica system were in reasonable agreement with experimental data. The contact angle increased with increasing pressure and decreasing temperature for this system. The contact angles of the dodecane+CO+water+silica system are relatively similar to those reported for the corresponding dodecane+water+silica system. The addition of CO increased the contact angle of the dodecane+water+silica system.
• #### Molecular Dynamics Modeling of Kaolinite Particle Associations

(The Journal of Physical Chemistry C, American Chemical Society, 2021-10-18) [Article]
We developed a new procedure for calculating finite-size kaolinite particles, their associations with complex surface chemistry, and the natural flexibility of sheets within a particle using a large-scale atomic/molecular massively parallel simulator. For the first time, all possible particle associations previously described in the literature were obtained using an atomic method. The structural configurations obtained were shifted face-face, angular edge-edge, corner-corner, and shifted face-face-face booklet associations. The simulations showed that if the initial angle between two interacting particles is less than 45°, the particles will form layer-to-layer aggregates. If the angle is larger than 60°, the particles will form an angular arrangement. The densities of kaolinite arrangements with dense and loose packings were evaluated as a function of the structure. The densest structures, as expected, were the layered structures, with four and two layers. The density of the shifted face-face packing was about the same density as the two. The face-face-face association showed lower density, and the angular edge-edge association showed a 3 times lower density than the densest, four-layer structure.
• #### Fully implicit two-phase VT-flash compositional flow simulation enhanced by multilayer nonlinear elimination

(Journal of Computational Physics, Elsevier BV, 2021-10) [Article]
VT flash calculation, with the variable specification of volume, temperature and mole numbers, is a newly-rising alternative to the conventional PT flash calculation for phase behavior modeling. Until now, VT flash is primarily used as standalone calculation to solve phase equilibria problems and few works apply it to compositional flow simulation. In this study, an improved two-phase VT-flash compositional flow algorithm is developed with the multilayer nonlinear elimination method. Compared to the preceding works, the robustness and efficiency of the new algorithm is significantly improved as the nonlinear elimination method implicitly removes locally large nonlinearities and thus restore large step length for Newton iterations. To enhance the computational efficiency, an adaptive time stepping control is used to adjust the timestep size. Moreover, unnecessary stability tests are bypassed using a modified shadow region method, which accommodates to substantial composition change under large time steps. A number of numerical examples are presented to demonstrate the robustness and efficiency of the proposed VT-flash compositional flow algorithm with multilayer nonlinear elimination. Even though the convergence issue is not fully resolved, which roots in the nondifferentiable equilibrium pressure at phase boundary, the occurrence of time refinements is significantly reduced with the help of multilayer nonlinear elimination. It is also found that multilayer nonlinear elimination generally increases the number of Newton iterations slightly but enlarges the timestep size significantly. Thus, the overall computational efficiency of the VT-flash compositional flow simulation is enhanced under the multilayer nonlinear elimination method.
• #### Overview of the Adsorption and Transport Properties of Water, Ions, Carbon Dioxide, and Methane in Swelling Clays

(ACS Earth and Space Chemistry, American Chemical Society (ACS), 2021-08-26) [Article]
It is quite challenging to understand the details of the complex mechanisms involved in the swelling processes of clays such as montmorillonite as a function of relative humidity (RH), and the adsorption and transport processes of CO2 and CH4 in swelling clays. Here we present a review of the molecular simulation results for the swelling clay systems which not only compare well with the experimental data but also provide deep insights into the details of these mechanisms. The presence of CO2 and CH4 hardly affects the distribution and mobility of the interlayer water and ions in these systems. Under all conditions, the stable basal d-spacing was mainly determined by the type of counterion present in the interlayer region and the amount of water in each hydration state was almost independent of the RH and the layer charge. The uptake of CH4 in the 1W state of the Na-clay was much smaller compared to that of CO2. The adsorbate mobility generally increased with increasing hydration/RH because of the associated swelling of the interlayer region. Interestingly, the uptake of CO2 in the high-charge clay was dramatically decreased, and the mobility of CO2 in each hydration state was almost independent of the type of cation. The preferential adsorption of CO2 over CH4 plays an important role in the diffusion processes. Such an understanding is important for the successful mitigation of climate change via storage of anthropogenic CO2 in geological formations.
• #### Dynamics of ion depletion in thin brine films

(Fuel, Elsevier BV, 2021-08-26) [Article]
Low-salinity water flooding (LSW) effects are generated by the reduction of ionic concentration of environment electrolytes to which thin brine films confined between oil and rock are exposed. We study the dynamics of ion depletion from thin brine films upon a reduction of environment electrolyte concentration using the Poisson-Nernst-Planck (PNP) model. Interestingly, the model predicts that the timescale of ion depletion is not prolonged, but slightly shortened, by charged oil and rock surfaces in comparison with the absence of surface charges. This phenomenon is also reproduced quantitatively using a reduced ion depletion model inspired by the membrane science literature, in which salt diffusion and Donnan equilibrium between brine film and environment are considered. Furthermore, the self-diffusion of ions confined between n-decane and negatively charged quartz surface is investigated via atomistic simulations. It is found that, on average, the diffusion of ions in nanometer-thin brine films is slowed down up to ~8 times compared to that in bulk, although the slowdown relevant to ion depletion in event of a salinity reduction in the environment is most likely only about 2–3 times. These results provide new, pore-scale insights into LSW processes. The reduced salt depletion model and molecular simulation of ion diffusion demonstrated here help to develop a multiscale, bottom-up modeling framework for predicting LSW processes.
• #### Spectral relaxation computation of electroconductive nanofluid convection flow from a moving surface with radiative flux and magnetic induction

(Journal of Computational Design and Engineering, Oxford University Press (OUP), 2021-07-20) [Article]
A theoretical model is developed for steady magnetohydrodynamic viscous flow resulting from a moving semi-infinite flat plate in an electrically conducting nanofluid. Thermal radiation and magnetic induction effects are included in addition to thermal convective boundary conditions. Buongiorno's two-component nanoscale model is deployed, which features Brownian motion and thermophoresis effects. The governing nonlinear boundary layer equations are converted to nonlinear ordinary differential equations by using suitable similarity transformations. The transformed system of differential equations is solved numerically, employing the spectral relaxation method (SRM) via the MATLAB R2018a software. SRM is a simple iteration scheme that does not require any evaluation of derivatives, perturbation, and linearization for solving a nonlinear system of equations. Effects of embedded parameters such as sheet velocity parameter$\lambda$, magnetic field parameter$\beta$, Prandtl number$Pr$, magnetic Prandtl number$Prm$, thermal radiation parameter$Rd$, Lewis number$Le$, Brownian motion parameter$Nb$, and thermophoresis parameter$Nt$ on velocity, induced magnetic field, temperature, and nanoparticle concentration profiles are investigated. The skin-friction results, local Nusselt number, and Sherwood number are also discussed for various values of governing physical parameters. To show the convergence rate against iteration, residual error analysis has also been performed. The flow is strongly decelerated, and magnetic induction is suppressed with greater magnetic body force parameter, whereas temperature is elevated due to extra work expended as heat in dragging the magnetic nanofluid. Temperatures are also boosted with increment in nanoscale thermophoresis parameter and radiative parameter, whereas they are reduced with higher wall velocity, Brownian motion, and Prandtl numbers. Both hydrodynamic and magnetic boundary layer thicknesses are reduced with greater reciprocal values of the magnetic Prandtl number Prm. Nanoparticle (concentration) boundary layer thickness is boosted with higher values of thermophoresis and Prandtl number, whereas it is diminished with increasing wall velocity, nanoscale Brownian motion parameter, radiative parameter, and Lewis number. The simulations are relevant to electroconductive nanomaterial processing.
• #### Bound-preserving inexact Newton algorithms on parallel computers for wormhole propagation in porous media

(Computers and Geotechnics, Elsevier BV, 2021-07-12) [Article]
Simulating wormhole propagation during reactive dissolution of carbonates is often significantly challenging, because of the high nonlinearity of the governing equations with complex fluid physics. High resolution grids are often required to represent the complex geological heterogeneity, which demands massively parallel computers with a large number of processors. Herein, we present a parallel and scalable simulator on parallel computers for the fully implicit solution of the wormhole propagation model. Our approach is based on a family of mixed finite element methods for the spatial discretization and the implicit Backward Euler scheme for the temporal integration, to handle the combination of complicated flow physics and high resolution grids in their full complexity. Moreover, the active-set reduced-space method, as a class of bound-preserving inexact Newton algorithms, is proposed for the resultant nonlinear system of equations, to guarantee the nonlinear consistency of the fully implicit discretization in a monolithic way and to ensure the boundedness requirement of the solution. Numerical experiments are presented to demonstrate the effectiveness and efficiency of the proposed simulator for a set of heterogeneous medium problems. Large-scale results are provided to show the scalability for reservoir simulation with hundreds of millions of unknowns by using several thousand processors.
• #### Modeling of Water Generation from Air Using Anhydrous Salts

(Energies, MDPI AG, 2021-06-25) [Article]
The atmosphere contains 3400 trillion gallons of water vapor, which would be enough to cover the entire Earth with a one-inch layer of water. As air humidity is available everywhere, it acts as an abundant renewable water reservoir, known as atmospheric water. The efficiency of an atmospheric water harvesting system depends on the sorption capacities of water-based absorption materials. Using anhydrous salts is an efficient process in capturing and delivering water from ambient air, especially under a condition of low relative humidity, as low as 15%. Many water-scarce countries, like Saudi Arabia, receive high annual solar radiation and have relatively high humidity levels. This study is focused on the simulation and modeling of the water absorption capacities of three anhydrous salts under different relative humidity environments: copper chloride (CuCl2), copper sulfate (CuSO4), and magnesium sulfate (MgSO4), to produce atmospheric drinking water in water-scarce regions. By using a mathematical model to simulate water absorption, this study attempts to compare and model the results of the current computed model with the laboratory experimental results under static and dynamic relative humidities. This paper also proposes a prototype of a system to produce atmospheric water using these anhydrous salts. A sensitivity analysis was also undertaken on these three selected salts to determine how the uniformity of their stratified structures, thicknesses, and porosities as applied in the mathematical model influence the results.
• #### A self-adaptive deep learning algorithm for intelligent natural gas pipeline control

(Energy Reports, Elsevier BV, 2021-06-15) [Article]
Natural gas has been recognized as a promising energy supply for modern society due to its relatively less air pollution in consumption, while pipeline transportation is preferred especially for long-distance transmissions. A simplified pipeline control scenario is proposed in this paper to deeply accelerate the management and decision process in pipeline dispatch, in which a direct relevance between compressor operations and the inlet flux at certain stations is established as the main dispatch logic. A deep neural network is designed with specific input and output features for this scenario and the hyper-parameters are carefully tuned for a better adaptability of this problem. The realistic operation data of two pipelines have been obtained and prepared for learning and testing. The proposed algorithm with the optimized network structure is proved to be effective and reliable in predicting the pipeline operation status, under both the normal operation conditions and abnormal situations. The successful definition of "ghost compressors" make this algorithm to be the first self-adaptive deep learning algorithm to assist natural gas pipeline intelligent control.
• #### Modeling and Simulation of Atmospheric Water Generation Unit Using Anhydrous Salts

(Springer International Publishing, 2021-06-09) [Conference Paper]
The atmosphere contains 3400 trillion gallons of water vapor, which would be enough to cover the entire earth in 1 inch of water. Air humidity is available everywhere, and it acts as a great alternative as a renewable reservoir of water known as atmospheric water. Atmospheric water harvesting system efficiency depends on the sorption capacity of water based on the adsorption phenomenon. Using anhydrous salts is an efficient process for capturing and delivering water from ambient air, especially at a low relative humidity as low as 15%. A lot of water-scarce countries like Saudi Arabia have much annual solar radiation and relatively high humidity. This study is focusing on modeling and simulating the water absorption and release of the anhydrous salt copper chloride (CuCl2 ) under different relative humidity to produce atmospheric drinking water in scarce regions.
• #### Investigation of the dynamics of immiscible displacement of a ganglion in capillaries

(Capillarity, Yandy Scientific Press, 2021-06-03) [Article]
In this work the problem of displacing a ganglion of a fluid by another immiscible one in capillaries is investigated. A modeling approach is developed to predict the location of the ganglion with time. The model describes two patterns; namely, when the ganglion totally exists inside the tube, and when the advancing interface of the ganglion has broken through the exit of the tube. The model is valid for the case in which the ganglion is wetting as well as when it is nonwetting to the wall of the tube. It also considers the situation in which both the advancing and the receding interfaces assume, generally, different contact angles. For the special case when the displacement process is quasistatic, both the receding and the advancing contact angles may be considered the same. Under these conditions, interfacial tension force plays no role and the ganglion moves as a plug inside the tube with a constant velocity. When the viscosity ratio between the invading fluid and the ganglion is one (i.e., both phases are having the same viscosity) the motion reduces to the Hagen-Poiseuille flow in pipes. Once the advancing interface breaks through the exit of the tube, interfacial tension starts to take part in the displacement process and the ganglion starts to accelerate or decelerate according to the viscosity ratio. When the ganglion is nonwetting, interfacial tension becomes in the direction of the flow and is opposite to the flow otherwise. The model accounts for external forces such as pressure and gravity in addition to capillarity. A computational fluid dynamics analysis of this system is conducted for both types of wettability scenarios and shows very good match with the results of the developed model during both the two modes of flow patterns. This builds confidence in the developed modeling approach. Other cases have also been explored to highlight the effects of other scenarios.
• #### Construction of a Minimum Energy Path for the VT Flash Model by the String Method Coupled with the Exponential Time Differencing Scheme

(COMMUNICATIONS IN COMPUTATIONAL PHYSICS, Global Science Press, 2021-06) [Article]
Flash calculation plays significant roles in petroleum and chemical industries. Since Michelsen proposed his milestone studies in 1982, through several decades of development, the current research interests on flash calculation have been shifted from accuracy to efficiency, but the ultimate goal remains the same; that is accurate determination of equilibrium phase amounts and compositions at a given condition. On the other hand, finding the transition route and its related saddle point is often crucial to understand the whole energy landscape of flash models, which would provide new insights for designing numerical algorithms or optimizing existing ones. In this study, an efficient numerical approach is developed by coupling the string method with the exponential time differencing (ETD) scheme to investigate the minimum energy paths and first-order saddle points of VT flash models with Peng-Robinson equation of state. As a promising alternative to the conventional approach, VT flash calculates phase equilibria under a new variable specification of volume and temperature. The Rosenbrock-type ETD scheme is used to reduce the computational difficulty caused by the high stiffness of the model systems. The proposed ETD-String method successfully calculates the minimum energy paths of single-component and two-component VT flash models with strong stiffness. Numerical results also show good feasibility and accuracy in calculation of equilibrium phase amounts and compositions.
• #### Numerical modeling on hydrate formation and evaluating the influencing factors of its heterogeneity in core-scale sandy sediment

(Journal of Natural Gas Science and Engineering, Elsevier BV, 2021-06) [Article]
Natural gas hydrate (NGH) has been regarded as a fossil fuel reserve for the future on account of its tremendous potential. The numerical modeling on NGH formation/dissociation mechanism contributes to better understanding its accumulation and distribution feature, and optimizing the development program. This paper aims to develop a new simulator for the NGH formation in the core-scale sandy sediments based on the computational fluids dynamic (CFD) methods. The mathematical model is established based on the kinetic reaction model of hydrate formation, the permeability reduction model by the NGH, model of heat and mass transfer in porous media. The hydrate formation model is programmed by C language, and used as a subroutine for Fluent software which is adopted to solve the governing equations of the multiphase flow. The simulator scheme is verified by comparison with the experiment and numerical simulation in literature. What's more, this study reproduces the same fluctuant tendency of temperature as the experiment during the 1.0 h–2.0 h for the first time. Different reaction surface models of NGH formation/dissociation are evaluated by the developed codes. The effects of the reaction surface of hydrate (RSH) model and the initial fluids distribution on the hydrate formation process are simulated and analyzed. The variation of the RSH in NGH formation/dissociation should be taken into consideration when modeling the hydrate re-formation in the exploitation of NGH. The initial distribution of water and gas has a great impact on the hydrate formation in the sealed reactor. The hydrate distribution is ununiform, even when assuming the water and methane are mixed uniformly in a homogeneous porous media. This study provides new insight for the parametric estimation of the RSH model in the hydrate formation and dissociation modeling.
• #### Lunar features detection for energy discovery via deep learning

(Applied Energy, Elsevier BV, 2021-05-19) [Article]
Because of the impending energy crisis and the environmental Impact of fossil fuels, researchers are actively looking for alternatives, such as Helium-3 on the Moon. Although it remains challenging to explore energies on the Moon due to the long physical distance, the lunar features, such as craters and rilles, can be the hotspots for such energy sources, according to recent studies. Thus, identifying lunar features, such as craters and rilles, can facilitate the discovery of Helium-3 on the Moon, which is enriched in such hotspots. However, previously, no computational method was developed to recognize the lunar features automatically for facilitating space energy discovery. In our research, we aim at developing the first deep learning method to identify multiple lunar features simultaneously for potential energy source discovery. Based on the state-of-the-art deep learning model, High Resolution Net, our model can efficiently extract semantic information and high-resolution spatial information from the input images, which ensures the performance for recognizing the lunar features. With a novel framework, our method can recognize multiple lunar features, such as craters and rilles, at the same time. We also used transfer learning to handle the data deficiency issue. With comprehensive experiments on three datasets, we show the effectiveness of the proposed method. All the datasets and codes are available online.
• #### Sorption and Diffusion of Methane, Carbon Dioxide, and Their Mixture in Amorphous Polyethylene at High Pressures and Temperatures

(Industrial & Engineering Chemistry Research, American Chemical Society (ACS), 2021-05-12) [Article]
Molecular dynamics (MD) simulations are performed to study the sorption and transport properties of CH4 and CO2 in amorphous polyethylene at temperatures from 350 to 600 K and pressures up to 500 bar. The uptake of CH4 and CO2 by polyethylene generally increased with increasing pressure and decreasing temperature. However, at high pressures, for example, the uptake of methane by polyethylene increases with temperature. The self-diffusion coefficients of methane and carbon dioxide generally increase with pressure. These results are, in general, consistent with the swelling behavior of the polymer. Interestingly, for the penetrants, the activation barrier of diffusion decreases with pressure. MD simulations are also carried out for the CH4/CO2 mixture in amorphous polyethylene. Here, the overall sorption and transport properties were similar to those reported for pure CH4 and pure CO2 in polyethylene. The sorption selectivity of CO2/CH4 decreases with increasing pressure and temperature and was mostly independent of the bulk mole fraction of methane. Importantly, at high pressures, the mobility of methane found here is higher than that of the corresponding pure methane in polyethylene and the opposite trend is observed in the case of carbon dioxide. These results might be due to the fact that the swelling of the polymer in the presence of carbon dioxide is significantly higher than that in the presence of methane, especially at high pressures. The diffusion and membrane selectivities of carbon dioxide/methane show a similar trend to the sorption selectivity data. Furthermore, the simulation data were in good agreement with the theoretical calculations based on the PC-SAFT equation of state.
• #### Improved IMPES Scheme for the Simulation of Incompressible Three-Phase Flows in Subsurface Porous Media

(Energies, MDPI AG, 2021-05-11) [Article]
In this work, an improved IMplicit Pressure and Explicit Saturation (IMPES) scheme is proposed to solve the coupled partial differential equations to simulate the three-phase flows in subsurface porous media. This scheme is the first IMPES algorithm for the three-phase flow problem that is locally mass conservative for all phases. The key technique of this novel scheme relies on a new formulation of the discrete pressure equation. Different from the conventional scheme, the discrete pressure equation in this work is obtained by adding together the discrete conservation equations of all phases, thus ensuring the consistency of the pressure equation with the three saturation equations at the discrete level. This consistency is important, but unfortunately it is not satisfied in the conventional IMPES schemes. In this paper, we address and fix an undesired and well-known consequence of this inconsistency in the conventional IMPES in that the computed saturations are conservative only for two phases in three-phase flows, but not for all three phases. Compared with the standard IMPES scheme, the improved IMPES scheme has the following advantages: firstly, the mass conservation of all the phases is preserved both locally and globally; secondly, it is unbiased toward all phases, i.e., no reference phases need to be chosen; thirdly, the upwind scheme is applied to the saturation of all phases instead of only the referenced phases; fourthly, numerical stability is greatly improved because of phase-wise conservation and unbiased treatment. Numerical experiments are also carried out to demonstrate the strength of the improved IMPES scheme.
• #### Characterization and microfabrication of natural porous rocks: from micro-CT imaging and digital rock modelling to micro-3D-printed rock analogs

(Journal of Petroleum Science and Engineering, Elsevier BV, 2021-04-18) [Article]
Tests on standard rock specimens with controlled and identical pore structure are critical to validating the analytical and numerical models. However, it is usually difficult to acquire two natural samples with the same internal structure for the destructive laboratory tests, for the sake of the heterogeneity of natural rock which is caused by the complex diagenetic processes. Three-dimensional (3D) printing technology provides an alternative approach to produce geometry-identical, features-controllable, and lab-testable analogs of natural rock from digital data in a faster and more cost-effective way. This paper presents a customized workflow of 3D-printed rock analogs from micro-CT images combining with digital rock modelling. Three types of natural rock specimens are imaged by micro-CT and processed as inputs for two types of 3D printing techniques. Rock analogs are printed at multiple magnifications from original CT volume in five curable resin materials. Petrophysical parameters of 3D-printed rock analogs are acquired through helium pycnometry (HP) and mercury intrusion porosimetry (MIP). The accuracy of 3D-printed rock analogs is evaluated by comparing the measured results with the benchmark data derived from the digital rock modelling. Both the advantages and the current challenges to reproduce the real pore structure of natural rock by the 3D-printed analogs are discussed. The results indicate that the gypsum-based printed analogs are prior to modelling the surface roughness and wettability properties to natural rock grains, while the resin-based printed analogs owe advantages on reproducing pore structure. As the first effort in literature, this study investigates the inherent relationship between digital rock and 3D-printed rock analogs via comprehensive comparison on petrophysical properties. The results approve that the 3D printing technique is a novel, feasible, and alternative approach for laboratory test to generate rock analogs from the digital model of the natural rock. However, it is still difficult to print the pore structure of the rock at the original dimension.