New stability estimates for an unfitted finite element method for two-phase Stokes problem

TL;DR

This paper develops new stability estimates for an unfitted finite element method applied to the two-phase Stokes problem, ensuring robustness with respect to viscosity jumps and providing theoretical and numerical validation.

Contribution

It introduces stability estimates for an unfitted finite element approach to the two-phase Stokes problem, with error bounds independent of viscosity variations.

Findings

Stability constants are independent of viscosity coefficients.

Finite element error estimates are robust against interface misalignment.

Numerical experiments confirm theoretical stability and accuracy.

Abstract

The paper addresses stability and finite element analysis of the stationary two-phase Stokes problem with a piecewise constant viscosity coefficient experiencing a jump across the interface between two fluid phases. We first prove a priori estimates for the individual terms of the Cauchy stress tensor with stability constants independent of the viscosity coefficient. Next, this stability result is extended to the approximation of the two-phase Stokes problem by a finite element method. In the method considered, the interface between the phases does not respect the underlying triangulation, which put the finite element method into the class of unfitted discretizations. The finite element error estimates are proved with constants independent of viscosity. Numerical experiments supporting the theoretical results are provided.

Tables

Table 1. Table 5.1. Example 1 with Ω − superscript Ω \Omega^{-} completely interior. Errors are shown for a fixed mesh, increasing values of ν + superscript 𝜈 \nu^{+} and decreasing values of ν − superscript 𝜈 \nu^{-} .

Parameters		$P_2$ – $P_0$		Mini-element
${\nu }^{-}$	${\nu }^{+}$	$𝚎\left(\text{u}\right)$	$𝚎\left(p\right)$	$𝚎\left(\text{u}\right)$	$𝚎\left(p\right)$
$1E-01$	$1E+01$	$0.01708$	$0.00448$	$0.06200$	$0.00409$
$1E-02$	$1E+02$	$0.01708$	$0.00448$	$0.06200$	$0.00409$
$1E-03$	$1E+03$	$0.01708$	$0.00448$	$0.06200$	$0.00409$
$1E-04$	$1E+04$	$0.01708$	$0.00448$	$0.06200$	$0.00409$

Table 2. Table 5.2. Example 1 with Ω + superscript Ω \Omega^{+} completely interior. Errors are shown for a fixed mesh, increasing values of ν + superscript 𝜈 \nu^{+} and decreasing values of ν − superscript 𝜈 \nu^{-} .

Parameters		$P_2$ – $P_0$		Mini-element
${\nu }^{-}$	${\nu }^{+}$	$𝚎\left(\text{u}\right)$	$𝚎\left(p\right)$	$𝚎\left(\text{u}\right)$	$𝚎\left(p\right)$
$1E-01$	$1E+01$	$0.01210$	$0.00282$	$0.06200$	$0.00401$
$1E-02$	$1E+02$	$0.01210$	$0.00282$	$0.06200$	$0.00401$
$1E-03$	$1E+03$	$0.01210$	$0.00282$	$0.06200$	$0.00401$
$1E-04$	$1E+04$	$0.01210$	$0.00282$	$0.06200$	$0.00401$

Table 3. Table 5.3. Example 2, errors for a sequence of uniform meshes, and fixed values of ν ± superscript 𝜈 plus-or-minus \nu^{\pm} . The solution ( u h , p h ) subscript u ℎ subscript 𝑝 ℎ (\text{\rm u}_{h},p_{h}) is approximated with P 2 subscript 𝑃 2 P_{2} – P 0 subscript 𝑃 0 P_{0} elements.

dofs	$𝚎\left(\text{u}\right)$	$𝚛\left(\text{u}\right)$	$𝚎\left(p\right)$	$𝚛\left(p\right)$	$𝚎(\text{u},p)$	$𝚛(\text{u},p)$
$815$	$5.9E-01$	$-$	$9.8E-01$	$-$	$1.1E+00$	$-$
$2867$	$2.6E-01$	$1.201$	$2.2E-01$	$2.157$	$3.4E-01$	$1.761$
$10867$	$1.1E-01$	$1.163$	$1.0E-01$	$1.106$	$1.5E-01$	$1.138$
$42199$	$5.2E-02$	$1.139$	$4.6E-02$	$1.154$	$6.9E-02$	$1.145$
$166303$	$2.5E-02$	$1.039$	$2.2E-02$	$1.035$	$3.4E-02$	$1.037$
$660243$	$1.3E-02$	$1.010$	$1.1E-02$	$0.998$	$1.7E-02$	$1.005$
$1481863$	$8.4E-03$	$1.004$	$7.5E-03$	$0.996$	$1.1E-02$	$1.000$

Table 4. Table 5.4. Example 2, errors for a sequence of uniform meshes, and fixed values of ν ± superscript 𝜈 plus-or-minus \nu^{\pm} . The solution ( u h , p h ) subscript u ℎ subscript 𝑝 ℎ (\text{\rm u}_{h},p_{h}) is approximated with the Mini-element.

dofs	$𝚎\left(\text{u}\right)$	$𝚛\left(\text{u}\right)$	$𝚎\left(p\right)$	$𝚛\left(p\right)$	$𝚎(\text{u},p)$	$𝚛(\text{u},p)$
$604$	$1.8E+00$	$-$	$9.5E-01$	$-$	$2.0E+00$	$-$
$2072$	$7.5E-01$	$1.226$	$2.0E-01$	$2.269$	$7.8E-01$	$1.362$
$7736$	$3.3E-01$	$1.189$	$6.0E-02$	$1.716$	$3.4E-01$	$1.213$
$29796$	$1.5E-01$	$1.121$	$2.0E-02$	$1.587$	$1.5E-01$	$1.132$
$116924$	$7.4E-02$	$1.043$	$7.0E-03$	$1.507$	$7.4E-02$	$1.049$
$463192$	$3.7E-02$	$1.012$	$2.5E-03$	$1.470$	$3.7E-02$	$1.015$
$1038836$	$2.4E-02$	$1.005$	$1.4E-03$	$1.460$	$2.4E-02$	$1.006$

Table 5. Table 5.5. Example 3, errors for a sequence of uniform meshes, and fixed values of ν ± superscript 𝜈 plus-or-minus \nu^{\pm} . The solution ( u h , p h ) subscript u ℎ subscript 𝑝 ℎ (\text{\rm u}_{h},p_{h}) is approximated with P 2 subscript 𝑃 2 P_{2} – P 0 subscript 𝑃 0 P_{0} elements.

dofs	$𝚎\left(\text{u}\right)$	$𝚛\left(\text{u}\right)$	$𝚎\left(p\right)$	$𝚛\left(p\right)$	$𝚎(\text{u},p)$	$𝚛(\text{u},p)$
$815$	$1.4E+00$	$-$	$5.8E+00$	$-$	$5.9E+00$	$-$
$2867$	$4.6E-01$	$1.568$	$2.9E-01$	$4.321$	$5.4E-01$	$3.445$
$10867$	$1.6E-01$	$1.484$	$8.1E-02$	$1.824$	$1.8E-01$	$1.564$
$42199$	$7.2E-02$	$1.192$	$3.3E-02$	$1.309$	$7.9E-02$	$1.213$
$166303$	$3.3E-02$	$1.130$	$1.3E-02$	$1.359$	$3.5E-02$	$1.164$
$660243$	$1.6E-02$	$1.033$	$6.0E-03$	$1.104$	$1.7E-02$	$1.042$
$1481863$	$1.1E-02$	$0.999$	$4.0E-03$	$0.962$	$1.1E-02$	$0.995$

Table 6. Table 5.6. Example 3, errors for a sequence of uniform meshes, and fixed values of ν ± superscript 𝜈 plus-or-minus \nu^{\pm} . The solution ( u h , p h ) subscript u ℎ subscript 𝑝 ℎ (\text{\rm u}_{h},p_{h}) is approximated with the Mini-element.

dofs	$𝚎\left(\text{u}\right)$	$𝚛\left(\text{u}\right)$	$𝚎\left(p\right)$	$𝚛\left(p\right)$	$𝚎(\text{u},p)$	$𝚛(\text{u},p)$
$604$	$3.0E+00$	$-$	$1.1E+01$	$-$	$1.1E+01$	$-$
$2072$	$1.5E+00$	$1.003$	$1.7E+00$	$2.639$	$2.3E+00$	$2.288$
$7736$	$5.9E-01$	$1.329$	$4.5E-01$	$1.927$	$7.4E-01$	$1.609$
$29796$	$2.4E-01$	$1.285$	$1.4E-01$	$1.717$	$2.8E-01$	$1.415$
$116924$	$1.1E-01$	$1.159$	$4.1E-02$	$1.732$	$1.2E-01$	$1.261$
$463192$	$5.2E-02$	$1.064$	$1.0E-02$	$2.014$	$5.3E-02$	$1.134$
$1038836$	$3.4E-02$	$1.021$	$4.8E-03$	$1.850$	$3.4E-02$	$1.043$

Table 7. Table 5.7. Example 4, errors for a sequence of uniform meshes, and fixed values of ν ± superscript 𝜈 plus-or-minus \nu^{\pm} . The solution ( u h , p h ) subscript u ℎ subscript 𝑝 ℎ (\text{\rm u}_{h},p_{h}) is approximated with P 2 subscript 𝑃 2 P_{2} – P 0 subscript 𝑃 0 P_{0} elements.

dofs	${𝚎}_d\left(\text{u}\right)$	${𝚛}_d\left(\text{u}\right)$	${𝚎}_{\infty }\left(\text{u}\right)$	${𝚛}_{\infty }\left(\text{u}\right)$	${𝚎}_{\infty }\left(p\right)$	${𝚛}_{\infty }\left(p\right)$
$1553$	$6.8E-02$	$-$	$3.7E-01$	$-$	$2.3E-01$	$-$
$5609$	$3.0E-02$	$1.183$	$1.7E-01$	$1.130$	$1.1E-01$	$1.045$
$21457$	$1.7E-02$	$0.847$	$1.1E-01$	$0.555$	$6.4E-02$	$0.828$
$83873$	$7.3E-03$	$1.184$	$4.7E-02$	$1.285$	$3.0E-02$	$1.079$
$331585$	$3.6E-03$	$1.024$	$2.1E-02$	$1.191$	$1.5E-02$	$1.026$
$1318473$	$1.7E-03$	$1.053$	$9.8E-03$	$1.066$	$7.4E-03$	$0.991$

Table 8. Table 5.8. Example 4, errors for a sequence of uniform meshes, and fixed values of ν ± superscript 𝜈 plus-or-minus \nu^{\pm} . The solution ( u h , p h ) subscript u ℎ subscript 𝑝 ℎ (\text{\rm u}_{h},p_{h}) is approximated with the Mini-element.

dofs	${𝚎}_d\left(\text{u}\right)$	${𝚛}_d\left(\text{u}\right)$	${𝚎}_{\infty }\left(\text{u}\right)$	${𝚛}_{\infty }\left(\text{u}\right)$	${𝚎}_{\infty }\left(p\right)$	${𝚛}_{\infty }\left(p\right)$
$1122$	$5.0E-01$	$-$	$1.7E+00$	$-$	$8.9E-01$	$-$
$3994$	$2.3E-01$	$1.124$	$9.4E-01$	$0.818$	$3.7E-01$	$1.277$
$15154$	$1.1E-01$	$1.055$	$4.8E-01$	$0.984$	$1.4E-01$	$1.432$
$58978$	$5.4E-02$	$1.025$	$2.4E-01$	$0.995$	$7.2E-02$	$0.916$
$232642$	$2.7E-02$	$1.013$	$1.2E-01$	$0.999$	$3.8E-02$	$0.905$
$923994$	$1.3E-02$	$1.006$	$6.0E-02$	$1.000$	$2.0E-02$	$0.964$
$2074098$	$8.9E-03$	$1.003$	$4.0E-02$	$1.000$	$1.3E-02$	$0.975$

Equations

$$\displaystyle-\text{\rm div}\big{(}\nu^{\pm}D(\text{\rm u})\big{)}+\nabla p^{\pm}$$

$$=g$$

$$[[\text{u}]]$$

$$[[\sigma (\text{u},p)\text{n}]]$$

$$=\text{0}$$

$$\Vert {\nu }^{\frac{1}{2}}D(\text{u}){\Vert }_{L^2(\Omega )}+\Vert {\nu }^{-\frac{1}{2}}p{\Vert }_{L^2(\Omega )}\le C\Vert {\nu }^{-\frac{1}{2}}\text{f}{\Vert }_{L^2(\Omega )},$$

$$\Vert {\nu }^{+}D({\text{u}}^{+}){\Vert }_{L^2({\Omega }^{+})}+\Vert p^{+}{\Vert }_{L^2({\Omega }^{+})}\le C(\Vert \text{f}{\Vert }_{L^2({\Omega }^{+})}+\sqrt{\frac{{\nu }^{+}}{{\nu }^{-}}}\Vert \text{f}{\Vert }_{L^2({\Omega }^{-})}).$$

$$\Vert \nu D(\text{u}){\Vert }_{L^2(\Omega )}+\Vert p{\Vert }_{L^2(\Omega )}\le C(\Vert \text{f}{\Vert }_{L^2(\Omega )}+\Vert \mathbf{\lambda }{\Vert }_{L^2(\Gamma )}+\Vert \nu g{\Vert }_{L^2(\Omega )}),$$

$$\begin{array}{rl}a(\text{u},\text{v})+b(\text{v},p)& =⟨\hat{\text{f}},\text{v}{⟩}_{-1}\forall \text{v}\in \text{V},\\ b(\text{u},q)& =-(g,q{)}_{\Omega }\forall q\in M,\end{array}$$

$$a(\text{u},\text{v}):=(\nu D(\text{u}),D(\text{v}){)}_{\Omega },b(\text{v},q):=-(\text{div}\text{v},q{)}_{\Omega },\text{and}\hat{\text{f}}\in {\text{V}}^{\ast }.$$

$$⟨\hat{\text{f}},\text{v}{⟩}_{-1}:={\int }_{\Omega }\text{f}\cdot \text{v}+{\int }_{\Gamma }\mathbf{\lambda }\cdot \text{v}.$$

$$\Vert \nu D(\text{u}){\Vert }_{\Omega }+\Vert p{\Vert }_{\Omega }\le C(\Vert \hat{\text{f}}{\Vert }_{-1}+\Vert \nu g{\Vert }_{\Omega }),$$

$$\Vert p^{\pm }-{\mathrm{avg}}_{{\Omega }^{\pm }}(p^{\pm }){\Vert }_{{\Omega }^{\pm }}\le C(\Vert {\nu }^{\pm }D({\text{u}}^{\pm }){\Vert }_{{\Omega }^{\pm }}+\Vert \hat{\text{f}}{\Vert }_{-1}),$$

$$\Vert p{\Vert }_{\Omega }\le C(\Vert \nu D(\text{u}){\Vert }_{\Omega }+\Vert \hat{\text{f}}{\Vert }_{-1}).$$

$$\begin{split}&\|p^{\pm}-\mathrm{avg}_{\Omega^{\pm}}(p^{\pm})\|_{\Omega^{\pm}}^{2}=(p^{\pm}-\mathrm{avg}_{\Omega^{\pm}}(p^{\pm}),\text{\rm div}\,\text{\rm v}^{\pm})_{\Omega^{\pm}}=(\text{\rm div}\,\text{\rm v}^{\pm},p^{\pm})_{\Omega}\\ &\qquad=\big{(}\nu^{\pm}D(\text{\rm u}^{\pm}),D(\text{\rm v}^{\pm})\big{)}_{\Omega^{\pm}}-\langle{\widehat{\text{\rm f}}},{\text{\rm v}^{\pm}}\rangle_{-1}\leq\big{(}\|\nu^{\pm}D(\text{\rm u}^{\pm})\|_{\Omega^{\pm}}+\|{\widehat{\text{\rm f}}}\|_{-1}\big{)}\|\text{\rm v}^{\pm}\|_{1,\Omega^{\pm}}\\ &\qquad\leq C\big{(}\|\nu^{\pm}D(\text{\rm u}^{\pm})\|_{\Omega^{\pm}}+\|{\widehat{\text{\rm f}}}\|_{-1}\big{)}\|p^{\pm}-\mathrm{avg}_{\Omega^{\pm}}(p^{\pm})\|_{\Omega^{\pm}},\end{split}$$

$$\|p\|_{\Omega}^{2}=(\text{\rm div}\,\text{\rm v},p)_{\Omega}=(\nu D(\text{\rm u}),D(\text{\rm v}))_{\Omega}-\langle{\widehat{\text{\rm f}}},\text{\rm v}\rangle_{-1}\leq C\|p\|_{\Omega}\big{(}\|\nu D(\text{\rm u})\|_{\Omega}+\|{\widehat{\text{\rm f}}}\|_{-1}\big{)},$$

$$\Vert {\nu }^{-}D({\text{u}}^{-}){\Vert }_{{\Omega }^{-}}\le C(\Vert \hat{\text{f}}{\Vert }_{-1}+\Vert \nu g{\Vert }_{\Omega }).$$

$$\Vert {\nu }^{\frac{1}{2}}D(\text{u}){\Vert }_{\Omega }^2=⟨\hat{\text{f}},\text{u}{⟩}_{-1}+(g,p{)}_{\Omega }.$$

$$⟨\hat{\text{f}},\text{u}{⟩}_{-1}\le C\Vert \hat{\text{f}}{\Vert }_{-1}\Vert D(\text{u}){\Vert }_{\Omega }\le \frac{C\Vert \hat{\text{f}}{\Vert }_{-1}}{\sqrt{{\nu }^{-}}}\Vert {\nu }^{\frac{1}{2}}D(\text{u}){\Vert }_{\Omega }$$

$$\begin{array}{rl}(g,p{)}_{\Omega }=& (g^{-},p^{-}-{\mathrm{avg}}_{{\Omega }^{-}}(p^{-}){)}_{{\Omega }^{-}}+(g^{-},{\mathrm{avg}}_{{\Omega }^{-}}(p^{-}){)}_{{\Omega }^{-}}\\ & +(g^{+},p^{+}-{\mathrm{avg}}_{{\Omega }^{+}}(p^{+}){)}_{{\Omega }^{+}}+(g^{+},{\mathrm{avg}}_{{\Omega }^{+}}(p^{+}){)}_{{\Omega }^{+}}.\end{array}$$

$$\begin{array}{rl}& (g^{-},p^{-}-{\mathrm{avg}}_{{\Omega }^{-}}(p^{-}){)}_{{\Omega }^{-}}+(g^{+},p^{+}-{\mathrm{avg}}_{{\Omega }^{+}}(p^{+}){)}_{{\Omega }^{+}}\\ & \le C(\sqrt{{\nu }^{-}}\Vert g^{-}{\Vert }_{{\Omega }^{-}}\Vert \sqrt{{\nu }^{-}}D({\text{u}}^{-}){\Vert }_{{\Omega }^{-}}+\sqrt{{\nu }^{+}}\Vert g^{+}{\Vert }_{{\Omega }^{+}}\Vert \sqrt{{\nu }^{+}}D({\text{u}}^{+}){\Vert }_{{\Omega }^{+}})+C\Vert g{\Vert }_{\Omega }\Vert \hat{\text{f}}{\Vert }_{-1}\\ & \le C(\Vert {\nu }^{\frac{1}{2}}g{\Vert }_{\Omega }\Vert {\nu }^{\frac{1}{2}}D(\text{u}){\Vert }_{\Omega }+\Vert g{\Vert }_{\Omega }\Vert \hat{\text{f}}{\Vert }_{-1}).\end{array}$$

$$\begin{array}{rl}& (g^{-},{\mathrm{avg}}_{{\Omega }^{-}}(p^{-}){)}_{{\Omega }^{-}}+(g^{+},{\mathrm{avg}}_{{\Omega }^{+}}(p^{+}){)}_{{\Omega }^{+}}\\ & =({\mathrm{avg}}_{{\Omega }^{-}}(g^{-}),{\mathrm{avg}}_{{\Omega }^{-}}(p^{-}){)}_{{\Omega }^{-}}+({\mathrm{avg}}_{{\Omega }^{+}}(g^{+}),{\mathrm{avg}}_{{\Omega }^{+}}(p^{+}){)}_{{\Omega }^{+}}\\ & =-\frac{|{\Omega }^{+}|}{|{\Omega }^{-}|}({\mathrm{avg}}_{{\Omega }^{+}}(g^{+}),{\mathrm{avg}}_{{\Omega }^{-}}(p^{-}){)}_{{\Omega }^{-}}+({\mathrm{avg}}_{{\Omega }^{+}}(g^{+}),{\mathrm{avg}}_{{\Omega }^{+}}(p^{+}){)}_{{\Omega }^{+}}.\end{array}$$

$$\begin{array}{rl}& (g^{-},{\mathrm{avg}}_{{\Omega }^{-}}(p^{-}){)}_{{\Omega }^{-}}+(g^{+},{\mathrm{avg}}_{{\Omega }^{+}}(p^{+}){)}_{{\Omega }^{+}}\le C\Vert g^{+}{\Vert }_{{\Omega }^{+}}\Vert p{\Vert }_{\Omega }\le C\Vert g^{+}{\Vert }_{{\Omega }^{+}}(\Vert \nu D(\text{u}){\Vert }_{\Omega }+\Vert \hat{\text{f}}{\Vert }_{-1})\\ & \le C\sqrt{{\nu }^{+}}\Vert g^{+}{\Vert }_{{\Omega }^{+}}\Vert {\nu }^{\frac{1}{2}}D(\text{u}){\Vert }_{\Omega }+C\Vert g{\Vert }_{\Omega }\Vert \hat{\text{f}}{\Vert }_{-1},\end{array}$$

$$(g,p{)}_{\Omega }\le C(\Vert {\nu }^{\frac{1}{2}}g{\Vert }_{\Omega }\Vert {\nu }^{\frac{1}{2}}D(\text{u}){\Vert }_{\Omega }+\Vert g{\Vert }_{\Omega }\Vert \hat{\text{f}}{\Vert }_{-1}).$$

$$\|\nu^{\frac{1}{2}}D(\text{\rm u})\|_{\Omega}\leq C\bigg{(}\frac{1}{\sqrt{\nu^{-}}}\|{\widehat{\text{\rm f}}}\|_{-1}+\|\nu^{\frac{1}{2}}g\|_{\Omega}\bigg{)}.$$

$$\Vert p^{-}-{\mathrm{avg}}_{{\Omega }^{-}}(p^{-}){\Vert }_{{\Omega }^{-}}\le C(\Vert \hat{\text{f}}{\Vert }_{-1}+\Vert \nu g{\Vert }_{\Omega }).$$

$$E\text{w}\in {\mathbf{H}}^1(\Omega ),E\text{w}=\text{w}\text{ in }{\Omega }^{+},\text{and}\Vert E\text{w}{\Vert }_{1,\Omega }\le C\Vert \text{w}{\Vert }_{1,{\Omega }^{+}},$$

$$D(\text{w})=\text{0},\text{and hence}\text{div}\text{w}=\text{tr}(D(\text{w}))=0.$$

$$\|\nu^{+}D(\text{\rm u}^{+})\|_{\Omega^{+}}\leq C\big{(}\|{\widehat{\text{\rm f}}}\|_{-1}+\|\nu g\|_{\Omega}\big{)}.$$

$$\text{v}=\{\begin{array}{ll}E{\text{u}}^{+},& \text{if }|\partial {\Omega }^{+}\cap \partial \Omega |>0,\\ E({\text{u}}^{+}-{\mathcal{P}}_{{\Omega }^{+}}^{RM}({\text{u}}^{+})),& \text{if }|\partial {\Omega }^{+}\cap \partial \Omega |=0.\end{array}$$

$$\Vert \text{v}{\Vert }_{1,\Omega }\le C\Vert D({\text{u}}^{+}){\Vert }_{{\Omega }^{+}}.$$

$${\nu }^{+}\Vert D({\text{u}}^{+}){\Vert }_{{\Omega }^{+}}^2=⟨\hat{\text{f}},\text{v}{⟩}_{-1}-({\nu }^{-}D({\text{u}}^{-}),D({\text{v}}^{-}){)}_{{\Omega }^{-}}+(\text{div}\text{v},p{)}_{\Omega },$$

Full text

Partially supported by NSF through the Division of Mathematical Sciences grant 1620100

Partially supported by NSF through the Division of Mathematical Sciences grant 1717516.

1. Introduction

We are interested in the analysis and a finite element method for the two-phase Stokes problem (also known in the literature as the Stokes interface problem). The system of equations is posed in a bounded Lipschitz domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, decomposed in two subdomains (phases) $\Omega^{\pm}$. The interface $\Gamma$ between two phases is a closed hypersurface immersed in $\Omega$, i.e., $\Gamma\subset\Omega$ and $\Gamma=\overline{\Omega^{+}}\cap\overline{\Omega^{-}}$. We assume $\Gamma$ is Lipschitz smooth. The Stokes interface problem reads as follows: Given a force field $\text{\rm f}\in L^{2}(\Omega)^{d}$, a source term $g\in L^{2}(\Omega)$, an interface force $\bm{\lambda}\in L^{2}(\Gamma)^{d}$, and viscosity coefficient $\nu^{\pm}$ constant and positive in each subdomain, find the fluid velocity u and the normalized kinematic pressure $p$ such that

[TABLE]

where $D(\text{\rm u}):=\frac{1}{2}(\nabla\text{\rm u}+(\nabla\text{\rm u})^{T})$ is the rate-of-strain tensor, $\sigma(\text{\rm u},p)=\nu\,D(\text{\rm u})-p\,\mathbb{I}$ is the Cauchy stress tensor, and n is a unit vector on $\Gamma$ pointing from $\Omega^{+}$ to $\Omega^{-}$. For any $f\in L^{1}(\Omega)$ we use notations $f^{\pm}$ for the restriction of $f$ on $\Omega^{\pm}$, i.e., $f^{\pm}=f|_{\Omega^{\pm}}$; same convention is used for vector functions. The jumps on the interface are then defined as $\llbracket{\sigma(\text{\rm u},p)\text{\rm n}}\rrbracket=\sigma(\text{\rm u}^{+},p^{+})\text{\rm n}-\sigma(\text{\rm u}^{-},p^{-})\text{\rm n}$ and $\llbracket{\text{\rm u}}\rrbracket=\text{\rm u}^{+}-\text{\rm u}^{-}$.

The studies of the Stokes interface problem are motivated by continuum models of two-phase flows. If the fluid is treated as Newtonian incompressible with immiscible phases separated by the sharp interface, then the system (1.1) is a reasonable model problem for the limit case of highly viscous fluid; see, e.g., [19, 6, 17, 18]. It also appears as an auxiliary problem in numerical simulations of two-phase incompressible flows [8]. According to the continuum surface force model, cf. [3], the effect of interfacial forces, such as the surface tension, are taken into account by using a localized force term at the interface, i.e., $\bm{\lambda}$ in (1.1).

Problem (1.1) is linear and a standard weak formulation (2.1) renders it as a saddle-point problem, thus yielding the well-posedness result and leading to Galerkin numerical methods; see, e.g., [7, 4]. This textbook analysis, however, does not provide an explicit information on the dependence of the stability and numerical errors estimates on the viscosity coefficient, in particular, on the ratio $\nu^{+}/\nu^{-}$, provided $\nu^{-}\leq\nu^{+}$. This robustness question becomes important if one addresses numerical stability of Galerkin methods, such as the finite element method, for the case of high variation in viscosity coefficient between two phases. The $\nu$-dependence of stability and finite element error estimates for (1.1) have been studied in the literature only recently; see [14, 13, 10, 11]. In those studies, stability and error analysis was done for the natural energy norm of the problem. In particular, under certain further assumptions on $\Omega^{\pm}$, the a priori estimate from [14] (proved there for $g=0$, $\bm{\lambda}=0$) reads

[TABLE]

with $C$ independent of $\nu$. Note that for single phase Stokes problem, a simple scaling argument provides uniform estimates for the quantities $\nu\,\text{\rm u}$ and $p$. Similar result does not follow from (1.2) for the velocity and pressure in each of the phases. For example, for $\nu^{+}\text{\rm u}^{+}$ and $p^{+}$ the estimate (1.2) yields

[TABLE]

We see that the right-hand side blows up for $\nu^{-}\to 0$. In the present paper, we prove the following stability result for the solution of (1.1):

[TABLE]

The improvement over (1.3) is clear: the re-scaled solution components, $\nu^{\pm}\text{\rm u}^{\pm}$ and $p^{\pm}$, enjoy uniform estimates in the corresponding subdomains, just as for the single phase problem. The estimate (1.4) can be also seen as the uniform estimate for the components of Cauchy stress tensor, an important quantity in practical fluid mechanics.

In the same spirit as (1.4) improves over the energy estimate (1.2), the finite element analysis developed in this paper extends the existing one by deriving robust in $\nu$ stability estimates and error estimates for the components of the finite element Cauchy stress tensor. Following [10], for the discretization of (1.1) we consider a geometrically unfitted finite element method known as Nitsche-XFEM or cutFEM. Geometrically unfitted methods use a fixed background mesh which does not respect the position of the interface. The main advantage of unfitted FEM is the relative ease of handling time-dependent domains, implicitly defined interfaces and problems with strong geometric deformations [2]. We prove uniform with respect to $\nu$ stability and error estimates for the unfitted FEM. These results hold for a family of bulk LBB-stable finite element Stokes pairs defined on the background mesh. These pairs include $P_{k+1}-P_{k}$, $k\geq 1$, and $P_{k+d}-P_{k}^{\rm disc}$ for $k\geq 0$, $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, and several other elements. We are able to accomplish this by combining ideas from the two papers [9, 5]. In [9] an unfitted FEM for the single phase Stokes problem was analyzed. We borrow some crucial inf-sup stability estimates from that paper. In [5] similar stability results were proved the Poisson interface problem. The chief tool was to use extension operators in Sobolev spaces. Similarly, here extension operators are essential, however, the pressure terms and the div-free condition add several new difficulties.

We organized the paper in five sections. In section 2 a notion of the weak solution is introduced and estimate (1.3) is proved. Section 3 describes the finite element method and proves the analogue of (1.3) for the finite element solution. In section 4 a $\nu$-independent optimal order error estimate is proved. Finally, a few illustrative results of numerical experiments are given in section 5.

2. A priori analysis for (1.1)

2.1. Preliminaries and problem setting

We introduce a variational formulation of (1.1) and several notations to be used throughout the paper. For an open set $\mathcal{O}\subset\mathbb{R}^{d}$ denote by $(\cdot,\cdot)_{\mathcal{O}}$ the $L^{2}$ inner product in $\mathcal{O}$, and by $\|\cdot\|_{\mathcal{O}}$ the corresponding norm. For the mixed variational formulation of (1.1), we set $\text{\rm V}:=[\mathrm{H}_{0}^{1}(\Omega)]^{d}$ for the space the vector field u belongs to, whereas for the pressure $p$ we set $M=L_{0}^{2}(\Omega)$, with $L_{0}^{2}(\mathcal{O})=\{p\in L^{2}(\mathcal{O}):(p,1)_{\mathcal{O}}=0\}$. We let $\|\cdot\|_{1,\mathcal{O}}$ denote the $H^{1}({\mathcal{O}})$-norm.

The norm of $\text{\rm V}^{*}$, the dual of V, is denoted by $\|\cdot\|_{-1}$, and $\langle\cdot,\cdot\rangle_{-1}$ denotes the pairing, with respect to the $L^{2}$-duality.

We consider the abstract mixed formulation: Find $(\text{\rm u},p)\in\text{\rm V}\times M$ such that

[TABLE]

where

[TABLE]

The problem (2.1) is the weak formulation of the Stokes interface equation (1.1) if we let

[TABLE]

2.2. Stability estimates for the weak solution

In this section, we analyze the variational formulation (2.1) of the Stokes interface problem (1.1). We are interested in the following stability result for the solution $(\text{\rm u},p)\in\text{\rm V}\times M$ of (2.1):

[TABLE]

with a constant $C>0$ independent of $\nu$ and depending only on $\Omega$ and $\Gamma$. A standard energy argument gives estimates for $\|\sqrt{\nu}\,D(\text{\rm u})\|_{\Omega}$ with the corresponding constant $C$ on the right-hand side dependent on $\nu$. Without any loss of generality we shall always assume $\nu^{-}\leq\nu^{+}$. The key to get the improved result (2.3) is using an energy argument to estimate $\|\nu^{-}D(\text{\rm u})\|_{\Omega^{-}}$, and then using an extension operator to estimate $\|\nu^{+}D(\text{\rm u})\|_{\Omega^{+}}$. This strategy is similar to the one taken in [5] for the Poisson interface problem. However, here the pressure term and the divergence condition require careful treatment, by a repetitive use of a continuous inf-sup condition.

In what follows, for an open set $\mathcal{O}$ and a function $q\in L^{2}(\mathcal{O})$, we denote its average over $\mathcal{O}$ by the expression $\mathrm{avg}_{\mathcal{O}}(q)=|\mathcal{O}|^{-1}(q,1)_{\mathcal{O}}$. We start by proving the following results.

Lemma 2.1**.**

Let $\text{\rm u}\in\text{\rm V}$ and $p\in M$ solve (2.1). Then there exists $C>0$, depending only on $\Omega$ and $\Gamma$, such that

[TABLE]

[TABLE]

Proof.

Since $p^{\pm}-\mathrm{avg}_{\Omega^{\pm}}(p^{\pm})\in L_{0}^{2}(\Omega^{\pm})$, the result from [1] ensures the existence of $\text{\rm v}^{\pm}\in\text{\rm H}_{0}^{1}(\Omega^{\pm})$ such that $\text{\rm div}\,\text{\rm v}^{\pm}=p^{\pm}-\mathrm{avg}_{\Omega^{\pm}}(p^{\pm})$ in $\Omega^{\pm}$, and $\|\text{\rm v}^{\pm}\|_{1,\Omega^{\pm}}\leq C\|p^{\pm}-\mathrm{avg}_{\Omega^{\pm}}(p^{\pm})\|_{\Omega}$, for some $C>0$ depending only on $\Omega$ and $\Gamma$. Therefore, using the first equation of (2.1) with $\text{\rm v}^{\pm}$ in $\Omega^{\pm}$ extended by zero on $\Omega^{\mp}$ and since $(\text{\rm div}\,\text{\rm v}^{\pm},1)_{\Omega^{\pm}}=0$, we get

[TABLE]

which gives (2.4). On the other hand, since $p\in M$, there exists $\text{\rm v}\in\text{\rm V}$ such that $\text{\rm div}\,\text{\rm v}=p$ in $\Omega$, and $\|\text{\rm v}\|_{1,\Omega}\leq C\|p\|_{\Omega}$, for some $C>0$, depending only on $\Omega^{\pm}$. Therefore, using the first equation of (2.1) with $\text{\rm v}\in\text{\rm V}$, we get

[TABLE]

which proves (2.5). ∎

In order to prove (2.3), we start by estimating $\nu^{-}D(\text{\rm u}^{-})$. This is an easier part of the desired estimate (2.3) since $\nu^{-}\leq\nu^{+}$.

Lemma 2.2**.**

Let $\text{\rm u}\in\text{\rm V}$ and $p\in M$ solve (2.1), then there exists $C>0$, depending only on $\Omega$, such that

[TABLE]

Proof.

The result follows from basic energy arguments and pressure estimates provided by Lemma 2.1. To see this, let us test the first equation of (2.1) with $\text{\rm v}=\text{\rm u}$ and the second equation with $q=p$ to get

[TABLE]

We start by estimating the first term on the right-hand side

[TABLE]

where we used Poincare’s inequality, Korn’s inequality and the fact that $0<\nu^{-}\leq\nu^{+}$. To estimate the second term, we use the following decomposition

[TABLE]

By (2.4) we have

[TABLE]

From $g\in M$ it follows that $(g,1)_{\Omega}=1$ and so $|\Omega^{+}|\mathrm{avg}_{\Omega^{+}}(p^{+})=-|\Omega^{-}|\mathrm{avg}_{\Omega^{-}}(p^{-})$. We employ this equality below to obtain

[TABLE]

Hence, we have

[TABLE]

where we used (2.5) and the fact that $\nu^{-}\leq\nu^{+}$. Thus, we have shown that

[TABLE]

Combining (2.7), (2.8), (2.9) and using that $\nu^{-}\leq\nu^{+}$ we arrive at

[TABLE]

This implies the result due to $\nu^{-}\leq\nu^{+}$. ∎

An immediate consequence of Lemma 2.2 and (2.4) is the desired pressure estimate in $\Omega^{-}$.

Lemma 2.3**.**

Let $\text{\rm u}\in\text{\rm V}$ and $p\in M$ solve (2.1), then there exists $C>0$, depending only on $\Omega$, such that

[TABLE]

It remains to show analogues bound for $\nu D(\text{\rm u})$ and $p$ on $\Omega^{+}$. To estimate $\nu D(\text{\rm u})$ in $\Omega^{+}$, we consider an extension operator. A detailed construction of this operator can be found in [16, Chapter VI, Section 3.3]. Let $E:\text{\rm H}^{1}(\Omega^{+})\to\text{\rm H}^{1}(\Omega)$ be the bounded extension operator from $\Omega^{+}$ to $\Omega$. That is, if $\text{\rm w}\in\mathbf{H}^{1}(\Omega^{+})$, then

[TABLE]

with a constant $C>0$, depending only on $\Omega^{+}$ and $\Omega$. Moreover, one can assume $E\text{\rm w}$ to vanish on $\partial\Omega$, i.e., $E\text{\rm w}\in\text{\rm V}$ for $\text{\rm w}\in\text{\rm V}$. Further, we note that if $\text{\rm w}(\text{\rm x}):=\text{\rm a}+\text{\rm B}\text{\rm x}\in\mathbb{RM}(\Omega)$, then

[TABLE]

Here $\mathbb{RM}(\mathcal{O})$ is the space of rigid body motions defined on $\mathcal{O}$. We denote by $\mathcal{P}^{RM}_{\mathcal{O}}(\text{\rm v})$ the $L^{2}$-orthogonal projection of $\text{\rm v}\in L^{2}(\mathcal{O})^{d}$ onto the subspace of rigid body motions.

Lemma 2.4**.**

Let $\text{\rm u}\in\text{\rm V}$ and $p\in M$ solve (2.1), then there exists $C>0$, depending only on $\Omega$ and $\Gamma$, such that

[TABLE]

Proof.

Define $\text{\rm v}\in\text{\rm V}$ by

[TABLE]

The boundeness of the extension operator in (2.10), Poincare’s and Korn’s inequalities imply

[TABLE]

The first equation of (2.1) with the above v gives

[TABLE]

where we used that $D(\text{\rm v}^{+})=D(\text{\rm u}^{+})$. We now bound each term on the right-hand side of (2.13). Using (2.12) we have

[TABLE]

Thanks to (2.12) and (2.6) we bound the second term on the right hand side of (2.13)

[TABLE]

For the third term we use decomposition:

[TABLE]

Using (2.3) and (2.12) we estimate

[TABLE]

For the next term in (2.16) we use $(\text{\rm div}\,\text{\rm v},1)_{\Omega}=0$ and the second equation in (1.1) to write

[TABLE]

Hence, with the help of (2.5) and (2.6) we get

[TABLE]

Similarly, for the last term in (2.16) we have,

[TABLE]

We combine the last two estimates to obtain the bound

[TABLE]

This estimate together with (2.15), (2.14), (2.13) leads to another bound

[TABLE]

This implies (2.11) after multiplication all through by $\nu^{+}$ and doing simple computations. ∎

Collecting (2.5), (2.6) and (2.11) we obtain the main result of this section.

Theorem 2.1**.**

Let $\text{\rm u}\in\text{\rm V}$ and $p\in M$ solve (2.1), then there exists $C>0$, depending only on $\Omega$ and $\Gamma$, such that

[TABLE]

The second inequality in (2.17) follows from the definition of the functional $\widehat{\text{\rm f}}$ in (2.2), the Poincaré inequality and the trace inequality, $\|\text{\rm v}\|_{\Gamma}\leq C\|\text{\rm v}\|_{1,\Omega}$.

3. The finite element method

3.1. Preliminaries and problem setting

For the discretization purpose, we assume that $\Omega$ is polygonal/polyhedral. Let $\{\mathcal{T}_{h}\}_{h>0}$ be an admissible family of triangulations of $\Omega$. We adopt the convention that the elements $T$ and edges $e$ are open sets, and use over-line symbol to refer to their closure. For each simplex $T\in\mathcal{T}_{h}$, let $h_{T}$ denote its diameter and define the global parameter of the triangulation by $h=\max_{T}h_{T}$. We assume that $\mathcal{T}_{h}$ is shape regular, that is, there exists $\kappa>0$ such that for every $T\in\mathcal{T}_{h}$, the radius $\rho_{T}$ of the inscribed sphere satisfies

[TABLE]

The sets of elements intersecting $\Omega^{\pm}$ and the set of elements cutting the interface $\Gamma$ are of interest. These are defined by

[TABLE]

We also define the sets of elements interior to each of subdomains $\Omega^{\pm}$, $\mathcal{T}_{h,i}^{\pm}=\{T\in\mathcal{T}_{h}^{\pm}:T\subset\Omega^{\pm}\}$. Finally, we let $\mathcal{E}_{h,i}^{\pm}$ be the collection of $d-1$, sub-simplexes of $\mathcal{T}_{h,i}^{\pm}$ (faces for $d=3$ and edges for $d=2$).

For $T\in\mathcal{T}_{h}^{\Gamma}$, we denote $T_{\Gamma}:=\overline{T}\cap\Gamma$. Under these definitions we define the $h$-dependent domains

[TABLE]

In particular, using the definition of the sets $\mathcal{T}_{h}^{+}$ and $\mathcal{T}_{h,i}^{-}$, we have that $\overline{\Omega}=\overline{\Omega_{h,i}^{-}}\cup\overline{\Omega_{h}^{+}}$. This fact will be useful when constructing a discrete extension operator. We also consider the layer of elements cut by the interface:

[TABLE]

and define the set of faces (edges) of $\mathcal{T}_{h}^{\Gamma}$ restricted to the interior of $\Omega_{h}^{\pm}$:

[TABLE]

For a piecewise smooth vector valued function v, the jump across an interior face $e=\mathrm{int}(\partial T_{1}\cap\partial T_{2})$ is defined by $\llbracket{\text{\rm v}}\rrbracket=\text{\rm v}|_{T_{1}}\cdot\text{\rm n}_{1}+\text{\rm v}|_{T_{2}}\cdot\text{\rm n}_{2}$, where $\text{\rm n}_{1}$ and $\text{\rm n}_{2}$ are the unit normal vectors to $e$, pointing outwards to $T_{1}$ and $T_{2}$, respectively. For a scalar function, we define $\llbracket{p}\rrbracket=p|_{T_{1}}\text{\rm n}_{1}+p|_{T_{2}}\text{\rm n}_{2}$.

The space of discontinuous and continuous finite element pressures are given by

[TABLE]

where integer $k_{p}\geq 0$ is a fixed polynomial degree. Throughout this paper, ${M_{h}^{\rm bulk}}=M_{h,\rm disc}$ for $k_{p}\geq 0$ or ${M_{h}^{\rm bulk}}=M_{h,\rm cont}$ for $k_{p}\geq 1$. We define $M_{h}^{\pm}:={M_{h}^{\rm bulk}}\cap L^{2}(\Omega_{h}^{\pm})$. Finally, our pressure space is given by

[TABLE]

Note that every element from $\mathcal{T}_{h}^{\Gamma}$ supports two finite element pressures corresponding to different phases. Only the restriction of these pressures to $\Omega^{+}$ or $\Omega^{+}$, respectively, makes sense as a numerical approximation of the true pressure solving the original problem (1.1). Same comment will be valid for the finite element velocity fields defined next. We consider the vector finite element space for $k\geq 1$,

[TABLE]

Next we consider a background velocity finite element space ${\text{\rm V}_{h}^{\rm bulk}}$ such that

[TABLE]

for some integers $s\geq k_{u}\geq 1$. Let $\text{\rm V}_{h}^{\pm}:={\text{\rm V}_{h}^{\rm bulk}}\cap\text{\rm H}^{1}(\Omega_{h}^{\pm})$. Finally our velocity space will be

[TABLE]

Also, we denote a generic element $\text{\rm v}_{h}\in\text{\rm V}_{h}$ by $\text{\rm v}_{h}:=(\text{\rm v}_{h}^{-},\text{\rm v}_{h}^{+})$.

Functions from $M_{h}$, $\text{\rm V}_{h}$ and their derivatives are multivalued in $\omega_{h}$, the overlap of $\Omega_{h}^{+}$ and $\Omega_{h}^{-}$. Below we use the $L^{2}(\Omega)$ norm notion for such functions to denote the norm of single-valued functions obtained by restricting $\Omega_{h}^{+}$-components on $\Omega^{+}\subset\Omega_{h}^{+}$ and $\Omega_{h}^{-}$-components on $\Omega^{-}\subset\Omega_{h}^{-}$. For example,

[TABLE]

for $p_{h}\in M_{h}$, $\text{\rm u}_{h}\in\text{\rm V}_{h}$, and so forth. The jump of a multivalued function over the interface is defined as the difference of components coming from $\Omega_{h}^{+}$ and $\Omega_{h}^{-}$, i.e. $\llbracket{\text{\rm v}_{h}}\rrbracket=\text{\rm v}_{h}^{+}-\text{\rm v}_{h}^{-}$ on $\Gamma$.

We consider a discrete norm $\|\cdot\|_{1,h}$, defined on $\text{\rm V}_{h}$, as follows: for all $\text{\rm v}_{h}\in\text{\rm V}_{h}$,

[TABLE]

where $\partial_{\text{\rm n}}^{\ell}$ is the $\ell-$th order normal derivative. We define a scaled norm

[TABLE]

and the augmented scaled norm

[TABLE]

We use the notation $\|\cdot\|_{\pm,h}$ to define a discrete norm on $M_{h}$, for all $q_{h}=(q_{h}^{-},q_{h}^{+})\in M_{h}$, by

[TABLE]

Spaces ${\text{\rm V}_{h}^{*}}$ and $M_{h}^{*}$ are dual to $\text{\rm V}_{h}$ and $M_{h}^{\pm}$ with respect to $\|\cdot\|_{1,h}$ and $\|\cdot\|_{\pm,h}$. For $\text{\rm F}_{h}\in{\text{\rm V}_{h}^{*}}$ and $G_{h}\in M_{h}^{*}$, we have by definition that

[TABLE]

where we agree that $\sup_{x}$ is taken over non-zero elements, if $x$ appears in the denominator. In this paper, we only consider bulk spaces ${\text{\rm V}_{h}^{\rm bulk}},{M_{h}^{\rm bulk}}$ which form inf-sup stable pairs. This is listed as an assumption.

Assumption 1**.**

There exists a constant $\xi>0$ such that

[TABLE]

We end this section by considering further assumptions on the mesh and on the pair of spaces $\{\text{\rm V}_{h},M_{h}\}$. These assumptions are essentially the ones made in [9]. For a generic set of elements $\mathcal{T}\subset\mathcal{T}_{h}$, denote $\omega(\mathcal{T})\subset\mathcal{T}_{h}$ the set of all tetrahedra having at least one vertex in $\mathcal{T}$.

Assumption 2**.**

For any $T\in\mathcal{T}_{h}^{\Gamma}$ we assume that the sets $W^{{\pm}}(T)=\mathcal{T}_{h,i}^{\pm}\cap\omega(\omega(T))$ are not empty.

We note that this assumption can be weaken by allowing in $W^{\pm}(T)$ neighbors of $T$ of degree $L$, with some finite and mesh independent $L\geq 2$.

Given $T\in\mathcal{T}_{h}^{\Gamma}$, we associate arbitrary but fixed $K_{T}^{{\pm}}\in W^{{\pm}}(T)$, which can be reached from $T$ by crossing faces in $\mathcal{E}_{h}^{\Gamma,\pm}$. More precisely, there exists simplexes $T=K_{1}^{{\pm}},K_{2}^{{\pm}},\ldots,K_{M}^{{\pm}}=K_{T}^{{\pm}}$ with $K_{j}^{{\pm}}\in\mathcal{T}_{h}^{\Gamma}$ for $j<M$. The number $M$ is uniformly bounded and only depends on the shape regularity of the mesh. Moreover, note that by (3.1) there exists a constant $c$ only depending on the shape regularity constant $\kappa$ such that

[TABLE]

For $T\in\mathcal{T}_{h,i}^{\pm}$, we set $K_{T}^{{\pm}}=T$.

Assumption 3**.**

Let $F\in\mathcal{E}_{h}^{\Gamma,{\pm}}$, with $F=\partial T_{1}\cap\partial T_{2}$. Assume $K_{T_{2}}^{\pm}$ can be reached from $K_{T_{1}}^{\pm}$ by crossing a finite, independent of $h$, number of faces of tetrahedra from $\mathcal{T}_{h,i}^{\pm}$.

The following assumption is also a type of inf-sup condition but restricted to interior elements in two phases, i.e. those lying inside $\Omega_{h,i}^{\pm}$.

Assumption 4**.**

There exists a constant $\beta>0$ such that

[TABLE]

where

[TABLE]

Examples of pair of spaces ${\text{\rm V}_{h}^{\rm bulk}}$–${M_{h}^{\rm bulk}}$ that satisfy the Assumption 4 can be found in [9, Section 6]; they include $P_{k+1}-P_{k}$, $k\geq 1$, and $P_{k+d}-P_{k}^{\rm disc}$ for $k\geq 0$, $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, and several other elements. In particular, if a pair ${\text{\rm V}_{h}^{\rm bulk}}$–${M_{h}^{\rm bulk}}$ satisfies these assumptions, we have the following result.

Theorem 3.1**.**

Suppose Assumptions 2–4 hold. There exists a constant $\theta>0$, independent of $h$ and $q$, and a constant $h_{0}>0$ such that for all $q\in L_{0}^{2}(\Omega_{h,i}^{\pm})\cap M_{h}^{\rm bulk}$ and $h\leq h_{0}$, we have

[TABLE]

Proof.

See [9, Theorem 1]. ∎

We will also need trace and inverse inequalities which can be found, for example, in [9].

Lemma 3.1**.**

Let $T\in\mathcal{T}_{h}$, then it holds

[TABLE]

with a constant $C$ independent of $T$ and on how $\Gamma$ intersects $T$.

3.2. Discrete Extension Operator

As in the continuous setting, a chief tool will be an extension operator. In this section we provide a discrete analogue of the extension operator in (2.10). Widlund [20] provided a discrete extension operator when the mesh fits the interface. For non-fitted meshes, as is the case here, a discrete extension operator can be found in [5] for piecewise-linear finite element functions and smooth interface.

In Lemma 3.2 below, we prove the result for unfitted meshes and finite elements of arbitrary degree that admit the existence of a local nodal basis. Let us make this assumption precise. The velocity bulk space is the space of vector functions, i.e. $\mathbf{V}_{h}^{\rm bulk}=\otimes_{j=1}^{d}V_{j,h}$. For each component space $V_{j,h}$ we assume that (i) there is a set of points (nodes) $\mathcal{N}(\mathcal{T}_{h})=\{\text{\rm y}_{1},\ldots,\text{\rm y}_{\ell}\}$ such that $v\in V_{j,h}$ is uniquely determined by $v(\text{\rm y}_{i})$ for $1\leq i\leq\ell$; (ii) for each $T$ there exists a local subset $\mathcal{N}(T)=\{\text{\rm y}\in\mathcal{N}(\mathcal{T}_{h})\,:\,\text{\rm y}\in\overline{T}\}$ such that $v|_{T}$ is uniquely determined by the values there; and (iii) if $\Phi:T\to\widehat{T}$ is the affine mapping to the reference simplex, then $\Phi(\mathcal{N}(T))$ is independent of $T\in\mathcal{T}_{h}$. The existence of the local nodal basis is, of course, standard if $\mathbf{V}_{h}^{\rm bulk}=\text{\rm W}_{h}^{k}$ for some integer $k\geq 1$. Moreover, we assume only Lipschitz regularity of the interface.

Lemma 3.2** (Finite element extension).**

Assume $\Gamma$ is Lipschitz, the meshes $\{\mathcal{T}_{h}\}$ are shape regular and satisfy Assumptions 2, and assume $\mathbf{V}_{h}^{\rm bulk}$ has a local nodal basis. There exists an extension operator $E_{h}:\text{\rm V}_{h}^{+}\to\text{\rm V}_{h}^{\rm bulk}$ with the following properties:

• a)

$E_{h}\text{\rm v}_{h}=\text{\rm v}_{h}$* on $\Omega_{h}^{+}$,*

• b)

There exists $C>0$, independent of $h$ and position of $\Gamma$ against the underlying mesh, such that for all $\text{\rm v}_{h}\in\text{\rm V}_{h}^{+}$,

[TABLE]

Proof.

The proof is given in the Appendix. ∎

3.3. The discrete variational formulation

In this section, we define the discrete counterpart of (2.1). The jumps over the interface are enforced weakly, and a term is added to enforce the symmetry of the bilinear form $a$. A discrete variational analogue of (2.1) is given by the problem of finding $(\text{\rm u}_{h},p_{h})\in\text{\rm V}_{h}\times M_{h}$ such that

[TABLE]

with bilinear forms given below. For all $\text{\rm u}_{h}=(\text{\rm u}_{h}^{-},\text{\rm u}_{h}^{+}),$ $\text{\rm v}_{h}=(\text{\rm v}_{h}^{-},\text{\rm v}_{h}^{+})\in\text{\rm V}_{h}$ and $p_{h}=(p_{h}^{-},p_{h}^{+})$, $q_{h}=(q_{h}^{-},q_{h}^{+})\in M_{h}$, we define

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]

Stabilization parameters $\gamma^{\pm}_{p}$, $\gamma^{\pm}_{\text{\rm u}}$ and $\gamma$ are all assumed to be independent of $\nu$, $h$ and position of $\Gamma$ against the underlying mesh. Parameter $\gamma$ needs to be large enough to provide the bilinear form $a(\cdot,\cdot)$ with coercivity. For the purpose of analysis, we set $\gamma^{\pm}_{p}=\gamma^{\pm}_{\text{\rm u}}=1$. In practice, these parameters can be tuned for better numerical performance (see section 5 for numerical examples) and the analysis below remains valid if all $\gamma^{\pm}_{p}$ and $\gamma^{\pm}_{\text{\rm u}}$ are $O(1)$ parameters.

The right-hand side $\text{\rm F}_{h}\in\text{\rm V}_{h}^{*}$ will be defined later on, and we assume $G_{h}\in M_{h}^{*}$ is given by

[TABLE]

It is straightforward to check that the norm of linear bounded functions $G_{h}^{\pm}$ can be expressed in terms of $M_{h}^{\pm}$ spaces. More precisely, it holds

[TABLE]

Also, we assume $G_{h}$ satisfies $G_{h}(1)=0$. In particular, this implies that if the second equation of (3.7) is satisfied by $q_{h}\in M_{h}$, then it is also satisfied by $q_{h}+c$ for any constant function $c$.

3.4. Well-posedness of the discrete scheme

We are now interested in a finite element counterpart of the a priori estimate (2.3). More precisely, for the solution $(\text{\rm u}_{h},p_{h})$ of (3.7), we shall prove the following stability result:

[TABLE]

with a constant $C>0$ independent of $\nu^{\pm}$, $h$ and the position of $\Gamma$ in the bulk mesh. The proof will largely follow the main steps made in section 2. We first prove specific estimates for discrete pressure in $\Omega$ and subdomains. The continuity of the bilinear forms $a_{h}$ and $b_{h}$, and the coercivity of $a_{h}$ in $\text{\rm V}_{h}$ will help us with energy estimates, which due to $\nu^{-}\leq\nu^{+}$ yield the desired control of $\nu D(\text{\rm u}_{h})$, pressure and stabilization terms in $\Omega^{-}$. To extend these estimates to $\Omega^{+}$, a crucial result about extension of finite element functions is stated (with its proof moved to the Appendix section). The result is then applied to gain control of finite element viscous stresses and pressure in $\Omega^{+}$. Define the natural energy norm for $a_{h}$ by

[TABLE]

We will need the following technical lemma which is essentially found in [12, Lemma 5.1]. The main difference is the use of Korn’s inequality and projections onto the rigid body motions instead of the Poincaré inequality and projection onto the constants.

Lemma 3.3**.**

There exists $C>0$, independent of $h$ and $\nu^{\pm}$, such that for every $q_{h}\in M_{h}$, it holds

[TABLE]

and for all $\text{\rm v}_{h}\in\text{\rm V}_{h}$,

[TABLE]

The following result discusses the continuity and coercivity of the bilinear form $a_{h}$. We omit the proof since they are by now standard and can easily be proved using Lemma 3.1. See, for example, similar results in [9].

Lemma 3.4**.**

Let $\text{\rm v}_{h},\,\text{\rm w}_{h}\in\text{\rm V}_{h}$ and $q_{h}\in M_{h}$. Then, there exists $C>0$ and $h_{0}$, independent of $h$ and $\nu^{\pm}$, $\text{\rm v}_{h}$ and $\text{\rm w}_{h}$, such that for all $h\leq h_{0}$ it holds

[TABLE]

[TABLE]

Additionally, for all $\text{\rm v}\in[\text{\rm H}^{s+1}(\Omega_{h}^{-})\times\text{\rm H}^{s+1}(\Omega_{h}^{+})]\times\text{\rm V}_{h}$ we have

[TABLE]

Finally, there exists $\alpha>0$, independent of $h$ and $\nu^{\pm}$, such that

[TABLE]

The following result is the discrete analogue of Lemma 2.1.

Lemma 3.5**.**

Let $\text{\rm u}_{h}\in\text{\rm V}_{h}$ and $p_{h}\in M_{h}$ solve (3.7). Then, there exists $C,h_{0}>0$, depending only on $\Omega$ and $\Gamma$, such that for $h\leq h_{0}$ it holds

[TABLE]

and

[TABLE]

Proof.

Let $q_{h}^{\pm}:=p_{h}^{\pm}-\mathrm{avg}_{\Omega_{h,i}^{\pm}}(p_{h}^{\pm})$. Noting that $J_{h}^{\pm}(p_{h}^{\pm},p_{h}^{\pm})=J_{h}(q_{h}^{\pm},q_{h}^{\pm})$ and employing (3.10) from Lemma 3.3, we get

[TABLE]

On the other hand, since $q_{h}^{\pm}\in L_{0}^{2}(\Omega_{h,i}^{\pm})$, Theorem 3.1 implies that there exist $\theta,h_{0}>0$, depending only on $\Omega$, such that for all $h\leq h_{0}$, it holds

[TABLE]

However, for $\text{\rm w}_{h}^{\pm}\in{\text{\rm V}_{h}^{\rm bulk}}\cap\text{\rm H}_{0}^{1}(\Omega_{h,i}^{\pm})$ we notice equalities

[TABLE]

where $\text{\rm v}_{h}=(\text{\rm w}_{h}^{-},0)$ or $\text{\rm v}_{h}=(0,\text{\rm w}_{h}^{+})$. Because $\text{\rm v}_{h}$ is supported on $\Omega_{h,i}^{\pm}$ we have

[TABLE]

With the help of the Cauchy-Schwarz inequality, inverse estimates and Korn’s inequality we obtain

[TABLE]

Using (3.18)–(3.19) in (3.17), we arrive at

[TABLE]

which leads to (3.14).

In order to prove (3.15), we consider $\alpha_{h}=p_{h}-q_{h}$ and observe that $\|\alpha_{h}\|_{L^{\infty}(\Omega)}\leq C\|p_{h}\|_{\Omega}$. Moreover, a simple calculation shows that

[TABLE]

Hence, using that $(p_{h}^{-},1)_{\Omega^{-}}+(p_{h}^{+},1)_{\Omega^{+}}=0$ and that ${|\Omega_{h}^{\pm}\setminus\Omega_{h,i}^{\pm}|}=|\Omega_{h}^{\pm}|-|\Omega_{h,i}^{\pm}|\leq C\,h$ we have

[TABLE]

In the case LABEL:when ${M_{h}^{\rm bulk}}=M_{h,\rm disc}$ we let $r_{h}\in{M_{h}^{\rm bulk}}$ to be the $L^{2}$ projection of $\alpha_{h}$ onto piecewise constants with respect to the mesh $\mathcal{T}_{h}$. In the case, ${M_{h}^{\rm bulk}}=M_{h,\rm cont}$ we let $r_{h}\in{M_{h}^{\rm bulk}}$ be the continuous piecewise linear function such that $r_{h}(x)=\alpha_{h}^{+}(x)$ if $x$ is a vertex and $x\in\Omega^{+}$, and $r_{h}(x)=\alpha_{h}^{-}(x)$ if $x$ is a vertex and $x\in\overline{\Omega^{-}}$. In either case, $r_{h}|_{\Omega_{h,i}^{\pm}}=\alpha_{h}^{\pm}$.

Recalling the notation $\omega_{h}=\Omega\setminus(\Omega_{h,i}^{-}\cup\Omega_{h,i}^{+})$, we then note that

[TABLE]

We let $\tilde{r}_{h}=r_{h}-\text{avg}_{\Omega}(r_{h})$. From (3.22) and (3.23) we have that

[TABLE]

Assumption 1 provides us with $\text{\rm v}_{h}\in{\text{\rm V}_{h}^{\rm bulk}}$ such that

[TABLE]

Let $\text{\rm w}_{h}\in\text{\rm V}_{h}$ be given by $\text{\rm w}_{h}=(\text{\rm v}_{h}|_{\Omega_{h}^{-}},\text{\rm v}_{h}|_{\Omega_{h}^{+}})$. It holds $\|\text{\rm w}_{h}\|_{1,h}\leq C\|\text{\rm v}_{h}\|_{1,\Omega}$. To verify the last inequality, we note that the first term in the definition of $\|\text{\rm w}_{h}\|_{1,h}$ vanish, while the second jump term can be estimated with the help of the finite element trace and inverse inequalities. Hence, we get

[TABLE]

where we used the first equation of (3.7). Thanks to (3.12) and (3.25) we have

[TABLE]

We then use the triangle inequality, (3.21), (3.24), and (3.26) to get

[TABLE]

The result now follows after taking $h_{0}$ small enough.

∎

In order to prove (3.8), we start by obtaining an estimate for $\|\nu^{-}D(\text{\rm u}_{h}^{-})\|_{\Omega^{-}}$ using the energy norm (3.9).

Lemma 3.6**.**

Let $(\text{\rm u}_{h},p_{h})\in V_{h}\times M_{h}$ be a solution of (3.7). Then, there exists $C>0$, independent of $h$ and $\nu^{\pm}$, such that

[TABLE]

Proof.

We use the first and second equation of (3.7) with $\text{\rm v}_{h}=\text{\rm u}_{h}$ and $q_{h}=p_{h}$, respectively, and the coercivity of $a_{h}$ for $\gamma$ large enough to get

[TABLE]

with some $\alpha>0$ independent of $\nu$ and $h$. By definition of the norms $\|\cdot\|_{{1,h}}$ and $\|\cdot\|_{\text{\rm V}_{h}}$, and since $\nu^{-}\leq\nu^{+}$, we get

[TABLE]

so that it only remains to estimate $|G_{h}(p_{h})|$. In order to do this, we will decompose the expression $G_{h}(p_{h})=G_{h}^{-}(p_{h}^{-})+G_{h}^{+}(p_{h}^{+})$ into three terms:

[TABLE]

where $\alpha_{h}^{\pm}:=\mathrm{avg}_{\Omega_{h,i}^{\pm}}(p_{h}^{\pm})$. Then, using (3.14) and $\nu^{-}\leq\nu^{+}$ we get

[TABLE]

On the other hand,

[TABLE]

where we used our assumption that $G_{h}^{+}(1)=-G_{h}^{-}(1)$ . Therefore, we have

[TABLE]

where we used that $\nu^{-}\leq\nu^{+}$. After re-scaling by $\frac{1}{\sqrt{\nu^{+}}}$ the estimate (3.15) yields for the last term on the right-hand side of (3.30) the bound,

[TABLE]

again due to $\nu^{-}\leq\nu^{+}$. Using now (3.29) and (3.30) to bound from above the right-hand side of (3.28) leads after some calculations to

[TABLE]

Therefore, we have

[TABLE]

The result now follows by multiplying both sides by $\nu^{-}$ and using that $\nu^{-}\leq\nu^{+}$. ∎

As an immediate consequence of Lemma 3.6 and (3.14), we have the following result.

Lemma 3.7**.**

Let $\text{\rm u}_{h}\in\text{\rm V}_{h}$ and $p_{h}\in M_{h}$ solve (3.7). Then, there exist $C,h_{0}>0$, depending only on $\Omega$ and $\Gamma$, such that for $h\leq h_{0}$,

[TABLE]

The next result is the discrete analogue of Lemma 2.4.

Lemma 3.8**.**

Let $\text{\rm u}_{h}\in\text{\rm V}_{h}$ and $p_{h}\in M_{h}$ solve (3.7). Then, there exist $C,h_{0}>0$, depending only on $\Omega$ and $\Gamma$, such that for $h\leq h_{0}$ it holds

[TABLE]

Proof.

We first note that the rigid body motions belong to the velocity finite element space , and using Lemma 3.2, we consider the discrete extension $E_{h}:\text{\rm V}_{h}^{+}\to\text{\rm V}_{h}^{\rm bulk}$ and define $\text{\rm w}_{h}\in\text{\rm V}_{h}^{\rm bulk}$ by

[TABLE]

Let $h_{0}>0$ be sufficiently small such that if $\Omega^{+}$ is the inclusion, then $\partial\Omega\cap\partial\Omega_{h}^{+}=\emptyset$. Then, if $\Omega^{+}$ is not the inclusion, we have by the boundedness of the extension $E_{h}$ and by Korn’s inequality that $\|\text{\rm w}_{h}\|_{1,\Omega}\leq\|\text{\rm u}_{h}^{+}\|_{1,\Omega_{h}^{+}}\leq C\|D(\text{\rm u}_{h}^{+})\|_{\Omega_{h}^{+}}$. On the other hand, if $\Omega^{+}$ is the inclusion, we have by the boundedness of the discrete extension $E_{h}$, Korn’s inequality and since

[TABLE]

that $\|\text{\rm w}_{h}\|_{1,\Omega}\leq C\|\text{\rm u}_{h}^{+}-\mathcal{P}_{\Omega_{h}^{+}}^{\mathrm{RM}}(\text{\rm u}_{h}^{+})\|_{1,\Omega_{h}^{+}}\leq C\|D(\text{\rm u}_{h}^{+})\|_{\Omega_{h}^{+}}$. In both cases, it holds

[TABLE]

We let $\text{\rm v}_{h}\in\text{\rm V}_{h}$ be given by $\text{\rm v}_{h}=(\text{\rm w}_{h}|_{\Omega_{h}^{-}},\text{\rm w}_{h}|_{\Omega_{h}^{+}})$. By the same arguments as in the proof of Lemma 3.5, one shows $\|\text{\rm v}_{h}\|_{1,h}\leq C\|\text{\rm w}_{h}\|_{1,\Omega}$. First equation from (3.7) yields

[TABLE]

where

[TABLE]

We bound every term on the right hand side of (3.35). With the help of (3.6) and (3.11), we get for the last term

[TABLE]

where

[TABLE]

consists of terms already estimated in Lemma 3.6. For the first term on the right-hand side of (3.35) we have

[TABLE]

It remains to estimate the second term on the right-hand side of (3.35). Noting that (3.33) implies $\text{\rm div}\,\text{\rm w}_{h}^{+}=\text{\rm div}\,\text{\rm u}_{h}^{+}$ and using $\Omega^{+}\subset\Omega^{+}_{h}$, we get

[TABLE]

Setting $\alpha_{h}^{-}=\text{avg}_{\Omega_{h,i}^{-}}(p_{h}^{-})$ and using that $(\text{\rm div}\,\text{\rm w}_{h},1)_{\Omega}=0$, we have

[TABLE]

We also have from the second equation of (3.7) with $q_{h}=(0,p_{h}^{+})$

[TABLE]

Now we see that

[TABLE]

We combine (3.38), (3.37), (3.36), and (3.35) to get

[TABLE]

where

[TABLE]

It remains to estimate the solution-dependent terms on the right-hand side of (3.39). Thanks to (3.15) we get

[TABLE]

From (3.34) we conclude that $\nu^{+}\|\text{\rm w}_{h}\|_{1,\Omega}\leq C\nu^{+}\|D(\text{\rm u}_{h}^{+})\|_{\Omega^{+}}\leq{C}\sqrt{L}$ and hence we have

[TABLE]

The result now follows after using (3.31) and (3.27). ∎

As an immediate consequence of Lemma 3.8 and (3.15), we have the following result.

Theorem 3.2**.**

Let $(\text{\rm u}_{h},p_{h})\in\text{\rm V}_{h}\times M_{h}$ be a solution of the finite element method (3.7). Then, there exists $C,h_{0}>0$, independent of $h$, position of $\Gamma$ in the mesh and $\nu^{\pm}$, such that for $h\leq h_{0}$, it holds

[TABLE]

4. Finite element error estimates

We use the stability results from the previous sections to obtain error estimates for both $\|\nu D(\text{\rm u}-\text{\rm u}_{h})\|_{\Omega}$ and $\|p-p_{h}\|_{\Omega}$. We assume that the solutions to the two-phase Stokes problem (1.1) sufficiently smooth in each subdomain. In particular, we assume that $\text{\rm u}^{\pm}\in\text{\rm H}^{s+1}(\Omega^{\pm})$ and $p^{\pm}\in\mathrm{H}^{k_{p}+1}(\Omega^{\pm})$. In such a case there exist extensions of u and $p$, from $\Omega^{\pm}$ to $\Omega^{\mp}$, denoted by $\text{\rm u}_{E}^{\pm}$ and $p_{E}^{\pm}$ and with the property $\text{\rm u}_{E}^{\pm}|_{\Omega^{\pm}}=\text{\rm u}^{\pm}$ and $p_{E}^{\pm}|_{\Omega^{\pm}}=p^{\pm}$, such that $\text{\rm u}_{E}^{\pm}\in\text{\rm H}^{s+1}(\Omega)$, $p_{E}^{\pm}\in\mathrm{H}^{k_{p}+1}(\Omega)$,

[TABLE]

where $C$ depends only on $\Omega$ and $\Gamma$. Further, we will identify $\text{\rm u}^{\pm}$ and $p^{\pm}$ with there extensions.

For the error analysis in addition to the augmented norm for the velocity (3.2) we also define

[TABLE]

Multiplying the first equation of (1.1) with $\text{\rm v}^{\pm}$, using that $\llbracket{\sigma(\text{\rm u},p)\text{\rm n}}\rrbracket=\bm{\lambda}$, and taking into account the choice of the weights for the average $\{\text{\rm v}\}_{\Gamma}$ in the definition of $a_{h}$, we define $\text{\rm F}_{h}$ by

[TABLE]

Then, we have the main result of this section.

Theorem 4.1**.**

Let $(\text{\rm u},p)\in\text{\rm V}\times M$ be a solution of (2.1) with $\text{\rm f}\in\text{\rm L}^{2}(\Omega)$, $\bm{\lambda}\in\text{\rm L}^{2}(\Gamma)$ and $g\in L_{0}^{2}(\Omega)$. Furthermore, let $(\text{\rm u}_{h},p_{h})$ be the approximation that solves (3.7) with $\text{\rm F}_{h}$ given by (4.1) and $G_{h}^{\pm}(q_{h}^{\pm})=-(g,q_{h}^{\pm})_{\Omega^{\pm}}$. Assume that $\text{\rm u}^{\pm}\in\text{\rm H}^{s+1}(\Omega^{\pm})$ and $p^{\pm}\in H^{k_{p}+1}(\Omega^{\pm})$. Then there exists a constant $C$ independent of $\nu$, $h$, u, $p$ such that

[TABLE]

Proof.

As discussed above we let $\text{\rm u}^{\pm}$ and $p^{\pm}$ be the extensions of the same functions to the entire $\Omega$. Then, by using (4.1), we one easily checks the consistency result:

[TABLE]

Hence, we have for an arbitrary $\text{\rm w}_{h}\in\text{\rm V}_{h}$ and $r_{h}\in Q_{h}$ the following consistency result:

[TABLE]

where

[TABLE]

and

[TABLE]

We can easily show, using for example (3.13), the following bound

[TABLE]

which implies

[TABLE]

Similarly,

[TABLE]

Hence,

[TABLE]

The result follows after applying Theorem 3.2 and triangle inequality. ∎

Using (3.3) and (3.4) with standard interpolation properties of finite element functions and the definition of the norms of the right-hand side of (4.2) the next results follows from the theorem.

Corollary 4.1**.**

Under the same assumptions as in Theorem 4.1, it holds:

[TABLE]

with a constant $C$ independent of $\nu$, $h$ and the position of the interface in the background mesh. The solution norms on the right-hand side of (4.5) are the norms in the broken Sobolev spaces ${H}^{\ell}(\Omega^{-})\times{H}^{\ell}(\Omega^{+})$,

[TABLE]

and similar for vector functions from ${H}^{\ell}(\Omega^{-})^{d}\times{H}^{\ell}(\Omega^{+})^{d}$.

5. Numerical experiments

Example 1

Consider the squared domain $\Omega:=(-1,1)\times(-1,1)$ and the embedded interface $\Gamma:x_{1}^{2}+x_{2}^{2}=R^{2}$ for $R=1/\sqrt{\pi}$. We define $\Omega^{-}=\{\text{\rm x}\in\mathbb{R}^{2}\,:\,|\text{\rm x}|<R\}$, and $\Omega^{+}=\Omega\setminus\overline{\Omega}^{+}$, and choose the data $\text{\rm f}\in\mathbf{L}^{2}(\Omega)$ so that the exact solution $(\text{\rm u},p)$ is given for all $\text{\rm x}=(x_{1},x_{2})$ by

[TABLE]

We observe that $\llbracket{\text{\rm u}}\rrbracket=\text{\rm 0}$ and $\llbracket{\nu D(\text{\rm u})}\rrbracket=0$, and since $p$ is continuous, we get $\llbracket{\sigma(\text{\rm u},p)\text{\rm n}}\rrbracket=\text{\rm 0}$. Also, u is divergence free. We will show that the errors $\|\nu\,D(\text{\rm u}-\text{\rm u}_{h})\|_{\Omega}$ and $\|p-p_{h}\|_{\Omega}$ are independent of $\nu$. In order to do this, we consider a uniform diagonal triangular decomposition of $\Omega$, and test the code with two pairs of spaces $\text{\rm V}_{h}^{\rm bulk}\times M_{h}^{\rm bulk}$ that satisfy Assumption 4. These are the Mini-element and the pair $P_{2}$ – $P_{0}$. We choose a mesh with $N=160092$ degrees of freedom for the Mini-element and $N=164697$ degrees of freedom for the pair $P_{2}$ – $P_{0}$.

We denote the errors by ${\tt e}(\text{\rm u}):=\|\nu D(\text{\rm u}-\text{\rm u}_{h})\|_{\Omega}$ and ${\tt e}(p):=\|p-p_{h}\|_{\Omega}$ and compute them experimentally, with decreasing values of $\nu^{-}$ and increasing values of $\nu^{+}$. As a good balance between stability and conditioning, we set the stabilization parameters to $\gamma=25$, $\gamma_{\text{\rm u}}^{\pm}=15$ and $\gamma_{p}^{\pm}=20$ for the Mini-element and $\gamma=20$, $\gamma_{\text{\rm u}}^{\pm}=10$ and $\gamma_{p}^{\pm}=15$ for the pair $P_{2}$ – $P_{0}$. The numerical results are summarized in Table 5.1.

From Table 5.1, we observe that the errors ${\tt e}(\text{\rm u})$ and ${\tt e}(p)$ remain unchanged for a fixed mesh when $\nu^{-}$ decreases and $\nu^{+}$ increases.

Now we switch the role of $\Omega^{+}$ and $\Omega^{-}$, and define $\Omega^{+}=\{\text{\rm x}\in\mathbb{R}^{2}\,:\,|\text{\rm x}|<R\}$, $\Omega^{-}=\Omega\setminus\overline{\Omega}^{+}$. Observe that $\Omega^{+}$ is now the inclusion. We choose the data $\text{\rm f}\in\text{\rm L}^{2}(\Omega)$ so that the exact solution $(\text{\rm u},p)$ is given by

[TABLE]

We choose a mesh with $N=160092$ degrees of freedom for the Mini-element and $N=164697$ degrees of freedom for the pair $P_{2}$ – $P_{0}$, and as a good balance between stability and conditioning we set the stabilization parameters to $\gamma=25$, $\gamma_{\text{\rm u}}^{\pm}=15$ and $\gamma_{p}^{\pm}=20$ for the Mini-element, and $\gamma=20$, $\gamma_{\text{\rm u}}^{\pm}=10$ and $\gamma_{p}^{\pm}=15$ for the pair $P_{2}$ – $P_{0}$, consider decreasing values of $\nu^{-}$ and increasing values of $\nu^{+}$, and summarize the results in Table 5.2.

Similarly, we observe that the errors ${\tt e}(\text{\rm u})$ and ${\tt e}(p)$ remain unchanged for a fixed mesh when $\nu^{-}$ decreases and $\nu^{+}$ increases.

Example 2

We consider the same exact solution $(\text{\rm u},p)$ given by (5.1), and the finite element errors

[TABLE]

We test the method for $P_{2}$ – $P_{0}$ bulk spaces and fixed viscosity and stabilization parameters $\nu^{-}=0.5$, $\nu^{+}=20$, $\gamma=20$, $\gamma_{\text{\rm u}}^{\pm}=10$ and $\gamma_{p}^{\pm}=15$. We consider a sequence of uniform triangular meshes with decreasing mesh size. The experimental rates of convergence are computed as

[TABLE]

where $\Phi_{j}$ is the corresponding error norm at mesh level $j$. The error norms and experimental rates are shown in Table 5.3. Also, we show plots of the approximate solution $(\text{\rm u}_{h},p_{h})$ in Figure 5.1.

We repeat the experiment using the Mini-element and the same set of parameters and meshes. The computed error norms and experimental rates are shown in Table 5.4.

Example 3

We consider the parameters $\nu^{-}=0.5$, $\nu^{+}=2$, and the exact solution $(\text{\rm u},p)$ given by

[TABLE]

We observe that $\llbracket{\text{\rm u}}\rrbracket=\text{\rm 0}$, and the jump $\llbracket{\sigma(\text{\rm u},p)\text{\rm n}}\rrbracket$ is non-zero and is given by

[TABLE]

Additionally, u is divergence-free and $(p,1)_{\Omega}=0$. We test the code with the stabilization parameters $\gamma=30$, $\gamma_{\text{\rm u}}^{\pm}=25$ and $\gamma_{p}^{\pm}=25$ for the Mini-element, $\gamma=30$, $\gamma_{\text{\rm u}}^{\pm}=25$ and $\gamma_{p}^{\pm}=20$ for the pair $P_{2}$ – $P_{0}$, decompose the domain by a sequence of uniform meshes with decreasing size, and summarize the results in the Tables 5.5 and 5.6.

We finish by showing plots of the computed finite element solution $(\text{\rm u}_{h},p_{h})$ in Figure 5.2.

Example 4

We consider the same exact solution $(\text{\rm u},p)$ given by (5.2) and repeat experiments from Example 3, this time with experimental errors for $\|\text{\rm div}(\text{\rm u}-\text{\rm u}_{h})\|_{\Omega}$, and also with respect to the $L^{\infty}$ norm for the (viscous part of the) stress tensor and pressure

[TABLE]

respectively. We test the method for $P_{2}$ – $P_{0}$ bulk spaces and fixed viscosity and stabilization parameters $\nu^{-}=0.5$, $\nu^{+}=20$, $\gamma=20$, $\gamma_{\text{\rm u}}^{\pm}=10$ and $\gamma_{p}^{\pm}=15$. We build a sequence of uniform triangular meshes with decreasing mesh size. The error norms and experimental rates are shown in Table 5.7.

We repeat the experiment using the Mini-element, the same viscosities $\nu^{\pm}$ and stabilization parameters $\gamma=30$, $\gamma_{\text{\rm u}}^{\pm}=15$ and $\gamma_{p}^{\pm}=5$. The computed error norms and experimental rates are shown in Table 5.8. Although, the error estimates in $L^{\infty}$ norm are not covered by the analysis of the paper, the numerical experiments demonstrate the first order of convergence.

Appendix A Finite element extension.

The goal of this section is to provide the $H^{1}$-bounded extension operator for finite elements of arbitrary degree, i.e. to prove Lemma 3.2. We build the desired extension for each velocity component independently. Let $V_{h}^{k}$ be the $H^{1}_{0}(\Omega)$-conformal FE space of degree $k$. Then $\mathbf{V}_{h}^{\rm bulk}=\otimes_{j=1}^{d}V_{j,h}$, and each $V_{j,h}$ satisfies $V_{h}^{1}\subset V_{j,h}\subset V_{h}^{r}$ for some integer $r>0$. Fix arbitrary $j\in\{1,\dots,d\}$. Using the definition of the local nodal basis, mapping $T\in\mathcal{T}_{h}$ to the reference simplex and the equivalence of norms in a space of finite dimension, one shows

[TABLE]

with some $C$ independent of $v$ and $T$.

The strategy will be to build an extension operator for piecewise linears first. Then to use that extension operator to build a general extension operator.

A.1. Extension operator for piecewise linears

Let $V_{h}^{1,+}=\{v|_{\Omega_{h}^{+}}:v\in V_{h}^{1}\}$. Consider $v\in V_{h}^{1,+}$ and let $Ev\in H_{0}^{1}(\Omega)$ be the Stein extension [16] such that

[TABLE]

We need $I^{\rm SZ}(Ev)$, the Scott-Zhang interpolant of $Ev$ onto $V_{h}^{1}$. The construction of $I^{\rm SZ}(Ev)\in V_{h}^{1}$ follows the standard procedure from [15]. However, some care is required to ensure that we are recovering the same $P^{1}$ function in all interior tetrahedra of $\Omega^{+}$,

[TABLE]

To provide (A.3), we exploit a freedom in choosing the Scott-Zhang interpolant pointed out in [15]: For every vertex y of $\mathcal{T}_{h}$ we need to associate either a $d$-dimensional simplex to y or a $d-1$-dimensional simplex. If y is a vertex for some $T\in\mathcal{T}_{h,i}^{+}$ and $\text{\rm y}\notin\partial\Omega$, then we associate one of these simplices from $\mathcal{T}_{h,i}^{+}$ with y. By the stability property of the Scott-Zhang interpolant we have

[TABLE]

Now we define a discrete extension operator for piecewise linear $v\in V_{h}^{1,+}$:

[TABLE]

Note that

[TABLE]

We decompose $\mathcal{T}_{h,i}^{-}=\mathcal{T}_{h,i}^{-,1}\cup\mathcal{T}_{h,i}^{-,\text{int}}$ where $\mathcal{T}_{h,i}^{-,1}=\{T\in\mathcal{T}_{h,i}^{-}:T\text{ has a vertex }\text{\rm y}\text{ such that }\text{\rm y}\in K,\text{ for some }K\in\mathcal{T}_{h}^{\Gamma}\}$ and $\mathcal{T}_{h,i}^{-,\text{int}}=\mathcal{T}_{h,i}^{-}\backslash\mathcal{T}_{h,i}^{-,1}$.

We then see that

[TABLE]

Therefore, we are left to bound $\sum_{T\in\mathcal{T}_{h,i}^{-,1}}\|\nabla E_{h}^{1}v\|_{T}^{2}$. We use the triangle inequality to get

[TABLE]

For ease of notation we set $w=E_{h}^{1}v-I^{\text{SZ}}Ev$. Now note that if $T\in T\in\mathcal{T}_{h,i}^{-,1}$ then $w$ vanishes on all vertices that do not belong to $\overline{\Omega_{h}^{+}}$. Hence, using inverse estimates we get

[TABLE]

Hence, we will have

[TABLE]

Recalling Assumption 2, if $T\in\mathcal{T}_{h}^{\Gamma}$ there exists $K_{T}\in\mathcal{T}_{h,i}^{+}$ such that $T=K_{1},K_{2},\ldots,K_{\ell}=K_{T}$ and $K_{j},K_{j+1}$ have a common $d-1$ face for $j<\ell$ and $K_{j}\subset\Omega_{h}^{+}$. The number $\ell\leq M$ where $M$ is uniformly bounded and only depends on the shape regularity of the mesh. Then we see from a simple scaling argument that

[TABLE]

Using (A.3) we have $w|_{K_{T}}\equiv 0$ since $K_{T}\in\mathcal{T}_{h,i}^{+}$ and so

[TABLE]

In the last inequality we used that by shape regularity $h_{K_{i}}\leq Ch_{T}$ where $C$ depends on $M$ and shape regularity constant. We then get

[TABLE]

In the last step we used (A.1). Therefore, combining all the inequalities above we obtain

[TABLE]

We hence, get that after using (A.2) and (A.4) that

[TABLE]

Finally, since $E_{h}^{1}v|_{\Omega_{h}^{+}}=v$ and $\overline{\Omega}=\overline{\Omega_{h,i}^{-}}\cup\overline{\Omega_{h}^{+}}$ as indicated in section 3.1, we get

[TABLE]

A.2. General Discrete Extension Operator

Building on the availability of $E_{h}^{1}$ we define the general extension operator. Let $V_{h}^{+}=\{v|_{\Omega_{h}^{+}}:v\in V_{j,h}\}$ and consider the subspace $W_{h}^{+}=\{v\in V_{h}^{+}:v(\text{\rm y})=0\text{ for all vertices }\text{\rm y}\text{ of }\mathcal{T}_{h}^{+}\}$. For $v\in W_{h}^{+}$ we consider the extension $Q_{h}v\in V_{j,h}$ by defining its nodal values as follows

[TABLE]

Then we can easily prove the following lemma.

Lemma A.1**.**

For $v\in W_{h}^{+}$ it holds

[TABLE]

where the constant $C_{1}$ only depends on the shape regularity of the mesh.

Proof.

Note that $Q_{h}v|_{\Omega_{h,i}^{-}}$ is supported on all simplices $T\in\mathcal{T}_{h}^{-}$ that have at least one edge belonging to $\partial\Omega_{h}^{+}$. Lets call this set of simplices $\omega_{h}^{-}$:

[TABLE]

and we also consider the set

[TABLE]

For each $T\in\omega_{h}^{-}$, let $\Xi(T):=\{\tau\in\mathcal{T}_{h}^{+}:\tau\text{ and }T\text{ share a common edge}\}$. Then we have due to the finite element inverse estimates:

[TABLE]

At the same time, with the help of (A.1) we have for each $T\in\omega_{h}^{+}$,

[TABLE]

For the second inequality we used that $Q_{h}v$ vanishes on all nodes except the nodal points of $T$ that belong $\partial T$ and belong to the boundary of $\partial\Omega_{h}^{+}$. Since $v\in W_{h}^{+}$ and so $v$ vanishes on all the vertices of such $\tau^{\prime}s$, we obtain

[TABLE]

Finally, applying inverse estimates gives

[TABLE]

Hence,

[TABLE]

∎

Let $I_{h}:C(\Omega_{h}^{+})\rightarrow V_{h}^{1,+}$ be the Lagrange interpolant. In section A.1 we defined a stable discrete extension operator $E_{h}^{1}:V_{h}^{1,+}(\Omega_{h}^{+})\to V_{h}^{1}\subset V_{j,h}$.

We finally can define our general discrete extension operator. For any $v\in V_{j,h}^{+}$, we define $E_{h}v\in V_{j,h}$ as follows

[TABLE]

Note that $v-I_{h}v\in W_{h}^{+}$ so indeed this definition makes sense.

We use the boundedness of $E_{h}^{1}$, (A.5), (A.6) and the stability of $I_{h}$ to obtain:

[TABLE]

with some $C$ independent of $h$ and the position of $\Gamma$ in the background mesh. Hence, we we have proven Lemma 3.2.