A numerical study of three-dimensional droplets spreading on chemically patterned surfaces

We study numerically the three-dimensional droplets spreading on physically flat chemically patterned surfaces with periodic squares separated by channels. Our model consists of the Navier-Stokes-Cahn-Hilliard equations with the generalized Navier boundary conditions. Stick-slip behavior and contact angle hysteresis are observed. Moreover, we also study the relationship between the effective advancing/receding angle and the two intrinsic angles of the surface patterns. By increasing the volume of droplet gradually, we find that the advancing contact line tends gradually to an equiangular octagon with the length ratio of the two adjacent sides equal to a fixed value that depends on the geometry of the pattern.


1.
Introduction. The well-regarded Young's equation [23] gives the following relationship γ SL + γ LG cos θ e = γ SG , (1) where γ SL , γ LG and γ SG are solid-liquid, liquid-gas and solid-gas surface tensions respectively. It defines the balance of forces caused by a wet drop on a flat and chemically homogeneous dry surface. The drop will form a spherical cap with Young's (or intrinsic) contact angle θ e defined above. This shape corresponds to the unique minimizer of the free energy of the system by neglecting gravity and three-phase line tension.
However, the real surfaces are neither chemically homogeneous nor perfectly flat. Wenzel's model [20] and Cassie's model [4] are widely used to predict the average contact angle in a sense. In this case, the free energy functional has multiple local minimizers due to the inhomogeneities, not just the global minimizer prescribed by Young's equation (1). It is believed that this results in the contact angle hysteresis.
The equilibrium state (metastable or stable) of a droplet depends not only on the state parameters (such as γ SL , γ LG , γ SG and the volume of the droplet) but also on the history of the droplet because of the contact line pinning. Therefore the contact angle between the liquid and the solid phase will exhibit a range of possible values. The largest (or smallest) one is referred to as the advancing (or receding) contact angle denoted by θ adv (or θ rec ). The contact angle hysteresis is defined as the difference of these two angles ∆θ = θ adv − θ rec .
There are two commonly used methods for measuring the contact angle hysteresis. One is referred to as the tilting base method. Once a droplet is placed on the horizontal surface, one can tilt the surface from 0 • to 90 • . The contact angle on the downhill side will increase while the contact angle on the uphill side will decrease. As the tilt increases, the downhill side is going to wet the surface while the uphill side is going to dewet the surface. The angles just prior to the drop moving represent the advancing and receding contact angles. The other is referred to as the add and remove volume method. One can add volume to the droplet slowly, the droplet is going to spread. The contact angle prior to the spreading is termed advancing contact angle. Similarly, one can also remove volume from the droplet slowly, the droplet is going to shrink. The contact angle prior to the shrinking is termed receding contact angle.
Contact angle hysteresis has been investigated for many years, both analytically and experimentally. In two dimensional case, there have been lots of numerical and analytic results for different models [10,13,2,8,19,22,24]. The stick-slip behavior in two dimensional droplet or channel flow with striped patterned surfaces were studied. There have been also some numerical results and approximate solutions for three dimensional problems [9,7,16,17,15,3,21,18,1,11,12]. For heterogeneous surfaces, Cassie's model has its own constraint. There have been several modified models [21,18,8,12]. In [1,21], the authors reported that the average contact angle of the global equilibrium state approaches Cassie angle with increasing drop volume and the contact line is, on a large scale, circular in shape. Larsen and Taboryski [12] proposed a Cassie-like law using triple phase boundary line fractions to model the advancing contact angles of droplets spreading on chemically heterogeneous surfaces of the same hexagonal pattern of defects. They reported that drops tend to take hexagonal shapes due to the pinning forces of the hydrophobic lattice structure.
We investigate the contact angle hysteresis and the contact line behavior of threedimensional droplets on chemically patterned surfaces with periodic squares separated by channels. We obtain the following results. (a) There is also stick-slip behavior of contact line. (b) For the square-channel like pattern with equal area fraction 1/2, numerical evidence shows that the advancing contact line tends to be an equiangular octagon with length ratio √ 2 of the two adjacent sides as the size of the droplet is much larger than the size of the pattern. (c) The effective advancing contact angle is between the Cassie angle and the larger one of the two intrinsic contact angles of channel and square, and the receding contact angle is between the Cassie angle and the smaller one. Moreover the effective advancing contact angle is linearly dependent on the intrinsic contact angle of channel for fixed smaller intrinsic contact angle of square; the effective receding contact angle is linearly dependent on the intrinsic contact angle of square for fixed larger intrinsic contact angle of channel. This is a variant of the two-dimensional case.
The rest of the paper is organized as follows. The phase filed model and the numerical method are presented in Section 2. In Section 3, we present the numerical simulations and make some analysis of the results. The paper is finally concluded in Section 4.
2. The phase field model and the numerical method. We study droplets spreading on physically flat surfaces with periodic square chemical patterns with different contact angles. The contact angle and contact line tend to be an interesting figure as the volume of the droplet increases. Here we use a phase field model with the generalized Navier boundary condition proposed in [14]. The coupled system of Cahn-Hilliard-Navier-Stokes equations is as follows, (4) Here φ is the phase function with value 1 for the droplet range and −1 for the outer range, and the value changes rapidly from 1 to −1 during the interfacial region, p is the pressure, σ v = η(∇v + ∇v T ) denotes the viscous part of the stress tensor, ρ, η are the fluid mass density and viscosity, which will be assumed to be constant in the next section, ρg ext is the external body force density (zero here), and M is the phenomenological mobility coefficient; µ = −K∆φ − rφ + uφ 3 is the chemical potential, and µ∇φ is the capillary force; K, r, u are the parameters that are related to the interface profile thickness ξ = K/r, the interfacial tension γ = 2 √ 2r 2 ξ/3u, and the two homogeneous equilibrium phases φ ± = ± r/u (= ±1 here).
The boundary conditions at top and bottom walls are v n = 0, ∂ n µ = 0, . Γ is a positive phenomenological parameter, θ e (x) is the static contact angle, and β is the slip coefficient. The subscripts n, τ stand for the normal and tangent components respectively. Eqs. (6,7) represent the impermeability condition and the generalized Navier boundary condition respectively.
There are six dimensionless parameters (1) L d = 3M γ/2 √ 2V L 2 , (2) R = ρV L/η, the Reynolds number, (3) B = 3γ/2 √ 2ηV , it is inversely proportional to the capil- an interpolation between two different wall-fluid interactions, (6) ε, it is the ratio between interface thickness ξ and characteristic length L. Using the above parameters, we can obtain the dimensionless system (10) Here µ = −ε∆φ − φ/ε + φ 3 /ε is the chemical potential. The boundary conditions at top and bottom walls are v n = 0, and periodic conditions at the lateral boundaries for all variables. Here L(φ, on position x is due to the chemical pattern. The equilibrium state (metastable or stable) corresponds to the minimizer of the free energy functional Efficient numerical methods are developed to solve the above dimensionless system of equations [6]. where where , and the lateral boundary conditions In order to solve Eq. (15) with the boundary condition (18), we decompose the original equation into a system Based on the Neumann boundary conditions at the lateral boundaries, fast Fourier transform (FFT) can be applied to solve this linear system directly. Once φ n+1 is solved, Navier-Stokes equations become Helmholtz equations for velocity updating and Poisson equation for pressure, which can be solved by fast Poisson solvers, e.g. FFT in a standard way.
3. Numerical simulations and analysis of the results. We consider the squarechannel like patterned physically flat surface shown in Figure 1. Squares of side a are separated by channels of width b with intrinsic equilibrium contact angle θ a for squares and θ b for channels. Next we sometimes call θ a angle1 and θ b angle2. By Cassie's law [4], one can compute the Cassie angle θ C where the two area fractions f a = ( a a+b ) 2 , f b = b(b+2a) (a+b) 2 . If we choose the particular choice of a and b such that a b = √ 2 + 1, then f a = f b = 0.5. In this case, the Cassie angle will be the same if we swap the two angles θ a and θ b . As is known, Cassie angle is an average of the local contact angles in a sense. When a droplet is placed on this patterned surface, it will not form a spherical cap any more because of the different contact angles for the squares and channels. One might simply expect that it will not change the behavior of droplets to swap the contact angles for the squares and channels. This is not true unless the strength of the inhomogeneity is weak to some extent [9,7] which implies negligible contact angle hysteresis. 3.1. Advancing and receding motion. We first consider the advancing motion of a droplet spreading on the two types of chemically patterned surfaces. Let θ a = 60 • , θ b = 110 • for type 1 surface, and θ a = 110 • , θ b = 60 • for type 2 surface. In two dimensional case [22,11], a droplet is placed on the chemically patterned stripes with two different intrinsic contact angles θ 1 < θ 2 . Assume that the static contact line is in the hydrophilic region initially. When adding volume to the droplet gradually, the contact line first moves outward until to the borderline of the hydrophilic and hydrophobic regions in order to keep the smaller contact angle θ 1 . This is the socalled slip motion. With increasing volume, the contact line pins to the borderline before it moves to the hydrophobic region until the contact angle attains the larger one θ 2 which results in the advancing angle. After that, the contact line slips in the hydrophobic region and the contact angle keeps constant θ 2 until the contact line moves to the borderline of hydrophobic and hydrophilic regions. It may lower the free energy by immediately spreading the droplet to the hydrophilic region with contact angle θ 1 or to the next chemical borderline, whichever is first reached.
In three dimensional case, we also observe the stick-slip motion of contact line. Distortion occurs in the contact line due to the chemical inhomogeneity which connects the hydrophilic and hydrophobic regions. The section of contact line in the hydrophilic region is convex while that in hydrophobic region is concave. The droplet no longer forms a strict spherical cap (shown in Figure 2). Therefore the local contact angle depends on the position of contact line. The macroscopic contact angle of the droplet is not well-defined. However, the interface is close to a sphere away from the base surface. We could define an effective contact angle θ ef by applying the popular method used in [11,5] cos θ ef = 1 − Hκ, where H is the height of the droplet, κ is the curvature at the top of the droplet.  Initially, a droplet of a spherical cap (shown in Figure 3) is placed on the patterned surface, with angle θ, base radius r and directed distance h from the base to the center of the sphere. The volume of this droplet is or V r (θ) = πr 3 (1 − cos θ) 2 (2 + cos θ) 3 sin 3 θ , θ ∈ (0, π).
V r (h) (or V r (θ)) is an increasing function of h (or θ).
In our simulations, we first choose ε = 0.01, ∆x = ∆y = 0.01, ∆t = 0.1∆x, (a, b) = ∆x (12,5) such that a b ≈ √ 2 + 1 and take the following parameters as used in [14,6]: L d = 5.0 × 10 −4 , R = 3, B = 12, V s = 500, l s = 0.0038. (29) The stability parameter s is taken to be s = 1.5 andα will be determined by the two intrinsic angles of the square and channel. The initial conditions are v(x, y, z, 0) = 0, φ(x, y, z, 0) = tanh where (x 0 , y 0 ) are the (x, y) coordinates of the center of initial spherical droplet and r, h are the parameters defined in Eq. (27). In Figure 4, we compare the results between ε = 0.01, ∆x = ∆y = 0.01, ∆t = 0.1∆x, (a, b) = ∆x (12,5) and ε = 0.005, ∆x = ∆y = 0.005, ∆t = 0.1∆x, (a, b) = ∆x (24,10). The equilibrium states are almost the same. Next it is enough to choose ε = 0.01. Figure 5 shows the advancing motion of droplets placed on type 1 surface with initial droplets centered at channel. The two intrinsic contact angles of squares and channels are 60 • , 110 • respectively. The horizontal line corresponds to the Cassie angle 85.5 • . The contact line pins in the horizontal and vertical directions before it jumps to the next hydrophilic square. For the effective contact angle, it reaches the largest one 98.0 • before the droplet spreads to the next hydrophilic square with effective contact angle 80.5 • . The jump of contact angle is 17.5 • . The contact line pinning is due to the free energy barrier of the system. In the diagonal directions, the contact line moves more easily. The second smaller jump of the contact angle is due to the contact line movement in these diagonal directions. The results we present in Figure 5 show the stick-slip motion of contact line on type 1 surface. Figure 6 shows the advancing motion of droplets placed on type 1 surface with initial droplets centered at square. From the simulation, we see that the results are very similar to that of the above case, with droplets centered at channel. The effective contact angle first increases to the largest one 98.0 • and then jumps to 80.6 • . Later on, a second smaller jump occurs. The contact line sticks to the outer edges of squares in the horizontal and vertical directions until it jumps to the next hydrophilic square. Next a smaller jump of contact line occurs in the diagonal directions. The effective advancing angle is 98.0 • for both cases with initial droplet centered at square and channel. Figure 7 shows the advancing motion of droplets placed on type 2 surface with initial droplets centered at channel. The two intrinsic contact angles of squares and channels are 110 • , 60 • respectively. The contact line is different with that in Figure 5. The effective advancing angle is 92.3 • which is smaller than that in Figure 5.
In the above cases, the initial droplets are placed symmetrically on the patterned surfaces. The contact line motions are different for surface 1 and surface 2. Moreover the effective advancing angles are between the Cassie angle and the larger one of the two intrinsic angles. For other asymmetric cases, the contact line of the equilibrium state may not display as a symmetric curve. It is much more complicated. We do not deal with these initial cases.
Removing volume from the droplet quasistatically gives the receding motion of contact line. Figure 8 and Figure

3.2.
Asymptotic behavior of advancing motion for large volume. We now study the asymptotic behavior of the effective contact angle and the contact line when the volume of droplet increases. By adding more volume to the droplet gradually, we find that in Figure 10 the largest contact angles before jump are almost the same and that in Figure 11 the contact line follows the same stick-slip behavior. The four pairs of effective contact angles before and after jump in Figure 10  In each period, the largest angles are almost the same, these common angle is referred to the effective advancing angle.
Continue adding volume to the droplet, we consider six more periods. The effective contact angle is shown in Figure 12, and contact line plot in Figure 15. Again, the largest angles before jump almost keep constant while the effective contact angles after jump increase with volume. The six more pairs of effective contact angles before and after jump in Figure 12 are Next we combine the above ten periods, and plot the dependence of effective contact angle on the cubic root of volume instead of the volume itself (shown in Figure 13). Now one can view the whole picture. Numerical evidence shows that the effective advancing angles hold the same while the jump decreases with volume. From Figure 13, one could guess that the jump of effective contact angle decreases to zero when the volume tends to infinity which is the same with that of two dimensional case. The following is an explanation of this guess.
Recall that Eq. (28) gives the dependence of volume on base radius and angle of a spherical cap, V r (θ) = πr 3 (1 − cos θ) 2 (2 + cos θ) 3 sin 3 θ , θ ∈ (0, π).  We could use the distance between the center of the base and the horizontal edge of the pinning contact line to approximate the base radius before and after jump r k,be = (k + 2)(a + b), r k,af = (k + 3)(a + b) in the k-th period. Notice that lim k→∞ r k,be r k,af = lim We conclude that lim k→∞ θ k,be θ k,af = 1, where θ k,be is the effective advancing contact angle θ adv which keeps constant. Therefore the effective contact angle just after the jump tends to the effective advancing contact angle, i.e. the jump of effective contact angle approaches zero when the volume tends to infinity. This property is the same with that of two dimensional case. From numerical simulations, Figure 14 shows that the cubic root of volume at jump forms an arithmetic sequence. The fitting equation of these points is V olume 1/3 = 0.250 P eriod + 0.524. ( This also shows that the droplet spreads periodically in another aspect. In Figure 15, the contact line shows some periodic property. Now we view the innermost two groups of contact lines as a whole. The innermost contact line pins to the horizontal and vertical edges which cover 6 periods of pattern. The horizontal and vertical edges are connected by the diagonal section which covers 4 squares. Due to symmetry, we only deal with the north and the northeastern edges. With a slight volume increase, the north pinning section holds while the central part in channel of the northeastern diagonal section jumps to cover the outer square partially. Next the diagonal part totally jumps to the outer squares. At the same time, the north pinning part covers 8 squares. With more volume added, the pinning section jumps to cover the same number 6 outer squares. This is the main jump after the weak jump in diagonal direction. While in the diagonal direction, only the central part section jumps to the outer 4 squares first, and then the diagonal section reproduces the picture between the fourth and fifth main jumps. In Figure 11, the contact line with the first 4 main jumps has no such property. This is because the size of droplet is not so much larger than the size of pattern. There is a competition between the pinning and the diagonal sections during the jump which means free energy release. Following the trend, before the (2k − 5)-th main jump (k = 3 represents the innermost pinning and jumping process in Figure 15), the diagonal section of contact line displays as covered by (k + 1) squares, intermediate weak jump -the central part jumping to the outer (k − 2) squares, and then covered by k squares. During the (2k − 5)-th main jump, the central part of diagonal section jumps to the outer (k − 1) squares. The diagonal section of contact line then reproduces the picture before the (2k − 5)-th main jump.
The length ratio of north and northeastern edges approaches Therefore, with equal area fraction 1/2, the contact line approximates to an equiangular octagon with length ratio √ 2 of pinning section and diagonal section when the size of droplet is much larger than the size of pattern. With other area fractions, it is highly probable that the contact line still approximates to an equiangular octagon. Here the pinning section of the contact line is almost a straight line due to the energy barrier, while the diagonal section of the contact line connecting the adjacent pinning sections displays as a curve. Therefore the equiangular octagon is only an approximate figure for the contact line of the advancing motion of droplet with large volume.
3.3. Linear dependence of advancing and receding angles on two intrinsic angles. Now we consider the dependence of effective advancing and receding contact angles on the two intrinsic angles. In two-dimensional case, the advancing angle is exactly the larger one of the two intrinsic angles, it has nothing to do with the smaller intrinsic angle; the receding angle is exactly the smaller one of the two intrinsic angles, it has nothing to do with the larger intrinsic angle. However in three-dimensional case we considered here, numerical evidence shows that the effective advancing/receding angle is between the Cassie angle and the larger/smaller intrinsic angle. Moreover the effective advancing angle is linearly dependent on the larger intrinsic angle of channel for fixed smaller intrinsic angle of square. The effective receding angle is linearly dependent on the smaller intrinsic angle of square for fixed larger intrinsic angle of channel. Figure 16 shows that for fixed intrinsic angle of square 60 • , the effective advancing angle is linear on the intrinsic angle of channel. We have done numerical simulations for the cases of angle of channel ranging from 60 • to 150 • with step 2.5 • . The effective advancing angle is linearly dependent on the intrinsic angle of channel. The linear fitting equation is If we take θ b = θ a = 60 • , it is the chemically homogeneous case. The advancing angle is exactly the intrinsic angle 60 • . The advancing angle predicted by the above fitting equation is 57.4 • . There is some difference. This is due to the error of the approximate approach to calculate the effective contact angle and the error of diffuse interface model to capture the interface.  Figure 17 shows that for fixed intrinsic angle of channel 120 • , the effective receding angle is linear on the intrinsic angle of square. We have done numerical simulations for the cases of angle of square ranging from 60 • to 120 • with step 2.5 • . The effective receding angle is linearly dependent on the intrinsic angle of square. The linear fitting equation is θ rec = 0.661 θ a + 39.5 • , θ a ∈ [60 • , 120 • ].
If we take θ a = θ b = 120 • , it is the chemically homogeneous case. The receding angle is exactly the intrinsic angle 120 • . The receding angle predicted by the above fitting equation is 118.8 • .

Conclusion.
We study three-dimensional droplets spreading on the chemically patterned surfaces with periodic squares separated by channels. Numerical simulations show that there are also stick-slip motions for three-dimensional problems.
In two-dimensional case, the advancing/receding angle is exactly the larger/smaller one of the two intrinsic contact angles of the chemically patterned stripes. In threedimensional case, we also find a linear dependence of the effective advancing (or receding) angle on the larger (or smaller) one of the two intrinsic angles with the other fixed. Moreover for droplet with large volume, the advancing contact line approaches an equiangular octagon with ratio √ 2 of the length of the two adjacent sides for the pattern with equal area fraction 1/2 as the size of the droplet is much larger than the size of the pattern.