Implementing Contact Angle Hysteresis in Moving Mesh-Based Two-Phase Flow Numerical Simulations.

Contact angle hysteresis is a common phenomenon in nature, which also plays an important role in industrial applications. A numerical model based on the moving mesh two-phase flow method is presented for modeling contact angle hysteresis. The implementation includes a displacement-based penalty method and a state variable method. The pinning, moving, and repinning of the contact lines can be simulated. This method is robust considering both two-dimensional and three-dimensional geometries. To further demonstrate the performance of this method, a fluid-solid interaction model with a cylinder fluctuating on a water surface considering contact angle hysteresis is demonstrated.

Introduction

The contact angle (CA) for a liquid droplet on a solid surface is commonly used to characterize the wettability of the solid surface.[1] The classical Young’s equation provides a unique CA when the three-phase system is in a static state with the lowest global surface energy. However, solid surfaces can hardly be perfectly homogeneous and smooth because most of them are covered with natural defects, rough structures, or artificial patterns, where the apparent contact angle is generally different from the equilibrium CA.[2−4] The CA when the interface is moving forward (called the advancing CA) is larger than that when moving backward (called the receding CA). The difference between the advancing CA and the receding CA is defined as contact angle hysteresis (CAH).[1,4−6] CAH is one of the most common phenomena in nature, which also plays an important role in industrial applications such as immersion lithography, fiber coatings, and inkjet printing.[2] Recently, many digital microfluidics[7−9] take advantage of surface wetting phenomena to achieve droplet manipulations and drug delivery. Numerical simulation is a powerful tool to help design these microfluidic devices. However, because of the complexity of CAH modeling, the CAH is usually ignored in these simulations for microfluidic designs.[10] Thus, developing numerical methods to model the CAH phenomenon is an important topic in computational fluid dynamics (CFD).

The simulations of droplet motion require a two-phase flow model, where the shape of the droplet interface is explicitly presented. The finite volume method or finite element method CFD solvers for these kinds of two-phase flow can be classified into two types. The first type is field-based methods including the volume of fluid (VOF) method,[11] level-set method,[12,13] or phase-field method,[14] which distinguishes the two-phase flow using an additional field variable. The second type is the moving mesh method,[15] which uses the deforming mesh to describe the shape of the interface. In some works,[16−19] the authors implement CAH using the first type of two-phase flow model. Fang et al. carried out a transient model by correcting boundary force balances through specifying the local CA and instantaneously updating the local angle values based on the variation of the volume fraction from previous time steps.[16] Dupont et al. carried out the numerical method based on the implementation of a “subgrid” description of the contact line that consists in imposing the apparent angle for static and moving contact lines.[17] The model reported by Ahmed et al. automatically locates the advancing and receding sections of the contact line, which enables the application of different contact angles at the advancing and receding fronts and therefore takes into account CAH.[18] In the model of He et al., a pair of pseudoline tensions in the receding and advancing states, respectively, are utilized to represent contact line interactions with the substrate.[19]

Here, we provide a method to implement CAH in a moving mesh-based two-phase flow CFD solver. With the moving mesh approach, the phase boundary is modeled as a geometrical surface separating two domains, with one phase at each side in the corresponding domains. The shape of the interface is more accurate compared to field-based methods. Forces and fluxes can be directly applied on the phase boundary. An advantage of the moving mesh method is that it can model one single phase while the other phase can be ignored. The drawback of the moving mesh model is that it cannot deal with topology changes.[15,20] This work is based on a commercial CFD solver provided by COMSOL Multiphysics (abbreviated as COMSOL). The “two-phase flow moving mesh” interface in COMSOL only provides the condition of a simple CA that follows the Young’s equation.[21,22] This work provides a method to implement CAH in COMSOL. The core idea of this method is based on the penalty method.[23,24] The advancing and receding CA are the input parameters in the simulations. The pinning, moving, and repinning of the contact line can be simulated. This method is robust in simulations considering both two-dimensional (2D) and three-dimensional (3D) geometries. To further demonstrate the performance of this method, we conduct a fluid–solid interaction (FSI) simulation, where a cylinder is fluctuating on a water surface, where the contact line exhibits CAH.

Results and Discussion

Numerical Models and Methods

The main characteristic of CAH is that the contact line is pinning at the initial position when the apparent CA is less than the advancing CA and larger than the receding CA. The contact lines begin to move as the apparent CA reaches advancing CA (θa) or receding CA (θr). The relation between the apparent CA and the displacement of the contact line is plotted in Figure a. We use a penalty method to describe the pinning characteristics. First, we use a 2D geometry to demonstrate this method. As shown in Figure b, a droplet is placed on a surface with a channel at the bottom to inject or remove liquid to or from the droplet. We define the input CA (θin) as a function of the displacement (d) of the contact point, i.e.,where θ0 is the initial CA and p is the penalty parameter being positive. The displacement (d) is defined aswhere ims is the unit surface binormal vector as shown in Figure b, (x0, y0) is the initial position of the contact point, and (x1, y1) is the current position of the contact point after moving. When the contact point is advancing, the value of d is positive and θin will become larger to prevent the contact point from advancing any further. When the contact point is receding, the value of d is negative and θin will become smaller to prevent the contact point from receding any further. When p is big enough, the value of d is very close to zero as the value of CA changes from θa to θr. Thus, eq makes the contact point pinning at the initial position (x0, y0). Based on eq , we construct an inequality relation to constrain the maximum and minimum values of θin.Then, θin will not become larger than θa or smaller than θr. The contact point will start to move as θin reaches θa or θr. We plot eq in Figure c with different values of p while the values of θ0, θa, and θr are 90, 100, and 80°, respectively. When p is big enough, the curve is nearly the same as in Figure a. Thus, eq can be used to simulate CAH.

Figure 1

Numerical method demonstration in 2D. (a) CA versus contact line displacement of CAH. (b) Modeling scheme with the black curve indicating the initial shape of the droplet and ims is the unit surface binormal vector. (c) Plots of eq with different p values. (d) Velocity distributions and shape deformations of the droplet. (e) Simulation results of contact point moving velocity and Mv corresponding to panel (d). (f) Simulation results of θ0 and θin corresponding to panel (d).

Some refs (2, 17, 25) reported that the apparent CA varies with the moving velocity of the contact point (vc). To include this effect, an additional velocity-dependent function, f (vc), can be added to describe this dynamic CA.The wetting behavior varies for different liquids and different substrates. It depends on the experimental observations whether the velocity-dependent CA requires to be considered in the simulations. In this work, we only use eq , where the velocity-dependent CA change is ignored.

When the CA reaches θa or θr, the contact point starts to move. Then, when the contact point stops, it will be pinned again in a new position. To describe the repinning process, the coordinates of the position (x0, y0) need to be updated. To simplify the formulas, we define a bool variable, Mv.

Mv equals 1 when the contact point is moving and 0 otherwise. Then, the position (x0, y0) is defined as a state variableWhen the contact point is pinning, the value of (x0, y0) keeps unchanged. When the contact point is moving, the pinning position is updated to be the current position. To avoid the eruptive change of θin, θ0 also needs to be defined as a state variable and is updated when Mv is 1.Using the above settings, we simulate such a process that the water is first drained out of the droplet and then injected into the droplet. The governing equations for the incompressible fluid flow arewhere u is the velocity, ρ is the density, ps is the pressure, μ is the dynamic viscosity, and t is the time. The equation for the droplet surface iswhere σ is the surface tension and n is the surface normal vector. The moving mesh provided by COMSOL will automatically deform the mesh to adapt the shape of the droplet. The numerical solver in COMSOL will apply a force at the contact point to maintain the apparent CA of the droplet to be the input CA. As for the example in Figure d, the initial CA is 60°. The advancing and receding CA is 70 and 50°, respectively. The initial distance between the two contact points is 5 mm. The material of the droplet is water, with its parameters taken from the built-in library in COMSOL. The instructions for implementing eqs 1–7 in COMSOL are provided in Figures S1–S3. The simulation results are shown in Figure d and Movie S1. Figure e shows the time-dependent plots of contact point moving velocity and the bool variable, Mv. Figure f shows the time-dependent plots of θ0 and θin. The contact point goes through four processes, namely, pinning, moving, repinning, and moving. In the first process, Mv is 0 and θ0 keeps unchanged while θin decreases. In the second process, Mv turns into 1 and θ0 turns into θr since the contact point moves. In this process, the state variables θ0 and (x0, y0) are updated in every numerical time step. As θin cannot be smaller than θr, the value of θ0 is always θr. The value of (x0, y0) is the output value of (x1, y1) in the last time step of the current process and is the input parameter for the simulation in the next time step. If the output value of (x1, y1) makes a negative value of d in the next time step, θ0 + pd will be less than θr. Then, Mv is still 1. If the output value of (x1, y1) makes a positive value of d in the next time step, θ0 + pd will be larger than θr, which means that the moving direction of the contact point is changed. Then, Mv is 0, and the state variable stops updating. Thus, the contact point is pinning again, which is described as the third process. In the fourth process, θ0 + pd is larger than θa, and the contact point is moving again. Thus, the CAH is successfully simulated.

3D CAH Simulations

Here, we demonstrate two 3D CAH examples. In 3D, the contact points become contact lines. Every contact point at the contact line has a value of CA. The implementing methods in 3D and 2D are similar. The only difference is that the coordinates of contact points and the vectors have three components in 3D.

The first example is a hemispheric droplet sliding down a ramp. The ramp has a gradient of 45°. The acceleration of gravity is 9.8 m/s2. The radius of the droplet is 1 mm. The initial CA is 90°. The density and viscosity of the liquid are 998.2 kg/m3 and 0.02 Pas, respectively. The surface tension is 0.01 N/m. The advancing and receding CAs are 110 and 80°, respectively. The initial shape of the droplet is hemispherical, and its deformation and movement are shown in Figure a and Movie S2. Figure b shows the time-dependent plots of the advancing CA and receding CA. The instructions for implementing CAH of this example in COMSOL are provided in Figures S4–S6. Figure c shows the profile change of the contact line over time. The results are qualitatively agreeable with experimental observations in the refs (18, 19).

Figure 2

Modeling of a sliding droplet: (a) Simulation results with the color presenting the velocity, (b) advancing CA and receding CA over time, and (c) contact line profiles over time.

The second example is the rising and falling of a liquid surface in a tube with an elliptical cross section. The long axis and short axis of the ellipse are 6 and 2.6 mm, respectively. The initial CA is 90°. The advancing and receding CAs are 110 and 70°, respectively. The density and viscosity of the liquid are 998.2 kg/m3 and 0.01 Pas, respectively. The surface tension is 0.07 N/m. The inlet and outlet velocity are both 4 mm/s. The instructions for implementing CAH of this example in COMSOL are provided in Figures S7–S9. The results are presented in Figure a and Movie S3. As shown in Figure a, we define two points, namely, A and B at the contact line. The time-dependent distributions of variable Mv, height of the contact line, and CA along the contact line from A to B are plotted in Figure b–d, respectively. The abscissa is the arc length and the ordinate is time, and the variable values are presented by colors. Because of the elliptical shape, it is easier to reach point B, advancing CA and receding CA. Thus, point B starts moving first and then point A finally starts to move. These two examples show that this method can be well applied to 3D cases.

Figure 3

Modeling of a liquid surface rising and falling in an ellipse tube: (a) Simulation results with the color presenting the height of the surface, (b) Mv distributions along arc AB, (c) surface height distributions along arc AB, and (d) apparent CA distributions along arc AB.

Fluid–Solid Interaction Simulation with CAH

To further demonstrate the performance of this method, we conduct an FSI simulation with CAH. Compared with a previous simulation, the difference here in the FSI simulation is that the solid surface also moves. The relative displacements need to be calculated. Here, we simulate a cylinder floating on the water surface (Figure a). The density of the cylinder is 600 kg/m3. The cylinder is pressed into the water at the beginning. After releasing, the cylinder is floating up and down. The initial CA is 90°, with the advancing and receding CAs of 100 and 85°, respectively. In the simulation, 2D axisymmetric space dimension is used. The cylinder is considered as a rigid body. The cylinder is subjected to the gravitational force, fluid force, and surface tension force. The acceleration of the cylinder (as) can be written aswhere g is the acceleration of gravity, T is the fluid force acting on the solid surface per area calculated by COMSOL, θin is the apparent CA, ms is the mass of the cylinder, and σ is the surface tension. The instructions for implementing eq in COMSOL are provided in Figure S10. The velocity of the cylinder (vs) is

Figure 4

FSI simulation: (a) Cylinder floating on the water surface, (b–d) simulation results with the color presenting the velocity, (e) simulation results of the velocities of the cylinder (vs) and the contact point (vp), and (f) simulation results of θ0 and θin.

The displacement of the cylinder (Ds) is

The instructions for implementing eqs 12 and 13 in COMSOL are provided in Figure S11. Then, Ds is used to control the deformation of the fluid domain. The instructions for implementing the deformation of the fluid domain in COMSOL are provided in Figure S12. In this model, the cylinder can only move in a vertical direction (z direction). We take an arbitrary point in the cylinder to trace its movement. The initial z coordinate of this point is zs0, and the current z coordinate is zs1. The displacement of the cylinder (ds) can be defined as

Then, the displacement of the contact point on the cylinder surface (dc) iswhere zc0 is the initial z coordinate of the contact point and zc1 is the current z coordinate of the contact point. Then, the input CA can be defined.To simulate the repinning of the contact point, zc0, θ0, and zs0 need to be defined as state variables.In FSI simulations, the displacements of the solid surface also need to be updated when depinning occurs. The instructions for implementing eqs 14–19 in COMSOL are provided in Figures S13–S15.

Using the above settings, the FSI with CAH is simulated. The results are shown in Figure b–d and Movie S4. Figure e shows the moving velocity of a cylinder (vs) and the contact point (vp). When the two curves are overlapping, the relative velocity between the cylinder and contact point is zero, indicating the pinning of the contact point at the cylinder surface. The time-dependent values of θ0 and θin are shown in Figure f. When the contact point is pinning, the values of θ0 and θin are different. From these results, we can see that FSI with CAH is successfully simulated.

Conclusions

In this paper, a method of implementing CAH in a moving mesh-based two-phase flow model has been presented and verified with several examples. This method is based on the displacement of the contact lines. It takes advantage of easily obtaining the coordinates of the contact lines in the moving mesh method. The penalty method is used to pin the contact line at a given position. The maximum and minimum values of the input CA are constrained to be the advancing and receding CAs. The repinning of the contact line is realized using the state variables to update the pinning positions. This method is robust when dealing with both 2D and 3D geometries and can be implemented in FSI simulations.

Some previous works[26,27] reported the velocity-based methods to describe the dynamic CA, where the CA is a function of the moving velocity of the contact line, which increases with the advancing velocity and decreases with the receding velocity. CAH can be mimicked by these methods. However, this method is difficult to model static CAH when the velocity of the contact line is zero. The velocity-dependent CA only gives one specific value when the velocity of the contact line is zero, which cannot describe the CA change when the contact line is pinning. Our displacement-based method can model static CAH when the contact line is pinning, and the velocity-based dynamic CA can also be incorporated in our method as shown in eq when the contact line is moving. From a numerical calculation point of view, penalty methods based on the contact line displacement are usually easier to solve and have better convergence than those based on the contact line velocity, as velocity is the derivative of displacement with respect to time. Moreover, although this method is implemented in COMSOL in this work, it is possible to be applied in other CFD solvers.

In addition, the reported methods[3,16−19] for CAH are unfriendly for researchers who are not familiar with the sophisticated CFD solver and programming, which cannot benefit many researchers from applying them in microfluidics. This work provides a method to simulate CAH through the commercial software COMSOL without any coding. With the detailed modeling setups provided in the Supporting Information, this method can be replicated by other researchers easily.