Simplify fluid flow by using symmetry property

Question

I am trying to describe the fluid flow within a radially symmetric tube and use the code as it can be found from Wolframs homepage as template (Stokes's equation).

I adapted the geometry and the inlet flow to my given system and get a reasonable solution which is nice (parabolic laminar velocity profile):

Now, with the goal to reduce computational demand and to use later on the common balance equations for radially symmetric systems I want to solve only the upper half of the system depicted above. By doing so the centre line of the tube, which is at r = 0.0024, is going to be the new x axis.

In principle this can be done by using the symmetry of the system and setting up Neumann boundary conditions at the bottom (r == 0):

$\frac{\partial u}{\partial r}\left ( r=0 \right )=\frac{\partial v}{\partial r}\left ( r=0 \right )=\frac{\partial p}{\partial r}\left ( r=0 \right )=0$

Unfortunately here I am at the point where I get stuck. By using the FE method NeumannValues need to be defined. According to the documentation NeumannValue specifies the value on the right side of this equation with u being the dependent variable:

$\vec{n}\left ( c\triangledown u+\alpha u-\gamma \right )=g-qu$

But how does it look like when I have PDEs with more than one dependent variable like in $\frac{\partial u}{\partial x}+\frac{\partial v}{\partial r}=0$ ?

Although I am not sure about the NeumannValue definition for a first test I set it to zero using this code:

Needs["NDSolve`FEM`"]

Ω = Rectangle[{0, 0}, {2.0016, 0.0024}]; 
RegionPlot[Ω, AspectRatio -> 0.5]

op = {D[μ*D[u[x, r], x], x] + D[μ*D[u[x, r], r], r] - 
 D[p[x, r], x], 
D[μ*D[v[x, r], x], x] + D[μ*D[v[x, r], r], r] - 
 D[p[x, r], r], D[u[x, r], x] + D[v[x, r], r]} /. μ -> 0.0001; 

Subscript[Γ, N] = {NeumannValue[0, r == 0 && x > 0], 
NeumannValue[0, r == 0 && x > 0], 
NeumannValue[0, r == 0 && x > 0]}; 

Subscript[Γ, 
D] = {DirichletCondition[u[x, r] == 2*0.065*(1 - (r/0.0024)^2), 
x == 0.], 
DirichletCondition[{u[x, 0.0024] == 0., v[x, 0.0024] == 0.}, 
Inequality[0, Less, x, LessEqual, 2.0016]], 
   DirichletCondition[p[x, r] == 0., x == 2.0016]}; 

mesh = ToElementMesh[Ω, 
MaxCellMeasure -> {"Length" -> 0.0012}]; 

{xVel, yVel, pressure} = 
NDSolveValue[{op == Subscript[Γ, N], 
Subscript[Γ, D]}, {u, v, p}, Element[{x, r}, mesh], 
Method -> {"FiniteElement", 
 "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]; 

ContourPlot[xVel[x, r], {x, 0, 2}, {r, 0, 0.0024}, 
AspectRatio -> 0.5, ColorFunction -> "TemperatureMap", 
PlotLegends -> Automatic]

But apparently this is wrong as there is a flux at the bottom border causing the axial velocity to be zero where it should be at its maximum:

Any help is appreciated, thanks.

Answer 1

No need for Neumann BCs at all. The independent variables got messed up a bit. This works:

Needs["NDSolve`FEM`"]

xL = 2.0016; yL = 0.0024;
Ω = Rectangle[{0, 0}, {xL, yL}];
RegionPlot[Ω, AspectRatio -> 0.5]

op = {D[μ*D[u[x, r], x], x] + D[μ*D[u[x, r], r], r] - 
     D[p[x, r], x], 
    D[μ*D[v[x, r], x], x] + D[μ*D[v[x, r], r], r] - 
     D[p[x, r], r], D[u[x, r], x] + D[v[x, r], r]} /. μ -> 
    0.0001;
Subscript[Γ, D] = {
   DirichletCondition[u[x, r] == 2*0.065*(1 - (r/yL)^2), x == 0.], 
   DirichletCondition[{u[x, r] == 0., v[x, r] == 0.}, 
    0 < x <= xL && r == yL], 
   DirichletCondition[p[x, r] == 0., x == xL]};

mesh = ToElementMesh[Ω, 
   MaxCellMeasure -> {"Length" -> 0.0012}];

{xVel, yVel, pressure} = 
  NDSolveValue[{op == {0, 0, 0}, Subscript[Γ, D]}, {u, 
    v, p}, Element[{x, r}, mesh], 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

ContourPlot[xVel[x, r], {x, 0, xL/100}, {r, 0, yL}, 
 AspectRatio -> 0.5, ColorFunction -> "TemperatureMap", 
 PlotLegends -> Automatic]