- Email: [email protected]
Micro-scale simulation of unidirectional capillary transport of wetting liquid through 3D fibrous porous media: Estimation of effective pore radii
International Journal of Multiphase Flow 77 (2015) 48–57
Contents lists available at ScienceDirect
International Journal of Multiphase Flow journal homepage: www.elsevier.com/locate/ijmultiphaseflow
Micro-scale simulation of unidirectional capillary transport of wetting liquid through 3D fibrous porous media: Estimation of effective pore radii Nikhil Kumar Palakurthi a,∗, Santosh Konangi a, Urmila Ghia a, Ken Comer b a b
Department of Mechanical and Materials Engineering, University of Cincinnati, 598 Rhodes Hall, PO Box 210072, Cincinnati, Ohio, OH 45221-0072, USA Procter and Gamble, Cincinnati, Ohio, OH 45201, USA
a r t i c l e
i n f o
Article history: Received 1 April 2015 Revised 15 July 2015 Accepted 30 July 2015 Available online 7 August 2015 Keywords: Fibrous porous media Capillary flow Lucas–Washburn equation Micro-scale simulation Effective pore radius Volume-of-fluid method
a b s t r a c t The effective and representative pore radii (R50 ) determined from the physical experiments such as capillary penetration and porosimetry were found to be different for porous materials like ceramic and particulate structures. These differences were attributed to contact angle variations, possibly induced by contaminations. In this paper, we studied the relationship between effective and representative pore radii using “numerical experiments”, in which the contact angle variations were eliminated by using a constant contact angle. Micro-scale simulations of capillary penetration of a wetting liquid (air–water configuration) through virtual fibrous media were performed using a finite-volume-based volume-of-fluid (VOF) method. From √ the simulations, it was found that unidirectional liquid penetration followed Lucas–Washburn kinetics (L∼ t) at the micro-scale, except during the initial stages that are dominated by the inertial forces. Even though the contact angle was held constant in the numerical simulations, the effective pore radii of the fibrous structures differed significantly from the representative pore radii (R50 ). The results suggest that a unique value for the characteristic pore radius does not exist for a porous medium, and the value of characteristic pore radius is strongly dependent on the experimental technique used to characterize the medium. Thus, the choice of experimental technique (e.g. porosimetry or capillary penetration experiment) to determine a characteristic pore radius should be based on the end application of the porous material. © 2015 Elsevier Ltd. All rights reserved.
Introduction Fibrous porous materials are used in a wide variety of applications as they are soft, structurally flexible, and have relatively high resistance to mechanical deformation (Ashari and Vahedi Tafreshi, 2009; Pradhan et al., 2012; Soltani et al., 2014). Capillary penetration of a liquid through a fibrous porous medium is an important process that occurs in many applications, including, but not limited to, printing, drug delivery patches, textile industries, sanitary wipes, and performance fabrics (sports apparel) (Ashari and Vahedi Tafreshi, 2009; Hong and Kim, 2007; Hyvaluoma et al., 2006; Pradhan et al., 2012; Simile, 2004). Understanding the kinetics of capillary penetration of liquid is essential for the characterization of fibrous media, and for the optimization and development of new products (Jaganathan et al., 2008b). Historically, unidirectional capillary transport of a liquid in porous media with a distinct liquid propagating front (i.e. saturated flow ∗
Corresponding author. Tel.: +513 410 6226; fax: +513 556 3390. E-mail address: [email protected], [email protected] Palakurthi). http://dx.doi.org/10.1016/j.ijmultiphaseflow.2015.07.010 0301-9322/© 2015 Elsevier Ltd. All rights reserved.
(N.K.
behind the wetting front) is described using the bundle-of-capillaries approach (Hyvaluoma et al., 2006; Lu et al., 2013; Zhmud et al., 2000). In this approach, a porous medium is assumed to consist of bundle of parallel capillary tubes of same radii, and the kinetics of liquid penetration through the medium is described by the Lucas–Washburn equation (Eq. 1).
L=
σ Rc cosθ t. 2μ
(1)
Here, L is the location of liquid propagating front, Rc is equivalent capillary radius of the porous medium, θ is the contact angle made by the wetting liquid on the solid fiber, t is the wetting time, σ is the surface tension between wetting and non-wetting phases, and μ is the dynamic viscosity of the wetting liquid. Derivation of the Lucas–Washburn equation is provided in Appendix A. The Lucas– Washburn equation has found applications in a variety of research areas, such as contact angle measurement, development of viscometers, and characterization of porous media (O’Loughlin et al., 2013; Xue et al., 2006). Additionally, a modified Lucas–Washburn equation (derivation is provided in Appendix B) that represents the viscous
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
losses using the Darcy law has been proposed in the literature.
L=
4σ kcosθ μϕ Rs
t,
(2)
where, Rs , k, and φ are static radius, permeability and porosity of the medium, respectively. It should be noted that the capillary penetration models described in Eqs. (1) and (2) are based on simplifying assumptions, and the applicability of these models to real porous media with tortuous pore space is still under investigation (Hyvaluoma et al., 2006; O’Loughlin et al., 2013). √ Lucas–Washburn kinetics (L∼ t): In the unidirectional capillary penetration experiments conducted at the macro-scale, the distance traveled by the wicking front or the height rise was found to be proportional to the square root of time (Fries and Quéré, 2010; Hyvaluoma et al., 2006; Masoodi, 2010). Typically, these experiments are carried out at length scales much larger than the structural inhomogeneities in a porous medium (Hyvaluoma et al., 2006). One of the interesting questions is to understand if the liquid imbibition at the micro-scale (length scales close to and even below√ those of structural variations) follows Lucas–Washburn kinetics (L∼ t) (Bousfield and Karles, 2004; Gane et al., 2004; Ridgway and Gane, 2002). If the liquid √ imbibition follows Lucas–Washburn kinetics at all scales, then L∼ t behavior might be intrinsic to the propagation of meniscii in a capillary penetration process. In the present effort, liquid penetration kinetics at the micro-scale (fiber-level) was studied by simulating unidirectional capillary transport of a wetting liquid through 3D virtual fibrous structures using a finite-volume-based volume-of-fluid (VOF) method, which relies only on conservation principles and allows for explicit description of geometry and interface dynamics (Berberovic´ et al., 2009; Deshpande et al., 2012; Ferrari and Lunati, 2013). Characterization of porous media: In applications where the rate of capillary penetration is important, characterization of porous media usually involves determination of an effective pore radius either equivalent capillary radius (Rc ) or static radius (Rs ), and an effective contact angle. However, the effective pore radius and effective contact angle of fibrous media cannot be independently determined from the results of a single capillary penetration (or wicking) experiment (Lavi et al., 2008; Marmur, 2003). Typically, this is achieved by performing two experiments (two-fluid method). A reference liquid, assumed to form a zero contact angle with the fibrous medium, is used to calculate the effective pore radius, and then the effective contact angle is calculated from experiments with the liquid of interest (Lavi et al., 2008). The pore space of fibrous media can also be characterized using porosimetry techniques, which give a distribution of pore-sizes in a fibrous network (Jaganathan et al., 2008a). A representative pore radius (R50 ), the pore radius at which 50% of the non-wetting liquid intrusion volume occurs, is typically reported as it is considered a well-accepted representative radius for a micro-structure with a distribution of pore features (Schoelkopf et al., 2002). It should be noted that the approach of using representative pore radius (R50 ) to characterize a porous medium with a distribution of pores appears to lack a solid scientific basis. Even though the goal of the characterization effort is to determine an unique characteristic pore radius, the effective pore radii obtained from capillary penetration experiments were typically found to be different from the representative pore radii obtained from pore-size distribution data (Hanžiˇc et al., 2010). Einset et al. (1996) studied liquid imbibition rates into particulate structures of carbon and silicon carbide, and found a discrepancy of 1–2 orders of magnitude between effective pore radius and representative pore radius. Li et al. (1994) conducted wicking experiments in ceramic structures using apolar, low-energy liquids, and reported an effective pore radius smaller than the representative pore radius by a factor of approximately two. Although there is no reason for the effective and representative pore
49
radii of a porous medium to be identical from a physics point of view, a single-valued characteristic pore radius was assumed in literature (Einset et al., 1996; Hanžiˇc et al., 2010; Li et al., 1994; Schoelkopf et al., 2002). The disparity between the values of effective and representative pore radii obtained from experiments were attributed to the variation in contact angle in experiments, possibly induced by contaminations, which to some extent are inevitable (Hamraoui and Nylander, 2002; Li et al., 1994; Schoelkopf et al., 2002). In contrast to the physical experiments, detailed numerical simulations of capillary penetration of a wetting liquid through porous media can be performed by keeping the contact angle constant (thereby eliminating the effects of contact angle variations) to study the relationship between effective and representative pore radii. To the best of the authors’ knowledge, no other study has been reported in the literature to numerically determine the effective pore radius of fibrous media. In the present effort, transient two-phase flow simulations of unidirectional capillary penetration of a wetting liquid (water) through 3D virtual fibrous structures media were performed with a constant contact angle using a finite-volume-based volume-of-fluid (VOF) method. Firstly, it will be investigated whether the liquid penetration through fibrous media follows Lucas–Washburn kinetics at the micro-scale. The penetration rates are then used to determine effective pore radii of the fibrous structures. The effective pore radii calculated from the micro-scale simulations are compared with the representative pore radii estimated from the pore-size distribution data. Methodology In the first part of this section, the construction of virtual 3D fibrous structures and determination of pore-size distributions are described. The flow conservation equations and the corresponding boundary conditions used in the micro-scale simulations are presented later, and finally, the generation of computational grids is discussed. Generation of virtual 3D fibrous structures In this study, virtual fibrous structures with different fiber orientations, namely unidirectional and isotropic structures, were constructed using GeoDict, a voxel-based commercial code (Math2Market GmbH, GeoDict2012R1, www.geodict.com). In GeoDict, the fiber orientation distribution is controlled by a density function p(ψ , ϕ ) defined in polar coordinates as
p(ψ , φ) =
1 4π
β sin ψ , (1 + (β 2 − 1)cos2 ψ)3/2
(3)
where β is an anisotropy parameter, and its value reflects the fiber orientation of the micro-structure generation, ψ ∈ [0,π ) and ϕ ∈ [0,2π ) are through-plane and in-plane angles, respectively (Ashari and Vahedi Tafreshi, 2009). The density function (Eq. 3) depends on through-plane angle only, and the fiber orientation is controlled by the value of β . For an isotropic structure, the value of β is equal to one. As the value of β increases, the fibers tend to become parallel to the in-plane. The degree of fiber alignment is described by the secondorder fiber orientation tensor (Soltani et al., 2014; Stylianopoulos et al., 2008),
=
1 li ltot
⎡
sin ψ cos2 φ × ⎣sin2 ψ sin φ cos φ cos ψ sin ψ cos φ 2
sin ψ sin φ cos φ 2 2 sin ψ sin φ cos ψ sin ψ sin φ 2
⎤
cos ψ sin ψ cos φ cos ψ sin ψ sin φ ⎦, cos2 ψ
(4)
where li is the length of the ith fiber, ltot is the total fiber length, and the sum is over all fibers within the domain. The values of diagonal components in the fiber orientation tensor (Eq. 4) are a measure
50
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
Fig. 1. Isometric view of virtual 3D fibrous structures considered in this study: (A) Unidirectional structure, and (B) Isotropic structure; Histograms of geometric and porosimetric pore-size distribution obtained using Full Morphology (FM) method for (C) Unidirectional structure, and (D) Isotropic structure.
of fiber alignment in the principal directions, while the off-diagonal components indicate alignment in a direction other than a principal direction (Soltani et al., 2014; Stylianopoulos et al., 2008). The trace of is always 1. For an unidirectional structure in which all the fibers are parallel to x-axis, xx = 1; yy = zz = 0, whereas for an isotropic structure, xx = yy = zz = 0.33. Unidirectional and isotropic fibrous structures were constructed using infinitely long cylindrical fibers with a diameter (d) of 50 μm in a domain of size 4 mm × 1 mm × 1 mm. The porosity of the fibrous structures was fixed at 0.4, and a resolution of 5 μm per voxel was used. The unidirectional and isotropic fibrous structures are shown in Fig. 1A and B, respectively. The virtual fibrous structures were exported from GeoDict as stereo-lithography (STL) files to create computational grids needed for the micro-scale simulations. It should be noted that the fibrous materials used in most of the industrial and consumer applications are highly porous (porosity > 0.65) as compared to the fibrous structures considered in this study (Hong and Kim, 2007; Soltani et al., 2014). The fibers were allowed to overlap and penetrate through each other during geometry creation. So, although cylindrical fibers were used to create the structures, the final structures contained objects with irregular shapes. However, it should be stated that each virtual structure (in their final form) has been characterized both by porosimetry (full morphology method) and capillary penetration numerical experiments to study the relationship between effective and representative pore radii. The use of low porous fibrous structures (without √ regularity in fiber shape) may affect or shift the inception of L∼ t behavior
(Schoelkopf et al., 2002), but it will not impact the overall liquid penetration kinetics at the micro-scale. Determination of pore-size distribution of the fibrous structures The full-morphology (FM) method implemented in GeoDict was used to determine pore-size distribution of the fibrous structures (Becker et al., 2007). The FM method assumes that the fibrous medium is connected at one end to a large imaginary reservoir of non-wetting fluid. The non-wetting fluid is forced into the pore-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. 5), can be intruded by the nonwetting phase at a given capillary pressure (Pc ),
Pc =
2σ cosθ . r
(5)
Erosion step: For each pressure increment, a spherical structural element of radius ‘r’, determined using the Young–Laplace equation at a given Pc , is used in a morphological operation of pore space. This erosion operation gives all possible center points in a domain where a given sphere could be fitted without intruding the fiber boundaries. Dilation step: In this step, the pore volume intruded by the nonwetting phase is obtained by dilating the structuring element, and the procedure is repeated by incrementing pressure in predetermined steps (Ashari and Vahedi Tafreshi, 2009; Becker et al., 2007). In the fluid-intrusion process, connectivity of the non-wetting phase to the
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
Fig. 2. (A) Schematic of the unidirectional capillary penetration experiment performed in this study; (B) An example of isometric view of the computational domain with a background hexahedral mesh.
reservoir was ignored, and the pore-size distribution thus obtained is called Geometric pore-size distribution (GPSD). In addition to GPSD, the porosimetric pore-size distribution (PPSD) of fibrous structures was determined by performing a connectivity analysis (before the dilation step) in which all the eroded spaces that are not connected to the non-wetting phase reservoirs were removed (Ashari and Vahedi Tafreshi, 2009; Becker et al., 2007). Fig. 1C and D show histograms of GPSD and PPSD for unidirectional and isotropic structures, respectively. It can be seen that the PPSDs are different from the GPSDs because of the requirement of connectivity of the non-wetting phase to the reservoir (connectivity check) in PPSDs before the dilation step. For a unidirectional structure (Fig. 1C), however, the PPSD and GPSD are similar, as the pore gaps are of constant diameters, and are connected throughout the domain. The PPSD of an isotropic structure has a narrower distribution, which could be attributed to the ink-bottle effect (Jaganathan et al., 2008a). The representative pore radii (R50 - median of the curves) of the fibrous structures were estimated from the GPSD and PPSDs. The R50 values obtained from the pore size distributions were compared with the effective pore radii calculated from the two-phase flow simulations in the results and discussion section. Flow field calculation
Single-phase flow: through-plane permeability calculation The flow field inside the inter-fiber void space of a fibrous structure was obtained by solving the steady, 3D incompressible Navier– Stokes equations (Eqs. 6 and 7) using a finite-volume-based computational fluid dynamics software (simpleFoam solver, OpenFOAM). The governing equations are as follows: Continuity (or mass conservation) equation: (6)
Momentum conservation equation:
U · ∇ U = −∇
P ρ
+ υ∇ 2U,
tions, at the inlet, a constant uniform velocity normal to the boundary was specified (Fig. 2B); The velocity at the inlet was chosen such that the value of Reynolds number, defined based on fiber diameter, was less than unity (Re << 1), resulting in a Stokes flow. At the outlet, a constant gauge pressure boundary condition with a fixed value of 0Pa was used. The symmetry boundary condition (normal velocity component and shear forces set to zero) was applied at the side boundaries, as the lateral flow is negligible as compared to the throughplane flow (Jaganathan et al., 2008b; Tahir and Tafreshi, 2009). The no-slip boundary condition was imposed on the surface of all the fibers. From the micro-scale simulation, the average pressure drop ( P) across the thickness of the fibrous medium was obtained, and the through-plane permeability of the virtual fibrous structure was calculated using Darcy law
V =−
P . μ z k
(8)
Here, k is the permeability constant, μ is the dynamic viscosity of the fluid, V is the velocity magnitude at the inlet boundary, and z is the thickness of the fibrous medium. The calculated permeability values were used to determine static radii (Rs ) of the fibrous structures. Two-phase flow: capillary penetration of a wetting liquid The unidirectional capillary penetration (neglecting gravity) of a wetting liquid through fibrous media was simulated using the VOF method, a fluid fraction-based interface-capturing technique used for modeling the flow of a system of two incompressible, immiscible fluids (Ashish Saha and Mitra, 2009; Berberovic´ et al., 2009). The VOF method treats a two-fluid system as a single-fluid system with spacedependent physical properties. In the conventional VOF method, a transport equation for the fluid fraction (Eq. 11) of one phase is solved simultaneously with the continuity (Eq. 9) and momentum equations (Eq. 10). The governing equations are: Continuity equation (mass conservation):
∇ · U = 0,
(9)
Momentum conservation equation:
This section describes the governing equations and boundary conditions for the fluid flow simulations. From Eq. 2, it can be seen that the permeability values of the fibrous structures are required to determine static radii (Rs ). The permeability values were calculated by performing single-phase flow calculations, and the capillary penetration of a wetting liquid through the fibrous structures was simulated using the VOF method. The flow simulations were performed using OpenFOAM, an open-source CFD code (OpenCFD Ltd., v 2.0, www.openfoam.org). Fig. 2A shows the schematic of a unidirectional capillary penetration numerical experiment performed in this study; an isometric view of a computational domain with background hexahedral mesh is shown in Fig. 2B.
∇ · U = 0,
51
(7)
where P is the pressure, U is the velocity vector, ρ and υ are the density and kinematic viscosity, respectively. For the boundary condi-
∂(ρU ) + ∇ · (ρUU ) = −∇ P + ∇ · μ(∇ U + ∇ U T ) + Fσ , ∂t
(10)
Fluid fraction advection equation:
∂α + ∇ · (U α) = 0. ∂t
(11)
Here, U is the velocity vector, α is the fluid fraction, Fσ is the surface tension force, ρ and μ are density and dynamic viscosity of the fluid, respectively. The fluid fraction (α ) is a characteristic function of one of the two fluids, and can take values within the range of 0 ≤ α ≤ 1; α = 0 corresponds to cells filled with the non-wetting phase, whereas α = 1 indicates cells containing the wetting phase. Gradients of the fluid fraction are encountered only in the region of the interface. The physical properties of the fluid are calculated as weighted averages based on the distribution of the fluid fraction function,
μ = μw α + μnw (1 − α),
(12)
ρ = ρw α + ρnw (1 − α).
(13)
Here, ρ w and ρ nw are the densities of the wetting and non-wetting phases, respectively. Whereas μw and μnw are the viscosities of the wetting and non-wetting phases, respectively. In OpenFOAM, a modified fluid fraction equation (or advection equation) based on a two-fluid approach is solved to ensure boundedness and conservation of the fluid fraction (α ); this is essential for obtaining a bounded solution for flow with fluids having high density
52
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
ratios (Berberovic´ et al., 2009; Deshpande et al., 2012). The modified fluid fraction equation is given as
∂α + ∇ · (U α) + ∇ · (Ur α(1 − α)) = 0.
∂t
(14)
CompressionTerm
Here, Ur is the compression velocity. The additional nonlinear convection term, referred to as the “compression term”, is active only near the interface and does not affect the solution outside the interface region. With an appropriate and a suitable choice of parameters, the compression terms ensures mass conservation of the fluid fraction, and guarantees a bounded and unbiased solution. In the momentum equation (Eq. 10), the surface tension force (Fσ ) at the fluid–fluid interface is computed using the continuum surface force (CSF) method (Brackbill et al., 1992). In this method, the force acting at the interface (which is a surface of discontinuity) is represented as a body or continuum force acting over the volume of the region occupied by the interface (0<α <1), and is formulated as (Hoang et al., 2013; Raeini et al., 2012).
Fσ = σ κ∇α ,
(15)
where ’κ ’ is the mean curvature of the free surface, and is determined from
∇α κ =∇· . |∇α|
(16)
In the CSF–VOF method, the interaction between the fluids and the fiber surface (wall adhesion) is included as a boundary condition on the surface normal vector at the solid wall.
nˆ = nˆ s cos θ + tˆs sin θ
(17)
Here, nˆ is the unit vector normal to the interface at the fiber surface, and nˆ s and tˆs are the unit vectors normal and tangential to the solid boundary, respectively; θ is the equilibrium or static contact angle at the triple-contact line (fluid–fluid interface and the solid wall). Two-phase flow simulations: initial and boundary conditions. The transient two-phase flow computations at micro-scale were performed with water (ρ = 1000 kg/m3 and μ = 1e-3Pa–s) and air (ρ = 1 kg/m3 and μ = 1.8 e–6 Pa–s) as the wetting and non-wetting phases, respectively. The surface tension between the wetting and non-wetting phases was specified to be 0.0725 N/m, and a fixed or static contact angle of 55° was used in the micro-scale simulations. At time t = 0 s, the domain was initialized in a manner that the fibers next to inlet (Fig. 2B) were in contact with the wetting phase (water). In order to simulate a capillary-driven flow, a constant gauge pressure condition (0 Pa) was specified at the domain inlet and outlet with standard velocity inlet/outlet boundary conditions. The fluid fraction value was set to unity at the inlet (wetting-phase reservoir), and a zero-gradient condition for fluid fraction was used at the outlet. The symmetry boundary condition (normal velocity and shear forces set to zero) was applied on the side boundaries of the domain, and the no-slip boundary condition was imposed on the surface of the fibers. Accuracy of the two-phase flow calculations performed using OpenFOAM were verified by studying the capillary transport between two infinite parallel plates for which analytical solution is available in literature.
Fig. 3. (A) Isometric view of a fibrous structure with a slice at the middle of z-axis (zmid ); Computational grid corresponding to 10 μm resolution at zmid for fibrous structures: (B) Unidirectional structure, and (C) Isotropic structure.
with a desired resolution and a cell aspect ratio of unity. 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. A grid convergence study was performed to ensure that the results were mesh-independent. Computational grids with a cell-aspect ratio close to one, which is desired for stable two-phase computations using VOF method, were constructed with five different grid resolutions ( x): 25 μm, 16.33 μm, 12.5 μm, 10 μm, 8 μm. An isometric view of a fibrous structure with a slice at the middle of z-axis (zmid ) is shown in Fig. 3A. Fig. 3B and C show the computational grids corresponding to 10 μm at zmid for the unidirectional and isotropic fibrous structures, respectively. 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. The computations were performed on the Procter and Gamble parallel-computing cluster consisting of 2.93 GHz Intel Xeon processors (X5670) with 2 Gigabytes of memory per processor. Transient two-phase flow calculations corresponding to 10 μm grid resolution for an unidirectional (1,911,145 cells) and isotropic structures (1,915,138 cells) were performed on 120 and 144 processors, respectively. These transient simulations for unidirectional and isotropic fibrous structures consumed 43 h and 186 h of wall clock time to simulate physical times of time 0.016 s and 0.044 s, respectively.
Computational grid generation Results and discussion The computational domains were discretized and hexahedraldominant computational grids were generated from the triangulated surface geometries (STL files) using foamProMesh (iconCFD, www.iconcfd.com). The tortuous fluid paths inside the fibrous structures necessitated the use of unstructured grids. As a first step in the meshing procedure, a background hexahedral mesh was created
The first part of this section deals with the kinetics of unidirectional capillary penetration of a wetting liquid. Next, the effective pore radii of fibrous structures calculated from the micro-scale simulations were compared with the representative radii obtained from the pore-size distribution data.
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
53
Fig. 4. Capillary transport between two parallel plates as a function of time: (A) Capillary penetration distance, and (B) Velocity of the wicking front.
Fig. 5. Effect of grid resolution on (A) Mass uptake per unit cross-sectional area, and (B) Mass flux for an unidirectional fibrous structure. Five different grid resolutions (࢞x = 25 μm, 16.33 μm, 12.5 μm, 10 μm, 8 μm) were considered for this study.
Capillary penetration of a wetting liquid: Lucas–Washburn kinetics Here, we first describe the results of capillary transport between two parallel plates, with the objective to verify our two-phase calculations. Later, the results of the grid-convergence study are explained, following which the mass uptake per unit cross-sectional area and mass flux values for different fibrous structures are presented. Finally, the kinetics of liquid penetration at the micro-scale is discussed in detail. Capillary transport between two parallel plates: verification of two-phase flow calculations The accuracy of the two-phase flow calculations was verified by simulating capillary transport between two infinite parallel plates, separated by a distance of 1 mm. The simulation was performed using an air-water configuration (with water as the wetting liquid); a static contact angle of 70° was fixed in the simulation. At time t = 0 s, the liquid meniscus position in the domain was initialized at 2 mm. At the inlet, the fluid fraction value was set to unity, and a zero-gradient condition for the fluid fraction was used at the outlet. Unidirectional capillary-driven flow was simulated by specifying a constant gauge pressure of 0 Pa at both inlet and outlet. Fig. 4 shows the capillary dynamics (average distance travelled and average velocity of the interface with time) obtained from the OpenFOAM computations along with the results of the Lucas– Washburn equation (Eq. 1), included for comparison. It can be seen (Fig. 4A) that the numerical results follow analytical results of the Lucas–Washburn equation, except during the initial stages which are dominated by inertial forces (Fries and Quéré, 2010; Lu et al., 2013). Fries and Quéré (2010) provided an analytical expression
(t1 = 2.1151 R2 ρ /μ = 1.41 s) to calculate the time required to reach the purely viscous stage, in which the Lucas–Washburn equation is valid. At time instant t1 = 1.41 s, the maximum difference between the numerical data and the results of the Lucas–Washburn equation is 1.43%, thereby indicating good agreement of the numerical data with the analytical results. Fig. 4B shows the average velocity of the wicking front predicted by the OpenFOAM calculations and the Lucas–Washburn equation. Here again, good agreement is observed between the numerical and analytical results in the purely viscous stage (after time t1 s). In the next section, the results of gridconvergence study are presented. Grid convergence study: micro-scale simulations Since two-phase flow computations are sensitive to grid resolution, the capillary penetration of a wetting liquid through an unidirectional fibrous structure was simulated on computational grids with different resolutions to estimate the effects of discretization. The mass uptake and mass flow rate of the wetting liquid were monitored during the simulations. Fig. 5 shows the mass uptake per unit crosssectional area and mass flux of the unidirectional fibrous structure for different grid resolutions as a function of time. It can be seen that, at a given time instant, mass uptake and mass flux decrease with an increase in grid resolution. Minimal variation was seen in the results of the fine grid (12.5 μm), finer or penultimate grid (10 μm), and finest grid (8 μm). At t = 0.015 s, difference in the mass uptake values predicted by the finest grid (8 μm) and finer or penultimate grid (10 μm) was 2%, and the difference in the mass uptake values predicted by finer grid (10 μm) and fine grid (12.5 μm) was 4.4%. Based on these results, the finer grid (10 μm) resolution was deemed sufficient, and was used in the subsequent micro-scale computations. In
54
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
Fig. 6. (A) Mass uptake per unit cross-sectional area, and (B) Mass flux as a function of time for unidirectional and isotropic fibrous structures.
Fig. 7. Fluid fraction contours on a representative slice at the middle of z-axis at t = 0.01 s for (A) Unidirectional fibrous structure, and (B) Isotropic fibrous structure.
the next section, the mass uptake and mass flux results for both the fibrous structures are presented. Capillary penetration of a wetting liquid through 3D fibrous structures: mass uptake and mass flux from micro-scale simulations The cumulative mass uptake per unit cross-sectional area and mass flux as a function of time for unidirectional and isotropic fibrous structures are shown in Fig. 6A and B, respectively. From Fig. 6, it can be seen that the cumulative mass uptake of the fibrous structures increases with time (Fig. 6A), while the mass flux of fibrous structures (Fig. 6B) increases initially, reaches a peak, and then decreases as the viscous effects in the flow become dominant. At any given time instant, as expected, the unidirectional structure has the highest mass uptake and mass flux as compared to an isotropic fibrous structure. Fig. 7A and B show the fluid fraction contours on a representative slice at the middle of z-axis at t = 0.01 s for the unidirectional and isotropic fibrous structures, respectively. It is evident that the wicking front propagates faster in a unidirectional structure as compared to an isotropic structure because of better connectivity between the inter-fiber spaces through which the fluid flows (Ashari et al., 2010; Soltani et al., 2014). Kinetics of liquid penetration at micro-scale The average distance (L) travelled by the wicking front with time was calculated from the results of the micro-scale simulations. Fig. 8 shows the location of the wicking front as a function of the square √ root of time ( t) for the fibrous structures. From Fig. 8, it can be
√ Fig. 8. Location (L) of liquid front as a function of square root of time ( t). A linear √ relationship was observed between L and t during the later stages of capillary transport.
seen that the average distance increases non-linearly during the early stages of liquid transport. However, during the later stages √ of liquid transport, a linear relationship is observed between L and t for both the fibrous structures. From these results, it can concluded that the kinetics of liquid penetration at the micro-scale√follow the behavior delineated by Lucas–Washburn equation (L∼ t), after the early time stages of capillary transport which are dominated by inertial forces (Hyvaluoma et al., 2006; Simile, 2004). Recent experimental works carried out at the nano-scale also support√ the validity of the Lucas–Washburn kinetics, which suggest that L∼ t behavior might be fundamental to the capillary penetration mechanism (Stukan et al., 2010). √ It is interesting to note that a viscous stage (L∼ t) occurs in unidirectional capillary transport through porous media at micro-scale, similar to the purely viscous flow stage observed in the capillary transport between two parallel plates. However, unlike the capillary √ transport in parallel plates or a single capillary, manifestation of L∼ t behavior in a porous medium does not mean that inertia plays no role in later stages. In a porous medium, there is an inertial contribution each time the liquid is accelerated (Schoelkopf et al., 2002). In this respect, the picture of a single capillary representing a porous substrate is misleading (Schoelkopf et al., 2002).
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
55
Table 1 Comparison of effective pore radii (Rc and Rs ) obtained from the micro-scale simulations with the representative radii (R50 ) calculated from the porosimetric (R50, P ) pore-size distributions (PSD). Radii
Pore-size distribution Geometric PSD Porosimetric PSD Two-phase flow simulations Bundle-of-capillaries Darcy law approach
Fibrous structure Unidirectional (k = 1.23 e–10 m2 )
Isotropic (k = 4.2 e–11 m2 )
R50, G R50, P
29.28 μm 29.26 μm
31.59 μm 22.98 μm
Rc Rs
47.47 μm 41.57 μm
10.18 μm 66.1 μm
Characterization of fibrous media: estimation of effective pore radius In the present study, the effective pore radius of a fibrous structure was determined using two approaches: bundle-of-capillaries approach (equivalent capillary radius) and Darcy law approach (static radius). Equivalent capillary radius (Rc ) using √ bundle-of-capillaries approach: From Fig. 8, the slopes of the L versus t curves, referred to as penetration rates or wicking coefficients, were determined for all three fibrous structures. Using the penetration rates, wetting phase properties, and contact angle, the equivalent capillary radii (Rc ) of the fibrous structures were calculated (Eq. 1). Once the effective pore radii were determined, the next step is to determine the static radii of the fibrous structures. Static radii (Rs ) using Darcy law approach: In order to determine Rs of a fibrous structure, we need permeability of the medium (Eq. 2) in addition to penetration rate, fluid properties, and porosity. The permeability values of the fibrous structures were calculated as described in the methodology section. The static radii of the fibrous structures were determined using Eq. 2. Table 1 lists the effective pore radii (Rc & Rs ) obtained from capillary penetration simulations, and the representative pore radii obtained from the GPSD (R50, G ) and PPSD (R50, P ) data. The throughplane permeability values of the fibrous structures are also listed in Table 1. The permeability of the unidirectional fibrous structure is greater than that of the isotropic fibrous structure, as the unidirectional structure offers less resistance to the fluid flow (Soltani et al., 2014; Tahir and Tafreshi, 2009; Wang et al., 2007). For the unidirectional fibrous structure, the R50, G and R50, P have approximately the same value, whereas for the isotropic fibrous structure, R50, P is smaller than R50, G because of the narrower PPSDs. From Eqs. 1 and 2, it can be seen that the average distance (L), and thus the mass uptake, of a fibrous structure is directly proportional to Rc , and inversely proportional to Rs . As expected, the unidirectional fibrous structure had the highest Rc and the lowest Rs . The Rs values for both the fibrous structures were greater than the R50 values obtained from the GPSDs and PPSDs. Next, on comparing R50 with Rc , we see that for a unidirectional structure, the values of R50 are smaller than that of Rc , whereas for the isotropic structure, the values of R50 are larger than Rc . The lack of meaningful correlation between the R50 and Rc values of the different fibrous structures can be attributed to the selective pore-filling that occurs because of the differences in the pore connectivity, which influences the liquid transport through the structure (Schoelkopf et al., 2002). As stated in introduction section, in all the simulations, the fluid properties and contact angle were kept constant, thereby removing any variations that may arise from the variation of contact angle, that are prevalent in physical experiments. Even though the Lucas– Washburn coefficients were kept constant in all the simulations, the Rc and Rs values determined from capillary penetration simulations were found to be different from R50 values obtained from pore-size distribution data (Hanžiˇc et al., 2010). Hence, from this study, it can be stated that the differences between effective and representative
pore radii of a fibrous medium do not arise from the contact angle variations. Conclusions In the present effort, direct numerical simulations were conducted to study the kinetics of liquid penetration at the microscale, and the relationship between the effective pore radii (Rc and Rs ) calculated from capillary penetration experiments and the representative radius (R50 ) estimated from pore-size distribution data (porosimetry). In the literature, for a given porous medium, different values of effective and representative pore radii were obtained, and the differences in pore radii were attributed to the contact angle variations that are prevalent in the experiments (possibly induced by contaminations). In contrast to physical experiments, numerical simulations allowed us to perform capillary penetration numerical “experiments” by keeping the contact angle constant, thus eliminating the effect of contact angle variations. Transient two-phase flow simulations of capillary penetration of a wetting liquid through 3D virtual fibrous media were performed with a constant contact angle using the finite-volume-based volume-of-fluid (VOF) method. The simulated kinetics of liquid penetration at the micro-scale fol√ lowed the Lucas–Washburn equation (L∼ t) during the later stages of the capillary transport. The deviation (non-linear behavior) from Lucas–Washburn kinetics during the early stages could be attributed to the inertial effects that are accounted for in the micro-scale simu√ lations. The inception of a purely viscous stage (L∼ t) suggests that the liquid propagation through fibrous structures during the later stages of capillary transport is controlled primarily by the viscous drag caused by the wetted fibers. Next, the effective pore radii of the fibrous structures were calculated using bundle-of-capillaries (equivalent capillary radius) and Darcy law (static radius) approaches from the rates of liquid penetration obtained from the micro-scale simulations. The equivalent capillary radii and static radii of the fibrous structures were compared with the representative pore radii obtained from the pore-size distributions. Even though the fluid properties and contact angle were kept constant in all the simulations, the effective pore radii (equivalent capillary radii and static radii) obtained from the simulations were found to be different from the representative radii. Thus, a unique value for the characteristic pore radius does not exist for a fibrous medium. The value of characteristic pore radius is strongly dependent on the experimental technique (e.g. porosimetry or capillary penetration experiment) used to characterize the medium. The selection of the characterization method should be based on the end application of the material. Acknowledgment This work is sponsored by the Procter & Gamble company through the University of Cincinnati Simulation Center. The authors appreciate the anonymous referee for his/her insightful comments and suggestions which have significantly improved the manuscript.
56
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
APPENDIX A. derivation of Lucas–Washburn equation Applying Newton’s second law on a liquid column rising in a circular tube of radius r, and assuming that flow is incompressible, viscous, laminar, and fully developed, one arrives at the 1D infiltration equation of motion (Fries and Quéré, 2010),
2
ρ LL + L
Inertia
2σ cos θ 8 − 2 μLL − r r
=
Capillary Pr essure
V iscousLosses
ρ gL
,
=
(19)
V iscousLosses
Capillary Pr essure
Integrating Eq. 19, one obtains Lucas–Washburn equation which predicts time evolution of length of the liquid column, and is given as
L=
rσ cos θ 1/2 t . 2μ
(20)
PenetrationRate
The Lucas–Washburn equation is also used to describe the kinetics of liquid penetration in a porous medium (Fries and Quéré, 2010; Masoodi, 2010); here, the cylinder radius r in Eq. 20 is replaced by the equivalent capillary radius of the porous medium Rc as shown in Eq. 21.
L=
σ Rc cosθ t. 2μ
(21)
APPENDIX B. derivation of Darcy Law based Lucas–Washburn equation for porous media The momentum equation of a liquid column rising in circular tube (Eq. 18) can be extended to a generic porous medium by replacing the Hagen–Poiseuille viscous losses using the Darcy law, and is given as (Fries and Quéré, 2010; Masoodi, 2010)
2 2σ cos θ ϕ ρ LL + L = − μLL − ρ gL , Rs k
HydrostaticForce Inertia
Capillary Pr essure
(22)
V iscousLosses
where Rs , k, and φ are static radius, permeability and porosity of the medium, respectively. It should be noted that the radius r in Eq. 18, used to calculate the capillary pressure in a circular tube, is replaced by the static radius Rs in Eq. 22 to emphasize that the porous structure does not consist of capillary tubes with a specific radius r; here, one assumes that the porous medium features a uniform pore size (Fries and Quéré, 2010; Masoodi, 2010). In Eq. 22, neglecting the influence of inertial and gravity forces, one obtains the balance between capillary pressure (driving force) and viscous drag,
2σ cos θ Rs
Capillary Pr essure
=
ϕ
μLL . k V iscousLosses
L=
4σ kcosθ μϕ Rs
t.
(24)
References
HydrostaticForce
8 μLL . 2 r
(18)
where L is the distance travelled by the advancing front, θ is the contact angle, σ is the surface tension of fluid–fluid interface, g is gravity, ρ and μ are density and viscosity of the wetting fluid, respectively. In Eq. 18, the capillary pressure (driving force) is balanced by the sum of inertial, viscous, and hydrostatic forces. Lucas and Washburn investigated the wicking dynamics in a circular tube neglecting the influence of inertial and gravity forces; here, the capillary pressure (driving force) is balanced by the viscous drag.
2σ cos θ r
Integrating Eq. 23, one obtains modified Lucas–Washburn equation which predicts time evolution of length of the liquid column in a porous medium (Fries and Quéré, 2010; Masoodi, 2010)
(23)
Ashari, A., Bucher, T., Tafreshi, H.V., Tahir, M., Rahman, M., 2010. Modeling fluid spread in thin fibrous sheets: effects of fiber orientation. Int. J. Heat Mass Trans. 53, 1750– 1758. Ashari, A., Vahedi Tafreshi, H., 2009. A two-scale modeling of motion-induced fluid release from thin fibrous porous media. Chem. Eng. Sci. 64, 2067–2075. Ashish Saha, A., Mitra, S.K., 2009. Effect of dynamic contact angle in a volume of fluid (VOF) model for a microfluidic capillary flow. J. Colloid Int. Sci. 339, 461–480. Becker, J., Wiegmann, A., Schulz, V., 2007. Design of fibrous filter media based on the simulation of pore size measures. In: Proceedings of the FILTECH 2007 conference. ´ E., van Hinsberg, N.P., Jakirlic, ´ S., Roisman, I.V., Tropea, C. , 2009. Drop imBerberovic, pact onto a liquid layer of finite thickness: dynamics of the cavity evolution. Phys. Rev. E 79, 036306. Bousfield, D.W., Karles, G., 2004. Penetration into three-dimensional complex porous structures. J. Coll. Int. Sci. 270, 396–405. Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surfacetension. J Comput. Phys. 100, 335–354. Deshpande, S.S., Anumolu, L., Trujillo, M.F. , 2012. Evaluating the performance of the two-phase flow solver interFoam. Comput. Sci. Discov. 5, 014016. Ferrari, A., Lunati, I., 2013. Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy. Adv. Water Res. 57, 19–31. Fries, N., Quéré, D., 2010. Capillary Transport Processes in Porous Materials-Experiment and Model. Cuvillier Verlag Göttingen Doctoral dissertation2010; ISBN 978-386955-507-2. Gane, P.A., Ridgway, C.J., Schoelkopf, J., 2004. Absorption rate and volume dependency on the complexity of porous network structures. Transp. Porous Media 54, 79–106. Hamraoui, A., Nylander, T., 2002. Analytical approach for the Lucas–Washburn equation. J. Coll. Interface Sci. 250, 415–421. Hanžiˇc, L., Kosec, L., Anžel, I., 2010. Capillary absorption in concrete and the Lucas– Washburn equation. Cement Concr. Compos. 32, 84–91. Hoang, D.A., van Steijn, V., Portela, L.M., Kreutzer, M.T., Kleijn, C.R., 2013. Benchmark numerical simulations of segmented two-phase flows in microchannels using the Volume of Fluid method. Comput. Fluids 86, 28–36. Hong, C.J., Kim, J.B., 2007. A study of comfort performance in cotton and polyester blended fabrics. I. vertical wicking behavior. Fibers Polym. 8, 218–224. Hyvaluoma, J., Raiskinmaki, P., Jasberg, A., Koponen, A., Kataja, M., Timonen, J. , 2006. Simulation of liquid penetration in paper. Phys. Rev. E Stat. Nonlin. Soft. Matt. Phys. 73, 036705. Jaganathan, S., Vahedi Tafreshi, H., Pourdeyhimi, B., 2008a. Modeling liquid porosimetry in modeled and imaged 3-D fibrous microstructures. J. Coll. Int. Sci. 326, 166– 175. Jaganathan, S., Vahedi Tafreshi, H., Pourdeyhimi, B., 2008b. A realistic approach for modeling permeability of fibrous media: 3-D imaging coupled with CFD simulation. Chem. Eng. Sci. 63, 244–252. Lavi, B., Marmur, A., Bachmann, J., 2008. Porous media characterization by the twoliquid method: effect of dynamic contact angle and inertia. Langmuir 24, 1918– 1923. Li, Z., Giese, R.F., Oss, C.J., Kerch, H.M., Burdette, H.E., 1994. Wicking technique for determination of pore size in ceramic materials. J. American Ceram. Soc. 77, 2220–2222. Lu, G., Wang, X.-D., Duan, Y.-Y., 2013. Study on initial stage of capillary rise dynamics. Coll. Surf. A: physicochem. Eng. Asp. 433, 95–103. Marmur, A., 2003. Kinetics of penetration into uniform porous media: testing the equivalent-capillary concept. Langmuir 19, 5956–5959. Masoodi, R., 2010. Modeling Imbibition of Liquids into Rigid and Swelling Porous Media. The University of Wisconsin–Milwaukee Doctoral dissertation. O’Loughlin, M., Wilk, K., Priest, C., Ralston, J., Popescu, M., 2013. Capillary rise dynamics of aqueous glycerol solutions in glass capillaries: a critical examination of the Washburn equation. J. Coll. Interface Sci. 411, 257–264. Pradhan, A.K., Das, D., Chattopadhyay, R., Singh, S., 2012. Effect of 3D fiber orientation distribution on transverse air permeability of fibrous porous media. Powder Tech. 221, 101–104. Raeini, A.Q., Blunt, M.J., Bijeljic, B., 2012. Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method. J. Comput. Phys. 231, 5653–5668. Ridgway, C.J., Gane, P.A., 2002. Dynamic absorption into simulated porous structures. Coll. Surf. A: physicochem. Eng. Aspects 206, 217–239. Schoelkopf, J., Gane, P.A., Ridgway, C.J., Matthews, G.P., 2002. Practical observation of deviation from Lucas–Washburn scaling in porous media. Coll. Surf. A: physicochem. Eng. Aspects 206, 445–454. Simile, C.B., 2004. Critical Evaluation of Wicking in Performance Fabrics. Georgia Institute of Technology Doctoral dissertation. Soltani, P., Johari, M.S., Zarrebini, M., 2014. Effect of 3D fiber orientation on permeability of realistic fibrous porous networks. Powder Tech. 254, 44–56. Stukan, M.R., Ligneul, P., Crawshaw, J.P., Boek, E.S., 2010. Spontaneous imbibition in nanopores of different roughness and wettability. Langmuir 26, 13342–13352.
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57 Stylianopoulos, T., Yeckel, A., Derby, J.J., Luo, X.-J., Shephard, M.S., Sander, E.A., Barocas, V.H. , 2008. Permeability calculations in three-dimensional isotropic and oriented fiber networks. Phys. Fluids 20, 123601. Tahir, M.A., Tafreshi, H.V., 2009. Influence of fiber orientation on the transverse permeability of fibrous media. Phys. Fluids 21, 083604–083605. Wang, Q., Maze, B., Tafreshi, H.V., Pourdeyhimi, B., 2007. Simulating through-plane permeability of fibrous materials with different fiber lengths. Modell. Simul. Mater. Sci. Eng. 15, 855–868.
57
Xue, H., Fang, Z., Yang, Y., Huang, J., Zhou, L., 2006. Contact angle determined by spontaneous dynamic capillary rises with hydrostatic effects: experiment and theory. Chem. Phys. Lett. 432, 326–330. Zhmud, B., Tiberg, F., Hallstensson, K., 2000. Dynamics of capillary rise. J. Coll. Interface Sci. 228, 263–269.
Contents lists available at ScienceDirect
International Journal of Multiphase Flow journal homepage: www.elsevier.com/locate/ijmultiphaseflow
Micro-scale simulation of unidirectional capillary transport of wetting liquid through 3D fibrous porous media: Estimation of effective pore radii Nikhil Kumar Palakurthi a,∗, Santosh Konangi a, Urmila Ghia a, Ken Comer b a b
Department of Mechanical and Materials Engineering, University of Cincinnati, 598 Rhodes Hall, PO Box 210072, Cincinnati, Ohio, OH 45221-0072, USA Procter and Gamble, Cincinnati, Ohio, OH 45201, USA
a r t i c l e
i n f o
Article history: Received 1 April 2015 Revised 15 July 2015 Accepted 30 July 2015 Available online 7 August 2015 Keywords: Fibrous porous media Capillary flow Lucas–Washburn equation Micro-scale simulation Effective pore radius Volume-of-fluid method
a b s t r a c t The effective and representative pore radii (R50 ) determined from the physical experiments such as capillary penetration and porosimetry were found to be different for porous materials like ceramic and particulate structures. These differences were attributed to contact angle variations, possibly induced by contaminations. In this paper, we studied the relationship between effective and representative pore radii using “numerical experiments”, in which the contact angle variations were eliminated by using a constant contact angle. Micro-scale simulations of capillary penetration of a wetting liquid (air–water configuration) through virtual fibrous media were performed using a finite-volume-based volume-of-fluid (VOF) method. From √ the simulations, it was found that unidirectional liquid penetration followed Lucas–Washburn kinetics (L∼ t) at the micro-scale, except during the initial stages that are dominated by the inertial forces. Even though the contact angle was held constant in the numerical simulations, the effective pore radii of the fibrous structures differed significantly from the representative pore radii (R50 ). The results suggest that a unique value for the characteristic pore radius does not exist for a porous medium, and the value of characteristic pore radius is strongly dependent on the experimental technique used to characterize the medium. Thus, the choice of experimental technique (e.g. porosimetry or capillary penetration experiment) to determine a characteristic pore radius should be based on the end application of the porous material. © 2015 Elsevier Ltd. All rights reserved.
Introduction Fibrous porous materials are used in a wide variety of applications as they are soft, structurally flexible, and have relatively high resistance to mechanical deformation (Ashari and Vahedi Tafreshi, 2009; Pradhan et al., 2012; Soltani et al., 2014). Capillary penetration of a liquid through a fibrous porous medium is an important process that occurs in many applications, including, but not limited to, printing, drug delivery patches, textile industries, sanitary wipes, and performance fabrics (sports apparel) (Ashari and Vahedi Tafreshi, 2009; Hong and Kim, 2007; Hyvaluoma et al., 2006; Pradhan et al., 2012; Simile, 2004). Understanding the kinetics of capillary penetration of liquid is essential for the characterization of fibrous media, and for the optimization and development of new products (Jaganathan et al., 2008b). Historically, unidirectional capillary transport of a liquid in porous media with a distinct liquid propagating front (i.e. saturated flow ∗
Corresponding author. Tel.: +513 410 6226; fax: +513 556 3390. E-mail address: [email protected], [email protected] Palakurthi). http://dx.doi.org/10.1016/j.ijmultiphaseflow.2015.07.010 0301-9322/© 2015 Elsevier Ltd. All rights reserved.
(N.K.
behind the wetting front) is described using the bundle-of-capillaries approach (Hyvaluoma et al., 2006; Lu et al., 2013; Zhmud et al., 2000). In this approach, a porous medium is assumed to consist of bundle of parallel capillary tubes of same radii, and the kinetics of liquid penetration through the medium is described by the Lucas–Washburn equation (Eq. 1).
L=
σ Rc cosθ t. 2μ
(1)
Here, L is the location of liquid propagating front, Rc is equivalent capillary radius of the porous medium, θ is the contact angle made by the wetting liquid on the solid fiber, t is the wetting time, σ is the surface tension between wetting and non-wetting phases, and μ is the dynamic viscosity of the wetting liquid. Derivation of the Lucas–Washburn equation is provided in Appendix A. The Lucas– Washburn equation has found applications in a variety of research areas, such as contact angle measurement, development of viscometers, and characterization of porous media (O’Loughlin et al., 2013; Xue et al., 2006). Additionally, a modified Lucas–Washburn equation (derivation is provided in Appendix B) that represents the viscous
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
losses using the Darcy law has been proposed in the literature.
L=
4σ kcosθ μϕ Rs
t,
(2)
where, Rs , k, and φ are static radius, permeability and porosity of the medium, respectively. It should be noted that the capillary penetration models described in Eqs. (1) and (2) are based on simplifying assumptions, and the applicability of these models to real porous media with tortuous pore space is still under investigation (Hyvaluoma et al., 2006; O’Loughlin et al., 2013). √ Lucas–Washburn kinetics (L∼ t): In the unidirectional capillary penetration experiments conducted at the macro-scale, the distance traveled by the wicking front or the height rise was found to be proportional to the square root of time (Fries and Quéré, 2010; Hyvaluoma et al., 2006; Masoodi, 2010). Typically, these experiments are carried out at length scales much larger than the structural inhomogeneities in a porous medium (Hyvaluoma et al., 2006). One of the interesting questions is to understand if the liquid imbibition at the micro-scale (length scales close to and even below√ those of structural variations) follows Lucas–Washburn kinetics (L∼ t) (Bousfield and Karles, 2004; Gane et al., 2004; Ridgway and Gane, 2002). If the liquid √ imbibition follows Lucas–Washburn kinetics at all scales, then L∼ t behavior might be intrinsic to the propagation of meniscii in a capillary penetration process. In the present effort, liquid penetration kinetics at the micro-scale (fiber-level) was studied by simulating unidirectional capillary transport of a wetting liquid through 3D virtual fibrous structures using a finite-volume-based volume-of-fluid (VOF) method, which relies only on conservation principles and allows for explicit description of geometry and interface dynamics (Berberovic´ et al., 2009; Deshpande et al., 2012; Ferrari and Lunati, 2013). Characterization of porous media: In applications where the rate of capillary penetration is important, characterization of porous media usually involves determination of an effective pore radius either equivalent capillary radius (Rc ) or static radius (Rs ), and an effective contact angle. However, the effective pore radius and effective contact angle of fibrous media cannot be independently determined from the results of a single capillary penetration (or wicking) experiment (Lavi et al., 2008; Marmur, 2003). Typically, this is achieved by performing two experiments (two-fluid method). A reference liquid, assumed to form a zero contact angle with the fibrous medium, is used to calculate the effective pore radius, and then the effective contact angle is calculated from experiments with the liquid of interest (Lavi et al., 2008). The pore space of fibrous media can also be characterized using porosimetry techniques, which give a distribution of pore-sizes in a fibrous network (Jaganathan et al., 2008a). A representative pore radius (R50 ), the pore radius at which 50% of the non-wetting liquid intrusion volume occurs, is typically reported as it is considered a well-accepted representative radius for a micro-structure with a distribution of pore features (Schoelkopf et al., 2002). It should be noted that the approach of using representative pore radius (R50 ) to characterize a porous medium with a distribution of pores appears to lack a solid scientific basis. Even though the goal of the characterization effort is to determine an unique characteristic pore radius, the effective pore radii obtained from capillary penetration experiments were typically found to be different from the representative pore radii obtained from pore-size distribution data (Hanžiˇc et al., 2010). Einset et al. (1996) studied liquid imbibition rates into particulate structures of carbon and silicon carbide, and found a discrepancy of 1–2 orders of magnitude between effective pore radius and representative pore radius. Li et al. (1994) conducted wicking experiments in ceramic structures using apolar, low-energy liquids, and reported an effective pore radius smaller than the representative pore radius by a factor of approximately two. Although there is no reason for the effective and representative pore
49
radii of a porous medium to be identical from a physics point of view, a single-valued characteristic pore radius was assumed in literature (Einset et al., 1996; Hanžiˇc et al., 2010; Li et al., 1994; Schoelkopf et al., 2002). The disparity between the values of effective and representative pore radii obtained from experiments were attributed to the variation in contact angle in experiments, possibly induced by contaminations, which to some extent are inevitable (Hamraoui and Nylander, 2002; Li et al., 1994; Schoelkopf et al., 2002). In contrast to the physical experiments, detailed numerical simulations of capillary penetration of a wetting liquid through porous media can be performed by keeping the contact angle constant (thereby eliminating the effects of contact angle variations) to study the relationship between effective and representative pore radii. To the best of the authors’ knowledge, no other study has been reported in the literature to numerically determine the effective pore radius of fibrous media. In the present effort, transient two-phase flow simulations of unidirectional capillary penetration of a wetting liquid (water) through 3D virtual fibrous structures media were performed with a constant contact angle using a finite-volume-based volume-of-fluid (VOF) method. Firstly, it will be investigated whether the liquid penetration through fibrous media follows Lucas–Washburn kinetics at the micro-scale. The penetration rates are then used to determine effective pore radii of the fibrous structures. The effective pore radii calculated from the micro-scale simulations are compared with the representative pore radii estimated from the pore-size distribution data. Methodology In the first part of this section, the construction of virtual 3D fibrous structures and determination of pore-size distributions are described. The flow conservation equations and the corresponding boundary conditions used in the micro-scale simulations are presented later, and finally, the generation of computational grids is discussed. Generation of virtual 3D fibrous structures In this study, virtual fibrous structures with different fiber orientations, namely unidirectional and isotropic structures, were constructed using GeoDict, a voxel-based commercial code (Math2Market GmbH, GeoDict2012R1, www.geodict.com). In GeoDict, the fiber orientation distribution is controlled by a density function p(ψ , ϕ ) defined in polar coordinates as
p(ψ , φ) =
1 4π
β sin ψ , (1 + (β 2 − 1)cos2 ψ)3/2
(3)
where β is an anisotropy parameter, and its value reflects the fiber orientation of the micro-structure generation, ψ ∈ [0,π ) and ϕ ∈ [0,2π ) are through-plane and in-plane angles, respectively (Ashari and Vahedi Tafreshi, 2009). The density function (Eq. 3) depends on through-plane angle only, and the fiber orientation is controlled by the value of β . For an isotropic structure, the value of β is equal to one. As the value of β increases, the fibers tend to become parallel to the in-plane. The degree of fiber alignment is described by the secondorder fiber orientation tensor (Soltani et al., 2014; Stylianopoulos et al., 2008),
=
1 li ltot
⎡
sin ψ cos2 φ × ⎣sin2 ψ sin φ cos φ cos ψ sin ψ cos φ 2
sin ψ sin φ cos φ 2 2 sin ψ sin φ cos ψ sin ψ sin φ 2
⎤
cos ψ sin ψ cos φ cos ψ sin ψ sin φ ⎦, cos2 ψ
(4)
where li is the length of the ith fiber, ltot is the total fiber length, and the sum is over all fibers within the domain. The values of diagonal components in the fiber orientation tensor (Eq. 4) are a measure
50
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
Fig. 1. Isometric view of virtual 3D fibrous structures considered in this study: (A) Unidirectional structure, and (B) Isotropic structure; Histograms of geometric and porosimetric pore-size distribution obtained using Full Morphology (FM) method for (C) Unidirectional structure, and (D) Isotropic structure.
of fiber alignment in the principal directions, while the off-diagonal components indicate alignment in a direction other than a principal direction (Soltani et al., 2014; Stylianopoulos et al., 2008). The trace of is always 1. For an unidirectional structure in which all the fibers are parallel to x-axis, xx = 1; yy = zz = 0, whereas for an isotropic structure, xx = yy = zz = 0.33. Unidirectional and isotropic fibrous structures were constructed using infinitely long cylindrical fibers with a diameter (d) of 50 μm in a domain of size 4 mm × 1 mm × 1 mm. The porosity of the fibrous structures was fixed at 0.4, and a resolution of 5 μm per voxel was used. The unidirectional and isotropic fibrous structures are shown in Fig. 1A and B, respectively. The virtual fibrous structures were exported from GeoDict as stereo-lithography (STL) files to create computational grids needed for the micro-scale simulations. It should be noted that the fibrous materials used in most of the industrial and consumer applications are highly porous (porosity > 0.65) as compared to the fibrous structures considered in this study (Hong and Kim, 2007; Soltani et al., 2014). The fibers were allowed to overlap and penetrate through each other during geometry creation. So, although cylindrical fibers were used to create the structures, the final structures contained objects with irregular shapes. However, it should be stated that each virtual structure (in their final form) has been characterized both by porosimetry (full morphology method) and capillary penetration numerical experiments to study the relationship between effective and representative pore radii. The use of low porous fibrous structures (without √ regularity in fiber shape) may affect or shift the inception of L∼ t behavior
(Schoelkopf et al., 2002), but it will not impact the overall liquid penetration kinetics at the micro-scale. Determination of pore-size distribution of the fibrous structures The full-morphology (FM) method implemented in GeoDict was used to determine pore-size distribution of the fibrous structures (Becker et al., 2007). The FM method assumes that the fibrous medium is connected at one end to a large imaginary reservoir of non-wetting fluid. The non-wetting fluid is forced into the pore-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. 5), can be intruded by the nonwetting phase at a given capillary pressure (Pc ),
Pc =
2σ cosθ . r
(5)
Erosion step: For each pressure increment, a spherical structural element of radius ‘r’, determined using the Young–Laplace equation at a given Pc , is used in a morphological operation of pore space. This erosion operation gives all possible center points in a domain where a given sphere could be fitted without intruding the fiber boundaries. Dilation step: In this step, the pore volume intruded by the nonwetting phase is obtained by dilating the structuring element, and the procedure is repeated by incrementing pressure in predetermined steps (Ashari and Vahedi Tafreshi, 2009; Becker et al., 2007). In the fluid-intrusion process, connectivity of the non-wetting phase to the
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
Fig. 2. (A) Schematic of the unidirectional capillary penetration experiment performed in this study; (B) An example of isometric view of the computational domain with a background hexahedral mesh.
reservoir was ignored, and the pore-size distribution thus obtained is called Geometric pore-size distribution (GPSD). In addition to GPSD, the porosimetric pore-size distribution (PPSD) of fibrous structures was determined by performing a connectivity analysis (before the dilation step) in which all the eroded spaces that are not connected to the non-wetting phase reservoirs were removed (Ashari and Vahedi Tafreshi, 2009; Becker et al., 2007). Fig. 1C and D show histograms of GPSD and PPSD for unidirectional and isotropic structures, respectively. It can be seen that the PPSDs are different from the GPSDs because of the requirement of connectivity of the non-wetting phase to the reservoir (connectivity check) in PPSDs before the dilation step. For a unidirectional structure (Fig. 1C), however, the PPSD and GPSD are similar, as the pore gaps are of constant diameters, and are connected throughout the domain. The PPSD of an isotropic structure has a narrower distribution, which could be attributed to the ink-bottle effect (Jaganathan et al., 2008a). The representative pore radii (R50 - median of the curves) of the fibrous structures were estimated from the GPSD and PPSDs. The R50 values obtained from the pore size distributions were compared with the effective pore radii calculated from the two-phase flow simulations in the results and discussion section. Flow field calculation
Single-phase flow: through-plane permeability calculation The flow field inside the inter-fiber void space of a fibrous structure was obtained by solving the steady, 3D incompressible Navier– Stokes equations (Eqs. 6 and 7) using a finite-volume-based computational fluid dynamics software (simpleFoam solver, OpenFOAM). The governing equations are as follows: Continuity (or mass conservation) equation: (6)
Momentum conservation equation:
U · ∇ U = −∇
P ρ
+ υ∇ 2U,
tions, at the inlet, a constant uniform velocity normal to the boundary was specified (Fig. 2B); The velocity at the inlet was chosen such that the value of Reynolds number, defined based on fiber diameter, was less than unity (Re << 1), resulting in a Stokes flow. At the outlet, a constant gauge pressure boundary condition with a fixed value of 0Pa was used. The symmetry boundary condition (normal velocity component and shear forces set to zero) was applied at the side boundaries, as the lateral flow is negligible as compared to the throughplane flow (Jaganathan et al., 2008b; Tahir and Tafreshi, 2009). The no-slip boundary condition was imposed on the surface of all the fibers. From the micro-scale simulation, the average pressure drop ( P) across the thickness of the fibrous medium was obtained, and the through-plane permeability of the virtual fibrous structure was calculated using Darcy law
V =−
P . μ z k
(8)
Here, k is the permeability constant, μ is the dynamic viscosity of the fluid, V is the velocity magnitude at the inlet boundary, and z is the thickness of the fibrous medium. The calculated permeability values were used to determine static radii (Rs ) of the fibrous structures. Two-phase flow: capillary penetration of a wetting liquid The unidirectional capillary penetration (neglecting gravity) of a wetting liquid through fibrous media was simulated using the VOF method, a fluid fraction-based interface-capturing technique used for modeling the flow of a system of two incompressible, immiscible fluids (Ashish Saha and Mitra, 2009; Berberovic´ et al., 2009). The VOF method treats a two-fluid system as a single-fluid system with spacedependent physical properties. In the conventional VOF method, a transport equation for the fluid fraction (Eq. 11) of one phase is solved simultaneously with the continuity (Eq. 9) and momentum equations (Eq. 10). The governing equations are: Continuity equation (mass conservation):
∇ · U = 0,
(9)
Momentum conservation equation:
This section describes the governing equations and boundary conditions for the fluid flow simulations. From Eq. 2, it can be seen that the permeability values of the fibrous structures are required to determine static radii (Rs ). The permeability values were calculated by performing single-phase flow calculations, and the capillary penetration of a wetting liquid through the fibrous structures was simulated using the VOF method. The flow simulations were performed using OpenFOAM, an open-source CFD code (OpenCFD Ltd., v 2.0, www.openfoam.org). Fig. 2A shows the schematic of a unidirectional capillary penetration numerical experiment performed in this study; an isometric view of a computational domain with background hexahedral mesh is shown in Fig. 2B.
∇ · U = 0,
51
(7)
where P is the pressure, U is the velocity vector, ρ and υ are the density and kinematic viscosity, respectively. For the boundary condi-
∂(ρU ) + ∇ · (ρUU ) = −∇ P + ∇ · μ(∇ U + ∇ U T ) + Fσ , ∂t
(10)
Fluid fraction advection equation:
∂α + ∇ · (U α) = 0. ∂t
(11)
Here, U is the velocity vector, α is the fluid fraction, Fσ is the surface tension force, ρ and μ are density and dynamic viscosity of the fluid, respectively. The fluid fraction (α ) is a characteristic function of one of the two fluids, and can take values within the range of 0 ≤ α ≤ 1; α = 0 corresponds to cells filled with the non-wetting phase, whereas α = 1 indicates cells containing the wetting phase. Gradients of the fluid fraction are encountered only in the region of the interface. The physical properties of the fluid are calculated as weighted averages based on the distribution of the fluid fraction function,
μ = μw α + μnw (1 − α),
(12)
ρ = ρw α + ρnw (1 − α).
(13)
Here, ρ w and ρ nw are the densities of the wetting and non-wetting phases, respectively. Whereas μw and μnw are the viscosities of the wetting and non-wetting phases, respectively. In OpenFOAM, a modified fluid fraction equation (or advection equation) based on a two-fluid approach is solved to ensure boundedness and conservation of the fluid fraction (α ); this is essential for obtaining a bounded solution for flow with fluids having high density
52
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
ratios (Berberovic´ et al., 2009; Deshpande et al., 2012). The modified fluid fraction equation is given as
∂α + ∇ · (U α) + ∇ · (Ur α(1 − α)) = 0.
∂t
(14)
CompressionTerm
Here, Ur is the compression velocity. The additional nonlinear convection term, referred to as the “compression term”, is active only near the interface and does not affect the solution outside the interface region. With an appropriate and a suitable choice of parameters, the compression terms ensures mass conservation of the fluid fraction, and guarantees a bounded and unbiased solution. In the momentum equation (Eq. 10), the surface tension force (Fσ ) at the fluid–fluid interface is computed using the continuum surface force (CSF) method (Brackbill et al., 1992). In this method, the force acting at the interface (which is a surface of discontinuity) is represented as a body or continuum force acting over the volume of the region occupied by the interface (0<α <1), and is formulated as (Hoang et al., 2013; Raeini et al., 2012).
Fσ = σ κ∇α ,
(15)
where ’κ ’ is the mean curvature of the free surface, and is determined from
∇α κ =∇· . |∇α|
(16)
In the CSF–VOF method, the interaction between the fluids and the fiber surface (wall adhesion) is included as a boundary condition on the surface normal vector at the solid wall.
nˆ = nˆ s cos θ + tˆs sin θ
(17)
Here, nˆ is the unit vector normal to the interface at the fiber surface, and nˆ s and tˆs are the unit vectors normal and tangential to the solid boundary, respectively; θ is the equilibrium or static contact angle at the triple-contact line (fluid–fluid interface and the solid wall). Two-phase flow simulations: initial and boundary conditions. The transient two-phase flow computations at micro-scale were performed with water (ρ = 1000 kg/m3 and μ = 1e-3Pa–s) and air (ρ = 1 kg/m3 and μ = 1.8 e–6 Pa–s) as the wetting and non-wetting phases, respectively. The surface tension between the wetting and non-wetting phases was specified to be 0.0725 N/m, and a fixed or static contact angle of 55° was used in the micro-scale simulations. At time t = 0 s, the domain was initialized in a manner that the fibers next to inlet (Fig. 2B) were in contact with the wetting phase (water). In order to simulate a capillary-driven flow, a constant gauge pressure condition (0 Pa) was specified at the domain inlet and outlet with standard velocity inlet/outlet boundary conditions. The fluid fraction value was set to unity at the inlet (wetting-phase reservoir), and a zero-gradient condition for fluid fraction was used at the outlet. The symmetry boundary condition (normal velocity and shear forces set to zero) was applied on the side boundaries of the domain, and the no-slip boundary condition was imposed on the surface of the fibers. Accuracy of the two-phase flow calculations performed using OpenFOAM were verified by studying the capillary transport between two infinite parallel plates for which analytical solution is available in literature.
Fig. 3. (A) Isometric view of a fibrous structure with a slice at the middle of z-axis (zmid ); Computational grid corresponding to 10 μm resolution at zmid for fibrous structures: (B) Unidirectional structure, and (C) Isotropic structure.
with a desired resolution and a cell aspect ratio of unity. 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. A grid convergence study was performed to ensure that the results were mesh-independent. Computational grids with a cell-aspect ratio close to one, which is desired for stable two-phase computations using VOF method, were constructed with five different grid resolutions ( x): 25 μm, 16.33 μm, 12.5 μm, 10 μm, 8 μm. An isometric view of a fibrous structure with a slice at the middle of z-axis (zmid ) is shown in Fig. 3A. Fig. 3B and C show the computational grids corresponding to 10 μm at zmid for the unidirectional and isotropic fibrous structures, respectively. 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. The computations were performed on the Procter and Gamble parallel-computing cluster consisting of 2.93 GHz Intel Xeon processors (X5670) with 2 Gigabytes of memory per processor. Transient two-phase flow calculations corresponding to 10 μm grid resolution for an unidirectional (1,911,145 cells) and isotropic structures (1,915,138 cells) were performed on 120 and 144 processors, respectively. These transient simulations for unidirectional and isotropic fibrous structures consumed 43 h and 186 h of wall clock time to simulate physical times of time 0.016 s and 0.044 s, respectively.
Computational grid generation Results and discussion The computational domains were discretized and hexahedraldominant computational grids were generated from the triangulated surface geometries (STL files) using foamProMesh (iconCFD, www.iconcfd.com). The tortuous fluid paths inside the fibrous structures necessitated the use of unstructured grids. As a first step in the meshing procedure, a background hexahedral mesh was created
The first part of this section deals with the kinetics of unidirectional capillary penetration of a wetting liquid. Next, the effective pore radii of fibrous structures calculated from the micro-scale simulations were compared with the representative radii obtained from the pore-size distribution data.
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
53
Fig. 4. Capillary transport between two parallel plates as a function of time: (A) Capillary penetration distance, and (B) Velocity of the wicking front.
Fig. 5. Effect of grid resolution on (A) Mass uptake per unit cross-sectional area, and (B) Mass flux for an unidirectional fibrous structure. Five different grid resolutions (࢞x = 25 μm, 16.33 μm, 12.5 μm, 10 μm, 8 μm) were considered for this study.
Capillary penetration of a wetting liquid: Lucas–Washburn kinetics Here, we first describe the results of capillary transport between two parallel plates, with the objective to verify our two-phase calculations. Later, the results of the grid-convergence study are explained, following which the mass uptake per unit cross-sectional area and mass flux values for different fibrous structures are presented. Finally, the kinetics of liquid penetration at the micro-scale is discussed in detail. Capillary transport between two parallel plates: verification of two-phase flow calculations The accuracy of the two-phase flow calculations was verified by simulating capillary transport between two infinite parallel plates, separated by a distance of 1 mm. The simulation was performed using an air-water configuration (with water as the wetting liquid); a static contact angle of 70° was fixed in the simulation. At time t = 0 s, the liquid meniscus position in the domain was initialized at 2 mm. At the inlet, the fluid fraction value was set to unity, and a zero-gradient condition for the fluid fraction was used at the outlet. Unidirectional capillary-driven flow was simulated by specifying a constant gauge pressure of 0 Pa at both inlet and outlet. Fig. 4 shows the capillary dynamics (average distance travelled and average velocity of the interface with time) obtained from the OpenFOAM computations along with the results of the Lucas– Washburn equation (Eq. 1), included for comparison. It can be seen (Fig. 4A) that the numerical results follow analytical results of the Lucas–Washburn equation, except during the initial stages which are dominated by inertial forces (Fries and Quéré, 2010; Lu et al., 2013). Fries and Quéré (2010) provided an analytical expression
(t1 = 2.1151 R2 ρ /μ = 1.41 s) to calculate the time required to reach the purely viscous stage, in which the Lucas–Washburn equation is valid. At time instant t1 = 1.41 s, the maximum difference between the numerical data and the results of the Lucas–Washburn equation is 1.43%, thereby indicating good agreement of the numerical data with the analytical results. Fig. 4B shows the average velocity of the wicking front predicted by the OpenFOAM calculations and the Lucas–Washburn equation. Here again, good agreement is observed between the numerical and analytical results in the purely viscous stage (after time t1 s). In the next section, the results of gridconvergence study are presented. Grid convergence study: micro-scale simulations Since two-phase flow computations are sensitive to grid resolution, the capillary penetration of a wetting liquid through an unidirectional fibrous structure was simulated on computational grids with different resolutions to estimate the effects of discretization. The mass uptake and mass flow rate of the wetting liquid were monitored during the simulations. Fig. 5 shows the mass uptake per unit crosssectional area and mass flux of the unidirectional fibrous structure for different grid resolutions as a function of time. It can be seen that, at a given time instant, mass uptake and mass flux decrease with an increase in grid resolution. Minimal variation was seen in the results of the fine grid (12.5 μm), finer or penultimate grid (10 μm), and finest grid (8 μm). At t = 0.015 s, difference in the mass uptake values predicted by the finest grid (8 μm) and finer or penultimate grid (10 μm) was 2%, and the difference in the mass uptake values predicted by finer grid (10 μm) and fine grid (12.5 μm) was 4.4%. Based on these results, the finer grid (10 μm) resolution was deemed sufficient, and was used in the subsequent micro-scale computations. In
54
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
Fig. 6. (A) Mass uptake per unit cross-sectional area, and (B) Mass flux as a function of time for unidirectional and isotropic fibrous structures.
Fig. 7. Fluid fraction contours on a representative slice at the middle of z-axis at t = 0.01 s for (A) Unidirectional fibrous structure, and (B) Isotropic fibrous structure.
the next section, the mass uptake and mass flux results for both the fibrous structures are presented. Capillary penetration of a wetting liquid through 3D fibrous structures: mass uptake and mass flux from micro-scale simulations The cumulative mass uptake per unit cross-sectional area and mass flux as a function of time for unidirectional and isotropic fibrous structures are shown in Fig. 6A and B, respectively. From Fig. 6, it can be seen that the cumulative mass uptake of the fibrous structures increases with time (Fig. 6A), while the mass flux of fibrous structures (Fig. 6B) increases initially, reaches a peak, and then decreases as the viscous effects in the flow become dominant. At any given time instant, as expected, the unidirectional structure has the highest mass uptake and mass flux as compared to an isotropic fibrous structure. Fig. 7A and B show the fluid fraction contours on a representative slice at the middle of z-axis at t = 0.01 s for the unidirectional and isotropic fibrous structures, respectively. It is evident that the wicking front propagates faster in a unidirectional structure as compared to an isotropic structure because of better connectivity between the inter-fiber spaces through which the fluid flows (Ashari et al., 2010; Soltani et al., 2014). Kinetics of liquid penetration at micro-scale The average distance (L) travelled by the wicking front with time was calculated from the results of the micro-scale simulations. Fig. 8 shows the location of the wicking front as a function of the square √ root of time ( t) for the fibrous structures. From Fig. 8, it can be
√ Fig. 8. Location (L) of liquid front as a function of square root of time ( t). A linear √ relationship was observed between L and t during the later stages of capillary transport.
seen that the average distance increases non-linearly during the early stages of liquid transport. However, during the later stages √ of liquid transport, a linear relationship is observed between L and t for both the fibrous structures. From these results, it can concluded that the kinetics of liquid penetration at the micro-scale√follow the behavior delineated by Lucas–Washburn equation (L∼ t), after the early time stages of capillary transport which are dominated by inertial forces (Hyvaluoma et al., 2006; Simile, 2004). Recent experimental works carried out at the nano-scale also support√ the validity of the Lucas–Washburn kinetics, which suggest that L∼ t behavior might be fundamental to the capillary penetration mechanism (Stukan et al., 2010). √ It is interesting to note that a viscous stage (L∼ t) occurs in unidirectional capillary transport through porous media at micro-scale, similar to the purely viscous flow stage observed in the capillary transport between two parallel plates. However, unlike the capillary √ transport in parallel plates or a single capillary, manifestation of L∼ t behavior in a porous medium does not mean that inertia plays no role in later stages. In a porous medium, there is an inertial contribution each time the liquid is accelerated (Schoelkopf et al., 2002). In this respect, the picture of a single capillary representing a porous substrate is misleading (Schoelkopf et al., 2002).
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
55
Table 1 Comparison of effective pore radii (Rc and Rs ) obtained from the micro-scale simulations with the representative radii (R50 ) calculated from the porosimetric (R50, P ) pore-size distributions (PSD). Radii
Pore-size distribution Geometric PSD Porosimetric PSD Two-phase flow simulations Bundle-of-capillaries Darcy law approach
Fibrous structure Unidirectional (k = 1.23 e–10 m2 )
Isotropic (k = 4.2 e–11 m2 )
R50, G R50, P
29.28 μm 29.26 μm
31.59 μm 22.98 μm
Rc Rs
47.47 μm 41.57 μm
10.18 μm 66.1 μm
Characterization of fibrous media: estimation of effective pore radius In the present study, the effective pore radius of a fibrous structure was determined using two approaches: bundle-of-capillaries approach (equivalent capillary radius) and Darcy law approach (static radius). Equivalent capillary radius (Rc ) using √ bundle-of-capillaries approach: From Fig. 8, the slopes of the L versus t curves, referred to as penetration rates or wicking coefficients, were determined for all three fibrous structures. Using the penetration rates, wetting phase properties, and contact angle, the equivalent capillary radii (Rc ) of the fibrous structures were calculated (Eq. 1). Once the effective pore radii were determined, the next step is to determine the static radii of the fibrous structures. Static radii (Rs ) using Darcy law approach: In order to determine Rs of a fibrous structure, we need permeability of the medium (Eq. 2) in addition to penetration rate, fluid properties, and porosity. The permeability values of the fibrous structures were calculated as described in the methodology section. The static radii of the fibrous structures were determined using Eq. 2. Table 1 lists the effective pore radii (Rc & Rs ) obtained from capillary penetration simulations, and the representative pore radii obtained from the GPSD (R50, G ) and PPSD (R50, P ) data. The throughplane permeability values of the fibrous structures are also listed in Table 1. The permeability of the unidirectional fibrous structure is greater than that of the isotropic fibrous structure, as the unidirectional structure offers less resistance to the fluid flow (Soltani et al., 2014; Tahir and Tafreshi, 2009; Wang et al., 2007). For the unidirectional fibrous structure, the R50, G and R50, P have approximately the same value, whereas for the isotropic fibrous structure, R50, P is smaller than R50, G because of the narrower PPSDs. From Eqs. 1 and 2, it can be seen that the average distance (L), and thus the mass uptake, of a fibrous structure is directly proportional to Rc , and inversely proportional to Rs . As expected, the unidirectional fibrous structure had the highest Rc and the lowest Rs . The Rs values for both the fibrous structures were greater than the R50 values obtained from the GPSDs and PPSDs. Next, on comparing R50 with Rc , we see that for a unidirectional structure, the values of R50 are smaller than that of Rc , whereas for the isotropic structure, the values of R50 are larger than Rc . The lack of meaningful correlation between the R50 and Rc values of the different fibrous structures can be attributed to the selective pore-filling that occurs because of the differences in the pore connectivity, which influences the liquid transport through the structure (Schoelkopf et al., 2002). As stated in introduction section, in all the simulations, the fluid properties and contact angle were kept constant, thereby removing any variations that may arise from the variation of contact angle, that are prevalent in physical experiments. Even though the Lucas– Washburn coefficients were kept constant in all the simulations, the Rc and Rs values determined from capillary penetration simulations were found to be different from R50 values obtained from pore-size distribution data (Hanžiˇc et al., 2010). Hence, from this study, it can be stated that the differences between effective and representative
pore radii of a fibrous medium do not arise from the contact angle variations. Conclusions In the present effort, direct numerical simulations were conducted to study the kinetics of liquid penetration at the microscale, and the relationship between the effective pore radii (Rc and Rs ) calculated from capillary penetration experiments and the representative radius (R50 ) estimated from pore-size distribution data (porosimetry). In the literature, for a given porous medium, different values of effective and representative pore radii were obtained, and the differences in pore radii were attributed to the contact angle variations that are prevalent in the experiments (possibly induced by contaminations). In contrast to physical experiments, numerical simulations allowed us to perform capillary penetration numerical “experiments” by keeping the contact angle constant, thus eliminating the effect of contact angle variations. Transient two-phase flow simulations of capillary penetration of a wetting liquid through 3D virtual fibrous media were performed with a constant contact angle using the finite-volume-based volume-of-fluid (VOF) method. The simulated kinetics of liquid penetration at the micro-scale fol√ lowed the Lucas–Washburn equation (L∼ t) during the later stages of the capillary transport. The deviation (non-linear behavior) from Lucas–Washburn kinetics during the early stages could be attributed to the inertial effects that are accounted for in the micro-scale simu√ lations. The inception of a purely viscous stage (L∼ t) suggests that the liquid propagation through fibrous structures during the later stages of capillary transport is controlled primarily by the viscous drag caused by the wetted fibers. Next, the effective pore radii of the fibrous structures were calculated using bundle-of-capillaries (equivalent capillary radius) and Darcy law (static radius) approaches from the rates of liquid penetration obtained from the micro-scale simulations. The equivalent capillary radii and static radii of the fibrous structures were compared with the representative pore radii obtained from the pore-size distributions. Even though the fluid properties and contact angle were kept constant in all the simulations, the effective pore radii (equivalent capillary radii and static radii) obtained from the simulations were found to be different from the representative radii. Thus, a unique value for the characteristic pore radius does not exist for a fibrous medium. The value of characteristic pore radius is strongly dependent on the experimental technique (e.g. porosimetry or capillary penetration experiment) used to characterize the medium. The selection of the characterization method should be based on the end application of the material. Acknowledgment This work is sponsored by the Procter & Gamble company through the University of Cincinnati Simulation Center. The authors appreciate the anonymous referee for his/her insightful comments and suggestions which have significantly improved the manuscript.
56
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57
APPENDIX A. derivation of Lucas–Washburn equation Applying Newton’s second law on a liquid column rising in a circular tube of radius r, and assuming that flow is incompressible, viscous, laminar, and fully developed, one arrives at the 1D infiltration equation of motion (Fries and Quéré, 2010),
2
ρ LL + L
Inertia
2σ cos θ 8 − 2 μLL − r r
=
Capillary Pr essure
V iscousLosses
ρ gL
,
=
(19)
V iscousLosses
Capillary Pr essure
Integrating Eq. 19, one obtains Lucas–Washburn equation which predicts time evolution of length of the liquid column, and is given as
L=
rσ cos θ 1/2 t . 2μ
(20)
PenetrationRate
The Lucas–Washburn equation is also used to describe the kinetics of liquid penetration in a porous medium (Fries and Quéré, 2010; Masoodi, 2010); here, the cylinder radius r in Eq. 20 is replaced by the equivalent capillary radius of the porous medium Rc as shown in Eq. 21.
L=
σ Rc cosθ t. 2μ
(21)
APPENDIX B. derivation of Darcy Law based Lucas–Washburn equation for porous media The momentum equation of a liquid column rising in circular tube (Eq. 18) can be extended to a generic porous medium by replacing the Hagen–Poiseuille viscous losses using the Darcy law, and is given as (Fries and Quéré, 2010; Masoodi, 2010)
2 2σ cos θ ϕ ρ LL + L = − μLL − ρ gL , Rs k
HydrostaticForce Inertia
Capillary Pr essure
(22)
V iscousLosses
where Rs , k, and φ are static radius, permeability and porosity of the medium, respectively. It should be noted that the radius r in Eq. 18, used to calculate the capillary pressure in a circular tube, is replaced by the static radius Rs in Eq. 22 to emphasize that the porous structure does not consist of capillary tubes with a specific radius r; here, one assumes that the porous medium features a uniform pore size (Fries and Quéré, 2010; Masoodi, 2010). In Eq. 22, neglecting the influence of inertial and gravity forces, one obtains the balance between capillary pressure (driving force) and viscous drag,
2σ cos θ Rs
Capillary Pr essure
=
ϕ
μLL . k V iscousLosses
L=
4σ kcosθ μϕ Rs
t.
(24)
References
HydrostaticForce
8 μLL . 2 r
(18)
where L is the distance travelled by the advancing front, θ is the contact angle, σ is the surface tension of fluid–fluid interface, g is gravity, ρ and μ are density and viscosity of the wetting fluid, respectively. In Eq. 18, the capillary pressure (driving force) is balanced by the sum of inertial, viscous, and hydrostatic forces. Lucas and Washburn investigated the wicking dynamics in a circular tube neglecting the influence of inertial and gravity forces; here, the capillary pressure (driving force) is balanced by the viscous drag.
2σ cos θ r
Integrating Eq. 23, one obtains modified Lucas–Washburn equation which predicts time evolution of length of the liquid column in a porous medium (Fries and Quéré, 2010; Masoodi, 2010)
(23)
Ashari, A., Bucher, T., Tafreshi, H.V., Tahir, M., Rahman, M., 2010. Modeling fluid spread in thin fibrous sheets: effects of fiber orientation. Int. J. Heat Mass Trans. 53, 1750– 1758. Ashari, A., Vahedi Tafreshi, H., 2009. A two-scale modeling of motion-induced fluid release from thin fibrous porous media. Chem. Eng. Sci. 64, 2067–2075. Ashish Saha, A., Mitra, S.K., 2009. Effect of dynamic contact angle in a volume of fluid (VOF) model for a microfluidic capillary flow. J. Colloid Int. Sci. 339, 461–480. Becker, J., Wiegmann, A., Schulz, V., 2007. Design of fibrous filter media based on the simulation of pore size measures. In: Proceedings of the FILTECH 2007 conference. ´ E., van Hinsberg, N.P., Jakirlic, ´ S., Roisman, I.V., Tropea, C. , 2009. Drop imBerberovic, pact onto a liquid layer of finite thickness: dynamics of the cavity evolution. Phys. Rev. E 79, 036306. Bousfield, D.W., Karles, G., 2004. Penetration into three-dimensional complex porous structures. J. Coll. Int. Sci. 270, 396–405. Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surfacetension. J Comput. Phys. 100, 335–354. Deshpande, S.S., Anumolu, L., Trujillo, M.F. , 2012. Evaluating the performance of the two-phase flow solver interFoam. Comput. Sci. Discov. 5, 014016. Ferrari, A., Lunati, I., 2013. Direct numerical simulations of interface dynamics to link capillary pressure and total surface energy. Adv. Water Res. 57, 19–31. Fries, N., Quéré, D., 2010. Capillary Transport Processes in Porous Materials-Experiment and Model. Cuvillier Verlag Göttingen Doctoral dissertation2010; ISBN 978-386955-507-2. Gane, P.A., Ridgway, C.J., Schoelkopf, J., 2004. Absorption rate and volume dependency on the complexity of porous network structures. Transp. Porous Media 54, 79–106. Hamraoui, A., Nylander, T., 2002. Analytical approach for the Lucas–Washburn equation. J. Coll. Interface Sci. 250, 415–421. Hanžiˇc, L., Kosec, L., Anžel, I., 2010. Capillary absorption in concrete and the Lucas– Washburn equation. Cement Concr. Compos. 32, 84–91. Hoang, D.A., van Steijn, V., Portela, L.M., Kreutzer, M.T., Kleijn, C.R., 2013. Benchmark numerical simulations of segmented two-phase flows in microchannels using the Volume of Fluid method. Comput. Fluids 86, 28–36. Hong, C.J., Kim, J.B., 2007. A study of comfort performance in cotton and polyester blended fabrics. I. vertical wicking behavior. Fibers Polym. 8, 218–224. Hyvaluoma, J., Raiskinmaki, P., Jasberg, A., Koponen, A., Kataja, M., Timonen, J. , 2006. Simulation of liquid penetration in paper. Phys. Rev. E Stat. Nonlin. Soft. Matt. Phys. 73, 036705. Jaganathan, S., Vahedi Tafreshi, H., Pourdeyhimi, B., 2008a. Modeling liquid porosimetry in modeled and imaged 3-D fibrous microstructures. J. Coll. Int. Sci. 326, 166– 175. Jaganathan, S., Vahedi Tafreshi, H., Pourdeyhimi, B., 2008b. A realistic approach for modeling permeability of fibrous media: 3-D imaging coupled with CFD simulation. Chem. Eng. Sci. 63, 244–252. Lavi, B., Marmur, A., Bachmann, J., 2008. Porous media characterization by the twoliquid method: effect of dynamic contact angle and inertia. Langmuir 24, 1918– 1923. Li, Z., Giese, R.F., Oss, C.J., Kerch, H.M., Burdette, H.E., 1994. Wicking technique for determination of pore size in ceramic materials. J. American Ceram. Soc. 77, 2220–2222. Lu, G., Wang, X.-D., Duan, Y.-Y., 2013. Study on initial stage of capillary rise dynamics. Coll. Surf. A: physicochem. Eng. Asp. 433, 95–103. Marmur, A., 2003. Kinetics of penetration into uniform porous media: testing the equivalent-capillary concept. Langmuir 19, 5956–5959. Masoodi, R., 2010. Modeling Imbibition of Liquids into Rigid and Swelling Porous Media. The University of Wisconsin–Milwaukee Doctoral dissertation. O’Loughlin, M., Wilk, K., Priest, C., Ralston, J., Popescu, M., 2013. Capillary rise dynamics of aqueous glycerol solutions in glass capillaries: a critical examination of the Washburn equation. J. Coll. Interface Sci. 411, 257–264. Pradhan, A.K., Das, D., Chattopadhyay, R., Singh, S., 2012. Effect of 3D fiber orientation distribution on transverse air permeability of fibrous porous media. Powder Tech. 221, 101–104. Raeini, A.Q., Blunt, M.J., Bijeljic, B., 2012. Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method. J. Comput. Phys. 231, 5653–5668. Ridgway, C.J., Gane, P.A., 2002. Dynamic absorption into simulated porous structures. Coll. Surf. A: physicochem. Eng. Aspects 206, 217–239. Schoelkopf, J., Gane, P.A., Ridgway, C.J., Matthews, G.P., 2002. Practical observation of deviation from Lucas–Washburn scaling in porous media. Coll. Surf. A: physicochem. Eng. Aspects 206, 445–454. Simile, C.B., 2004. Critical Evaluation of Wicking in Performance Fabrics. Georgia Institute of Technology Doctoral dissertation. Soltani, P., Johari, M.S., Zarrebini, M., 2014. Effect of 3D fiber orientation on permeability of realistic fibrous porous networks. Powder Tech. 254, 44–56. Stukan, M.R., Ligneul, P., Crawshaw, J.P., Boek, E.S., 2010. Spontaneous imbibition in nanopores of different roughness and wettability. Langmuir 26, 13342–13352.
N.K. Palakurthi et al. / International Journal of Multiphase Flow 77 (2015) 48–57 Stylianopoulos, T., Yeckel, A., Derby, J.J., Luo, X.-J., Shephard, M.S., Sander, E.A., Barocas, V.H. , 2008. Permeability calculations in three-dimensional isotropic and oriented fiber networks. Phys. Fluids 20, 123601. Tahir, M.A., Tafreshi, H.V., 2009. Influence of fiber orientation on the transverse permeability of fibrous media. Phys. Fluids 21, 083604–083605. Wang, Q., Maze, B., Tafreshi, H.V., Pourdeyhimi, B., 2007. Simulating through-plane permeability of fibrous materials with different fiber lengths. Modell. Simul. Mater. Sci. Eng. 15, 855–868.
57
Xue, H., Fang, Z., Yang, Y., Huang, J., Zhou, L., 2006. Contact angle determined by spontaneous dynamic capillary rises with hydrostatic effects: experiment and theory. Chem. Phys. Lett. 432, 326–330. Zhmud, B., Tiberg, F., Hallstensson, K., 2000. Dynamics of capillary rise. J. Coll. Interface Sci. 228, 263–269.