lc = 0.03;
rc = 0.01;
xp = 0.01;
c = 0.005;
rp = rc - c;
lp = lc - xp;
Subscript[T, 0] = 300;
Subscript[[Eta], 0] = 1.846*10^-5;
Subscript[P, 1] = 6*10^5 ;
Subscript[P, 0] = 10^5;
Subscript[c, P] = 1004.9;
Subscript[c, [Nu]] = 717.8;
Subscript[R, 0] = Subscript[c, P] - Subscript[c, [Nu]];
[CapitalOmega] = RegionDifference[
   Rectangle[{0, 0}, {lc, rc}],
   Rectangle[{xp, 0}, {xp + lp, rp}]];

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

eqn1 = D[[Rho][x, r]*Subscript[[Nu], x][x, r], x] +
    D[r*[Rho][x, r]*Subscript[[Nu], r][x, r], r]/r == 0 ;
eqn2 = D[[Rho][x, r]*Subscript[[Nu], x][x, r]^2 +
      Subscript[R, 0] [Rho][x, r]*T[x, r], x] +
    D[r*([Rho][x, r]*Subscript[[Nu], x][x, r]*
          Subscript[[Nu], r][x, r] +
         Subscript[[Eta], 0]*D[Subscript[[Nu], x][x, r], r]), r]/
     r == 0 ;
eqn3 = D[[Rho][x, r]*Subscript[[Nu], x][x, r]*
       Subscript[[Nu], r][x, r] +
      Subscript[[Eta], 0]*D[Subscript[[Nu], r][x, r], x], x] +
    D[r*([Rho][x, r]*Subscript[[Nu], r][x, r]^2 +
         Subscript[R, 0] [Rho][x, r]*T[x, r]), r]/r == 0;
eqn4 = Subscript[
     c, [Nu]]*[Rho][x,
      r]*(Subscript[[Nu], x][x, r]*D[T[x, r], x] +
       Subscript[[Nu], r][x, r]*D[T[x, r], r]) +
    Subscript[R, 0]*[Rho][x, r]*
     T[x, r]*(D[Subscript[[Nu], x][x, r], x] +
       D[r*Subscript[[Nu], r][x, r], x]/r) + (2*
        D[Subscript[[Nu], x][x, r], x]^2 +
       2*D[Subscript[[Nu], r][x, r],
          r]^2 + (D[Subscript[[Nu], x][x, r], r] +

          D[Subscript[[Nu], r][x, r],
           x])^2 - ((D[Subscript[[Nu], x][x, r], x] +
            D[r*Subscript[[Nu], r][x, r], x]/r)^2)*2/3)*
     Subscript[[Eta], 0] == 0;

bc1 = Subscript[R, 0] [Rho][0, r]*Subscript[T, 0] == Subscript[P, 1]
bc2 = Subscript[R, 0] [Rho][lc, r]*Subscript[T, 0] == Subscript[P, 0]
bc3 = DirichletCondition[{Subscript[[Nu], r][x, 0] == 0,
   D[Subscript[[Nu], r][x, r], r] == 0,
   D[Subscript[[Nu], x][x, r], r] == 0, D[[Rho][x, r], r] == 0,
   D[T[x, r], r] == 0}, r == 0 && (0 <= x <= xp  )]
bc4 = DirichletCondition[{Subscript[[Nu], r][x, r] == 0,
    Subscript[[Nu], x][x, r] ==
     0}, (0 <= r <= rp && x == xp ) || (r == rp &&
      xp <= x <= xp + lp)  || (r == rc && 0 <= x <= lc) ] == 0
bcs = {bc1, bc2, bc3, bc4};

{[Nu]xsol, [Nu]rsol, Tsol, [Rho]sol} =
  NDSolveValue[{eqns, , bcs}, {Subscript[[Nu], x], Subscript[[Nu],
    r], [Rho], T}, {x, r} [Element] mesh,
   Method -> {"FiniteElement",
     "InterpolationOrder" -> {Subscript[[Nu], x] -> 2,
       Subscript[[Nu], r] -> 2, [Rho] -> 1, T -> 1},
     "IntegrationOrder" -> 5}];