An atomistic–continuum hybrid scheme for numerical simulation of droplet spreading on a solid surface

Abstract

We present an atomistic–continuum hybrid method to investigate spreading dynamics of drops on solid surfaces. The Navier–Stokes equations are solved by the finite-volume method in a continuum domain comprised of the main body of the drop, and atomistic molecular dynamics simulations are used in a particle domain in the vicinity of the contact line. The spatial coupling between the continuum and particle domains is achieved through constrained dynamics of flux continuities in an overlap domain.

Appendix

The computational domain is fixed as a unit square by using the transformation

$$\xi=\frac{x}{X},\quad \eta=\frac{y}{h} .$$

(16)

The transformation metrics for this coordinate transformation are given by

$$\xi_{,x}=\frac{1}{X}; \quad \eta_{,x}=-\frac{h'}{h}\eta$$

(17)

$$\eta_{,y}=\frac{1}{h} ; \quad \eta_{,xx}=-\frac{h''}{h}\eta+\frac{2{h'}^2}{h^2}\eta$$

(18)

$$h'=\frac{1}{X}h_{,\xi}; \quad h''=\frac{1}{X^2}h_{,\xi\xi}$$

(19)

where the subscript,“x” denotes differentiation with respect to x. The Navier–Stokes equations are transformed to the following equations:

$$\begin{aligned} &u_{,t}+u\left(\xi_{,x} u_{,\xi}+\eta_{,x}u_{,\eta}\right)+v\eta_{,y}u_{,\eta} = -\left(\xi_{,x} p_{,\xi}+\eta_{,x} p_{,\eta}\right)\\ &+ \frac{1}{Re}\left[\xi_{,x}^2 u_{,\xi \xi}+2\xi_{,x}\eta_{,x} u_{,\xi \eta}+\eta_{,xx} u_{,\eta}+\left(\eta_{,x}^2+\eta_{,y}^2\right) u_{,\eta \eta}\right] \end{aligned}$$

$$\begin{aligned} &v_{,t}+u\left(\xi_{,x} v_{,\xi}+\eta_{,x} v_{,\eta}\right)+v\eta_{,y} v_{,\eta}= -\eta_{,y} p_{,\eta}\\ &+\frac{1}{Re}\left[\xi_{,x}^2 v_{,\xi\xi}+2\xi_{,x}\eta_{,x} v_{,\xi \eta}+\eta_{,xx} v_{,\eta} +\left(\eta_{,x}^2+\eta_{,y}^2\right) v_{,\eta\eta}\right] \end{aligned}$$

Similarly, the boundary conditions are transformed as follows:

$$u=v=0;\quad p_{,\eta}=0\quad at \quad \eta=0$$

(20)

$$u=0; \quad v_{,\xi}=0;\quad p_{,\xi}=0 \quad at \quad \xi=0$$

(21)

$$\begin{aligned} &Re p+\frac{2h'}{1+{h'}^{2}}\left(\eta_{,y} u_{,\eta}+\xi_{,x} v_{,\xi}+\eta_{,x} v_{,\eta}\right)-\eta_{,y} v_{,\eta}\\ &-\left(\xi_{,x} u_{,\xi}+\eta_{,x} u_{,\eta}\right){h'}^{2} +h{''}\left(1+{h}^{'2}\right)^{-\frac{3}{2}}=0 \quad at \quad \eta=1 \end{aligned}$$

(22)

$$\begin{aligned} &\left(\eta_{,y} u_{,\eta}+\xi_{,x} v_{,\xi}+\eta_{,x} v_{,\eta}\right)\left(1-{h'}^2\right)\\ &-2\left(\xi_{,x} u_{,\xi}+\eta_{,x} u_{,\eta}-\eta_{,y} v_{,\eta}\right)h{'}=0 \quad at \quad \eta=1 \end{aligned}$$

(23)

$$\frac{\partial h}{\partial t} +uh'-v=0 \quad at\quad \eta=1$$

(24)

In these equations, the Reynolds number Re ≡ ρ UR/μ is defined based on the characteristic velocity U≡ γ/μ.

To solve these equations, the computational domain is divided into cells of size $\varDelta\xi=1/I$ by $\varDelta\eta=1/K$, with the coordinates of the cell centers given by $\xi_i=(i-1/2)\varDelta\xi$ (i = 1, ..., I), and $\eta_k=(k-1/2)\varDelta\eta (k=1,\ldots, K)$. A staggered computational grid is used wherein the pressure nodes are positioned at the center of the cells, the x- velocity nodes at the midpoint of the right faces, and the y- velocity nodes at the midpoint of the top faces. Spatial discretization of the transformed equations using a second-order accurate finite difference approximation results in a nonlinear coupled set of equations which were solved using the following time-splitting algorithm [34, 35]:

$$\frac{{\bf u}^{*}-{\bf u}^{n}}{\varDelta t_c}=-{\bf u}^{n}\cdot\nabla {\bf u}^{n}+\frac{\nabla^{2}{\bf u}^{n}}{Re}$$

(25)

$$\nabla^2 p^{n+1}=\frac{1}{\varDelta t_c} \nabla\cdot {\bf u}^{*}$$

(26)

$${\bf u}^{n+1}={\bf u}^{*}-\varDelta t_c\nabla p^{n+1}$$

(27)

The transformed equations resulting from the splitting algorithm take the following form:

$$\begin{aligned} &u^*=u+\varDelta t_c -u\left(\xi_{,x} u_{,\xi}+\eta_{,x} u_{,\eta}\right)-v \eta_{,y} u_{,\eta}\\ &+\frac{1}{Re}\left[\xi_{,x}^{2} u_{,\xi\xi}+2\xi_{,x}\eta_{,x} u_{,\xi \eta}+\eta_{,xx} u_{,\eta}+\left(\eta_{,x}^2+\eta_{,y}^2\right) u_{,\eta\eta}\right] \end{aligned}$$

(28)

$$\begin{aligned} &v^*=v+\varDelta t_c -u\left(\xi_{,x} v_{,\xi}+\eta_{,x} v_{,\eta}\right)-v \eta_{,y} v_{,\eta}\\ &+\frac{1}{Re}\left[\xi_{,x}^{2} v_{,\xi\xi}+2\xi_{,x}\eta_{,x} v_{,\xi \eta}+\eta_{,xx} v_{,\eta}+\left(\eta_{,x}^2+\eta_{,y}^2\right) v_{,\eta\eta}\right] \end{aligned}$$

(29)

$$\begin{aligned} &\xi_{,x}^2 p_{,\xi \xi}+2\xi_{,x}\eta_{,x} p_{,\xi\eta}+\eta_{,xx} p_{,\eta}+\left(\eta_{,x}^2+\eta_{,y}^2\right) p_{,\eta \eta}\\ &=\frac{1}{\varDelta t_c}\left(\xi_{,x} u^*_{,\xi}+\eta_{,x} u^*_{,\eta}+\eta_{,y} v^*_{,\eta}\right) \end{aligned}$$

(30)

$$u=u^*-\varDelta t_c\left(\xi_{,x} p_{,\xi}+\eta_{,x} p_{,\eta}\right)$$

(31)

$$v=v^*-\varDelta t_c\left(\eta_{,y} p_{,\eta}\right)$$

(32)

During each time step, the intermediate velocity field u ^* is calculated from Eq. (25), and used to obtain the pressure at the next time step through an iterative solution of Eq. (26). Subsequently, the new pressure distribution is used in Eq. (27) to advance the velocity field to the next time step. Once the new velocity field, u ⁿ⁺¹, is determined, the interface position $h(x,t+\varDelta t_c)$ is updated using the kinematic condition (Eq. (8)).