I am trying to solve the heat equation for a system made of a cylinder (hot) that is suddenly immersed in a cooling medium (cold). I am doing this using a 1D approximation with single differential equation (in cylindrical coordinates) and physical properties that vary in space. The code is:
diamcyl = 0.800;
T1 = 140;
T2 = -20;
k2 = 0.58; (* W/mK *)
rho2 = 1000; (* kg/m3 *)
Cp2 = 4200; (* J/KgK *)
k1 = 0.128; (*W/(mK)*)
rho1 = 800; (* kg/m3 *)
Cp1 = 1670; (*J/(KgK)*)
(*properties as a function of space*)
rho[x_] := (rho1 + (rho2 - rho1)* UnitStep[x - diamcyl/2]);
k[x_] := 10^6*(k1 + (k2 - k1)* UnitStep[x - diamcyl/2]);
Cp[x_] := (Cp1 + (Cp2 - Cp1)* UnitStep[x - diamcyl/2]);
heateq = 1/x*D[x*k[x]*D[u[x, t], x], x] == rho[x]*Cp[x]*D[u[x, t], t];
u0[x_] := (T1 + (T2 - T1)*UnitStep[x - diamcyl/2]);
solm20 = First[
NDSolve[{heateq, u[x, 0] == u0[x],
u[10, t] == T2, (D[u[x, t], x] /. x -> 0.001) == 0},
u, {x, 0.001, 2}, {t, 0, 5},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 1000}}]];
The solution transient temperature profile that I obtain looks smooth and is continuous over the boundary of the cylinder. However, the heat flux exhibits a discontinuity over this boundary.
Show[Table[
Plot[u[x, t] /. solm20, {x, 0.001, 2},
PlotRange -> {{0, 1}, {T2, T1}}, PlotRangePadding -> {None, 10},
Prolog -> {LightGray, Rectangle[{0, -100}, {diamcyl/2, 200}]},
GridLines -> {{diamcyl/2}, {T2}},
FrameLabel -> {"distance [mm]", "temperature [\[Degree]C]"}], {t,
0, 2, 0.05}],
Plot[u0[x], {x, 0.01, 2}, PlotRange -> All, Exclusions -> None,
PlotStyle -> Red]]
Show[Table[
Plot[Evaluate[k[x]*D[u[x, t], x] /. solm20], {x, 0.01, 2},
PlotRange -> {{0, 1}, All}, PlotRangePadding -> {None, 10},
GridLines -> {{diamcyl/2}, {T2}},
FrameLabel -> {"distance [mm]", "temperature [\[Degree]C]"}], {t,
0.01, 2, 0.05}]]
How to obtain a solution with continuous heat flux? Or, I could also ask it in this way: how to impose a Neumann condition over the boundary of the cylinder?
