Hi, I wrote the code to simulate the flow in a pipe. The model uses the Navier-Stokes equation and Mathematica solves them by means of finite element. At a first stage I use the cartesian coordinates and the code show good results, in the case of cylindrical coordinate there is something wrong. I report the first part of the code where the geometry are defined, the equation are written and also the first iteration. Please, tell me if there is something wrong.

Thank you.

\[CapitalOmega] = ImplicitRegion[True, {{r, Rin, Rout}, {z, 0, L}}];
Needs["NDSolve`FEM`"]
mesh = ToElementMesh[\[CapitalOmega], 
MaxCellMeasure -> {"Length" -> .02}, 
"MaxBoundaryCellMeasure" -> 0.0035, 
"MeshElementType" -> TriangleElement];
mesh["Wireframe"]



\[Mu] = 0.003;
R = 286;
T = 300;
eqc = 1/r \!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\((r*u[r, z])\)\) + \!\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\((v[r, z])\)\);
eqm1 = u[r, z]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(u[r, z]\)\) + v[r, z]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\(u[r, z]\)\) + 
   1/\[Rho][r, z]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(p[r, z]\)\) - \[Mu] (\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\((u[r, z])\)\)\) + 1/r*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\((u[r, z])\)\) + \!\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\((u[r, z])\)\)\) - u[r, z]/r^2);
eqm2 = u[r, z]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(v[r, z]\)\) + v[r, z] \!\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\(v[r, z]\)\) + 
   1/\[Rho][r, z]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\(p[r, z]\)\) - \[Mu] (\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\((v[r, z])\)\)\) + 1/r*\!\(
\*SubscriptBox[\(\[PartialD]\), \(r\)]\((v[r, z])\)\) + \!\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\(
\*SubscriptBox[\(\[PartialD]\), \(z\)]\((v[r, z])\)\)\));
NSOperator = {eqc, eqm1, eqm2};

bcs = {DirichletCondition[{u[r, z] == 0, 
     v[r, z] == 10}, {z == 0 || z == L}],
   DirichletCondition[{u[r, z] == 0, 
     v[r, z] == 0}, {r == Rout || r == Rin}], 
   DirichletCondition[p[r, z] == 101325, r == Rout && z == 0]};

StokesOperator = 
  NSOperator /. {u[r, z] -> 0, v[r, z] -> 0, \[Rho][r, z] -> 1.18, 
\!\(\*SuperscriptBox[\(\[Rho]\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[r, z] -> 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[r, z]/(R*T), 
\!\(\*SuperscriptBox[\(\[Rho]\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[r, z] -> 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[r, z]/(R*T)};
pde = StokesOperator == {0, 0, 0};
{u0, v0, p0} = 
  NDSolveValue[{pde, bcs}, {u, v, p}, {r, z} \[Element] mesh, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, 
     "IntegrationOrder" -> 5}];


ContourPlot[u0[r, z], {r, z} \[Element] mesh, 
 AspectRatio -> Automatic, PlotRange -> All, 
 ColorFunction -> ColorData["TemperatureMap"], Contours -> 10, 
 ImageSize -> Large]