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.
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.
$\endgroup$DirichletCondotion[]
is not relevant sinceCross-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. $\endgroup$