1
$\begingroup$

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", 
    0, \[Rho]*Cp*Velocity.Inactive[Grad][T, X]];
  Q = If[Source === "NoSource", 0, Source];
  If[FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
  If[VectorQ[a], a = DiagonalMatrix[a]];
  (*Note the-sign in the operator*)
  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["NDSolve`FEM`"]
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]] 
$\endgroup$
5
  • $\begingroup$ Your heat capacity for polystyrene probably should be 1250 in SI units, but there still appears to be issues with the generated mesh. $\endgroup$ – Tim Laska Aug 3 '20 at 18:02
  • 2
    $\begingroup$ 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. $\endgroup$ – Tim Laska Aug 3 '20 at 19:22
  • $\begingroup$ 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 $\endgroup$ – kmulc Aug 3 '20 at 20:59
  • $\begingroup$ 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. $\endgroup$ – Tim Laska Aug 3 '20 at 23:31
  • $\begingroup$ I got the blue parts at the beginning to go away by changing the mesh size! $\endgroup$ – kmulc Aug 4 '20 at 14:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.