Meniscus behaviors and capillary pressures in capillary channels having various

Yicun Tang ,Jingchun Min ,*,Xuan Zhang ,Guiling Liu

1 Department of Engineering Mechanics,Tsinghua University,Beijing 100084,China

2 Department of Energy and Power Engineering,Tsinghua University,Beijing 100084,China

Keywords:Capillary phenomenon Polygonal channel Meniscus shape Capillary pressure

ABSTRACT A numerical study has been conducted to simulate the liquid/gas interface(meniscus)behaviors and capillary pressures in various capillary channels using the volume of fluid(VOF)method.Calculations are performed for four channels whose cross-sectional shapes are circle,regular hexagon,square and equilateral triangle and for four solid/liquid contact angles of 30°,60°,120°and 150°.No calculation is needed for the contact angle of 90°because the liquid/gas interface in this case can be thought to be a plane surface.In the calculations,the liquid/gas interface in each channel is assumed to have a flat surface at the initial time,it changes towards its due shape thereafter,which is induced by the combined action of the surface tension and contact angle.After experiencing a period of damped oscillation,it stabilizes at a certain geometry.The interface dynamics and capillary pressures are compared among different channels under three categories including the equal inscribed circle radius,equal area,and equal circumscribed circle radius.The capillary pressure in the circular channel obtained from the simulation agrees well with that given by the Young-Laplace equation,supporting the reliability of the numerical model.The channels with equal inscribed circle radius yield the closest capillary pressures,while those with equal circumscribed circle radius give the most scattered capillary pressures,with those with equal area living in between.A correlation is developed to calculate the equivalent radius of a polygonal channel,which can be used to compute the capillary pressure in such a channel by combination with the Young-Laplace equation.

1.Introduction

Capillary phenomena exist in many fields such as hydrology,petroleum industry,and polymer electrolyte fuel cells(PEFCs)area,and the capillary action has a broad range of applications[1-5].Early studies on capillarity focused mainly on the capillary rise,e.g.,Green and Ampt[6]studied the permeability of soil to air and water,and developed a theoretical model for calculating the capillary rise in a thin tube,they assumed the porous soil to consist of a bundle of capillary tubes which were irregular in area,length,direction and shape,and derived the liquid velocity equation based on Poiseuille's law.Lucas[7]and Washburn[8]further developed this theory independently,they derived the relation between the meniscus height and time,which is often called the Lucas-Washburn equation,based on the momentum balance with ignorance of the gravitational and kinetic effects.Efforts were then made to justify the Lucas-Washburn equation by including the missing or neglected terms such as the surface free energy,contact angles,and gravity and kinetic effects,which may have significant in fluences on the capillary behaviors.Masoodi et al.[9]derived the governing equations for capillary rise in a vertical tube using the energy balance with consideration of the kinetic,gravity,and viscous effects.They normalized and numerically solved the governing equations and reported that the obtained results could reflect the liquid rise and oscillation in the tube.Fries and Dreyer[10]developed an analyticalsolution for the capillary rise of liquid in a cylindrical tube.Their solution included the gravity term and predicted a liquid rise which lasted a longer time than that calculated by the classic Lucas-Washburn equation.Lago and Araujo[11]proposed a model describing the capillary pressure evolution in a regular assembly of spheres based on a quasi-static advance of meniscus with a piston-like motion.The results generated by the model were very close to those obtained by fitting the data to the Washburn equation.Siebold etal.[12]experimentally investigated the effects of dynamic contact angle on capillary rise and proposed various approaches to determine the real constant term in the calculation of particle size.Fries[13]systematically studied the capillary transport phenomena as well as their in fluencing factors including the viscous,gravity and evaporation affects.Zhmud et al.[14]reviewed the classical theory of capillarity and investigated the capillary rises of surfactant solutions,in which a simplified relation was developed for capillary rise dynamics in the case of strong depletion of the interfacial region and a qualitative agreement with the experiment was obtained.

The pores and throats in porous media are often assumed to consist of capillary tubes with different radii and cross-sectional shapes,and a pore network model can thus be used to simulate the transport phenomena in porous media[15-18].Prat[15]reported that the pore shape and contact angle significantly in fluenced the rate of isothermal drying of porous materials and a better agreementwith the experiment could be achieved when film effect was taken into account in modeling of the pore network.Segura and Toledo[16]investigated the effects of gravity,pore shape and pore size distribution on saturation and transport characteristics during isothermal drying process.Their results suggested that the pore network dimensionality,pore size skewness and variance,as well as gravity effect all had pronounced in fluences on saturation distribution.The liquid menisci and capillary pressures in porous structures are important to understand the phase equilibrium and phase interface dynamic behaviors in porous media.Wong et al.[17]proposed a method to solve the Young-Laplace equation indirectly by including conjoining/disjoining thin film forces,they experimentally investigated the meniscus shape in a square glass capillary and reported that the calculation agreed well with the experiment.Mason and Morrow[18]theoretically explored the methods for evaluating the meniscus curvatures in capillary channels with a uniform cross-section.The conditions of pore geometry and contact angles which gave rise to wedge menisci are discussed and presented by examples including liquid menisci in tubes of polygonal sections.The capillary-induced flow can also be described using the Leverett approach[19,20],and the relative Leverett function may explicitly include the surface wetting effects[19],which are important to many applications such as vapor condensation,dehumidifying heat exchangers,etc.[21-24].The driving force forliquid movementwas the capillary pressure caused by the difference in the energies of dry and wet surfaces in the tubes[25].

The Young-Laplace equation can be used to calculate the capillary pressure(capillary force)if the meniscus curvature is known.However,such a curvature is usually unknown for liquid meniscus in a noncircular channelbecause the meniscus geometry in this case is nota sphericalcap.Different approaches can be used to simulate the multiphase flow in pore spaces,they include the level-set method,lattice Boltzmann method,smooth particle hydrodynamics,etc.[26-28].These approaches allow to obtain accurate results for any geometry whereas the complicated governing equations and computational cost restrict their applicability.So,it would be helpful if a correlation relating the capillary pressure and irregular pore shape can be provided to calculate the capillary pressure based on the pore shape.

In the present research,the geometries of liquid menisci in capillary channels with various cross-sectional shapes are simulated numerically using the volume of fluid method,and the capillary pressures,which are defined as the pressure difference between the liquid and gas fluids,are determined based on the meniscus geometries.Simulationsare conducted for four channels whose cross-sectionalshapes are circle,regular hexagon,square and equilateral triangle,and for four wetting angles of the liquid to channel wall whose values are 30°,60°,120°and 150°.In the simulations,the liquid/gas interface in each channel is assumed to have a flat surface initially,it changes thereafter and stabilizes eventually at a certain shape.The capillary pressures are compared among different channels under three conditions including equal cross-sectional area and equal internally and externally tangent circle radii.A correlation is provided to calculate the equivalent radius of a polygonal channel,which can be used to compute the capillary pressure in such a channel by combining it with the Young-Laplace equation.

2.Theoretical Model

2.1.Physical model

Consider the liquid/gas interfaces in four capillary channels whose cross-sectional shapes are circle,regular hexagon,square and equilateral triangle,respectively,asillustrated in Fig.1.Each channelhas an inscribed circle radius of r0and a heightof H,ithas a sealed bottomand an open top exposed to the ambientair.Allchannels are partly filled with liquid water thathas a flatsurface in contactwith the air and a heightof h atthe initial time.The surface then change towards a curved one due to the surface tension effect and the inherent contact angle between the water and channel wall(θ).The surface undergoes an oscillation process before it eventually stabilizes at its due shape.The capillary pressure(capillary force),which is defined as the pressure difference between the two sides of the air/water interface,can be obtained from the static pressures of the air and water fluids.The coordinate system adopted in this work is shown in the figure.

2.2.Mathematic model

The simulation considers channels in small scale so the gravitational effect can be ignored.Further,it considers the transient fluid dynamics during the formation of the free liquid/gas interface in each channel,no thermal effect is present.The mass conservation equation can thus be represented by where ρ is the density,and u,v and w denote the velocities in the x,y and z directions.Correspondingly,the momentum equation can be expressed as follows:

Fig.1.Liquid/gas interface in capillary tubes whose cross-sectional shapes are:(a)circle,(b)regular hexagon,(c)square,and(d)equilateral triangle.

x-direction

y-direction

z-direction

where p is the static pressure and F stands for the body force term.

The volume of fluid(VOF)model is often used to track the free surfaces between immiscible fluids or the interfaces between fluid and deformable structures[29,30].The model introduces a parameter α to express the volume fraction of each fluid in element cells,it uses αlandαgto express the volume fractions ofliquid and gas,αl=1 indicates that the cell is full of liquid while αl=0 means that the cell contains no liquid,and 0＜αl＜1 implies that the cell contains both liquid and gas whose interface is a free surface that needs to be further identified.It is usually assumed that the liquid/gas interface appears at a volume fraction of 0.5 as in Refs.[31,32],and its tracking can be realized by solving the continuity equations of fluid phases.The volume fraction equation for liquid and gas phases can be represented by

where the subscripts l and g express the liquid and gas phases,respectively.Obviously,the volume fractions of the liquid and gas satisfy

The physical properties and other parameters all take the volumeweighted average values,e.g.,the density in each cell can be obtained from

By solving the governing equations for the volume fractions of different fluids at each time step,the free surface can be reconstructed based on the calculated volume fractions.

Tracking of an interface between immiscible fluids can be implemented using the geo-reconstruct scheme[33],which is an accurate but computationally expensive scheme.The interface can be obtained by connecting the line segments in adjacent cells as depicted by Fig.2.For cell containing a free surface,its normal vectorcan be represented by

which defines the included angle of the free surface by

Normalization of it gives

The free surface is determined by the normalized inclination angle and volume fraction in each cell,and the whole interface can be achieved by connecting all line segments in the relevant cells.

The continuum surface force(CSF)model[34]is used as the source term in the momentum equation in order to include the surface tension effect in the results generated by the VOF model.For a curved interface(meniscus)between liquid and gas fluids,the pressure difference between the two sides of the meniscus is determined by the surface tension coefficientσand the surface curvature measured by two principal radii in orthogonal directions R1and R2:

Fig.2.Interface construction through the geo-reconstruct scheme.

where pland pgare the pressures of liquid and gas fluids.The surface curvature can be calculated from the volume fraction gradientperpendicular to the surface.Ifisused to expressthe surface normal,i.e.the gradientof the volume fraction function α,then

The curvature κ is defined as the divergence of the unit normal^n,i.e.

The surface tension can be linked to the pressure difference across the interface,which can be converted into the volume force term based on the divergence theorem and is used as the source term in the momentum equation.It has a form of

For a liquid/gas system,κl+ κg=0,∇αl+ ∇αg=0,Eq.(12)thus becomes

where ρ is the volume-averaged density,given by Eq.(5).

The option to specify a wall adhesion angle in conjunction with the surface tension model is available in the VOF model[34].Instead of imposing this boundary condition at the wall itself,the contact angle specified to representthe wall wettability is used to adjust the interface normal in element cells near the wall,and the interface normal at cells next to the wall can be expressed as

As stated early,all channels have open tops exposed to the ambient fluid,which has a pressure of po,it also serves as a reference pressure in this research.The boundary conditions can thus be expressed as:at the channel walls:u=v=w=0,=0;at the gas domain:p=po.

3.Numerical Method

The preceding setofequations is numerically solved with the help of the CFD code Fluent.The liquid/gas interface is tracked by the georeconstruct scheme embedded in the volume of fluid(VOF)model.For the balance of the pressure gradient in the momentum equation,the implicit body forces are taken into account in the VOF model to mprove the convergence of the calculations.The pressure implicit split operator(PISO)algorithm is adopted for faster convergence by using a two-step modification method which allows a larger time step.The least squares cell-based gradient evaluation method is used for the spatial discretization.The pressure staggering option method is employed for the pressure terms,with the second order upwind scheme being used for the momentum terms.A double precision representation of real numbers is adopted to reduce round off errors.

4.Results and Discussion

4.1.Meniscus and capillary pressure in a circular channel and model validation

Simulations are first conducted to calculate the meniscus geometry and capillary pressure in a two-dimensional(2D)circular channel containing air and water fluids for four liquid/solid contact angles of θ=30°,60°,120°and 150°.As shown by Fig.3,the water/air interface is assumed to be flat at the initial time.It begins to change when the calculation starts because the value specified for the contact angle differs from that of the initial angle which is 90°(corresponding to the flat surface),the water/air interface eventually reaches a curved shape,with the contactangle becoming equal to the specified one.Table 1 presents the values taken for some key variables included in the capillary model in the simulations,in which the geometry creation and grid generation are completed in Gambit,and the subsequent calculations are implemented in Fluent,with the time step being set as 10-6s.A grid independence investigation is carried out for three grid systems of 20×80,40×160,80×320,the difference in the calculated capillary pressure between the first two grid systems is 0.48%and that between the last two grid systems is 0.37%,the grid system of 40×160 is eventually adopted for the simulations.

Fig.3.Meniscus shape evolution in a 2D circular channel.

Table 1 Values taken for some key variables in the calculations

Fig.4.Evolution of liquid/gas interface shape in a circular channel for contact angle of 30°.

Fig.4 illustrates the evolution of liquid/gas interface geometry in a circular channel whose radius is r0=0.5 mm for the contact angle of θ =30°.Since the initial contact angle is set as 90°,the liquid climbs up along the channel wall to make the wetting angle change towards the specified value of 30°,as a result,the neighboring liquid surface caves in,forming a w-shape liquid surface profile at an early stage and a U-shape profile afterwards,and the deepest concave appears at about t=1.0 ms.The lowest point turns to rise thereafter,and the centralpartbecomes nearly flatat t=1.8 ms,with the whole liquid surface becoming almost unchanged when the time reaches t=7.0 ms.Fig.5 is similar to Fig.4,itis for the contactangle ofθ =150°.Since the inherent contact angle is specified as θ =150°,which is larger than the initial contact angle of 90°,the liquid climbs down along the channel wall,forming a convex profile.The liquid surface stabilizes when the time attains to t=15.0 ms.It is important to note that what Figs.4 and 5 show are not a convergence process of the numerical calculations but a real physical process,which means that if one can initially keep a liquid surface flat and then let it go,the surface would experience a change in shape as shown by these figures.

Fig.6 presents the evolutions of the positions ofliquid surface center for five contact angles of θ =30°,60°,90°,120°and 150°.For θ =90°,the liquid surface is flat and its center remains stationary throughout,no calculation is needed in this case.For the other contact angles,the liquid surface center always experiences a period of oscillation with attenuating amplitudes and eventually stabilizes at a position differing from that at the beginning.For θ =30°and 60°,the channel wall is hydrophilic and the liquid surface centerstabilizes ata position lowerthan the origin,with a smaller contact angle yielding a larger deviation.For θ=120°and 150°,the channelwallis hydrophobic and the liquid surface center stabilizes at a position higher than the origin,with a larger contact angle giving a larger deviation.Meanwhile,the duration of oscillation tends to increase with increasing contact angle,it is approximately 7 ms for θ =30°and 15 ms for θ =150°,this is why the last picture in Fig.4 is for t=7.0 ms and that in Fig.5 is for t=15.0 ms.

Fig.7 depicts the distributions of static pressure along the channel axis at the time when the water/air interface becomes stable for five contact angles of θ=30°,60°,90°,120°and 150°,and Table 2 compares the capillary pressures obtained from the simulation and those given by the Young-Laplace equation for the same contact angles.The gas pressure serves as the reference pressure and is presented as zero in Fig.7.For θ =90°,the liquid surface is flat,so the liquid and gas phases possess the same pressure;for the other contact angles,a jump in pressure exists at the liquid/gas interface,and the capillary pressure can be obtained by calculating the pressure difference between the liquid and gas phases as defined by Eq.(9).The calculation results are presented in Table 2,which compares the capillary pressures obtained from the simulation and those from the Young-Laplace equation.For liquid meniscus in a circular channel,the two principal radii of the meniscus are equal to each other,and Eq.(9)i.e.the Young-Laplace equation therefore becomes

where r0is the radius of the circular channel.Table 2 indicates that the capillary pressures given by the simulation agree well with those by the Young-Laplace equation,they agree within 1.4%,supporting the reliability of our numerical procedure.

Fig.5.Evolution of liquid/gas interface shape in a circular channel for contact angle of 150°.

4.2.Effects of the cross-sectional shape on the meniscus and capillary pressure

Fig.6.Evolutions of water/air interface center position in a circular channel for various contact angles.

Fig.7.Distributions of static pressure along channel axis for various contact angles.

Table 2 Capillary pressure comparison between the simulation and Young-Laplace equation for circular channel for various contact angles

Simulations are then performed to calculate the meniscus shapes and capillary pressures in channels whose cross-sectional shapes are circle,regular hexagon,square and equilateral triangle.A grid system containing about0.3 million hexahedralmeshes is adopted to discretize the 3D models for all cases,which is chosen based on the results of grid independence investigations.Concretely speaking,the grid number is 299460 for the circular channel,289800 for the hexagonal channel,300000 for the square channel,and 340740 for the triangular channel.Taking the triangular channel with 60°contact angle as an example,a grid independence study is conducted for three grid systems of 163920,340740 and 771360,the difference in the calculated capillary pressure between the first two grid systems is 1.4%and that between the last two grid systems is 0.34%,the grid system of 340740 is eventually used for the simulations.Meanwhile,the values taken for some key variables in the calculations are the same as those presented in Table 1.Since the channel scale has a significant in fluence on the capillary behaviors,comparisons are made among different channels under three categories including:(1)equal inscribed circle radius,(2)equal area,and(3)equal circumscribed circle radius,which are called scale A,scale B,and scale C,respectively.Fig.8 compares the evolutions of the water/air interfaces in channels having circular,hexagonal,square and triangular cross-sectional shapes whose inscribed circle radii are all the same as r0=0.5 mm(belonging to the scale A case)for contact angle of θ =60°.The figure shows that the water/air interface undergoes a damped oscillation and eventually stabilizes at a geometry differing from a perfect spherical crown in polygonal channels,this is why the Young-Laplace equation cannot be directly used to compute the capillary pressure in such channels.The time length of the oscillation sorts in the order of the triangular,square,hexagonal and circular channels,implying that the larger the number of the polygon edges,the longer the duration of the oscillation.

Fig.9 compares the capillary pressures in channels having circular,hexagonal,square and triangularcross-sectionalshapes forvarious contact angles under scales A,B,and C.For all cases,the circular channel serves as a baseline channel whose radius is r=r0=0.5 mm in which r0is the inscribed circle radius of a polygonal channel.The figure shows that the capillary pressures in these channels are closest to one another under scale A but scatter most seriously under scale C,meaning that the inscribed circle radius is the most appropriate dimension in evaluating the capillary pressure in a non-circular polygonal channel.

Fig.8.Evolutions of water/air interfaces in channels whose cross-sectional shapes are:(a)circle,(b)regular hexagon,(c)square,and(d)equilateral triangle,for contact angle of 60°.

Fig.9.Comparisons of capillary pressures in channels having circular,hexagonal,square and triangular cross-sectional shapes for various contact angles under scales A,B,and C.

Although scale A is superior to scales B and C in evaluating the capillary pressure in a polygonal channel,it still cannot yield equal capillary pressure in differently shaped channels,as seen in Fig.9.For this reason,a correlation is developed based on the Fig.9 data to calculate the equivalent radius(r*)of a polygonal channel,i.e.

Fig.10.Comparison of capillary pressures in various channels obtained from the simulations and from Eqs.(15)and(16)for various contact angles.

where r0is the inscribed circle radius of a polygonal channel,N is the number of polygon edges,and the sign “±”takes positive for θ≤ 90°and negative for θ N 90°.Eq.(16)is valid for regular polygonal channels overa channelsurface contactangle range from30 to 150°.The capillary pressure in a polygonalchannelcan be calculated by combining Eq.(16)with Eq.(15).For example,for θ =60°,Eq.(16)yields r*/r0=1 for the circular channel,r*/r0=1.052 for the hexagonal channel,r*/r0=1.083 for the square channel,and r*/r0=1.121 for the triangular channel.For r0=0.5 mm,the equivalent radius is r*=0.5 mm for the circular channel,r*=0.526 mm for the hexagonal channel,r*=0.542 mm for the square channel,and r*=0.561 mm for the triangular channel.The capillary pressures in these channels differ(see Fig.9)because their equivalent radii are different.Fig.10 compares the capillary pressures in various channels obtained from the simulations and from Eqs.(15)and(16)for various contact angles.The figure shows that the capillary pressures given by Eqs.(15)and(16)agree with those generated by the simulation within 1.5%,supporting the effectiveness ofthe combination ofEqs.(15)and(16)in evaluating the capillary pressure in a polygonal channel.

5.Conclusions

The meniscus behaviors and capillary pressures in capillary channels whose cross-sectional shapes are circle,regular hexagon,square and equilateral triangle are numerically investigated with the help of the Fluent software,and the following conclusions can be drawn:

·The capillary behaviors in circular and polygonal channels can be simulated numerically and the capillary pressure can be obtained from the fluid static pressure distribution.

·The capillary pressure in a circular channel obtained from the simulation agrees well with that given by the Young-Laplace equation,supporting the effectiveness of the numerical approach.

·When the liquid/gasinterface in a capillary channelsetsoutfrom a flat surface,it undergoes an oscillation with attenuating amplitude and eventually stabilizes at a geometry which is determined by the channel shape and liquid/solid contact angle.

·The capillary channels with equal inscribed circle radius yield the closestcapillary pressures while those with equalcircumscribed circle radius give the most scattered capillary pressures,with those with equal area in between.

·The capillary pressure in a polygonal channel can be calculated by combining the Young-Laplace equation with an equivalent radius correlation,which can be developed based on the simulation results obtained for various channels and various contact angles.

Nomenclature

F volume forces

H channel height,m

h liquid column height,m

N number of edges of polygonal channel

^n unit normal vector

^nwunit vector normal to the wall

p pressure,Pa

p0atmospheric pressure,Pa

pccapillary pressure,Pa

R1,R2principal radii in orthogonal directions,m

r radius of circular channel,m

r* equivalent radius of polygonal channel,m

r0inscribed circle radius,m

T0system temperature,K

t time,s

^twunit vector tangential to the wall

u,v,w velocities in x,y and z directions,m·s-1

x,y,z coordinate system,m

α volume fraction

θ contact angle,(°)

ĸ curvature,m-1

μ dynamic viscosity,kg·m·s-1

ξ angle of free surface in cell,(°)

ρ density,kg·m-3

Subscripts

a air

g gas phase

l liquid phase