This is a follow-up to a previous question (see here). We would like to solve the two-dimensional Stokes equations using the FEM package in Mathematica, when we prescribe traction boundary conditions. Let's formulate a concrete example for such a two-dimensional flow problem. Assume that the domain is an annulus with an outer radius $b=1$ and an inner radius $a=1/5$, and generate a mesh:
<< NDSolve`FEM`
mesh = ToElementMesh[Annulus[{0, 0}, {1/5, 1}],
"RegionHoles" -> {0, 0},"MaxCellMeasure" -> 0.005, "MaxBoundaryCellMeasure" -> 0.05];
and display the mesh and the corresponding boundary element markers:
GraphicsGrid[{{mesh["Wireframe"], mesh["Wireframe"["MeshElement" -> "BoundaryElements","MeshElementMarkerStyle" -> Black,"MeshElementStyle" -> {Red,Blue}]]}}]
This allows us to see that (by default) the inner boundary is identified by ElementMarker == 1 and the outer boundary is given by ElementMarker == 2 which is useful to know when setting up the boundary conditions.
Let's setup the boundary conditions. I would like to prescribe the physical traction at the inner boundary; namely, the stress vector at that boundary is given by $$\left.\boldsymbol{n}\cdot\boldsymbol{\sigma}\right|_{r=a} = \boldsymbol{\hat{r}}\,N(\theta) + \boldsymbol{\hat{\theta}}\,S(\theta),$$
where for simplicity we take the normal traction to be $N(\theta)=0$, whereas the the tangential component is given by $S(\theta)=\cos(2\theta)$ (similar to a squirmer model). Here, $\boldsymbol{\hat{r}}=\boldsymbol{\hat{\imath}}\cos\theta +\boldsymbol{\hat{\jmath}}\sin\theta $ is the radial unit vector, and the tangent unit vector $\boldsymbol{\hat{\theta}}=-\boldsymbol{\hat{\imath}}\sin\theta +\boldsymbol{\hat{\jmath}}\cos\theta$, with $\boldsymbol{\hat{\imath}}$ and $\boldsymbol{\hat{\jmath}}$ as the Cartesian unit vectors (along $x$ and $y$ respectively). On the other boundary, we set both the traction and the velocity to be zero: $$\left.\boldsymbol{n}\cdot\boldsymbol{\sigma}\right|_{r=b} = \boldsymbol{0}; \qquad\left.\boldsymbol{u}\right|_{r=b}=0.$$ Note that the velocity is not specified at the inner disk, only the traction, $\boldsymbol{n}\cdot\boldsymbol{\sigma}$, is prescribed. Here, $\boldsymbol{\sigma}$ is the fluid stress tensor, that is, $$\boldsymbol{\sigma} = -p\,\mathsf{I}+\mu\left[(\boldsymbol{\nabla}\boldsymbol{u})+(\boldsymbol{\nabla}\boldsymbol{u})^{\mathsf{T}}\right],$$ which can be written (just to be clear) in the component form as follows: $\sigma_{ij} = -p\hspace{1pt}\delta_{ij} + \mu\left(\partial_i u_j+\partial_j u_i\right)$, with $u_i$ as the Cartesian components of the velocity field $\boldsymbol{u}$, and $p$ being the fluid pressure. Stokes equations are given by the force balance equation (note that is a vector equation), $$\boldsymbol{\nabla}\cdot\boldsymbol{\sigma} = \boldsymbol{0},$$ and the incompressibility condition, $$\boldsymbol{\nabla}\cdot\boldsymbol{u} = 0.$$ The latter equation is typically used to simplify the former equation to $-\mu\hspace{1pt}\Delta\boldsymbol{u}+\boldsymbol{\nabla}p = \boldsymbol{0}$. Although the form of this equation is identical to $\boldsymbol{\nabla}\cdot\boldsymbol{\sigma} = \boldsymbol{0}$, their respective Neumann conditions are not (this is briefly discussed here). Thus, in order to get the physical traction as a Neumann boundary condition, the initial form of the force-balance equation must be retained. Let's rewrite this in Cartesian form (using Mathematica):
pdes = Div[-p[x, y] IdentityMatrix[2]
+\[Mu] Grad[{u[x, y], v[x, y]}, {x, y}]
+\[Mu] Transpose[Grad[{u[x, y], v[x, y]}, {x, y}]], {x,y}]//FullSimplify
where $\boldsymbol{u} = u \boldsymbol{\hat{\imath}} + v\boldsymbol{\hat{\jmath}} = (u,v)$. Hence, along the $x$-direction, we have $$-\frac{\partial}{\partial x}\hspace{1pt}p(x,y) + \mu\left(\frac{\partial^2}{\partial y^2}\hspace{1pt}u(x,y)+\frac{\partial^2}{\partial x\hspace{1pt}\partial y}\hspace{1pt}v(x,y)+2\frac{\partial^2}{\partial x^2}\hspace{1pt}u(x,y)\right) = 0,$$ whereas, along the $y$-direction, we have
$$-\frac{\partial}{\partial y}\hspace{1pt}p(x,y) + \mu\left(\frac{\partial^2}{\partial x^2}\hspace{1pt}v(x,y)+\frac{\partial^2}{\partial x\hspace{1pt}\partial y}\hspace{1pt}u(x,y)+2\frac{\partial^2}{\partial y^2}\hspace{1pt}v(x,y)\right) = 0.$$ Note that the mixed derivatives can be removed by using the incompressibility condition, that is,
$$\frac{\partial}{\partial x}\hspace{1pt}u(x,y)+\frac{\partial}{\partial y}\hspace{1pt}v(x,y)=0,$$ and to retrieve the standard form of Stokes equations. In Mathematica, we would write these equations in their Inactive form as follows:
ipde1 = Inactive[Div][Inactive[Plus][
{{0, 0}, {-\[Mu], 0}}.Inactive[Grad][v[x, y], {x, y}],
{{-2 \[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][u[x, y], {x, y}],
Inactive[Times][{1, 0}, p[x, y]]
], {x, y}]
ipde2 = Inactive[Div][Inactive[Plus][
{{0, -\[Mu]}, {0, 0}}.Inactive[Grad][u[x, y], {x, y}],
{{-\[Mu], 0}, {0, -2 \[Mu]}}.Inactive[Grad][v[x, y], {x, y}],
Inactive[Times][{0, 1}, p[x, y]]
], {x, y}]
Notice that the activated form of these equations Activate[{ipde1,ipde2}]=={0,0} is equivalent to those in pdes=={0,0}; the minus sign change is because we want to express the equations in Mathematica's canonical coefficient form of PDEs: 
Now, using the code written by @user21 (see the answer here) we can define the domain:
nr = ToNumericalRegion[mesh];
vd = NDSolve`VariableData[{"DependentVariables","Space"} -> {{u, v, p}, {x, y}}];
sd = NDSolve`SolutionData[{"Space"} -> {nr}];
mdata = InitializePDEMethodData[vd, sd, Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];
In addition, to setting up the pressure in the conservative convection form, we also need to include the other diffusion coefficients that arise from the mixed terms, as seen in the inactivated form of the Stokes equations:
cdata = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {
{{{-2 \[Mu], 0}, {0, -\[Mu]}}, {{0, 0}, {-\[Mu], 0}}, 0},
{{{0, -\[Mu]}, {0, 0}}, {{- \[Mu], 0}, {0, -2 \[Mu]}}, 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}}
}
];
Lastly, set up the boundary conditions (where t denotes the angle $\theta$):
bcdata = InitializeBoundaryConditions[vd, sd, {
{DirichletCondition[u[x, y] == 0, ElementMarker == 2],
NeumannValue[0, ElementMarker == 2]
NeumannValue[-(-Sin[t])*(Cos[2t])/.{t -> ArcTan[x,y]}, ElementMarker == 1]},
{DirichletCondition[v[x, y] == 0, ElementMarker == 2],
NeumannValue[0, ElementMarker == 2]
NeumannValue[-(+Cos[t])*(Cos[2t])/.{t -> ArcTan[x,y]}, ElementMarker == 1]},
{}
}];
As mentioned in the beginning, ElementMarker == 2 identifies the outer boundary elements, while ElementMarker == 1 gives us the inner disk. Note that (-Sin[t]) and (Cos[t]) terms are just the $x$- and $y$-components of the unit vector $\boldsymbol{\hat{\theta}}$. The extra minus sign accounts for the overall sign change of the equations, with the Neumann value signifying the negative of the prescribed physical traction (which in case is only tangential and varies as $\cos(2\theta)$, that is, Cos[2t] term).
Now, process and numerically solve:
sdNew = PDESolve[cdata, bcdata, vd, sd, mdata];
{xVel, yVel, pressure} = ProcessPDESolutions[mdata, sdNew]
Plot the numerical results for the velocity as a StreamPlot and the fluid pressure as DensityPlot:
Show[{
DensityPlot[pressure[x,y],{x,y} ∈ mesh,PlotRange -> All,
ColorFunction -> "LightTemperatureMap", BoundaryStyle -> Directive[Gray],
PlotLegends -> Placed[BarLegend[Automatic, LegendMarkerSize -> 400,
LegendLabel -> Style["p(x,y)", 18, FontFamily->"Times"],
LabelStyle -> Directive[Black, FontSize -> 12, FontFamily -> "Times"]],Right]],
StreamPlot[{xVel[x,y],yVel[x,y]},{x,y} ∈ mesh, PlotRange->All,
StreamColorFunction -> Automatic, StreamPoints -> Fine,
RegionBoundaryStyle -> None, RegionFillingStyle -> None],
mesh["Wireframe"["MeshElementStyle"->EdgeForm[Directive[Black,Opacity[0.03]]]]]
}, Frame -> True, FrameStyle -> Directive[Black, FontSize -> 12, FontFamily -> "Times"],
FrameLabel -> {Style["x",FontSize->18],Style["y",18]}, ImageSize -> 480]
Looks promising! Now, let's check whether this numerical solution agrees with the analytical solution. The radial symmetry of our domain allows us to solve this problem in exact form. I won't spend too much time on how one obtains such a solution, but I would be more than happy if someone else could also check my results (see below). The idea is to rewrite the equations and to separately solve for the stream function and the pressure. Namely, by taking the divergence of the Stokes equations we end up with $$\Delta p =0.$$ The pressure is simply a harmonic function within the domain. By taking instead the (2D) curl of the Stokes equation, and also defining the stream function $\Psi$ to be given by $u=\partial_y\Psi$ and $v=-\partial_x\Psi$, we find that $$\Delta^2 \Psi=0$$ which means that $\Psi$ solves the biharmonic equation. We work in polar coordinates $(r,\theta)$, where we write that $\boldsymbol{u}=v_r\boldsymbol{\hat{r}}+v_\theta\boldsymbol{\hat{\theta}}$, with $v_r = \frac{1}{r}\frac{\partial}{\partial\theta}\Psi(r,\theta)$ and $v_\theta = -\frac{\partial}{\partial r}\Psi(r,\theta)$. Moreover, the stress vector can now be rewritten as follows: $$\boldsymbol{n}\cdot\boldsymbol{\sigma} = \left[-p(r,\theta)+2\mu\left(\frac{1}{r}\frac{\partial^2\Psi}{\partial\theta\partial r}-\frac{1}{r^2}\frac{\partial\Psi}{\partial\theta}\right)\right]\boldsymbol{\hat{r}} +\mu\left(\frac{1}{r^2}\frac{\partial^2\Psi}{\partial \theta^2}+\frac{1}{r}\frac{\partial\Psi}{\partial r}-\frac{\partial^2\Psi}{\partial r^2}\right)\boldsymbol{\hat{\theta}}.$$ We can solve the Laplace equation and the biharmonic equation by means of separation of variables. Knowing that the tangential traction at $r=a$ is $\cos(2\theta)$, we find that the radial velocity to be: $$v_r(r,\theta) = \frac{a^4(b^2-r^2)^3\sin(2\theta)}{6(b^2-a^2)(a^4+b^4)r^3\mu},$$ the azimuthal velocity is given by $$v_\theta(r,\theta) = -\frac{a^4(b^6-3b^2r^4+2r^6)\cos(2\theta)}{6(b^2-a^2)(a^4+b^4)r^3\mu},$$ and pressure $$p(r,\theta)=-\frac{a^2(b^4-r^4)\sin(2\theta)}{(a^4+b^4)r^2}.$$
Let's plot these fields:
So, it seems that there is a major disagreement between the analytical solution and the numerical results from before! I'm not entirely sure what went wrong here. Perhaps I'm still not implementing correctly the Stokes equations in Mathematica. Can you please help me with this?
{{0, 0}, {-\[Mu], 0}}should that not be{{0, -\[Mu]}, {0, 0}}? It seems that the issue is with the inner boundary. $\endgroup$BoundaryUnitNormalto model the normal/tangent. WithNeumannValue[-Indexed[BoundaryUnitNormal[x, y]*(Cos[2 ArcTan[x, y]]), 2], ElementMarker == 1]andNeumannValue[ Indexed[BoundaryUnitNormal[x, y]*(Cos[2 ArcTan[x, y]]), 1], ElementMarker == 1]for the x and y component, respectively. $\endgroup$