Solving Navier-Stokes equations for a steady-state compressible viscous flow in a 2D axisymmetric step

step.m

lc = 0.03;
rc = 0.01;
xp = 0.01;
c = 0.005;
rp = rc - c;
lp = lc - xp;
T0 = 300;
eta0 = 1.846*10^-5;
P1 = 6*10^5 ;
P0 = 10^5;
cP = 1004.9;
cnu = 717.8;
R0 = cP - cnu;

Omega = RegionDifference[Rectangle[{0, 0}, {lc, rc}], 
   Rectangle[{xp, 0}, {xp + lp, rp}]];

Needs["NDSolve`FEM`"];
mesh = ToElementMesh[Omega, "MaxBoundaryCellMeasure" -> 0.00001, 
   MaxCellMeasure -> {"Length" -> 0.0008}, 
   "MeshElementConstraint" -> 20, MeshQualityGoal -> "Maximal"];



eqn1 = D[rho[x, r]*nux[x, r], x] + D[r*rho[x, r]*nur[x, r], r]/r == 0;
eqn2 = D[rho[x, r]*nux[x, r]^2 + R0*rho[x, r]*T[x, r], x] + 
    D[r*(rho[x, r]*nux[x, r]*nur[x, r] + eta0*D[nux[x, r], r]), r]/
     r == 0;
eqn3 = D[rho[x, r]*nux[x, r]*nur[x, r] + eta0*D[nur[x, r], x], x] + 
    D[r*(rho[x, r]*nur[x, r]^2 + R0*rho[x, r]*T[x, r]), r]/r == 0;
eqn4 = cnu*
     rho[x, r]*(nux[x, r]*D[T[x, r], x] + nur[x, r]*D[T[x, r], r]) + 
    R0*rho[x, r]*
     T[x, r]*(D[nux[x, r], x] + 
       D[r*nur[x, r], x]/r) + (2*D[nux[x, r], x]^2 + 
       2*D[nur[x, r], r]^2 + (D[nux[x, r], r] + 
          D[nur[x, r], 
           x])^2 - ((D[nux[x, r], x] + D[r*nur[x, r], x]/r)^2)*2/3)*
     eta0 == 0;
eqns = {eqn1, eqn2, eqn3, eqn4};


bc1 = R0* rho[0, r]*T0 == P1;
bc2 = R0 *rho[lc, r]*T0 == P0;
bc3 = DirichletCondition[{nur[x, 0] == 0, D[nur[x, r], r] == 0, 
    D[nux[x, r], r] == 0, D[rho[x, r], r] == 0, D[T[x, r], r] == 0}, 
   r == 0 && (0 <= x <= xp)];
bc4 = DirichletCondition[{nur[x, r] == 0, 
    nux[x, r] == 
     0}, (0 <= r <= rp && x == xp) || (r == rp && 
      xp <= x <= xp + lp) || (r == rc && 0 <= x <= lc)];
bcs = {bc1, bc2, bc3, bc4};

{nuxsol, nursol, rhosol, Tsol} = 
 NDSolveValue[{eqns, bcs}, {nux, nur, rho, T}, {x, r} \[Element] mesh,
   Method -> {"FiniteElement", 
    "InterpolationOrder" -> {nux -> 2, nur -> 2, rho -> 1, T -> 1}, 
    "IntegrationOrder" -> 5}]