Technical Note pubs.acs.org/ac
Time-Dependent Model for Fluid Flow in Porous Materials with Multiple Pore Sizes Brian M. Cummins,† Rukesh Chinthapatla, Frances S. Ligler, and Glenn M. Walker* Department of Biomedical Engineering, University of North Carolina at Chapel Hill and North Carolina State University, Raleigh, North Carolina 27695, United States S Supporting Information *
ABSTRACT: An understanding of fluid transport through porous materials is critical for the development of lateral flow assays and analytical devices based on paper microfluidics. Models of fluid transport within porous materials often assume a single capillary pressure and permeability value for the material, implying that the material comprises a single pore size and that the porous material is fully saturated behind the visible wetted front. As a result, current models can lead to inaccuracies when modeling transport over long distances and/or times. A new transport model is presented that incorporates a range of pore sizes to more accurately predict the capillary transport of fluid in porous materials. The model effectively predicts the time-dependent saturation of rectangular strips of Whatman filter no. 1 paper using the manufacturer’s data, published pore-size distribution measurements, and the fluid’s properties.
P
The Lucas-Washburn model also assumes complete saturation behind the wetted front. However, in many of the porous materials used for analytical applications there is a broad pore size distribution, and others have proven that the porous material behind the visible wetted front is not fully saturated.22,23 Modeling volumetric flow rates of liquid into paper with the Lucas-Washburn equation leads to errors between the predicted absorbed volume and measured absorbed volume. Bico and Quere introduced a simple two pore-size model to explain this discrepancy, but porous media typically contains a large range of pore sizes.23 Here, we extend previous modeling work by developing and experimentally validating a model that predicts the volume of liquid imbibed by a porous media with an arbitrary number of pore sizes. The model can be adapted to incorporate resistances upstream of the porous material, porous materials with varying cross-sectional areas, and even varying porous materials placed in series. In particular, the model will be useful for predicting fluid transport within porous media over long distances or transport times (e.g., paper microfluidics) or in cases where multiple porous materials are put in series (e.g., lateral flow tests or paper pumps attached to microfluidic channels).
orous materials are found in a wide variety of products, ranging from moisture-wicking clothing to tea biscuits.1 These materials are particularly important to the analytical chemistry field, where their ability to wick fluid forms the basis of most point-of-use analytical devices including lateral flow tests and paper microfluidics.2−9 Paper, in particular, has received considerable attention recently as an attractive substrate for creating analytical devices because of its innate ability to control fluid flow.10−16 The initial wetting process, or imbibition, into these materials is key to their functionality, and a better understanding of the process through improved mathematical models will yield new applications. The development of computational models for fluid imbibition in porous media is an active area of research.17,18 While these models provide insight into imbibition at micrometer length scales, they are not well-suited for modeling the meso- and macroscale volumetric flow rates associated with paper-based microfluidic pumps19,20 and paper fluidic devices. Instead, the analytical Lucas-Washburn model (eq 1) is often used to relate the position of the visible wetted front, L, to the square root of time, t.21 L2 =
γ cos(θ)r t 2μ
(1)
The Lucas-Washburn model assumes the porous media is a bundle of parallel capillary tubes that have the same radius, r, and the liquid has surface tension γ, contact angle θ, and dynamic viscosity μ (Figure S-1). These values can only be identified empirically by fitting experimental data to eq 1. © XXXX American Chemical Society
Received: November 28, 2016 Accepted: March 15, 2017
A
DOI: 10.1021/acs.analchem.6b04717 Anal. Chem. XXXX, XXX, XXX−XXX
Technical Note
Analytical Chemistry
■
The pressure difference, ΔP, is generated by the capillary pressure, Pc, which drives fluid into a porous material and is determined by the surface tension of the fluid, the contact angle of the fluid with the material, and the inner radius of the capillaries or pores. This negative pressure pulls the fluid toward unwetted regions of the paper and can be calculated with the equation
MODEL FOR TRANSPORT THROUGH POROUS MATERIAL Figure 1A illustrates the wetting of a strip of paper as modeled by Darcy’s law. Capillaries of uniform pore size are filled with
Pc =
2γ cos(θ ) r
(4)
By setting Q = φA dL/dt, combining eqs 2−4, and integrating, a form of the Lucas-Washburn equation results (eq 5). When C = φ/8, eq 5 equals the common form of the Lucas-Washburn equation (eq 1). L2 =
liquid at each cross-section of the wetted paper and the wetted front exists at position L1(t). Darcy’s law assumes uniform pore size and complete wetting behind the liquid front. Darcy’s law can be written24,25 κAΔP μL
(2)
where Q is the volumetric flow rate of a fluid through a porous material, A is the cross-sectional area of all pores within the porous material, and ΔP is the pressure drop across the porous material. The permeability, κ, is a measure of how easily fluid can flow through a given porous material and can be approximated with the equation κ = Cr 2
(5)
Figure 1B shows how liquid wets paper in the new model. In this example, three pore sizes exist within the paper, each with its own wetted front position: L1(t), L2(t), and L3(t). However, any number of pore sizes can be used. At the first cross-section, behind the first wetted front, all pores are wetted. At the second cross-section, behind the second wetted front, only the smaller two pore sizes are wetted and so on. Each pore size accounts for a certain percentage of the total porosity of the material, and the pore size distribution is assumed to be constant throughout the entire porous material. Figure 1C shows the fluidic circuit model of the simplified porous material with the three pore sizes: 1 (small), 2 (medium), and 3 (large). Each pore size has its own wetted front position, Li(t). Since smaller pores pull fluid from larger pores due to their more-negative capillary pressures, the fluid front within smaller pore sizes advances ahead of the fluid front within larger pore sizes (i.e., L1(t) > L3(t)). Voltage sources model the capillary pressures generated at the boundary of the fluid front of a pore size, (e.g., Pc,1 represents the capillary pressure of the smallest pore size group at its fluid front, L1(t)). This model assumes that there is no pressure gradient across the thickness or width of the porous material. Each resistor, Ri,j(t), is the hydraulic resistance due to pore size i between two wetted fronts of different pore size, j. For example, the resistor R1,2(t) is the resistance of the smallest pores between the wetted fronts of the medium and large pores. Because each resistor is a function of its wetted length, the resistances will increase over time as fluid flows into the porous material. In this model, the dimensionless permeability coefficient, C, (eq 3) is assumed to be the same for all pore sizes. Each pore size is assumed to be fully saturated behind its wetted front and completely empty at positions ahead of its wetted front. As a result, Darcy’s law can be applied independently to each pore size. The length from the fluid source to the wetted front of each pore size group can then be calculated as a function of time by solving the circuit diagram of Figure 1C. The volumetric flow rates of fluid into the dry regions of the small, medium, and large pores are Q1, Q2, and Q3, respectively. The velocity of each wetted front position Li(t) is the volumetric flow rate divided by the pore size’s cross-sectional void area (Ai = A × Pcti), where Pcti is the percentage of the total void volume that the pore size represents. The flow rates can be calculated with the equations
Figure 1. Schematic of fluid transport through a rectangular strip of porous material. Planes indicate positions of cross-sectional areas during transport of a perfectly wetting fluid in a porous material. (A) Cross-sectional areas for a strip of material with a single effective pore size show complete saturation behind the visible wetted front, L1. (B) The cross-sectional areas of a material with three different pore sizes show a position-dependent saturation behind the visible wetted front, L1. The wetted front positions for the two larger pore sizes, L2 and L3, are also shown. (C) Circuit model of fluid flow in a porous material with three pore sizes. Qt is the total volumetric flow rate into the porous material. Q1, Q2, and Q3 are the volumetric flow rates into the small, medium, and large pores, respectively.
Q=
4Cγ cos(θ )r t μφ
(3)
where C is a dimensionless permeability coefficient of the material. This coefficient accounts for the tortuosity of the flow path, the porosity of the material, φ, and the shapes of the pores. Porosity is the ratio of the porous material’s void volume to its total volume. B
DOI: 10.1021/acs.analchem.6b04717 Anal. Chem. XXXX, XXX, XXX−XXX
Technical Note
Analytical Chemistry Q1 =
dL1 P − P2 A1φ = 1 dt R1,1
(6)
Q2 =
P2 − P3 dL 2 − Q1 A 2φ = 2 dt (∑1 R i ,2−1)−1
(7)
Q3 =
dL 3 P3 A3φ = − 3 dt (∑1 R i ,3−1)−1
accounted for the hydrostatic pressure in addition to the capillary pressures of each pore size group. Constants used in the model are listed in Table S-1. Transport and Saturation Study. After determining θ and C, horizontal imbibition studies were performed to validate the predictions of the MATLAB model. Rectangular strips (labeled A−G) of Whatman no. 1 filter paper (2 cm × 21 cm) were laminated in pouches, leaving the last 2 cm open to the atmosphere. After lamination, the last 1 cm of laminated paper was cut from each strip to allow trapped air to escape during imbibition. Graduations were marked in 1 cm increments from the start of the lamination, producing a total of 18 segments. Strips B−G were placed on a horizontal surface with the protruding paper tucked behind a 90° stainless steel corner brace in an empty trough. The trough was then filled with dyed water to a height that made the length of paper between the top of the fluid and the start of lamination 1 cm. A plastic sheet was placed on top of the setup to minimize evaporation, and a timelapse video (Movie S-1) was taken from above to monitor the wicking into the paper strips. Strip A served as a 0 min strip (i.e., not dipped in water). Strips B, C, D, E, F, and G were removed at 5, 10, 20, 40, 80, and 160 min, respectively. Immediately after removing a strip, it was cut into segments at the graduated markings along the length of the strip and the segments were weighed to determine the weight of the segment at that time point. Strip A served as a control strip for 0 min and was weighed in the same manner after other strips were finished. The actual dimensions of the segments were later measured with a micrometer. Using this data, the prewetted weight and void volume of each segment were calculated. The difference in the measured weight of a wetted segment and the calculated weight of that segment was assumed to be due to the water that had wetted that segment. The saturation of a segment was determined by comparing the corresponding volume of this amount of water to the expected void volume of that segment.
2
∑ Qi i=1
(8)
For the material in Figure 1C, eqs 6−8 can be solved numerically to track the positions of the wetted fronts for each pore size: L1(t), L2(t), and L3(t). The visible wetted front that is commonly reported for porous materials when using the LucasWashburn equation is L1(t). Once the wetted front of a pore size group reaches the end of the porous material, there is no additional volumetric capacity for that pore size. As a result, the pressure at this position quickly equilibrates with the capillary pressure of next larger pore size and the voltage source for the filled pore size disappears from the circuit. The percent saturation of the entire paper strip can be calculated by adding up the volumes contained in each pore group at any given time. Note that the model predicts eq 5 (Lucas-Washburn flow) if a single pore size is used (i.e., r1 = r and Pct1 = 1.0). The circuit in Figure 1C also predicts that the total flow rate into the porous material, QT, can be solved by knowing the position of the wetted front for the group with the largest pores (eq 9) because the porous material is fully saturated behind this position. In practice, it is difficult to measure the position where full saturation exists because the last liquid front is not visible. However, this feature of the model could be useful in identifying properties of the material and/or have implications for the design of certain materials that control flow-rate. 3
QT =
■
∑ Qi = i=1
P3 3 (∑1 R i ,3−1)−1
■
(9)
RESULTS AND DISCUSSION Vertical Imbibition. While it is possible to calculate the contact angle, θ, by waiting for the fluid to reach its final height, it would be difficult to ensure equilibration, and the process could take weeks. Instead, the height of the visible wetted front was measured at various times over several days, and this data was compared with predictions from the numerical models that were generated by iteratively choosing values for θ and C. Evaporation and the changes in thickness due to the overlap of strips were considered to have a negligible effect on the timedependent position of the front and were not included in the model. An error was calculated for each modeled curve by comparing it to the experimental data, and a least-squares approach was used to identify the θ and C combination, 75° and 0.0063, respectively, that generated a modeled curve that least differed from the experimental curve (Figure 2). This calculated contact angle is close to the contact angle measured by Hong and Kim using a similar setup (θ = 83°).27 The constants were then used to predict the time-dependent saturation of horizontal strips of paper. The dynamic hydrostatic pressure in the vertical orientation makes it possible to identify θ and C. As the wetted fronts of each pore size climb, the driving pressures between fronts change due to the increasing hydrostatic pressures. For a fluid with a known density, the hydrostatic pressure provides a known value and allows θ and C to be estimated via curve
EXPERIMENTAL METHODS Pore Size Calculation. Using the pore size distribution data for Whatman no. 1 filter paper, obtained by mercury porosimetry and corrected for compression (Figure S-2A),26 we extracted ten groups of pore radii linearly spaced from 1 to 19 μm (Figure S-2B). Interpolation was used to identify the percentage of total pore volume that each of these groups represented. Vertical Imbibition. Accurate prediction of the timedependent saturation of a rectangular strip of a porous material requires knowing the material’s porosity, width, thickness, length, dimensionless permeability coefficient, pore size distribution and the fluid viscosity, density, surface tension, and contact angle with the matrix. For Whatman no. 1 paper, all variables are known except the contact angle with the porous matrix, θ, and the permeability coefficient, C. By measuring the vertical height of a fluid column in Whatman no. 1 paper over time and calculating the pressure required to support it, enough information can be obtained to solve for θ and C. Equations 6−8 were expanded to model 10 pore size groups to account for the pore size distribution shown in Figure S-2. The equations were solved in MATLAB to predict the vertical imbibition (i.e., height vs time) of water into the rectangular strip, producing estimates for θ and C. The MATLAB model C
DOI: 10.1021/acs.analchem.6b04717 Anal. Chem. XXXX, XXX, XXX−XXX
Technical Note
Analytical Chemistry
Figure 2. Schematic of the experimental setup for the vertical imbibition study (inset). The experimentally determined position of the visible wetted front in a 2 cm × 1.2 m rectangular strip of laminated Whatman no. 1 filter paper (red squares) was fit to a curve using the values θ = 75° and C = 0.0063 (black curve) which has an R2 = 0.99. The zero position indicates the location of the fluid source.
fitting. There is no such unique solution to imbibition in the horizontal plane because all combinations of θ and C that have the same cos(θ)C value predict the identical imbibition. Only by incorporating a known hydrostatic pressure can individual θ and C values be found. The drawbacks of the vertical strip study are the length of the study time (days to weeks) and the amount of material required for the study. The hydrostatic pressure of the fluid in the paper must be on the order of the capillary pressure of the paper or there will be multiple combinations of θ and C that predict a curve similar to the experimental data. In addition, many specialty papers are sold in relatively small pieces which makes assembling ∼1 m lengths a challenge. Transport and Saturation Study. The time-dependent position of the visible wetted front was measured from timelapse video (Movie S-1) as it moved through the horizontal strips (Figure 3A). Using the θ and C that were identified in the vertical imbibition study and the discretized pore size distribution (Figure S-2B) in the 10-pore size MATLAB model, the model accurately predicts the position of the visible wetted front (Figure 3A). Figure 3A can be fit with a Lucas-Washburn-like curve (L2 = Bt, B = 3.53 × 10−6, R2 = 1). However, the coefficient B assumes a uniform pore size within the paper and cannot be predicted from material properties. In contrast, eq 6 assumes multiple pore sizes within the paper and can be used to predict the position of the visible wetted front using known quantities. The saturation profiles (percent saturation vs position) of paper strips filling with water are shown in Figure 3B. Strips B− G, sacrificed at 5, 10, 20, 40, 80, and 160 min, respectively, represent a snapshot of the saturation profile as water imbibes into the paper. The curves are the predicted saturation from the 10-pore size MATLAB model. Experimental results from the horizontal study show that the porous material is not fully saturated behind the visible wetted front (L1). At every time point, the position of the wetted front of the smallest pore size, L1, leads while the position of the wetted front of the largest pore size, L10, remains closest to the fluid reservoir. These relative positions of wetted fronts L1−L10 remain the same until L1 reaches the edge of the paper, and L2 becomes the leading wetted front. During imbibition, a saturation gradient always exists from the fluid source (100%
Figure 3. (A) Position of the visible wetted front in a horizontal imbibition study using the modeled parameters (black line) and the experimentally determined position of the visible wetted front (points). The additional lines show the predicted positions of the wetted fronts, Li, of other pore size groups. (B) Time-dependent saturation of paper strips. Computed results are shown by solid lines and experimental data are shown by points. Position x = 0 mm is the location of the fluid source.
at x = 0) to the position of the leading wetted front (0% at x = Li). The multiple pore size model described here effectively predicts the characteristics of the saturation profile and closely matches the measured saturation gradients at the time points tested (Figure 3B). Each modeled saturation profile comprises 10 vertical line segments, and the x-axis position of each segment is the predicted wetted front position of that pore size. The vertical length of each line segment is fixed and is proportional to the fraction of total porosity of that pore size. For continuous pore size distributions (e.g., in filter paper), the number of vertical line segments approaches infinity and the modeled curve becomes smooth. Because the relative horizontal positions of the vertical line segments within the saturation profile stay constant during imbibition, it should be possible to predict the volume of fluid that has been absorbed by simply measuring the position of the visible wetted front. Compared to the common assumption that 100% of the void volume is filled behind the wetted front, experimental results showed an average of 87% of the void volume is filled in Whatman no. 1 paper when the wetted front reaches the terminal edge of the paper (Figure S-3). The saturation is less than 100% because only the smallest pore size is completely saturated, while the larger pore sizes have yet to fill. A theoretical paper with a monodisperse pore size would D
DOI: 10.1021/acs.analchem.6b04717 Anal. Chem. XXXX, XXX, XXX−XXX
Technical Note
Analytical Chemistry Present Address
have 100% saturation. The MATLAB model predicted that 80% of the void volume behind the visible wetted front should be filled, showing good agreement with experimental results. In contrast, the Lucas-Washburn model predicts 100% saturation. While the multiple pore size model is a better predictor of the saturation profile than assuming the front is fully saturated behind the visible wetted front, there are slight differences between the experimental data and the model. These differences can be explained by experimental error in the cutting and weighing process and by continued movement of fluid within the strip in the short time between removing it from the fluid source and cutting it into segments. Incorrect assumptions about the permeability (i.e., the use of a constant C for all pore sizes), the lack of a pressure gradient across the thickness of the material, and the discretized pore size distribution could contribute to differences between the model and measured data as well. The model also assumes that the contact angle and pore size distribution is constant for the duration of the experiment and that the pore volume is completely void at the start of the experiment.
†
B.M.C.: Abbott Laboratories, Diagnostics Division, 1921 Hurd Drive, Irving, Texas 75038, United States. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported in part by the Chancellor’s Innovation Fund and funds from the Ross Lampe Distinguished Chair at North Carolina State University.
■
■
CONCLUSIONS The model presented here uses 10 pore sizes to predict timedependent fluid transport into Whatman no. 1 filter paper, although any porous material and any number of pore sizes can be modeled. Experiments prove that a heterogeneous porous material is neither fully saturated nor constantly saturated behind the visible wetted front. By incorporating the material properties and the pore size distribution into the model, the time-dependent position of the visible wetted front, the saturation gradient behind the visible wetted front, and the total volume absorbed can be predicted. Information gleaned from the model can translate measurable, physical properties of porous materials into predictable behaviors for fluid transport, expediting the development of lateral flow assays and analytical devices based on paper microfluidics. While the model developed here has immediate applications to the field of paper-fluidics and lateral flow assays, we envision it ultimately being applied to any field where porous materials are of interest, such as clothing, consumer products, and soil science.
■
ASSOCIATED CONTENT
* Supporting Information S
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.analchem.6b04717. Original data for pore size distribution and extracted size distribution, predicted and measured front positions for vertical study, and table of constants used in MATLAB model (PDF) Time lapsed video of the horizontal imbibition study with the six identical strips (MOV)
■
REFERENCES
(1) Fisher, L. Nature 1999, 397, 469. (2) Martinez, A. W.; Phillips, S. T.; Butte, M. J.; Whitesides, G. M. Angew. Chem., Int. Ed. 2007, 46, 1318−1320. (3) Byrnes, S.; Thiessen, G.; Fu, E. Bioanalysis 2013, 5, 2821−2836. (4) Yamada, K.; Henares, T. G.; Suzuki, K.; Citterio, D. Angew. Chem., Int. Ed. 2015, 54, 5294−5310. (5) Parolo, C.; Merkoci, A. Chem. Soc. Rev. 2013, 42, 450−457. (6) Yetisen, A. K.; Akram, M. S.; Lowe, C. R. Lab Chip 2013, 13, 2210−2251. (7) Nash, M. A.; Hoffman, J. M.; Stevens, D. Y.; Hoffman, A. S.; Stayton, P. S.; Yager, P. Lab Chip 2010, 10, 2279−2282. (8) Ballerini, D. R.; Li, X.; Shen, W. Microfluid. Nanofluid. 2012, 13, 769−787. (9) Li, X.; Ballerini, D. R.; Shen, W. Biomicrofluidics 2012, 6, 011301. (10) Zhong, Z. W.; Wang, Z. P.; Huang, G. X. D. Microsyst. Technol. 2012, 18, 649−659. (11) Koo, C. K. W.; He, F.; Nugen, S. R. Analyst 2013, 138, 4998− 5004. (12) Lu, Y.; Shi, W. W.; Qin, J. H.; Lin, B. C. Anal. Chem. 2010, 82, 329−335. (13) Liu, H.; Crooks, R. M. J. Am. Chem. Soc. 2011, 133, 17564− 17566. (14) Dornelas, K. L.; Dossi, N.; Piccin, E. Anal. Chim. Acta 2015, 858, 82−90. (15) Cassano, C. L.; Fan, Z. H. Microfluid. Nanofluid. 2013, 15, 173− 181. (16) Carrilho, E.; Martinez, A. W.; Whitesides, G. M. Anal. Chem. 2009, 81, 7091−7095. (17) Hyväluoma, J.; Raiskinmäki, P.; Jäsberg, A.; Koponen, A.; Kataja, M.; Timonen, J. Phys. Rev. E 2006, 73, 036705. (18) Jaganathan, S.; Tafreshi, H.; Pourdeyhimi, B. J. Colloid Interface Sci. 2008, 326, 166−175. (19) Cummins, B.; Chinthapatla, R.; Lenin, B.; Ligler, F.; Walker, G. Technology 2017, 1. (20) Wang, X.; Hagen, J. A.; Papautsky, I. Biomicrofluidics 2013, 7, 014107. (21) Washburn, E. W. Phys. Rev. 1921, 17, 273−283. (22) Williams, R. J. Colloid Interface Sci. 1981, 79, 287−288. (23) Bico, J.; Quéré, D. Europhys. Lett. 2003, 61, 348−353. (24) Whitaker, S. Transp. Porous Media 1986, 1, 3−25. (25) Kauffman, P.; Fu, E.; Lutz, B.; Yager, P. Lab Chip 2010, 10, 2614−2617. (26) Gribble, C. M.; Matthews, G. P.; Laudone, G. M.; Turner, A.; Ridgway, C. J.; Schoelkopf, J.; Gane, P. A. C. Chem. Eng. Sci. 2011, 66, 3701−3709. (27) Hong, S.; Kim, W. Microfluid. Nanofluid. 2015, 19, 845−853.
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Phone: 919-513-4390. Fax: 919513-3814. ORCID
Glenn M. Walker: 0000-0003-1980-7990 E
DOI: 10.1021/acs.analchem.6b04717 Anal. Chem. XXXX, XXX, XXX−XXX
