SEPTEhlBER

2001

813

Modeling Capillary Flow in Complex Geometries DILIPRAJAGOPALAN AND ARUNP. ANEJA E. I. di Poiit de Neriiotirs and Coiiipaiiy, IVibiiiiigtoii, Dela,r.nre 19880, U.S.A.

JEAN-MARIE MARCHAL Flrreiit hie., B-1348 Loiivaiii-In-Nernte, Belgiirin ABSTRACT We discuss mathematical models of capillary flow in complex geometries representative of the void spaces formed between fibers in a textile yarn. Moisture transport in textile yarns and fabrics is an important factor affecting physiological comfort. We have extended an existing analytical model for capillary flow in circular tubes to more complex geometries. We validate this model using detailed computational fluid dynamics simulations of this flow. These models are used to understand the effect of geometric and material parameters on moisture transport. In vertical wicking in a bundle of filaments, the model predicts that as the nonroundness of the filaments increases, or the void area between the filaments decreases, the maximum liquid height increases while the initial rate of penetration decreases.

Moisture transport in textile fabrics 'is one of the critical factors affecting physiological comfort. Fabrics that rapidly transport moisture away from the human body make wearers feel more comfortable by keeping them dry. This enhanced moisture transport may also help wearers feel cool as the body provides latent heat to evaporate sweat at an enhanced rate. The comfort afforded by textile fabrics can be improved by understanding the key geometric and material parameters that contribute to moisture transport. Mathematical modeling of surface-tension-driven flow in yarns and fabrics could provide a way to develop such an understanding. Capillary flow in yarns and fabrics has been extensively studied (e.g., [5-7]), and the subject was reviewed by Kissa [9]. Fabrics are typically constructed by knitting or weaving textile yarns, which are essentially bundles of several fibers or filaments. For movement of liquid in a fabric, the liquid must wet the fabric surface before being transported by capillary action through the fabric pores formed between fibers and yarns. This capillary action is determined by the interaction of the liquid and the fabric material, by liquid properties such as viscosity and surface tension, and by the geometric structure of the pores. The size and shape of the fibers, as well as their alignment, determine the geometry of the void spaces o r pores through which the liquid is transported. However, the complexity of a fabric structure makes it impossible to predict pore structure, and very difficult to arrive at a Texrile Rex 1. 71(9), 813-821 (2001)

detailed structure experimentally. Furthermore, movement of liquid through the pores can cause shifting of fibers and changes in the pore structure. For some materials, the fibers can absorb liquid and swell considerably, thereby changing the pore structure even more. Thus, a detailed mathematical o r computational model of capillary flow in fabric structures is not a reasonable goal. The most important phenomenon in fabric wicking is the motion of liquid in the void spaces between fibers in a yarn [6]. The larger pores between yarns are not as important in long range motion of liquid. Thus, a study of wicking in yarns should provide a way to understand the role of geometric and material parameters in fabric wicking. Mathematical models of surface-tension-driven flow in cylindrical capillaries have long been used to interpret results of wicking in yarns [ 5 . 6, 81. Modeling in this area dates back to the work of Lucas [ 131 and Washburn [IS], who independently derived the result that the extent of liquid penetration into a cylindrical capillary is proportional to the square root of time. This resu1t;typically referred to as Lucas-Washburn kinetics, is derived by balancing viscous and capillary forces, neglecting gravitational and inertial forces. Hence, this cequation is strictly applicable to capillary flow in horizontal geometries. The early modeling work on this problem is reviewed by Levine et al. [ 121. They also proposed their own model of this flow, which accounts for entrance effects of the 0010-5175/$2.00

TEXTILE RESEARCH JOURNAL

814 liquid into the capillary, and two-dimensional flow effects near the advancing meniscus. They concluded that neither effect is very important, provided the capillary size is 6 (lo-‘) m or less and the penetration length is larger than the diameter of the capillary. More recently, Reed and Wilson [16] derived a mathematical model of liquid rise in vertical cylindrical capillaries using a macroscopic force balance approach. The predictions of their model are identical to those of Levine et al. [12] once the penetration length into the capillary is larger than the capillary diameter. The cylindrical capillary geometry used in these models makes it difficult to understand the role of fiber shape and size on wicking in yarns. Even with circular fiber shapes, the cross-sectional shape of the pores between fibers is far from circular. The pore shapes become more complex when the yarn is made of nonround fibers. The goal of this work is to develop a mathematical model of surface-tension-drivkn flow in capillaries of arbitrary cross section. In the next section, we discuss the simplified model developed by Reed and Wilson [16] and derive an extension of this model to flow channels of arbitrary crgss section. In the third section, w e discuss a computational fluid dynamics (CFD) model of this flow, which is solved using Polyflow, a commercial CFD tool. In the last section, we compare results from the simplified model and the CFD calculations for both circular capillaries and channels of arbitrary cross section.

Analytical Flow Modeling Reed and Wilson [I61 studied capillary rise when a circular tube of radius r is dipped in a reservoir of liquid of density p, viscosity 7, and surface tension T (see Figure 1). If the liquid wets the capillary, i.e., the contact angle 8 is less than go”, the liquid will rise into the capillary. This problem and other related flows have received considerable attention in the past. Reed and Wilson wrote a “macroscopic” force balance for this flow that balanced the upward force of surface tension against the downward forces of viscous drag, inertia, and gravity. Under conditions where inertia is not important, they solved this force balance to derive the following expression for the time t required by the liquid to reach a height x above the reservoir:

In Equation 1, .yo is the length of the capillary tube dipped into the reservoir (see Figure 1). At infinite time, the liquid reaches a final height x, when surface

tension forces are balanced by gravity, and this is given by x,=-

2T cos 6 rPg

,

where g is the magnitude of gravitational acceleration. Reed and Wilson compared the predictions of this equation to experimental data for capillary rise in circular tubes and found good agreement.

FIGURE 1. Capillary rise schematic.

W e are interested in generalizing Reed and Wilson’s approach to predict liquid rise in fiber bundles of specified void fraction, fiber shape and size, etc. A first step towards that goal is to extend their approach to an idealized geometric system more representative of fiber bundles than flow in a round capillary. For example, we could consider axial flow in the void spaces formed between a regular array of very long cylinders, a cross section of which is shown in Figure 2. For the purpose of capillary rise modeling, we characterize the geometry by the total perimeter C available for the liquid to wet and the total cross-sectional area A available for liquid flow. For the example shown above, C is simply the total circumference of all the cylinders, and A is the crosssectional area of the spaces between the cylinders. C and A are determined from the radius of the cylinder, and the center-to-center distance of the cylinders, which is a measure of the void fraction or degree of packing of the assembly. Following Reed and Wilson’s approach [16], we analyze the surface-tension-driven flow of liquid in the void spaces between a bundle of cylinders after it has been dipped a distance xo into a reservoir. At time t ,

SEMEhlBER

815

2001

Void Spaces for Flow

am

FIGURE 2. Cross section of array of cylinders.

when the liquid has reached a height x above the reservoir level, the force balance on the liquid column is given by

The viscous drag force F , exerted by the wall on the liquid is given by the product of the wall shear stress and the wall area. Applying this result to the problem of surface-tensiondriven liquid rise between the array of cylinders, the average velocity u is given by dxlclt, the wall area is C( x + x0), and the viscous drag exerted by the cylinders on the fluid is ( q h C 2 / 8 A ) ( x+ x,)dx/dt, which is the term used in the force balance equation (3). The shape factor A depends only on the shape of the conduit and not on its size. Shape factors are available in the literature for many standard shapes, and some of these are catalogued by Miller [14]. For arbitrarily shaped conduits not previously studied in the literature, the shape factor can be computed by numerically solving the equations for laminar, pressure-driven, axial flow in the conduit. This approach is discussed in in the last section for the case of axial flow in between a regular array of cylinders. At long times, the upward force due to surface tension is balanced by the gravitational force, and the meniscus reaches is steady-state position given by x, =

In Equation 3, the first term is the surface tension force, the second is the gravitational force, the third is the viscous drag force, the fourth is the inertial force, and A is a shape factor discussed below. The shape factor (A) in the viscous drag term of Equation 3 arises from a generalization of Reed and Wilson’s derivation to noncircular conduits. Laminar flow in nonround conduits is discussed by Denn [3] (section 8.7). Combining the expression for the Fanning friction factor c f ) in such flows with the definition of the friction factor, we obtain the following expression for the pressure drop per unit length: (4) where u is the axial velocity averaged over the cross section, and D , is the hydraulic diameter defined as 4A D H = c

.

From the overall force balance for axial flow relating the average wall shear stress T,, to the pressure drop, we obtain the following expression for the wall shear stress:

CT COS

e

(7)

%4

In the case of wicking in the void spaces in textile yarns, the size of the channels is so small that inertial forces can be neglected, and the force balance can be converted to a nonlinear first-order differential equation,

by means of the following substitution: (9)

Equation 8 is subject to the initial condition that x = 0 for t = 0, and the solution is

The velocity of the liquid front v = dx/dt is given by v = - dx =-

dt

8pg

qA

C

( , X ? )

,

( 1 1)

and the volumetric flow rate is given by Q = vA. The theoretical analysis reveals that the maximum height x, increases as the wetted perimeter of the channel C increases. In contrast, x, decreases as the contact angle 0 or the cross-sectional area A increases. In the context of filament bundles, this implies that selecting a

TEXTILE RESEARCH JOURNAL

816

filament shape that maximizes the ratio CIA of the channels formed between filaments leads to maximum penetration I,. The capacity of the bundle to hold liquid is defined as x,A, and it is apparent from Equation 7 that this quantity will increase as C increases. As for the rate of liquid rise into the interfilament .channel, Equation 11 shows that the initial velocity at time t = 0 is given by v(t =

8T cos 8 A 0 ) = ~qhr, c .

Thus, the initial velocity of penetration increases as the ratio AIC increases. An interesting prediction of the model is that the maximum height and initial velocity show opposite trends with respect to the ratio C/A. This implies that the maximum height can only be increased at the expense of initial velocity, or vice versa. Reed and Wilson's simple model and the extension derived in this section take a macroscopic view of the fluid dynamics and ignore details of the flo&-near the advancing meniscus and the entrance to the capillary. Nevertheless, Reed and Wilson demonstrated that their model agrees well with experimental data. In the'next section, we will present a more detailed computational fluid dynamics model of this flow and use it to validate the analytical model discussed in this section.

Computational Flow Modeling A detailed computational fluid dynamics (CFD) model of isothermal, surface-tension-driven flow consists of mass and momentum conservation equations and boundary conditions, which are solved by a suitable numerical method. Such a model includes details of the flow field near the advancing meniscus and the entrance to the geometry, which are neglected in the simplified analytical model discussed in the previous section. For isothermal, creeping flow of a Newtonian liquid, the mass and momentum conservation laws are

V.v=O

(13)

and -Pp

define a unit cell by taking advantage of symmetry lines. As shown in Figure 3, the cross section of the unit cell ABCD is contained within three symmetry lines and an arc of the cylinder. The coordinate system is centered on the cylinder, and the flow direction (s)is normal to plane of the cross section.

+ qV2v + pg = 0

,

(14)

where V is the gradient operator, v is the velocity vector, p is the pressure, and g is the gravitational acceleration vector. The inertial terms have been neglected in the momentum equation, because their contribution is unimportant in typical filament bundle systems. These equations are solved in a flow geometry that represents the system of interest. For example, to solve a flow problem for the system shown in Figure 2, we can

FIGURE 3. Flow geometry for computation.

The CFD model is completed by specifying boundary and initial conditions. At the three symmetry planes, we set the tangential stress and normal velocity to zero. At the cylinder wall, we apply the no-slip condition and set all three velocity components to zero. At the bottom of the fiber assembly ( x = -xo), we set the normal force equal to the hydrostatic head of the reservoir above that height, i e . , to pgx,. The surface tension effects are incorporated at the top boundary, which is a free surface that rises vertically into the array of cylinders. At this surface, we set the tangential stress to zero and apply the Young-Laplace equation, which relates the normal stress difference across the interface to the surface tension coefficient and the mean curvature of the interface 'X:

nn : (PI - ~ { ( V V+) ( V V ) ~+) 2XT = 0 ,

(15)

where n is the unit normal vector to the free surface, I is the identity tensor, and (Vv) is the velocity gradient tensor. Equation 15 is integrated by parts using Green's theorem on Riemann surfaces [17] in two dimensions (2D) and in three dimensions (3D), so that only first derivatives of positions are required. This integration introduces a boundary term that corresponds to the static contact angle (8) between the fluid and the wall. Thus, the contact angle appears as a contribution to the momentum equation. At time t = 0, the fluid velocity and pressure are set to zero, and the free surface is assumed to be flat and located at I = 0.

SEMEhlBER

817

2001

The governing equations are solved with a commercial (versions 3.5.4 and 3.6), which uses the finite element method. Flow inside circular capillaries is simulated using both 2Dlaxisymmetric versions of the code as well as the full 3D version. Flow between arrays of cylinders is simulated using the 3D code. The spatial derivatives of the governing equations are discretized by the finite element method, and the partial differential equations are converted into a set of differential algebraic equations. Quadratic interpolation for velocity and linear interpolation for pressure are used in the 2D/axisymmetric simulations. In the 3D case, the mini-element [4] (linear bubble velocity and constant pressure) is used. The time integration involves a predictor corrector scheme with an explicit Euler predictor step and an implicit Euler corrector step. Step size is controlled automatically to maintain discretization errors below a specified error tolerance. Details of the time stepping scheme are available in the Polyflow user's manual. The resulting set of nonlinear algebraic equations is solved by a Newton-Raphson iteration. Polyflow is the only conimercial CFD code that enables this problem to be solved with this numerical technique, in which all the equations are coupled together and ,solved simultaneously by an implicit time stepping scheme. The advantage of this -method is that the time stepping scheme is not subject to stability limits on the step size, and much larger step sizes are possible than with an explicit time integration. This translates into lower execution times for this flow. Polyflow employs a boundaryfitted deforming mesh to solve this problem. The spines method [ 101is used as the remeshing technique to move the internal nodes of the mesh in response to upward motion of the free surface. The numerical solution yields the free surface shape, velocities, and pressure as a function of time, and the integration is carried out until steady state is reached. Polyflow allows a direct specification of a contact angle in the 2Dlaxisymmetric simulation. In the 3D simulations, the contact angle is specified by adding a vertical force contribution related to the vertical angle between the free surface and the wall. The magnitude of the force contribution to the momentum equation is equal to the product of the surface tension coefficient and the cosine of the vertical angle. This force points in the direction of the unit outward tangent vector to the free surface at the dynamic contact line formed by the intersection of the free surface and the cylinder wall. Also, at contact line, the two components of the velocity vector in the plane of the cross section are set to zero. The no-slip boundary condition on the cylinder wall poses a special problem for the finite-element simulation at the dynamic contact line. Details of the fluid dynamics near this three-phase (airlliquidlsolid) interface are not

CFD tool-Polyflow

+

well understood. The problem arises because when we impose the no slip condition, it is not possible for this contact line to move upwards. The common way to resolve this difficulty in finite-element simulations is to allow partial slip along the solid wall [lo]. The strict condition of zero tangential velocity on all or part of the solid wall is replaced by the Navier slip law, which states that the flux of momentum tangent to the wall is proportional to the local slip velocity according to Pt

- v = t n : -q{(Vv)

+ (Vv)?

,

(16)

where t, n are the local unit tangent and normal vectors to the solid wall, and /3 is known as the slip coefficient. Perfect slip is obtained when p = 0, no slip when p + m. In this work, the slip law is applied along the entire solid wall. The appropriate slip coefficient value must usually be determined by some independent means, for example, experimental observations. In this work, we select the value of p such that the simulation results match Reed and Wilson's analytical solution for circular capillaries. All other simulations use the same value of dimensionless slip coefficient (see the next section for details).

Results and Discussion AXISYhlXETRIC SlXlULATlONS '

W e first ran a set of Polyflow simulations to determine the appropriate value of the slip coefficient @ As discussed in the previous section, we adjusted the value of /3 until the 2Dlaxisymmetric simulations matched Reed and Wilson's analytical results for circular tubes. We then used the dimensionless value of p (see below) determined in this manner for all remaining Polyflow simulations. We used the following parameter values for this simulation: circular tube radius r = 2 X lo-' m, surface tension T = 7.2 X N/m, contact angle 8 = 60", fluid density p = 1000 kg/m3. fluid viscosity q = 1 X Pa s, magnitude of gravitational acceleration g = 9.82 m/s2, and initial dip length into the reservoir of xo = 1 X m. W e used a low value of surface tension to limit the final steady-state height of the meniscus, but calculations can be made for any value of surface tension. In Figure 4, we compare the simulation results to the analytical result (Equation I ) for slip coefficient values of 10, 100, 1000, and 10000 kg/s/m2. We report the liquid height above the base of the circular tube. Thus, the initial liquid height is equal to the dip length .yo, and the height above the reservoir is obtained by subtracting the initial dip length from the y-axis value. We use this shifted axis for all the figures in this section.

TEXTILE RESEARCH JOURNAL

818

Analytical

0.0

0.1

0.2

0.3

o

0.4

0.5

FIGURE 4. Effect of slip coefficient on Polyflow results.

The Polyflow simulation matches the analytical solution best for p = 100 kg/s/m2.-Translating this result to determine the appropriate slip coefficient value for other simulations requires some care. When the Navier slip law, Equation 16, is made dimensionless with the fluid viscosity 9 and an appropriate length scale L, we find that a dimensionless slip coefficient is defined as pL/q. For flow in a circular tube, the natural choice for the length scale is the tube radius I-, and the dimensionless slip coefficient value for the simulation result shown in Figure 4 is 20. The appropriate way to determine the dimensionless slip coefficient value (p) for other simulations in circular capillaries is to maintain the value of the dimensionless slip coefficient (PI-Iq)fixed at 20. We verified this by running numerical simulations for a different fluid viscosity and tube radius and adjusting the dimensional slip coefficient value appropriately. The numerical results were in excellent agreement with the analytical curve for liquid rise. We will discuss the appropriate definition of the characteristic length scale I, for geometries other than circular capillaries below. We also conducted mesh and time step refinement studies to verify numerical convergence of the circular capillary simulations, and we tested the simulations against the analytical model by varying material parameters such as the surface tension and contact angle. We found that the axisymmetric Polyflow model correctly predicts liquid rise in a circular capillary.

THREE-DIMENSIONAL SIX~ULATIONS The 3D numerical model allows US to predict liquid rise in more complex geometries such as the void

8

spaces formed by a regular array of cylinders. We discussed the modifications that must be made to the Polyflow model in the previoub section. Before applying this model to complex geometries of interest, we validated the model by using it to solve liquid rise in a circular capillary. We applied the 3D model to a 45" slice of the circular tube and found close agreement between the results and the analytical curve for circular tubes. We now discuss application of the 3D Polyflow model to liquid rise between regular arrays of cylinders. We compare the Polyflow results to the predictions of the simple analytical model derived earlier in the paper. As discussed there, before using the analytical model, we must estimate the value of the shape factor h for the particular geometry of interest. This shape factor depends only on the shape of the flow geometry, not on its size, and it can be calculated by solving the equations for fully developed, laminar, pressure-driven axial flow in the geometry of interest. We solved these equations numerically using Polyflow to obtain the pressure drop per unit length necessary to drive a specified flow rate through the cross section. We solved the problem for the case of a regular square array of cylinders over a range of void fractions, and used the results along with Equation 4 to calculate the values of the shape factor A. The shape factor as a function of void fraction for this geometry is shown in Figure 5, and it varies by almost an order of magnitude over the range of void fractions examined.

4

60 -

-

40 -

-

0 20

40

*

60

80

Void Area ( X ) FIGURE5. Shape factor for axial flow between a square array of cylinders.

100

SEPTEMBER 2001

819

We now discuss application of the 3D Polyfloiv model to solve the problem of liquid rise in the void space between a square array of cylinders. Due to symmetry considerations (see the section on computational flow modeling), we can restrict our attention to the unit cell ABCD in the cross section of Figure 2. We modeled the case of cylinder radius I- = 2 X l 0-4 m and a 50% void fraction, which sets the distance of the point C from the center of the cylinder (L,) according to the following equation:

rr

L = -2 \ ( 1 - F " )

'

0.006

0.005

-E

0.004

x

0.003

-

0.002

(17)

where F,, is the void fraction. The other parameter values used in the simulation were T = 7.2 X N/m, 0 = 50°, p = 1000 kg/m3, q = 1 X Pa s, g = 9.82 m/s2, and xo = 1 X m. As discussed before, we must specify the appropriate value of the slip coefficient so that the dimensionless slip coefficient is held Constant. This requires appropriate definition of the characteristic length scale for this geometry. iVe used the ratio of the area of the unit cell to the length of the cylinder arc as the length scale, because such a measure ' accounts for changes in the shear stresses due to changes in both the cylinder radius and the void area or cylinder spacing. Based on this definition, the length scale for this case equals 1 X lo-' m, and the appropriate slip coefficient value is /3 = 100 kg/s/m2. The Polyflow result is compared to the predictions of the analytical solution in Figure 6, and once again, there is excellent agreement. A snapshot of the finite-element mesh and the free surface shape is shown from the simulation above for t = 6.1 X lo-' s in Figure 7. The effect of void fraction on liquid rise for constant cylinder size is demonstrated in Figure 8. When the void fraction of the square array of r = 2 X lo-' m cylinders is increased from 50% to 70%, the maximum height attained decreases. All other parameters are the same as in Figure 6, and the characteristic length scale value for this geometry is used to specify the appropriate value of /3 = 42.8 kg/s/ni2 so as to maintain the dimensionless slip coefficient constant at 20. Once again, the 3D Polyflow simulation agrees well with the analytical result (shown by the solid line). The result in Figure 6 for 50% void fraction is also shown for comparison, and the height above the reservoir in the 50% void case is 2.34 times the height above the reservoir for the 70% void fraction case. This ratio is exactly equal to the ratio of the cross-sectional area in the 70% case to the 50% case, as one would expect. Examination of Equation 7 for the final steady state reveals that the height above the reservoir is inversely proportional to the cross-sectional area

Polyflouv AnalyLical

0.001

0.0

0

-

0.5

t

1.o

(4

FIGUKE 6. Comparison of Polyflow and analytical results for a 50% void square array of cylinders.

Symmetry Plane

A

/

\

FIGURE 7. Finite element mesh and free surface shape at t = 6.1 X lo-' 5.

of the void space, provided all other factors are held constant as in this comparison. The effect of changing the void fraction from 50% to 70% on liquid rise at early times is shown in Figure 9. In contrast to the final height, the liquid rises higher in the 70% case at early times and the profiles cross over at some intermediate time. Thus, it is possible to arrive at misleading conclusions by examining the liquid rise profiles at intermediate times. The comparison of the liquid rise profiles for the two cases above shows that increasing the void area at constant cylinder radius decreases the final height attained, but in'creases the initial rate of penetration into the array of cylinders. Finally, the effect of perimeter on liquid rise is shown in Figure 10. We considered the case of flow between a regular array of elliptical cylinders with 70% void frac-

TEXTILE RESEARCH JOURNAL above the reservoir for the flow in the array of elliptical cylinders. This is confirmed by the liquid rise profiles shown in Figure 10, in which the profile for the ellipse is compared to that of the circular cylinder. The Polyflow simulations for the elliptical cylinders also agree well with the analytical result shown by the solid line.

tion. We arbitrarily set the ratio of major to minor axes of the ellipse to 2. We set the length of the semi-minor axis to b = 1.4 14 X 10-‘ m, so that the ellipse had the same cross-sectional area as the cylinders discussed above with radii r = 2 X lo-‘ m. Based on the characteristic length for this geometry, the slip coefficient p was set to 46.6 kg/s/m2. All other parameters were the same as those in Figure 8. The elliptical cylinder defined above has an 8.8% higher perimeter than the circular cylinders of the previous calculations. Based on Equation 7, this should translate to an 8.8% higher final liquid height

We have presented analytical and computational models for surface-tension-driven flow in complex geometries. Reed and Wilson’s model [ 161 for surface-tensiondriven liquid rise in a circular capillary forms the starting point for our work. We have discussed computational fluid dynamics (CFD) simulations of this flow using Polyflow, and validated the computational model against Reed and Wilson’s results. We have then extended their analysis to more complex geometries, e.g., flow in the void spaces between a regular array of cylinders. This geometry is a model system for flow in filament bundles. We have also applied the Polyflow simulations to these more complex geometries and confirmed the simplified theoretical analysis we have derived. This extension to Reed and Wilson’s equation for more complex geometries forms the basis of our understanding of flow in filament blfndles. That equation provides a simple way to understand the effect of changing either geometric or material parameters on the wicking process. For example, the model predicts that at constant void area, increasing the perimeter of the filament increases the maximum height attained by the liquid. Con-

82 1

SEPTEMBER 2001 versely, increasing the void area at constant perimeter decreases the final height, but increases the initial rate of liquid penetration. The model also illustrates the effect of changing material parameters such as the contact angle on moisture transport.

Literature Cited 1. Adler, hl. hf., and Walsh, W. K., Mechanisms of Transient hloisture Transport Between Fabrics, Textile Res. J. 5 4 5 ) .

334-342 (1984). 2. Burgeni, A. A., and Kapur, C., Capillary Sorption Equilibria in Fiber hlasses, Textile Res. J. 37(5), 356-366 ( 1967). 3. Denn, hl. hl., “Process Fluid hlechanics,” Prentice Hall, Englewood Cliffs. NJ, 1980. 4. Fortin, M., Old and New Finite Elements for Incompressible Flows, hit. J. M O I LMeth. Fluids 1, 347-363 (1981). 5. Hollies, N. R., Kaessinger, hl. M., and Bogaty, H.. Water Transport hlechanisms in Textile hlaterials, Part I: The Role of Yarn Roughness in Capillary-type Pevetration, Textile Res. J. 26( I I), 829-835 (1956). 6. Hollies, N. R., Kaessinger, M. )I., Watson, B. S., and Bogaty. H., Water Transport Mechanisms in Textile hlaterials, Part 11: Capillary Type Penetration in Yarns and Fabrics, Te-rtile Res. J. 27( I), 8-1 3 (1957 J. 7. Hsieh, Y. L., Liquid Transport in Fabric Structures, Te.rtile Res. J. 65(5), 299-307 (1995). 8. Kamath, Y. K., Hornby, S. B., Weigmann. H.-D., and Wilde, M.F., Wicking of Spin Finishes and Related Liquids into Continuous Filament Yarns, Te.xtile Res. J. 61(1), 33-40 ( 1 993). 9. Kissa, E., Wetting and Wicking, Textile Res. J. 66(10), 660-668 (1996).

10. Kistler, S. F., and Scriven. L. E., Coating Flows, in, “Computational Analysis of Polymer Processing,” J. R. A. Pearson and S. M. Richardson, Eds., Apblied Science Publish-

ers, NY, 1983, pp. 243-299. 1 I. Laughlin, R. D., and Davies, J. E., Some Aspects of Capillary Adsorption in Fibrous Textile Wicking, Texrile Res. J. 31(10), 903-910 (1961). 12. Levine, S., Lowndes, J., Watson, E. J., and Neale, G., A Theory of Capillary Rise of a Liquid in a Vertical Cylindrical Tube and in a Parallel Plate Channel, J. Colloid Irite@ce Sci. 73(1). 136-151 (1980). 13. Lucas, R., Ueber das zeitgesetz des kapillaren aufstiegs von flussigkeiten, Kolloid Z. 23(1), 15-22 (1918). 14. Miller, C.. Predicting Non-Newtonian Flow Behavior in Ducts of Unusual Cross Section, hid. D i g . CImr. Friiidmz. 11(4), 524-528 (1972). 5. Minor, F. W., Schwartz, A. M., Wilkow, E. A., and Buckles, L. C., The Migration of Liquids in Textile Assemblies, Part 11: The Wicking of Liquids in Yarns. Te.rtilc Res. J. 31(12), 931-939 (1961). 6. Reed, C. M.,and Wilson, N., The Fundamentals of Absorbency of Fibres, Textile Structures and Polymers, I: The Rate of Rise of a Liquid in Glass Capillaries, J. Plry. D Appl. PIJJS.26(9). 1378-1381 (1993). 17. Ruschak, K. J., A Three-dimensional Linear Stability Analysis for Two-dimensional Free Boundary Flows by the Finite Element Method, Corrrp. Fluids 11(4), 391-401 (1 983). 18. Washburn, E. W., The Dynamics of Capillary Flow, PIJJS. Rev. 17(3), 273-283 (1921). 19. Woodcock, A. H., hloisture Transfer in Textile Systems, Part 11, Textile Res. J. 32(9), 719-722 (1962). Afaiiirscript received October 5. 2000: nccepprrd Decrritber 17. 2000.