Effective boundary conditions for dynamic contact angle hysteresis on chemically inhomogeneous surfaces

Consider a region $\Omega$ near the contact line as in Figure 1. For simplicity, we suppose $\Omega$ is a two-dimensional domain. The boundary of $\Omega$ is composed of two parts, the lower solid boundary $\Gamma_{S}$ and the outer boundary $\Gamma_{O}$. On $\Gamma_{S}$, there exists a moving contact line, which is an intersection point between the solid boundary and the two-phase flow interface $\Gamma$. Then $\Omega$ represents an open system near the contact point. In the following, we will derive a sharp interface continuum model for moving contact lines. We adopt the Onsager principle to derive the model. The method has been used to develop the generalized Navier slip boundary condition for a diffuse interface model in Qian et al. (2006). In the derivation, we ignore the gravitational force and the inertial effect, which can be added simply if needed.

To use the Onsager principle, we first compute the rate of change of the total energy in the system. The total free energy consists of three interface energies,

$$\mathcal{E}=\int_{\Gamma_{SA}}\gamma_{SA}ds+\int_{\Gamma_{SB}}\gamma_{SB}ds+\int_{\Gamma}\gamma ds,$$

(2.1)

where $\Gamma_{SA}$ and $\Gamma_{SB}$ are respectively the parts of solid boundary in contact with the fluid $A$ and fluid $B$, $\gamma_{SA}$ and $\gamma_{SB}$ are corresponding surface tensions, and $\gamma$ is the surface tension of the two-phase flow interface. Direct calculations give

$\displaystyle\dot{\mathcal{E}}$

$\displaystyle=\gamma(\cos\theta_{d}-\cos\theta_{Y})v_{ct}+\gamma\int_{\Gamma}v_{n}\kappa ds.$

(2.2)

Here $\theta_{d}$ is the dynamic contact angle with respect to the fluid $A$, $\theta_{Y}$ is Young’s angle, $v_{ct}$ is the velocity of the contact line, $v_{n}=\mathbf{v}\cdot\mathbf{n}$ is the normal velocity, and $\kappa$ is the curvature of the interface. In the derivation, we have also used the Young’s equation $\gamma\cos\theta_{Y}=\gamma_{SB}-\gamma_{SA}$.

The energy dissipation function, which is defined as half of the total energy dissipation rate in the system, can be written as

$$\Phi=\int_{\Omega_{A}}\frac{\mu_{A}}{2}|\nabla\mathbf{v}|^{2}dx+\int_{\Omega_{B}}\frac{\mu_{B}}{2}|\nabla\mathbf{v}|^{2}dx+\int_{\Gamma_{SA}}\frac{\beta_{A}}{2}v_{\tau}^{2}dx+\int_{\Gamma_{SB}}\frac{\beta_{B}}{2}v_{\tau}^{2}dx+\frac{\xi}{2}v_{ct}^{2},$$

(2.3)

where $\Omega_{A}$ and $\Omega_{B}$ are the regions in $\Omega$ occupied by fluid A and fluid B, respectively, $\mathbf{v}$ is the corresponding velocity field, $v_{\tau}$ is the slip velocity on the solid boundary, $\mu_{A}$ and $\mu_{B}$ are viscosity parameters, $\beta_{A}$ and $\beta_{B}$ are phenomenological slip coefficients, and $\xi$ is the friction coefficient of the contact line. The normal velocity on the solid boundary $\Gamma_{S}$ is zero.

Since the system is an open system, we need also consider the work to the outer fluids at the boundary $\Gamma_{O}$. It is given by

$$\dot{\mathcal{E}}^{*}=-\int_{\Gamma_{O}}\ \mathbf{F}_{ext}\cdot\mathbf{v}ds=-\int_{\Gamma_{OA}}\mathbf{F}_{ext}\cdot\mathbf{v}_{A}ds-\int_{\Gamma_{OB}}\mathbf{F}_{ext}\cdot\mathbf{v}_{B}ds,$$

(2.4)

where $\Gamma_{OA}$ and $\Gamma_{OB}$ are the respective open boundary in contact with fluid A and fluid B.

With the above definitions, the Onsager principle states that (Doi, 2013) the dynamical equation of the system is given by minimizing the Onsager-Machlup action defined as follows

$$\mathcal{O}=\dot{\mathcal{E}}+\dot{\mathcal{E}}^{*}+\Phi,$$

(2.5)

under the constraint of incompressibility condition

$$\nabla\cdot\mathbf{v}=0.$$

(2.6)

Introduce a Lagrangian multiplier $p(x)$ in $\Omega$. We minimize the following modified functional with respect to the velocity

$$\mathcal{R}_{\lambda}=\dot{\mathcal{E}}+\dot{\mathcal{E}}^{*}+\Phi-\int_{\Omega_{A}\cup\Omega_{B}}p\nabla\cdot\mathbf{v}dx.$$

(2.7)

Direct calculation for the first variation of $\mathcal{R}_{\lambda}$ gives

	$\displaystyle\delta\mathcal{R}_{\lambda}=$	$\displaystyle\int_{\Gamma}\gamma\kappa\delta v_{n}ds+\gamma(\cos\theta_{a}-\cos\theta_{Y})\delta v_{ct}-\int_{\Gamma_{O}}\mathbf{F}_{ext}\cdot\delta\mathbf{v}ds+\int_{\Omega_{A}}\mu_{A}\nabla\mathbf{v}\cdot\nabla\delta\mathbf{v}dx$
		$\displaystyle+\int_{\Omega_{B}}\mu_{B}\nabla\mathbf{v}\cdot\nabla\delta\mathbf{v}dx+\int_{\Gamma_{SA}}\beta_{A}v_{\tau}\delta v_{\tau}ds+\int_{\Gamma_{SB}}\beta_{B}v_{\tau}\delta v_{\tau}ds+\xi v_{ct}\delta v_{ct}$
		$\displaystyle-\int_{\Omega_{A}\cup\Omega_{B}}p\nabla\cdot\delta\mathbf{v}dx.$

Integration by parts gives

	$\displaystyle\int_{\Omega_{A}}\mu_{A}\nabla\mathbf{v}\cdot\nabla\delta\mathbf{v}dx=\int_{\partial\Omega_{A}}\mu_{A}(\mathbf{n}\cdot\nabla)\mathbf{v}\cdot\delta\mathbf{v}dx-\int_{\Omega_{A}}\mu_{A}\Delta\mathbf{v}\cdot\delta\mathbf{v},$
	$\displaystyle-\int_{\Omega_{A}}p\nabla\cdot\delta\mathbf{v}dx=-\int_{\partial\Omega_{A}}p\mathbf{n}\cdot\delta\mathbf{v}dx+\int_{\Omega_{A}}\nabla p\cdot\delta\mathbf{v}dx.$

Similar calculations can be done in $\Omega_{B}$. Combing all this calculations, we immediately derive the corresponding Euler-Lagrange equation for (2.7),

$$\left\{\begin{array}[]{ll}-\mu(x)\Delta\mathbf{v}+\nabla p=0&\\ \nabla\cdot\mathbf{v}=0&\end{array}\right.\qquad\qquad\hbox{in }\Omega,$$

(2.8)

coupled with the interface condition on $\Gamma$

$$[\mathbf{v}]=0,\quad\big{[}\mu(\mathbf{n}\cdot\nabla)\mathbf{v}-p\mathbf{n}\big{]}=\gamma\kappa\mathbf{n}\qquad\qquad\hbox{on }\Gamma,$$

(2.9)

the Navier slip boundary condition on $\Gamma_{S}$

$$v_{n}=0,\quad\beta(x)v_{\tau}=-\mu\frac{\partial v_{\tau}}{\partial n}\qquad\qquad\hbox{on }\Gamma_{S}$$

(2.10)

and the condition for the moving contact line

$$\xi v_{ct}=-\gamma(\cos\theta_{d}-\cos\theta_{Y}).$$

(2.11)

Here

$$\mu(x)=\left\{\begin{array}[]{ll}\mu_{A},&\hbox{if }x\in\Omega_{A};\\ \mu_{B},&\hbox{if }x\in\Omega_{B};\\ \end{array}\right.\quad\hbox{and}\quad\beta(x)=\left\{\begin{array}[]{ll}\beta_{A},&\hbox{if }x\in\Gamma_{SA};\\ \beta_{B},&\hbox{if }x\in\Gamma_{SB}.\\ \end{array}\right.$$

On the outer boundary $\Gamma_{O}$, we also have a relation for the external force that

$$\mathbf{F}_{ext}=\mu(\mathbf{n}\cdot\nabla)\mathbf{v}-p\mathbf{n}.$$

The boundary condition (2.11) for moving contact lines are exactly the model derived in Ren & E (2007). Specifically, when $\xi=0$, the condition is reduced to

$$\theta_{d}=\theta_{Y}.$$

(2.12)

This implies that the microscopic dynamic contact angle is equal to the Young’s angle. Together with the Navier-Slip boundary condition (2.10), this is the model used in the asymptotic analysis in Cox (1986). Therefore, the above analysis indicates that some well-known continuum models for moving contact lines can be derived in a variational framework of the Onsager principle. In the following, we will use the variational principle to derive a reduced model for the apparent dynamic contact angle.

It is known that the two-phase flow has multiscale property near moving contact lines. In general the macroscopic contact angle is different from the microscopic one; see for example Cox (1986), where a dynamic boundary condition for the apparent contact angle was derived by matched asymptotic analysis. In the following, we derive a similar boundary condition by using the Onsager principle as an approximation tool. The derivation is much simpler but still captures the essential dynamics of the apparent contact angle.

As shown in Figure 2, we separate the domain near the contact line in three different scales. The macroscopic region $\Omega$ has a characteristic length $L$. The microscopic region is of molecular scale with a characteristic length $l$. In the microscopic scale, the interaction between the liquid molecules and the solid molecules induces a friction of the contact line (Guo et al., 2013). The mesoscopic region has a characteristic length $D\ll L$, but still much larger than the molecular scale $l$. In this region, the no-slip boundary condition is still a good approximation. The contact angle $\theta_{a}$ represents the apparent angle in macroscopic point of view.

We are interested in the contact line dynamics in the mesoscopic region. For this purpose, we analyze the system by using the Onsager principle as an approximation tool. We make the ansatz that the interface between the two fluids in this region can be approximated well by a straight line $\tilde{\Gamma}$ which has a tilting angle equal to the macroscopic contact angle $\theta_{a}$. With this assumption, we can use the Onsager principle to derive a relation between the apparent contact angle $\theta_{a}$ and the contact line motion as follows.

We first calculate the rate of change of the total free energy in the system. Similar to (2.1), the total approximate free energy $\tilde{\mathcal{E}}$ is composed of three surface energies. The changing rate of the total interface energies with respect to the motion of the contact line is given by

$$\dot{\tilde{\mathcal{E}}}=\gamma(\cos\theta_{a}-\cos\theta_{Y})v_{ct}+\int_{\tilde{\Gamma}}\gamma\kappa v_{n}ds=\gamma(\cos\theta_{a}-\cos\theta_{Y})v_{ct}.$$

(2.13)

Here the second term disappears since the curvature is zero for a straight line. To use the Onsager principle for the open system, we need consider the work to the exterior region on the outer boundary, $\dot{\tilde{\mathcal{E}}}^{*}=-\int_{\tilde{\Gamma}_{out}}\mathbf{F}_{ext}\cdot\mathbf{v}$. This is a higher order term in comparison with $\dot{\tilde{\mathcal{E}}}$, since it is of order $|\mathbf{F}_{ext}|v_{ct}D$ with $D\ll L=O(1)$. We will ignore this term in the following calculations.

We now compute the energy dissipations in the wedge region as shown in Figure 2 (right plot). The total energy dissipations is calculated approximately as

$$\tilde{\Psi}={\xi}v_{ct}^{2}+\int_{\tilde{\Omega}_{A}}\mu_{A}|\nabla\mathbf{v}|^{2}dx+\int_{\tilde{\Omega}_{B}}\mu_{B}|\nabla\mathbf{v}|^{2}dx,$$

(2.14)

In the above formula, $\tilde{\Omega}_{A}$ and $\tilde{\Omega}_{B}$ are the domains corresponding to liquids $A$ and $B$ in the mesoscopic region. The contact line friction term ${\xi}v_{ct}^{2}$ is due to the dissipations in the microscopic region. The velocity $\mathbf{v}$ can be obtained by solving the Stokes equations in wedge regions assuming the interface moving in a steady state. By careful calculations (see Appendix A), we can compute the energy dissipation function approximately

$$\tilde{\Phi}=\frac{1}{2}\tilde{\Psi}\approx\frac{1}{2}{\xi}v_{ct}^{2}+\frac{\mu_{A}|\ln\zeta|}{2\mathcal{F}(\theta_{a},\lambda)}v_{ct}^{2}\approx\frac{1}{2}\frac{\mu_{A}|\ln\zeta|}{\mathcal{F}(\theta_{a},\lambda)}v_{ct}^{2},$$

(2.15)

where the dimensionless parameter $\zeta=D/l$ is the ratio between the mesoscopic size and the microscopic size, and

		$\displaystyle\mathcal{F}(\theta_{a},\lambda)$
	$\displaystyle=$	$\displaystyle\frac{\lambda(\theta_{a}^{2}-\sin^{2}\theta_{a})(\pi-\theta_{a}+\sin\theta_{a}\cos\theta_{a})+((\pi-\theta_{a})^{2}-\sin^{2}\theta_{a})(\theta_{a}-\sin\theta_{a}\cos\theta_{a})}{2\sin^{2}\theta_{a}\Big{(}\lambda^{2}(\theta_{a}^{2}-\sin^{2}\theta_{a})+2\lambda(\sin^{2}\theta_{a}+\theta_{a}(\pi-\theta_{a}))+((\pi-\theta_{a})^{2}-\sin^{2}\theta_{a})\Big{)}},$		(2.16)

with $\lambda=\frac{\mu_{B}}{\mu_{A}}$. In the last approximation in (2.15), we have assumed that $\xi$ is much smaller than the viscous friction coefficient ${\mu_{A}|\ln\zeta|}/{\mathcal{F}(\theta_{a},\lambda)}$. Actually, the friction coefficient $\xi$ is measured directly in experiments in Guo et al. (2013), which is given by $\xi\approx(0.8\pm 0.2)\mu$. Notice that $\ln\zeta$ is generally chosen as a constant of order $10$ (In de Gennes et al. (2003), it is chosen as $13.6$). Direct computations also show that $1/\mathcal{F}(\theta_{a},\lambda)=O(1)$. Therefore, $\xi$ can be regarded a small perturbation to the viscous friction coefficient and can be ignored.

By using the Onsager principle, we minimize the Rayleighian $\mathcal{R}=\tilde{\Phi}+\dot{\tilde{\mathcal{E}}}$ with resepct to $v_{ct}$. This leads to an equation for $v_{ct}$

$$\frac{\mu_{A}|\ln\zeta|}{\mathcal{F}(\theta_{a},\lambda)}v_{ct}=-\gamma(\cos\theta_{a}-\cos\theta_{Y}).$$

(2.17)

This equation basically means that the local dissipative force (due to the viscous dissipation) is balanced by the effective unbalanced Young force, which is the dominant term in the driving force in the mesoscopic region. It gives a relation between the contact line velocity and the apparent contact angle $\theta_{a}$. The relation (2.17) can be rewritten in a dimensionless form

$$|\ln\zeta|Ca=-\mathcal{F}(\theta_{a},\lambda)(\cos\theta_{a}-\cos\theta_{Y}),$$

(2.18)

where $Ca={\mu_{A}v_{ct}}/{\gamma}$ is the capillary number. This is a Cox-type formula for the apparent contact angle and the capillary number of the contact line.

This equation (2.18) is consistent with Cox’s formula in leading order when $Ca$ is small. This will be discussed as follows. We recall the Cox’s formula for small $Ca$,

$$|\ln\zeta|Ca=\mathcal{K}(\theta_{a},\lambda)-\mathcal{K}(\theta_{Y},\lambda).$$

(2.19)

where

	$\displaystyle\mathcal{K}(\theta,\lambda)$
	$\displaystyle=\int_{0}^{\theta}\frac{\lambda(\beta^{2}-\sin^{2}\beta)(\pi-\beta+\sin\beta\cos\beta)+((\pi-\beta)^{2}-\sin^{2}\beta)(\beta-\sin\beta\cos\beta)}{2\sin\beta\Big{(}\lambda^{2}(\beta^{2}-\sin^{2}\beta)+2\lambda(\sin^{2}\beta+\beta(\pi-\beta)+(\pi-\beta)^{2}-\sin^{2}\beta)\Big{)}}\,\mathrm{d}\beta,$
	$\displaystyle=\int_{0}^{\theta}\mathcal{F}(\beta,\lambda)\sin\beta\,\mathrm{d}\beta.$

A first order approximation of Cox’s formula leads to

	$\displaystyle|\ln\zeta|Ca=$	$\displaystyle\int_{\theta_{Y}}^{\theta_{a}}\mathcal{F}(\beta,\lambda)\sin\beta\mathrm{d}\beta$
	$\displaystyle\approx$	$\displaystyle-\mathcal{F}(\theta_{a},\lambda)(\cos\theta_{a}-\cos\theta_{Y})+O((\theta_{a}-\theta_{Y})^{2}).$

This implies that Cox’s formula is consistent with our analysis by using the Onsager principle in leading order. The difference between the two formula is illustrated in Figure 3. We can see that their difference is small for small capillary number case (e.g. $Ca\leq 0.01$). When $Ca$ is large, the linear approximation of the interface in the mesoscopic region is not that accurate any more. Then there is a larger difference between the equation (2.18) and Cox’s formula.

The equation (2.18) can be regarded as a coarse-graining boundary condition for the two-phase Navier-Stokes equation in the macroscopic region $\Omega$. It gives a relation between the apparent contact angle and the contact line velocity. The results are similar to that in de Gennes et al. (2003). One can use (2.18) instead of the equation (2.11) to avoid resolving the mesoscopic region by very fine meshes in numerical simulations.

In some problems with small characteristic size, the energy dissipations in the mesoscopic region may dominate. Then one can derive a reduced model for such problems. One typical example is the spreading of a small droplet on hydrophilic surfaces which gives the so-called Tanner’s law (Tanner, 1979). Other examples can be found in de Gennes et al. (2003); Chan et al. (2020); Xu & Wang (2020). In the following, we will introduce two typical examples. Then we show that they can be described by a unified reduced model, which will be used to study the dynamic contact angle hysteresis in the next section.

Example 1. A moving contact line problem in a micro-channel.

In the example, we consider a contact line problem in a microfluidics, which is similar to that considered in Joanny & De Gennes (1984). As shown in Figure 4, suppose that the two walls of a channel are moving in a velocity $u_{wall}$. There is a bar on the left side of the channel to keep the liquid from moving out. Suppose that the height $h_{0}$ of the channel is smaller than the capillary length. Then we could assume that the interface keeps almost circular if the velocity $u_{wall}$ is small. Then the position of the contact point is fully determined by the dynamic contact angle since the volume of the liquid is conserved.

Suppose the volume of the liquid is $V_{0}$. We suppose the left boundary of the liquid domain is $x=0$ and the $x$-coordinate of the contact point is $x_{ct}$. Suppose the apparent contact angle is $\theta_{a}$. Then the volume of liquid is calculated by

	$\displaystyle V_{0}$	$\displaystyle=h_{0}x_{ct}+\frac{h_{0}^{2}}{4}\frac{\theta_{a}-{\pi}/{2}-\sin(\theta_{a}-{\pi}/{2})\cos(\theta_{a}-{\pi}/{2})}{\sin^{2}(\theta_{a}-\pi/2)}$
		$\displaystyle=h_{0}x_{ct}+\frac{h_{0}^{2}}{8}\frac{(2\theta_{a}-{\pi}+\sin(2\theta_{a}))}{\cos^{2}\theta_{a}}.$

This equation gives a relation between ${x}_{ct}$ and the apparent contact angle $\theta_{a}$. We do time derivative for the above equation and notice that $\dot{V}_{0}=0$. Then we have

$$\dot{x}_{ct}=-\frac{h_{0}}{2}\frac{\cos\theta_{a}+(\theta_{a}-\frac{\pi}{2})\sin\theta_{a}}{\cos^{3}\theta_{a}}\dot{\theta}_{a}.$$

(2.20)

On the other hand, since the relative velocity of the contact line with respect to the two walls is $\dot{x}_{ct}-u_{wall}$, the equation (2.18) implies

$$\frac{|\ln\zeta|\mu}{\gamma}(\dot{x}_{ct}-u_{wall})=-\mathcal{F}_{1}(\theta_{a})(\cos\theta_{a}-\cos\theta_{Y}),$$

(2.21)

where

$$\mathcal{F}_{1}(\theta_{a})=\mathcal{F}(\theta_{a},0)=\frac{(\theta_{a}-\sin\theta_{a}\cos\theta_{a})}{2\sin^{2}\theta_{a}}.$$

(2.22)

The two equations (2.20) and (2.21) compose a complete system of ordinary differential equations (ODEs) for the apparent contact angle and the contact point position. Denote

$$\mathcal{G}_{1}(\theta)=-\frac{\cos^{3}\theta_{a}}{\cos\theta_{a}+(\theta_{a}-\frac{\pi}{2})\sin\theta_{a}},$$

then the ODE system is written as

$$\left\{\begin{array}[]{l}\dot{\theta}_{a}=\frac{2\mathcal{G}_{1}(\theta)}{h_{0}}\Big{(}-\frac{\gamma}{|\ln\zeta|\mu}\mathcal{F}_{1}(\theta_{a})(\cos\theta_{a}-\cos\theta_{Y})+u_{wall}\Big{)},\\ \dot{x}_{ct}=-\frac{\gamma}{|\ln\zeta|\mu}\mathcal{F}_{1}(\theta_{a})(\cos\theta_{a}-\cos\theta_{Y})+u_{wall}.\end{array}\right.$$

(2.23)

Example 2. A capillary problem along a moving thin fiber.

Motivated by the recent experiments in Guan et al. (2016a), we consider a capillary problem along a moving fiber. As shown in Figure 5, we suppose a fiber is inserted in a liquid. When the fiber moves up and down, the contact line will recede and advance accordingly. We assume that the radius $r_{0}$ of the fiber is much smaller than the capillary length $l_{c}$. By the Young-Laplace equation, the radial symmetric liquid-vapor interface can be described by the following equation (see de Gennes et al. (2003))

$$x=H(r)=h-r_{0}\cos\theta_{a}\ln\frac{r+\sqrt{r^{2}-r_{0}^{2}\cos^{2}\theta_{a}}}{r_{0}\cos\theta_{a}},\qquad r\geq r_{0}.$$

(2.24)

Here we assume the upper direction is the positive $x$-direction. There are two parameters $h$ and $\theta_{a}$ in this equation. By the definition of the capillary length, we can assume that the interface intersects with the horizontal surface $x=0$ at $r=l_{c}$. Then we have

$$H(l_{c})=0.$$

(2.25)

This gives a relation between $\theta_{a}$ and $h$, which is

$$h=r_{0}\cos\theta_{a}\ln\frac{l_{c}+\sqrt{l_{c}^{2}-r_{0}^{2}\cos^{2}\theta_{a}}}{r_{0}\cos\theta_{a}}.$$

(2.26)

Notice that the $x$-coordinate of the contact line is given by

$$x_{ct}:=H(r_{0})=h-r_{0}\cos\theta_{a}\ln\frac{1+\sin\theta_{a}}{\cos\theta_{a}}=r_{0}\cos\theta_{a}\ln\frac{l_{c}+\sqrt{l_{c}^{2}-r_{0}^{2}\cos^{2}\theta_{a}}}{r_{0}(1+\sin\theta_{a})}.$$

Direct calculation gives

$$\dot{x}_{ct}=r_{0}\mathcal{G}^{-1}_{2}(\theta_{a})\dot{\theta}_{a},$$

(2.27)

where

$$\mathcal{G}^{-1}_{2}(\theta_{a})=r_{0}^{-1}\frac{\partial x_{ct}}{\partial\theta_{a}}\approx-\Big{(}\sin\theta_{a}\ln\frac{2l_{c}}{r_{0}(1+\cos\theta_{a})}+1-\sin\theta_{a}\Big{)}.$$

The equation (2.27) gives a relation between the time derivative of the contact line position and that of the dynamic contact angle.

In this example, since the relative velocity of the contact line with respect to the fiber surface is $\dot{x}_{ct}-u_{wall}$, the equation (2.18) turns to be

$$|\ln\zeta|\frac{\mu}{\gamma}(\dot{x}_{ct}-u_{wall})=-\mathcal{F}_{1}(\theta_{a})(\cos\theta_{a}-\cos\theta_{Y}).$$

(2.28)

The two equations compose a complete system for the capillary rising problem. They can be rewritten as

$$\left\{\begin{array}[]{l}\dot{\theta}_{a}=\frac{\mathcal{G}_{2}(\theta_{a})}{r_{0}}\Big{(}-\frac{\gamma}{|\ln\zeta|\mu}\mathcal{F}_{1}(\theta_{a})(\cos\theta_{a}-\cos\theta_{Y})+u_{wall}\Big{)},\\ \dot{x}_{ct}=-\frac{\gamma}{|\ln\zeta|\mu}\mathcal{F}_{1}(\theta_{a})(\cos\theta_{a}-\cos\theta_{Y})+u_{wall}.\end{array}\right.$$

(2.29)

This is the formula given in Xu & Wang (2020).

We can see that the ordinary differential systems (2.23) and (2.29) in the two examples have the same structure. They can be written as a unified form as follows,

$$\left\{\begin{array}[]{l}\dot{\theta}_{a}=\frac{\mathcal{G}(\theta_{a})}{l_{0}}\Big{(}-\frac{\gamma}{|\ln\zeta|\mu}\mathcal{F}_{1}(\theta_{a})(\cos\theta_{a}-\cos\theta_{Y})+u_{wall}\Big{)},\\ \dot{x}_{ct}=-\frac{\gamma}{|\ln\zeta|\mu}\mathcal{F}_{1}(\theta_{a})(\cos\theta_{a}-\cos\theta_{Y})+u_{wall}.\end{array}\right.$$

(2.30)

The second equation of (2.30) is due to the condition (2.18) where $\mathcal{F}_{1}(\theta_{a})$ is given in (2.22). Since $\mathcal{F}_{1}(\theta_{a})=\mathcal{F}(\theta_{a},0)$, it corresponds to the case when $\frac{\mu_{B}}{\mu_{A}}=0$ (see in (2.16)), i.e. a liquid-vapor system. In the first equation of (2.30), $l_{0}$ is a characteristic length and the function $\mathcal{G}(\theta_{a})$ comes from the geometric setup of a problem. Both $l_{0}$ and $\mathcal{G}(\theta_{a})$ may change for different problems. Hereinafter, we will use the equation (2.30) as a model problem to study the contact angle hysteresis for a liquid-vapor system on chemically patterned surface.

We first introduce some properties satisfied by the dynamic factor $f$ and $g$:

$$f(\theta)\geq 0,f(0)=0;\quad-M\leq g(\theta)\leq-m,$$

for some positive numbers $m$ and $M$. In addition, we also assume that $f$ is monotonically increasing in the interval $[0,\pi]$. These conditions are quite general and are satisfied by the two examples in the previous section.

We assume that the substrate has periodic chemical patterns with period $\varepsilon$ so that Young’s angle at the position $x$ satisfies $\hat{\theta}_{Y}(x)=\varphi(\frac{x}{\varepsilon})$, where $\varphi$ is a periodic continuously differentiable function with period 1. The minimum and maximum of $\varphi$ are $\theta_{A}$ and $\theta_{B}$ respectively with $0<\theta_{A}<\theta_{B}<\pi$. One example of such functions is

$$\varphi(z)=\frac{\theta_{A}+\theta_{B}}{2}+\frac{\theta_{B}-\theta_{A}}{2}\sin(2\pi z).$$

(3.3)

In the following, we will use the method of averaging to derive the effective dynamics of the contact angle and contact line position. This method has been studied in details (e.g. in E (2011); Pavliotis & Stuart (2008)) and has been widely used in obtaining the effective boundary conditions (e.g., Miskis & Davis (1994))

We introduce the fast spatial variable $y=\frac{\hat{x}_{ct}}{\varepsilon}$ and fast temporal variable $\tau=\frac{t}{\varepsilon}$. Consider the multiple scale asymptotic expansions

$$\theta_{a}=\theta_{0}(t,\tau)+\varepsilon\theta_{1}(t,\tau)+\cdots,\qquad y=y_{0}(t,\tau)+\varepsilon y_{1}(t,\tau)+\cdots.$$

(3.4)

Then the system of equations (3.2) can be rewritten as

		$\displaystyle\frac{1}{\varepsilon}\frac{\partial\theta_{0}}{\partial\tau}+(\frac{\partial\theta_{0}}{\partial t}+\frac{\partial\theta_{1}}{\partial\tau})+\varepsilon(\frac{\partial\theta_{1}}{\partial t}+\frac{\partial\theta_{2}}{\partial\tau})+\cdots$
	$\displaystyle=$	$\displaystyle g(\theta_{0}+\varepsilon\theta_{1}+\cdots)\Big{(}f(\theta_{0}+\varepsilon\theta_{1}+\cdots)(\cos\varphi(y_{0}+\varepsilon y_{1}+\cdots)-\cos(\theta_{0}+\varepsilon\theta_{1}+\cdots))+v\Big{)},$
		$\displaystyle\frac{1}{\varepsilon}\frac{\partial y_{0}}{\partial\tau}+(\frac{\partial y_{0}}{\partial t}+\frac{\partial y_{1}}{\partial\tau})+\varepsilon(\frac{\partial y_{1}}{\partial t}+\frac{\partial y_{2}}{\partial\tau})+\cdots$
	$\displaystyle=$	$\displaystyle\frac{1}{\varepsilon}f(\theta_{0}+\varepsilon\theta_{1}+\cdots)(\cos\varphi(y_{0}+\varepsilon y_{1}+\cdots)-\cos(\theta_{0}+\varepsilon\theta_{1}+\cdots)).$

First order equations in $O(\frac{1}{\varepsilon})$. We have the two equations in the fast time scale:

	$\displaystyle\frac{\partial\theta_{0}}{\partial\tau}=0,$		(3.6)
	$\displaystyle\frac{\partial y_{0}}{\partial\tau}=f(\theta_{0})(\cos\varphi(y_{0})-\cos\theta_{0}).$		(3.7)

The first equation implies that $\theta_{0}=\theta_{0}(t)$ is indeed a slow variable that does not depend on the fast time scale $\tau$. Then for given $\theta_{0}$, the second equation is a simple ordinary differential equation for $y_{0}$ with respect to $\tau$, that can be easily solved. We discuss its solution in two different cases for the parameter $\theta_{0}$.

In the first case when $\theta_{0}\in[\theta_{A},\theta_{B}]$, for any given initial value for $y_{0}$, the solution of (3.7) approaches to some equilibrium value $y_{0,\infty}$, which satisfies $\varphi(y_{0,\infty})=\theta_{0}$. By the phase plane analysis (see Figure 6), we know that there are three types of equilibrium points: $y_{0,\infty}$ is asymptotically stable, semi-stable, or unstable if $\varphi^{\prime}(y_{0,\infty})>0$, $\varphi^{\prime}(y_{0,\infty})=0$ or $\varphi^{\prime}(y_{0,\infty})<0$ respectively. Every solution path must be attracted to either an asymptotically stable point or a semi-stable point, regardless of the initial value of $y_{0}$. For instance, if the initial value $y_{0}^{init}$ satisfies $\varphi(y_{0}^{init})<\theta_{0}$, then $y_{0}$ increases monotonically with $\tau$ until it reaches the nearest equilibrium point $y_{0,\infty}$ which is larger than $y_{0}^{init}$. This equilibrium point must satisfy $\varphi^{\prime}(y_{0,\infty})\geqslant 0$ and thus becomes asymptotically stable or semi-stable. If the initial value $y_{0}^{init}$ satisfies $\varphi(y_{0}^{init})>\theta_{0}$, then $y_{0}$ decreases monotonically with $\tau$ until it reaches the nearest equilibrium point $y_{0,\infty}$ that is smaller than $y_{0}^{init}$. This equilibrium point must also satisfy $\varphi^{\prime}(y_{0,\infty})\geqslant 0$ and thus becomes asymptotically stable or semi-stable.

In the second case that $\theta_{0}>\theta_{B}$ or $\theta_{0}<\theta_{A}$, $\cos\varphi(y_{0})-\cos\theta_{0}$ is either always positive or always negative. As a result, $y_{0}$ is monotonic in $\tau$ and diverges to $\pm\infty$ either increasingly or decreasingly. Moreover, using the method of separable variables, $y_{0}$ can be solved and represented implicitly as

$$Q(y_{0}(t,\tau),\theta_{0})=Q(y_{0}(t,0),\theta_{0})+\tau,$$

(3.8)

where $Q(z,\phi)=\int\frac{\mathrm{d}z}{f(\phi)(\cos\varphi(z)-\cos\phi)}$ is monotonically increasing or decreasing in $z$ and thus can be inverted to give $y_{0}(t,\tau)$.

Second order equations in $O(1)$. We consider the $O(1)$ order equation for $\theta$, which is given by

$$\frac{\partial\theta_{0}}{\partial t}+\frac{\partial\theta_{1}}{\partial\tau}=g(\theta_{0})\Big{(}f(\theta_{0})(\cos\varphi(y_{0})-\cos\theta_{0})+v\Big{)}.$$

(3.9)

We can assume that $\theta_{1}$ is of sublinear growth in $\tau$, i.e., there exists a constant $\alpha\in[0,1)$ such that $|\theta_{1}(t,\tau)|\leqslant C|\tau|^{\alpha}$ for some constant $C$ (see also E (2011)). Define the time averaging operator $<\cdot>_{\tau}$ as

$$<F>_{\tau}=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}F(\tau)\mathrm{d}\tau.$$

Then applying this time averaging operator to equation (3.9) gives rise to

$$\frac{\partial\theta_{0}}{\partial t}=g(\theta_{0})\Big{(}<f(\theta_{0})(\cos\varphi(y_{0})-\cos\theta_{0})>_{\tau}+v\Big{)},$$

(3.10)

where we have used the sublinearity of $\theta_{1}$ to eliminate $<\frac{\partial\theta_{1}}{\partial\tau}>_{\tau}$.

From the previous analysis for the leading order equations, $y_{0}$ is monotonic in $\tau$. Thus we can simplify the averaged term by using (3.7),

$\displaystyle<f(\theta_{0})(\cos\varphi(y_{0})-\cos\theta_{0})>_{\tau}=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}\frac{\partial y_{0}}{\partial\tau}\mathrm{d}\tau=\lim\limits_{T\rightarrow\infty}\frac{y_{0}(t,T)-y_{0}(t,0)}{T}.$

We discuss the limit on the right side for different cases on $\theta_{0}$. In the first case that $\theta_{0}\in[\theta_{A},\theta_{B}]$, this limit is zero since every solution path $y(t,\cdot)$ converges to an asymptotically stable or semi-stable equilibrium point (by the analysis for the leading order equations). In the second case that $\theta_{0}>\theta_{B}$ or $\theta_{0}<\theta_{A}$, we shall show that this limit is $\Big{(}\int_{0}^{1}\frac{\mathrm{d}z}{f(\theta_{0})(\cos\varphi(z)-\cos\theta_{0})}\Big{)}^{-1}$, the harmonic average of $f(\theta_{0})(\cos\varphi(z)-\cos\theta_{0})$ over $z\in[0,1]$.

The argument for the second case is as follows. Let

$$b=\int_{0}^{1}\frac{\mathrm{d}z}{f(\theta_{0})(\cos\varphi(z)-\cos\theta_{0})}=\frac{1}{f(\theta_{0})}\int_{0}^{1}\frac{\mathrm{d}z}{(\cos\varphi(z)-\cos\theta_{0})}.$$

From (3.8) we know that if $\tau=nb$ for any integer $n$, then

$$y_{0}(t,\tau)-y_{0}(t,0)=n.$$

Writing $T=nb+a$ with $n=\lfloor T/b\rfloor$ (the largest integer no greater than $T/b$) and $0\leqslant a<b$, and letting $n\rightarrow\infty$, we have $\lim\limits_{T\rightarrow\infty}\frac{y_{0}(t,T)-y_{0}(t,0)}{T}=b^{-1}$. This implies that $<f(\theta_{0})(\cos\varphi(y_{0})-\cos\theta_{0})>_{\tau}$ is exactly equal to the harmonic average of $f(\theta_{0})(\cos\varphi(z)-\cos\theta_{0})$ over a period $[0,1]$.