# Stokes equations in 2D with traction boundary conditions

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:

<< NDSolveFEM
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
+\[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 = NDSolveVariableData[{"DependentVariables","Space"} -> {{u, v, p}, {x, y}}];
sd = NDSolveSolutionData[{"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?

• The color maps look fine except for the streamplots. Perhaps you used Cartesian coordinates in the 1st one, and Polar coordinates for your 2nd plot? Apr 12, 2021 at 16:43
• @Ferca I've doubled checked that. I used Cartesian coordinates in both cases. Pressure fields are different too. It's probably quite hard to see that in those plots. My money is on the numerical results being wrong; possibly in the way I've implemented the equations Apr 12, 2021 at 19:12
• You say: "Note that the mixed derivatives can be removed by..." but it seems to me that you do not actually remove them. Is this correct? Also you use {{0, 0}, {-\[Mu], 0}} should that not be {{0, -\[Mu]}, {0, 0}}? It seems that the issue is with the inner boundary. Apr 13, 2021 at 6:28
• You solution, does not that have 0 normal velocity at the inner boundary. Should it not? Apr 13, 2021 at 7:46
• Another way to model this is to use BoundaryUnitNormal to model the normal/tangent. With NeumannValue[-Indexed[BoundaryUnitNormal[x, y]*(Cos[2 ArcTan[x, y]]), 2], ElementMarker == 1] and NeumannValue[ Indexed[BoundaryUnitNormal[x, y]*(Cos[2 ArcTan[x, y]]), 1], ElementMarker == 1] for the x and y component, respectively. Apr 13, 2021 at 7:47

As mentioned user21 for setting traction boundary conditions we need to know normal and tangent to the surface which is subjected to loading. On the inner surface ($$r=a$$) for the unit tangent vector we have $$\vec{\tau}=(-y/a,x/a)$$ that corresponds to anticlockwise direction. On the inner surface the normal component of traction vector is zero, whereas shear stress $$S(\theta)=cos{2\theta}=cos^2{\theta}-sin^2{\theta}=(x/a)^2-(y/a)^2$$. On the outer surface ($$r=b$$) no slip conditions $$u=v=0$$ are applied. Such mixed BC can be implemented by the next way.

Needs["NDSolveFEM"];
(*geometry of computational domain*)
a = 1/5;
b = 1;
(*dynamic viscosity of the fluid*)
mu = 1;
(*FE mesh generation*)

(*uniform triangular mesh*)

mesh1 = ToElementMesh[Annulus[{0, 0}, {a, b}],
"RegionHoles" -> {0, 0}, "MaxCellMeasure" -> 0.001,
"MaxBoundaryCellMeasure" -> 0.05];

(*adaptive mesh refinement near the inner surface*)

mesh2 = ToElementMesh[Annulus[{0, 0}, {a, b}],
"RegionHoles" -> {0, 0},
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.000025 (0.1 + 20 Norm[Mean[vertices]])]];

(*selection of FE mesh for calculation *)
mesh = mesh2;
(* mesh=mesh1 *)


We will use FE programming routines

Clear[vd, sd, methodData, discretePDE, initCoeffs, InitBCs,discreteBCs]
vd = NDSolveVariableData[{"DependentVariables",
"Space"} -> {{u, v, p}, {x, y}}];
sd = NDSolveSolutionData[{"Space"} -> {mesh}];
methodData =
InitializePDEMethodData[vd, sd,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];


Define the coefficients of PDE

initCoeffs = InitializePDECoefficients[vd, sd,
"DiffusionCoefficients" -> {
{{{-2 mu,0}, {0, -mu}}, {{0, 0}, {-mu, 0}}, {{0, 0}, {0, 0}}},
{{{0,-mu}, {0, 0}}, {{-mu, 0}, {0, -2 mu}}, {{0, 0}, {0, 0}}},
{0, 0, 0}
},

"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}}
}

];


Calculation of stiffness matrix and load vector

discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
split = Span @@@Transpose[{Most[# + 1], Rest[#]} &[
methodData["IncidentOffsets"]]];


Implementing of the boundary conditions

GammaDu = DirichletCondition[u[x, y] == 0., ElementMarker == 2];
GammaNu = NeumannValue[(-y/a)*((x/a)^2 - (y/a)^2), ElementMarker == 1];

GammaDv = DirichletCondition[v[x, y] == 0., ElementMarker == 2];
GammaNv = NeumannValue[(x/a)*((x/a)^2 - (y/a)^2), ElementMarker == 1];

};

InitBCs =
InitializeBoundaryConditions[vd,
discreteBCs = DiscretizeBoundaryConditions[InitBCs, methodData, sd];



Solution of the linear system of equations

solution = LinearSolve[stiffness, load, "Method" -> "Pardiso"]


Construction of interpolation functions for velocity and pressure

ufun = ElementMeshInterpolation[{mesh}, solution[[split[]]]]
vfun = ElementMeshInterpolation[{mesh}, solution[[split[]]]]
pfun = ElementMeshInterpolation[{MeshOrderAlteration[mesh, 1]},
solution[[split[]]]]


Visualization of pressure and velocity fields

Show[{DensityPlot[pfun[x, y], {x, y} \[Element] 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[{ufun[x, y], vfun[x, y]}, {x, y} \[Element] 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] Velocity components in polar CS can be obtained by

(*radial component*)
Vr[r_, th_] := Cos[th]*ufun[Cos[th]*r, Sin[th]*r] +
Sin[th]*vfun[Cos[th]*r, Sin[th]*r]
(*azimuthal component*)
Vt[r_, th_] := -Sin[th]*ufun[Cos[th]*r, Sin[th]*r] +
Cos[th]*vfun[Cos[th]*r, Sin[th]*r]


Distribution of $$V_{r}(r,\theta)$$ along line $$\theta=const$$ will be the next Plot[{Vr[a, th], 0.02 Sin[2 th]}, {th, 0, 2 Pi},
PlotRange -> All, Frame -> True, FrameLabel -> {"\[Theta]", "Vr"},
PlotStyle -> {{RGBColor[0, 0, 1],
Thickness[0.01]}, {RGBColor[1, 0, 0]}}, BaseStyle -> 18,
FrameStyle -> RGBColor[0, 0, 0], ImageSize -> 700,
PlotLabel -> "radial velocity on inner surface",
GridLines -> {{Pi/4, Pi/4 + Pi/2, Pi + Pi/4, 3 Pi/2 + Pi/4}, None},
GridLinesStyle -> {{Thickness[0.003], Dashed}, {}},
PlotLegends -> {"Vr(a,\[Theta])", "0.02sin(2\[Theta])"}]

Plot[{Vt[a, th], 0.054 Cos[2 th]}, {th, 0, 2 Pi},
PlotRange -> All, Frame -> True, FrameLabel -> {"\[Theta]", "Vt"},
PlotStyle -> {{RGBColor[0, 0, 1],
Thickness[0.01]}, {RGBColor[1, 0, 0]}}, BaseStyle -> 18,
FrameStyle -> RGBColor[0, 0, 0], ImageSize -> 700,
PlotLabel -> "tangential velocity on inner surface",
GridLines -> {{Pi/4, Pi/4 + Pi/2, Pi + Pi/4, 3 Pi/2 + Pi/4}, None},
GridLinesStyle -> {{Thickness[0.003], Dashed}, {}},
PlotLegends -> {"Vt(a,\[Theta])", "0.054cos(2\[Theta])"}]


For $$V_{r}(a,\theta)$$ and $$V_{\theta}(a,\theta)$$ on inner surface we have  As expected the maximum of azimuthal velocity on inner surface corresponds to the maximum values of shear stress

Calculation of the shear and normal stresses on inner surface

Under the given velocity and pressure fields let's calculate shear and normal stresses on inner surface by

$$S_{num}=-\sigma_{r\theta}\Big |_{r=a}=-\mu\left(\frac{1}{r}\frac{\partial V_r}{\partial \theta}+\frac{\partial V_{\theta}}{\partial r}-\frac{V_{\theta}}{r}\right)=-\mu\left( \frac{\partial x}{\partial r}\frac{\partial V_{\theta}}{\partial x}+\frac{\partial y}{\partial r}\frac{\partial V_{\theta}}{\partial y} +\frac{1}{r}(\frac{\partial x}{\partial \theta}\frac{\partial V_r}{\partial x}+\frac{\partial y}{\partial \theta}\frac{\partial V_r}{\partial y})-\frac{V_{\theta}}{r} \right)$$

$$N_{num}=\sigma_{rr}\Big |_{r=a}=-p+2\mu\frac{\partial V_r}{\partial r}=-p+2\mu\left( \frac{\partial x}{\partial r}\frac{\partial V_r}{\partial x}+\frac{\partial y}{\partial r}\frac{\partial V_r}{\partial y} \right)$$

Clear[Vr, Vt];
Vr[x_, y_] :=
x/Sqrt[x^2 + y^2]*ufun[x, y] + y/Sqrt[x^2 + y^2]*vfun[x, y]
Vt[x_, y_] := -(y/Sqrt[x^2 + y^2])*ufun[x, y] +
x/Sqrt[x^2 + y^2]*vfun[x, y]

VrInt = ElementMeshInterpolation[{mesh}, VrData];
VtInt = ElementMeshInterpolation[{mesh}, VtData];

ShStress[x_, y_] := Module[{r},
r = Sqrt[x^2 +
y^2]; -mu*(x/
r*(Derivative[1, 0][VtInt][x, y] +
Derivative[0, 1][VrInt][x, y]) +
y/r*(Derivative[0, 1][VtInt][x, y] -
Derivative[1, 0][VrInt][x, y]) - VtInt[x, y]/r)
]

NormStress[x_, y_] := Module[{r},

r = Sqrt[x^2 + y^2];

pfun[x, y] -
2 mu*(x*Derivative[1, 0][VrInt][x, y] +
y*Derivative[0, 1][VrInt][x, y])/r
]


It is interesting to compare calculated surface stresses with those we applied as a boundary conditions in NeumanValue.

datSnum =
Table[{th, ShStress[a*Cos[th], a*Sin[th]]}, {th, 0, 2 Pi,
2 Pi/64}];
datNnum =
Table[{th, NormStress[a*Cos[th], a*Sin[th]]}, {th, 0, 2 Pi,
2 Pi/64}];
datS = Table[{th, Cos[2 th]}, {th, 0, 2 Pi, 2 Pi/64}];

ListLinePlot[{datS, datSnum, datNnum},
PlotLegends -> {"cos(2\[Theta])",
"shear stress from numerical solution",
"normal stress from numerical solution"}, Frame -> True,
BaseStyle -> 18, PlotStyle -> AbsoluteThickness,
FrameLabel -> {"\[Theta], rad", "stress"}, ImageSize -> 400,
FrameStyle -> RGBColor[0, 0, 0]] Mesh refinement near the inner surface leads to diminishing of the discrepancy between calculated and applied in BC surface forces. It can be concluded that traction BC are implemented correctly.

• Great. I think the diffusion coefficient should be: { {{{-2 mu, 0}, {0, -mu}}, {{0, -mu}, {0, 0}}, 0}, {{{0, 0}, {-mu, 0}}, {{-mu, 0}, {0, -2 mu}}, 0}, {0, 0, 0} Apr 15, 2021 at 11:18
• Your numerical solution seems similar to my numerical solution (at least, up to a minus sign!). Indeed, at the inner disk, I expect that the tangential velocity to be $V_\theta \propto \cos(2\theta)$ and the radial velocity to be $V_r\propto \sin(2\theta)$. Though I still find it rather strange that the radial velocity is not monotonically decreasing with $r$ in both our numerical results (there shouldn't be a bump in the radial profile!). So, I don't think that $\cos(2\theta)$ in the argument of NeumannValue corresponds to the actual tangential traction. Apr 15, 2021 at 21:54
• @user21 Thanks. The array for diffusion coefficients you wrote and those from my answer give the same results if $\mu=const$ Apr 16, 2021 at 11:20
• @AlexR Of course numerical results should give $V_{\theta}\propto \cos{2\theta}$ and $V_{r}\propto \sin{2\theta}$ as you expected. Let you see updated unswer. Apr 16, 2021 at 11:43
• @AlexR It is nescessary to calculate $-\mu\frac{\partial V_{\theta}}{\partial r}$ on inner surface to be shure that NeumanValue adiquately approximates shear stress. Apr 16, 2021 at 12:01

I just wanted to point out two things: 1) You can use top level code for the Stokes operator, you just need to split it in a sum of Divs, and 2) You can use BoundaryUnitNormal for the normal and tangent computation. While computing the normal and tangent is possible to do symbolically here, it might not be in other cases. The BoundaryUnitNormal comes in handy.

Here is the code:

Needs["NDSolveFEM"]
pars = Association[a -> 1/5, b -> 1, mu -> 1];
region = Annulus[{0, 0}, {a, b}] /. pars;
mesh = ToElementMesh[region, "RegionHoles" -> {0, 0},
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.000025 (0.1 + 20 Norm[Mean[vertices]])]];


The Stokes flow as a sum of Inactive Divs

op = {Inactive[Div][{{-2*mu, 0}, {0, -mu}} .
Inactive[Grad][u[x, y], {x, y}], {x, y}] +
Inactive[Div][{{0, 0}, {-mu, 0}} .
Inactive[Grad][v[x, y], {x, y}], {x, y}] +
Inactive[Div][{1, 0}*p[x, y], {x, y}],
Inactive[Div][{{0, -mu}, {0, 0}} .
Inactive[Grad][u[x, y], {x, y}], {x, y}] +
Inactive[Div][{{-mu, 0}, {0, -2*mu}} .
Inactive[Grad][v[x, y], {x, y}], {x, y}] +
Inactive[Div][{0, 1}*p[x, y], {x, y}],
Derivative[1, 0][u][x, y] + Derivative[0, 1][v][x, y]
} /. pars;


Boundary conditions:

GammaWall =
DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, ElementMarker == 2];


Here is the place where we use the BoundaryUnitNormal and a 'Cross' product to compute the tangent and extract the u and v parts, respectively:

GammaU =
NeumannValue[
Indexed[Cross[-BoundaryUnitNormal[x, y]], 1]*((x/a)^2 - (y/a)^2),
ElementMarker == 1];
GammaV =
NeumannValue[
Indexed[Cross[-BoundaryUnitNormal[x, y]], 2]*((x/a)^2 - (y/a)^2),
ElementMarker == 1];


Solving this is now much simpler, because we can directly call NDSolve:

{ufun, vfun, pfun} =
NDSolveValue[{op == {GammaU, GammaV, 0} /. pars, GammaWall}, {u, v,
p}, {x, y} \[Element] mesh,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}},
DependentVariables -> {u, v, p}];


With out specifying a boundary condition for the pressure, the pressure value will be floating and NDSolve will give a warning that not enough boundary conditions are specified.

Visualization:

Show[{DensityPlot[pfun[x, y], {x, y} \[Element] mesh, Sequence[
ColorFunction -> "LightTemperatureMap", PlotLegends -> Placed[
BarLegend[Automatic, LegendLabel -> "p(x,y)"], Right],
PlotRange -> All]],
StreamPlot[{ufun[x, y], vfun[x, y]}, {x, y} \[Element] mesh,
Sequence[
StreamPoints -> Fine, RegionFillingStyle -> None]]}, Sequence[
FrameLabel -> {"x", "y"}, ImageSize -> Medium]] 