- Email: [email protected]
Prediction of capillary pressure-saturation relationship for primary drainage in a 3D fibrous porous medium using volume-of-fluid method
European Journal of Mechanics / B Fluids 67 (2018) 357–365
Contents lists available at ScienceDirect
European Journal of Mechanics / B Fluids journal homepage: www.elsevier.com/locate/ejmflu
Prediction of capillary pressure-saturation relationship for primary drainage in a 3D fibrous porous medium using volume-of-fluid method Nikhil Kumar Palakurthi a, *, Santosh Konangi a , Aravind Kishore a , Ken Comer b , Urmila Ghia a a b
Department of Mechanical and Materials Engineering, University of Cincinnati, Cincinnati, OH 45221-0072, USA The Procter and Gamble Company, Cincinnati, OH 45201, USA
article
info
Article history: Received 9 January 2017 Received in revised form 26 September 2017 Accepted 9 October 2017
Keywords: Fibrous media Volume-of-fluid method Full-morphology method Capillary-pressure saturation Microscale simulation Quasi-static drainage
a b s t r a c t The unsaturated flow through fibrous porous media at the macroscale is typically described using the Richards equation which requires constitutive relations for capillary pressure and relative permeability as a function of liquid saturation. In literature, these constitutive relations are typically estimated using reduced-order modeling approaches such as pore-network (PN) models, full-morphology (FM) method etc. In this paper, we determine quasi-static capillary pressure-saturation relationship (Pc -S) for primary drainage in a 3D isotropic fibrous medium by performing direct numerical simulations (at microscale) with the volume-of-fluid (VOF) method, which takes into account the detailed description of the pore structure and interface dynamics. As a first step, the accuracy of the VOF method was verified by simulating a test case of quasi-static drainage between two cylinders. Next, a grid-convergence study was performed to determine the optimal grid resolution required for the 3D simulations of quasi-static drainage. Using this grid, the Pc -S relationship for primary drainage in an isotropic fibrous structure was calculated from the direct simulations, and the results were compared with the Pc -S estimated using fullmorphology method, which is a quasi-static reduced-order geometric approach. In the VOF method, no modeling assumptions are made on the geometry of the pore structure and the liquid propagation, and thus, the results from direct simulations can be used to gain insights into the limitations of the reducedorder models. © 2017 Elsevier Masson SAS. All rights reserved.
1. Introduction Fibrous porous materials are utilized in a variety of applications such as sanitary products, performance apparel and drug delivery patches, etc., owing to their softness, high mechanical strength and unique properties in absorbing, storing, and releasing fluids [1,2]. As engineered fibrous materials become increasingly sophisticated, a fundamental understanding of the physical processes involved in unsaturated flow through fibrous materials is imperative for successful product design and development [3,4]. Traditionally, unsaturated flow through porous media at the macroscale is described using Richards equation [5]. Richards equation for unsaturated flow:
*
Correspondence to: 598 Rhodes Hall, PO Box 210072, Cincinnati, OH 452210072, USA. E-mail address: [email protected] (N.K. Palakurthi). https://doi.org/10.1016/j.euromechflu.2017.10.008 0997-7546/© 2017 Elsevier Masson SAS. All rights reserved.
For two-phase flow through porous media, modified Darcy’s law (Buckingham–Darcy law) in the absence of gravity is formulated as
( q=−
kkr (S)
µ
)
∇ Pc (S).
(1)
Here, q is the specific discharge or volumetric flow rate per unit area perpendicular to flow direction, k is the intrinsic permeability, kr is the relative permeability (a dimensionless measure of the effective permeability of that phase), µ is the dynamic viscosity of the fluid, Pc is capillary pressure which is the pressure difference across the interface between two fluids, and S is the wettingphase saturation. Combining the modified Darcy equation with the continuity equation (mass balance) results in the well-known Richards equation [5], which can be solved to obtain medium’s saturation as a function of time and space,
φ∂ S dt
[ = −∇.
kkr (S)
µ
] ∇ Pc
(2)
358
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Here, φ is the porosity of the medium. It is important to note that both relative permeability and capillary pressure are dependent on liquid saturation (S) of a porous medium. In order to predict liquid transport at the macroscale, Richards equation requires two constitutive relations: capillary pressure (Pc ) and relative permeability (kr ) as a function of liquid saturation (S) [4]. These constitutive relations are usually determined from capillary desaturation experiments, which are cumbersome to perform. The desaturation experiment starts with a porous medium that is fully saturated with the wetting phase. The wetting phase is drained out of porous medium in a quasi-static manner by either injecting a non-wetting phase or by incrementally increasing the pressure of the non-wetting phase reservoir. This process is referred to as primary drainage. Following primary drainage, the direction of saturation change is reversed wherein the wettingphase is imbibed into the system by displacing the non-wetting phase — this process is referred to as primary imbibition. Subsequently, above drainage-imbibition process is repeated to determine capillary pressure-saturation relationship during secondary drainage, secondary imbibition followed by the intermediate scanning curves. Pore-scale Modeling: In order to estimate the constitutive relations needed for the macroscale models, various microscale modeling approaches such as pore-network (PN) models, full-morphology (FM), lattice Boltzmann (LB), and smoothed particle hydrodynamics (SPH) methods have been used in the last two decades [6–13]. These numerical approaches differ in the amount of input data required, effort of data processing, and prediction potential [12]. Among these microscale modeling methods, PN modeling is the preferred method in literature as it is computationally less expensive, and thereby allows for modeling of larger domain sizes than possible with the other microscale approaches. However, in the PN method, the pore-space is idealized using primitive geometries and the fluid flow propagation is governed by simple set of rules [12]. Unlike the PN method, the FM method allows for a detailed representation of the porous geometry, and morphological openings are used to determine the flow (or open) paths inside the pore-space [14]. One disadvantage of the FM method is that it neglects the complex interplay between the various forces encountered during dynamic liquid propagation processes in the porous media [6,11,12,14]. In recent times, rapid increase in computing power and the development of robust numerical algorithms have enabled detailed numerical simulations of fluid flow through porous micro-structures using SPH and LB methods. The SPH and LB methods allow for sub-pore resolution and honor conservation equations in the appropriate limit [15–17]. However, the surface tension and contact angle are commonly implemented as special interaction forces acting between different particles and lattice nodes in SPH and LB methods, respectively and require calibration in some cases [18]. Also, porous structures like fibrous media pose an additional challenge to SPH and LBM methods, as capturing tortuous flow paths accurately would require using a large number of particles or a very fine grid. An alternative to the PN, FM, SPH and LB microscale approaches is to perform direct numerical simulations using the finite-volume based Volume-of-Fluid (VOF) method. The VOF method relies only on conservation principles and solves a system of equations to model the flow field evolution of two or more immiscible fluids [18–22]. Unlike the other microscale modeling approaches, VOF method allows for explicit description of interface dynamics by directly including the effects of surface tension at the interface, and the wall adhesion between the wetting liquid and solid fibers (in the form of contact angle) [18,19]. Since no modeling assumptions are made on the geometry of the pore structure and flow propagation, direct numerical simulations with VOF method are computationally challenging but will provide valuable insight into
the multi-phase flow dynamics at the porescale. The results from direct simulations can be used to understand the limitations of the reduced-order models (such as pore-network models, etc.), and for their calibration; verified reduced-order models can be used to study physical problems with larger length and time scales. To the best of authors’ knowledge, no study has been reported using the VOF method to determine capillary pressure-saturation (Pc -S) relationship in a fibrous porous medium. In this paper, we present the quasi-static Pc -S relationship for primary drainage in a 3D virtual isotropic fibrous structure determined using the finite-volume based VOF method. Two-phase flow simulations at the microscale were performed using the VOF method employed in OpenFOAM, an open-source CFD code (OpenCFD Ltd., v 2.0, www.openfoam.org). The accuracy of the two-phase calculations was established by simulating a test case of quasi-static drainage between two cylinders, and the numerical results were compared with the analytical solution published in the literature. The quasi-static Pc -S relationship for primary drainage in an isotropic fibrous structure was determined from the direct numerical simulation with VOF method. The results from VOF method were compared with the Pc -S relationship estimated using FM method (a reduced-order approach). 2. Methodology The first part of this section presents the details of the construction of a 3D isotropic fibrous structure, following which we describe the conservation equations of the VOF method, and the boundary conditions used in the two-phase flow simulations. Finally, the computational grid generation is presented. 2.1. Generation of 3D virtual isotropic fibrous structure A 3D virtual isotropic fibrous structure with a domain size of 2 mm ×2 mm × 1 mm was constructed using the commercial software GeoDict, a voxel-based code (Math2Market GmbH, GeoDict2012R1, www.geodict.com). The fiber diameter and porosity were kept constant at 50 µm and 0.8, respectively. In GeoDict, the fiber orientation distribution was controlled by a density function p(ψ ) defined in polar coordinates as, p(ψ ) =
1 4π
(
β sin ψ (1 + (β 2 − 1)cos2 ψ )3/2
) (3)
where, β is the anisotropy parameter which controls the fiber orientation of the micro-structure generation, and ψ ∈ [0, π ) is the through-plane angle [14]. Due to assumed isotropy in the xyplane, the density function in Eq. (3) depends on through-plane angle only [14], and the fiber orientation is controlled by the value of β ; for an isotropic structure, β is equal to unity. Fig. 1 shows the isotropic structure constructed in this study. The virtual isotropic structure constructed in GeoDict software was exported as a stereo-lithography (STL) file to create computational grids (described in Section 2.3) needed for the two-phase flow simulations. 2.2. Two-phase flow simulations: Volume-of-Fluid (VOF) method Transient 3D simulation of quasi-static drainage of a wetting liquid through fibrous structure was conducted using the VOF method. The VOF method is a fluid fraction-based interfacecapturing technique used for modeling the flow of a system of two or more incompressible, immiscible fluids [19,23,24]. The method is based on a whole-domain formulation that treats a two-fluid system as a single-fluid system with space-dependent physical properties. In the conventional VOF method, the continuity equation (Eq. (4)), momentum conservation equations (Eq. (5)), and a
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Fig. 1. Isometric view of the 3D virtual isotropic fibrous structure considered in this study.
transport equation for the fluid fraction of one phase (Eq. (6)) are solved simultaneously. The governing equations are as follows: Continuity equation (mass conservation) : ∇ · U = 0,
∂ (ρ U) + ∇ · (ρ UU) Momentum equation : ∂t T = −∇ P + ∇ · µ(∇ U + ∇ U ) + fσ , Fluid-fraction advection equation :
∂F + ∇ · (UF ) = 0. ∂t
(4)
(6)
Here, U is the velocity vector, F is the fluid-fraction, fσ is the surface tension force, ρ and µ are density and dynamic viscosity of the fluid, respectively. The fluid-fraction (F ) lies between 0 and 1 (0 ≤ F ≤ 1); F = 0 corresponds to cells filled with the non-wetting phase, whereas F = 1 indicates cells containing the wetting phase. Gradients of the fluid fraction are encountered only in the region of the interface between the two phases. The physical properties of the fluids are calculated as weighted averages based on the distribution of the fluid fraction function,
µ = µw F + µnw (1 − F ),
(7)
ρ = ρw F + ρnw (1 − F ).
(8)
Here, µw and µnw are the viscosities of the wetting and nonwetting phases, respectively, and ρw and ρnw are the densities of the wetting and non-wetting phases, respectively. The surface tension force (fσ ) at the fluid–fluid interface in the momentum conservation equation (Eq. (5)) is computed using the continuum surface force (CSF) method [25]. In the CSF method [21,22], the surface tension force acting at the interface is represented as a body force or volume force over the region occupied by the interface (0 < F < 1), and is formulated as fσ = σ κ∇ F ,
(9)
Here, σ is the surface tension between the wetting and nonwetting phases, κ is the mean curvature of the free surface which represents the magnitude of the interface normal flux at a specific face of the cell, and is formulated as
( ) κ = ∇ · nˆ c
Fig. 2. Saturated fibrous structure connected to imaginary non-wetting and wetting phase reservoirs located at the top and bottom boundaries, respectively.
In the VOF method, the interaction between the fluids and the solid surface (wall adhesion) is included as a boundary condition on the surface normal vector at the solid wall. nˆ c = nˆ s cos θ + tˆs sin θ
(5)
(10)
Here, nˆ c is the unit vector normal to the fluid interface calculated from the fluid fraction field as (∇ F ) nˆ c = (11) (|∇ F |)
359
(12)
Here, nˆ s and tˆs are the unit vectors normal and tangential to the solid boundary, respectively; θ is the equilibrium or static contact angle. The surface normal vector is corrected at the contact line to satisfy the specified static contact angle boundary condition following which the interface curvature (Eq. (10)) is calculated using this corrected surface normal vector. In OpenFOAM, a modified fluid fraction advection equation based on a two-fluid approach is solved to ensure boundedness and conservation of the fluid fraction field (F ), which is critical for two-phase flows with high density ratios [18,20,23]. Readers are referred to Deshpande et al. [20] for additional details of the modified fluid fraction advection equation implemented in OpenFOAM. Initial and boundary conditions of quasi-static drainage simulation: The quasi-static drainage computations were performed with water (ρ = 1000 kg/m3 and µ = 1e-3 Pa.s) and air (ρ = 1 kg/m3 and µ = 1.8e-5 Pa.s) as wetting and non-wetting phases, respectively. The surface tension between the wetting and non-wetting phases was specified as σ = 0.0725 N/m. A fixed or static contact angle of θ = 55 degrees was used; the dynamic variation of contact angle with contact line velocity was ignored. At time t = 0 s, the inter-fiber pore space is fully saturated with the wetting liquid (water). Fig. 2 shows a saturated isotropic structure connected to imaginary reservoirs of non-wetting and wetting phases on the top and bottom boundaries, respectively. A fixed gauge pressure of 0 Pa was specified at the outlet (wetting phase reservoir) boundary. The non-wetting fluid was forced into the pore-space by incrementally increasing the pressure at the inlet (non-wetting phase reservoir) in pre-determined steps, starting from 500 Pa. Symmetry boundary condition (normal velocity and shear forces set to zero) was applied on the lateral boundaries of the domain, and a no-slip boundary condition was imposed on the surface of the fibers. The effect of gravity is negligible at the porescale and was not included in our simulations. To perform a quasi-static drainage simulation, air flow (non-wetting phase) through the outlet boundary (wetting phase reservoir) must be prevented [26]. This is because when the non-wetting phase breaks through the wetting-phase reservoir boundary (during a primary drainage), a least resistant pathway for non-wetting phase is formed. At this point, any further increase in non-wetting phase reservoir pressure forces the non-wetting phase to flow through the least resistant path, and no additional Pc S data can be determined. To avoid this, in quasi-static desaturation
360
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Fig. 3. (a) Isometric view of a fibrous structure with a slice at the middle of y-axis (ymid ), and (b) Computational grid corresponding to 10 µm resolution at ymid .
experiment, a frit (with pore-size smaller than the porous medium of interest) is used which makes it possible to measure full the Pc S curve. In our numerical study, we employed a similar approach and used a virtual frit (array of parallel cylinders of diameter 20 µm separated by a distance of 10 µm, as shown in Fig. 3a) at the outlet (wetting phase reservoir) of our computational domain. The back pressure created by the virtual frit prevented the non-wetting phase (air) from escaping the outlet (wetting-phase reservoir) in our simulation. 2.3. Computational grid generation One of the main challenges with porescale simulations using conventional CFD methods is the generation of a computational grid which adequately resolves the tortuous inter-fiber void spaces, and offers the feasibility to perform two-phase flow simulations. In this paper, hexahedral-dominant computational grids were generated for the 3D virtual isotropic fibrous structure using foamProMesh (iconCFD, www.iconcfd.com). The tortuous fluid paths inside the fibrous structure necessitated the use of unstructured grids. As a first step in the meshing procedure, a background hexahedral mesh was created using the desired resolution maintaining a cell-aspect ratio of unity. Subsequently, the unwanted cells inside the solid region (of the fibers) were removed, and a boundary-fitted grid was obtained by iteratively refining and snapping the cells to the fiber surfaces. Fig. 3b shows the computational grid corresponding to 10 µm resolution at a slice in the middle of the y-axis (shown in Fig. 3a). The results of gridconvergence study are presented in Section 3.2. 3. Results and discussion In this section, we first present the results of the verification case (quasi-static drainage between two cylinders) which was studied to check the accuracy of the two-phase flow calculations. Next, we describe the results of our grid convergence study for the microscale simulations. Subsequently, we present the results of quasi-static drainage simulation performed using VOF method in a 3D isotropic fibrous structure. Finally, the Pc -S relationship determined from VOF method is compared with the Pc -S relationship estimated using FM method, and discussed in detail. 3.1. Quasi-static drainage of wetting liquid between two horizontal parallel cylinders: Verification of two-phase flow calculations using VOF method In literature, Princen [27] studied the equilibrium configuration of limited amounts of liquid in horizontal assemblies of parallel cylinders, and proposed analytical solutions for radius of curvature and capillary pressure. We adopted Princen’s work [27] and modified the problem setup to simulate quasi-static drainage of a
wetting liquid between two cylinders. A schematic of the cross section of a liquid column between two infinite horizontal cylinders of radii (r), separated by a distance 2d, is shown in Fig. 4. For this configuration, the radius of curvature (Eq. (13)) of meniscus and capillary pressure (Eq. (14)) were formulated by Princen [27] as
[ R= Pc =
d + r − r cos α
] (13)
cos(α + θ )
σ R
=
σ r
[
cos(α + θ ) d/r + 1 − cos α
] (14)
Here, σ is the surface tension between wetting and non-wetting phases, θ is the contact angle made by wetting liquid on solid cylinders, α is the angle between the line connecting the centers of the cylinders and the radius to the solid–liquid–air interface, R is the curvature of the interface, and Pc is the difference between non-wetting and wetting phase reservoir pressures. In the present study, the radii of the cylinders (r) and the separation distance (2d) were fixed at 1 mm and 0.2 mm, respectively, and a static contact angle (θ ) of 55 degrees was used. The quasi-static drainage was simulated for air–water configuration (properties presented in Section 2.2). For a capillary pressure of Pc = 0 Pa (starting point of the drainage process), alpha angle is α = 35 degrees. For a given geometric configuration (fixed d/r ratio) of parallel cylinders and contact angle (θ ), there exists a critical value for alpha αcr (and thus a critical capillary pressure) beyond which a constant contact angle cannot be maintained; at this point, capillary pressure at the interface will not be sufficient to oppose any disturbances created at the interface resulting in complete drainage of the wetting liquid. Critical alpha (αcr ) can be determined by differentiating Eq. (14) with respect to alpha (α ) and c = 0), resulting in equating it to zero ( dP dα
αcr = sin−1
(
sin θ d/ r + 1
)
−θ
(15)
The critical alpha (from Eq. (15)) and critical capillary pressure (from Eq. (14)) values calculated are αcr = −6.87 degrees and 448.36 Pa, respectively. Hence, as capillary pressure increases from 0 Pa to 448.36 Pa, alpha angle (α ) decreases from 35 degrees to −6.87 degrees, at which point the wetting phase will be completely drained out of the system. Any value of capillary pressure greater than 448.36 Pa will instantly drain out the wetting phase between two cylinders. To evaluate the accuracy of the two-phase flow calculations and contact angle model implemented in OpenFOAM, we simulated the drainage of wetting phase between two cylinders in a quasistatic manner. Fig. 4b shows the computational domain for the verification problem with boundary conditions used in the numerical simulation. Since the domain is symmetric, computations were carried out only on one half of the domain, and symmetry BC was applied at the mid-plane between the cylinders. No-slip boundary condition was imposed on the cylinder surface. The gauge pressure
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
361
Fig. 4. Liquid drainage between two horizontal parallel cylinders. (a) Schematic of cross section of two infinite horizontal cylinders of radii (r) separated by a distance 2d, (b) Computational domain with boundary conditions.
at the lower boundary (water reservoir) was fixed at 0 Pa. A grid convergence study was performed for this verification case using four grid resolutions (d/6, d/12, d/18, d/24; d = 0.1 mm) to ensure the results are grid-independent (not included in the manuscript); the difference between numerical results and analytical solution was found to be less than 2%. From the grid-convergence study, a computational grid with a grid-resolution of d/12 (8.3 µm) was used for the quasi-static drainage simulation. Quasi-static drainage of wetting liquid was simulated by incrementally increasing the pressure of the top boundary (air reservoir) from 0 Pa to 448.36 Pa (critical Pc for the configuration considered) with a relaxation phase between successive pressure steps to allow the system to reach equilibrium as shown in Fig. 5a. From the numerical results (wetting-phase saturation), alpha angle was estimated, and plotted as a function of time in Fig. 5b. It can be observed that, for each pressure step, the alpha angle decreases and then reached an equilibrium value (indicated by gray circles). Fig. 5c shows the numerical (equilibrium values from Fig. 5b) and analytical values of alpha for a range of capillary pressures (Fig. 5a) considered in this study; the analytical values of alpha angle for a given capillary pressure or vice-versa can be determined from Eq. (14). It can be seen that for all the Pc values, the alpha values determined from the simulation agree well with the analytical values, except near critical alpha (αcr ), a point at which wetting liquid completely drains out of the system. In the numerical simulation, the wetting phase drainage through the cylinders occurred slightly before the analytical prediction of critical capillary of 448.36 Pa. This is because, the analytical solution is valid only at the equilibrium, and do not take into account the dynamic aspects of the interface propagation. In a quasi-static simulation performed using the VOF method, each pressure step is akin to a dynamic simulation (with applied pressure) during which wetting liquid drains out until the applied pressure and inertial effects (dynamic pressure) are balanced by the capillary pressure at the interface (contact angle) and viscous forces on the cylinder; at equilibrium (end of each pressure step), the inertial forces become negligible. Near the critical alpha angle, the dynamic effects (captured in simulations) resulted in an early drainage of wetting liquid through the cylinders. The good agreement between numerical and analytical results establishes confidence in the predictive capabilities and accuracy of the VOF method used in this study. In the next section, we present the results of the grid-convergence study of microscale simulations of quasi-static drainage in 3D isotropic fibrous structure.
3.2. Grid-convergence study: Drainage simulation using volume-offluid (VOF) method Since the two-phase computations using VOF method are sensitive to the grid resolution [18–20], a grid convergence study was performed to ensure that the results were grid-independent. In order to estimate the effects of discretization, drainage of wetting liquid through an isotropic fibrous structure was simulated using computational grids with different resolutions (50 µm, 25 µm, 16.33 µm, 12.5 µm, 10 µm, and 8 µm). A cell-aspect ratio close to one, which is desired for stable two-phase computations, was used. A grid refinement factor of at least 1.25 was maintained between successive grids. An isometric scaling (N ∼ ∆x−3 ) was observed between the number of cells (N) and grid resolution (∆x), indicating uniform internal refinement. For these simulations, the fluid was initialized (t = 0 s) such that inter-fiber space was completely saturated with the wetting liquid. For this case, a constant pressure differential was applied in the simulation by specifying gauge pressure values of 900 Pa and 0 Pa at the inlet (non-wetting phase reservoir) and outlet (wetting phase reservoir), respectively. Fig. 6 shows the variation in wetting phase saturation as a function of time for different grid resolutions. It was observed that the liquid saturation decreases and reaches an equilibrium value. From Fig. 6, it can be noticed that the variation in saturation levels was found to be minimal between finest grid (8 µm), finer or penultimate grids (10 µm), and fine grid (12.5 µm). At t = 0.01 s, the difference between the saturation levels predicted by finest grid (8 µm) and finer or penultimate grid (10 µm) was 0.009%, and the difference between finest grid (8 µm) and fine grid (12.5 µm) was 1.54%. Based on these results, the finer grid resolution (10 µm) was deemed sufficient, and was used for the quasi-static drainage numerical simulation. 3.3. Quasi-static drainage numerical simulation 3.3.1. Volume-of-fluid (VOF) method At the beginning (t = 0 s) of the quasi-static drainage simulation, the inter-fiber space was fully saturated with the wetting liquid (water). The pressure of the wetting phase reservoir was fixed at 0 Pa during the simulation. As described in Section 2.2, the quasi-static drainage was simulated by increasing the pressure of the non-wetting phase reservoir incrementally. Fig. 7a shows the pressure applied at the inlet (non-wetting phase boundary) or pressure differential between non-wetting and wetting phase reservoirs as a function of time. The wetting phase saturation in
362
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Fig. 5. (a) Pressure applied at the top or non-wetting phase boundary as a function of time, (b) Transient variation of alpha angle in simulations; equilibrium alpha values determined for each pressure step are indicated by gray circles, and (c) Comparison of numerically determined alpha values with analytical alpha values for a range of capillary pressure values considered in this study.
Fig. 8 shows the fluid fraction contours at a slice in the middle of y-plane (Fig. 3a) for different wetting phase saturations (a) S = 1, (b) S = 0.83, (c) S = 0.47, and (d) S = 0.28. It can be observed that, with the increase in capillary pressure, the non-wetting phase intrudes into the complex inter-fiber space displacing the wetting phase. In the VOF method, no modeling assumptions are made on the liquid drainage process, and at each pressure step, the non-wetting phase intrusion process is dictated by the smaller inter-fiber gaps, pore connectivity, and complex interplay between applied pressure differential, inertial, viscous and capillary forces. The P-c -S relationship determined from the VOF method was compared with the results from the FM method (implemented in GeoDict software), a commonly used quasi-static reduced-order geometric approach explained in the next Section 3.3.2. Fig. 6. Effect of grid resolution on wetting phase saturation in drainage simulation performed using volume-of-fluid (VOF) method; transient variation of wetting phase saturation was shown for six different grid resolutions ∆x = 50 µm, 25 µm, 16.33 µm, 12.5 µm, 10 µm, 8 µm.
the inter-fiber space was monitored during the simulations and is shown in Fig. 7b. For each pressure increment, the wetting phase saturation decreases initially, and then reaches an equilibrium value (indicated by gray circles in Fig. 7b). The equilibrium values of wetting phase saturation obtained at the end of each pressure step were used to plot the quasi-static Pc -S relationship. We present the capillary pressure saturation relationship (Pc -S) determined from the VOF method in Section 3.4.
3.3.2. Full-morphology method The full-morphology (FM) method implemented in GeoDict software was used to simulate the quasi-static drainage in a 3D isotropic fibrous structure for air–water configuration (properties presented in Section 2.2). In FM method, the non-wetting fluid was forced into the saturated inter-fiber space by incrementally increasing the pressure of the non-wetting phase reservoir. It is assumed that only certain pore radii (r), as described by the Young– Laplace equation (Eq. (16)), can be intruded by the non-wetting phase at a given capillary pressure (Pc ), Pc =
2σ cosθ r
(16)
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
363
Fig. 7. (a) Pressure applied at the inlet or non-wetting phase boundary as a function of time; a gauge pressure of 0 Pa was applied at the outlet or wetting phase boundary. (b) Transient variation of wetting saturation inside the isotropic fiber structure.
Fig. 8. Fluid-fraction contours for different wetting phase saturations on a slice at the middle of y-axis of the model; (a) S = 1, (b) S = 0.83, (c) S = 0.47, and (d) S = 0.28.
Briefly, the two steps in the FM method are as follows: Erosion step: For each pressure increment, a spherical structural element of radius r, determined using the Young–Laplace equation at a given Pc , was used in a morphological operation of the porespace. This erosion operation renders all possible center points in a domain where a given sphere could be fitted without intruding the fiber boundaries. The erosion operation was performed without any consideration to the connectivity of the fluids to their respective reservoirs [14]. In order to determine Pc -S relationship, a connectivity analysis was performed to remove all those eroded spaces which are not connected to the fluid reservoirs. Dilation step: Following the connectivity analysis, the pore volume occupied by the wetting phase was determined by dilating the structuring element, and the above procedure was repeated by incrementing pressure in predetermined steps [14,28]. Fig. 9 shows intermediate stages of the quasi-static drainage simulation: (a) Pc = 475 Pa, S = 1, (b) Pc = 1232 Pa, S = 0.62, (c) Pc = 1584 Pa, S = 0.3, and (d) Pc = 3696 Pa, S = 0.04. In Fig. 9, the wetting phase (water) is shown in dark gray, the fibers are show in light gray and the non-wetting phase is not colored. It can be seen that the non-wetting phase intrudes the inter-fiber space with the increase in capillary pressure, thereby displacing the wetting phase (water) from the structure. The results of the quasi-static Pc -S relationship
determined from the FM method are presented along with the results from VOF method in the next section. 3.4. Comparison of capillary pressure-saturation (Pc -S) relationship determined using volume-of-fluid (VOF) method and full-morphology (FM) method Fig. 10 shows the quasi-static Pc -S relationship for primary drainage of a 3D isotropic fibrous medium determined using the VOF method along with the results from the FM method, which is a simplified porescale technique. It can be seen that, for 0.15 < S < 0.8, a reasonably good agreement was observed between the Pc -S relationships predicted by the VOF and FM methods. However, for S > 0.9, the FM method predicted high entry pressures as compared to the VOF method (Fig. 10); this is because the VOF method takes into account the dynamic aspects of liquid propagation whereas FM method does not account for the inertial effects which could assist in the drainage of wetting phase. Also, for Pc > 2500 Pa, the FM method predicted slightly lower residual wetting phase saturation as compared to the VOF method (Fig. 10). This could be attributed to the lack of contact line physics (and surface tension force of the interface) in the FM method. At lower saturations, this becomes especially important as liquid blobs could be trapped in small pore spaces, and their effective
364
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Fig. 9. Intermediate stages of quasi-static drainage simulation performed using full-morphology method: (a) P c = 475 Pa, S = 1, (b) P c = 1232 Pa, S = 0.62, (c) P c = 1584 Pa, S = 0.3, and (d) P c = 3696 Pa, S = 0.04. The solid fibers are shown in light gray color, the wetting phase (water) is shown in dark gray color.
capillary pressure could be greater than predicted by the Young– Laplace equation (with a spherical structural element of radius r) used in FM method. Even though the FM method uses simpler set of rules (morphological openings) for liquid propagation in the simulations, the Pc -S determined from the FM method compared reasonably well with the results from direct simulations performed using VOF method (in which no modeling assumptions are made on the liquid propagation). The results indicate that the simplified equations in the FM method reasonably represent multi-phase flow dynamics of quasi-static primary drainage process. In this fashion, the numerical results from direct simulations such as VOF method can be used to verify the accuracy of the modeling assumptions in the reduced-order porescale techniques.
pressure-saturation (Pc -S) relationship for primary drainage in a 3D isotropic fibrous medium. The accuracy of the two-phase VOF calculations was verified by simulating a test case of quasi-static drainage between two infinite parallel cylinders; the numerical results agreed well with the analytical solution. Next, a grid convergence study was performed for quasi-static drainage process in a 3D isotropic fibrous structure, and an optimal grid needed for the microscale simulation was determined. Using this grid, the quasistatic drainage simulation was conducted, and Pc -S relationship was calculated. The Pc -S relationship determined from the direct simulations using the VOF method was compared with full-morphology (FM) method, a popularly used reduced-order approach [4,14]. Good agreement was observed between the results of VOF and FM methods for wetting phase saturation levels of S = 0.15–0.8. The results show that the FM method which relies on quasi-static geometric approach, a relatively simpler physics as compared to VOF method, may be sufficient for estimating Pc -S relationship for a quasi-static primary drainage process in a fibrous medium. Unlike reduced-order modeling approaches (e.g. pore-network models), the VOF method allows for a detailed description of the porous structure, and honors fluid flow conservation equations. The liquid propagation in the VOF method is governed by complex interplay between applied pressure differential, viscous and capillary forces at the microscale. Since no modeling assumptions are made on the geometry of the pore structure and liquid transport, the direct simulations with VOF method are computationally challenging to perform. However, these direct simulations (VOF, SPH, LBM etc.) provide reliable results that can be used for verification and calibration of reduced-order models, which can be used to study physical problems over large length and time scales. The direct simulations can also be used to gain valuable insights into the effect of various porescale processes on macroscale properties under both quasi-static and dynamic conditions.
4. Conclusions
Acknowledgments
In the present study, we performed direct numerical simulations at the microscale using the finite-volume based volumeof-fluid (VOF) method to determine the quasi-static capillary
This work is sponsored by The Procter & Gamble Company through the University of Cincinnati Simulation Center. The computations were performed on a parallel computer cluster at the
Fig. 10. Comparison of capillary pressure-saturation relationship for primary drainage in a 3D isotropic fibrous structure determined using the volume-of-fluid (VOF) and full-morphology (FM) methods.
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Procter and Gamble Company consisting of 2.93 GHz Intel Xeon processors (X5670) with 2 Gbytes of memory per processor. The authors thank Dr. Thomas Baer from the Procter & Gamble Company for his valuable suggestions. The authors also thank the anonymous referees for their insightful comments and suggestions which have significantly improved the manuscript. References [1] A. Ashari, T. Bucher, H. Vahedi Tafreshi, A semi-analytical model for simulating fluid transport in multi-layered fibrous sheets made up of solid and porous fibers, Comput. Mater. Sci. 50 (2) (2010) 378–390. [2] C.J. Hong, J.B. Kim, A study of comfort performance in cotton and polyester blended fabrics. I. Vertical wicking behavior, Fibers Polym. 8 (2) (2007) 218– 224. [3] M. Amiot, M. Lewandowski, P. Leite, M. Thomas, A. Perwuelz, An evaluation of fiber orientation and organization in nonwoven fabrics by tensile, air permeability and compression measurements, J. Mater. Sci. 49 (1) (2014) 52–61. [4] A. Ashari, T. Bucher, H.V. Tafreshi, M. Tahir, M. Rahman, Modeling fluid spread in thin fibrous sheets effects of fiber orientation, Int. J. Heat Mass Transfer 53 (9) (2010) 1750–1758. [5] L.A. Richards, Capillary conduction of liquids through porous mediums, Physics 1 (5) (1931) 318–333. [6] M. Hilpert, C.T. Miller, Pore-morphology-based simulation of drainage in totally wetting porous media, Adv. Water Resour. 24 (3) (2001) 243–255. [7] V. Joekar-Niasar, S.M. Hassanizadeh, H. Dahle, Non-equilibrium effects in capillarity and interfacial area in two-phase flow: dynamic pore-network modelling, J. Fluid Mech. 655 (1) (2010) 38–71. [8] V. Joekar-Niasar, M. Prodanović, D. Wildenschild, S. Hassanizadeh, Network model investigation of interfacial area, capillary pressure and saturation relationships in granular porous media, Water Resour. Res. 46 (6) (2010). [9] P.P. Mukherjee, Q. Kang, C.-Y. Wang, Pore-scale modeling of two-phase transport in polymer electrolyte fuel cells—progress and perspective, Energy Environ. Sci. 4 (2) (2011) 346–369. [10] M.L. Porter, M.G. Schaap, D. Wildenschild, Lattice-Boltzmann simulations of the capillary pressure–saturation–interfacial area relationship for porous media, Adv. Water Resour. 32 (11) (2009) 1632–1640. [11] V.P. Schulz, J. Becker, A. Wiegmann, P.P. Mukherjee, C.-Y. Wang, Modeling of two-phase behavior in the gas diffusion medium of PEFCs via full morphology approach, J. Electrochem. Soc. 154 (4) (2007) B419–B426. [12] H.-J. Vogel, J. Tölke, V. Schulz, M. Krafczyk, K. Roth, Comparison of a latticeBoltzmann model a full-morphology model and a pore network model for determining capillary pressure–saturation relationships, Vadose Zone J. 4 (2) (2005) 380–388.
365
[13] M.S. Riasi, N.K. Palakurthi, C. Montemagno, L. Yeghiazarian, A feasibility study of the pore topology method (PTM), A medial surface-based approach to multiphase flow simulation in porous media, Transp. Porous Media 1 (2016) 1–21. [14] A. Ashari, H. Vahedi Tafreshi, A two-scale modeling of motion-induced fluid release from thin fibrous porous media, Chem. Eng. Sci. 64 (9) (2009) 2067– 2075. [15] A. Nabovati, E.W. Llewellin, A. Sousa, A general model for the permeability of fibrous porous media based on fluid flow simulations using the lattice Boltzmann method, Composites 40 (6) (2009) 860–869. [16] A.M. Tartakovsky, P. Meakin, Pore scale modeling of immiscible and miscible fluid flows using smoothed particle hydrodynamics, Adv. Water Resour. 29 (10) (2006) 1464–1478. [17] Y. Zhu, P.J. Fox, J.P. Morris, A pore-scale numerical model for flow through porous media, Int. J. Numer. Anal. Methods Geomech. 23 (9) (1999) 881–904. [18] A. Ferrari, I. Lunati, Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy, Adv. Water Resour. 57 (2013) 19– 31. [19] A. Ashish Saha, S.K. Mitra, Effect of dynamic contact angle in a volume of fluid (VOF) model for a microfluidic capillary flow, J. Colloid Interface Sci. 339 (2) (2009) 461–480. [20] S.S. Deshpande, L. Anumolu, M.F. Trujillo, Evaluating the performance of the two-phase flow solver interfoam, Comput. Sci. Discov. 5 (1) (2012) 014016. [21] D.A. Hoang, V. van Steijn, L.M. Portela, M.T. Kreutzer, C.R. Kleijn, Benchmark numerical simulations of segmented two-phase flows in microchannels using the volume of fluid method, Comput. & Fluids 86 (2013) 28–36. [22] A.Q. Raeini, M.J. Blunt, B. Bijeljic, Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method, J. Comput. Phys. 231 (17) (2012) 5653–5668. [23] E. Berberović, N.P. van Hinsberg, S. Jakirlić, I.V. Roisman, C. Tropea, Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution, Phys. Rev. E 79 (3) (2009) 036306. [24] N.K. Palakurthi, S. Konangi, U. Ghia, K. Comer, Micro-scale simulation of unidirectional capillary transport of wetting liquid through 3D fibrous porous media: Estimation of effective pore radii, Int. J. Multiph. Flow. 77 (2015) 48–57. [25] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface-tension, J. Comput. Phys. 100 (2) (1992) 335–354. [26] B. Ahrenholz, J. Tölke, P. Lehmann, A. Peters, A. Kaestner, M. Krafczyk, W. Durner, Prediction of capillary hysteresis in a porous material using latticeBoltzmann methods and comparison to experimental data and a morphological pore network model, Adv. Water Resour. 31 (9) (2008) 1151–1173. [27] H. Princen, Capillary phenomena in assemblies of parallel cylinders III. Liquid columns between horizontal parallel cylinders, J. Colloid Interface Sci. 34 (2) (1970) 171–184. [28] J. Becker, A. Wiegmann, V. Schulz, Design of fibrous filter media based on the simulation of pore size measures, in: Proceedings of the FILTECH 2007 conference, 2007.
Contents lists available at ScienceDirect
European Journal of Mechanics / B Fluids journal homepage: www.elsevier.com/locate/ejmflu
Prediction of capillary pressure-saturation relationship for primary drainage in a 3D fibrous porous medium using volume-of-fluid method Nikhil Kumar Palakurthi a, *, Santosh Konangi a , Aravind Kishore a , Ken Comer b , Urmila Ghia a a b
Department of Mechanical and Materials Engineering, University of Cincinnati, Cincinnati, OH 45221-0072, USA The Procter and Gamble Company, Cincinnati, OH 45201, USA
article
info
Article history: Received 9 January 2017 Received in revised form 26 September 2017 Accepted 9 October 2017
Keywords: Fibrous media Volume-of-fluid method Full-morphology method Capillary-pressure saturation Microscale simulation Quasi-static drainage
a b s t r a c t The unsaturated flow through fibrous porous media at the macroscale is typically described using the Richards equation which requires constitutive relations for capillary pressure and relative permeability as a function of liquid saturation. In literature, these constitutive relations are typically estimated using reduced-order modeling approaches such as pore-network (PN) models, full-morphology (FM) method etc. In this paper, we determine quasi-static capillary pressure-saturation relationship (Pc -S) for primary drainage in a 3D isotropic fibrous medium by performing direct numerical simulations (at microscale) with the volume-of-fluid (VOF) method, which takes into account the detailed description of the pore structure and interface dynamics. As a first step, the accuracy of the VOF method was verified by simulating a test case of quasi-static drainage between two cylinders. Next, a grid-convergence study was performed to determine the optimal grid resolution required for the 3D simulations of quasi-static drainage. Using this grid, the Pc -S relationship for primary drainage in an isotropic fibrous structure was calculated from the direct simulations, and the results were compared with the Pc -S estimated using fullmorphology method, which is a quasi-static reduced-order geometric approach. In the VOF method, no modeling assumptions are made on the geometry of the pore structure and the liquid propagation, and thus, the results from direct simulations can be used to gain insights into the limitations of the reducedorder models. © 2017 Elsevier Masson SAS. All rights reserved.
1. Introduction Fibrous porous materials are utilized in a variety of applications such as sanitary products, performance apparel and drug delivery patches, etc., owing to their softness, high mechanical strength and unique properties in absorbing, storing, and releasing fluids [1,2]. As engineered fibrous materials become increasingly sophisticated, a fundamental understanding of the physical processes involved in unsaturated flow through fibrous materials is imperative for successful product design and development [3,4]. Traditionally, unsaturated flow through porous media at the macroscale is described using Richards equation [5]. Richards equation for unsaturated flow:
*
Correspondence to: 598 Rhodes Hall, PO Box 210072, Cincinnati, OH 452210072, USA. E-mail address: [email protected] (N.K. Palakurthi). https://doi.org/10.1016/j.euromechflu.2017.10.008 0997-7546/© 2017 Elsevier Masson SAS. All rights reserved.
For two-phase flow through porous media, modified Darcy’s law (Buckingham–Darcy law) in the absence of gravity is formulated as
( q=−
kkr (S)
µ
)
∇ Pc (S).
(1)
Here, q is the specific discharge or volumetric flow rate per unit area perpendicular to flow direction, k is the intrinsic permeability, kr is the relative permeability (a dimensionless measure of the effective permeability of that phase), µ is the dynamic viscosity of the fluid, Pc is capillary pressure which is the pressure difference across the interface between two fluids, and S is the wettingphase saturation. Combining the modified Darcy equation with the continuity equation (mass balance) results in the well-known Richards equation [5], which can be solved to obtain medium’s saturation as a function of time and space,
φ∂ S dt
[ = −∇.
kkr (S)
µ
] ∇ Pc
(2)
358
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Here, φ is the porosity of the medium. It is important to note that both relative permeability and capillary pressure are dependent on liquid saturation (S) of a porous medium. In order to predict liquid transport at the macroscale, Richards equation requires two constitutive relations: capillary pressure (Pc ) and relative permeability (kr ) as a function of liquid saturation (S) [4]. These constitutive relations are usually determined from capillary desaturation experiments, which are cumbersome to perform. The desaturation experiment starts with a porous medium that is fully saturated with the wetting phase. The wetting phase is drained out of porous medium in a quasi-static manner by either injecting a non-wetting phase or by incrementally increasing the pressure of the non-wetting phase reservoir. This process is referred to as primary drainage. Following primary drainage, the direction of saturation change is reversed wherein the wettingphase is imbibed into the system by displacing the non-wetting phase — this process is referred to as primary imbibition. Subsequently, above drainage-imbibition process is repeated to determine capillary pressure-saturation relationship during secondary drainage, secondary imbibition followed by the intermediate scanning curves. Pore-scale Modeling: In order to estimate the constitutive relations needed for the macroscale models, various microscale modeling approaches such as pore-network (PN) models, full-morphology (FM), lattice Boltzmann (LB), and smoothed particle hydrodynamics (SPH) methods have been used in the last two decades [6–13]. These numerical approaches differ in the amount of input data required, effort of data processing, and prediction potential [12]. Among these microscale modeling methods, PN modeling is the preferred method in literature as it is computationally less expensive, and thereby allows for modeling of larger domain sizes than possible with the other microscale approaches. However, in the PN method, the pore-space is idealized using primitive geometries and the fluid flow propagation is governed by simple set of rules [12]. Unlike the PN method, the FM method allows for a detailed representation of the porous geometry, and morphological openings are used to determine the flow (or open) paths inside the pore-space [14]. One disadvantage of the FM method is that it neglects the complex interplay between the various forces encountered during dynamic liquid propagation processes in the porous media [6,11,12,14]. In recent times, rapid increase in computing power and the development of robust numerical algorithms have enabled detailed numerical simulations of fluid flow through porous micro-structures using SPH and LB methods. The SPH and LB methods allow for sub-pore resolution and honor conservation equations in the appropriate limit [15–17]. However, the surface tension and contact angle are commonly implemented as special interaction forces acting between different particles and lattice nodes in SPH and LB methods, respectively and require calibration in some cases [18]. Also, porous structures like fibrous media pose an additional challenge to SPH and LBM methods, as capturing tortuous flow paths accurately would require using a large number of particles or a very fine grid. An alternative to the PN, FM, SPH and LB microscale approaches is to perform direct numerical simulations using the finite-volume based Volume-of-Fluid (VOF) method. The VOF method relies only on conservation principles and solves a system of equations to model the flow field evolution of two or more immiscible fluids [18–22]. Unlike the other microscale modeling approaches, VOF method allows for explicit description of interface dynamics by directly including the effects of surface tension at the interface, and the wall adhesion between the wetting liquid and solid fibers (in the form of contact angle) [18,19]. Since no modeling assumptions are made on the geometry of the pore structure and flow propagation, direct numerical simulations with VOF method are computationally challenging but will provide valuable insight into
the multi-phase flow dynamics at the porescale. The results from direct simulations can be used to understand the limitations of the reduced-order models (such as pore-network models, etc.), and for their calibration; verified reduced-order models can be used to study physical problems with larger length and time scales. To the best of authors’ knowledge, no study has been reported using the VOF method to determine capillary pressure-saturation (Pc -S) relationship in a fibrous porous medium. In this paper, we present the quasi-static Pc -S relationship for primary drainage in a 3D virtual isotropic fibrous structure determined using the finite-volume based VOF method. Two-phase flow simulations at the microscale were performed using the VOF method employed in OpenFOAM, an open-source CFD code (OpenCFD Ltd., v 2.0, www.openfoam.org). The accuracy of the two-phase calculations was established by simulating a test case of quasi-static drainage between two cylinders, and the numerical results were compared with the analytical solution published in the literature. The quasi-static Pc -S relationship for primary drainage in an isotropic fibrous structure was determined from the direct numerical simulation with VOF method. The results from VOF method were compared with the Pc -S relationship estimated using FM method (a reduced-order approach). 2. Methodology The first part of this section presents the details of the construction of a 3D isotropic fibrous structure, following which we describe the conservation equations of the VOF method, and the boundary conditions used in the two-phase flow simulations. Finally, the computational grid generation is presented. 2.1. Generation of 3D virtual isotropic fibrous structure A 3D virtual isotropic fibrous structure with a domain size of 2 mm ×2 mm × 1 mm was constructed using the commercial software GeoDict, a voxel-based code (Math2Market GmbH, GeoDict2012R1, www.geodict.com). The fiber diameter and porosity were kept constant at 50 µm and 0.8, respectively. In GeoDict, the fiber orientation distribution was controlled by a density function p(ψ ) defined in polar coordinates as, p(ψ ) =
1 4π
(
β sin ψ (1 + (β 2 − 1)cos2 ψ )3/2
) (3)
where, β is the anisotropy parameter which controls the fiber orientation of the micro-structure generation, and ψ ∈ [0, π ) is the through-plane angle [14]. Due to assumed isotropy in the xyplane, the density function in Eq. (3) depends on through-plane angle only [14], and the fiber orientation is controlled by the value of β ; for an isotropic structure, β is equal to unity. Fig. 1 shows the isotropic structure constructed in this study. The virtual isotropic structure constructed in GeoDict software was exported as a stereo-lithography (STL) file to create computational grids (described in Section 2.3) needed for the two-phase flow simulations. 2.2. Two-phase flow simulations: Volume-of-Fluid (VOF) method Transient 3D simulation of quasi-static drainage of a wetting liquid through fibrous structure was conducted using the VOF method. The VOF method is a fluid fraction-based interfacecapturing technique used for modeling the flow of a system of two or more incompressible, immiscible fluids [19,23,24]. The method is based on a whole-domain formulation that treats a two-fluid system as a single-fluid system with space-dependent physical properties. In the conventional VOF method, the continuity equation (Eq. (4)), momentum conservation equations (Eq. (5)), and a
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Fig. 1. Isometric view of the 3D virtual isotropic fibrous structure considered in this study.
transport equation for the fluid fraction of one phase (Eq. (6)) are solved simultaneously. The governing equations are as follows: Continuity equation (mass conservation) : ∇ · U = 0,
∂ (ρ U) + ∇ · (ρ UU) Momentum equation : ∂t T = −∇ P + ∇ · µ(∇ U + ∇ U ) + fσ , Fluid-fraction advection equation :
∂F + ∇ · (UF ) = 0. ∂t
(4)
(6)
Here, U is the velocity vector, F is the fluid-fraction, fσ is the surface tension force, ρ and µ are density and dynamic viscosity of the fluid, respectively. The fluid-fraction (F ) lies between 0 and 1 (0 ≤ F ≤ 1); F = 0 corresponds to cells filled with the non-wetting phase, whereas F = 1 indicates cells containing the wetting phase. Gradients of the fluid fraction are encountered only in the region of the interface between the two phases. The physical properties of the fluids are calculated as weighted averages based on the distribution of the fluid fraction function,
µ = µw F + µnw (1 − F ),
(7)
ρ = ρw F + ρnw (1 − F ).
(8)
Here, µw and µnw are the viscosities of the wetting and nonwetting phases, respectively, and ρw and ρnw are the densities of the wetting and non-wetting phases, respectively. The surface tension force (fσ ) at the fluid–fluid interface in the momentum conservation equation (Eq. (5)) is computed using the continuum surface force (CSF) method [25]. In the CSF method [21,22], the surface tension force acting at the interface is represented as a body force or volume force over the region occupied by the interface (0 < F < 1), and is formulated as fσ = σ κ∇ F ,
(9)
Here, σ is the surface tension between the wetting and nonwetting phases, κ is the mean curvature of the free surface which represents the magnitude of the interface normal flux at a specific face of the cell, and is formulated as
( ) κ = ∇ · nˆ c
Fig. 2. Saturated fibrous structure connected to imaginary non-wetting and wetting phase reservoirs located at the top and bottom boundaries, respectively.
In the VOF method, the interaction between the fluids and the solid surface (wall adhesion) is included as a boundary condition on the surface normal vector at the solid wall. nˆ c = nˆ s cos θ + tˆs sin θ
(5)
(10)
Here, nˆ c is the unit vector normal to the fluid interface calculated from the fluid fraction field as (∇ F ) nˆ c = (11) (|∇ F |)
359
(12)
Here, nˆ s and tˆs are the unit vectors normal and tangential to the solid boundary, respectively; θ is the equilibrium or static contact angle. The surface normal vector is corrected at the contact line to satisfy the specified static contact angle boundary condition following which the interface curvature (Eq. (10)) is calculated using this corrected surface normal vector. In OpenFOAM, a modified fluid fraction advection equation based on a two-fluid approach is solved to ensure boundedness and conservation of the fluid fraction field (F ), which is critical for two-phase flows with high density ratios [18,20,23]. Readers are referred to Deshpande et al. [20] for additional details of the modified fluid fraction advection equation implemented in OpenFOAM. Initial and boundary conditions of quasi-static drainage simulation: The quasi-static drainage computations were performed with water (ρ = 1000 kg/m3 and µ = 1e-3 Pa.s) and air (ρ = 1 kg/m3 and µ = 1.8e-5 Pa.s) as wetting and non-wetting phases, respectively. The surface tension between the wetting and non-wetting phases was specified as σ = 0.0725 N/m. A fixed or static contact angle of θ = 55 degrees was used; the dynamic variation of contact angle with contact line velocity was ignored. At time t = 0 s, the inter-fiber pore space is fully saturated with the wetting liquid (water). Fig. 2 shows a saturated isotropic structure connected to imaginary reservoirs of non-wetting and wetting phases on the top and bottom boundaries, respectively. A fixed gauge pressure of 0 Pa was specified at the outlet (wetting phase reservoir) boundary. The non-wetting fluid was forced into the pore-space by incrementally increasing the pressure at the inlet (non-wetting phase reservoir) in pre-determined steps, starting from 500 Pa. Symmetry boundary condition (normal velocity and shear forces set to zero) was applied on the lateral boundaries of the domain, and a no-slip boundary condition was imposed on the surface of the fibers. The effect of gravity is negligible at the porescale and was not included in our simulations. To perform a quasi-static drainage simulation, air flow (non-wetting phase) through the outlet boundary (wetting phase reservoir) must be prevented [26]. This is because when the non-wetting phase breaks through the wetting-phase reservoir boundary (during a primary drainage), a least resistant pathway for non-wetting phase is formed. At this point, any further increase in non-wetting phase reservoir pressure forces the non-wetting phase to flow through the least resistant path, and no additional Pc S data can be determined. To avoid this, in quasi-static desaturation
360
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Fig. 3. (a) Isometric view of a fibrous structure with a slice at the middle of y-axis (ymid ), and (b) Computational grid corresponding to 10 µm resolution at ymid .
experiment, a frit (with pore-size smaller than the porous medium of interest) is used which makes it possible to measure full the Pc S curve. In our numerical study, we employed a similar approach and used a virtual frit (array of parallel cylinders of diameter 20 µm separated by a distance of 10 µm, as shown in Fig. 3a) at the outlet (wetting phase reservoir) of our computational domain. The back pressure created by the virtual frit prevented the non-wetting phase (air) from escaping the outlet (wetting-phase reservoir) in our simulation. 2.3. Computational grid generation One of the main challenges with porescale simulations using conventional CFD methods is the generation of a computational grid which adequately resolves the tortuous inter-fiber void spaces, and offers the feasibility to perform two-phase flow simulations. In this paper, hexahedral-dominant computational grids were generated for the 3D virtual isotropic fibrous structure using foamProMesh (iconCFD, www.iconcfd.com). The tortuous fluid paths inside the fibrous structure necessitated the use of unstructured grids. As a first step in the meshing procedure, a background hexahedral mesh was created using the desired resolution maintaining a cell-aspect ratio of unity. Subsequently, the unwanted cells inside the solid region (of the fibers) were removed, and a boundary-fitted grid was obtained by iteratively refining and snapping the cells to the fiber surfaces. Fig. 3b shows the computational grid corresponding to 10 µm resolution at a slice in the middle of the y-axis (shown in Fig. 3a). The results of gridconvergence study are presented in Section 3.2. 3. Results and discussion In this section, we first present the results of the verification case (quasi-static drainage between two cylinders) which was studied to check the accuracy of the two-phase flow calculations. Next, we describe the results of our grid convergence study for the microscale simulations. Subsequently, we present the results of quasi-static drainage simulation performed using VOF method in a 3D isotropic fibrous structure. Finally, the Pc -S relationship determined from VOF method is compared with the Pc -S relationship estimated using FM method, and discussed in detail. 3.1. Quasi-static drainage of wetting liquid between two horizontal parallel cylinders: Verification of two-phase flow calculations using VOF method In literature, Princen [27] studied the equilibrium configuration of limited amounts of liquid in horizontal assemblies of parallel cylinders, and proposed analytical solutions for radius of curvature and capillary pressure. We adopted Princen’s work [27] and modified the problem setup to simulate quasi-static drainage of a
wetting liquid between two cylinders. A schematic of the cross section of a liquid column between two infinite horizontal cylinders of radii (r), separated by a distance 2d, is shown in Fig. 4. For this configuration, the radius of curvature (Eq. (13)) of meniscus and capillary pressure (Eq. (14)) were formulated by Princen [27] as
[ R= Pc =
d + r − r cos α
] (13)
cos(α + θ )
σ R
=
σ r
[
cos(α + θ ) d/r + 1 − cos α
] (14)
Here, σ is the surface tension between wetting and non-wetting phases, θ is the contact angle made by wetting liquid on solid cylinders, α is the angle between the line connecting the centers of the cylinders and the radius to the solid–liquid–air interface, R is the curvature of the interface, and Pc is the difference between non-wetting and wetting phase reservoir pressures. In the present study, the radii of the cylinders (r) and the separation distance (2d) were fixed at 1 mm and 0.2 mm, respectively, and a static contact angle (θ ) of 55 degrees was used. The quasi-static drainage was simulated for air–water configuration (properties presented in Section 2.2). For a capillary pressure of Pc = 0 Pa (starting point of the drainage process), alpha angle is α = 35 degrees. For a given geometric configuration (fixed d/r ratio) of parallel cylinders and contact angle (θ ), there exists a critical value for alpha αcr (and thus a critical capillary pressure) beyond which a constant contact angle cannot be maintained; at this point, capillary pressure at the interface will not be sufficient to oppose any disturbances created at the interface resulting in complete drainage of the wetting liquid. Critical alpha (αcr ) can be determined by differentiating Eq. (14) with respect to alpha (α ) and c = 0), resulting in equating it to zero ( dP dα
αcr = sin−1
(
sin θ d/ r + 1
)
−θ
(15)
The critical alpha (from Eq. (15)) and critical capillary pressure (from Eq. (14)) values calculated are αcr = −6.87 degrees and 448.36 Pa, respectively. Hence, as capillary pressure increases from 0 Pa to 448.36 Pa, alpha angle (α ) decreases from 35 degrees to −6.87 degrees, at which point the wetting phase will be completely drained out of the system. Any value of capillary pressure greater than 448.36 Pa will instantly drain out the wetting phase between two cylinders. To evaluate the accuracy of the two-phase flow calculations and contact angle model implemented in OpenFOAM, we simulated the drainage of wetting phase between two cylinders in a quasistatic manner. Fig. 4b shows the computational domain for the verification problem with boundary conditions used in the numerical simulation. Since the domain is symmetric, computations were carried out only on one half of the domain, and symmetry BC was applied at the mid-plane between the cylinders. No-slip boundary condition was imposed on the cylinder surface. The gauge pressure
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
361
Fig. 4. Liquid drainage between two horizontal parallel cylinders. (a) Schematic of cross section of two infinite horizontal cylinders of radii (r) separated by a distance 2d, (b) Computational domain with boundary conditions.
at the lower boundary (water reservoir) was fixed at 0 Pa. A grid convergence study was performed for this verification case using four grid resolutions (d/6, d/12, d/18, d/24; d = 0.1 mm) to ensure the results are grid-independent (not included in the manuscript); the difference between numerical results and analytical solution was found to be less than 2%. From the grid-convergence study, a computational grid with a grid-resolution of d/12 (8.3 µm) was used for the quasi-static drainage simulation. Quasi-static drainage of wetting liquid was simulated by incrementally increasing the pressure of the top boundary (air reservoir) from 0 Pa to 448.36 Pa (critical Pc for the configuration considered) with a relaxation phase between successive pressure steps to allow the system to reach equilibrium as shown in Fig. 5a. From the numerical results (wetting-phase saturation), alpha angle was estimated, and plotted as a function of time in Fig. 5b. It can be observed that, for each pressure step, the alpha angle decreases and then reached an equilibrium value (indicated by gray circles). Fig. 5c shows the numerical (equilibrium values from Fig. 5b) and analytical values of alpha for a range of capillary pressures (Fig. 5a) considered in this study; the analytical values of alpha angle for a given capillary pressure or vice-versa can be determined from Eq. (14). It can be seen that for all the Pc values, the alpha values determined from the simulation agree well with the analytical values, except near critical alpha (αcr ), a point at which wetting liquid completely drains out of the system. In the numerical simulation, the wetting phase drainage through the cylinders occurred slightly before the analytical prediction of critical capillary of 448.36 Pa. This is because, the analytical solution is valid only at the equilibrium, and do not take into account the dynamic aspects of the interface propagation. In a quasi-static simulation performed using the VOF method, each pressure step is akin to a dynamic simulation (with applied pressure) during which wetting liquid drains out until the applied pressure and inertial effects (dynamic pressure) are balanced by the capillary pressure at the interface (contact angle) and viscous forces on the cylinder; at equilibrium (end of each pressure step), the inertial forces become negligible. Near the critical alpha angle, the dynamic effects (captured in simulations) resulted in an early drainage of wetting liquid through the cylinders. The good agreement between numerical and analytical results establishes confidence in the predictive capabilities and accuracy of the VOF method used in this study. In the next section, we present the results of the grid-convergence study of microscale simulations of quasi-static drainage in 3D isotropic fibrous structure.
3.2. Grid-convergence study: Drainage simulation using volume-offluid (VOF) method Since the two-phase computations using VOF method are sensitive to the grid resolution [18–20], a grid convergence study was performed to ensure that the results were grid-independent. In order to estimate the effects of discretization, drainage of wetting liquid through an isotropic fibrous structure was simulated using computational grids with different resolutions (50 µm, 25 µm, 16.33 µm, 12.5 µm, 10 µm, and 8 µm). A cell-aspect ratio close to one, which is desired for stable two-phase computations, was used. A grid refinement factor of at least 1.25 was maintained between successive grids. An isometric scaling (N ∼ ∆x−3 ) was observed between the number of cells (N) and grid resolution (∆x), indicating uniform internal refinement. For these simulations, the fluid was initialized (t = 0 s) such that inter-fiber space was completely saturated with the wetting liquid. For this case, a constant pressure differential was applied in the simulation by specifying gauge pressure values of 900 Pa and 0 Pa at the inlet (non-wetting phase reservoir) and outlet (wetting phase reservoir), respectively. Fig. 6 shows the variation in wetting phase saturation as a function of time for different grid resolutions. It was observed that the liquid saturation decreases and reaches an equilibrium value. From Fig. 6, it can be noticed that the variation in saturation levels was found to be minimal between finest grid (8 µm), finer or penultimate grids (10 µm), and fine grid (12.5 µm). At t = 0.01 s, the difference between the saturation levels predicted by finest grid (8 µm) and finer or penultimate grid (10 µm) was 0.009%, and the difference between finest grid (8 µm) and fine grid (12.5 µm) was 1.54%. Based on these results, the finer grid resolution (10 µm) was deemed sufficient, and was used for the quasi-static drainage numerical simulation. 3.3. Quasi-static drainage numerical simulation 3.3.1. Volume-of-fluid (VOF) method At the beginning (t = 0 s) of the quasi-static drainage simulation, the inter-fiber space was fully saturated with the wetting liquid (water). The pressure of the wetting phase reservoir was fixed at 0 Pa during the simulation. As described in Section 2.2, the quasi-static drainage was simulated by increasing the pressure of the non-wetting phase reservoir incrementally. Fig. 7a shows the pressure applied at the inlet (non-wetting phase boundary) or pressure differential between non-wetting and wetting phase reservoirs as a function of time. The wetting phase saturation in
362
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Fig. 5. (a) Pressure applied at the top or non-wetting phase boundary as a function of time, (b) Transient variation of alpha angle in simulations; equilibrium alpha values determined for each pressure step are indicated by gray circles, and (c) Comparison of numerically determined alpha values with analytical alpha values for a range of capillary pressure values considered in this study.
Fig. 8 shows the fluid fraction contours at a slice in the middle of y-plane (Fig. 3a) for different wetting phase saturations (a) S = 1, (b) S = 0.83, (c) S = 0.47, and (d) S = 0.28. It can be observed that, with the increase in capillary pressure, the non-wetting phase intrudes into the complex inter-fiber space displacing the wetting phase. In the VOF method, no modeling assumptions are made on the liquid drainage process, and at each pressure step, the non-wetting phase intrusion process is dictated by the smaller inter-fiber gaps, pore connectivity, and complex interplay between applied pressure differential, inertial, viscous and capillary forces. The P-c -S relationship determined from the VOF method was compared with the results from the FM method (implemented in GeoDict software), a commonly used quasi-static reduced-order geometric approach explained in the next Section 3.3.2. Fig. 6. Effect of grid resolution on wetting phase saturation in drainage simulation performed using volume-of-fluid (VOF) method; transient variation of wetting phase saturation was shown for six different grid resolutions ∆x = 50 µm, 25 µm, 16.33 µm, 12.5 µm, 10 µm, 8 µm.
the inter-fiber space was monitored during the simulations and is shown in Fig. 7b. For each pressure increment, the wetting phase saturation decreases initially, and then reaches an equilibrium value (indicated by gray circles in Fig. 7b). The equilibrium values of wetting phase saturation obtained at the end of each pressure step were used to plot the quasi-static Pc -S relationship. We present the capillary pressure saturation relationship (Pc -S) determined from the VOF method in Section 3.4.
3.3.2. Full-morphology method The full-morphology (FM) method implemented in GeoDict software was used to simulate the quasi-static drainage in a 3D isotropic fibrous structure for air–water configuration (properties presented in Section 2.2). In FM method, the non-wetting fluid was forced into the saturated inter-fiber space by incrementally increasing the pressure of the non-wetting phase reservoir. It is assumed that only certain pore radii (r), as described by the Young– Laplace equation (Eq. (16)), can be intruded by the non-wetting phase at a given capillary pressure (Pc ), Pc =
2σ cosθ r
(16)
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
363
Fig. 7. (a) Pressure applied at the inlet or non-wetting phase boundary as a function of time; a gauge pressure of 0 Pa was applied at the outlet or wetting phase boundary. (b) Transient variation of wetting saturation inside the isotropic fiber structure.
Fig. 8. Fluid-fraction contours for different wetting phase saturations on a slice at the middle of y-axis of the model; (a) S = 1, (b) S = 0.83, (c) S = 0.47, and (d) S = 0.28.
Briefly, the two steps in the FM method are as follows: Erosion step: For each pressure increment, a spherical structural element of radius r, determined using the Young–Laplace equation at a given Pc , was used in a morphological operation of the porespace. This erosion operation renders all possible center points in a domain where a given sphere could be fitted without intruding the fiber boundaries. The erosion operation was performed without any consideration to the connectivity of the fluids to their respective reservoirs [14]. In order to determine Pc -S relationship, a connectivity analysis was performed to remove all those eroded spaces which are not connected to the fluid reservoirs. Dilation step: Following the connectivity analysis, the pore volume occupied by the wetting phase was determined by dilating the structuring element, and the above procedure was repeated by incrementing pressure in predetermined steps [14,28]. Fig. 9 shows intermediate stages of the quasi-static drainage simulation: (a) Pc = 475 Pa, S = 1, (b) Pc = 1232 Pa, S = 0.62, (c) Pc = 1584 Pa, S = 0.3, and (d) Pc = 3696 Pa, S = 0.04. In Fig. 9, the wetting phase (water) is shown in dark gray, the fibers are show in light gray and the non-wetting phase is not colored. It can be seen that the non-wetting phase intrudes the inter-fiber space with the increase in capillary pressure, thereby displacing the wetting phase (water) from the structure. The results of the quasi-static Pc -S relationship
determined from the FM method are presented along with the results from VOF method in the next section. 3.4. Comparison of capillary pressure-saturation (Pc -S) relationship determined using volume-of-fluid (VOF) method and full-morphology (FM) method Fig. 10 shows the quasi-static Pc -S relationship for primary drainage of a 3D isotropic fibrous medium determined using the VOF method along with the results from the FM method, which is a simplified porescale technique. It can be seen that, for 0.15 < S < 0.8, a reasonably good agreement was observed between the Pc -S relationships predicted by the VOF and FM methods. However, for S > 0.9, the FM method predicted high entry pressures as compared to the VOF method (Fig. 10); this is because the VOF method takes into account the dynamic aspects of liquid propagation whereas FM method does not account for the inertial effects which could assist in the drainage of wetting phase. Also, for Pc > 2500 Pa, the FM method predicted slightly lower residual wetting phase saturation as compared to the VOF method (Fig. 10). This could be attributed to the lack of contact line physics (and surface tension force of the interface) in the FM method. At lower saturations, this becomes especially important as liquid blobs could be trapped in small pore spaces, and their effective
364
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Fig. 9. Intermediate stages of quasi-static drainage simulation performed using full-morphology method: (a) P c = 475 Pa, S = 1, (b) P c = 1232 Pa, S = 0.62, (c) P c = 1584 Pa, S = 0.3, and (d) P c = 3696 Pa, S = 0.04. The solid fibers are shown in light gray color, the wetting phase (water) is shown in dark gray color.
capillary pressure could be greater than predicted by the Young– Laplace equation (with a spherical structural element of radius r) used in FM method. Even though the FM method uses simpler set of rules (morphological openings) for liquid propagation in the simulations, the Pc -S determined from the FM method compared reasonably well with the results from direct simulations performed using VOF method (in which no modeling assumptions are made on the liquid propagation). The results indicate that the simplified equations in the FM method reasonably represent multi-phase flow dynamics of quasi-static primary drainage process. In this fashion, the numerical results from direct simulations such as VOF method can be used to verify the accuracy of the modeling assumptions in the reduced-order porescale techniques.
pressure-saturation (Pc -S) relationship for primary drainage in a 3D isotropic fibrous medium. The accuracy of the two-phase VOF calculations was verified by simulating a test case of quasi-static drainage between two infinite parallel cylinders; the numerical results agreed well with the analytical solution. Next, a grid convergence study was performed for quasi-static drainage process in a 3D isotropic fibrous structure, and an optimal grid needed for the microscale simulation was determined. Using this grid, the quasistatic drainage simulation was conducted, and Pc -S relationship was calculated. The Pc -S relationship determined from the direct simulations using the VOF method was compared with full-morphology (FM) method, a popularly used reduced-order approach [4,14]. Good agreement was observed between the results of VOF and FM methods for wetting phase saturation levels of S = 0.15–0.8. The results show that the FM method which relies on quasi-static geometric approach, a relatively simpler physics as compared to VOF method, may be sufficient for estimating Pc -S relationship for a quasi-static primary drainage process in a fibrous medium. Unlike reduced-order modeling approaches (e.g. pore-network models), the VOF method allows for a detailed description of the porous structure, and honors fluid flow conservation equations. The liquid propagation in the VOF method is governed by complex interplay between applied pressure differential, viscous and capillary forces at the microscale. Since no modeling assumptions are made on the geometry of the pore structure and liquid transport, the direct simulations with VOF method are computationally challenging to perform. However, these direct simulations (VOF, SPH, LBM etc.) provide reliable results that can be used for verification and calibration of reduced-order models, which can be used to study physical problems over large length and time scales. The direct simulations can also be used to gain valuable insights into the effect of various porescale processes on macroscale properties under both quasi-static and dynamic conditions.
4. Conclusions
Acknowledgments
In the present study, we performed direct numerical simulations at the microscale using the finite-volume based volumeof-fluid (VOF) method to determine the quasi-static capillary
This work is sponsored by The Procter & Gamble Company through the University of Cincinnati Simulation Center. The computations were performed on a parallel computer cluster at the
Fig. 10. Comparison of capillary pressure-saturation relationship for primary drainage in a 3D isotropic fibrous structure determined using the volume-of-fluid (VOF) and full-morphology (FM) methods.
N.K. Palakurthi et al. / European Journal of Mechanics / B Fluids 67 (2018) 357–365
Procter and Gamble Company consisting of 2.93 GHz Intel Xeon processors (X5670) with 2 Gbytes of memory per processor. The authors thank Dr. Thomas Baer from the Procter & Gamble Company for his valuable suggestions. The authors also thank the anonymous referees for their insightful comments and suggestions which have significantly improved the manuscript. References [1] A. Ashari, T. Bucher, H. Vahedi Tafreshi, A semi-analytical model for simulating fluid transport in multi-layered fibrous sheets made up of solid and porous fibers, Comput. Mater. Sci. 50 (2) (2010) 378–390. [2] C.J. Hong, J.B. Kim, A study of comfort performance in cotton and polyester blended fabrics. I. Vertical wicking behavior, Fibers Polym. 8 (2) (2007) 218– 224. [3] M. Amiot, M. Lewandowski, P. Leite, M. Thomas, A. Perwuelz, An evaluation of fiber orientation and organization in nonwoven fabrics by tensile, air permeability and compression measurements, J. Mater. Sci. 49 (1) (2014) 52–61. [4] A. Ashari, T. Bucher, H.V. Tafreshi, M. Tahir, M. Rahman, Modeling fluid spread in thin fibrous sheets effects of fiber orientation, Int. J. Heat Mass Transfer 53 (9) (2010) 1750–1758. [5] L.A. Richards, Capillary conduction of liquids through porous mediums, Physics 1 (5) (1931) 318–333. [6] M. Hilpert, C.T. Miller, Pore-morphology-based simulation of drainage in totally wetting porous media, Adv. Water Resour. 24 (3) (2001) 243–255. [7] V. Joekar-Niasar, S.M. Hassanizadeh, H. Dahle, Non-equilibrium effects in capillarity and interfacial area in two-phase flow: dynamic pore-network modelling, J. Fluid Mech. 655 (1) (2010) 38–71. [8] V. Joekar-Niasar, M. Prodanović, D. Wildenschild, S. Hassanizadeh, Network model investigation of interfacial area, capillary pressure and saturation relationships in granular porous media, Water Resour. Res. 46 (6) (2010). [9] P.P. Mukherjee, Q. Kang, C.-Y. Wang, Pore-scale modeling of two-phase transport in polymer electrolyte fuel cells—progress and perspective, Energy Environ. Sci. 4 (2) (2011) 346–369. [10] M.L. Porter, M.G. Schaap, D. Wildenschild, Lattice-Boltzmann simulations of the capillary pressure–saturation–interfacial area relationship for porous media, Adv. Water Resour. 32 (11) (2009) 1632–1640. [11] V.P. Schulz, J. Becker, A. Wiegmann, P.P. Mukherjee, C.-Y. Wang, Modeling of two-phase behavior in the gas diffusion medium of PEFCs via full morphology approach, J. Electrochem. Soc. 154 (4) (2007) B419–B426. [12] H.-J. Vogel, J. Tölke, V. Schulz, M. Krafczyk, K. Roth, Comparison of a latticeBoltzmann model a full-morphology model and a pore network model for determining capillary pressure–saturation relationships, Vadose Zone J. 4 (2) (2005) 380–388.
365
[13] M.S. Riasi, N.K. Palakurthi, C. Montemagno, L. Yeghiazarian, A feasibility study of the pore topology method (PTM), A medial surface-based approach to multiphase flow simulation in porous media, Transp. Porous Media 1 (2016) 1–21. [14] A. Ashari, H. Vahedi Tafreshi, A two-scale modeling of motion-induced fluid release from thin fibrous porous media, Chem. Eng. Sci. 64 (9) (2009) 2067– 2075. [15] A. Nabovati, E.W. Llewellin, A. Sousa, A general model for the permeability of fibrous porous media based on fluid flow simulations using the lattice Boltzmann method, Composites 40 (6) (2009) 860–869. [16] A.M. Tartakovsky, P. Meakin, Pore scale modeling of immiscible and miscible fluid flows using smoothed particle hydrodynamics, Adv. Water Resour. 29 (10) (2006) 1464–1478. [17] Y. Zhu, P.J. Fox, J.P. Morris, A pore-scale numerical model for flow through porous media, Int. J. Numer. Anal. Methods Geomech. 23 (9) (1999) 881–904. [18] A. Ferrari, I. Lunati, Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy, Adv. Water Resour. 57 (2013) 19– 31. [19] A. Ashish Saha, S.K. Mitra, Effect of dynamic contact angle in a volume of fluid (VOF) model for a microfluidic capillary flow, J. Colloid Interface Sci. 339 (2) (2009) 461–480. [20] S.S. Deshpande, L. Anumolu, M.F. Trujillo, Evaluating the performance of the two-phase flow solver interfoam, Comput. Sci. Discov. 5 (1) (2012) 014016. [21] D.A. Hoang, V. van Steijn, L.M. Portela, M.T. Kreutzer, C.R. Kleijn, Benchmark numerical simulations of segmented two-phase flows in microchannels using the volume of fluid method, Comput. & Fluids 86 (2013) 28–36. [22] A.Q. Raeini, M.J. Blunt, B. Bijeljic, Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method, J. Comput. Phys. 231 (17) (2012) 5653–5668. [23] E. Berberović, N.P. van Hinsberg, S. Jakirlić, I.V. Roisman, C. Tropea, Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution, Phys. Rev. E 79 (3) (2009) 036306. [24] N.K. Palakurthi, S. Konangi, U. Ghia, K. Comer, Micro-scale simulation of unidirectional capillary transport of wetting liquid through 3D fibrous porous media: Estimation of effective pore radii, Int. J. Multiph. Flow. 77 (2015) 48–57. [25] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface-tension, J. Comput. Phys. 100 (2) (1992) 335–354. [26] B. Ahrenholz, J. Tölke, P. Lehmann, A. Peters, A. Kaestner, M. Krafczyk, W. Durner, Prediction of capillary hysteresis in a porous material using latticeBoltzmann methods and comparison to experimental data and a morphological pore network model, Adv. Water Resour. 31 (9) (2008) 1151–1173. [27] H. Princen, Capillary phenomena in assemblies of parallel cylinders III. Liquid columns between horizontal parallel cylinders, J. Colloid Interface Sci. 34 (2) (1970) 171–184. [28] J. Becker, A. Wiegmann, V. Schulz, Design of fibrous filter media based on the simulation of pore size measures, in: Proceedings of the FILTECH 2007 conference, 2007.