Pore-scale Modelling of Gravity-driven Drainage in Disordered Porous Media
PPore-scale Modelling of Gravity-driven Drainage in Disordered Porous Media
Guanzhe Cui , Mingchao Liu
1, 2 , Weijing Dai and Yixiang Gan School of Civil Engineering, The University of Sydney, NSW 2006, Australia Department of Engineering Mechanics, CNMM & AML, Tsinghua University, Beijing 100084, China * Corresponding author, email: [email protected]
Abstract:
Multiphase flow through a porous medium involves complex interactions between gravity, wettability and capillarity during drainage process. In contrast to these factors, the effect of pore distribution on liquid retention is less understood. In particular, the quantitative correlation between the fluid displacement and level of disorder has not yet been established. In this work, we employ direct numerical simulation by solving the Navier-Stokes equations and using volume of fluid method to track the liquid-liquid interface during drainage in disordered porous media. The disorder of pore configuration is characterized by an improved index to capture small microstructural perturbation, which is pivotal for fluid displacement in porous media. Then, we focus on the residual volume and morphological characteristics of saturated zones after drainage and compare the effect of disorder under different wettability (i.e., the contact angle) and gravity (characterized by a modified Bond number) conditions. Pore-scale simulations reveal that the highly-disordered porous medium is favourable to improve liquid retention and provide various morphologies of entrapped saturated zones. Furthermore, the disorder index has a positive correlation to the characteristic curve index ( n ) in van Genuchten equation, controlling the shape of the retention characteristic curves. It is expected that the findings will benefit to a broad range of industrial applications involving drainage processes in porous media, e.g., drying, carbon sequestration, and underground water remediation. Keywords:
Multiphase flow, disordered porous media, gravity-driven drainage, Bond number, wettability.
Highlights:
1. The geometrical features can be quantitatively characterized by the newly defined disorder index πΌ " . 2. With increasing disorder, the formation of local fluid clusters dominates the morphological distribution and promotes higher total residual saturation. Quantitative correlations are derived for describing the fluid displacement with the disorder index. 3. For the small Bond number region (i.e., Bo<1.0), the relationship between the residual saturation and Bond number comply with van Genuchten relation. While for the high Bond number regime, adsorbed liquid dominates for the residual volume. INTRODUCTION
Drainage in porous media, a typical multiphase flow problem, plays a significant role in many environmental applications and industrial processes, including carbon capture and storage (CCS), oil and gas extraction, geological hazards, pharmaceutics, and food processes (Parker et al., 1987; Olivella et al., 1994; Bandara et al., 2011; Colombo and Fairweather, 2015; Zhao et al., 2016; Rabbani et al., 2017). As a wetting (e.g., liquid) phase is withdrawn and replaced by another non-wetting (e.g., gas) phase, the flow phenomena and entrapped saturated zone distribution are dominated by complex interactions among different constituents, including gas, liquid and solid, and attracting interests from different research areas related to unsaturated soil mechanics, groundwater remediation, and physics of porous media (MΓ©heust et al., 2002; Islam et al., 2014; Ghanbarian et al., 2015). The multiphase interactions combining several factors (including intrinsic topological features, gravity, capillary, and wettability) determine macroscopic drainage properties, such as the residual saturation of wetting phase and the temporal/spatial distributions of saturated zones (Yang et al., 1988; Herring et al., 2016; Li et al., 2017). The majority of previous experiments and numerical models of drainage and injection are focusing on the macroscopic parameters, e.g., porosity, permeability, system size and aspect ratio (Succi et al., 1989; Toussaint et al., 2005; Babadagli et al., 2015; Moura et al., 2015; Rognmo et al., 2017; March et al., 2018), in which the spatial configuration, pore connectivity and their influences on liquid retention are ignored (Prat, 1995; Lin et al., 2018; Yekta et al., 2018). However, for most natural and synthetic porous media, disordered microstructures are dominance (Anguy et al., 2001; Woo et al., 2004). These disorder features not only affect the mechanical properties of porous media (Laubie et al., 2017a; Laubie et al., 2017b), but also hinder the fundamental understanding of drainage processes (Holtzman, 2016; Fantinel et al., 2017). Further studies are required to evaluate the degree of microstructural disorder and bridge the microscopic information to macroscopic effective properties of porous media. Displacement efficiency of two-phase flow in disordered porous media is important in many applications (Dias and Payatakes, 1986). It has been found that the efficiency of immiscible fluids is influenced by the microstructure in those media (Zhao et al., 2016; Hu et al., 2018). For example, low disorder of porous geometry is found to be advantageous for fluid displacement in partially saturated porous media; on the other hand, in terms of liquid mixing and reaction, highly-disordered porous materials is desirable (Holtzman, 2016). These applications indicate that the disordered microstructure and its effects on liquid-liquid displacement need to be quantitatively assessed. However, the quantitative and systematic studies about the effects of disordered microstructure on the multiphase flow problems, particularly on liquid retention during drainage processes, remain relatively scarce. In the past few years, increasing attention has been focused on the effects of disordered geometry on immiscible two-phase flow (Ferrari et al., 2015; Holtzman, 2016; Rabbani et al., 2017). The ordered and disordered porous media have clear differences in the behaviours and mechanisms of multiphase fluid displacement processes (Pak et al., 2015). With an increase of the disorder degree, the displacement stability is weakened, which results in a lower degree of saturation than that in relatively more homogeneous pore network (Liu et al., 2015). Ferrari et al. (2015) studied numerically the effect of different disordered microstructure controlled by a standard deviation of solid phase size on the unstable invasion structure in Hele-Shaw cells consisting of a large number of cylindrical obstacles. The quantitative assessment on the impact on the disorder of pore distribution keeps relatively unexplored. In this paper, we employed the direct numerical simulation (DNS) by using the volume of fluid (VOF) method based on Eulerian algorithms to capture the detailed flow dynamics within the porous media, which describes the distribution of two phases as a fluid function with different properties. Then, quasi-two-dimensional (2D) simulations of gravity-driven drainage are performed on a broad range of scenarios and the temporal and spatial information of saturated clusters is exacted and analysed. Our objectives are twofold: firstly, we investigate the reliability of an introduced disorder index to characterize a quantitative relation between the degree of microstructural disorder and retention characteristics; secondly, we demonstrate the effect of gravitational force and wettability on drainage consequence under different disorder degree conditions. NUMERICAL MODELING
In this section, we introduce the numerical framework for generating and characterizing the disordered porous media, via an improved disorder index using Voronoi tessellations of the pore space, and modelling the gravity-driven drainage processes employing VOF method in
OpenFOAM . Disorder index, π° v To quantitatively investigate the influence of the degree of microstructural disorder in porous media on the drainage process, we explore a meaningful measure of the disorder of the pore space, as well as generating samples with the controlled disorder quantity. Here a model 2D porous medium containing mono-sized circular discs is focused, seen in Fig. 1(a). Random samples can be achieved by deviating the discs from the original regular arrangement, e.g., a hexagonal or square lattice. The porous medium is in form of π disk-shaped solid phases of radius π that is embedded in a rectangular plate of size πΏ β Γ πΏ ) with the periodic boundary condition along the π₯ -axis. The centers of discs are moved with a random vector π + ,,,,,β ( .π + ,,,,,β. β€πΏ ). The disordered porous media contain only non-overlapping discs at a constant global porosity, π , are constructed and the corresponding degree of disorder can be controlled by varying πΏ . The movement method of discs allows a larger range of configuration from the completely ordered system, πΏ = 0 , to highly-disordered systems by applying the moving steps iteratively. Fig. 1.
Schematic diagrams of (a) a typical computational domain used for the drainage process containing solid (grey), liquid (red) and gas (blue) phases, (b) Voronoi tessellations on the pore space and three characteristic liquid zones (saturated cluster, adsorbed droplet and liquid bridge), (c) VOF method with the volume fraction of wetting liquid inside each square grid.
Since the porosity π as a global parameter is determined by the number and size of discs in the domain, the next high-order parameter is a measurement of the fluctuations of local porosity π , e.g., πΌ = [β©π βͺ β π ] < 8β as the index of disorder proposed by Laubie et al. (2017a). Note that the average operator β©ββͺ is volumetric average. The local porosity π is evaluated at every point i in a square background mesh of the domain with a mesh size of π = (ππ /π)
The comparison between disordered indices πΌ (original definition) and πΌ D (improved index). Inserts (a)-(d) indicate the different background meshes used for calculating the local porosity: (a) πΌ =0 and (b) πΌ D =0.02, (c) πΌ =0.043 and (d) πΌ D =0.037, where (a-b) and (c-d) are two cases with identical pore structure. Fig. 2 shows the comparison between two disorder indices, πΌ and πΌ D . When the ordered primitive structure is slightly disturbed, e.g., πΏ β€ π 2β β π , πΌ = 0 , while πΌ D varies from 0 to 0.03. The reason is that, when the circular discs only randomly move within the square meshes, the local porosity π equals to mean global porosity π , as shown in Fig. 2(a)-(d) of the two identical microstructures. This indicates that πΌ cannot be used to characterize the small variation in a unit cell (as mentioned as βRandomness Case Iβ in Laubie et al. (2017b)); on the contrary, the improved definition using Voronoi tessellation is sensitive to such small disturbances. However, this limited movement has a quite significant impact on the drainage processes, which will be demonstrated later in Section 3. As a result, for investigating the pore-scale multiphase flow in this paper, we employed the πΌ D based on Voronoi tessellations to evaluate the degree of microstructural disorder in porous media. Moreover, the resulting ranges of πΌ and πΌ D also differ from πΌ β [0,0.1] and πΌ D β [0,0.06] , respectively. Volume of Fluid
In traditional computational fluid dynamics, during the procedure of solving the Navier-Stokes equations, it is extremely difficult and computationally expensive to solve moving interfacial boundary condition problems (Horgue et al., 2013). However, the VOF method employs treatments of two immiscible and isothermal, incompressible liquid phases as a single-fluid phase with different properties including density and viscosity. It utilizes an additional interfacial force replacing the jump condition at the liquid-liquid interface (Owkes and Desjardins, 2014). In VOF method, the process of drainage in porous media can be described as a wetting phase (βdefending fluidβ) replaced by a non-wetting phase (βinvading fluidβ) (Ferrari and Lunati 2013). The spatial distribution of two immiscible fluid phases is described as the phase indicator for fluids, shown in Fig.1 (c), as πΌ = O 0 1 in the non-wetting ( nw ) fluid , in the wetting ( w ) fluid . (2) According to the dynamic law of fluid phase, the immiscible fluids can be expressed by the conversation of mass and momentum as a set of Navier-Stokes equations: π» β π =0 , (3) RSπRT + π» β (πππ) = βπ»π + π» β (2ππ¬) + π [ , (4) where π is the fluid velocity, π is the pressure, E = <8 ( β u+ βπ’ ^ ) is the strain rate tensor, π and π are the density and viscosity considering multiple fluid phases, respectively, as: π(π) = πΌπ β + (1 β πΌ)π aβ , (5) π(π) = πΌπ b + (1 β πΌ)π cb , (6) In Eq. (4), the last term represents the effect of Laplace pressure at the interface and is defined as: π d = 2ππππΏ h , (7) where π is the interfacial tension between the two fluids, πΏ ^ is a unit function that is 0 in the whole domain except at the interface, π is the normal to the interface, and π is the mean curvature at the interface: π = <8 π» β π = <8 π» β i jkβjkβ m . (8) In this work, the no-slip condition is adopted, i.e., the vectorial component of velocity is zero on the solid boundary. To reproduce the wetting behaviour at the triple line region, the boundary condition in VOF equals to Youngβs law as defined as (Ferrari and Lunati, 2013): π| o ᴦ = π π¬ sinπ + π π¬ cosπ , (10) where π is the normal to two fluids interface at the solid surface, π π¬ is the unit tangent toward the wetting phase, π π¬ is the unit normal pointing into the solid, and π is the contact angle. When the above equations are used to compute the whole-domain velocity field and describe the evolution of the fluid-fluid displacement, the advection equation should be included as RkRT + π» β (πΌπ) + π» β (π y πΌ(1 β πΌ)) = 0 , (11) where π y ( π y = π < β π ) is the relative velocity between the two fluids. Since OpenFOAM use the surface compression method to implement the VOF method, the second term on the left hand side is an artificial compression term applied to prevent the numerical diffusion which smears the sharp interface (Boyce et al., 2016).
Simulation setup and materials parameters
The numerical simulations are implemented in an open-source CFD platform,
OpenFOAM , with a multiphase flow solver interFoam . This solver employs Finite Volume (FV) schemes for the discretization of partial differential equations to calculate multiple immiscible fluids phases. Inside the pore space, the classical Navier-Stokes equations are solved by VOF method tracking the temporal and spatial evolution of fluid-fluid interface. A mesh generator, snappyHexMesh , is used to automatically generate complex meshes of hexahedral and split-hexahedral cells from triangulated surface geometry in the domain. In order to construct a quasi-2D model, the out-of-plane thickness of the sample is particularly set to be much smaller than the in-plane edges in x and y axis. The initial background hexagonal mesh is iteratively refined to better reproduce the cylindrical surface with the meshes. To investigate the morphological characteristics of entrapped wetting fluid, including absorbed liquid, liquid bridge and saturated cluster, as shown in Fig. 1(b), we simulated drainage under atmospheric condition in a rectangular domain (Fig. 1(a)) with inlet on top, outlet on bottom and periodic boundary condition in the horizontal direction. The geometrical parameters of computational models are shown in Table 1. The dynamics of fluid-fluid displacement, morphology and distribution of wetted zone during gravity-driven drainage depend on the properties of wetting and non-wetting phases including density, viscosity and surface tension, which are selected based on the parameters of water and air (see Table 1).
Table 1.
Geometrical parameters of the quasi-two-dimensional porous media and properties of two immiscible fluids used in numerical simulations.
Property
Value
Radius of discs, π (mm) 1 Contact angle, π ( Β° ) ~150 Β° Viscosity of wetting phase (m /s) 5.93Γ10 -7 Viscosity ratio, π b /π cb (-) 0.38 Density of wetting phase, π β (kg/m ) 997 Density of non-wetting phase, π a (kg/m ) 11.69 Density difference, βπ (kg/m ) 985.31 Void ratio, e (-) 2.65 Gravity acceleration, g (m/s ) 9.81 Surface tension, πΎ (N/m) 0.064 Table 2.
A sensitivity study on the domain size based on the same discretization ( =27).
Small Medium Large Dimensions (mm ) 35 Γ Number of discs (-) 100 400 1600 Number of meshes (-) 825,547 2,975,769 11,990,042 Porosity (-) 0.704 0.705 0.704 Degree of residual saturation 17.73% 16.90% 14.94%
Table 3.
A sensitivity study on discretization based on the medium-sized domain.
Discretisation,
16 27 32 41 Number of meshes (-) 1,080,298 2.975,769 5,268,631 6,662,617 Porosity (-) 0.704 0.705 0.705 0.708 Degree of residual saturation 15.86% 16.90% 17.11% 16.93%
Before the detailed simulation of drainage processes, we have performed sensitivity studies on the mesh discretization and domain size to identify the cost-efficient mesh density and representative model size. A compromise between the level of discretization and size of computational domains has to be found. To this end, we tested several cases including three computational domains and four levels of discretization, defined by the ratio of disc diameter to typical mesh cell size (i.e., =16, 27, 32, 41). The benchmark is defined in terms of the residual saturation at the end of gravitational drainage process, and all cases are conducted under the same porosity and pore structure, with π =30 Β° and Bond number of 0.36. Based on these results (shown in Table 2 and Table 3), we decided to select the case with medium domain size of 70 Γ with 400 solid discs and discretisation of =27 (meshes), which ensures that differences with the results obtained with the finest mesh discretization and largest domain size are within 3%. RESULTS ANALYSIS
By employing the proposed numerical model in the last section, we can investigate the gravity-driven drainage processes. More specifically, in this section, the effects of degree of disorder, gravity and wettability on the residual volume and morphology distribution of wetting phase will be discussed.
Effects of disorder index, π° π― Fig. 3.
Simulation results of four different samples with increasing disorder index: (a) πΌ D =0, (b) πΌ D = 0.0056, (c) πΌ D =0.015 and (d) πΌ D =0.055. (e) Liquid-holding capacity (the residual volume of wetting phase normalized by the solid volume as a function of disorder index πΌ D . The scale shows the range of liquid-holding capacity varying from 0.41 to 0.68. We first demonstrate the reliability of disorder index πΌ D by characterizing the effect of disorder degree of pore topology on the liquid-holding capacity defined as the residual volume of wetting phase normalized by the solid volume. Note that, in order to exclude upper and lower boundary effects in the calculation, the volume of residual liquid within the top one row and bottom two rows of discs are removed during image processing. Fig. 3(a)-(d) show the spatial distributions of liquid zones (i.e., the entrapped saturated zones) in four samples with different pore topologies after drainage, which correspond to disorder index: πΌ D =0 (ordered distribution), πΌ D =0.0056, πΌ D =0.015 and πΌ D =0.055 (highly-disordered topology), respectively. The morphological characteristics of entrapped saturated zones are significantly different in the four samples. Such variations in the morphology of the entrapped saturated zones at pore-scale will be analysed in detail in Section 3.2. In the macroscopic scale, the numerical results demonstrate that the degree of disorder has a monotonously increasing relationship with the volume of entrapped saturated zones, as shown in Fig. 3(e), in which the solid line is for a guiding purpose. It is obvious that πΌ D varying in the range from 0 to 0.02 has a more significant impact on the residual volume than the rest range. Absorbed droplets are rapidly connected into liquid bridges and even saturated clusters with the increasing degree of disorder within this range, which are more efficient to enhance the liquid-holding capacity than the other πΌ D range that mainly contains clusters. Note that πΌ d in the corresponding range keeps nearly zero when πΌ D vary from 0 to 0.02 (see Figure 2), failing to capture this transition region. These results indicate that the parameter πΌ D is effective to quantitatively evaluate the influence of the topological disorder on the liquid holding capacity. In addition, the increase of the disorder index is unfavourable for drainage efficiency in porous media. Distribution of residual saturation cluster
In this section, the study aims to find a relationship between the disorder degree and the morphological characteristics of entrapped saturated zones after drainage. The spatial information of the saturated clusters is extracted by image processing. In statistics, the area of liquid bridges range from 0.8 mm to 4 mm , and other zones of smaller or larger areas correspond to absorbed droplets and saturated clusters, respectively, with typical examples shown in Fig. 1(b). In the post-processing, the extremely small liquid volume with a size smaller than one-tenth of a solid disc area ( π mm ) volume is neglected. The probability density distribution of saturated zones volume after drainage processes are plotted in Fig. 4. The distributions of zone area using a uniform bin-width of 0.4 mm are extracted from four samples with different disorder degree (i.e., πΌ D =0, πΌ D =0.0056, πΌ D =0.015 and πΌ D =0.055), corresponding to the pore topologies in Fig. 3(a)-(d), respectively. Fig. 4(a) shows that, when the solid discs are assigned regularly ( πΌ D =0), the maximum probability of cluster area is possessed by the bin from 0.4 mm to 0.8 mm , which is mostly consisted of the absorbed droplets. In Fig. 4(b) and Fig. 4(c), the span of the size distribution of saturated zones becomes wider with increasing disorder of pore structure. Specifically, the three kinds of wetting phase morphology (absorbed droplet, liquid bridge and saturated cluster) simultaneously exist at the final drainage stages. When disorder index arrives the approximately maximum value ( πΌ D =0.055), the saturated zones in the highly-disordered porous media have the maximum cluster area reaching 40 mm as well as the minimum number of absorbed droplets and liquid bridges (see Fig. 4(d)). These results demonstrate that the increase of disorder contributes to the aggregation of solid discs, which is advantageous for the merging of absorbed droplets into liquid bridges and further into residual saturated clusters. Fig. 4.
Probability density function (PDF) and lognormal distribution of saturated zones area and maps of morphological characteristics with the four disorder parameters: (a) πΌ D =0, π =-0.061, π =0.59; (b) πΌ D =0.0056, π =0.24, π =0.85; (c) πΌ D =0.015, π =0.38, π =1.04; (d) πΌ D =0.055, π =0.66, π =1.20. Note that the different colours in the inserted morphological character map represent the disconnected wetting phase zones. The dataset can be characterized by employing the lognormal distribution function, which has been employed in interpreting experimentally observed saturation clusters (Li et al. 2017) as π(π|π, π ) = <(cid:136)(8(cid:137)(cid:138) (cid:139) )(cid:140) ππ₯π (cid:142)β ((cid:143)a(cid:140)(cid:144)(cid:145)) (cid:139) (cid:139) (cid:146) . (12) According to the maximum likelihood estimation method, two determinative parameters, i.e., the scale parameter ( π ) and shape parameter ( π ), can be obtained by fitting the probability density distribution of liquid zone areas. It is clearly seen from Fig. 4(a)-(d) that the mean π and variance π have a positive relationship with the disorder index, indicating the increase of the mean and variation of residual liquid areas of high degree of disorder by enhancing the formation of saturated clusters. The reason is that a large πΌ D improves the solid phase aggregation, individual liquid bridges and absorbed droplets merging and eventually leading to the formation of large saturated zones in those dense regions. This further illustrates why the liquid-holding capacity is gradually enhanced with the increase of disorder degree and how the disorder of porous media statistically varies the total volume and morphological characteristics of residual wetting phase after the gravity-driven drainage process. Liquid retention under varying Bond number
During the gravity-driven drainage process, the interplay between the gravity and capillary forces has a significant effect on the fluid displacement behaviours. In order to quantify the relative magnitude of these two forces, a dimensionless parameter, i.e., the Bond number, is introduced in the previous work (MΓ©heust et al., 2002; Moebius and Or, 2014) as π΅ (cid:148) = βS(cid:149)(cid:150) (cid:139) (cid:151) , (13) where π is the characteristic length, βπ density difference between wetting and non-wetting phases, π the gravity acceleration, and πΎ the surface tension. The characteristic length is usually defined as the throat size (Blunt and Scher, 1995) or the radius of pore (Birovljev et al., 1991; LΓΈvoll et al., 2005). However, because the disorder degree introduces non-uniform solid disc separations, the throat size is not easy to quantify. Here we adopt the uniform radius of solid discs as the characteristic length. In order to consider the variation due to the void ratio, we introduce the volumetric ratio between the pore and solid phase into the traditional form Eq. (13) and obtain the modified Bond number: π΅ (cid:148)(cid:154) = βS(cid:149)(cid:155) (cid:139) (cid:151) β π , (14) where R is the radius of the solid disc and π = π (cid:157) /π (cid:158) is the void ratio, in which π (cid:157) and π (cid:158) are the total volume of pore and solid phase, respectively. To demonstrate the validity of the modified Bond number, the parameters including the density difference, the gravity acceleration, the radius of solid disc, the void ratio are individually adjusted with the surface tension to keep the B o' the same at 0.37. Fig. 5 demonstrates that the degree of saturation is similar when the B o' is kept the same for given contact angle, across a range of wettability. In the following analyses, we will focus on the combined effects among the modified Bond number, disorder and wettability. Fig. 5.
Degree of residual saturation as a function of wettabilities at B o' =0.37 with a different combination of density difference ( βπ , kg/m ), gravity acceleration ( π , m/s ), disc radius ( π , mm), void ratio ( π ) and surface tension ( πΎ , N/m). To elucidate the combined effect of disorder degree and the modified Bond number, B o' , simulations are performed for the four disordered samples as shown in Fig. 3(a)-(d) with a broad range of B o' from 0.011 to 7.70. The liquid-holding capacity of the disordered samples is plotted in Fig. 6 as a function of B o' . The results demonstrate that, based on the influence of disorder, the drainage process can be divided into two categories, the capillary regime and adsorbed regime. The B o' threshold separating these two categories is approximately 1.0, beyond which the gravitational force overtakes capillary force in drainage process, resulting in the disappearance of saturation clusters and capillary bridges. When B o' is less than the threshold, capillarity dominates the fluid-fluid displacement. Based on comparison between the ordered and the highly disordered medium (see Fig. 6(a) and 6(b)), higher disorder can dramatically improves the liquid-holding capacity of porous media with the same B o' ; in contrast to the only small absorbed droplets in the ordered system, high disorder leads to large entrapped saturated zones, which demonstrates that disorder degree has a significant effect on the capillary interaction during drainage. In addition, it can be clearly seen that there is a morphological transformation from saturated zones (Fig. 6(a)) to liquid bridge (Fig. 6(c)) with the increasing gravitational force. In the second regime, i.e. B o' >1.0, the gravitational force has overcome the capillary counterpart. Fig. 6(d) indicates that only the absorbed droplets exist on the surface of solid discs instead of the saturated patch and liquid bridges between discs in Fig 6(b) and 6(c), respectively. The absorbed volume gradually declines with the increase of Bond number. Since the absorbed liquid phase is only proportional to the total surface area (length in 2D) of the solid phase, there is no observed difference for ordered and disordered pore spaces, exhibiting similar relationship as shown in Fig. 6. Fig. 6.
Liquid-holding capacity as a function of Bond number with four different disorder indices and maps of wetted cluster distribution at contact angle ΞΈ =30 Β° : (a) πΌ D =0.055, Bo ' =0.073, (b) πΌ D =0, Bo ' =0.073, (c) πΌ D =0, Bo ' =1.10, (d) πΌ D =0.027, Bo ' =7.65. The van Genuchten equation originally describes the relationship between the degree of saturation and matric suction in the soil-water characteristics curves (SWCCs) (Zhou et al., 2016). As an analogue, here the residual volume is governed by the interaction between the capillary and gravitational forces (characterized by Bond number). Thus, we replaced the suction term in the original format with the Bond number οΌ π y = <(cid:144)(cid:160) ‘’ Β£<β( B o' /k) Β₯ Ζ Β§Β€Β§/Β₯ + π y' , (15) where π y is the degree of saturation, π ya the adsorbed saturation, πΌ the scaling parameter positively proportional to the air-entry value, and π is the characteristic curve index controlling the slope of SWCC (Fredlund et al., 2012). Fig. 7 shows that the liquid retention curve can be quantitatively predicted by the van Genuchten equation (i.e., Eq. (15)). Furthermore, we conclude that with the same porosity and solid disc diameter, the increase of disorder index πΌ D decreases the curve index n . During the curve fitting process, we used πΌ =0.055 and π βΉ(cid:150) =0.20, since the former was found to be insensitive to πΌ D and latter can be measured at B o' = 1 . Thus, the potential relationship between the disorder of pore space and liquid characteristic curves can be evaluated quantitatively by the application of πΌ D under the framework of unsaturated soil mechanics. Fig. 7.
Degree of residual saturation as a function Bond number fitted by van Genuchten equation (with πΌ =0.055, π βΉ(cid:150) =0.20): (a) πΌ D =0, n =16.07, (b) πΌ D =0.0056, n =6.80, (c) πΌ D =0.015, n =3.75, (d) πΌ D =0.055, n =2.54. Surface wettability
It is well known that the wettability of porous media can dramatically change the liquid-liquid displacement behavior during both drainage and imbibition processes (Zhao et al., 2016), but whether the disorder alters such influence of wettability has yet been explored. To this end, a series of simulations were performed on the disordered porous samples with the wettability of solid surface varying over a broad range from 0 Β° to 150 Β° at the same Bond number of 0.37. The results are plotted in Fig. 8. It is clearly seen that, with the increase of disorder index πΌ D , the volume of liquid retention is gradually elevated at the same wettabilities, which indicates the influence of the wettability on retention can be controlled by adjusting the disorder index (see Fig. 8(a) and 8(b)). Meanwhile, as for varying wettability, it is relatively more efficient to alter the residual volume when the disorder index is within πΌ D β (0-0.02), than the larger degree regime, πΌ D β (0.02-0.06). Due to the disorder-induced reinforcement of liquid-holding capacity, the maximum and minimum relative difference, compared to the ordered configuration at the same wettability, are 108% at the hydrophilic ( ΞΈ =60 Β° ) and 716% at hydrophobic ( ΞΈ =15 ) samples, respectively. Figure. 8.
Liquid-holding capacity as a function of wettability for the different disorder degree and morphological patterns of drainage results with different wettability conditions as inserts: (a) πΌ D =0.055, π =10 Β° (b) πΌ D =0, π =10 Β° ; and (c) πΌ D =0.055, π =110 Β° . In addition, when the surface wettability of porous media varies from hydrophilic to hydrophobic, wetting phase begins to flow out more and the liquid-holding capacity drastically decreases until approximately the inflection point (i.e., ΞΈ~ Β° ). During this transition, in terms of high disorder degree, the morphological characteristics of final drainage clearly transforms from a variety of liquid clusters to only entrapped droplets between the circular solid phase (see Fig. 8(c)). Such entrapment can be noticeably strengthened in high disorder samples and eventually keeps the residual volume of wetting phase nonzero. When the contact angle is within 0 Β° to 110 Β° , the effect of wettability on drainage efficiency is more obvious. DISUCUSSIONS
Our numerical results have shown that the disorder parameter πΌ D is reliable to quantitatively characterize the influence of disorder degree on liquid retention in porous media after the gravity-driven drainage process. The residual volume has a monotonously increasing relationship with the disorder index under the same porosity and radius of solid discs. It demonstrates that highly-disordered media is disadvantageous for drainage efficiency, however, the efficiency can be more sensitively controlled within the range of disorder index πΌ D (0-0.02) than the larger range (0.02-0.06). For representing the morphological characteristics, we employed the lognormal function to fit the probability density distribution of saturated areas. The scale parameters π and shape parameter π are positively correlated with disorder index πΌ D , which indicates that enhancing disorder is advantageous for the formation of large fluid clusters. Due to the reinforcement of solid phase aggregation, a large quantity of absorbed droplets and liquid bridges can be gradually connected into liquid clusters. In addition, a modified Bond number, including the void ratio, is introduced. Based on the magnitude of gravitational and capillary forces, the drainage process is divided into two categories and the transition point of is approximately B o' =1.0. When the gravitational force dominates the drainage process (B o' >1.0), only absorbed droplets exist on the solid surface without the presence of saturated clusters and liquid bridges, resulting in almost the same liquid-holding capacity for different disorder degree; on the other hand, when capillary force is larger than gravitational force (B o' <1.0), the relationship between the degree of saturation and Bond number can be well represented by van Genuchten equation, in which disorder index πΌ D has a positive relationship with the characteristic curve index n . Finally, the study of fluid drainage characteristics with various surface wettability demonstrates that the effect of microstructural disorder on residual volume is clearly weakened with the increase of hydrophobicity. In addition, when the surface wettability vary from hydrophilic to hydrophobic, morphological characteristics of highly-disordered topology have an obvious transformation from various saturated zones to only entrapped droplets. In summary, this paper provides a fresh perspective into the significant applications of two-phase flows, in which displacement efficiency and morphological feature of wetting phase are crucial. For example, decreasing disorder degree can control the drainage processes, desirable for improving the displacement efficiency of water during drying of food and soil aeration processes, and accelerating the recovery of contaminants and CO geosequestration. On the other hand, increasing disorder enlarges the interfacial area, which is advantageous for the subsurface remediation, fluid-gas mixing, and reaction enhancement in microfluidics.
5. CONCLUSIONS
Based on Voronoi tessellations of the pore space, the disorder index πΌ D has been introduced to describe the drainage characteristics in disordered porous media. With the help of this index, we studied the gravity-driven drainage with different pore configurations exhibiting a broad range of disorder degree. The results demonstrate that πΌ D is able to quantitatively characterize the influence of the disorder of porous media on the residual volume and morphological characteristics of residual saturated zones. Among other factors (capillarity, gravity, wettability, porosity and characteristic size), this study highlights the correlations between the disorder index to the following drainage behaviour: (1) residual saturation, (2) morphological features and (3) statistical distribution of residual wetting phase. The results presented in this paper contribute to fundamentally understand the physical principles behind multiphase flow in disordered porous media and to improve applications related to control and optimize the drainage processes. ACKNOWLEDGEMENTS
This work was financially supported by Australian Research Council (Projects DP170102886) and The University of Sydney SOAR Fellowship.
REFERENCES
Anguy, Y., et al. (2001). "On the ability of a class of random models to portray the structural features of real, observed, porous media in relation to fluid flow." Cement and Concrete Composites (2-3): 313-330. Babadagli, T., et al. (2015). "Effect of surface roughness and lithology on the waterβgas and waterβoil relative permeability ratios of oil-wet single fractures." International Journal of Multiphase Flow : 68-81. Bandara, U. C., et al. (2011). "Pore-scale study of capillary trapping mechanism during CO2 injection in geological formations." International Journal of Greenhouse Gas Control (6): 1566-1577. Birovljev, A., et al. (1991). "Gravity invasion percolation in two dimensions: Experiment and simulation." Physical Review Letters (5): 584. Blunt, M. J. and H. Scher (1995). "Pore-level modeling of wetting." Physical Review E (6): 6387. Boyce, C., et al. (2016). "Intrusion of a Liquid Droplet into a Powder under Gravity." Langmuir (34): 8631-8640. Colombo, M. and M. Fairweather (2015). "Multiphase turbulence in bubbly flows: RANS simulations." International Journal of Multiphase Flow : 222-243. Dias, M. M. and A. C. Payatakes (1986). "Network models for two-phase flow in porous media Part 1. Immiscible microdisplacement of non-wetting fluids." Journal of Fluid Mechanics : 305-336. Fantinel, P., et al. (2017). "Drying in a microfluidic chip: experiments and simulations." Scientific Reports (1): 15572. Ferrari, A., et al. (2015). "Challenges in modeling unstable two β phase flow experiments in porous micromodels." Water Resources Research (3): 1381-1400. Ferrari, A. and I. Lunati (2013). "Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy." Advances in Water Resources : 19-31. Fredlund, D. G., et al. (2012). Unsaturated soil mechanics in engineering practice, John Wiley & Sons. Ghanbarian, B., et al. (2015). "Gas and solute diffusion in partially saturated porous media: Percolation theory and effective medium approximation compared with lattice Boltzmann simulations." Journal of Geophysical Research: Solid Earth (1): 182-190. Greenshields, C. J. (2015). "Openfoam user guide." OpenFOAM Foundation Ltd, version (1). Herring, A. L., et al. (2016). "Impact of wettability alteration on 3D nonwetting phase trapping and transport." International Journal of Greenhouse Gas Control : 175-186. Holtzman, R. (2016). "Effects of pore-scale disorder on fluid displacement in partially-wettable porous media." Scientific Reports : 36221. Horgue, P., et al. (2013). "Experimental and numerical study of two-phase flows in arrays of cylinders." Chemical Engineering Science : 335-345. Hu, R., et al. (2018). "Wettability and flow rate impacts on immiscible displacement: A theoretical model." Geophysical Research Letters (7): 3077-3086. Islam, A., et al. (2014). "Characterization of the crossover from capillary invasion to viscous fingering to fracturing during drainage in a vertical 2D porous medium." International Journal of Multiphase Flow : 279-291. Laubie, H., et al. (2017b). "Disorder-induced stiffness degradation of highly disordered porous materials." Journal of the Mechanics and Physics of Solids : 207-228. Laubie, H., et al. (2017a). "Stress transmission and failure in disordered porous media." Physical Review Letters (7): 075501. Li, S., et al. (2017). The distribution of saturated clusters in wetted granular materials. EPJ Web of Conferences, EDP Sciences. Lin, J., et al. (2018). "Experimental investigation of N 2 injection to enhance gas drainage in CO 2-rich low permeable seam." Fuel : 665-674. Liu, H., et al. (2015). "Lattice Boltzmann simulation of immiscible fluid displacement in porous media: Homogeneous versus heterogeneous pore network." Physics of Fluids (5): 052103. LΓΈvoll, G., et al. (2005). "Competition of gravity, capillary and viscous forces during drainage in a two-dimensional porous medium, a pore scale study." Energy (6): 861-872. March, R., et al. (2018). "Assessment of CO2 Storage Potential in Naturally Fractured Reservoirs With Dual β Porosity Models." Water Resources Research. MΓ©heust, Y., et al. (2002). "Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects." Physical Review E (5): 051603. Moebius, F. and D. Or (2014). "Pore scale dynamics underlying the motion of drainage fronts in porous media." Water Resources Research (11): 8441-8457. Moura, M., et al. (2015). "Impact of sample geometry on the measurement of pressure β saturation curves: Experiments and simulations." Water Resources Research (11): 8900-8926. Olivella, S., et al. (1994). "Nonisothermal multiphase flow of brine and gas through saline media." Transport in Porous Media (3): 271-293. Owkes, M. and O. Desjardins (2014). "A computational framework for conservative, three-dimensional, unsplit, geometric transport with application to the volume-of-fluid (VOF) method." Journal of Computational Physics : 587-612. Pak, T., et al. (2015). "Droplet fragmentation: 3D imaging of a previously unidentified pore-scale process during multiphase flow in porous media." Proceedings of the National Academy of Sciences (7): 1947-1952. Parker, J., et al. (1987). "A parametric model for constitutive properties governing multiphase flow in porous media." Water Resources Research (4): 618-624. Prat, M. (1995). "Isothermal drying on non-hygroscopic capillary-porous materials as an invasion percolation process." International Journal of Multiphase Flow (5): 875-892. Rabbani, H. S., et al. (2017). "New insights on the complex dynamics of two-phase flow in porous media under intermediate-wet conditions." Scientific Reports (1): 4584. Rognmo, A., et al. (2017). "Nanotechnology for improved CO2 utilization in CCS: Laboratory study of CO2-foam flow and silica nanoparticle retention in porous media." International Journal of Greenhouse Gas Control : 113-118. Succi, S., et al. (1989). "Three-dimensional flows in complex geometries with the lattice Boltzmann method." EPL (Europhysics Letters) (5): 433. Toussaint, R., et al. (2005). "Influence of pore-scale disorder on viscous fingering during drainage." EPL (Europhysics Letters) (4): 583. Woo, H.-J., et al. (2004). "Modeling desorption of fluids from disordered mesoporous materials." Langmuir (11): 4743-4747. Yang, Y.-W., et al. (1988). "Capillary flow phenomena and wettability in porous media: I. Static characteristics." Journal of Colloid and Interface Science (1): 24-34. Yekta, A., et al. (2018). "Determination of HydrogenβWater Relative Permeability and Capillary Pressure in Sandstone: Application to Underground Hydrogen Injection in Sedimentary Formations." Transport in Porous Media: 1-24. Zhao, B., et al. (2016). "Wettability control on multiphase flow in patterned microfluidics." Proceedings of the National Academy of Sciences (37): 10251-10256. Zhou, A., et al. (2016). "Capillary water retention curve and shear strength of unsaturated soils." Canadian Geotechnical Journal53