I have not checked @Alex Trounev's answer, but this answer shows that there is good agreement between _Mathematica_ and COMSOL Multiphysics.

Since you have a variety of thicknesses, I create a little routine so that I could mesh each region with the same number of elements (100 each).

    Needs["NDSolve`FEM`"]
    (* User Supplied Parameters *)
    g = {0.25, 0.114, 0.04};(*thickness*)
    gw = {0}~Join~Accumulate[g];
    λ = {8, 1.8, 44};
    ρ = {3100, 2100, 7800};
    cp = {1050, 1100, 540};
    (* Create a Multiregion Mesh *)
    ClearAll[seg, appendCrdRight]
    seg[thick_, nelm_, marker_] := Module[{crd, inc, marks},
      crd = Subdivide[0, thick, nelm];
      inc = Partition[Range[crd // Length], 2, 1];
      marks = ConstantArray[marker, inc // Length];
      <|"c" -> crd, "i" -> inc, "m" -> marks|>
      ]
    appendCrdRight[a1_, a2_] := Module[{crd, inc, marks, len, lcrd},
      len = a1["c"] // Length;
      lcrd = a1["c"] // Last;
      inc = Join[a1["i"], a2["i"] + len - 1];
      crd = Join[a1["c"], Rest[a2["c"] + lcrd]];
      marks = Join[a1["m"], a2["m"]];
      <|"c" -> crd, "i" -> inc, "m" -> marks|>]
    a = Fold[appendCrdRight, MapIndexed[seg[#1, 100, First[#2]] &, g]];
    mesh = ToElementMesh["Coordinates" -> Partition[a["c"], 1], 
       "MeshElements" -> {LineElement[a["i"], a["m"]]}, 
       "BoundaryElements" -> {PointElement[{{1}, {a["c"] // Length}}, {1, 
           2}]}];
    Show[mesh["Wireframe"["MeshElementStyle" -> {Red, Green, Blue}]], 
     PlotRange -> {-0.01, 0.01}]

[![Multiregion Mesh][1]][1]

Now, we can set up the PDE system and solve it on our newly created mesh.

    σ = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
    Trob = 1700;
    Tamb = 297;
    h = 10;
    ε = 0.85;
    rhocp = Evaluate[
       Piecewise[{{ρ[[1]] cp[[1]], gw[[1]] <= x <= gw[[2]]},
         {ρ[[2]] cp[[2]], gw[[2]] <= x <= gw[[3]]},
         {ρ[[3]] cp[[3]], gw[[3]] <= x <= gw[[4]]}}]];
    k = Evaluate[Piecewise[{{λ[[1]], gw[[1]] <= x <= gw[[2]]},
         {λ[[2]], gw[[2]] <= x <= gw[[3]]},
         {λ[[3]], gw[[3]] <= x <= gw[[4]]}}]];
    bc1 = DirichletCondition[T[t, x] == Trob, x == 0];
    bc2conv = NeumannValue[h*(Tamb - T[t, x]), x == Last@gw];
    bc2rad = NeumannValue[ε*σ*(Tamb^4 - T[t, x]^4), 
       x == Last@gw];
    ic1 = T[0, x] == Tamb;
    op = Inactive[Div][{{-k}}.Inactive[Grad][T[t, x], {x}], {x}] + 
       rhocp*Derivative[1, 0][T][t, x];
    pde = op == bc2conv + bc2rad;
    sol = NDSolveValue[{pde, bc1, ic1}, 
       T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];

The model I set up in COMSOL Multiphysics (v 5.5) shows similar results to those shown in the OP.

[![COMSOL Results][2]][2]

For comparison purposes, I extracted the temperature data at each phase boundary point in COMSOL.

[![Temperature Data at Phase Boundaries][3]][3]

I exported these data to compare versus the _Mathematica_ solution.

    data = {{0, 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, 
        10000, 11000, 12000, 13000, 14000, 15000, 16000, 17000, 18000, 
        19000, 20000, 21000, 22000, 23000, 24000, 25000, 26000, 27000, 
        28000, 29000, 30000, 31000, 32000, 33000, 34000, 35000, 
        36000}, {1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 
        1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 
        1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 1700, 
        1700, 1700, 1700, 1700, 1700, 1700}, {297, 297.9169787`, 
        320.0562147`, 374.4552427`, 444.9013611`, 517.6131837`, 
        587.4876631`, 652.6604327`, 712.3644688`, 766.9603206`, 
        816.5391802`, 861.866491`, 902.8730203`, 940.4564489`, 
        974.5556695`, 1005.867455`, 1034.417079`, 1060.665637`, 
        1084.866141`, 1107.411419`, 1128.099762`, 1146.931167`, 
        1164.637928`, 1180.832645`, 1195.499525`, 1208.917884`, 
        1221.536363`, 1233.003818`, 1243.320249`, 1252.972747`, 
        1261.872597`, 1269.909554`, 1277.155111`, 1284.007597`, 
        1290.216067`, 1295.780522`, 1300.901468`}, {297, 297.0000101`, 
        297.0108185`, 297.2403045`, 298.3422144`, 301.3296677`, 
        306.8304462`, 315.0786727`, 326.0187665`, 339.0198185`, 
        353.9950315`, 370.1369655`, 387.5159699`, 405.1722292`, 
        423.1836315`, 440.8382141`, 458.14222`, 474.6735528`, 
        490.3439464`, 504.9171794`, 518.5145476`, 531.1360512`, 
        542.7808248`, 553.4493263`, 563.1415743`, 571.9455027`, 
        580.0023514`, 587.2015743`, 593.5431713`, 599.3724133`, 
        604.6264161`, 609.2270331`, 613.2390417`, 617.0233547`, 
        620.3526001`, 623.2267777`, 625.8287217`}, {297, 297.0000065`, 
        297.0084849`, 297.2058139`, 298.1991325`, 300.9831864`, 
        306.2034638`, 314.1201414`, 324.7019404`, 337.3400768`, 
        351.9481631`, 367.722907`, 384.7337123`, 402.0228897`, 
        419.6676093`, 436.9560503`, 453.8952359`, 470.0643493`, 
        485.3780489`, 499.6031165`, 512.8593059`, 525.1466173`, 
        536.4765686`, 546.8430665`, 556.2458626`, 564.7760878`, 
        572.5801167`, 579.5433842`, 585.6658902`, 591.2927421`, 
        596.3610853`, 600.7928104`, 604.6517643`, 608.293677`, 
        611.4944415`, 614.2540579`, 616.7511966`}};
    Show[Plot[Evaluate[sol[t, #] & /@ gw], {t, 0, 36000}], 
     ListPlot[data[[2 ;; -1]], DataRange -> {0, 36000}]]

[![COMSOL Mathematica Comparison][4]][4]

As you can see, there is very little difference between COMSOL (dots) and _Mathematica_ (solid lines).


  [1]: https://i.sstatic.net/MRxDC.png
  [2]: https://i.sstatic.net/CBfXO.png
  [3]: https://i.sstatic.net/uaDCh.png
  [4]: https://i.sstatic.net/LwowR.png