Finite elements for Stokes with traction boundary conditions

Question

Suppose we are given the Stokes equations with Neumann conditions on part of the boundary:

$-\nabla\cdot\boldsymbol{\sigma} = \mathbf{f}, \quad \text{and} \quad \nabla\cdot \mathbf{u} = 0 \quad \text{in} \quad \Omega$

Here $\boldsymbol{\sigma} := \nu(\nabla \mathbf{u} + (\nabla \mathbf{u})^T) - p I$, and $\mathbf{u}=[u_1,\dots,u_d]^T$. Additionally, we assume the boundary conditions of Dirichlet and Neumann-type

$\mathbf{u} = 0 \quad \text{on} \quad \partial\Omega\backslash\Gamma, \quad \text{and} \quad \boldsymbol{\sigma}\cdot\mathbf{n} = \mathbf{t}\quad \text{on} \quad \Gamma$,

where $\mathbf{n}$ denotes the outward pointing normal vector. Assuming that $\nu$ is constant in the whole of $\Omega$, we can rewrite the momentum equation as

$-\nu\Delta\mathbf{u} + \nabla p = \mathbf{f}, \quad \text{and} \quad \nabla\cdot \mathbf{u} = 0 \quad \text{in} \quad \Omega$

by using the incompressibility $\nabla \cdot \mathbf{u} = 0$ and some elementary manipulations.

What puzzles me now is the following:

By the usual procedure, we obtain the following weak form: find $(\mathbf{u},p) \in \mathbf{H}^1_{0,\partial\Omega\backslash\Gamma}(\Omega)\times L^2(\Omega)$ such that

$\int_\Omega (\nu\nabla\mathbf{u} : \nabla\mathbf{v} - p\nabla\cdot\mathbf{v} - q\nabla\cdot\mathbf{u}) \,{\rm d}x - \int_\Gamma (\nu\nabla\mathbf{u}-pI)\cdot\mathbf{n}\cdot\mathbf{v}\,{\rm d}s = \int_\Omega \mathbf{f}\cdot\mathbf{v}\,{\rm d}x$

for all test functions $(\mathbf{v},q) \in \mathbf{H}^1_{0,\partial\Omega\backslash\Gamma}(\Omega)\times L^2(\Omega)$. At this point, one would normally replace the remaining boundary terms by the Neumann boundary conditions and put them to the right hand side. However, when considering the physical traction boundary conditions given above, we clearly have that

$\mathbf{t} = \boldsymbol{\sigma}\cdot\mathbf{n} = (\nu(\nabla \mathbf{u} + (\nabla \mathbf{u})^T)\cdot\mathbf{n} - p\mathbf{n} \ne \nu\nabla \mathbf{u}\cdot\mathbf{n} - p\mathbf{n}$

unless for some reason $(\nabla \mathbf{u})^T\cdot\mathbf{n} = 0$ of course. This now leads to the weak form: find $(\mathbf{u},p) \in \mathbf{H}^1_{0,\partial\Omega\backslash\Gamma}(\Omega)\times L^2(\Omega)$ such that

$\int_\Omega (\nu\nabla\mathbf{u} : \nabla\mathbf{v} - p\nabla\cdot\mathbf{v} - q\nabla\cdot\mathbf{u}) \,{\rm d}x + \int_\Gamma \nu(\nabla\mathbf{u})^T\cdot\mathbf{n}\cdot\mathbf{v}\,{\rm d}s = \int_\Omega \mathbf{f}\cdot\mathbf{v}\,{\rm d}x - \int_\Gamma \mathbf{t}\cdot\mathbf{v}\,{\rm d}s$.

for all test functions $(\mathbf{v},q) \in \mathbf{H}^1_{0,\partial\Omega\backslash\Gamma}(\Omega)\times L^2(\Omega)$. However, what I am really interested in is the problem in finite dimensional spaces $\mathbf{V}_h \times Q_h \subset \mathbf{H}^1_{0,\partial\Omega\backslash\Gamma}(\Omega)\times L^2(\Omega)$ for the velocity and pressure, respectively. Let us assume that these are LBB stable. Following the standard procedure for showing stability, one would still have to show that:

$a(\mathbf{u},\mathbf{v}) := \int_\Omega \nu\nabla\mathbf{u} : \nabla\mathbf{v} \,{\rm d}x + \int_\Gamma \nu (\nabla\mathbf{u})^T\cdot\mathbf{n}\cdot\mathbf{v}\,{\rm d}s$.

is uniformly coercive on the kernel space $\mathbf{\tilde V}_h := \{ \mathbf{v}_h \in \mathbf{V}_h : b(\mathbf{v}_h,q_h) = 0, \forall q\in Q_h \}.$ This is something, which I could not find in the literature, and did not manage to show myself. Perhaps some stabilization terms or sophisticated analysis techniques might be needed.

If one wants to discretize the above problem by the finite element method, one can of course always discretize the alternative form

${\tilde a}(\mathbf{u},\mathbf{v}) := \int_\Omega \tfrac{\nu}{2}(\nabla\mathbf{u} + (\nabla\mathbf{u})^T) : (\nabla\mathbf{v} + (\nabla\mathbf{v})^T)\,{\rm d}x$

where the traction boundaries are no problem to incorporate during the derivation.

However, for constant viscosities it might sometimes be advantageous to not implement the above form with the cross-derivatives involved. Reasons may be that the effort to implement this in some legacy code might be too large or efficiency concerns due to the additional couplings, etc.

tl;dr: Most people who use the Laplacian-Stokes only consider pseudo-traction boundary conditions, i.e., they impose $\nu\nabla \mathbf{u}\cdot\mathbf{n} - p\mathbf{n} = \mathbf{t}$, where it assumed that $(\nabla \mathbf{u})^T\cdot\mathbf{n} \approx 0$. I kind of doubt that there is always a physical justification for this, and rather expect that this is mostly done for convenience. However, I would be surprised if nobody ever investigated the use of the Laplacian form with full traction conditions. I was not very successful yet with my literature study and hope that somebody might have additional resources or some more information for me.

What I have found so far is:

Laplace form of Navier-Stokes equations: A safe path or wrong way? by A. Limache and S. Idelsohn (2006)

This paper raises some similar questions.

Answer 1

You're learning that the two forms (with $-\Delta u$ and with $-2\nabla \cdot \varepsilon(u)$) are not equivalent. They lead to different boundary terms in the bilinear form after integration by parts. In the first case, you get a term involving $n\cdot \nabla u$, in the latter $2n \cdot \varepsilon(u)$, and the "natural" boundary conditions you can enforce are, respectively, the one where either $\nu n\cdot \nabla u-pn=g$ or $2\nu n \cdot \varepsilon(u)-pn=g$. The respectively other one is, simply, not "natural" for each formulation.

From a physical viewpoint, my belief is that the "correct" formulation should be the one with the $-2\nabla \cdot \varepsilon(u)$ term, as it is the one that comes out of modeling flow and also because that's the boundary conditions that make sense. The other formulation may be equivalent when you only have Dirichlet boundary conditions, but it is not equivalent whenever you have natural boundary conditions.

An alternative point of view is to ask what "natural" boundary conditions are anyway. What does it really mean to say that at an outflow segment of the boundary there is no traction on the fluid. To assume that there are no forces on the fluid at a particular point is equivalent to assuming certain things about the flow field past the boundary -- for example, with any boundary condition we may prescribe at such a boundary, we make an assumption that the flow will proceed mostly unobstructed for a long distance, but what exactly is it that we're trying to say with that? When we use no-traction boundary conditions $n \cdot \sigma = 0$ at these places, then this is mostly for convenience because we can not put into formulas what exactly it is that we want to say about what we think is happening outside the domain. But if this is so, whether we say that $n \cdot (\nu\nabla u - pI)=0$ or that $n \cdot (2\nu [\nabla u + \nabla u^T] - pI)=0$, either is about equally arbitrary and simply expresses slightly different notions about the kinds of forces we believe to be zero at this boundary. Given that either is arbitrary, we may as well choose whichever is more convenient in a given situation, given which form we've chosen for the (Navier-)Stokes equations.

Answer 2

In the meantime we did some more research into this issue and documented our findings in:

• M. Huber, U. Rüde, C. Waluga, B. Wohlmuth. Highly sparse surface couplings for subdomain-wise isoviscous Stokes finite element discretizations, Journal of Scientific Computing, 2017 (article, preprint).

In summary, there are a couple of points that can be made:

• Yes, it is possible to combine the Laplacian formulation with some boundary modification to tackle Stokes interface problems and/or traction-type boundary conditions with reduced cross-couplings between the velocity components.
• To make this work we either need finite element spaces (V,Q) for velocity and pressure, respectively, where div(V)=Q, or we need to add some suitable pressure stabilization (we discuss and analyze one choice). The problem with (higher order) spaces that do not require such a stabilization is mostly inf-sup stability and this is a problem of the spaces themselves, not of our modification. For instance in 2D one would have to go to special meshes or use a P4-P3disc-pairing to avoid stability-related issues. In 3D the question of stability of such spaces is much harder and to my best knowledge a long-standing open question (see references in the paper).
• The modification adds boundary integral terms at the interface/boundary that introduce only couplings in tangential direction for H1 conforming velocity spaces.
• The modified Laplacian approach can also be exploited in matrix-free codes in such a way that only stencils at the interfaces need to be replaced and no integration over the boundary facets is required. Combined with suitable matrix-free solvers this can speed-up computations on large-scale architectures with only minor changes to the code.