Incorporated @user21's comment.

Source Link

edited Jun 3, 2020 at 6:39

16.4k
1
35
58

(* 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};
σ = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
Trob = 1700;
Tamb = 297;
h = 10;
ε = 0.85;
bmesh = ToBoundaryMesh["Coordinates" -> Partition[gw, 1], 
  "BoundaryElements" -> {PointElement[{{1}, {2}, {3}, {4}}]}]; nrEle \
= 100; pt = Partition[gw, 2, 1]; mesh = 
 ToElementMesh[bmesh, 
  "RegionMarker" -> 
   Transpose[{Partition[(Mean /@ pt), 1], {1, 2, 3}, 
     Abs[Subtract @@@ pt]/nrEle}]]
Show[mesh["Wireframe"["MeshElementStyle" -> {Red, Green, Blue}]], 
 PlotRange -> {-0.01, 0.01}]
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;
(* Coefficient Form *)
op = Inactive[Div][{{-k}}.Inactive[Grad][T[t, x], {x}], {x}] + 
   rhocp*Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Tcoef = NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Coefficient form withAlex's thermal"Basic diffusivityForm" *)
op = Inactive[Div][{{-k/rhocp}}.Inactive[Grad][T[trhocp*D[T[t, x], {x}], {x}] +t] 
 - k D[T[t, Derivative[1x], 0][T][tx, x];
pde = op == bc2conv + bc2rad;
TalphainsideTbasic = 
  NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh,
    StartingStepSize -> 0.01];
(* Alex'sCoefficient "Basicform Form"with thermal diffusivity *)
opbc2conv = rhocp*D[T[tNeumannValue[h*(Tamb - T[t, x])/rhocp, t]x == Last@gw];
bc2rad = NeumannValue[ε*σ*(Tamb^4 - kT[t, D[T[tx]^4)/
     rhocp, x == Last@gw];
op = Inactive[Div][{{-k/rhocp}}.Inactive[Grad][T[t, x], {x}], {x}] + 
   Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
TbasicTalphainside = 
  NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh,
    StartingStepSize -> 0.01];
(* Plot Alex's "Basic Form" *)
Show[Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}], 
 ListPlot[data[[2 ;; -1]], DataRange -> {0, 36000}]]
(* Comparison of Methods *)
Show[Plot[Evaluate[Tcoef[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> ConstantArray[{Opacity[0.2], Thickness[0.015]}, 4]], 
 Plot[Evaluate[Talphainside[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> Dashed], 
 Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> DotDashed]]

(* 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};
σ = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
Trob = 1700;
Tamb = 297;
h = 10;
ε = 0.85;
bmesh = ToBoundaryMesh["Coordinates" -> Partition[gw, 1], 
  "BoundaryElements" -> {PointElement[{{1}, {2}, {3}, {4}}]}]; nrEle \
= 100; pt = Partition[gw, 2, 1]; mesh = 
 ToElementMesh[bmesh, 
  "RegionMarker" -> 
   Transpose[{Partition[(Mean /@ pt), 1], {1, 2, 3}, 
     Abs[Subtract @@@ pt]/nrEle}]]
Show[mesh["Wireframe"["MeshElementStyle" -> {Red, Green, Blue}]], 
 PlotRange -> {-0.01, 0.01}]
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;
(* Coefficient Form *)
op = Inactive[Div][{{-k}}.Inactive[Grad][T[t, x], {x}], {x}] + 
   rhocp*Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Tcoef = NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Coefficient form with thermal diffusivity *)
op = Inactive[Div][{{-k/rhocp}}.Inactive[Grad][T[t, x], {x}], {x}] + 
    Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Talphainside = 
  NDSolveValue[{pde, bc1, ic1}, T, {t, 0, 36000}, {x} ∈ mesh,
    StartingStepSize -> 0.01];
(* Alex's "Basic Form" *)
op = rhocp*D[T[t, x], t] - k D[T[t, x], x, x];
pde = op == bc2conv + bc2rad;
Tbasic = NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Plot Alex's "Basic Form" *)
Show[Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}], 
 ListPlot[data[[2 ;; -1]], DataRange -> {0, 36000}]]
(* Comparison of Methods *)
Show[Plot[Evaluate[Tcoef[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> ConstantArray[{Opacity[0.2], Thickness[0.015]}, 4]], 
 Plot[Evaluate[Talphainside[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> Dashed], 
 Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> DotDashed]]

(* 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};
σ = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
Trob = 1700;
Tamb = 297;
h = 10;
ε = 0.85;
bmesh = ToBoundaryMesh["Coordinates" -> Partition[gw, 1], 
  "BoundaryElements" -> {PointElement[{{1}, {2}, {3}, {4}}]}]; nrEle \
= 100; pt = Partition[gw, 2, 1]; mesh = 
 ToElementMesh[bmesh, 
  "RegionMarker" -> 
   Transpose[{Partition[(Mean /@ pt), 1], {1, 2, 3}, 
     Abs[Subtract @@@ pt]/nrEle}]]
Show[mesh["Wireframe"["MeshElementStyle" -> {Red, Green, Blue}]], 
 PlotRange -> {-0.01, 0.01}]
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;
(* Coefficient Form *)
op = Inactive[Div][{{-k}}.Inactive[Grad][T[t, x], {x}], {x}] + 
   rhocp*Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Tcoef = NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Alex's "Basic Form" *)
op = rhocp*D[T[t, x], t] - k D[T[t, x], x, x];
pde = op == bc2conv + bc2rad;
Tbasic = NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Coefficient form with thermal diffusivity *)
bc2conv = NeumannValue[h*(Tamb - T[t, x])/rhocp, x == Last@gw];
bc2rad = NeumannValue[ε*σ*(Tamb^4 - T[t, x]^4)/
     rhocp, x == Last@gw];
op = Inactive[Div][{{-k/rhocp}}.Inactive[Grad][T[t, x], {x}], {x}] + 
   Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Talphainside = 
  NDSolveValue[{pde, bc1, ic1}, T, {t, 0, 36000}, {x} ∈ mesh,
    StartingStepSize -> 0.01];
(* Plot Alex's "Basic Form" *)
Show[Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}], 
 ListPlot[data[[2 ;; -1]], DataRange -> {0, 36000}]]
(* Comparison of Methods *)
Show[Plot[Evaluate[Tcoef[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> ConstantArray[{Opacity[0.2], Thickness[0.015]}, 4]], 
 Plot[Evaluate[Talphainside[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> Dashed], 
 Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> DotDashed]]

Added Alex's request to compare the "basic form" to COMSOL.

Source Link

edited Jun 3, 2020 at 2:44

Tim Laska

16.4k
1
35
58

Update to Include the Basic Form

@AlexTrounev requested a comparison of the basic form to COMSOL as defined by:

$$\rho {{\hat C}_p}\frac{{\partial T}}{{\partial t}} - \lambda \frac{{{\partial ^2}T}}{{\partial {x^2}}} = 0$$

To use the FEM method, I recommend to cast your equations into coefficient form as shown FEM Tutorial.

$$\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$

I find it easier to make comparisons of commercial solver (such as COMSOL) results to Mathematica results.

As shown with the following workflow, the Alex's basic form also matches COMSOL quite closely. I also included a case where I tried to thermal diffusivity in coefficient form and it fails to match COMSOL. Finally, it may be interesting to note that COMSOL's Laplace Equation Interface does not contain a Laplacian, rather:

$$\nabla \cdot \left( { - \nabla u} \right) = 0$$

(* 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};
σ = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
Trob = 1700;
Tamb = 297;
h = 10;
ε = 0.85;
bmesh = ToBoundaryMesh["Coordinates" -> Partition[gw, 1], 
  "BoundaryElements" -> {PointElement[{{1}, {2}, {3}, {4}}]}]; nrEle \
= 100; pt = Partition[gw, 2, 1]; mesh = 
 ToElementMesh[bmesh, 
  "RegionMarker" -> 
   Transpose[{Partition[(Mean /@ pt), 1], {1, 2, 3}, 
     Abs[Subtract @@@ pt]/nrEle}]]
Show[mesh["Wireframe"["MeshElementStyle" -> {Red, Green, Blue}]], 
 PlotRange -> {-0.01, 0.01}]
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;
(* Coefficient Form *)
op = Inactive[Div][{{-k}}.Inactive[Grad][T[t, x], {x}], {x}] + 
   rhocp*Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Tcoef = NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Coefficient form with thermal diffusivity *)
op = Inactive[Div][{{-k/rhocp}}.Inactive[Grad][T[t, x], {x}], {x}] + 
   Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Talphainside = 
  NDSolveValue[{pde, bc1, ic1}, T, {t, 0, 36000}, {x} ∈ mesh,
    StartingStepSize -> 0.01];
(* Alex's "Basic Form" *)
op = rhocp*D[T[t, x], t] - k D[T[t, x], x, x];
pde = op == bc2conv + bc2rad;
Tbasic = NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Plot Alex's "Basic Form" *)
Show[Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}], 
 ListPlot[data[[2 ;; -1]], DataRange -> {0, 36000}]]
(* Comparison of Methods *)
Show[Plot[Evaluate[Tcoef[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> ConstantArray[{Opacity[0.2], Thickness[0.015]}, 4]], 
 Plot[Evaluate[Talphainside[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> Dashed], 
 Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> DotDashed]]

Update to Include the Basic Form

@AlexTrounev requested a comparison of the basic form to COMSOL as defined by:

$$\rho {{\hat C}_p}\frac{{\partial T}}{{\partial t}} - \lambda \frac{{{\partial ^2}T}}{{\partial {x^2}}} = 0$$

To use the FEM method, I recommend to cast your equations into coefficient form as shown FEM Tutorial.

$$\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$

I find it easier to make comparisons of commercial solver (such as COMSOL) results to Mathematica results.

As shown with the following workflow, the Alex's basic form also matches COMSOL quite closely. I also included a case where I tried to thermal diffusivity in coefficient form and it fails to match COMSOL. Finally, it may be interesting to note that COMSOL's Laplace Equation Interface does not contain a Laplacian, rather:

$$\nabla \cdot \left( { - \nabla u} \right) = 0$$

(* 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};
σ = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
Trob = 1700;
Tamb = 297;
h = 10;
ε = 0.85;
bmesh = ToBoundaryMesh["Coordinates" -> Partition[gw, 1], 
  "BoundaryElements" -> {PointElement[{{1}, {2}, {3}, {4}}]}]; nrEle \
= 100; pt = Partition[gw, 2, 1]; mesh = 
 ToElementMesh[bmesh, 
  "RegionMarker" -> 
   Transpose[{Partition[(Mean /@ pt), 1], {1, 2, 3}, 
     Abs[Subtract @@@ pt]/nrEle}]]
Show[mesh["Wireframe"["MeshElementStyle" -> {Red, Green, Blue}]], 
 PlotRange -> {-0.01, 0.01}]
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;
(* Coefficient Form *)
op = Inactive[Div][{{-k}}.Inactive[Grad][T[t, x], {x}], {x}] + 
   rhocp*Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Tcoef = NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Coefficient form with thermal diffusivity *)
op = Inactive[Div][{{-k/rhocp}}.Inactive[Grad][T[t, x], {x}], {x}] + 
   Derivative[1, 0][T][t, x];
pde = op == bc2conv + bc2rad;
Talphainside = 
  NDSolveValue[{pde, bc1, ic1}, T, {t, 0, 36000}, {x} ∈ mesh,
    StartingStepSize -> 0.01];
(* Alex's "Basic Form" *)
op = rhocp*D[T[t, x], t] - k D[T[t, x], x, x];
pde = op == bc2conv + bc2rad;
Tbasic = NDSolveValue[{pde, bc1, ic1}, 
   T, {t, 0, 36000}, {x} ∈ mesh, StartingStepSize -> 0.01];
(* Plot Alex's "Basic Form" *)
Show[Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}], 
 ListPlot[data[[2 ;; -1]], DataRange -> {0, 36000}]]
(* Comparison of Methods *)
Show[Plot[Evaluate[Tcoef[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> ConstantArray[{Opacity[0.2], Thickness[0.015]}, 4]], 
 Plot[Evaluate[Talphainside[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> Dashed], 
 Plot[Evaluate[Tbasic[t, #] & /@ gw], {t, 0, 36000}, 
  PlotStyle -> DotDashed]]

Source Link

created Jun 1, 2020 at 3:35

Tim Laska

16.4k
1
35
58

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}]

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.

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

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}]]

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

Stack Exchange Network

Return to Answer

Update to Include the Basic Form

Update to Include the Basic Form