BarLegend inconsistent with DensityPlot colors

Question

In this post, 2D transient heat equation solution, @Alex Trounev helped me, in setting up a 2D simulation of a heat transfer problem where multiple materials are involved. The code is here:

Needs["NDSolve`FEM`"]


(*mesh*)
fibdiam = 0.06;   
fibrad = fibdiam/2;
time1=0.8;

reg1 = Disk[{0, 0}, fibrad];   
reg2 = RegionDifference[
   Rectangle[{-2*fibrad, 0}, {2*fibrad, 2*fibrad}], reg1];  
reg = RegionUnion[Rectangle[{-2*fibrad, 0}, {2*fibrad, 2*fibrad}], 
   reg1];  
mesh1 = ToElementMesh[reg, MaxCellMeasure -> fibdiam/10000];

(*material parameters*)
kf = 0.33; rhof = 940; Cpf = 2100;
kr = 0.22; rhor = 1300; Cpr = 1300;

k[x_, y_] := (10^6)*If[Element[{x, y}, reg2], kf, kr];
rho[x_, y_] := If[Element[{x, y}, reg2], rhof, rhor];
cp[x_, y_] := If[Element[{x, y}, reg2], Cpf, Cpr];

(*initial condition*)

TfibNW = 40;
TairNW = 40;
Tmelt = 220;

T01[x_, y_] := If[Element[{x, y}, reg2], Tmelt,If[Element[{x, y}, reg1], TfibNW, TairNW]]

(*solving PDE*)

eq1 = cp[x, y] rho[x, y] D[T1[x, y, t], t] - 
   Div[k[x, y] Grad[T1[x, y, t], {x, y}], {x, y}];

ic = T1[x, y, 0] == T01[x, y];
bc1 = NeumannValue[(10^-3)*h*(T1[x, y, t] - TairNW)/kr, y <= 0];
bc2 = NeumannValue[(10^-3)*h*(T1[x, y, t] - TairNW)/kf, y >= 2*fibrad];

h = 10^3;

 sol1 = 
 NDSolve[{eq1 == bc1 + bc2, ic}, T1, 
  Element[{x, y}, mesh1], {t, 0, time1}];


(*2D plotting*)

Table[
 
 Show[DensityPlot[Evaluate[T1[x, y, t] /. sol1[[1]]], 
   Element[{x, y}, mesh1],
   ColorFunction -> "TemperatureMap",
   PlotLegends -> BarLegend[{Automatic, {0, 250}}],
   PlotRange -> {All, All, {0, 250}},
   AspectRatio -> Automatic,
   PlotPoints -> 100,
   MaxRecursion -> 2,
   ImageSize -> 250,
   PlotLabel -> Row[{"t = ", t}]],
  
  Graphics[Circle[{0, 0}, fibrad]]],
 
 {t, 0, time1, time1/5}]

This leads to the following 2D graphs:

What I noticed is that the scale bar and the actual values of T1[x, y, t] /. sol1[[1]]] are inconsistent. This is obvious when I plot a cross section of the DensityPlot that is temperature at a fixed x coordinate (eg x==0) as function of time.

Is this a bug? How to fix it?

For 1D temperature profiles vertical through the center of the disk:

(*1D Temp profiles*)

Show[Table[
  Plot[Evaluate[T1[0, y, t] /. sol1[[1]]], {y, -fibrad, 2*fibrad},
   PlotRange -> {{-2*fibrad, 2*fibrad}, All},
   GridLines -> {{-fibrad, fibrad}, None},
   GridLinesStyle -> Red,
   Frame -> True,
   Axes -> None],
  
  {t, 0, time1, time1/5}]]

the temperature, at long times (close to time1) should be 160 in all the mesh points, but this is not consistent with the corresponding 2D image (3rd image, 2nd row).

Answer 1

Actually, I found some help in this post: Why does PlotRange in DensityPlot have no effect. A way out is to rescale the ColorFunction:

Table[Show[
  DensityPlot[Evaluate[T1[x, y, t] /. sol1[[1]]], 
   Element[{x, y}, mesh1],
   PlotRange -> {Automatic, Automatic, {20, 230}},
   PlotLegends -> BarLegend[{Automatic, {20, 230}}],
   ColorFunctionScaling -> False,
   ColorFunction -> (ColorData["TemperatureMap"][Rescale[#, {20, 230}]] &),
   AspectRatio -> Automatic,
   PlotPoints -> 100,
   MaxRecursion -> 2,
   ImageSize -> 250,
   PlotLabel -> Row[{"t = ", t}]],
  
  Graphics[Circle[{0, 0}, fibrad]]], {t, 0, time1, time1/5}]

Other views are also appreciated.