# Moving Contact Line meeting

### 26 November 2021

## Advancements in the characterization of contact line friction

### Contact line friction models

Setup: shear droplet (like in the MD-PF-VoF paper):


#### Wall relaxation (Yue & Feng 2011)
Obtained from the wetting boundary condition of the Cahn-Hilliard-Navier-Stokes model with the following assumptions:
- Sharp interface limit
- Small capillary number
- No-slip

$$
\color{red}{\mu_f^*(\theta)}Ca = \cos{\theta_0} -\cos{\theta} \; , \quad \mu_f^* = \frac{2\sqrt{2}}{3}\hat{\mu}_f^*\sin{\theta} \; ,
$$

being $Ca=\mu U/\gamma$, where $U$ is the _absolute_ contact line speed, $\mu$ is the fluid's viscosity and $\gamma$ is the interface's surface tension; $\hat{\mu}_f^*=\mu_f/\mu$ is the non-dimentional baseline c.l. friction parameter to be estimated.

This model predicts the line friction parameter to be a function of the absolute contact angle; however, the dependence is too weak to capture asymmetries between wetting and de-wetting:


If we fit advancing and receding contact lines separately, we obtain a different baseline contact line friction value:


| $\theta_0$ | Advancing $\hat{\mu}_f^*$ | | Receding $\hat{\mu}_f^*$ |
| ---- | ----- |-| ----- |
| $127^\circ$ | $0.79\pm0.019$ | $\color{red}{\boldsymbol{>}}$ | $0.093\pm0.030$ |
| $95^\circ$ | $3.21\pm0.059$ | $\color{green}{\boldsymbol{\sim}}$ | $3.20\pm0.13$ |
| $69^\circ$ | $7.20\pm0.24$ | $\color{green}{\boldsymbol{\sim}}$ | $7.48\pm0.38$ |
| $38^\circ$ | $7.48\pm1.04$ | $\color{blue}{\boldsymbol{<}}$ | $20.5\pm1.64$ |

It would be tempting to assert that **hydrophilic surfaces are easier to wet rather than de-wet, and vice-versa for hydrophobic surfaces**. However, we know that friction is not dictated by the liquid-surface interaction energy itself, but rather by energy barriers between adsorption sites (or equivalently, by the corrugation of the surface energy landscape).

#### Molecular Kinetic Theory (Blake & Haynes 1969)


Obtained by appying the reaction kinetics formalism to the adsorption-desorption process of a single layer of molecules on a substrate. The full formula contains many parameters that are related to the molecular details of the contact line motion and to the liquid-surface interaction energy. We can define the following angle deviation:

$$
\delta\cos\theta = \cos{\theta_0}-\cos{\theta}
$$

and expanding until order 3 the following simplified formula is obtained:

$$
\color{red}{\mu_f^*(\delta\cos\theta)}Ca = \delta\cos\theta \; , \quad \mu_f^*(\delta\cos\theta) = \frac{\hat{\mu}_f^*}{1+\beta(\delta\cos\theta)^2} \; ,
$$

where $\hat{\mu}_f^*$ and $\beta$ are two shape parameters to be obtaine by fitting experimental or simulation data.

Contact line friction is a **symmetric function of the difference** between the cosine of the equilibrium contact angle and the cosine of the dynamic contact angle. Therefore, this model too is unable to capture asymmetric wetting/de-wetting.


| $\theta_0$ | | $127^\circ$ | $95^\circ$ | $69^\circ$ | $38^\circ$ |
| - | - | - | - | - | - |
| $\hat{\mu}_f^*$ | | $0.433\pm0.14$ | $2.436\pm0.048$ | $5.870\pm0.33$ | $4.782\pm0.42$ |

However the overall trend is similar: **a higher equilibrium contact angle generally corresponds to a lower c.l. friction**.

The baseline contact line friction value corresponds to the inverse slope of the linearized model close to $\theta=\theta_0$, $Ca=0$.

#### Two-layers model (Johansson & Hess 2018)


This model assumes friction to be the product of energy barriers in the motion of molecules in the _first and second_ liquid layers closest to the solid wall. The steeper the contact angle (assume $\theta<\pi/2$), the easier is for liquid molecules at the interface to be transported to the c.l. region; we thus have an energy barrier that scales with the contact angle:

$$
\color{red}{\mu_f^*(\theta)}Ca = \delta\cos\theta \; , \quad \mu_f^*(\theta) = \hat{\mu}_f^*\cdot e^{\color{red}{\Delta E(\theta)}} \; , \quad \Delta E(\theta) = a\Bigg(\frac{1}{2}\sin\theta+\cos\theta\Bigg)^2 \; ,
$$

where now the parameters to estimate are again the baseline friction $\hat{\mu}_f^*$ and the scaling coefficient of the energy barrier $a$.

Not only this model is **strongly nonlinear**, but also **strongly dependent on the absolute contact angle**. Therefore there is hope this model could capture asymmetric wetting/de-wetting automatically.


However when I look at the numbers, they don't seem to make much sense:

| $\theta_0$ | | $127^\circ$ | $95^\circ$ | $69^\circ$ | $38^\circ$ |
| - | - | - | - | - | - |
| $\hat{\mu}_f^*$ | | $0.35\pm1.51$ | $2.70\pm0.05$ | $6.97\pm0.30$ | $0.002\pm0.003$ |
| $a$ | | $1.51\pm2.15$ | $0.19\pm0.05$ | $-0.19\pm0.06$ | $6.90\pm1.28$ |

**To-do list**
- Understand how $\hat{\mu}_f^*$ and $a$ should vary when the wall is more hydrophylic (resp. hydrophobic).
- In other words: can the model be refined to add energy barriers not only for jump/rolling but also for wetting and de-wetting?

### Agreement between forced and spontaneous wetting

Let's consider now the situation of a droplet spreading spontaneously:


We can plot the measurement of $Ca$ and $\theta$ obtained from droplet spreading/retracting simulations (marked $\circ$) and from droplets under shear conditions (marked $\diamond$).


Visually, they agree very well. But does it prove spontaneous and forced wetting can be both captured by the same contact line friction parameters? In order to prove it we follow the procedure below:
- Estimate the contact line friction by fitting one (or more) model to droplet spreading measurements;
- Evaluate the model with the contact angle from shear droplet simulations;
- Predict the capillary number and compare it with the one imposed by the wall motion (assuming no-slip).

$$
\widehat{Ca} = \frac{\cos\theta_0-\cos\theta_{shear}}{\mu_f^*(\theta_{shear})} \; .
$$

Evaluating the line friction model can be understood by looking at the 'rollercoaster plots' below. The equilibrium contact angle is $\theta_0\simeq69^\circ$.


The imposed capillary number (marked $\square$) is plotted alongside the predicted ones using the two-layer model and MKT (marked $\circ$ and $\diamond$).


The error bars corresponds to the uncertainty on the dynamic contact angle ($\delta\theta\mapsto\delta\mu$), **not** the modelling uncertainty, which is still to be estimated. The imposed $Ca$ is almost always within the error bars of both models (expected, since there should be no asymmetry for $q_3$)

**To-do list**
- Repeat the same analysis for all other substrate silica charge values $(q_1,q_2,q_4)$.
- Estimate the error in the parameter's models (inclided $\theta_0$).

### Characterization of friction through the classification of molecular motion

Can we connect the motion of single molecules to the displacement of the contact line?

How is the three-phase contact-line even defined?


We use a **clustering** procedure ([DBSCAN](https://en.wikipedia.org/wiki/DBSCAN)) to distinguish between the molecules in the liquid phase and the one absorbed on the substrate ahead (Gauss-Pomeau film?); after that, we construct a **concave hull** ([alpha shape](https://en.wikipedia.org/wiki/Alpha_shape)), which allows to define the contact line position and displacement as a function of the system's width at a given time frame:

$$
x_{cl} = x_{cl}(y,t) \; \implies \; \Delta x_{cl}(y,t) = x_{cl}(y,t) - x_{cl}(y,t-\tau) \; .
$$

[Video of clustering + alpha shape](https://drive.google.com/file/d/1r1kFj82SXIDlzXCYY3u-GFqQ83bB-G3V/view?usp=sharing)

[Video of contact line displacement](https://drive.google.com/file/d/1eNsh1QIqys2JE2EILal3PfghOP52PBGr/view?usp=sharing)

The determination of the contact line motion is hence sensitive to three arbitrary parameters:
- The network distance used by DBSCAN, $d$;
- The 'sharpness' of the alpha shape, $\alpha$;
- The time lag between two consecutive positions, $\tau$.

We map the motion of molecules to a displacement to compare to the contact line motion:

$$
(x_i(t),y_i(t)) \; \mapsto \; x_i(y_i,t) \; \implies \; \Delta x_i(y_i^*,t) = x_i(y_i,t) - x_i(y_i,t-\tau) \; , \quad y_i^* =y_i(t) \; ,
$$

and compare it with the contact line displacement, finding the molecules whose motion correlates more with the contact line one.

A molecule is considered a direct contributor to contact line motion (_true positive_) if:
- The molecule is not too far from the contact line;
- The molecule has moved far enough in $\tau$, so that it's motion is probably not just a thermal fluctuation;
- The molecule hasn't moved across the periodic boundary (difficult to deal with);
- The molecule displacement has the same sign as the c.l. one (sanity check).
If one of the above is not true, the motion is discarded as a _false positive_.

The most arbitrary parameters in this procedure are $\{d,\alpha,\tau\}$, so we perform a **numerical scanning** of the three to find the combination that minimizes the fraction of false positives or (possibly and) maximizes the number of true positives.


The results shown are for an equilibrium contact line having contact angle $\theta_0=69^\circ$. The clustering distance $d$ does not seem to influence the results too much, so we take a value close to the **characteristic distance of h-bonds** ($d=0.375$nm); on the other hand, **the sharper the contact line, the more true positives are found** (makes sense) and $\tau=1$ps optimizes both the fraction of false positives and the number of true positives.

For $\theta_0=69^\circ$, at equilibrium (no moving wall: c.l. fluctuates but do not move on average) over 3500 frames (3.5ns):
- events motion forth = 2949
- events motion back = 3054

Doesn't seem too bad.

**To-do list**
- Add the wall velocity **bias** $\delta x_b = U_w\tau$ to the classification scheme and apply it to a **moving wall** situation.
- Classify molecular motion: rolls/jumps (Johannson, Hess 2018), motion involving more than one molecule at the time?

## New research directions

### Streamlines

Since we have quite a lot of MD data, let's take a long-time average of the momentum flow at the contact line and try to obtain streamlines.

The following figures present a zoomed window of $\sim$ 12.37nm x 6.55nm; streamlines are obtained via Paraview's stream tracer instrument.

$\theta_0=127^\circ$, $Ca=0.6$, receding c.l.

$\theta_0=95^\circ$, $Ca=0.2$, receding c.l.

$\theta_0=69^\circ$, $Ca=0.1$, receding c.l.

Streamlines are the level sets of the stream function $\psi$, which in 2D can be obtained from the following elliptic equation:

$$
\nabla^2\psi = -\omega = -(\nabla\times\boldsymbol{u})\cdot\boldsymbol{\hat{k}} \; .
$$

Unfortunately $\boldsymbol{u}=(u,v)$ is too irregular and we cannot compute its curl. However, we can cast the stream function equation into its integral (weak) formulation; considering a sufficiently regular test function $\varphi\in\mathcal{H}^1_0(\Omega)$ (right?):

$$
\int_{\Omega} \varphi\nabla^2\psi = -\int_{\Omega}\nabla\times\boldsymbol{u}\cdot\boldsymbol{\hat{k}}\varphi \; ,
$$

$$
\int_{\partial\Omega}\varphi\nabla\psi\cdot\boldsymbol{\hat{n}} - \int_{\Omega}\nabla\varphi\cdot\nabla\psi = \int_{\Omega}u\frac{\partial\varphi}{\partial z} - \int_{\Omega}v\frac{\partial\varphi}{\partial x} - \oint_{\partial\Omega}(\boldsymbol{u}\times\boldsymbol{\hat{n}})\cdot\boldsymbol{\hat{k}}\varphi \; ,
$$

discard contour terms (choice of Sobolev's space):

$$
-\int_{\Omega}\nabla\varphi\cdot\nabla\psi = \int_{\Omega}u\frac{\partial\varphi}{\partial z} - \int_{\Omega}v\frac{\partial\varphi}{\partial x}
$$

and solve it with FreeFem++ taking the $(u,v)$ from MD (and with a Diriclet condition on some part of the boundary to ensure the problem has only one solution).

**To-do list**
- Obtain the streamfunction by solving the Laplace equation (solving the equation in its weak formulation).
- Compare to theoretical models:
- Huh & Scriven 1971
- Fricke et al. 2019
- Shikhmurzaev (?)

### System of two immiscible fluids

Biphase system of water and cyclohexane (completely immiscible), taken from the [GROMACS tutorial](http://www.mdtutorials.com/gmx/biphasic/index.html); I added the silica walls (structurally the same as the previous projects, different force field).


**Advantages**
Velocity is well-defined in the non-water phase also; hence:
- Better streamlines;
- Possibility of testing GNBC boundary condition from MD (Quian et al 2003): $\beta v_x^{slip}+[\eta\partial_n v_x](0)+\widetilde{\sigma}^Y_{nx}(0)=0$, difficult if velocity in the other phase is not defined:
- Maybe the equivalent PF won't have streamlines crossing and thus we can focus on the localized slip at the contact line?

**To-do**
- Implement pressure coupling with position restraints;
- Re-calibrate the contact angle with the new force field and the new fluid.

### Droplet's equilibrium, mobility and evaporation over random surfaces

Equilibrium properties of droples over non-homogeneous surfaces are still not well understood ([How Wenzel and Cassie Were Wrong](https://pubs.acs.org/doi/full/10.1021/la062634a)).

The main physical aspects to model in random surfaces are:
- The amount of defects (ratio between chemical species);
- The size of defects (in proportion to the size of the droplet);
- The proxmity of defects to the three-phase contact line.

The first two aspects can be captured by a Ising-like model, where hydrophobic/hydrophilic sites correspond to the orientiation of what would be magnetic spins.

$$
H = -\sum_{i,j} J_{ij}\sigma_i\sigma_j - \mu\sum_j h_j\sigma_j
$$

Example:

We can produce a substrate of high-$q$ / low-$q$ quadrupoles mimicking the output of a Ising model MC simulation.
Blue quadrupoles are hydrophilic and yellow ones are hydrophobic.


**Questions**
- What is the equilibrium contact angle? Is it really a function of the proportion between high-$q$ and low-$q$ alone? Or is more related to the contact line metastability?
- How does defects' size and density influence c.l. shape and near-wall interface curvature?
- How is droplet mobility related to defects' size and density?
- How defects' size and density influence heat tranfer and evaporation?