# Stokes equations in 2D with Neumann conditions

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["NDSolveFEM"]
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.

• Do you have a literature reference for this, perhaps even an expected solution? Mar 24, 2021 at 11:25
• No, it is impossible to implement mixed boundary conditions in the current version 12.2. In this case we have message NDSolveValue::fembdcc: Cross-coupling of dependent variables in DirichletCondition[(u x)/Sqrt[x^2+y^2]-(v y)/Sqrt[x^2+y^2]==G[x,y],True] is not supported in this version. Mar 24, 2021 at 20:31
• @AlexTrounev, the error message you show is for DirichletConditions. I do not think that this is relevant here. Mar 25, 2021 at 8:20
• @user21 Yes you are absolutely right, DirichletCondotion[] is not relevant since Cross-coupling of dependent variables in DirichletCondition[] is not supported in this version. But it can be relevant if somebody from developer will think about it. Mar 25, 2021 at 10:18
• @Alex R Did you take this example from Fenics Project on fenicsproject.org/olddocs/dolfin/1.5.0/python/demo/documented/… ? Mar 25, 2021 at 11:41

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 = NDSolveVariableData[{"DependentVariables",
"Space"} -> {{u, v, p}, {x, y}}];
sd = NDSolveSolutionData[{"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 NDSolveFEMBoundaryUnitNormal (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.

• Very nice answer! But how we can implement bc with G=1; H={1,1}? Mar 25, 2021 at 10:34
• @AlexTrounev, I think it is better to discuss this with a more concrete example. Mar 25, 2021 at 10:51
• Well, let take $\vec{u}=0, -p \vec {n}+ \frac {\partial \vec {u}}{\partial \vec {n}}=0$ like in this paper mat.fsv.cvut.cz/nales/preprints/preprinty/2009/ulozeneclanky/… Mar 25, 2021 at 11:09
• @AlexTrounev, That is DirichletCondition[{u[x,y]==0,v[x,y]==0}, walls] and Neumann 0 for the remainder. The BC 1 from the paper is probably not correct. It misses a factor nu. Mar 25, 2021 at 13:34
• And in fact equation 5 has the factor. Mar 25, 2021 at 13:43