# Couple a PDE and ODE in NDSolve

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.

• Did I answer your question? Commented Aug 5, 2016 at 14:33
• @Young : Yes, your answer works for me, thank you. Commented Aug 8, 2016 at 8:35

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]


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]


• 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? Commented Aug 4, 2016 at 1:01
• @xzczd The x is ignored because the the reservoir response is 0-D Commented Aug 4, 2016 at 1:02
• @xzczd Said another way, T0 is only valid at x=0 Commented Aug 4, 2016 at 1:14
• @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) Commented Aug 4, 2016 at 2:50
• 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.) Commented Aug 4, 2016 at 3:48