capillary_filling_pf

Solved with COMSOL Multiphysics 4.2 © 2 0 1 1 C O F i l l i n g o f a C ap i l l a r y Chann e l— Pha s e F i e l d Introduction Surface tension and wall adhesive forces are often used to transport fluid through microchannels in MEMS devices or to measure the transport and position of small amounts of fluid using micropipettes. Multiphase flow through a porous medium and droplets on solid walls are other examples where wall adhesion and surface tension strongly influence the dynamics of the flow. This example studies a narrow vertical cylinder placed on top of a reservoir filled with water. Because of wall adhesion and surface tension at the air/water interface, water rises through the channel. The model calculates the velocity field, the pressure field, and the shape and position of the water surface. This example demonstrates how to model the filling of a capillary channel using two physics interfaces available in the Microfluidics Module. You can use either the Laminar Two-Phase Flow, Level Set or the Laminar Two-Phase Flow, Phase Field interface. The Level Set interface uses a reinitialized level set method to represent the fluid interface between the air and the water. The Phase Field interface, on the other hand, uses a Cahn-Hilliard equation, including a chemical potential to represent a diffuse interface separating the two phases. Both interfaces use the Navier-Stokes equations to describe the momentum transport, including surface tension, and conservation of mass. This model uses the phase field method ; for a version using the level set method, see Filling of a Capillary Channel—Level Set. Model Definition The model consists of a capillary channel of radius 0.15 mm attached to a water reservoir. Water can flow freely into the reservoir. Because both the channel and the reservoir are cylindrical, you can use the axisymmetric geometry illustrated in Figure 1. Initially, the thin cylinder is filled with air. Wall adhesion causes water to creep up along the cylinder boundaries. The deformation of the water surface induces surface tension at the air/water interface, which in turn creates a pressure jump across the interface. The pressure variations cause water and air to move upward. The fluids continue to rise M S O L 1 | F I L L I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D Solved with COMSOL Multiphysics 4.2 2 | F I L L I until the capillary forces are balanced by the gravity force building up as the water rises in the channel. In the present example, the capillary forces dominate over gravity throughout the simulation. Consequently, the interface moves upwards during the entire simulation. W al l Sy m m et ry li ne , r = 0 Outflow, p=0 Inflow, hydrostatic pressure Wall Wall Figure 1: Axisymmetric geometry with boundary conditions. R E P R E S E N T A T I O N A N D C O N V E C T I O N O F T H E F L U I D I N T E R F A C E Level Set Method The Laminar Two-Phase Flow, Level Set interface automatically sets up the equations for the convection of the interface. The fluid interface is represented by the 0.5 contour of the level set function  . In air  0= and in water  1= . The level set function can thus be thought of as the volume fraction of water. The transport of the fluid interface separating the two phases is given by The  parameter determines the thickness of the interface. When stabilization is used for the level set equation, you can typically use an interface thickness of hc2, where hc is the characteristic mesh size in the region passed by the interface. The  parameter determines the amount of reinitialization. A suitable value for  is the maximum  t------ u +   – 1 –   ----------   = N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D © 2 0 1 1 C O M S O L Solved with COMSOL Multiphysics 4.2 © 2 0 1 1 C O velocity magnitude occurring in the model. The density and viscosity are calculated in the following way: Due to these definitions, the density and viscosity vary smoothly across the fluid interface. The delta function is approximated by and the interface normal is calculated from Phase Field Method In the Laminar Two-Phase Flow, Phase Field interface the two-phase flow dynamics is governed by a Cahn-Hilliard equation. The equation tracks a diffuse interface separating the immiscible phases. The diffuse interface is defined as the region where the dimensionless phase field variable  goes from 1 to 1. When solved in COMSOL Multiphysics, the Cahn-Hilliard equation is split up into two equations where u is the fluid velocity (m/s),  is the mobility (m3·s/kg),  is the mixing energy density (N) and  (m) is the interface thickness parameter. The  variable is referred to as the phase field help variable. The following equation relates the mixing energy density and the interface thickness to the surface tension coefficient: You can typically set the interface thickness parameter to hc2, where hc is the characteristic mesh size in the region passed by the interface. The mobility parameter  determines the time scale of the Cahn-Hilliard diffusion and must be chosen judiciously. It must be large enough to retain a constant interfacial thickness but small enough so that the convective terms are not overly damped. The default value, 2,  air water air– +=  air water air– +=  6  1 –  = n ----------= t  u +   2 ------=   2 – 2 1– +=  2 2 3 --------------= M S O L 3 | F I L L I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D Solved with COMSOL Multiphysics 4.2 4 | F I L L is usually a good initial guess. This model uses a higher mobility to obtain the correct pressure variation over the interface. In the Phase Field interface, the volume fractions of the individual fluids are In the present model water is defined as Fluid 1 and air as Fluid 2. The density (kg/ m3) and the viscosity (Pa·s) of the mixture are defined to vary smoothly over the interface by letting where the single phase water properties are denoted w and the air properties air. M A S S A N D M O M E N T U M TR A N S P O R T The Navier-Stokes equations describe the transport of mass and momentum for fluids of constant density. In order to account for capillary effects, it is crucial to include surface tension in the model. The Navier-Stokes equations are then Here,  denotes the density (kg/m3),  equals the dynamic viscosity (Ns/m2), u represents the velocity (m/s), p denotes the pressure (Pa), and g is the gravity vector (m/s2). Fst is the surface tension force acting at the air/water interface. Surface Tension In the Level Set interface the surface tension force is computed as Here, I is the identity matrix, n is the interface normal,  equals the surface tension coefficient (N/m), and  equals a Dirac delta function that is nonzero only at the fluid interface. When you use the finite element method to solve the Navier-Stokes equations, you multiply the equations by test functions and then integrate over the computational domain. If you use integration by parts, you can move derivatives of T Vf1 1 – 2 ------------,= Vf2 1 + 2 ------------=  w air w– Vf2+=  w air w– Vf2+=  t u  u  u+  pI–  u u T+ +  Fst g+ +=  u 0= Fst  T= T  I nnT – = I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D © 2 0 1 1 C O M S O L Solved with COMSOL Multiphysics 4.2 © 2 0 1 1 C O to the test functions. This is used in the Laminar Two-Phase Flow interface and results in an integral over the computational domain plus a boundary integral of the form (1) where  is the contact angle (see Figure 2). If you apply a no-slip boundary condition, the boundary term vanishes because test(u)0 on that boundary, and you cannot specify the contact angle. Instead, the interface remains fixed on the wall. However, if you allow a small amount of slip, it is possible to specify the contact angle. The Wetted wall boundary condition adds the term given by Equation 1 and consequently makes it possible to set the contact angle. In the Phase Field interface, the diffuse interface representation makes it possible to compute the surface tension by where  is the phase field parameter, and G is the chemical potential (J/m3) As seen above, the phase field surface tension is computed as a distributed force over the interface using only  and the gradient of the phase field variable. This computation avoids using the surface normal and the surface curvature, which are troublesome to represent numerically. I N I T I A L C O N D I T I O N S Initially, the reservoir is filled with water and the capillary channel is filled with air. The initial velocity is zero. B O U N D A R Y C O N D I T I O N S Inlet The hydrostatic pressure, pgz, gives the pressure at the inflow boundary. Only water enters through the inlet, so the level set function (that is, the volume fraction of water) is 1 here. test u   nwall n cos –   Sd   Fst G= G  2–   2 1–  2 -----------------------+  2 -----= = M S O L 5 | F I L L I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D Solved with COMSOL Multiphysics 4.2 6 | F I L L Outlet At the outlet, the pressure is equal to zero, that is, equal to the pressure at the top of the inflow boundary. Because it is an outflow boundary, you do not have to set any condition on the level set function. Walls The Wetted wall boundary condition is suitable for solid walls in contact with a fluid interface. It sets the velocity component normal to the wall to zero; that is, and adds a frictional boundary force Here,  is the slip length. The boundary condition also allows you to specify the contact angle , that is, the angle between the wall and the fluid interface (see Figure 2). In this example, the contact angle is 67.5 and the slip length equals the mesh element size, h. Air Water Wall  Figure 2: Definition of the contact angle . Results and Discussion The initial development of the fluid interface is shown in Figure 3. During this stage the surface changes drastically in order for it to obtain the prescribed contact angle with the wall. When this is achieved, the surface tension imposed by the surface curvature begins to pull water up through the vertical cylinder. Due to the instantaneous start, the surface oscillates slightly during the rise. u nwall 0= Ffr  ---– u= I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D © 2 0 1 1 C O M S O L Solved with COMSOL Multiphysics 4.2 © 2 0 1 1 C O Figure 3: Snapshots of the position of the interface during the first 0.15 ms. Level Set (left) and Phase Field (right) model results. Figure 4 shows the interface and the velocity field at three different times following the initial stage. After about 0.6 ms the shape of the water surface remains approximately constant and forms a rising concave meniscus. Comparing the velocity field in the Level Set and the Phase Field models, the Level Set results display a small velocity near M S O L 7 | F I L L I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D Solved with COMSOL Multiphysics 4.2 8 | F I L L the wall/interface contact point, something that is not present in the Phase Field results. This is due to a difference in the Wetted wall boundary condition. The Level Set interface requires a wall slip length for the interface to move along the wall. As shown in Figure 4, the imposed slip velocity at the wall is small. In the Phase Field model a slip length is not necessary and the fluid velocity is truly zero on the wall. t = 0.4 ms t = 0.6 mst = 0.2 ms Figure 4: Interface and velocity field at different times. Level Set (top) and Phase Field (bottom) model results. Figure 5 shows surface plots of the pressure at t0.6 ms. At the fluid interface there is a pressure jump of roughly 300 Pa. The jump is caused by the surface tension, and forces the water and air to rise through the vertical cylinder. I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D © 2 0 1 1 C O M S O L Solved with COMSOL Multiphysics 4.2 © 2 0 1 1 C O Figure 5: Pressure at t = 0.6 ms. Level Set (top) and Phase Field (bottom) model results. You can easily calculate the position of the interface/wall contact point by integrating the level set function along the thin cylinder wall. Figure 6 shows the position of the contact point as a function of time. The slight oscillations of the water surface noted above is seen here also in the contact point plot. The contact plots from the Level Set and Phase Field models compare very well, except for two minor points. The surface oscillation is a bit more pronounced in the Level Set model, and the surface end point M S O L 9 | F I L L I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D Solved with COMSOL Multiphysics 4.2 10 | F I L is somewhat higher up in this case. Both these differences are small and are most likely related to the different implementations of the Wetted wall boundary condition. Figure 6: Position of the interface/wall contact point as a function of time. Level Set (top) and Phase Field (bottom) model result. The velocity is approximately constant after t = 0.6 ms. Finally, you can verify the obtained contact angle. It is defined by cos nTnwall. In this case, the normal to the wall is nwaller. The contact angle is thus acos nr, where nr is the radial component of the interface normal. Due to the slight oscillations of the surface, the contact angle varies during the rise. As Figure 7 shows, at t0.6 ms the contact angle is  for both the Level Set and the Phase Field models. Both results L I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D © 2 0 1 1 C O M S O L Solved with COMSOL Multiphysics 4.2 © 2 0 1 1 C O are close to the imposed contact angle of 381.18rad 67.5. The contact angle further approaches the imposed value if the mesh is refined. Figure 7: Plot of acos(nr). At the wall, this gives the contact angle. In the Level Set model (top) the wall angle is  and in the Phase Field model (bottom) it is 67.8. Notes About the COMSOL Implementation The model is straightforward to set up using either the Laminar Two-Phase Flow, Level Set or the Laminar Two-Phase Flow, Phase Field interface. At walls in contact with the fluid interface, you can use the Wetted wall boundary condition. M S O L 11 | F I L L I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D Solved with COMSOL Multiphysics 4.2 12 | F I L The simulation procedure consists of two steps. First the phase field and level set functions are initialized, then the time dependent calculation starts. This is automatically set up by the software. You only need to specify appropriate times for the intialization step and the time-dependent analysis. Model Library path: Microfluidics_Module/Two-Phase_Flow/ capillary_filling_pf Modeling Instructions M O D E L W I Z A R D 1 Go to the Model Wizard window. 2 Click the 2D axisymmetric button. 3 Click Next. 4 In the Add physics tree, select Fluid Flow>Multiphase Flow>Two-Phase Flow, Phase Field>Laminar Two-Phase Flow, Phase Field (tpf). 5 Click Next. 6 In the Studies tree, select Preset Studies>Transient with Initialization. 7 Click Finish. G E O M E T R Y 1 1 In the Model Builder window, click Model 1>Geometry 1. 2 Go to the Settings window for Geometry. 3 Locate the Units section. From the Length unit list, select mm. Rectangle 1 1 Right-click Model 1>Geometry 1 and choose Rectangle. 2 Go to the Settings window for Rectangle. 3 Locate the Size section. In the Width edit field, type 0.3. 4 In the Height edit field, type 0.15. 5 Locate the Position section. In the z edit field, type -0.15. 6 Click the Build Selected button. L I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D © 2 0 1 1 C O M S O L Solved with COMSOL Multiphysics 4.2 © 2 0 1 1 C O Rectangle 2 1 In the Model Builder window, right-click Geometry 1 and choose Rectangle. 2 Go to the Settings window for Rectangle. 3 Locate the Size section. In the Width edit field, type 0.15. 4 In the Height edit field, type 0.5. 5 Click the Build Selected button. 6 Click the Zoom Extents button on the Graphics toolbar. Form Union In the Model Builder window, right-click Form Union and choose Build Selected. M A T E R I A L S 1 In the Model Builder window, right-click Model 1>Materials and choose Open Material Browser. 2 Go to the Material Browser window. 3 Locate the Materials section. In the Materials tree, select Built-In>Air. 4 Right-click and choose Add Material to Model from the menu. 5 In the Model Builder window, right-click Materials and choose Open Material Browser. 6 Go to the Material Browser window. 7 Locate the Materials section. In the Materials tree, select Built-In>Water, liquid. 8 Right-click and choose Add Material to Model from the menu. L A M I N A R TW O - P H A S E F L O W, P H A S E F I E L D Fluid Properties 1 1 In the Model Builder window, expand the Model 1>Laminar Two-Phase Flow, Phase Field node, then click Fluid Properties 1. 2 Go to the Settings window for Fluid Properties. 3 Locate the Fluid 1 Properties section. From the Fluid 1 list, select Air. 4 Locate the Fluid 2 Properties section. From the Fluid 2 list, select Water, liquid. 5 Locate the Phase Field Parameters section. In the  edit field, type 50. 6 In the pf edit field, type 6.5e-6. Initial Values 2 1 In the Model Builder window, right-click Laminar Two-Phase Flow, Phase Field and choose Initial Values. M S O L 13 | F I L L I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D Solved with COMSOL Multiphysics 4.2 14 | F I L 2 Select Domain 1 only. 3 Go to the Settings window for Initial Values. 4 Locate the Initial Values section. Click the Fluid 2 button. Initial Interface 1 1 In the Model Builder window, click Initial Interface 1. 2 Select Boundary 4 only. Inlet 1 1 In the Model Builder window, right-click Laminar Two-Phase Flow, Phase Field and choose Inlet. 2 Select Boundary 8 only. 3 Go to the Settings window for Inlet. 4 Locate the Boundary Condition section. From the Boundary condition list, select Pressure, no viscous stress. 5 In the Vf edit field, type 1. 6 Locate the Pressure, No Viscous Stress section. In the p0 edit field, type -tpf.rho2*z*g_const. This is the hydrostatic pressure at the inlet. Outlet 1 1 In the Model Builder window, right-click Laminar Two-Phase Flow, Phase Field and choose Outlet. 2 Select Boundary 5 only. Wall 2 1 In the Model Builder window, right-click Laminar Two-Phase Flow, Phase Field and choose Wall. 2 Select Boundaries 6 and 7 only. 3 Go to the Settings window for Wall. 4 Locate the Boundary Condition section. From the Boundary condition list, select Wetted wall. 5 In the w edit field, type (3*pi/8)[rad]. Gravity 1 1 In the Model Builder window, right-click Laminar Two-Phase Flow, Phase Field and choose Gravity. 2 Go to the Settings window for Gravity. L I N G O F A C A P I L L A R Y C H A N N E L — P H A S E F I E L D © 2 0 1 1 C O M S O L Solved with COMSOL Multiphysics 4.2 © 2 0 1 1 C O 3 Locate the Domain Selection section. From the Selection list, select All domains. D E F I N I T I O N S Next, define a variable for the contact angle. Variables 1 1 In the Model Builder window, right-click Model 1>Definitions and choose Variables. 2 Go to the Settings window for Variables. 3 Locate the Variables section. In the Variables table, enter the followin