COMSOL 6.4 - Capillary Filling

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 multiphysics coupling features available in the CFD Module. You can use either the Two-Phase Flow, Level Set or the Two-Phase Flow, Phase Field multiphysics coupling features. The Level Set interface uses a reinitialized level set method to represent the fluid interface between the air and the water. The Phase Field in Fluids interface, on the other hand, uses a Cahn–Hilliard equation, including a chemical potential to represent a diffuse interface separating the two phases. The Navier–Stokes equations are used to describe the momentum transport and the conservation of mass.

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 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 upward during the entire simulation.

Figure 1: Axisymmetric geometry with boundary conditions.

Representation and Convection of the Fluid Interface

Level Set Method

The 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 velocity magnitude occurring in the model. The multiphysics coupling feature defines the density and viscosity according to:

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 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 value γ = ε2 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 multiphysics coupling feature defines the density (kg/m3) and the viscosity (Pa·s) of the mixture to vary smoothly over the interface by letting

where the single phase water properties are denoted w and the air properties air.

Mass and Momentum Transport

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, n is the interface normal, σ is the surface tension coefficient (N/m), κ = −∇ ⋅ n is the curvature, and δ equals a Dirac delta function that is nonzero only at the fluid interface.

The following boundary force is added to enforce the contact angle:

(1)

where θ is the contact angle (see Figure 2). If you apply a no slip boundary condition, the velocity vanishes 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 coupling feature 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.

Initial Conditions

Initially, the reservoir is filled with water and the capillary channel is filled with air. The initial velocity is zero.

Boundary Conditions

Inlet

The hydrostatic pressure, p = ρgz, gives the pressure at the inflow boundary. The pressure boundary condition automatically compensates for hydrostatic pressure so the actual value for the pressure is set to zero. Only water enters through the inlet, so the volume fraction of water is 1 here.

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.

Wetted Walls

The Wetted Wall feature is suitable for solid walls in contact with a fluid-fluid interface. In the Level Set case, it is introduced as a multiphysics coupling feature that 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.

Figure 2: Definition of the contact angle θ.

In the Phase Field simulation, one can combine the Wall feature under the Laminar Flow physics interface imposing a no slip condition for the velocity with the Wetted Wall feature under the Phase Field physics interface imposing 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.

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 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 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.

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.

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 are 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 endpoint 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 approximately 69° for the Level Set model and approximately 68° for the Phase Field model. Both results 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 approximately 69° and in the Phase Field model (bottom) it is approximately 68°.

Notes About the COMSOL Implementation

The model is straightforward to set up using either the Level Set or the Phase Field interface. At walls in contact with the fluid interface, you can use the Wetted Wall coupling feature for the Level Set case.

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 COMSOL Multiphysics. You only need to specify appropriate times for the time-dependent analysis.

CFD_Module/Multiphase_Flow/capillary_filling_ls

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > > .

Click .

Click

In the tree, select > .

Click

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.3.

In the text field, type 0.15.

Locate the section. In the text field, type -0.15.

Click

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.15.

In the text field, type 0.5.

Click

Form Union (fin)

In the window, click .

In the window for , click

Multiphysics

Two-Phase Flow, Level Set 1 (tpf1)

In the window, under > click .

In the window for , locate the section.

Click

In the window, click .

Locate the section. Select the checkbox.

Materials

Phase 1 (mpmat1.phase1)

In the window, under > > click .

In the window for , locate the section.

Click

. This button is found when expanding the options next to the list.

Add Material to Phase 1 (mpmat1.phase1)

Go to the window.

In the tree, select > .

Click .

Materials

Phase 2 (mpmat1.phase2)

In the window, under > > click .

In the window for , locate the section.

Click

. This button is found when expanding the options next to the list.

Add Material to Phase 2 (mpmat1.phase2)

Go to the window.

In the tree, select > .

Click .

Multiphysics

Wetted Wall 1 (ww1)

In the window, under > click .

Select Boundaries 6 and 7 only.

In the window for , locate the section.

In the θw text field, type (3*pi/8)[rad].

Level Set (ls)

Level Set Model 1

In the window, under > click .

In the window for , locate the section.

In the εls text field, type 5e-6.

Initial Values, Fluid 2

In the window, click .

Select Domain 1 only.

Add gravity volume force contribution by enabling gravity in the interface.

Laminar Flow (spf)

In the window, under click .

In the window for , locate the section.

Select the checkbox.

Inlet 1

In the toolbar, click

and choose .

Select Boundary 8 only.

In the window for , locate the section.

From the list, choose .

Level Set (ls)

In the window, under click .

In the toolbar, click

and choose .

Select Boundary 8 only.

In the window for , locate the section.

From the list, choose φ.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 5 only.

Laminar Flow (spf)

In the window, under click .

In the toolbar, click

and choose .

Select Boundary 5 only.

In the window for , locate the section.

Clear the checkbox.

Next, define a variable for the contact angle using the interface’s variable for the r component of the interface normal, tpf.intnormr.

Initial Values 2

In the toolbar, click

and choose .

Select Domain 2 only.

In the window for , locate the section.

Clear the checkbox.

To set the correct hydrostatic pressure, the “Compensate for hydrostatic pressure” option is applied at the inlet and in the initial values node for the water-filled domain 1, but not in the initial values node for the air-filled domain 2. A warning will appear in the solver, but it can be ignored.

Definitions

Variables 1

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
theta	(acos(tpf1.intnormr))[1/deg]		Contact angle expression

Meniscus Position

In the toolbar, click

and choose .

In the window for , type Meniscus Position in the text field.

In the text field, type z_pos.

Locate the section. Clear the checkbox.

Locate the section. Click

Select Boundary 6 only.

Click in the upper-right corner of the section. From the menu, choose > > .

Locate the section.

Select the checkbox. In the associated text field, type Contact point position.

Locate the section. From the list, choose .

Mesh 1

Mapped 1

In the toolbar, click

Size

In the window, click .

In the window for , locate the section.

From the list, choose .

Click

Study 1

Step 2: Time Dependent

In the window, under click .

In the window for , locate the section.

In the text field, type range(0,0.25e-4,1e-3).

Solution 1 (sol1)

In the toolbar, click

The automatic time stepping can have difficulties to start when the pressure initially tries to adjust to the surface tension. A manually specified time step is therefore a more appealing option in this case.

In the window, expand the node, then click .

In the window for , click to expand the section.

From the list, choose .

In the text field, type 5e-6.

In the text field, type 1e-3.

In the toolbar, click

The fifth default plot group shows the volume fraction of air. While the position of the air/water interface appears clearly, you can obtain an even sharper interface by plotting the 0.5 level of the same quantity using a filled contour plot, as in Figure 3.

Results

Contour 1

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Select the checkbox.

From the list, choose .

Surface 1

In the window, right-click and choose .

Click to confirm.

Volume Fraction of Fluid 1 (ls)

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Click the button in the toolbar, then zoom in on the lower part of the capillary. Compare the resulting plot with that in the upper-left panel of Figure 3.

Reproduce the remaining plots on the left in Figure 3 by plotting the solution for the time values 5e-5, 1e-4, and 1.5e-4.

Velocity (spf)

The second default plot shows a surface plot of the velocity magnitude of the fluids combined with one contour lines to identify the air/water-interface. This plot can be changed to reproduce the combined velocity-field arrows and air/water-interface plot shown in Figure 4.

Surface

In the window, expand the node.

Right-click and choose .

Click to confirm.

Arrow Surface 1

In the window, right-click and choose .

In the window for , locate the section.

Find the subsection. In the text field, type 30.

Contour 1