17
$\begingroup$

I would like to solve an example of non-stationary heat transfer with a coupled PDE and ODE. Let's assume that we have 1 dimensional bar of length $L$ with uniform initial temperature. The right end is isolated (no heat flux), and the left end is attached to a (small) isolated reservoir of hot water. The water in the reservoir is well mixed, and we assume it has 0 spatial dimensions. On the boundary between the reservoir and the bar, heat is exchanged by convection, with a known constant coefficient. After enough time has passed, we get uniform temperature in the bar and reservoir, somewhere between both initial temperatures.

Question: How to solve this example with NDSolve command or some similar elegant piece of code?

I have a problem with coupling a PDE ($T[x,t]$) that describes heat transfer inside the bar, and the ODE ($T0[t]$) that describes cooling of the reservoir. It seems that NDSolve can't handle two equations with different number of variables. This is my solution so far, where I solve both equations over the whole domain, but I omit the diffusion term for the reservoir. I think qualitatively the temperature evolution in the bar is appropriate, but the values are wrong.

L = 1.; (* length of domain *)
end = 20.; (* end time *)
k = 2; (* convection coefficient *)
m0 = 100.; (* mass of reservoir *)
Tinit = 0. ;(* initial temperature of the bar *)
T0init = 100.; (* initial temperature of the fluid *)

(* Tsol ... temperature in the bar, T0sol ... temperature on the boundary *)
(* Constants,such as heat conduction coefficient inside the bar, cross-section,etc are omitted for simplicity (assumed to be 1).*)
{Tsol, T0sol} =
NDSolveValue[{
  D[T[t, x], t] - D[T[t, x], x, x] == NeumannValue[-k (T0[t, x] - T[t, x]), x == 0] + NeumannValue[0., x == L]
   , m0*D[T0[t, x], t] == NeumannValue[k (T0[t, x] - T[t, x]), x == 0]
   , T[0, x] == Tinit
   , T0[0, x] == T0init
   }, {T, T0}
  , {t, 0, end}
  , {x, 0, L}
  , Method -> {"MethodOfLines", "SpatialDiscretization" -> "FiniteElement"}
  ]

Manipulate[
 Plot[{Tsol[t, x], T0sol[t, x]}, {x, 0, L}
, PlotRange -> {Automatic, {-1, 110}}
, AxesLabel -> {"x", "Temperature"}
, Epilog -> {Red, PointSize[0.015], Point[{0, T0sol[t, 0]}]}
, PlotLegends -> {"bar", "reservoir"}]
, {t, 0, end}]

Bonus question: How can I extend the solution for another dimension? For example, we have a rectangular domain (the bar), isolated on the right side, with fixed temperature on top and bottom and 1D fluid flow on the left side (advection equation), where heat is exchanged by convection.

$\endgroup$
  • $\begingroup$ Did I answer your question? $\endgroup$ – Young Aug 5 '16 at 14:33
  • $\begingroup$ @Young : Yes, your answer works for me, thank you. $\endgroup$ – Pinti Aug 8 '16 at 8:35
18
+100
$\begingroup$

Modified Newton Cooling Law implementation and reduced FEA cell length

Clear["Global`*"]

len = 1;(*length of domain*)
end = 5;(*end time*)
h = 2;(*convection coefficient*)
cpmR = 9;(*heat capacity/mass ratio*)
Tinit = 0;(*initial temperature of the bar*)
T0init = 100;(*initial temperature of the fluid*)

(*Tsol=temperature in the bar,T0sol=temperature on the boundary*)
(*Constants,such as heat conduction coefficient inside the
bar,cross-section,etc are omitted for simplicity (assumed to be 1)*)

{Tsol, T0sol} = NDSolveValue[{
    D[T[t, x], t] - D[T[t, x], x, x] == 
    NeumannValue[h (T[t, x] - T0[t, x]), x == 0] + NeumannValue[0, x == len],
    cpmR D[T0[t, x], t] == -h (T0[t, x] - T[t, x]),
   T[0, x] == Tinit, T0[0, x] == T0init},
  {T, T0}, {t, 0, end}, {x, 0, len},
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 1*^-4}}}}];

result = Table[Plot[{Tsol[t, x], T0sol[t, 0]}, {x, 0, len},
    PlotRange -> {{0, 1}, {0, T0init}}, 
    FrameLabel -> {"x", "Temperature"}, 
    PlotLabels -> {"Bar", "Reservoir"},
    PlotLabel -> Row@{"Heat Capacity/Mass Ratio = ", N[cpmR], "     ", 
                      "Time = ", PaddedForm[N[t, 3], {3, 2}]},
    BaseStyle -> {FontFamily -> Times, FontSize -> 12}, Frame -> True],
   {t, 0, end, end/100}];
ListAnimate[result]

enter image description here

enter image description here

2D Answer:

Clear["Global`*"]

len = 1;(*length of domain*)
wid = 1;(*width of domain*)

end = 3;(*end time*)

h = 2;(*convection coefficient*)
cpmR = 9;(*heat capacity/mass ratio*)

Tinit = 0;(*initial temperature of the bar*)
T0init = 100;(*initial temperature of the fluid*)

Note: Constants, such as thermal conductivity, heat capacity, mass, cross-section, etc are all omitted for simplicity (assumed to be 1)

{Tsol2D, T0sol2D} = NDSolveValue[{
    D[T[t, x, y], t] - (D[T[t, x, y], x, x] + D[T[t, x, y], y, y]) == 
     NeumannValue[-h (T[t, x, y] - T0[t, x, y]), x == 0] + 
      NeumannValue[0, x == len],
    cpmR D[T0[t, x, y], t] == -h (T0[t, x, y] - T[t, x, y]),
    T[0, x, y] == Tinit, T0[0, x, y] == T0init},
   {T, T0}, {t, 0, end}, {x, 0, len}, {y, 0, wid},
   Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}];

result2D = 
  Table[Show[
    Plot3D[Tsol2D[t, x, y], {x, 0, len}, {y, 0, wid}, 
     PlotRange -> {{0, len}, {0, wid}, {0, T0init}},
     ColorFunction -> Function[{x, y, z}, ColorData["TemperatureMap"][z/T0init]], 
     ColorFunctionScaling -> False,
     MeshFunctions -> {#3 &}, Mesh -> {Range[5, T0init, 5]}, 
     AxesLabel -> Automatic, 
     PlotLabel ->  Row@{"Heat Capacity/Mass Ratio = ", N[cpmR], "     ", 
                        "Time = ", PaddedForm[N[t, 3], {3, 2}]},
     BaseStyle -> {FontFamily -> Times, FontSize -> 12}], 
    Graphics3D[{Red, 
      Line[{{0, 0, T0sol2D[t, 0, 0]}, {0, 1, T0sol2D[t, 0, 0]}}]}]],
   {t, 1/10, end, end/100}];
ListAnimate[result2D, AnimationRate -> 3]

enter image description here

$\endgroup$
  • $\begingroup$ The reservoir is just attached to the left end, are you sure the equation mR D[T0[t, x], t] == -h (T0[t, x] - T[t, x]) is correct? $\endgroup$ – xzczd Aug 4 '16 at 1:01
  • 1
    $\begingroup$ @xzczd The x is ignored because the the reservoir response is 0-D $\endgroup$ – Young Aug 4 '16 at 1:02
  • 1
    $\begingroup$ @xzczd Said another way, T0 is only valid at x=0 $\endgroup$ – Young Aug 4 '16 at 1:14
  • 2
    $\begingroup$ @xzczd I respectfully disagree ... the boundary condition is being applied to a 0-D space and and Q=cp m dT ... cpR here is a relative term for mass and cp ... so cpR dT/dt = q = -k(T0-T) $\endgroup$ – Young Aug 4 '16 at 2:50
  • 1
    $\begingroup$ You're right. (Seems that I'm a little tired…) Still, I think there should be a k before D[T[t, x], x, x]. (This mistake is inherited from OP's code, of course.) $\endgroup$ – xzczd Aug 4 '16 at 3:48

Your Answer

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

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