1
$\begingroup$

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:

enter image description here

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

enter image description here

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).

$\endgroup$

1 Answer 1

1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.