# Why is there no heat being transfered in my program?

I have been working on this code and am new to this coding language. I have been trying to model the transfer of heat through a hot plastic being cooled in a water bath. For some reason I have not been able to see the plastic being cooled with my code. I have provided the code below, thank you for helping!

Immediately below is the code I used to setup the heat transfer equation and the properties of the water and the plastic.

ClearAll[HeatTransferModel]
HeatTransferModel[T_, X_List, k_, \[Rho]_, Cp_, Velocity_, Source_] :=
Module[{V, Q, a = k},
V = If[Velocity === "NoFlow",
Q = If[Source === "NoSource", 0, Source];
If[FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
If[VectorQ[a], a = DiagonalMatrix[a]];
a = PiecewiseExpand[Piecewise[{{-a, True}}]];
Inactive[Div][a.Inactive[Grad][T, X], X] + V - Q]
TimeHeatTransferModel[T_, TimeVar_, X_List, k_, \[Rho]_, Cp_,
Velocity_, Source_] := \[Rho]*Cp*D[T, {TimeVar, 1}] +
HeatTransferModel[T, X, k, \[Rho], Cp, Velocity, Source]
(* material properties *)
Subscript[\[Rho], water] =
QuantityMagnitude[ThermodynamicData["Water", "Density"]];
Subscript[Cp, water] =
QuantityMagnitude[
ThermodynamicData["Water", "IsobaricHeatCapacity"]];
Subscript[k, water] =
QuantityMagnitude[
ThermodynamicData["Water", "ThermalConductivity"]];

Subscript[\[Rho], polystyrene] = 1045;
Subscript[Cp, polystyrene] = 1.25;
Subscript[k, polystyrene] = 0.14;


This is where I created the mesh, I was helped out in another post.

(* mesh *)
Needs["NDSolveFEM"]
bathx = 30;
bathy = 30;
plastic = Rectangle[{10, 10}, {20, 20}];
waterbath =
RegionDifference[Rectangle[{0, 0}, {bathx, bathy}], plastic];
waterbathCoordinate = {1, 1};
plasticCoordinate = {15, 15};
markerColors = {Orange, Blue};
markerSpecification = {{plasticCoordinate, 1,
0.25}, {waterbathCoordinate, 2, 1.}};
mesh = ToElementMesh[waterbath, "RegionHoles" -> None,
"RegionMarker" -> markerSpecification]
GraphicsRow[{mesh[
"Wireframe"[
"MeshElementStyle" ->
Map[Directive[FaceForm[#], EdgeForm[]] &, markerColors]]],
mesh["Wireframe"[
"MeshElementStyle" -> Map[FaceForm[#] &, markerColors]]]}]



Here I begin to set up the initial conditions.

(* parameters and initial conditions *)
parameters = {\[Rho] ->
Piecewise[{{Subscript[\[Rho], polystyrene],
ElementMarker == 1}, {Subscript[\[Rho], water],
ElementMarker == 2}}],
Cp -> Piecewise[{{Subscript[Cp, polystyrene],
ElementMarker == 1}, {Subscript[Cp, water],
ElementMarker == 2}}],
k -> Piecewise[{{Subscript[k, polystyrene],
ElementMarker == 1}, {Subscript[k, water],
ElementMarker == 2}}]};

Subscript[T, cold] = 20;
Subscript[T, hot] = 200;
ic = {T[0, x, y] ==
Piecewise[{{Subscript[T, hot],
ElementMarker == 1}, {Subscript[T, cold],
ElementMarker == 2}}]};
(* air *)
Subscript[T, amb] = 20;
h = 500;
Subscript[\[CapitalGamma], convective] =
NeumannValue[h*(Subscript[T, amb] - T[t, x, y]), y == 30];
(* solving *)
tend = 30;
pde = {TimeHeatTransferModel[T[t, x, y], t, {x, y}, k, \[Rho], Cp,
"NoFlow", "NoSource"] == Subscript[\[CapitalGamma], convective],
ic} /. parameters;
measure =
AbsoluteTiming[
MaxMemoryUsed[
Monitor[Tfun =
NDSolveValue[pde, T, {t, 0, tend}, {x, y} \[Element] mesh,
EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])],
monitor]]/(1024.^2)];
Print["Time -> ", measure[[1]], "\nMemory -> ", measure[[2]]]


Now I am creating the visualization. I think below is where I have done something wrong, as another problem I have is getting blue spots inside my center rectangle.


(* visualization *)
TRange = {Subscript[T, cold], Subscript[T, hot]};
legendBar =
BarLegend[{"TemperatureMap", TRange},
Sequence[20,
LegendLabel -> Text[Style["[\[Degree]C]", Opacity[0.6]]]]];
options = {PlotRange -> TRange,
Sequence[ColorFunction -> ColorData[{"TemperatureMap", TRange}],
ContourStyle -> Opacity[0.1], ColorFunctionScaling -> False,
Contours -> 10, PlotPoints -> 41,
FrameLabel -> {"x", "y"},
PlotLabel -> Text[Style["Temperature Field: T(t,x,y)", 18]],
AspectRatio -> Automatic, ImageSize -> Medium,
Prolog -> {ColorData[{"TemperatureMap", TRange}][
Subscript[T, amb]],
Rectangle[Scaled[{0, 0}], Scaled[{1, 1}]]}]};

frames = Table[
Labeled[Show[
Legended[
ContourPlot[Tfun[t, x, y], {x, y} \[Element] mesh,
Evaluate[options]], legendBar]],
Text[Style[Row[{"Time", " = ", N[t, 5], " [s]"}], 14,
Opacity[0.6]]]], {t, PowerRange[0.1, tend, 1.25]}];
frames = Rasterize[#1, "Image", ImageResolution -> 80] & /@ frames;
ListAnimate[Sequence[frames, SaveDefinitions -> True]]

• Your heat capacity for polystyrene probably should be 1250 in SI units, but there still appears to be issues with the generated mesh. – Tim Laska Aug 3 '20 at 18:02
• You should also note that you are simulating the cooling of a 10mx10m plastic block for 30 seconds. Nothing in likely to occur on that time scale for that large of a block. You should reduce the size of your block or dramatically increase tend. – Tim Laska Aug 3 '20 at 19:22
• Thanks. I changed the dimensions to smaller and it works very well now. I am still having the problem where when the visualization starts, there are blue parts inside the plastic – kmulc Aug 3 '20 at 20:59
• There are some very steep gradients when the simulation starts. This can be a problem with transient simulations. You may need substantial refinement at the boundaries to mitigate the ringing. Alternatively, you could choose a start time after the start-up transient. – Tim Laska Aug 3 '20 at 23:31
• I got the blue parts at the beginning to go away by changing the mesh size! – kmulc Aug 4 '20 at 14:10