A Two-Phase SPH Model for Dynamic Contact Angles Including Fluid–Solid Interactions at the Contact Line

Abstract

The description of wetting phenomena on the continuum scale is a challenging problem, since intermolecular interactions, like van der Waals forces between liquid and solid, alter the flow field at the contact line. Recently, these effects were included in the smoothed particle hydrodynamics method by introducing a contact line force (CLF) on the continuum scale. This physically based contact line force model is employed here to simulate two-phase flow in a wide range of wetting dynamics parametrized by capillary number. In particular, dynamic contact angles at various capillary number values are calculated by CLF model and compared to measured values. We find that there is significant disagreement between simulated and measured results, specially at low wetting speeds. It is indeed expected that most of the driving force is dissipated to overcome strong liquid–solid interactions, which are not adequately accounted for in the existing CLF model. Therefore, we have extended that model to account for stick-slip (SSL) behavior of the contact line caused by solid–fluid interactions. The new SSL model results in dynamic contact angle values that are in good agreement with experimental data for the full range of wetting dynamics.

Appendices

Appendix A

Dynamic contact angles $\theta _\mathrm{D}$ are usually modeled as a as a feature of viscous dissipation (de Gennes 1985). The term $\mu \Delta \mathbf {v}$ in the Navier–Stokes equation (Eq. 1) is the source of the viscous energy dissipation $\varPhi$. Thus, the effect of stick-slip phenomena in our model can be expressed as follows. At a high capillary number ($Ca > 0.1$ in the example given in Fig. 5), there is almost no viscosity increase and it becomes almost constant ($\mu _\mathrm{SSL} \approx \mu$). Since $\Delta \mathbf {v}$ monotonically increases as the contact line velocity increases $\left( \frac{\hbox {d}\Delta \mathbf {v}}{\hbox {d}v_\mathrm{cl}} > 0 \right)$, the viscous dissipation is also increasing $\left( \frac{\hbox {d}\varPhi }{\hbox {d}v_\mathrm{cl}} > 0 \right)$, because $\mu$ is almost constant. We refer to this behavior as steady contact line movement in Fig. 15. Now, with strongly decreasing viscosity and with increasing dynamics at slow wetting processes $\left( \frac{\hbox {d}\mu _\mathrm{SSL}}{\hbox {d}v_\mathrm{cl}} < 0 \right)$, it can happen that the viscous dissipation decreases, while the contact line accelerates $\frac{\hbox {d}\varPhi }{\hbox {d}v_\mathrm{cl}} < 0$. This property of our model is in accordance with results from molecular-dynamics simulation of Thompson and Robbins (1990), where the frictional force is reduced in the transition from the ordered state, when the fluid molecules are sticking to the solid, to the unordered sliding state. Such a flow is unstable and will adjust itself toward one of the two possible stable states with $\frac{\hbox {d}\varPhi }{\hbox {d}v_\mathrm{cl}} > 0$. This stick-slip behavior is qualitatively shown in Fig. 15.

Fig. 15

Schematic plots of the viscous dissipation in the vicinity of a contact line in the transition between stick-slip wetting regime and a steady contact line movement. The dashed lines 1 and 2 represent the viscous dissipation with a constant high viscosity $\mu _\mathrm{SSL} = \mu (1+\alpha )$ (for $Ca \rightarrow 0$) and constant low viscosity $\mu _\mathrm{SSL} = \mu$, respectively. The solid line is a qualitative plot when the new stick-slip model is applied

Appendix B

An appropriate test case to analyze the accuracy of the applied SPH scheme is the two-dimensional Taylor–Green flow. This test case is well discussed in the literature for the SPH method (Oger et al. 2016). We consider a Reynolds number $Re = \frac{l\,u_\mathrm{max}}{\nu } = 100$ with characteristic length $l = 1m$, the maximum initial velocity $u_\mathrm{max} = 1\frac{m}{s}$, periodic boundary conditions and a resolution $l/\Delta x$ between 25 and 100. The initial divergence-free velocity field is given as:

$$\begin{aligned} \begin{aligned} v_x^*&= \sin (2\pi x^*)\cos (2\pi y^*)\\ v_y^*&= -\cos (2\pi x^*)\sin (2\pi y^*) \end{aligned} \end{aligned}$$

(47)

The dimensionless velocities are $v_x^*= \frac{v_x}{v_\mathrm{max}}$ and $v_y^*= \frac{v_y}{v_\mathrm{max}}$, respectively. The dimensionless coordinates are $x^* = \frac{x}{l}$ and $y^* = \frac{y}{l}$, respectively. The analytic solution of the dimensionless kinetic energy of the system over time is given as:

$$\begin{aligned} E_\mathrm{kin}^* = \frac{E_\mathrm{kin}}{E_\mathrm{kin}^0} = e^{\left( -16\,\pi ^2\,\nu \,t\right) } \end{aligned}$$

(48)

with the initial kinetic energy $E_\mathrm{kin}^0$ and the kinematic viscosity $\nu$. The analytic solution for the local pressure field is given as:

$$\begin{aligned} p^*(x,y) = \frac{1}{2} \left[ \cos (4\pi x^*) + \cos (4\,\pi y^*)\right] e^{(-16 \pi ^2/Re)t^*}. \end{aligned}$$

(49)

In Fig. 16, we exemplary show the velocity and pressure distribution for $Re = 100$ at $t^* = 0.5$ with a resolution of $L/\Delta x = 100$. In Fig. 17a, the convergence for the kinetic energy decay predicted by the proposed ISPH scheme is shown and compared against the analytic solution for $Re = 100$ for resolutions $L/\Delta x = 25$, $L/\Delta x = 50$ and $L/\Delta x = 100$. In Fig. 17b, the comparison between simulated dimensionless pressures $p^*$ and analytic solution is shown on the line between the points $P_0 = (0|0)$ and $P_1 = (1|1)$. The line is also plotted in Fig. 16b.

Fig. 16

Taylor–Green flow: Particle snapshots for $t^*=0.5$, a dimensionless velocity field $v^*$, b dimensionless pressure field $p^*$

Fig. 17

Taylor–Green flow: convergence of the kinetic energy decay at $Re=100$ obtained with proposed ISPH scheme