A Parallel Finite Element Method for 3D Two-Phase Moving Contact Line Problems in Complex Domains

Abstract

Moving contact line problem plays an important role in fluid-fluid interface motion on solid surfaces. The problem can be described by a phase-field model consisting of the coupled Cahn–Hilliard and Navier–Stokes equations with the generalized Navier boundary condition (GNBC). Accurate simulation of the interface and contact line motion requires very fine meshes, and the computation in 3D is even more challenging. Thus, the use of high performance computers and scalable parallel algorithms are indispensable. In this paper, we generalize the GNBC to surfaces with complex geometry and introduce a finite element method on unstructured 3D meshes with a semi-implicit time integration scheme. A highly parallel solution strategy using different solvers for different components of the discretization is presented. More precisely, we apply a restricted additive Schwarz preconditioned GMRES method to solve the systems arising from implicit discretization of the Cahn–Hilliard equation and the velocity equation, and an algebraic multigrid preconditioned CG method to solve the pressure Poisson system. Numerical experiments show that the strategy is efficient and scalable for 3D problems with complex geometry and on a supercomputer with a large number of processors.

Appendix: The Structures of Element Stiffness Matrix and Element Right-Hand Side of the Proposed Discretization Schemes

In this appendix we provide some details for the structures of element stiffness matrix and element right-hand side in Step 1, 3, and 4 of the proposed discretization schemes.

In Step 1, we denote $\{\varphi _i\}^{N_{w_h}}_{i=1}$ as the piecewise linear finite element basis functions for $W_h$ with $N_{w_h}$ as its dimension within an element. The structures of the element stiffness matrix and element right-hand side arising from (2.8)–(2.9) are as follows:

$$\begin{aligned} \begin{pmatrix}K^n_{\phi \phi } &{} K^n_{\phi \mu }\\ K^n_{\mu \phi } &{} K^n_{\mu \mu }\end{pmatrix} ~\text {and}~ \begin{pmatrix}F_\phi ^n\\ F_\mu ^n\end{pmatrix}, \end{aligned}$$

(5.1)

where for $1\le i,j \le N_{w_h}$,

$$\begin{aligned} K^n_{\phi \phi }(i,j)&=\frac{1}{\delta t}({\varphi _i},\varphi _j),\\ K^n_{\phi \mu }(i,j)&={{\mathcal {L}}}_d(\nabla \varphi _i,\nabla \varphi _j),\\ K^n_{\mu \phi }(i,j)&=-\epsilon (\nabla \varphi _i,\nabla \varphi _j)-\frac{s}{\epsilon }(\varphi _i,\varphi _j)-\left( \frac{1}{{{\mathcal {V}}}_s\delta t}+\tilde{\alpha }\right) \langle \varphi _i,\varphi _j \rangle _{\varGamma _w},\\ K^n_{\mu \mu }(i,j)&=(\varphi _i, \varphi _j),\\ F_{\phi }^n(i)&=\frac{1}{\delta t}(\phi _h^n,\varphi _i)-(\mathbf{u}^n_h\cdot \nabla \phi _h^{n},\varphi _i),\\ F_{\mu }^n(i)&=\frac{1}{\epsilon }\left( (\phi _h^n)^3-(1+s)(\phi _h^n), \varphi _i\right) \\&\quad +\,\left\langle \frac{1}{\mathcal {V}_s}\left( {{- \phi _h^n }\over {\delta t}}+u^n_{\tau _1,h}\partial _{\tau _1}\phi _{h}^{n}+u^n_{\tau _2,h}\partial _{\tau _2}\phi _{h}^{n}\right) ,\varphi _i\right\rangle _{\varGamma _w}\\&\quad +\,\left\langle \left( Q\left( \phi ^n_h\right) -\tilde{\alpha }\phi _h^n\right) ,\varphi _i\right\rangle _{\varGamma _w}. \end{aligned}$$

In Step 3, we denote $\{\psi _i\}^{N_{u_h}}_{i=1}$ as the piecewise linear finite element basis functions for $\mathbf U _h$ with $N_{u_h}$ as its dimension for each component within an element. The structures of the element stiffness matrix and element right-hand side arising from (2.13) are as follows:

$$\begin{aligned} \begin{pmatrix}K^n_{u_xu_x} &{} K^n_{u_xu_y}&{} K^n_{u_xu_z} \\ K^n_{u_yu_x} &{} K^n_{u_yu_y} &{} K^n_{u_yu_z} \\ K^n_{u_zu_x} &{} K^n_{u_zu_y} &{} K^n_{u_zu_z} \end{pmatrix} ~\text {and}~ \begin{pmatrix}F_{u_x}^n\\ F_{u_y}^n\\ F_{u_z}^n\end{pmatrix}, \end{aligned}$$

(5.2)

where for $1\le i,j \le N_{u_h}$,

$$\begin{aligned} K^n_{u_xu_x}(i,j)&=\,Re\left( \psi _i,\left( \frac{{\frac{1}{2}}(\rho _h^{n+1}+\rho _h^n)}{\delta t} +\rho _h^{n+1}(\mathbf{u}^n_h\cdot \nabla ) +{\frac{1}{2}}\left( \nabla \cdot (\rho _h^{n+1}\mathbf{u}^n_h)\right) \right) \psi _j\right) \\&\quad +\,(\nabla \psi _i,\eta _h^{n+1}\nabla \psi _j)+(\partial _x\psi _i,\eta _h^{n+1}\partial _x\psi _j)\\&\quad +\,\left\langle \psi _i{\tau _{1,x}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{1,x} \right\rangle _{\varGamma _w}\\&\quad +\,\left\langle \psi _i{\tau _{2,x}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{2,x} \right\rangle _{\varGamma _w},\\ K^n_{u_xu_y}(i,j)&=(\partial _y\psi _i,\eta _h^{n+1}\partial _x\psi _j)+\left\langle \psi _i{\tau _{1,x}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{1,y} \right\rangle _{\varGamma _w}\\&\quad +\,\left\langle \psi _i{\tau _{2,x}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{2,y} \right\rangle _{\varGamma _w},\\ K^n_{u_xu_z}(i,j)&=\,(\partial _z\psi _i,\eta _h^{n+1}\partial _x\psi _j)+\left\langle \psi _i{\tau _{1,x}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{1,z} \right\rangle _{\varGamma _w}\\&\quad +\,\left\langle \psi _i{\tau _{2,x}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{2,z} \right\rangle _{\varGamma _w},\\ K^n_{u_yu_x}(i,j)&=\,(\partial _x\psi _i,\eta _h^{n+1}\partial _y\psi _j)+ \left\langle \psi _i{\tau _{1,y}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{1,x} \right\rangle _{\varGamma _w}\\&\quad +\, \left\langle \psi _i{\tau _{2,y}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{2,x} \right\rangle _{\varGamma _w},\\ K^n_{u_yu_y}(i,j)&=\,Re\left( \psi _i,\left( \frac{{\frac{1}{2}}(\rho _h^{n+1}+\rho _h^n)}{\delta t} +\rho _h^{n+1}(\mathbf{u}^n_h\cdot \nabla ) +{\frac{1}{2}}\left( \nabla \cdot (\rho _h^{n+1}\mathbf{u}^n_h)\right) \right) \psi _j\right) \\&\quad +\,(\nabla \psi _i,\eta _h^{n+1}\nabla \psi _j)+(\partial _y\psi _i,\eta _h^{n+1}\partial _y\psi _j)\\&\quad +\,\left\langle \psi _i{\tau _{1,y}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{1,y} \right\rangle _{\varGamma _w}\\&\quad +\,\left\langle \psi _i{\tau _{2,y}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{2,y} \right\rangle _{\varGamma _w}, \end{aligned}$$

$$\begin{aligned} K^n_{u_yu_z}(i,j)&=\,(\partial _z\psi _i,\eta _h^{n+1}\partial _y\psi _j)+\left\langle \psi _i{\tau _{1,y}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{1,z} \right\rangle _{\varGamma _w}\\&\quad +\, \left\langle \psi _i{\tau _{2,y}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{2,z} \right\rangle _{\varGamma _w},\\ K^n_{u_zu_x}(i,j)&=\,(\partial _x\psi _i,\eta _h^{n+1}\partial _z\psi _j)+\left\langle \psi _i{\tau _{1,z}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{1,x} \right\rangle _{\varGamma _w}\\&\quad +\, \left\langle \psi _i{\tau _{2,z}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{2,x} \right\rangle _{\varGamma _w},\\ K^n_{u_zu_y}(i,j)&=(\partial _y\psi _i,\eta _h^{n+1}\partial _z\psi _j)+\left\langle \psi _i{\tau _{1,z}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{1,y} \right\rangle _{\varGamma _w}\\&\quad +\, \left\langle \psi _i{\tau _{2,z}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{2,y} \right\rangle _{\varGamma _w},\\ K^n_{u_zu_z}(i,j)&=\,Re\left( \psi _i,\left( \frac{{\frac{1}{2}}(\rho _h^{n+1}+\rho _h^n)}{\delta t} +\rho _h^{n+1}(\mathbf{u}^n_h\cdot \nabla ) +{\frac{1}{2}}\left( \nabla \cdot (\rho _h^{n+1}\mathbf{u}^n_h)\right) \right) \psi _j\right) \\&\quad +\,(\nabla \psi _i,\eta _h^{n+1}\nabla \psi _j)+(\partial _z\psi _i,\eta _h^{n+1}\partial _z\psi _j)\\&\quad +\,\left\langle \psi _i{\tau _{1,z}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{1,z} \right\rangle _{\varGamma _w}\\&\quad +\,\left\langle \psi _i{\tau _{2,z}},\eta _h^{n+1}\left( {{\mathcal {L}}}_s{l_s}_h^{n+1}\right) ^{-1}\psi _j\tau _{2,z} \right\rangle _{\varGamma _w},\\ F_{u_x}^n(i)&=\,\frac{Re}{\delta t}(\rho _h^nu_{x,h}^n,\psi _i) + {{\mathcal {B}}}(\mu _h^{n+1}\partial _x\phi _h^{n+1},\psi _i)-(2\partial _xp_h^n-\partial _xp_h^{n-1},\psi _i)\\&\quad +\,\mathcal {B}\left\langle \left( \epsilon \partial _n\phi _h^{n+1}+Q\left( \phi ^{n+1}_h\right) +\tilde{\alpha }\left( \phi ^{n+1}_h-\phi ^n_h\right) \right) \partial _{\tau _1}\phi _h^{n+1},\psi _i\tau _{1,x}\right\rangle _{\varGamma _w}\\&\quad +\,\mathcal {B}\left\langle \left( \epsilon \partial _n\phi _h^{n+1}+Q\left( \phi ^{n+1}_h\right) +\tilde{\alpha }\left( \phi ^{n+1}_h-\phi ^n_h\right) \right) \partial _{\tau _2}\phi _h^{n+1},\psi _i\tau _{2,x}\right\rangle _{\varGamma _w},\\ F_{u_y}^n(i)&=\,\frac{Re}{\delta t}(\rho _h^nu_{y,h}^n,\psi _i) + {{\mathcal {B}}}(\mu _h^{n+1}\partial _y\phi _h^{n+1},\psi _i)-(2\partial _yp_h^n-\partial _yp_h^{n-1},\psi _i)\\&\quad +\,\mathcal {B}\left\langle \left( \epsilon \partial _n\phi _h^{n+1}+Q\left( \phi ^{n+1}_h\right) +\tilde{\alpha }\left( \phi ^{n+1}_h-\phi ^n_h\right) \right) \partial _{\tau _1}\phi _h^{n+1},\psi _i\tau _{1,y}\right\rangle _{\varGamma _w}\\&\quad +\,\mathcal {B}\left\langle \left( \epsilon \partial _n\phi _h^{n+1}+Q\left( \phi ^{n+1}_h\right) +\tilde{\alpha }\left( \phi ^{n+1}_h-\phi ^n_h\right) \right) \partial _{\tau _2}\phi _h^{n+1},\psi _i\tau _{2,y}\right\rangle _{\varGamma _w},\\ F_{u_z}^n(i)&=\,\frac{Re}{\delta t}(\rho _h^nu_{z,h}^n,\psi _i) + {{\mathcal {B}}}(\mu _h^{n+1}\partial _z\phi _h^{n+1},\psi _i)-(2\partial _zp_h^n-\partial _zp_h^{n-1},\psi _i)\\&\quad +\,\mathcal {B}\left\langle \left( \epsilon \partial _n\phi _h^{n+1}+Q\left( \phi ^{n+1}_h\right) +\tilde{\alpha }\left( \phi ^{n+1}_h-\phi ^n_h\right) \right) \partial _{\tau _1}\phi _h^{n+1},\psi _i\tau _{1,z}\right\rangle _{\varGamma _w}\\&\quad +\,\mathcal {B}\left\langle \left( \epsilon \partial _n\phi _h^{n+1}+Q\left( \phi ^{n+1}_h\right) +\tilde{\alpha }\left( \phi ^{n+1}_h-\phi ^n_h\right) \right) \partial _{\tau _2}\phi _h^{n+1},\psi _i\tau _{2,z}\right\rangle _{\varGamma _w}. \end{aligned}$$

In Step 4, we denote $\{\chi _i\}^{N_{q_h}}_{i=1}$ as the piecewise linear finite element basis functions for $M_h$ with $N_{q_h}$ as its dimension within an element. the element stiffness matrix $K^n_{p}$ and element right-hand side $F_p^n$ arising from (2.14) are

$$\begin{aligned}&K^n_{p}(i,j)=(\nabla \chi _i,\nabla \chi _j), \\&F_p^n(i)=-\frac{\bar{\rho }}{\delta t}Re(\nabla \cdot \mathbf{u}_h^{n+1},\chi _i)+ (\nabla p_h^n,\nabla \chi _i). \end{aligned}$$

Here $1\le i,j \le N_{q_h}$.