Stokes equations in 2D with Neumann conditions

Question

I would like to solve the 2D Stokes equations within a unit disk, say $\Omega$, by using the finite element method (FEM) as it is implemented in NDSolve (by loading the finite element package). The examples shown in the FEM documentation only deal with Dirichlet boundary conditions (see Fluid flow section). Specifically, I want to numerically find $\boldsymbol{u}=(u,v)$ that solves $-\mu \nabla^2\boldsymbol{u}+\nabla p =0$ and $\nabla\cdot\boldsymbol{u}=0$, everywhere within the disk $\Omega$. Here $u$ is the fluid velocity in the $x$ direction and $v$ represents the velocity in the $y$ direction, whist $p$ stands for the pressure in the fluid. On the boundary, $\partial\Omega$, both the normal velocity and the tangential component of the traction are given by$$\left.\boldsymbol{n}\cdot\boldsymbol{u}\right|_{\partial\Omega}=G$$ $$(\mathsf{I}-\boldsymbol{n}\otimes\boldsymbol{n})\left[\boldsymbol{n}\cdot(\mu \nabla\boldsymbol{u}-p\hspace{1pt}\mathsf{I})\right]_{\partial\Omega} = H$$ where $\boldsymbol{n}$ is the unit normal vector and $\mathsf{I}$ is the identity matrix; $G$ and $H$ are prescribed functions.

Here's a starting setup:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[Disk[]];
op = {
Inactive[Div][{{-μ, 0}, {0, -μ}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] + Derivative[1, 0][p][x, y], 
Inactive[Div][{{-μ, 0}, {0, -μ}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] + Derivative[0, 1][p][x, y], 
Derivative[1, 0][u][x, y] + Derivative[0, 1][v][x, y]
} /. {μ -> 1};
NDSolveValue[op=={0,0,0}, {u, v, p}, {x, y} ∈ mesh]]

How can one include the above boundary conditions using the NeumannValue function on the $u$ and $v$ equations; for instance, the strategy

NDSolveValue[op=={NeumannValue[1,True],NeumannValue[1,True],0}, {u, v, p}, {x, y} ∈ mesh]]

does not quite work as this only provides the fluxes on $u$ and $v$ (namely, the velocity gradients at the boundary) ignoring the pressure term. Can this even be implemented with the FEM NSolve package?

Edit: Stokes equations are linear, and in the case of a unit disk, as I described the problem above, one could actually solve for this in an analytical form, by using a stream function (that satisfies $u=\partial\Psi/\partial y$ and $v=-\partial\Psi/\partial x$). However, I'm not particularly interested in such radially symmetric solutions (although this would be useful as a check), since I would like to extend this to a nontrivial boundary, and thus FEM is a must. So, I am wondering whether the FEM package in Mathematica can deal with such traction boundary conditions. Here's a short lecture that covers this topic. NDSolve seems to be doing the weak formulation internally and leaves very little options to the user.

Answer 1

Here is how I'd approach that. The first thing to realize is that you equation is not in the correct form. You have a non-conservative form of the stokes equation. Your boundary condition, however, suggests that you need a conservative form (You can read up on that in the MassTransport tutorial).

Unfortunately, the parser does not handle that in a straight forward manner so we resort to the low level functions (See finite element programming).

Set up the numerical region, variables, and solution data:

nr = ToNumericalRegion[\[CapitalOmega]];
vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u, v, p}, {x, y}}];
sd = NDSolve`SolutionData[{"Space"} -> {nr}];

Set up the method data:

mdata = InitializePDEMethodData[vd, sd, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

Now, comes the key point. The pressure needs to be in conservative form:

cdata = InitializePDECoefficients[vd, sd,
   "DiffusionCoefficients" -> {
      {-mu IdentityMatrix[2], 0, 0},
      {0, -mu IdentityMatrix[2], 0},
      {0, 0, 0}
      } /. mu -> 1,
   "ConservativeConvectionCoefficients" -> {
     {{0, 0}, {0, 0}, {1, 0}},
     {{0, 0}, {0, 0}, {0, 1}},
     {{0, 0}, {0, 0}, {0, 0}}
     },
   "ConvectionCoefficients" -> {
     {{0, 0}, {0, 0}, {0, 0}},
     {{0, 0}, {0, 0}, {0, 0}},
     {{1, 0}, {0, 1}, {0, 0}}
     }
   ];

Note, how the pressure is in the conservative convection form, while the divergence of the velocity is in the convection coefficient.

Now, you can use NeumannValue to reason about $(\mu \nabla\boldsymbol{u}-p)$ on the boundary. The other thing you are likely going to need is NDSolve`FEM`BoundaryUnitNormal (Also shown in the mass transport tutorial).

Unfortunately, you did not really provide an example to work on so I can not really show how to set up the BCs. The remainder would then be set up like this:

bcdata = InitializeBoundaryConditions[vd, sd, {
    {},
    {},
    {}
    }];
sdNew = PDESolve[cdata, bcdata, vd, sd, mdata];
{xVel, yVel, pressure} = ProcessPDESolutions[mdata, sdNew]

Would be nice to see a result or an updated question.