A Nitsche cut finite element method for the Oseen problem with general Navier boundary conditions

TL;DR

This paper introduces a Nitsche-based cut finite element method for the Oseen problem that robustly handles generalized Navier boundary conditions, including zero slip length, ensuring stability and optimal error estimates across all flow regimes.

Contribution

A novel Nitsche-based weak imposition method for generalized Navier conditions that remains stable and well-posed for all slip lengths, including zero, improving upon existing approaches.

Findings

Method is stable for all slip lengths including zero.

Numerical tests confirm robustness and stability across flow regimes.

Achieves L2-optimal velocity error estimates.

Abstract

In this work a Nitsche-based imposition of generalized Navier conditions on cut meshes for the Oseen problem is presented. Other methods from literature dealing with the generalized Navier condition impose this condition by means of substituting the tangential Robin condition in a classical Galerkin way. These methods work fine for a large slip length coefficient but lead to conditioning and stability issues when it approaches zero. We introduce a novel method for the weak imposition of the generalized Navier condition which remains well-posed and stable for arbitrary choice of slip length, including zero. The method proposed here builds on the formulation done by [1]. They impose a Robin condition for the Poisson problem by means of Nitsche's method for an arbitrary combination of the Dirichlet and Neumann parts of the condition. The analysis conducted for the proposed method is done in a similar fashion as in [2], but is done here for a more general type of boundary condition. The analysis proves stability for all flow regimes and all choices of slip lengths. Also an L2-optimal estimate for the velocity error is shown, which was not conducted in the previously mentioned work. A numerical example is carried out for varying slip lengths to verify the robustness and stability of the method with respect to the choice of slip length. Even though proofs and formulations are presented for the more general case of an unfitted grid method, they can easily be reduced to the simpler case of a boundary-fitted grid with the removal of the ghost-penalty stabilization terms.

Equations

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}-\nabla \cdot (2\nu \mathbf{ϵ}(\mathbf{u}))+\nabla p$$

$$\nabla \cdot \mathbf{u}$$

$$\left(\mathbf{u}-\mathbf{g}\right){\mathbf{P}}^n$$

$$\left(\epsilon (2\nu \mathbf{ϵ}(\mathbf{u})\mathbf{n}-\mathbf{h})+\nu (\mathbf{u}-\mathbf{g})\right){\mathbf{P}}^t$$

$$2\nu \mathbf{ϵ}(\mathbf{u})\mathbf{n}{\mathbf{P}}^t=\mathbf{h}{\mathbf{P}}^t-\frac{\nu }{\epsilon }\left(\mathbf{u}-\mathbf{g}\right){\mathbf{P}}^t\text{ on }\partial \Omega .$$

$$\sigma \mathbf{u}+\mathbf{\beta }\cdot \nabla \mathbf{u}-\nabla \cdot (2\nu \mathbf{ϵ}(\mathbf{u}))+\nabla p$$

$$\nabla \cdot \mathbf{u}$$

$$\left(\mathbf{u}-\mathbf{g}\right){\mathbf{P}}^n$$

$$\left(\epsilon (2\nu \mathbf{ϵ}(\mathbf{u})\mathbf{n}-\mathbf{h})+\nu (\mathbf{u}-\mathbf{g})\right){\mathbf{P}}^t$$

$$\mathbf{g}\in [H^{3/2}(\Gamma ){]}^d\text{ and }\mathbf{h}\in [H^{1/2}(\Gamma ){]}^d,$$

$${\mathcal{V}}_{\mathbf{g}}$$

$${\mathcal{V}}_0^n$$

$$\mathcal{Q}$$

$$a(\mathbf{u},\mathbf{v})+b(p,\mathbf{v})-b(q,\mathbf{u})+⟨\frac{\nu }{ϵ}\mathbf{u}{\mathbf{P}}^t,\mathbf{v}{⟩}_{\Gamma }=l(\mathbf{v})+⟨(\mathbf{h}+\frac{\nu }{ϵ}\mathbf{g}){\mathbf{P}}^t,\mathbf{v}{⟩}_{\Gamma }\forall (\mathbf{v},q)\in {\mathcal{V}}_0^n\times \mathcal{Q},$$

$$a(\mathbf{u},\mathbf{v})$$

$$b(p,\mathbf{v})$$

$$l(\mathbf{v})$$

$${\mathcal{T}}_h:=\{T\in {\hat{\mathcal{T}}}_h:T\cap \Omega \ne \varnothing \},$$

$${\mathcal{T}}_{\Gamma }:=\{T\in {\mathcal{T}}_h:T\cap \Gamma \ne \varnothing \}.$$

$${\mathcal{F}}_{\Gamma }:=\{F\in {\mathcal{F}}_i:T_F^{+}\cap \Gamma \ne \varnothing \vee T_F^{-}\cap \Gamma \ne \varnothing \}.$$

$${\mathcal{X}}_h$$

$${\mathcal{V}}_h=[{\mathcal{X}}_h{]}^d,{\mathcal{Q}}_h={\mathcal{X}}_h,{\mathcal{W}}_h={\mathcal{V}}_h\times {\mathcal{Q}}_h$$

$$-⟨(2\nu \mathbf{ϵ}({\mathbf{u}}_h)\mathbf{n}){\mathbf{P}}^t,{\mathbf{v}}_h{⟩}_{\Gamma }=⟨(-\mathbf{h}+{\epsilon }^{-1}\nu (\mathbf{u}-\mathbf{g})){\mathbf{P}}^t,{\mathbf{v}}_h{⟩}_{\Gamma }.$$

$$A_h(U_h,V_h)+S_h(U_h,V_h)+G_h(U_h,V_h)$$

$$A_h(U_h,V_h)$$

$$a_h({\mathbf{u}}_h,{\mathbf{v}}_h)$$

$$-{\zeta }_u⟨{\mathbf{u}}_h\cdot \mathbf{n},(2\nu \mathbf{ϵ}({\mathbf{v}}_h)\mathbf{n})\cdot \mathbf{n}{⟩}_{\Gamma }+⟨\frac{\nu }{{\gamma }^{n}h}{\mathbf{u}}_h\cdot \mathbf{n},{\mathbf{v}}_h\cdot \mathbf{n}{⟩}_{\Gamma }$$

$$+⟨\frac{{\varphi }_u}{{\gamma }^{n}h}{\mathbf{u}}_h\cdot \mathbf{n},{\mathbf{v}}_h\cdot \mathbf{n}{⟩}_{\Gamma }-⟨(\mathbf{\beta }\cdot \mathbf{n}){\mathbf{u}}_h,{\mathbf{v}}_h{⟩}_{{\Gamma }_{\mathrm{in}}}$$

$$+⟨\frac{1}{\epsilon +{\gamma }^{t}h}\epsilon (2\nu \mathbf{ϵ}({\mathbf{u}}_h)\mathbf{n}){\mathbf{P}}^t,{\mathbf{v}}_h{⟩}_{\Gamma }+⟨\frac{1}{\epsilon +{\gamma }^{t}h}\nu {\mathbf{u}}_h{\mathbf{P}}^t,{\mathbf{v}}_h{⟩}_{\Gamma }$$

$$-{\zeta }_u⟨\frac{{\gamma }^{t}h}{\epsilon +{\gamma }^{t}h}\epsilon (2\nu \mathbf{ϵ}({\mathbf{u}}_h)\mathbf{n}){\mathbf{P}}^t,2\mathbf{ϵ}({\mathbf{v}}_h)\mathbf{n}{⟩}_{\Gamma }-{\zeta }_u⟨\frac{{\gamma }^{t}h}{\epsilon +{\gamma }^{t}h}\nu {\mathbf{u}}_h{\mathbf{P}}^t,2\mathbf{ϵ}({\mathbf{v}}_h)\mathbf{n}{⟩}_{\Gamma },$$

$$b_h(p_h,{\mathbf{v}}_h)$$

$$L_h(V_h)$$

$$-{\zeta }_u⟨\mathbf{g}\cdot \mathbf{n},(2\nu \mathbf{ϵ}({\mathbf{v}}_h)\mathbf{n})\cdot \mathbf{n}{⟩}_{\Gamma }+⟨\frac{\nu }{{\gamma }^{n}h}\mathbf{g}\cdot \mathbf{n},{\mathbf{v}}_h\cdot \mathbf{n}{⟩}_{\Gamma }$$

$$-⟨\mathbf{g}\cdot \mathbf{n},q_h{⟩}_{\Gamma }+⟨\frac{{\varphi }_u}{{\gamma }^{n}h}\mathbf{g}\cdot \mathbf{n},{\mathbf{v}}_h\cdot \mathbf{n}{⟩}_{\Gamma }-⟨(\mathbf{\beta }\cdot \mathbf{n})\mathbf{g},{\mathbf{v}}_h{⟩}_{{\Gamma }_{\mathrm{in}}}$$

$$+⟨\frac{1}{\epsilon +{\gamma }^{t}h}\epsilon \mathbf{h}{\mathbf{P}}^t,{\mathbf{v}}_h{⟩}_{\Gamma }+⟨\frac{1}{\epsilon +{\gamma }^{t}h}\nu \mathbf{g}{\mathbf{P}}^t,{\mathbf{v}}_h{⟩}_{\Gamma }$$

$$-{\zeta }_u⟨\frac{{\gamma }^{t}h}{\epsilon +{\gamma }^{t}h}\epsilon \mathbf{h}{\mathbf{P}}^t,2\mathbf{ϵ}({\mathbf{v}}_h)\mathbf{n}{⟩}_{\Gamma }-{\zeta }_u⟨\frac{{\gamma }^{t}h}{\epsilon +{\gamma }^{t}h}\nu \mathbf{g}{\mathbf{P}}^t,2\mathbf{ϵ}({\mathbf{v}}_h)\mathbf{n}{⟩}_{\Gamma },$$

$${\varphi }_{u,T}:=\nu +\Vert \mathbf{\beta }{\Vert }_{0,\infty ,T}h+\sigma h^2.$$

$$a_h({\mathbf{u}}_h,{\mathbf{v}}_h)$$

$$-{\zeta }_u⟨{\mathbf{u}}_h\cdot \mathbf{n},(2\nu \mathbf{ϵ}({\mathbf{v}}_h)\mathbf{n})\cdot \mathbf{n}{⟩}_{\Gamma }+⟨\frac{\nu }{{\gamma }^{n}h}{\mathbf{u}}_h\cdot \mathbf{n},{\mathbf{v}}_h\cdot \mathbf{n}{⟩}_{\Gamma }$$

$$+⟨\frac{{\varphi }_u}{{\gamma }^{n}h}{\mathbf{u}}_h\cdot \mathbf{n},{\mathbf{v}}_h\cdot \mathbf{n}{⟩}_{\Gamma }-⟨(\mathbf{\beta }\cdot \mathbf{n}){\mathbf{u}}_h,{\mathbf{v}}_h{⟩}_{{\Gamma }_{\mathrm{in}}}$$

$$-{\zeta }_u⟨{\gamma }^{t}h(2\nu \mathbf{ϵ}({\mathbf{u}}_h)\mathbf{n}){\mathbf{P}}^t,2\mathbf{ϵ}({\mathbf{v}}_h)\mathbf{n}{⟩}_{\Gamma },$$

$$L_h(V_h)$$

$$-{\zeta }_u⟨\mathbf{g}\cdot \mathbf{n},(2\nu \mathbf{ϵ}({\mathbf{v}}_h)\mathbf{n})\cdot \mathbf{n}{⟩}_{\Gamma }+⟨\frac{\nu }{{\gamma }^{n}h}\mathbf{g}\cdot \mathbf{n},{\mathbf{v}}_h\cdot \mathbf{n}{⟩}_{\Gamma }$$

Full text

††journal: journal

url]www.lnm.mw.tum.de/staff/magnus-winter

url]www.lnm.mw.tum.de/staff/benedikt-schott

url]www.andremassing.com

url]www.lnm.mw.tum.de/staff/wall

1 Introduction

For an incompressible Newtonian fluid in $\Omega\subset\mathbb{R}^{d}$, $d=2,3$, which is a bounded, connected domain with the boundary $\partial\Omega$, the following Navier-Stokes equations with a general Navier boundary condition are valid,

[TABLE]

where ${\boldsymbol{\epsilon}}(\boldsymbol{u}):=\frac{1}{2}(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{T})$ and $\boldsymbol{u},\boldsymbol{g},\boldsymbol{h}:\left[0,T\right]\times\Omega\rightarrow\mathbb{R}^{d}$ and $\nu,p:\left[0,T\right]\times\Omega\rightarrow\mathbb{R}$. Here $\boldsymbol{u}$ is the velocity and $p$ and $\nu$ the pressure and kinematic viscosity of the fluid respectively. The terms $\boldsymbol{g}$ and $\boldsymbol{h}$ are the velocity and traction prescribed at the boundary $\partial\Omega$. The projection matrices ($\boldsymbol{P}^{n},\boldsymbol{P}^{t}\in\mathbb{R}^{d\times d}$) are constructed from the outward pointing unit normal $\boldsymbol{n}$ of the boundary $\partial\Omega$. The projection in normal direction is defined as $\boldsymbol{P}^{n}:=\boldsymbol{n}\otimes\boldsymbol{n}$ and the projection onto the tangential plane of the surface as $\boldsymbol{P}^{t}:=\boldsymbol{I}-\boldsymbol{n}\otimes\boldsymbol{n}$, where $\boldsymbol{I}$ is the $d\times d$ identity matrix.

The boundary condition applied to the problem (1.1)-(1.2) is the general Navier boundary condition. It prescribes a Dirichlet condition in the normal direction (1.3), where the normal velocity of the fluid needs to be the same as the normal velocity prescribed at the boundary. In the tangential plane (1.4), a Robin condition is imposed, which is a condition comprised of a linear combination of a Dirichlet and Neumann condition. The slip length parameter $\varepsilon:\partial\Omega\rightarrow\mathbb{R}_{0}^{+}\cup\{\infty\}$ determines the influence of the two parts in the tangential direction. In the limiting case where $\varepsilon\rightarrow 0$, the formulation reduces to the classic Dirichlet condition $\boldsymbol{u}=\boldsymbol{g}\text{ on }\partial\Omega$; and in the case where $\varepsilon\rightarrow\infty$, a combination of a Dirichlet condition in normal direction $\boldsymbol{u}\boldsymbol{P}^{n}=\boldsymbol{g}\boldsymbol{P}^{n}$ and a Neumann condition in tangential direction $2\nu{\boldsymbol{\epsilon}}(\boldsymbol{u})\boldsymbol{n}\boldsymbol{P}^{t}=\boldsymbol{h}\boldsymbol{P}^{t}$ is enforced.

This type of boundary condition was first proposed by Navier [3] with $\boldsymbol{g}=\boldsymbol{h}=\bf 0$ as a boundary condition for incompressible fluids at solid walls. However, in the majority of cases, it has been verified that the slip length is negligibly small and as such the no-slip condition, which is a pure Dirichlet type condition with $\boldsymbol{u}=\boldsymbol{g}$, is normally used to describe the motion of fluids. Nevertheless, there are some cases where a Navier slip type model has merit, such as for flow over super-hydrophobic surfaces [4] and in the modeling of rough surfaces [5]. It can also be shown that for contact between smooth rigid bodies in an incompressible fluid, the no-slip condition is not a feasible boundary condition [6]. One option to alleviate this issue is to employ a Navier slip boundary condition at the boundaries of the rigid bodies [7, 8]. Furthermore, in modeling the motion of contact lines for multiphase flows, it is well known that the use of the no-slip boundary condition yields non-physical infinite dissipation in the vicinity of the contact line [9]. To overcome this issue, an alternative is to model the contact line motion by using a Navier slip type boundary condition [10, 11].

The work done on how to impose the Navier boundary condition can be divided based on how it enforces the no-penetration condition (1.3), i.e. the condition in normal direction. In the work done by Verfürth, a strong imposition is proposed in [12], and by means of a weak imposition through Lagrange multipliers in [13, 14]. The constraint in normal direction can also be enforced weakly by means of Nitsche’s method [15] or by a penalty method [16]. In all the previously mentioned works, the tangential condition (1.4) is treated by substituting the traction in the variational formulation at the boundary with the following,

[TABLE]

This method of imposing the boundary condition leads to numerical difficulties when $\varepsilon\rightarrow 0$ and is not defined for a Dirichlet condition $\varepsilon=0$. Thus, the slip length can not be chosen arbitrarily, and it is difficult to predict at what value numerical issues will arise. The same problem can be observed in the Poisson equation with Robin boundary conditions, for which a solution was presented by Juntunen and Stenberg [1]. In their formulation, the Robin condition is imposed weakly by means of Nitsche’s method, which permits any choice of $\varepsilon$. Based on this method, the present work proposes a formulation to extend these results to the Navier-Stokes equations.

A few examples have already been mentioned, such as contact between rigid bodies in an incompressible fluid and moving contact lines in multi-phase flows, where large deformation and topological changes occur and where the Navier slip boundary condition has merit. To simulate these types of problems, an unfitted mesh instead of a boundary-fitted mesh method is advantageous. This is due to the fact that the boundary-fitted mesh methods, which require Arbitrary-Lagrangian-Eulerian (ALE) based mesh moving algorithms, will break down and necessitate expensive remeshing [17] under those conditions.

Unfitted mesh methods have already successfully been used in simulating a variety of different problems. These include single-phase [18] and multi-phase flows [19, 20, 21, 22], which were conducted by using an implicit level set to describe the interface between the fluids and an unfitted finite element approach for solving the governing equations. For fluid-structure interaction problems with large deformations of the structure [23, 24, 25, 26], and in the case of contact of submerged bodies [27], a fixed-grid Eulerian approach for the fluid and a Lagrangian description for the structure have been applied successfully. Unfitted mesh methods also offer the possibility to use embedding meshes [28, 29, 30, 31] to simplify meshing or to improve accuracy in certain regions of the mesh.

The Cut Finite Element Method (CutFEM) is the unfitted mesh method of choice in this work. This method stems from the eXtended Finite Element Method [32, 33]; for an overview of CutFEM see, e.g., [34]. To prescribe boundary conditions on non-fitted domain boundaries, a common choice is to impose them weakly. This can be done for instance with Nitsche’s method [35]. Stability and a priori error estimates have been derived for this setup for the Poisson [36], Stokes [37] and recently for the Oseen problem [2] with weakly imposed Dirichlet boundary conditions. Another advantage of Nitsche’s method is its flexibility and straightforward extensibility to coupled multiphysics problems [26], where the strength of the imposition of the boundary conditions is regulated by the resolution of the physical variables of the computational mesh.

The Oseen equations, which are addressed in this work, can be viewed as a linearized version of the Navier-Stokes equations when time-stepping methods are applied. Since this work is also investigating the stability and a priori error properties of an unfitted Nitsche’s method for the Oseen problem, but with a more general type of boundary condition, a large portion of the analysis relies on the work done in [2]. As the proposed method uses equal order continuous finite elements for the velocity and pressure, stabilization of the discretized equations is a necessity. The method chosen here is the continuous interior penalty (CIP) method by Burman et al. [38]; see [39] for an overview of different stabilization methods applicable to the Oseen problem. Further stabilization is necessary at the boundaries as these are cut out from the background mesh. If such stabilization is not applied, the method is not necessarily stable for pathological cut cases. Our choice of stabilization at the boundary is the ghost-penalty method introduced by Burman [40] for the Poisson equation and adapted to the Oseen and Navier-Stokes equations by [2] and [18].

The contributions of this work can be summarized as follows: A weakly imposed general Navier boundary condition is introduced. It consists of a normal and a tangential component, which are separated by the use of projection operators at the boundary. In the normal direction a Dirichlet condition is imposed by means of Nitsche’s method. In the tangential plane a Robin condition is also imposed by a Nitsche’s method, which has the advantage over the classical Galerkin substitution (1.5) in that it remains stable and well-posed for all slip lengths $\varepsilon$. The formulation in the tangential plane is inspired by the work of [1], where a similar problem was solved for the Robin boundary condition for the Poisson equation. A numerical analysis is conducted to prove that our formulation is stable and has optimal convergence behavior. The analysis builds upon the work done in [2], where stability and optimality of the error estimates for the Oseen equations with CutFEM was done for a weak imposition of a Dirichlet boundary condition. Here, however, the analysis is carried out for a more general case of the boundary condition.

The paper is organized as follows. Section 2 introduces the Oseen problem, and the necessary spaces and assumptions are specified; also the variational formulation of the problem is presented. In Section 3 the definitions for the cut finite element spaces are introduced. The general Navier boundary condition imposed by Nitsche’s method is presented and its differences and advantages are highlighted against alternative methods. The stabilized weak form of the Oseen equations, discretized by means of the CutFEM and the necessary stabilizations of the domain through the CIP method and the cut boundaries by means of ghost-penalty stabilization, are also explained in this section. In Section 4 basic approximation properties and interpolation operators and norms are introduced. Sections 5 and 6 are dedicated to proving the inf-sup stability and the optimal a priori error of the proposed method. In Section 7 the method is applied to a numerical example and a convergence study is conducted to verify the theoretical results. The final Section 8 gives a summary and conclusion of this work.

2 The Oseen Problem with General Navier Boundary Conditions

2.1 Basic Notation on Function Spaces

We define our domain as $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, where $\Omega$ is either a bounded convex domain or a plane bounded domain with Lipschitz and piecewise $C^{2}$ boundary with convex angles. The standard Sobolev spaces are denoted by $W^{m,q}(U)$ for $U\in\{\Omega,\Gamma\}$ where $0\leqslant m<\infty$ and $1\leqslant q\leqslant\infty$ with associated norms $\|\cdot\|_{m,q,U}$. We write $H^{m}(U)=W^{m,2}(U)$ with norm $\|\cdot\|_{m,U}$. For the associated inner products we write $(\cdot,\cdot)_{m,U}$ for measurable subsets $U\subseteq\mathbb{R}^{d}$ and $\langle\cdot,\cdot\rangle_{m,U}$ for subsets $U\subseteq\mathbb{R}^{d-1}$. Occasionally, we write $(\cdot,\cdot)_{U}$, $\langle\cdot,\cdot\rangle_{U}$ and $\|\cdot\|_{U}$ if $m=0$. Fractional Sobolev trace spaces $[H^{s-\frac{1}{2}}(Y)]^{d}$ at subsets $Y\subset\mathbb{R}^{d-1}$, which in practice will be a part of the domain boundary $\Gamma$, denote the set of boundary traces of all $\mathbb{R}^{d}$-valued functions in $[H^{s}(\Omega)]^{d}$. For the Oseen problem, we make use of the specific function spaces $H_{0}(\nabla\cdot;\Omega)\subset[L^{2}(U)]^{d}$, which denotes the space of divergence-free functions, and $L^{2}_{0}(\Omega)$, which denotes the function space consisting of functions in $L^{2}(\Omega)$ with zero average on $\Omega$. To shorten the presentation, for broken norms and inner products which involve a collection of geometric entities $\mathcal{P}_{h}$, i.e. finite elements or facets, we write $\|\cdot\|_{\mathcal{P}_{h}}^{2}=\sum_{P\in\mathcal{P}_{h}}\|\cdot\|_{P}^{2}$ whenever $\|\cdot\|_{P}$ is well-defined.

2.2 Problem Formulation

Considering the non-linear Navier-Stokes equations (1.1)–(1.2), after applying a time discretization method and a linearization step, many solution algorithms can be reduced to solving a sequence of auxiliary problems of Oseen type for the velocity field $\boldsymbol{u}:\Omega\rightarrow\mathbb{R}^{d}$ and the pressure field $p:\Omega\rightarrow\mathbb{R}$:

[TABLE]

where ${\boldsymbol{\epsilon}}(\boldsymbol{u}):=\frac{1}{2}(\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{T})$ denotes the strain rate tensor, $\boldsymbol{\beta}\in[W^{1,\infty}(\Omega)]^{d}\cap H_{0}(\nabla\cdot;\Omega)$ the given divergence-free advective velocity field and $\boldsymbol{f}\in[L^{2}(\Omega)]^{d}$ the body force. The assumption is made that the reaction coefficient $\sigma$ and the viscosity $\nu$ are positive real-valued constants. For the boundary conditions (2.3)–(2.4) the functions are defined as

[TABLE]

and the normal and tangential projection matrices are constructed from the outward pointing unit normal $\boldsymbol{n}$ of the boundary $\Gamma$ as $\boldsymbol{P}^{n}:=\boldsymbol{n}\otimes\boldsymbol{n}$ and $\boldsymbol{P}^{t}:=\boldsymbol{I}-\boldsymbol{n}\otimes\boldsymbol{n}$, where $\boldsymbol{I}$ is the $d\times d$ identity matrix. Utilizing a non-negative slip length function $\varepsilon:\Gamma\rightarrow\mathbb{R}_{0}^{+}\cup\{\infty\}$, the classical full Dirichlet boundary conditions can be recovered by setting $\varepsilon=0$, which for $\nu>0$ states that $\boldsymbol{u}=\boldsymbol{g}\text{ on }\Gamma_{D}:=\{\boldsymbol{x}\in\Gamma\text{ with }\varepsilon=0\}\subseteq{\Gamma}$; however, the tangential condition (2.4) does not contribute in the Darcy limit, i.e. for $\nu=0$, and as such, the normal condition (2.3) is sufficient to define a Dirichlet boundary condition. By choosing $\varepsilon=\infty$ and $\boldsymbol{h}=\bf 0$ the boundary condition reduces to a full-slip condition as the tangential velocity on $\Gamma$ is not constrained anymore. The part of $\Gamma$ with $0\leqslant\varepsilon<\infty$ is denoted by $\Gamma^{\varepsilon}$ and is assumed being non-vanishing, i.e. $\operatorname{meas}(\Gamma^{\varepsilon})>0$, to avoid not well-posed pure slip-boundary problems. It is further assumed that the inflow boundary $\Gamma_{\mathrm{in}}:=\{\boldsymbol{x}\in\Gamma\leavevmode\nobreak\ |\leavevmode\nobreak\ \boldsymbol{\beta}\cdot\boldsymbol{n}<0\}$ is part of the Dirichlet boundary, i.e. $\Gamma_{\mathrm{in}}\subseteq\Gamma_{D}$. In the discrete setting, the constant pressure mode needs to be filtered out during the solution procedure, as it is determined only up to a constant.

2.3 Variational Formulation

Let us denote the velocity and pressure space by

[TABLE]

The corresponding weak formulation of the Oseen problem (2.1)–(2.4) is to find the velocity and the pressure field $(\boldsymbol{u},p)\in\mathcal{V}_{\boldsymbol{g}}\times\mathcal{Q}$ such that

[TABLE]

where

[TABLE]

For the well-posedness and solvability of the continuous Oseen problem (2.9), we refer the reader to, e.g. [41, 42, 14].

3 A Stabilized Nitsche-type Cut Finite Element Method for the Oseen Problem

This section is devoted to the presentation of our cut finite element method. After defining suitable cut finite element function spaces, we address the weak imposition of general Navier boundary conditions for the Oseen problem by a Nitsche-type method. A possible stabilization technique to overcome issues of classical finite element approximations is recalled and we discuss how the resulting discrete formulation can be extended to non-boundary-fitted approximation with the help of boundary-zone ghost-penalty stabilizations. Finally, we summarize our cut finite element formulation and introduce suitable norms for the numerical stability and a priori error analysis.

3.1 Computational Meshes and Cut Finite Element Spaces

While for standard finite element methods the computational mesh is fitted to the boundary or an interpolatory approximation is given, in cut finite element methods the boundary is allowed to intersect the mesh. For the sake of simplicity, in this work, we assume a sequence of quasi-uniform meshes $\widehat{\mathcal{T}}_{h}=\{T\}$, each consisting of shape-regular finite elements $T$ with mesh size parameter $h$ and covering the physical domain $\Omega$. For each background mesh $\widehat{\mathcal{T}}_{h}$, the finite element solution is then approximated on an active part of the background mesh

[TABLE]

consisting of all elements in $\widehat{\mathcal{T}}_{h}$ which intersect the physical domain $\Omega$. The possibly enlarged domain which is covered by the union of all elements $T\in\mathcal{T}_{h}$ is denoted by $\Omega^{\ast}_{h}$. A mesh $\mathcal{T}_{h}$ is called a fitted mesh if $\overline{\Omega}=\overline{\Omega^{\ast}_{h}}$ and an unfitted mesh if $\overline{\Omega}\subsetneq\overline{\Omega^{\ast}_{h}}$. The subset of all elements of the respective active mesh which are located in the vicinity of the boundary $\Gamma$ are denoted as

[TABLE]

For stabilization purposes we need to define the set of all facets by $\mathcal{F}_{h}$, where $\mathcal{F}_{i}$ is the set of all interior facets $F$ which are shared by exactly two elements, denoted by $T^{+}_{F}$ and $T^{-}_{F}$. The notation $\mathcal{F}_{\Gamma}$ is used for the set of all interior facets belonging to elements intersected by the boundary $\Gamma$,

[TABLE]

Figure  3.1 summarizes the notation.

By analogy to other preceding works on cut finite element methods [43, 36, 44, 45, 2], we assume that the domain and the boundary are reasonably resolved by each quasi-uniform mesh $\mathcal{T}_{h}$. This entails that the intersection between the facet $F\in\mathcal{F}_{i}$ and the boundary $\Gamma$ is simply connected. Furthermore, there needs to exist a plane $S_{T}$ with a piecewise smooth parametrization $\Phi:S_{T}\cap T\rightarrow\Gamma\cap T$ for each element $T$ intersected by $\Gamma$. Finally an assumption is made that for each element $T\in\mathcal{T}_{\Gamma}$ there exists $T^{\prime}\in\mathcal{T}_{h}\setminus\mathcal{T}_{\Gamma}$ such that the sequence $\{T\}_{j=1}^{N}$ with $T_{1}=T,\,T_{N}=T^{\prime}$ and $T_{j}\cap T_{j+1}\in\mathcal{F}_{i},\;j=1,\ldots,N-1$ is at most $N$ elements.

Let $\mathcal{X}_{h}$ be the finite element spaces consisting of continuous piecewise polynomials of order $k$ for a given mesh $\mathcal{T}_{h}$

[TABLE]

For the discrete cut finite element approximation of the solution to the Oseen problem, we use equal-order interpolations for velocity and pressure, where

[TABLE]

are the discrete velocity space, the discrete pressure space and the total approximation space.

3.2 Weak Imposition of Generalized Navier Boundary Conditions

In this section, we present the imposition of generalized Navier boundary conditions for the Oseen problem. As previously mentioned, this boundary condition consists of a Dirichlet condition in normal direction to the boundary and a Robin condition in the tangential plane of the boundary. Since the boundary-normal constraint (2.3) is of Dirichlet-type, either a strong or a weak imposition can be used. Nevertheless, in the case of an unfitted mesh, a weak imposition is recommended. In the case when the slip-length coefficient is sufficiently large, the tangential constraint (2.4) can be imposed by direct substitution of the boundary terms, which stems from integrating the viscous bulk term in (2.1) by parts, i.e. by substituting

[TABLE]

Such a formulation provides a reasonable and stable method for large slip-length coefficients $\varepsilon$. However, it is well known that for $\varepsilon\rightarrow 0$, i.e. when we approach a Dirichlet condition in the tangential direction, the conditioning of the formulation deteriorates and, as a result, does not provide an accurate and stable formulation, see discussion by Juntunen and Stenberg [1]. To overcome this issue in the case of the Poisson problem, Juntunen and Stenberg proposed a Nitsche’s method of imposing a Robin boundary condition, which provides stability, optimal a priori error estimates and a bounded conditioning w.r.t. small choices of $\varepsilon\in[0,\infty]$.

Next, we present our novel Nitsche-type method for the Oseen problem with general Navier boundary conditions. It combines the advantages of the Nitsche-type method by Juntunen and Stenberg [1] for generalized Robin-type conditions with the Nitsche formulations developed by Burman et al. [38] and Massing et al. [2] for incompressible low- and high-Reynolds-number flow with Dirichlet-type constraints.

The unfitted Nitsche-type finite element formulation for the Oseen problem is to find $U_{h}=(\boldsymbol{u}_{h},p_{h})\in\mathcal{W}_{h}$ such that for all $V_{h}=(\boldsymbol{v}_{h},q_{h})\in\mathcal{W}_{h}$

[TABLE]

where

[TABLE]

where $\phi_{u}:\Omega_{h}^{*}\rightarrow\mathbb{R}$ is defined on each element as

[TABLE]

The stabilization operators $S_{h},G_{h}$ will be specified later in Sections 3.3 and 3.4. Starting from the weak formulation (2.9), it can be seen that due to the weak constraint enforcement the standard consistency boundary terms remain in the momentum equation (3.9) and (3.14). Equivalent to a Nitsche formulation for a pure Dirichlet boundary condition, these standard consistency terms on the boundary $\Gamma$ are potential sources for instabilities and are analyzed further in Sections 5 and 6. The constraints being imposed in wall-normal and tangential directions, i.e. (2.3) and (2.4), are enforced by adding additional consistent boundary terms.

Wall-normal constraint: The enforcement of the boundary-normal constraint follows the standard Nitsche-technique by adding an adjoint (in-)consistent viscous term (depending on $\zeta_{u}\in\{-1,1\})$ and a consistent and optimal convergent viscous symmetric penalty term (3.10) with appropriate right-hand-side terms (3.16). Choosing an adjoint-consistent method ($\zeta_{u}=1$), the discrete formulation requires the penalty parameter $0<\gamma^{n}\leqslant C$ to be chosen small enough, where the constant $C$ depends on the shape and polynomial order of the finite element, cf. the trace estimate (4.1) in Section 4. Furthermore, this choice allows to deduce optimal convergence for the velocity $L^{2}$-error as stated in Theorem 6.5. Even though an adjoint-inconsistent formulation ($\zeta_{u}=-1$) enjoys improved inf-sup stability for any $0<\gamma^{n}<\infty$ and thereby ensures optimal convergence w.r.t. an energy norm, optimality for the velocity $L^{2}$-error is not guaranteed in this case anymore. For an analysis of penalty-free Nitsche methods, i.e. for $\gamma_{n}=\infty$, the interested reader is referred to, e.g. [46, 47]. Note that even for vanishing viscosity the wall-normal constraints needs to be enforced. Thus, to ensure inf-sup stability and optimal convergence for all flow regimes, an adjoint-consistent pressure term (3.14) and a symmetric penalty term (3.11), which accounts for the different flow regimes as reflected by the definition of the piecewise constant stabilization scaling function $\phi_{u}$, are added to the formulation. The respective right-hand-side terms (3.17) again guarantee consistency of the method. Note that additionally at inflow boundaries $\Gamma_{\mathrm{in}}\subseteq\Gamma_{D}$ where it holds that $\varepsilon=0$ and $\boldsymbol{\beta}\cdot\boldsymbol{n}<0$, Dirichlet boundary conditions need to be imposed in all spatial directions as stated by the additional advective inflow stabilization terms in (3.11) and (3.17).

Wall-tangential constraint: The enforcement of the tangential boundary condition follows the technique introduced by Juntunen and Stenberg [1]. By adding terms where the tangential condition (2.4) is tested with $\boldsymbol{v}_{h}$ and ${\boldsymbol{\epsilon}}(\boldsymbol{v}_{h})$, consistency is ensured. Thereby, choosing the weights of these added terms as presented, guarantees coercivity, optimal a priori error estimates and bounded system conditioning w.r.t. the choice of the slip-length coefficient $\varepsilon\in[0,\infty]$. Similar to the viscous wall-normal constraint enforcement, the choice between an adjoint-consistent ($\zeta_{u}=1$) and an adjoint-inconsistent ($\zeta_{u}=-1$) formulation poses equivalent restrictions to $\gamma^{t}$, i.e. $0<\gamma^{t}\leqslant C$ for $\zeta_{u}=1$ (with $C$ stemming from estimate (4.1)) and $0<\gamma^{t}<\infty$ for $\zeta_{u}=-1$.

Remark 3.1**.**

In the limiting case of $\varepsilon\rightarrow 0$, it can easily be seen that the formulation (3.7) ends up the same as already presented in [2] for the imposition of weak Dirichlet boundary conditions by means of Nitsche’s method.

Remark 3.2**.**

For the case where $\varepsilon\rightarrow\infty$ and $\gamma^{t}h\ll\varepsilon$ the formulation does not return the familiar imposition of a Neumann condition in the tangential plane. For clarity, we explicitly state the terms $a_{h}(\boldsymbol{u}_{h},\boldsymbol{v}_{h})$ and $L_{h}(V_{h})$ here,

[TABLE]

This limiting case gives a standard Neumann condition as in (3.26) with the addition of the terms (3.24) and (3.28). These terms are, however, added consistently to the formulation and as such do not ruin its validity. Nevertheless, these terms are of importance as they impose limitations on the choice of $\gamma^{t}$. The conditioning of this formulation becomes increasingly bad as $\gamma^{t}$ is increasing. From this it is clear that the parameter is effectively limited from above even in the case of an adjoint-inconsistent formulation, i.e. $\zeta_{u}=-1$, which is proven to be inf-sup stable for $0\leqslant\gamma^{t}<\infty$ in Section 5. Another observation is that in the other limiting case where $\gamma^{t}\rightarrow 0$ we end up with a classical standard Galerkin imposition of a Neumann boundary condition in the tangential plane, as the terms (3.24) and (3.28) disappear.

3.3 The Continuous Interior Penalty (CIP) Stabilizations

It is well-known that a finite element based discretization of formulation (2.9) needs to be stabilized to allow for equal-order interpolation spaces $\mathcal{W}_{h}=\mathcal{V}_{h}\times\mathcal{Q}_{h}$ due to its saddle-point structure, see e.g. [48], and to suppress spurious oscillations in the numerical solution in the case of convection-dominant flow. For a detailed overview of different stabilization techniques, see e.g. the overview article [39].

Following the presentation of the cut finite element formulation from [2], we apply continuous interior penalty (CIP) stabilization operators originally introduced by Burman et al. [38] and adapted by Massing et al. [2] such that

[TABLE]

consists of symmetric stabilization terms, which penalize the jump of the velocity and pressure gradients over interior element facets $F\in\mathcal{F}_{i}$. The stabilization operators are defined by

[TABLE]

where for any, possibly vector-valued, piecewise discontinuous function $\phi$ on the computational mesh $\mathcal{T}_{h}$, we denote the jump and average over an interior facet $F\in\mathcal{F}_{i}$ with $\llbracket\phi\rrbracket:=(\phi_{F}^{+}-\phi_{F}^{-})$ and $\phi_{F}:=\tfrac{1}{2}(\phi_{F}^{+}+\phi_{F}^{-})$ where $\phi^{\pm}(\boldsymbol{x})=\lim_{t\to 0^{+}}\phi(\boldsymbol{x}\pm t\boldsymbol{n}_{F})$ for some chosen normal unit vector $\boldsymbol{n}_{F}$ on $F$. The element-wise constant stabilization parameter $\phi_{u}$ is as defined in (3.20) and $\phi_{\beta},\phi_{p}$ are given as

[TABLE]

Throughout this work, we use the notation $a\lesssim b$ for $a\leqslant Cb$ for some positive generic constant $C$ which varies with the context, however, is always independent of the mesh size $h$ and the intersection of the mesh $\mathcal{T}_{h}$ by $\Gamma$.

For a (Lipschitz)-continuous $\boldsymbol{\beta}\in[W^{1,\infty}(\Omega)]^{d}\subseteq[C^{0,1}(\Omega)]^{d}$, in the forthcoming numerical analysis we assume a piecewise constant approximation satisfying

[TABLE]

Furthermore, in the unfitted mesh case we assume that there exists an extension $\boldsymbol{\beta}^{\ast}\in[W^{1,\infty}(\Omega^{\ast})]^{d}$ from $\Omega$ to $\Omega^{\ast}$ satisfying $\|\boldsymbol{\beta}^{\ast}\|_{1,\infty,\Omega^{\ast}}\lesssim\|\boldsymbol{\beta}\|_{1,\infty,\Omega}$. Similar to the preceding works by [2, 38], we assume that the flow field $\boldsymbol{\beta}$ is sufficiently resolved by the mesh such that $\forall\,T\in\mathcal{T}_{h}$

[TABLE]

where $\omega(T)$ denotes a local patch of elements neighboring $T$. As a result, the piecewise constant stabilization parameters are comparable locally in a neighborhood of elements, i.e.

[TABLE]

Remark 3.3**.**

Note that for facets $F\in\mathcal{F}_{i}$, which are intersected by the boundary $\Gamma$, the inner products of all CIP stabilization operators (3.30)–(3.32) have to be evaluated along the entire cut facets.

3.4 Stabilizing Cut Elements – The Role of Ghost Penalties

To strengthen the stability properties of the discrete formulation in the boundary zone for non-boundary-fitted meshes $\mathcal{T}_{h}$, additional measures are required. So-called ghost-penalty stabilizations, comprised in the operator $G_{h}$, are active in the boundary zone and augment the stabilized bilinear form $A_{h}+S_{h}$ to account for small cut elements $|T\cap\Omega|\ll|T|,\;T\in\mathcal{T}_{h}$, in the vicinity of the boundary $\Gamma$. For detailed elaborations on this concept see, e.g., the works [40, 36, 44, 2].

Ghost-penalty stabilizations extend the stability and approximation properties of the discrete scheme to the entire active background mesh, i.e. from $\Omega$ to $\Omega^{\ast}_{h}$, and thus give control of the discrete velocity and pressure solutions, where they have no physical significance. As a major advantage, the resulting scheme has improved optimality properties with highly reduced sensitivity of the errors and guaranteed uniformly bounded conditioning irrespective of how the boundary $\Gamma$ intersects the underlying mesh $\mathcal{T}_{h}$. For this purpose, the different terms need to be designed such that inf-sup stability is guaranteed for all flow regimes and (weak) consistency and optimality of the numerical scheme is maintained.

In the preceding work by Massing et al. [2], a set of different ghost-penalty terms have been developed to stabilize the Oseen equations. The suggested terms consist of CIP-type jump penalties of polynomial order $k$ for velocity and pressure and are recalled in the following

[TABLE]

with

[TABLE]

where the $j$-th normal derivative $\partial_{\boldsymbol{n}}^{j}v$ is given by $\partial_{\boldsymbol{n}}^{j}v=\sum_{|\alpha|=j}D^{\alpha}v(\boldsymbol{x})\boldsymbol{n}^{\alpha}$ for multi-index $\alpha=(\alpha_{1},\ldots,\alpha_{d})$, $|\alpha|=\sum_{i}\alpha_{i}$ and $\boldsymbol{n}^{\alpha}=n_{1}^{\alpha_{1}}n_{2}^{\alpha_{2}}\cdots n_{d}^{\alpha_{d}}$. In contrast to the related CIP stabilizations, these terms need to control the entire discrete polynomial, this entails all the higher-order normal derivatives contained in the adjacent elements polynomials. However, they are only evaluated along facets in the vicinity of the boundary, i.e. $\forall\,F\in\mathcal{F}_{\Gamma}$. The ghost-penalty terms (3.40)–(3.42) extend control over the respective instabilities to the enlarged domain $\Omega_{h}^{\ast}$, whereas the reactive and viscous ghost-penalties (3.38)–(3.39) extend control over the scaled $L^{2}$- and $H^{1}$-norms in the following fashion:

Proposition 3.4**.**

Let $\Omega$, $\Omega^{\ast}_{h}$ and $\mathcal{F}_{\Gamma}$ be defined as in Section 3.1 and the ghost penalty operators be given as in (3.38)–(3.42). Then for scalar functions $p_{h}\in\mathcal{Q}_{h}$ as well as for vector-valued equivalents $\boldsymbol{u}_{h}\in\mathcal{V}_{h}$ the following estimates hold

[TABLE]

Furthermore, let the scaling functions $\phi_{\beta},\phi_{u}$ and $\phi_{p}$ be defined as in (3.20) and (3.33) and let $\boldsymbol{\beta}_{h}^{0}\in[\mathcal{X}_{h}^{\mathrm{dc},0}]^{d}$ be a piecewise constant approximation to $\boldsymbol{\beta}$ on $\mathcal{T}_{h}$, which satisfies the approximation properties specified in (3.34). Then the following estimates hold

[TABLE]

*with the non-dimensional scaling functions $\omega_{h}$ and $\Phi$ given as in [2] by *

[TABLE]

Note that the hidden constants in (3.45)–(3.48) depend only on the shape-regularity and the polynomial order, but not on the mesh or the location of $\Gamma$ within $\mathcal{T}_{h}$.

Proof. For detailed proofs of these estimates, the reader is referred to the preceding work by Massing et al. [2] (cf. Lemma 5.4, Corollary 5.5 and 5.7 and Lemma 5.8) and the references therein.

Remark 3.5**.**

Note that the function $\boldsymbol{\beta}$ occurring in the stabilization operators $s_{\beta}$ (3.30) and $g_{\beta}$ (3.40) can be replaced by a proper continuous interpolation $\boldsymbol{\beta}_{h}$ satisfying the assumption specified in (3.34), without changing the final results of the stability and a priori estimates. When solving the incompressible Navier-Stokes equations, $\boldsymbol{\beta}_{h}$ is the finite element approximation of $\boldsymbol{u}$ from a previous time or iteration step.

Remark 3.6**.**

As suggested in [2], the convection and the incompressibility related stabilization forms $s_{\beta},s_{u}$ and $g_{\beta},g_{u}$, i.e. (3.30), (3.31) and (3.40), (3.41), can be replaced by single operators

[TABLE]

with $\overline{\phi}_{\beta}=\|\boldsymbol{\beta}\|_{0,\infty,F}^{2}\phi_{\beta}$, which simplifies the implementation of the purposed method.

3.5 Final Discrete Formulation and Norms

The full stabilized cut finite element method for the Oseen problem with generalized Navier boundary conditions then reads: find $U_{h}=(\boldsymbol{u}_{h},p_{h})\in\mathcal{W}_{h}$ such that $\forall\,V_{h}=(\boldsymbol{v}_{h},q_{h})\in\mathcal{W}_{h}$

[TABLE]

where $A_{h},L_{h}$ are the discrete operators including the weak constraint enforcement of the boundary conditions (see Section 3.2, (3.8)–(3.19)), $S_{h}(\cdot,\cdot)$ is the CIP operator to balance instabilities in the interior of the physical domain (see Section 3.3, (3.29)–(3.32)) and $G_{h}(\cdot,\cdot)$ the ghost-penalty operator, which extends stability control to the boundary zone when non-boundary-fitted meshes are used for the approximation (see Section 3.4, (3.37)–(3.42)).

Following the works by Massing et al. [2] on a related cut finite element method with Dirichlet boundary conditions and the work by Juntunen and Stenberg [1] on a Nitsche-type method for Robin-type constraints, for the subsequent numerical stability and convergence analysis we introduce the following (semi-)norms according to our cut finite element method (3.52). For functions $U=(\boldsymbol{u},p)$ with $\boldsymbol{u}\in H^{1}(\Omega)$ and $p\in L^{2}(\Omega)$ we define the mesh-dependent energy norms related to our Nitsche-type formulation

[TABLE]

For the stability analysis we will also utilize the augmented natural energy norm for $U_{h}\in\mathcal{V}_{h}\times\mathcal{Q}_{h}$

[TABLE]

and a semi-norm which is defined as

[TABLE]

Note that it holds that $\Phi\lesssim\phi_{u}^{-1}$, which allows us to create a lower bound of the locally scaled pressure norms. Moreover, $C_{P}$ denotes the so-called Poincaré constant as defined in (4.3) in Section 4, which scales as the diameter of $\Omega$.

For the a priori error analysis, additional control over boundary fluxes on $\Gamma$ and the divergence of the velocity is desired. For this purpose, for sufficiently regular functions $U=(\boldsymbol{u},p)\in H^{2}(\Omega)\times H^{1}(\Omega)$ we define

[TABLE]

From the inverse estimate (4.1) (see Section 4) and the norm equivalences from Proposition 3.4, discrete functions satisfy

[TABLE]

4 Preliminary Estimates

In this section, we collect important and useful estimates, which will be used frequently throughout the stability and a priori error analysis of the proposed cut finite element method (3.52) in Sections 5 and 6. First, we collect trace inequalities and inverse estimates and comment on the Korn-type inequality for the strain-rate tensor in combination with Nitsche boundary terms to enforce a generalized Navier boundary condition. Finally, we introduce suitable interpolation error estimates, which will be used to establish a priori error estimates.

4.1 Useful Inequalities and Estimates

For discrete functions $v_{h}\in\mathcal{X}_{h}$ the following well-known generalized inverse and trace inequalities hold for elements $T$ which are arbitrarily intersected by the boundary $\Gamma$:

[TABLE]

For functions $v\in H^{1}(\Omega^{\ast}_{h})$, the following trace inequalities are valid

[TABLE]

as shown in [43, 49]. Finally, assuming that $\operatorname{meas}(\Gamma^{\varepsilon})>0$, we recall the generalized Poincaré inequality for $[H^{1}(\Omega)]^{d}$ functions with non-vanishing boundary trace from [50]

[TABLE]

In the following, a generalized Korn-type inequality is presented and it is shown how to control the kernel of the strain-rate-deformation tensor ${\boldsymbol{\epsilon}}(\cdot)$ with the help of Nitsche-type penalty terms.

Proposition 4.1**.**

Let a semi-norm on $[H^{1}(\Omega)]^{d}$ be defined as $\|\boldsymbol{u}\|_{\Gamma^{\varepsilon}}^{2}:=\int_{\Gamma^{\varepsilon}}{\boldsymbol{u}^{2}}\,\mathrm{d}s$, where $\|\varepsilon\|_{\infty,\Gamma^{\varepsilon}}\leqslant c_{\varepsilon}<\infty$. Let us assume that $\operatorname{meas}(\Gamma^{\varepsilon})>0$, then $\|\cdot\|_{\Gamma^{\varepsilon}}$ defines a norm on the space of rigid body motions

[TABLE]

where $S_{d}$ is the space of anti-symmetric $d\times d$ matrices, i.e. $\boldsymbol{W}=-\boldsymbol{W}^{T}$. Note, that $RM(\Omega)$ is the kernel of the symmetric strain-rate-deformation tensor ${\boldsymbol{\epsilon}}(\cdot)$ in $[H^{1}(\Omega)]^{d}$ with $d\in\{2,3\}$.

Proof. For a more explicit description of why $RM(\Omega)$ is the kernel of the tensor ${\boldsymbol{\epsilon}}(\cdot)$ in $[H^{1}(\Omega)]^{d}$, the reader is referred to, e.g. [51, 52]. The claim regarding the norm property follows directly from the fact that $\operatorname{meas}(\Gamma^{\varepsilon})>0$ and from the uniqueness of the trivial solution $(\boldsymbol{c},\boldsymbol{W})=(\boldsymbol{0},\boldsymbol{0})\in\mathbb{R}^{d}\times S_{d}$ of the linear system $\boldsymbol{c}+\boldsymbol{W}\boldsymbol{x}=\boldsymbol{0}$ on $\Gamma^{\varepsilon}$.

Theorem 4.2** (Korn-type inequality).**

Let $\operatorname{meas}(\Gamma^{\varepsilon})>0$ and $\|\cdot\|_{\Gamma^{\varepsilon}}$ be as defined in Proposition 4.1, then there exists a constant $C_{K}$ such that

[TABLE]

Proof. The proof follows the technique proposed in [47] and [52] for Korn-type inequalities and is therefore only sketched in the following. Let $[\hat{H}^{1}(\Omega)]^{d}:=\{\boldsymbol{v}\in[H^{1}(\Omega)]^{d}:|\int_{\Omega}\nabla\times\boldsymbol{v}\,\mathrm{d}x|=0\wedge|\int_{\Omega}\boldsymbol{v}\,\mathrm{d}x|=0\}$, then it holds $[H^{1}(\Omega)]^{d}=[\hat{H}^{1}(\Omega)]^{d}\oplus RM(\Omega)$ such that for each $\boldsymbol{u}\in[H^{1}(\Omega)]^{d}$ there exists a unique pair $(\boldsymbol{z},\boldsymbol{w})\in[\hat{H}^{1}(\Omega)]^{d}\oplus RM(\Omega)$ with $\boldsymbol{u}=\boldsymbol{z}+\boldsymbol{w}$ satisfying $\|\boldsymbol{z}\|_{1,\Omega}+\|\boldsymbol{w}\|_{1,\Omega}\lesssim\|\boldsymbol{u}\|_{1,\Omega}$ (open mapping theorem by Lax [53]). Following [47] the claim can be proven via contradiction. Assuming non-boundedness in (4.5), there exists a sequence $\{\boldsymbol{u}_{n}\}\subseteq[H^{1}(\Omega)]^{d}$ with $\|\boldsymbol{u}_{n}\|_{1,\Omega}=1$ and $\|{\boldsymbol{\epsilon}}(\boldsymbol{u}_{n})\|_{\Omega}+\|\boldsymbol{u}_{n}\|_{\Gamma^{\varepsilon}}<\frac{1}{n}$. Splitting $\boldsymbol{u}_{n}=\boldsymbol{z}_{n}+\boldsymbol{w}_{n}\in[\hat{H}^{1}(\Omega)]^{d}\oplus RM(\Omega)$ for each $n$ it holds that $\|{\boldsymbol{\epsilon}}(\boldsymbol{z}_{n})\|_{\Omega}=\|{\boldsymbol{\epsilon}}(\boldsymbol{u}_{n})\|_{\Omega}<\frac{1}{n}$, since ${\boldsymbol{\epsilon}}(\boldsymbol{w}_{n})=\boldsymbol{0}$. Applying the second Korn’s inequality, see e.g. [54], we obtain

[TABLE]

Since $\boldsymbol{z}_{n}+\boldsymbol{w}_{n}$ is a bounded sequence and $RM(\Omega)$ is finite dimensional, $\{\boldsymbol{w}_{n}\}$ is a bounded sequence and as such there exists a convergent subsequence for which holds $\boldsymbol{u}=\lim_{n_{k}\to\infty}{\boldsymbol{w}_{n_{k}}}$ with $\boldsymbol{u}\in RM(\Omega)$ since $\lim_{n_{k}\to\infty}\boldsymbol{z}_{n_{k}}=\boldsymbol{0}$ in $[H^{1}(\Omega)]^{d}$. Note that by assumption $\|\boldsymbol{u}\|_{\Gamma^{\varepsilon}}=0$ and $\|\boldsymbol{u}\|_{1,\Omega}=1$. Since $\|\cdot\|_{\Gamma^{\varepsilon}}$ defines a norm on $RM(\Omega)$ whenever $\operatorname{meas}(\Gamma^{\varepsilon})>0$ (see Proposition 4.1), it follows that $\boldsymbol{u}=\boldsymbol{0}\in RM(\Omega)$ and therefore $\boldsymbol{u}=\boldsymbol{0}\in[H^{1}(\Omega)]^{d}$, which contradicts $\|\boldsymbol{u}\|_{1,\Omega}=1$. As a result, the assumption of unboundedness of (4.5) is refuted, which proves the claim.

The subsequent corollary states that the Nitsche penalty terms occurring in our Nitsche-type cut finite element formulation (3.52) are sufficient to control the rigid body motions in $[H^{1}(\Omega)]^{d}$, which remain undetermined by the strain-rate-deformation tensor ${\boldsymbol{\epsilon}}(\cdot)$.

Corollary 4.3**.**

Let the domain be bounded, i.e. $\operatorname{diam}(\Omega)<\infty$, the Nitsche penalty parameters $\gamma^{t},\gamma^{n}<\infty$ and assume that normal and tangential velocities are constrained by Nitsche-type penalty terms on a non-vanishing part of the boundary $\Gamma^{\varepsilon}$, i.e. $\operatorname{meas}(\Gamma^{\varepsilon})>0$, then

[TABLE]

Proof. The proof follows immediately from applying the Korn-type inequality deduced in Theorem 4.2, the definition of $\Gamma^{\epsilon}$ and the specified assumptions in the corollary

[TABLE]

4.2 Interpolation Operators

Since the finite element approximation space is defined on the enlarged domain $\Omega_{h}^{\ast}$, we first comment on the construction of an appropriate interpolation operator $L^{2}(\Omega)\to\mathcal{X}_{h}$. From [55] it is well known that for the Sobolev spaces $W^{m,q}(\Omega)$, $0\leqslant m<\infty$, $1\leqslant q\leqslant\infty$, a linear extension operator can be defined

[TABLE]

and we write $v^{\ast}:=Ev$. Following the analysis provided by Massing et al. [2], let $\pi_{h}$ denote the Clément operator, see for instance [56], then for $u\in H^{s}(\Omega)$ we define its “fictitious domain” extension $\pi_{h}^{\ast}:H^{s}(\Omega)\to\mathcal{X}_{h}$ by $\pi_{h}^{\ast}u:=\pi_{h}(u^{\ast})$. For some fixed Lipschitz-domain $\Omega^{\ast}$ satisfying $\Omega^{\ast}_{h}\subseteq\Omega^{\ast}$ for $h\lesssim 1$, then the following interpolation estimates hold for functions $v\in H^{r}(\Omega^{\ast})$ and its fictitious domain variant

[TABLE]

owing to the boundedness of the extension operator (4.11), with $s,t\in\mathbb{N}$, $m=\min\{r,k+1\}$, $k$ the interpolation order of $\mathcal{X}_{h}$ and $\omega(T)$ the set of elements in $\mathcal{T}_{h}$ sharing at least one vertex with $T$. Throughout the analysis we write $\boldsymbol{\pi}_{h}^{\ast}$ for the Clément interpolant of vector-valued functions $\boldsymbol{v}$ and $\Pi_{h}^{\ast}$ for functions in a product space.

5 Continuity and Stability Estimates

In this section, we establish stability properties of the proposed cut finite element method (3.52). Our presentation is closely related to the analysis presented in the preceding work by Massing et al. [2] and mainly differs in the analysis of the boundary terms related to the weak enforcement of general Navier boundary conditions. The major stability result relies on a modified coercivity estimate for the total bilinear form $A_{h}+S_{h}+G_{h}$. Since fluid instabilities and the extension of the finite element method to unfitted meshes can be treated almost equivalently to Massing et al. [2], we only recall the required statements from the latter publication without presenting proofs.

After deriving continuity estimates for parts of the stabilized bilinear form in Lemma 5.1, we start by proving a coercivity estimate for our stabilized formulation w.r.t. the semi-norm $|U_{h}|_{h}$ on $\mathcal{W}_{h}$ in Lemma 5.2. We then recall how to recover required control over three additional (semi-)norms

[TABLE]

in Lemmas 5.3, 5.4 and 5.6 with the help of CIP and ghost-penalty stabilizations as introduced in Sections 3.3 and 3.4 using techniques provided in [57, 2]. By combining the above mentioned estimates, a global inf-sup stability on $\mathcal{W}_{h}$ is shown in Theorem 5.7.

Lemma 5.1** (Continuity Estimates).**

For an arbitrary choice of functions $\boldsymbol{u},\boldsymbol{v}\in[H^{2}(\Omega)]^{d}$, $p\in H^{1}(\Omega)$ and discrete functions $\boldsymbol{u}_{h},\boldsymbol{v}_{h}\in\mathcal{V}_{h}$, $p_{h}\in\mathcal{Q}_{h}$ the following continuity estimates hold provided that $\gamma^{n},\gamma^{t}<\infty$

[TABLE]

Proof. We start by proving the continuity estimate for $a_{h}$ with neglected advective bulk term as in (5.2). The proof is straightforward as it follows directly from applying the Cauchy-Schwarz inequality for each term. For the bulk terms we have

[TABLE]

It remains to estimate the boundary terms. By analogy, the estimate for all symmetric viscous and advective inflow Nitsche penalty terms follows directly from applying Cauchy-Schwarz inequality and the definition of $|||\cdot|||_{\ast}$. The non-symmetric boundary terms in $a_{h}$ can be estimated as

[TABLE]

All remaining boundary terms in $a_{h}$ can be estimated analogously, which yields

[TABLE]

and as a consequence (5.2). Estimate (5.3) can easily be deduced by integrating the advective bulk term by parts and applying the Cauchy-Schwarz inequality for the resulting boundary term

[TABLE]

Setting $\boldsymbol{v}=\boldsymbol{u}=\boldsymbol{0}$ in (5.14)–(5.16) and applying relation (3.64), i.e. $|\mspace{-1.0mu}|\mspace{-1.0mu}|\boldsymbol{v}_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{\ast}\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|\boldsymbol{v}_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}$, gives

[TABLE]

which leads to (5.4) and (5.5).

Similar estimates can be established for $b_{h}$ by applying the Cauchy-Schwarz inequality

[TABLE]

which proves (5.6). Choosing $\boldsymbol{v}=\boldsymbol{0}$ and using (3.64) yields estimate (5.7).

The next lemma shows how the Nitsche-type boundary terms, which have been consistently added to enforce the general Navier boundary conditions, guarantee stability of the bilinear form with respect to a semi-norm on $\mathcal{W}_{h}$ and how ghost-penalty terms extend stability to the enlarged domain $\Omega_{h}^{\ast}$.

Lemma 5.2** (Coercivity estimate).**

The stabilized cut finite element formulation is coercive, i.e.

[TABLE]

holds whenever the CIP and ghost penalty stability parameters $\gamma_{\nu},\gamma_{\sigma},\gamma_{\beta},\gamma_{u},\gamma_{p}$ are chosen large enough. For $\zeta_{u}=1$, the Nitsche penalty parameters need to be chosen small enough, i.e. $0<\gamma^{n},\gamma^{t}\leqslant C<\infty$, where $C$ depends on the shape and polynomial order of the elements (see inverse estimate (4.1)), however, it is independent of how the boundary intersects the element. By contrast, for $\zeta_{u}=-1$, the bilinear form is coercive for any choice $0<\gamma^{n},\gamma^{t}<\infty$. Note that $\operatorname{meas}{(\Gamma^{\varepsilon})}>0$ is assumed.

Proof. Starting from the definition of $A_{h}$, see (3.8), we have

[TABLE]

Integration by parts for the advective term together with continuity of $\boldsymbol{\beta}$ yields

[TABLE]

Using the assumption $\nabla\cdot\boldsymbol{\beta}=0$, with the help of the advective inflow control imposed on $\Gamma_{\mathrm{in}}\subseteq\Gamma_{D}\subseteq\Gamma$, the boundary control can be recovered at the entire boundary $\Gamma$ as

[TABLE]

For the viscous boundary terms, the following estimates can be established. By applying a $\delta$-scaled Cauchy-Schwarz inequality, a trace inequality (4.1) and the viscous norm equivalence from Proposition 3.4 for cut meshes, the following estimate can be established

[TABLE]

By analogy, the following boundary term from (5.23) can be estimated as

[TABLE]

The last term in (5.23) can be bounded by applying a trace inequality, such that

[TABLE]

From the Corollary 4.3 a part of the remaining viscous terms in (5.23) can be estimated as

[TABLE]

and by combining the previous inequalities (5.23), (5.25), (5.27), (5.28), (5.29) with (4.1), the following estimate can be done

[TABLE]

At this stage it is trivial to see that if an adjoint-inconsistent approach, i.e. $\zeta_{u}=-1$, is chosen, stability is guaranteed for any positive choice of $\gamma^{n},\gamma^{t}<\infty$. More delicate is to show stability in the adjoint-consistent case, i.e. $\zeta_{u}=1$, which is considered from here on. From the last two rows we have an upper bound with $\delta_{1},\delta_{2}\leqslant c<1$. The second row can be rewritten as

[TABLE]

Choosing, for instance, $\delta_{1}=\delta_{2}=1/2$, stability of the adjoint-consistent Nitsche-type method can be guaranteed by bounding the stabilization parameters $\gamma^{n},\gamma^{t}$ by $\gamma^{n}\leqslant\delta_{1}/8=1/16$ and $\gamma^{t}\leqslant\delta_{2}/8=1/16$. Note that the parameters also need to be lower bounded, i.e. $\gamma^{n},\gamma^{t}>0$, not to lose contributions of the bulk measures contained in the energy-type semi-norm $|U_{h}|_{h}$ and to prevent ill-conditioning. Finally, the claim follows by adding all other stabilization terms in $S_{h}$ and $G_{h}$, not already included in the analysis.

The next three lemmas are recalled from the work by Massing et al. [2] and show how the CIP, the ghost-penalty and the Nitsche-related stabilization terms can be used to recover the missing velocity and pressure semi-norm parts which are included in $|\mspace{-1.0mu}|\mspace{-1.0mu}|V_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}$.

First, we start by recovering control over the incompressibility constraint with the help of the stabilization terms $s_{u}$ and $g_{u}$. This is required for low Reynolds-numbers when the viscous $H^{1}$-semi-norm control $\|\nu^{\frac{1}{2}}{\boldsymbol{\epsilon}}(\boldsymbol{u}_{h})\|_{\Omega}$ vanishes.

Lemma 5.3**.**

There is a constant $c_{1}>0$ such that for each $\boldsymbol{u}_{h}\in\mathcal{V}_{h}$ there exists a $q_{h}\in\mathcal{Q}_{h}$ satisfying

[TABLE]

and the stability estimate

[TABLE]

with $\phi\in\{\phi_{u},\Phi^{-1}\}$ whenever the CIP and ghost-penalty parameter $\gamma_{u}$ is chosen large.

Proof. A detailed proof can be found in Massing et al. [2] (cf. Lemma 6.3).

The next lemma shows how additional control over a mixed semi-norm of the form $\|\phi_{\beta}^{\frac{1}{2}}(\boldsymbol{\beta}\cdot\nabla\boldsymbol{u}_{h}+\nabla p_{h})\|_{\Omega}^{2}$ can be recovered with the help of the CIP stabilizations $s_{\beta},s_{p}$ and the ghost penalty stabilizations.

Lemma 5.4**.**

There is a constant $c_{2}>0$ such that for each $U_{h}=(\boldsymbol{u}_{h},p_{h})\in\mathcal{W}_{h}$ we can construct a $\boldsymbol{v}_{h}\in V_{h}$ satisfying

[TABLE]

and the stability estimate

[TABLE]

whenever the stability parameters $\gamma_{\beta},\gamma_{p},\gamma_{\nu},\gamma_{\sigma}$ are chosen sufficiently large and the Nitsche penalty parameter $\gamma^{n}$ is chosen sufficiently small.

Proof. A detailed proof of this estimate has been presented in Massing et al. [2] (cf. Lemma 6.4).

Next, we recall two estimates from [2] which will be useful in deriving the final stability estimate in Theorem 5.7.

Lemma 5.5**.**

Let $\boldsymbol{u}_{h},\boldsymbol{v}_{h}\in\mathcal{V}_{h}$, then the following estimates hold

[TABLE]

[TABLE]

Proof. A detailed proof can be found in Massing et al. [2] (cf. Lemma 6.5).

Finally, a stabilized inf-sup condition for the operator $b_{h}(p_{h},\boldsymbol{v}_{h})$ holds, provided that the pressure stabilization operator $s_{p}$ is added, which gives the desired $L^{2}$-pressure norm control.

Lemma 5.6**.**

There is a constant $c_{3}>0$ such that for each pressure $p_{h}\in\mathcal{Q}_{h}$ there exists a velocity field $\boldsymbol{v}_{h}\in\mathcal{V}_{h}$ satisfying

[TABLE]

and the stability estimate

[TABLE]

whenever the stability parameters $\gamma_{\nu},\gamma_{\sigma},\gamma_{\beta},\gamma_{u},\gamma_{p}$ are chosen large enough.

Proof. A detailed proof of this modified inf-sup condition has been presented in Massing et al. [2] (cf. Lemma 6.6).

Collecting the previous lemmas, the final inf-sup stability estimate of our stabilized cut finite element method $A_{h}+S_{h}+G_{h}$, see (3.52), can be stated with respect to the energy norm $|\mspace{-1.0mu}|\mspace{-1.0mu}|U_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}^{2}$. The subsequent theorem ensures existence and uniqueness of a discrete velocity and pressure solution.

Theorem 5.7**.**

Let $U_{h}=(\boldsymbol{u}_{h},p_{h})\in\mathcal{W}_{h}$. Then, under the assumptions of Lemma 5.2 on the Nitsche penalty parameter $\gamma^{t}$ depending on $\zeta_{u}\in\{-1,1\}$, on $\gamma^{n}$ being sufficiently small as stated in Lemmas 5.2 and 5.4 and that CIP and ghost penalty stability parameters $\gamma_{\nu},\gamma_{\sigma},\gamma_{\beta},\gamma_{u},\gamma_{p}$ are chosen large enough (see Lemmas 5.2–5.6), the cut finite element method (3.52) is inf-sup stable

[TABLE]

Note that the hidden stability constant is independent of the mesh size $h$, the slip-length coefficient $\varepsilon\in[o,\infty]$ and in case of unfitted meshes independent of the position of the boundary relative to the background mesh.

Proof. Since the proof follows the procedure proposed by Massing et al. [2], the major steps are only sketched in the following. For a given $U_{h}\in\mathcal{W}_{h}$ we construct a test function $V_{h}\in\mathcal{W}_{h}$ based on estimates derived in Lemmas 5.2, 5.3, 5.4 and 5.6.

Step 1. Choosing the test function $V_{h}^{1}:=(\boldsymbol{0},q_{h}^{1})$ with $q_{h}^{1}$ from Lemma 5.3 to recover divergence control and utilizing its stability bounds (5.33) and (5.34) we have

[TABLE]

for any positive $\gamma^{n}\leqslant C<\infty$. This shows the need for the boundary-normal Nitsche penalty term at $\Gamma$ for convective and reactive dominant flows as well as the need for the pressure and divergence CIP and ghost penalty terms.

Step 2. Inserting $V_{h}^{2}:=(\boldsymbol{v}_{h}^{2},0)$ with $\boldsymbol{v}_{h}^{2}$, from Lemma 5.4 into the formulation $A_{h}+S_{h}+G_{h}$, followed by applying the continuity estimate (5.4) for $a_{h}$ and the stability bounds (5.36) and (5.37) yields the following

[TABLE]

Step 3. The $L^{2}$-pressure norm term can be constructed by testing with $V_{h}^{3}:=(\boldsymbol{v}_{h}^{3},0)$, where $\boldsymbol{v}_{h}^{3}$ is now chosen as in Lemma 5.6. Utilizing the continuity estimate (5.5) for $a_{h}$ and making use of the estimates (5.39), (5.40), and the stability bound (5.42) for the chosen test function allows us to deduce that

[TABLE]

Step 4. From the coercivity estimate in Lemma 5.2 we obtain the positive semi-norm term $|U_{h}|_{h}$ with $V_{h}^{4}:=U_{h}$

[TABLE]

Step 5. It remains to combine Step 1–4 by choosing $\delta$ sufficiently small and defining the final test function

[TABLE]

for a given $U_{h}$. Choosing $\eta>0$ sufficiently small for some $2\eta\leavevmode\nobreak\ \sim(C_{1}(\delta)+C_{2}(\delta)+C_{3}(\delta))^{-1}$ allows us to gain control over all desired norm parts in $|\mspace{-1.0mu}|\mspace{-1.0mu}|U_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}$ and at the same time to absorb the defective $|U_{h}|_{h}$-contribution, which stem from testing with $V_{h}^{1},V_{h}^{2},V_{h}^{3}$. Consequently,

[TABLE]

The inf-sup stability estimate can be concluded as

[TABLE]

after dividing by $|\mspace{-1.0mu}|\mspace{-1.0mu}|V_{h}^{5}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}$ and choosing the supremum over $V_{h}\in\mathcal{W}_{h}\backslash\{0\}$. To conclude the proof, it remains to prove that $|\mspace{-1.0mu}|\mspace{-1.0mu}|V_{h}^{5}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|U_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}$. Note that

[TABLE]

Thanks to the stability estimate (5.37) and norm definitions (3.59) and (3.60) it holds that

[TABLE]

Similarly, the stability bound (5.34) for $q_{h}^{1}$ implies that $|\mspace{-1.0mu}|\mspace{-1.0mu}|V_{h}^{1}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}^{2}=\Phi\|q_{h}^{1}\|_{\Omega}^{2}+|q_{h}^{1}|_{h}^{2}\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|U_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}^{2}$ and from (5.42) it is obtained that $|\mspace{-1.0mu}|\mspace{-1.0mu}|V_{h}^{3}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}^{2}=|\mspace{-1.0mu}|\mspace{-1.0mu}|\boldsymbol{v}_{h}^{3}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}^{2}+(1+\omega_{h})^{-1}\|\phi_{\beta}^{\frac{1}{2}}\boldsymbol{\beta}\cdot\nabla\boldsymbol{v}_{h}^{3}\|_{\Omega}^{2}+\|\phi_{u}^{\frac{1}{2}}\nabla\cdot\boldsymbol{v}_{h}^{3}\|_{\Omega}^{2}\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|U_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}^{2}$. As a result it holds

[TABLE]

which concludes the proof.

6 A Priori Error Estimates

This section is devoted to the a priori error analysis of our cut finite element method. We derive an energy-norm error estimate for the discrete velocity and pressure solution in Theorem 6.4 and use a standard duality technique to establish optimal error convergence in the velocity $L^{2}$-norm for the adjoint-consistent formulation with dominating viscous flow effects in Theorem 6.5.

6.1 Consistency and Interpolation Error Estimates

We start by showing that the discrete formulation (3.52) satisfies a weakened form of the Galerkin orthogonality.

Lemma 6.1** (Weakened Galerkin Orthogonality).**

Assume that the solution $U=(\boldsymbol{u},p)$ of the formulation (2.1)–(2.4) is sufficiently regular, i.e. in $[H^{2}(\Omega)]^{d}\times H^{1}(\Omega)$, and let $U_{h}=(\boldsymbol{u}_{h},p_{h})\in\mathcal{V}_{h}\times\mathcal{Q}_{h}$ be the cut finite element solution to the discrete weak formulation (3.52). Then, the error $U-U_{h}$ satisfies a weak Galerkin orthogonality property

[TABLE]

Proof. The proof follows standard techniques. Multiplying the problem formulation (2.1)–(2.2) with test functions $V_{h}=(\boldsymbol{v}_{h},q_{h})\in\mathcal{V}_{h}\times\mathcal{Q}_{h}$, integrating over $\Omega$, performing integration by parts and using the fact that all additional Nitsche-related terms vanish for $U$ satisfying the boundary condition, yields

[TABLE]

where the last equation holds for the discrete solution $U_{h}\in\mathcal{V}_{h}\times\mathcal{Q}_{h}$. The combination $\textbf{I}-\textbf{II}$ yields the claim.

The subsequent lemma recalls the well known weak consistency property of continuous interior penalty and ghost penalty stabilization operators $S_{h},G_{h}$, and as such they do not deteriorate the optimality of our cut finite element scheme.

Lemma 6.2** (Weak Consistency).**

For all functions $(\boldsymbol{u},p)\in[H^{r}(\Omega)]^{d}\times H^{s}(\Omega)$ there holds

[TABLE]

where $k$ is the polynomial degree of the respective Clément interpolants $\Pi^{\ast}_{h}U=(\boldsymbol{\pi}^{\ast}_{h}\boldsymbol{u},\pi^{\ast}_{h}p)$ for the velocity and pressure and $r_{u}:=\min\{r,k+1\}$ and $s_{p}:=\min\{s,k+1\}$.

Proof. A detailed proof of the weak consistency of the proposed continuous interior penalty and ghost penalty operators $S_{h}$ and $G_{h}$ has been given in [2] (see Lemma 7.2) and the references therein.

The next lemma ensures that the interpolation error between a continuous solution and its Clément interpolation converges with optimal rates.

Lemma 6.3** (Interpolation Estimates).**

Assume that $(\boldsymbol{u},p)\in[H^{r}(\Omega)]^{d}\times H^{s}(\Omega)$ and let $r_{u}:=\min\{r,k+1\}\geqslant 2$, $s_{p}:=\min\{s,k+1\}\geqslant 1$ where $k$ is the polynomial degree of the approximation spaces for the velocity and pressure. Then

[TABLE]

with $\phi\in\{\phi_{u},\Phi^{-1}\}$ in the pressure estimate.

Proof. The proof follows the techniques proposed in [2] (see Lemma 7.4). By applying the interpolation estimate (4.13), the viscous, the reactive and the divergence bulk error measures can be estimated as

[TABLE]

Considering the normal and tangential boundary semi-norms with $\gamma^{t},\gamma^{n}\geqslant c>0$ and $\varepsilon\geqslant 0$, the claim follows by combining the interpolation estimate (4.13), the trace inequality (4.2) and the definition of the scaling function $\phi_{u}$ such that

[TABLE]

The remaining boundary terms can be estimated in a similar fashion,

[TABLE]

Collecting all estimates and noting that $\Phi\lesssim\phi_{u}^{-1}$ yields the claim.

6.2 A Priori Error Estimates

Subsequently, the main a priori estimate for the velocity and pressure errors w.r.t a natural energy norm is stated.

Theorem 6.4** (Energy norm error estimate).**

Assume that the continuous solution of the Oseen problem (2.1)–(2.4) resides in $U=(\boldsymbol{u},p)\in[H^{r}(\Omega)]^{d}\times H^{s}(\Omega)$ and let $U_{h}=(\boldsymbol{u}_{h},p_{h})\in\mathcal{V}_{h}\times\mathcal{Q}_{h}$ be the discrete solution of problem (3.52). Let the energy type norm be defined as in (3.63), then

[TABLE]

where $r_{u}:=\min\{r,k+1\}\geqslant 2$ and $s_{p}:=\min\{s,k+1\}\geqslant 1$. Note that the hidden constants are independent of $h$, are bounded with respect to the slip length coefficient $\varepsilon\in[0,\infty]$ and, owing to the ghost penalty stabilization terms $G_{h}$, independent of how the boundary intersects the mesh $\mathcal{T}_{h}$. The (hidden) scaling functions $\omega_{h},\Phi$ are as defined in (3.49).

Proof. We first split the total discretization error into a discrete error and an interpolation part by applying (3.66) and the norm definitions in Section 3.5

[TABLE]

As the term $|\mspace{-1.0mu}|\mspace{-1.0mu}|U-\Pi^{\ast}_{h}U|\mspace{-1.0mu}|\mspace{-1.0mu}|_{\ast}$ is readily estimated from the interpolation estimates (6.5) and (6.6), only the discrete error $|\mspace{-1.0mu}|\mspace{-1.0mu}|\Pi^{\ast}_{h}U-U_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}$ is considered from here on. From the inf-sup condition (5.43) there exists a $|\mspace{-1.0mu}|\mspace{-1.0mu}|V_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}=1$ such that

[TABLE]

where the last step follows from applying the weak Galerkin orthogonality (6.1). By applying a Cauchy-Schwarz inequality on the stabilization terms $S_{h}(\Pi^{\ast}_{h}U,V_{h})$ and $G_{h}(\Pi^{\ast}_{h}U,V_{h})$ and utilizing the results from Lemma 6.2 the following estimate can be made

[TABLE]

The term $A_{h}(\Pi^{\ast}_{h}U-U,V_{h})$ can be estimated by integrating $b_{h}(q_{h},\boldsymbol{\pi}^{\ast}_{h}\boldsymbol{u}-\boldsymbol{u})$ by parts and applying the continuity estimates (5.3) and (5.7) to $a_{h}(\boldsymbol{\pi}^{\ast}_{h}\boldsymbol{u}-\boldsymbol{u},\boldsymbol{v}_{h})$ and $b_{h}(\pi^{\ast}_{h}p-p,\boldsymbol{v}_{h})$ respectively, such that

[TABLE]

Thanks to the interpolation estimate Lemma 6.3 and keeping in mind that $|\mspace{-1.0mu}|\mspace{-1.0mu}|\boldsymbol{v}_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}+\|\phi_{u}^{\frac{1}{2}}\nabla\cdot\boldsymbol{v}_{h}\|_{\Omega}\lesssim|\mspace{-1.0mu}|\mspace{-1.0mu}|V_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{h}=1$, it only remains to estimate the last term. Its estimate follows readily from applying a Cauchy-Schwarz inequality, using (4.12) and the definition of $\phi_{\beta}$ to arrive at

[TABLE]

which concludes the proof of the a priori error estimate (6.18).

6.3 $L^{2}$-Optimal Estimate for Flows with Large Viscosity

We now proceed to deduce an optimal $L^{2}$-error estimate for the velocity of an adjoint-consistent ($\zeta_{u}=1$) formulation for flows dominated by viscous forces, i.e. from here on we assume

[TABLE]

The proof follows the standard Aubin–Nitsche duality argument, see e.g. [58, 59], and requires the established optimal energy-type error estimate from Theorem 6.4. A similar estimate for the Oseen problem has been derived for boundary-fitted meshes by Burman et al. [38]. We introduce the following dual problem to (2.1)–(2.4): find adjoint velocity and pressure $(\boldsymbol{w},r)$ such that

[TABLE]

for which we assume additional elliptic regularity and that the solution belongs to $(\boldsymbol{w},r)\in[H^{2}(\Omega)]^{d}\times H^{1}(\Omega)$ so that it satisfies

[TABLE]

provided that the boundaries are sufficiently smooth, see e.g. [41, 60, 61]. Note that due to the homogeneous boundary conditions in (6.31)–(6.32) the estimate (6.33) is independent of the boundary data. However, the dimensionless constant $C$ depends on the physical parameters $\nu,\sigma,\boldsymbol{\beta}$ and the domain $\Omega$. By choosing $\tilde{\boldsymbol{f}}=\boldsymbol{u}-\boldsymbol{u}_{h}\in L^{2}(\Omega)$ as the right hand side of the dual momentum equation, the desired error quantity bounds the dual solution $(\boldsymbol{w},r)$ in (6.33). Note that the dual advective velocity is set as the negative advective velocity field of the primal problem (2.1)–(2.4). As a result, inflow and outflow parts of the boundary of primal and dual problems swap, respectively. Similarly, the dual pressure solution $r$ changes the sign compared to the pressure solution $p$ of the primal problem. For further explanations on the dual problem, see e.g. the textbook [61]. Furthermore, it is assumed that $\boldsymbol{\beta}\cdot\boldsymbol{n}=0$ on $\Gamma\setminus\Gamma_{\mathrm{D}}$, which simplifies (6.32) in the following.

Theorem 6.5** (Velocity $L^{2}$-error estimate).**

Let $U=(\boldsymbol{u},p)\in[H^{r}(\Omega)]^{d}\times H^{s}(\Omega)$ be the continuous solution to the Oseen problem (2.1)–(2.4) and $U_{h}=(\boldsymbol{u}_{h},p_{h})\in\mathcal{V}_{h}\times\mathcal{Q}_{h}$ be the discrete solution of problem (3.52). Under previously specified assumptions, for an adjoint-consistent Nitsche-type formulation ($\zeta_{u}=1$) we have

[TABLE]

provided that viscous forces dominate the flow as assumed in (6.28) with $r_{u}:=\min\{r,k+1\}$ and $s_{p}:=\min\{s,k+1\}$. Note that the hidden constant is independent of $h$ and independent of how the boundary intersects the mesh $\mathcal{T}_{h}$. However, the constant depends on the physical parameters.

Proof. We need to estimate the desired velocity $L^{2}$-error $\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0,\Omega}$. For this purpose, we choose $\tilde{\boldsymbol{f}}=\boldsymbol{u}-\boldsymbol{u}_{h}$ as right hand side of the dual momentum equation (6.29) and multiply (6.29)–(6.30) with test functions $\boldsymbol{v}:=\boldsymbol{u}-\boldsymbol{u}_{h}$ and $-q:=-(p-p_{h})$, respectively. After integrating by parts, using $\nabla\cdot\boldsymbol{\beta}=0$ and the relation

[TABLE]

as $\boldsymbol{w}=\bf 0$ on $\Gamma_{D}$ and $\boldsymbol{\beta}\cdot\boldsymbol{n}=0$ on $\Gamma\setminus\Gamma_{\mathrm{D}}$, we obtain

[TABLE]

Using the boundary conditions for the normal and tangential directions (6.31)–(6.32) with the assumption that $\boldsymbol{\beta}\cdot\boldsymbol{n}=0$ on $\Gamma\setminus\Gamma_{\mathrm{D}}$, which then reduce to $\boldsymbol{w}\cdot\boldsymbol{n}=0$ on $\Gamma$ and $(\varepsilon(2\nu{\boldsymbol{\epsilon}}(\boldsymbol{w})\boldsymbol{n})+\nu\boldsymbol{w})\boldsymbol{P}^{t}=\boldsymbol{0}$ on $\Gamma$, the following terms can be consistently added to (6.38)

[TABLE]

Continuing in (6.38), we note that $b(q,\boldsymbol{w})+I=b_{h}(q,\boldsymbol{w})$. The viscous term can be split into directional parts such that

[TABLE]

Collecting all terms from (6.38)–(6.43) and defining $W:=(\boldsymbol{w},r)$, the $L^{2}$-error can be expressed as

[TABLE]

i.e. it can be expressed in terms of the discrete bilinear operator (3.8) associated to the primal problem (2.1)–(2.4), where the sufficiently smooth solution $W$ of the dual problem now takes the role of the test function.

Note that under the assumption of dominating viscous forces (6.28), it holds for the advective term occurring in the continuity estimate of $a_{h}$ in (5.2) that

[TABLE]

where the interpolation estimate (4.13) for the Clément interpolant was used in the last step. The $L^{2}$-velocity error can then be further estimated by using the weak Galerkin orthogonality from Lemma 6.1 with $\Pi^{\ast}_{h}W\in\mathcal{V}_{h}\times\mathcal{Q}_{h}$, the continuity of $a_{h},b_{h}$ provided in Lemma 5.1 (equations (5.2) and (5.6)) and by applying a Cauchy Schwarz inequality to the remaining stabilization operators

[TABLE]

In line (6.52) we use the energy-norm a priori error estimate from Theorem 6.4 in (6.18) and (6.20) for $(\boldsymbol{u},p)\in[H^{r_{u}}(\Omega)]^{d}\times H^{s_{p}}(\Omega)$ under the assumption of dominant viscous effects (6.28). Note that for $\phi\in\{\phi_{u},\Phi^{-1}\}$ it holds that $\nu/\phi\lesssim 1$. Thanks to the interpolation estimate Lemma 6.3, in (6.52) we gain the desired additional power of $h$ for $(\boldsymbol{w},r)\in[H^{2}(\Omega)]^{d}\times H^{1}(\Omega)$ by estimating

[TABLE]

together with the boundedness (6.33) of the solution to the dual problem (6.29)–(6.32). Note that for $S_{h}$ and $G_{h}$ the weak consistency estimates from Lemma 6.2 hold. Finally, the claim follows after dividing by $\|\boldsymbol{u}-\boldsymbol{u}_{h}\|_{0,\Omega}$ in (6.54).

7 Numerical Example

To ensure the validity of the proposed method and the theoretical results presented, a numerical example is conducted. The error estimates obtained in the a priori error analysis, as summarized in the Theorems 6.4 and 6.5, are validated by the results of a 2D box-flow case.

The stabilization parameters used for the following simulations are taken by large from Massing et al. [2] and Schott and Wall [18] and are repeated here for completeness. The CIP-stabilization terms (3.30)–(3.32) are set as $\gamma_{\beta}=\gamma_{p}=0.01$ and $\gamma_{u}=0.05\gamma_{\beta}$, as suggested in [62]. For the convective, pressure and velocity ghost penalty stabilizations (3.40)–(3.42) the same parameters as for the CIP-stabilization are used. In the case of the viscous (3.39) and (pseudo-) reactive (3.38) ghost-penalty a value of $\gamma_{\nu}=0.05$ and $\gamma_{\sigma}=0.005$ is prescribed, respectively. The second order terms of the ghost-penalties (i.e. for $j=2$ in the equations (3.38),(3.39) and (3.42) and $j=1$ in (3.40) and (3.41)) are scaled by an extra $0.05$, as without it too strong enforcement of these terms were observed, ruining the solution. Higher order terms of the ghost-penalties, i.e. $j>2$ are neglected as their influence on the solution and stability is negligible for linear and quadratic elements. Note that for simulations with higher-order approximations, i.e. $k>1$, the simplified variant $\overline{g}_{\beta}$ (3.51) for the convective and divergence ghost penalty terms $g_{\beta},g_{u}$ (see (3.40) and (3.41)) and the related continuous interior penalty stabilizations $s_{\beta},s_{u}$ (see (3.30) and (3.31)) are used.

The different flow regimes appearing in $\phi_{u},\phi_{\beta},\phi_{p}$ (see (3.20) and (3.33)) are weighted as $\nu+c_{u}(\|\boldsymbol{\beta}\|_{0,\infty,T}h)+c_{\sigma}(\sigma h^{2})$ with $c_{u}=1/6$ and $c_{\sigma}=1/12$ as suggested in [21]. If nothing else is mentioned, then the simulations are conducted with an adjoint-consistent Nitsche’s method (i.e. $\zeta_{u}=1$) and with penalty parameters of $1/\gamma^{n}=1/\gamma^{t}=10.0$. Furthermore, as there is a need to integrate cut elements of non-regular shapes, standard integration rules can not be applied. To overcome this issue, the integration rules proposed by Sudhakar et al. [63] are used.

All simulations presented in this publication have been performed using the parallel finite element software environment “Bavarian Advanced Computational Initiative” (BACI), see [64].

7.1 Problem Setup – 2D Box Flow

This numerical example is inspired from the example done for Stokes flow in [15], but made somewhat more complex and extended to the Oseen equations. We use the technique of manufactured solution to create the example. A divergence-free velocity field ($\nabla\cdot\boldsymbol{u}=0$) and a suitable pressure field are chosen which, when put into the Oseen equations with $\boldsymbol{\beta}=\boldsymbol{u}$, generate an associated volumetric body force $\boldsymbol{f}$. By providing the volumetric body force and appropriate boundary conditions the solution is known and error studies can be carried out for the example. The domain $\Omega$ used here is a square $\{(x,y),-1<x<1,-1<y<1\}$ and the following choice for the solution field is

[TABLE]

as visualized in Figure 7.1. The geometry here is chosen such that we have no geometric approximation error by the meshes, thus avoiding issues stemming from the Babuška paradox [65]. It can easily be seen that the chosen analytical velocity field ($\boldsymbol{u}=(u_{1},u_{2})$) satisfies $\nabla\cdot\boldsymbol{u}=0$. The volume force field $\boldsymbol{f}$ is chosen in accordance with the given fields (7.1)–(7.3) such that the Oseen equations (2.1)–(2.2) are satisfied. In the discrete case, the force field $\boldsymbol{f}_{h}$ and the advective velocity $\boldsymbol{\beta}_{h}$ are given by the nodal interpolations of their continuous counterparts.

The boundary conditions applied to this problem are the same as introduced in the Oseen problem (2.3)–(2.4). The velocity at the boundary is prescribed as $\boldsymbol{g}=\boldsymbol{u}$, and the traction at the boundary as $\boldsymbol{h}=2\nu{\boldsymbol{\epsilon}}(\boldsymbol{u})$. It should be mentioned that $\boldsymbol{u}\boldsymbol{P}^{n}=\boldsymbol{0}$ on $\Gamma$, which satisfies the conditions assumed throughout the analysis. Furthermore, with this choice of $\boldsymbol{g}$ and $\boldsymbol{h}$, the prescribed solution is independent of the choice of slip length $\varepsilon$. The boundary condition is imposed on $\Gamma$ by the method (3.52) introduced in Section 3.2.

The background mesh $\widehat{\mathcal{T}}_{h}$ covers a rectangular domain $[-1.6,1.6]^{2}$ which is rotated by the angle $\theta=0.25\pi$ around the origin to ensure a non-trivial computational mesh $\mathcal{T}_{h}$, see Figure 7.2. It consists of $N\times N$ equally sized square elements equipped with either bi-linear $\mathcal{Q}^{1}$ or bi-quadratic $\mathcal{Q}^{2}$ equal-order approximations for velocity and pressure, i.e. $\mathcal{V}^{k}_{h}\times\mathcal{Q}^{k}_{h}$ where $k\in\{1,2\}$. As a consequence the element length is $h=3.2/N$. The parameters for the Oseen problem are chosen as $\nu=1.0$ and $\sigma=1.0$. In the case of convection dominated flow (i.e. for small $\nu$), the tangential components disappear as they scale with $\nu$, and effectively the same formulation as presented and studied in [2] is retained. For this reason, we focus on investigating viscous flows here. As the pressure is only defined up to a constant, the constant pressure mode is filtered out by imposing $\int_{\Omega}p_{h}dx=0$ during the solution process.

7.2 Mesh Refinement Study

To verify the results from the a priori error analysis, mesh refinement studies are conducted for both linear and quadratic interpolations. In the case of the linear $\mathcal{Q}^{1}$ elements, a series of mesh sizes are chosen with $N\in[8;512]$ and for the quadratic $\mathcal{Q}^{2}$ case $N\in[8;224]$. To demonstrate the robustness to the choice of slip length, the convergence studies are conducted with a choice of three different $\varepsilon\in[10^{-10};1.0;10^{10}]$.

In Figure 7.3 the $L^{2}$-errors of $\boldsymbol{u}_{h}$, $\nabla\boldsymbol{u}_{h}$ and $p_{h}$ are presented for the linear $\mathcal{Q}^{1}$ approximations evaluated in $\Omega$ and on $\Gamma$. The convergence plots verify the optimal convergence $\mathcal{O}(h^{2})$ for the velocity $L^{2}$-error (see Theorem 6.5) for all choices of $\varepsilon$. Furthermore, the errors for the different choices of the various slip lengths remain of comparable size throughout the mesh refinement. The convergence of the velocity gradient and pressure, as seen from the error at the boundary, converges with at least first order, in accordance with theory stated in Theorem 6.4.

For the quadratic $\mathcal{Q}^{2}$ interpolation case the same error norms as for the linear case are studied, as shown in Figure 7.4. Similar conclusions can be drawn as for the linear case. An optimal convergence $\mathcal{O}(h^{3})$ for the velocity $L^{2}$-error and an order of at least $\mathcal{O}(h^{2})$ for the velocity gradient and pressure are observed for all $\varepsilon$, in accordance with theory. This verifies the applicability of the method for higher order elements as well. Also notable is that the quadratically interpolated elements demonstrate less sensitivity to the cut position than the linear elements.

7.3 Nitsche Stabilization Parameter Study

To observe the behavior of the method for different choices of the Nitsche parameter, a sensitivity study is conducted. For the adjoint consistent case the study is done with $\mathcal{Q}^{1}$ elements and for the adjoint inconsistent case for both $\mathcal{Q}^{1}$ and $\mathcal{Q}^{2}$ elements. In the case of linear elements a $512\times 512$ mesh is used and for the quadratic case a $224\times 224$ mesh. Both the tangential and normal stabilization parameters are varied at the same time, i.e. $\gamma=\gamma^{t}=\gamma^{n}$, and studies are conducted for three different choices of slip length $\varepsilon=\{10^{-10};1.0;10^{10}\}$.

In Figure 7.5 the results for the adjoint-consistent formulation ($\zeta_{u}=1$) can be seen. As is verified from the inf-sup stability analysis in Section 5, the error diverges as $1/\gamma$ becomes too small, and in this case the stability limit is reached at around $1/\gamma\approx 4.0$. Notable from the figure is the different behavior of the velocity error at the interface for varying $\varepsilon$. In the case of $\varepsilon=10^{-10}$ the error is slowly decreasing for a choice of larger $1/\gamma$, whereas for the cases with $\varepsilon=1.0$ and $\varepsilon=10^{10}$ the error is strongly increasing. This discrepancy can be explained by the fact that in the latter cases the penalty effect on the normal constraint is way stronger than for the tangential part and thus allows for worsening of the imposition of the tangential condition while better enforcing the normal constraint for a large penalty parameter $1/\gamma$. In contrast, in the former case the velocity is enforced in both normal and tangential direction, thus imposing the velocity well at the boundary. However, note that at the same time larger $1/\gamma$ influences the velocity gradient and pressure at the boundary negatively, as expected.

The introduced method also works for the adjoint-inconsistent case ($\zeta_{u}=-1$). This enables the choice of smaller values of $1/\gamma$, as predicted from theory. However, the choice of both a small $1/\gamma^{t}$ and a large $\varepsilon$ leads to conditioning issues, see Remark 3.2. This can be seen in Figures 7.6 and 7.7 where for both the linear and quadratic case with $\varepsilon=10^{10}$ small values of the stabilization parameter lead to the case where the used solver will not converge to the desired tolerance. This occurs for the linear case with $1/\gamma<10^{-3}$ and for the quadratic case with $1/\gamma<10^{-4}$. For the other choices of $\varepsilon$, the error remains stable for smaller choices of $1/\gamma$ confirming the results from theory which states that $1/\gamma>0$ is a sufficient choice for stability. In the quadratic case, see Figure 7.7, similar results as for the linear case are observed. However, the two observed minima for the pressure error are not present. This indicates that this phenomenon stems from the linear elements and is not an inherent property of the proposed method. Another observation is the fact that the error for the quadratic approximations remains stable longer before a large stabilization parameter $1/\gamma$ effects the solution negatively.

From the graphs $1/\gamma\approx 10$ gives a clear minimum for the pressure error for both the linear and quadratic approximations. However, in the quadratic case the velocity error is smaller for a choice of $1/\gamma\approx 35$. From these observations a good choice of the stabilization parameter should be around $1/\gamma\in(10,50)$. This proposed interval agrees well with the observed results in Schott and Wall [18].

7.4 Slip Length Sensitivity Study

To demonstrate the robustness of the method to the choice of slip length, errors are evaluated for a series of $\varepsilon\in[10^{-10},10^{10}]$ for linear $\mathcal{Q}^{1}$ elements. A comparison is done with the substitution method (3.6) explained in Section 3.2 to emphasize the advantages of our method. As can be seen in Figure 7.8, when imposing the general Navier condition by means of our Nitsche’s method the error remains almost constant and the difference between the errors of the limiting cases (i.e. $\varepsilon\rightarrow 0$ or $\varepsilon\rightarrow\infty$) are small. Even though the analytic solution is independent of $\varepsilon$ some difference between the limiting cases are expected for the numerical simulations. Furthermore, the advantage of the proposed method to the substitution method is clear from this Figure. The substitution method starts to produce noticeably larger errors for $\varepsilon<10^{-5}$, and for smaller choices of slip length consecutively worse results are observed up until our linear solver could not solve the system for $\varepsilon<10^{-8}$. These results are expected for the substitution method (3.6) as the conditioning of the system becomes increasingly bad when the slip length approaches zero, as discussed in Section 3.2. Furthermore, it is worth to be noted that in the limiting case of $\varepsilon\rightarrow\infty$ the errors between the two methods are of comparable size and the additional consistent terms (3.24) and (3.28) do not deteriorate the error.

8 Conclusions

In this work, a novel method for the Oseen problem with a general Navier slip boundary condition is introduced. This boundary condition is best explained as a Robin condition in the tangential plane and a Dirichlet condition in the normal direction of the boundary. The proposed method imposes the boundary condition weakly by means of a Nitsche’s method for both the tangential and normal parts of the condition. It remains well-posed and stable for all choices of slip lengths and for both low and high Reynolds numbers. Furthermore, the presented formulation is a first step for building a more general formulation for the imposition of boundary conditions for the Oseen and Navier-Stokes equations since in the tangential direction it allows for any linear combination of Dirichlet and Neumann conditions to be set within the same framework, which previously was not possible. This method is presented for the case of unfitted grids and an equal-order interpolated cut finite element method is used for the discretization. In the analysis, it is shown that the proposed formulation remains stable irrespective of where the boundary intersects the background mesh. The theoretical findings in this work remain valid also for the simpler boundary-fitted grid case and can thus readily be applied for this case as well.

To show the validity of the proposed formulation, inf-sup stability is shown. Furthermore, an a priori error analysis for an energy-type norm and for the $L^{2}$–error of the velocity are conducted for the adjoint-consistent formulation. A numerical example corroborates the findings, where optimal order of convergence was observed for the error norms of both linear and quadratic approximations. Also observed in a numerical example is that our proposed method performs better for small slip lengths compared to the classical substitution method, and for larger values of the slip length it performs just as well. This is in agreement with theory, as our method is well-posed in both the Dirichlet and Neumann limits for the tangential condition, whereas the classical substitution method is not defined in the Dirichlet limit, which coincides with the choice of a small slip length.

As the proposed method enables the possibility to use both no-slip and slip with the same method, it is a promising approach for a number of future applications. These include for instance fluid-structure-contact interaction as pure no-slip walls may lead to a not well-posed problem in the region where submerged bodies are close to each other. Another promising application currently investigated by the authors is the use of a localized Navier-slip condition in the modeling of contact line motion for multi-phase flows.

Acknowledgement

This work is supported by the International Graduate School of Science and Engineering (IGSSE) of the Technical University of Munich, Germany, under project 6.02 and the Erasmus Mundus Joint Doctorate SEED project (European Commission). The authors would also like to thank Antonio Huerta for the discussions had about the formulation.