Spontaneous shrinkage of droplet on wetting surface in phase-field model

TL;DR

This paper investigates the spontaneous shrinking of droplets on wetting surfaces using phase-field modeling, revealing a critical radius dependent on contact angle and other parameters, below which droplets vanish.

Contribution

It introduces a theoretical prediction of the critical radius for droplet disappearance and verifies it through numerical simulations across various conditions.

Findings

Existence of a critical radius for droplet disappearance.

Critical radius depends on contact angle, domain size, and interface width.

Numerical results confirm the theoretical predictions.

Abstract

Phase field theory is widely used to model multi-phase flows. A drop can shrink or grow spontaneously due to the redistribution of interface and bulk energies to minimize the system energy. In this paper, the spontaneous behaviour of a drop on a flat surface is investigated. It is found that there exists a critical radius dependent on the contact angle, the domain size and the interface width, below which the droplet will eventually disappear. In particular, the critical radius can be very large when the contact angle is hydrophilic. The theoretical prediction of the critical radius is verified numerically by simulating a drop on a surface with various contact angles, the domain sizes and the interface widths.

Equations

$$\frac{\partial \varphi }{\partial t}=\nabla \cdot (M\nabla \mu ),$$

$$F={\int }_{\Omega }\left[f(\varphi )+\frac{\kappa }{2}|\nabla \varphi {|}^2\right]d\Omega +{\int }_{\partial \Omega }{f}_w(\varphi )d\partial \Omega ,$$

$$f_w(\varphi )=-\sigma \cos \theta \frac{\varphi (3-{\varphi }^2)}{4}+\frac{{\sigma }_{w1}+{\sigma }_{w2}}{2},$$

$$\mu =f^{\prime }(\varphi )-\kappa \Delta \varphi .$$

$$\frac{d}{dt}{\int }_{\Omega }\varphi d\Omega ={\int }_{\partial \Omega }M\nabla \mu \cdot \mathbf{n}d\Omega =0,$$

$$\frac{dF(\varphi )}{dt}$$

$${\mathbf{u}}_{\partial \Omega }$$

$$\mathbf{n}\cdot \nabla \mu {|}_{\partial \Omega }$$

$$(\kappa \mathbf{n}\cdot \nabla \varphi +f_w^{\prime }(\varphi )){|}_{\partial \Omega }$$

$$\delta \varphi =\frac{2\pi (R_0^3-R^3){\theta }_c}{3V},$$

$$\begin{array}{rl}F(R)& =\sigma S+{\sigma }_{w1}{d}_1+{\sigma }_{w2}(D-d_1)+\beta ({\varphi }_1^2-1{)}^2{V}_1+\beta ({\varphi }_2^2-1{)}^2(V-V_1)\\ & \simeq \sigma S+({\sigma }_{w1}-{\sigma }_{w2})d_1+{\sigma }_{w2}D+4\beta V\delta {\varphi }^2,\end{array}$$

$$F(r)=2\sigma \pi r(1-\cos \theta )-\sigma \pi r\cos \theta {\sin }^2\theta +{\sigma }_{w2}D+\frac{16\beta {\pi }^2(R_0^2+r^3-2r^{3/2}{R}_0^3){\theta }_c^2}{9V},$$

$$\frac{\partial F}{\partial r}=2\sigma \pi (1-\cos \theta )-\pi \sigma \cos \theta {\sin }^2\theta +\frac{16\beta {\pi }^2(3r^2-3\sqrt{r}{R}_0^3){\theta }_c^2}{9V},$$

$$\frac{{\partial }^{2}F}{\partial r^2}=\frac{16\beta {\pi }^2\left(6r-3R_0^3/2/\sqrt{r}\right){\theta }^2}{9V},$$

$$\frac{{\partial }^{3}F}{\partial r^3}=\frac{16\beta {\pi }^2\left(6+\frac{3R_0^3}{4r^{3/2}}\right){\theta }^2}{9V}.$$

$$R_c={\left(\frac{4}{{\cos }^3\theta -3\cos \theta +2}\right)}^{1/4}{\left(\frac{2^{-4/3}VW}{3\pi }\right)}^{1/4},$$

$$R_c={\left(\frac{\pi }{\theta -\sin \theta \cos \theta }\right)}^{1/3}{\left(\frac{\sqrt{3}VW}{16\pi }\right)}^{1/3}.$$

Full text

††preprint: APS/123-QED

In recent years, the phase-field model has received extensive attention in the fields of natural sciences and engineering applications, such as solidification and crystal growth Karma (2001); Rojas et al. (2015); Nestler et al. (2005), crack propagation Ambati et al. (2015); Pons and Karma (2010), multiphase flows Badalassi et al. (2003); Shen and Yang (2015) and moving contact line Ding and Spelt (2007); Yue and Feng (2011). Within this framework, the interface between two-phase fluids is treated as a thin layer of several grid sizes instead of a sharp interface. The interfacial profile can be identified by a suitably defined phase field variable, which is usually described by certain diffusive models such as the Cahn-Hilliard equation or the Allen-Cahn equation. A striking advantage of the phase-field model is free of explicitly tracking the interface between two-phase fluids, which leads to the phase field method easy to implement for interfacial dynamics in complex flows Badalassi et al. (2003); Teigen et al. (2011).

However, it is found that the size of a circular bubble can shrink to zero due to the dissipative mechanism in the phase field theory framework Liu and Shen (2003). Yue et al. Yue et al. (2004) noticed that the interface can shrink slightly even when the interface of a drop is initialized by the exact hypertangential profile. Similarly, Lee et al. Lee and Fischer (2006) observed that the radius of a droplet decreases first and then reaches a steady-state value in simulations. To understand this behavior, Yue et al. Yue et al. (2007) made a theoretical analysis of the Cahn-Hilliard model and found there existed a critical drop radius dependent on the domain size and interface width, below which the droplet will eventually disappear. Following the same procedure of Yue et al, zheng et al. Zheng et al. (2014) demonstrated that there is a similar critical radius for a pseudo-Vander Waals fluid in the single-fluid diffuse interface method. However, in both works, no solid walls are involved and the result is inapplicable to systems with solid bodies. In deed, it was found numerically that a drop on a solid surface can shrink with different rates for varying wetting properties Zheng et al. (2014). As the wetting of liquid on the solid surface is of great importance in defining the flow characteristics Bonn et al. (2009); Sohrabi et al. (2011), it is important to study the critical behavior of a droplet on a solid surface. In this paper we will present an such analysis in view of the Cahn-Hillard equation.

The Cahn-Hillard equation for an immiscible two-phase system can be written as

[TABLE]

where $\phi$ is the phase-field variable, i.e., $\phi_{1}=1$ in fluid $1$ and $\phi_{2}=-1$ in fluid $2$, $M$ is the mobility which is assumed to be constant in the present study and $\mu$ is the chemical potential defined from the mixing energy in the whole domain. For a two-phase system with a solid surface, The total mixing energy combines a bulk contribution from the Ginzburg-Landau free energy and a surface energy contribution from the substrate Cahn and Hilliard (1958),

[TABLE]

where $f(\phi)=\beta(\phi^{2}-1)^{2}$ is a double well potential, $\kappa$ is a constant and $\partial\Omega$ is the boundary of computational domain $\Omega$. Both $\beta$ and $\kappa$ are related to the interface tension $\sigma$ and the interface thickness $W$, i.e., $\beta=3\sigma/(4W)$, $\kappa=3/(8\sigma W)$. The wall energy $f_{w}$ can be given by Yue et al. (2010),

[TABLE]

where $\sigma_{w1}$ ($\sigma_{w2}$) is the fluid-solid interfacial tension between the wall and fluid 1 (fluid 2), $\theta$ is the static contact angle. As shown in Fig. 1, the contact angle satisfies the Young’s equation, i.e., $\sigma_{w2}-\sigma_{w1}=\sigma\cos\theta$,. The chemical potential $\mu$ is defined by

[TABLE]

For the Cahn-Hilliard equation, it can be easily identified that the total mass is conserved,

[TABLE]

and the total energy of the system is non-increasing over time,

[TABLE]

with the following boundary conditions,

[TABLE]

where $\bm{n}$ is the unit normal vector pointing into the wall. These two properties play important roles in analyzing the critical radius of a drop on a flat solid surface.

We consider a general case where a drop with initial radius $R_{0}$ is in contact with a flat solid surface in a domain with size $L\times L\times L$, as shown in Fig. 2. Following the same assumption of Ref.Yue et al. (2007), there is an equal and uniform shift $\delta\phi$ inside and outside of the drop when the drop radius shrink to $R$ from $R_{0}$. In addition, we further assume that the contact angle remains unchanged and the shape of the drop maintains a spherical cap during shrinking. With these assumptions, we can estimate the critical radius of the drop on the solid surface.

From Fig. 2, it can be seen that the volume of spherical cap is $V_{1}=\pi(3R-H)H^{2}/3$ with $H=(1-\cos\theta)R$ being the height of the spherical cap. According to the conservation property of the phase field variable, the shift of $\phi$ is

[TABLE]

where $\theta_{c}=(\cos\theta)^{3}-3\cos\theta+2$ and $V=L^{3}$ is the volume of the domain. Although $R_{0}^{3}-R^{3}$ is finite here, it can be confirmed that $\delta\phi\sim r_{c}^{3}/V\sim W^{3/4}/V^{1/4}\ll 1$. By neglecting cubic and quartic terms in $\delta\phi$, the system free energy $F(R)$ can be approximated by

[TABLE]

where $S=2\pi RH$ is the area of the surface of spherical cap, $d_{1}=\pi(\sin\theta R)^{2}$ is the area between fluid $1$ and the solid surface and $D$ is the total interface area between the fluids and the solid surface. As the surface energy of the drop in the three-dimensional (3D) space is proportional to area (for two-dimensional case, the surface energy is proportional to perimeter), the above equation can be rewritten as,

[TABLE]

where $r=R^{2}$. As shown in Fig. 3, the critical radius of the drop should correspond to an inflection point of the system energy, below which the energy decreases continuously. Thus we explicitly give the first to third derivatives of $F(r)$ with respect to $r$,

[TABLE]

[TABLE]

[TABLE]

Based on Eq. (15), a inflexion point $R_{i}$ can be determined from $\frac{\partial^{2}F}{\partial r^{2}}=0$, leading to $R_{i}=2^{2/3}R_{0}^{2}/4$. A vanishing drop requires $\frac{\partial F}{\partial r}|_{R=R_{i}}=0$. As a result, the critical radius is given by

[TABLE]

suggesting that the drop with $R_{0}<R_{c}$ will disappear eventually due to energy minization. The resulting formula is suitable to all contact angles (i.e., $\theta\in(0,\pi]$) although the contact angle in the above derivation is less than $\pi/2$. In particular, as $\theta=\pi$, the solid surface is completely non-wetting, the wall effect can be ignored, and the formula of the critical radius can be reduced to the result for a droplet without walls in Yue et al. (2007). Similarly, the critical radius of a two-dimensional (2D) drop can be determined as

[TABLE]

The above formula with $\phi=\pi$ is consistent with the one without walls in Ref. Yue et al. (2007) as well. The critical radius ($R_{c}$) profile of a drop with the products of the domain sizes and the interface widths $VW$ and contact angles $\theta$ based on Eq. (16) is shown in Fig 4. It can be observed that the critical radius is largely determined by the domain size and the contact angle.

We now demonstrate numerically the theoretical critical radius of a drop on the solid surface with various contact angles varying from $\pi/6$ to $\pi$ using the lattice Boltzmann method Zhang et al. (2019). A uniform cartesian grid is used and all simulation parameters are expressed in lattice units. First, we verify the relationship between the critical radius and the domain size as well as interface width. In the simulations, the domain is $L\times L$ for 2D simulations and $L\times L\times L/2$ for $3D$ simulations. The parameters $\sigma$ is fixed at $0.05$. Figure. 5 shows the results as $V=200$ and $300$ with $W=4$ and $6$. It can be seen that the numerical results are consistent with the theoretical predictions. Then, we test the theoretical formulas for both 2D drops and 3D drops with various contact angles. The results of 2D droplets are shown in Fig. 6. Again, good agreements with the theoretical predictions are observed.

To further illustrate the effect of contact angle on critical radius, a drop of radius $R=33$ in contact with the solid surface with $\theta=80^{\circ}$ or $100^{\circ}$ is simulated in a domain size of $300\times 300$. According to Eq. (17), the critical radiuses for $\theta=80^{\circ}$ and $100^{\circ}$ are $36$ and $31$, respectively. Thus, the drop with $\theta=80^{\circ}$ will vanish while the drop with $\theta=100^{\circ}$ will be maintained. The shrinkage processes for both contact angles are shown in Fig. 7. The contours of $\phi=0$ are plotted to indicate the interface locations. It can be observed that the drop with $\theta=80^{\circ}$ gradually vanishes while the drop with $\theta=100^{\circ}$ is close to a steady-state shape, which is in line with the theoretical prediction. As the Cahn-Hilliard dynamics tends to minimize the system energy, the equilibrium state of a drop should correspond to the minimum energy of the system. Thus, we also plot the energy curves of the above two drops, as shown in Fig. 8. As seen, the energy of the drop with $\theta=80^{\circ}$ finally decreases to zero while the total energy of the drop with $\theta=100^{\circ}$ approaches a non-zero stable value, which confirms the theoretical predictions Eq. (17) and Eq.(16).

In this paper, we carried out a theoretical analysis of the spontaneous shrinkage of a drop on a flat surface with a wetting boundary condition. It is shown that there exists a critical radius for the drop related to the domain size, interface width and the contact angle. All drops below the critical radius will eventually vanish. The critical radius could be rather large for a small contact angle for a given domain and interface width. Although drop shrinkage is very slow, it may be important for some practical applications with wetting phenomena, such as prediction of the relative permeability in porous media with the Cahn-Hilliard model .

We thank Pentau Yue for helpful discussions. This study was supported by the National Science Foundation of China(51836003).