kf = 0.000654(*[kW/mK]*)
ks = 0.00025(*[kW/mK]*)
Prandtlf[T_] :=
 2.193*(10^(-7))*(T^4) - 6.062*(10^(-5))*(T^3) + 0.006767*(T^2) -
  0.3959*T + 12.65
Muf[T_] :=
 2.72*(10^(-11))*(T^4) - 7.411*(10^(-09))*(T^3) +
  8.259*(10^(-07))*(T^2) - 4.969*(10^(-05))*T + 0.001726(*[Ns/m^2] *)
ReynoldsG[T_] := ((0.5*983.3)/(0.00152053*0.71))*(0.0015/(Muf[T]))
NusseltP1[T_] :=
 3.22*((ReynoldsG[T]*Prandtlf[T])^(1/3)) +
  0.117*(ReynoldsG[T]^0.8)*Prandtlf[T]^0.4
hp1[T_] := NusseltP1[T]*kf/0.5;

cps[T_] := (161.88/57.89)*(1.766*
   Exp[-((T - 58.41)/0.6049)^2] + 1.925*Exp[-((T - 57.34)/1.365)^2] +
  2.038*Exp[-((T - 32.38)/20.57)^2] +
  0.6765*Exp[-((T - 39.4)/2.519)^2] +
  0.6589*Exp[-((T - 30.2)/16.44)^2] +
  0.9032*Exp[-((T - 41.79)/50.14)^2])

eps = 0.0000000001
sol = NDSolveValue[{(r + eps)*Rhos *cps[T[r, t]]*D[T[r, t], t] -
     ks* D[r*(D[T[r, t], r]), r] ==
    NeumannValue[0, r == 0] -
     NeumannValue[hp1[T[r, t]]*(T[r, t] - 60), r == 0.0015],
   T[r, 0] == 20}, T, {r, 0, 0.0015}, {t, 0, 100},
  Method -> {"MethodOfLines",
    "SpatialDiscretization" -> "FiniteElement"}]

Plot[{sol[0.0005, t], sol[0, t], sol[0.0013, t]}, {t, 0, 10},
 PlotRange -> All, AxesLabel -> Automatic,
 PlotLegends -> "Expressions", PlotStyle -> {Blue, Red, Green}]
Plot3D[sol[r, t], {r, 0, 0.0015}, {t, 0, 10}, PlotRange -> All,
 AxesLabel -> Automatic, PlotLegends -> "Expressions"]