Can Matlab give me an analytic solution for transient Heat diffusion?

I want an analytic solution to the transient Heat Diffusion equation in cylindrical coordinates. My equation is

dT/dt = d^2T/dr^2 +1/r dT/dr + d^2T/dz^2

This problem is cylindrically symmetric, so the theta term is zero.

Is there a way to solve this in MatLab? My prof tells me MatLab should be able to do this with the symbolic toolkit, but all I can find on it concerns ODEs.

The boundary conditions for this problem are T(0,z,t) = g(z), T(r,z,0)=0

[Shaihroz Khan]

Listen I don't know but u should better try to visualize your problem in fluent rather than wasting your time matlab, coz for that u have to program the code and you wont find it here.

Im also having the same problem but in momentum transfer for getting velocity profile in a cylinder with velocity being a function of r and z and time.

[Walter Roberson]

Shaihroz, this was from two years ago, and discussion stalled because the poster was unable to define the desired behaviour along the axis of the cylinder. If a behaviour had been defined we might have been able to program it.

Answer 1

[Walter Roberson]

Older versions of MATLAB used the Maple symbolic package, which has a PDE solver. I have not encountered a PDE solver for the MuPAD symbolic engine as yet, as partly because I do not have that Toolbox.

Maple believes that the only solution is if g(z) = 0, with T(r,z,t) = 0.

Maple is able to do a partial resolution if you omit the T(0,z,t)=g(z) initial condition. That is partly because solving PDE typically depends critically on details of all of the functions, but it is also because it encounters a singularity at r=0. The formula has a 1/r which at r=0 would give "undefined" leaving the rest of the formula undefined. When that first initial condition is omitted,

T(r, z, t) = -_C5 * _C1 * (_C3 * exp(_c[2]^(1/2) * z - t * RootOf(BesselJ(0, _Z * r) * BesselY(1, _Z * r) - BesselJ(1, _Z * r) * BesselY(0, _Z * r))^2 + t  * _c[2]) + _C4 * exp(-_c[2]^(1/2) * z - t * RootOf(BesselJ(0, _Z * r) * BesselY(1, _Z * r) - BesselJ(1, _Z * r)  * BesselY(0, _Z * r))^2 + t * _c[2])) * (-BesselY(1, RootOf(BesselJ(0, _Z * r) * BesselY(1, _Z * r) - BesselJ(1, _Z * r) * BesselY(0, _Z * r)) * r) * BesselJ(0, csgn(RootOf(BesselJ(0, _Z * r) * BesselY(1, _Z*r) - BesselJ(1, _Z*r) * BesselY(0, _Z*r))) * RootOf(BesselJ(0, _Z*r) * BesselY(1, _Z * r) - BesselJ(1, _Z*r) * BesselY(0, _Z*r)) * r) + BesselJ(1, RootOf(BesselJ(0, _Z*r) * BesselY(1, _Z * r) - BesselJ(1, _Z*r) * BesselY(0, _Z*r)) * r) * BesselY(0, csgn(RootOf(BesselJ(0, _Z*r) * BesselY(1, _Z * r) - BesselJ(1, _Z*r) * BesselY(0, _Z*r))) * RootOf(BesselJ(0, _Z*r) * BesselY(1, _Z*r) - BesselJ(1, _Z * r) * BesselY(0, _Z*r)) * r)) / BesselY(1, RootOf(BesselJ(0, _Z*r) * BesselY(1, _Z*r) - BesselJ(1, _Z*r) * BesselY(0, _Z*r)) * r)

This is Maple notation.

_C1, _c[2], _C3, _C4, and _C5 should all be read as the names of constants of integration that need to be resolved according to further initial conditions.

MuPAD uses lower-case names for the Bessel functions, such as bessely where Maple uses BesselY

The form RootOf(f(_Z)) in Maple translates into MuPAD as RootOf(f(z1), z1) . It stands for the set of values of the variable _Z such that f(_Z) is 0 -- the "roots" of the expression.

At r=0, BesselY(1, _Z*r) would be BesselY(1, 0) which is undefined. lim x => -0 BesselY(1, x) is +infinity and lim x <= +0 BesselY(1,x) is -infinity so there is an unresolvable singularity.

Perhaps T(0,z,t) = g(z) is not really a "boundary condition" but rather a definition overriding the behavior at r = 0 ?

[Joshua]

What do you mean by "a definition overriding the behavior at r = 0" ?

[Walter Roberson]

Suppose you define a function

y = 6/x

Then as you approach x = 0 from above, y goes to +infinity, but as you approach x = 0 from below, y goes to -infinity. So in this example function, y is well defined for every real number except x = 0 exactly, where it must be undefined, as it has a fundamental singularity there.

As opposed to the example

y = 6 / abs(x)

then as x approaches 0 from above , y goes to +infinity, an as x approaches 0 from below, y goes the +infinity. Right at x = 0 it has a singularity, but it is a "removable singularity", because if you were to _define y(0) as +infinity then there would be continuity as you flowed upward toward y = infinity as you approach x = 0 from below, then you stay at infinity at x = 0 exactly, then you come down from infinity as you get further positive away from x = 0.

Thus, if what you were to do would be to define your function as:

t(r,z,t) = piecewise( r == 0, -infinity, otherwise dT/dt = d^2T/dr^2 +1/r dT/dr + d^2T/dz^2 )

then the function would be well-defined at every point because at what would otherwise be the plane of singularity at r = 0 you have defined a result.

The difficulty with doing this is that you would still have a singularity at r = 0 because although the function is well defined for every r = -delta_R, delta_R > 0, and is well defined for every r = 0, and is well defined for every r = +delta_R, delta_R = -infinity, there is no continuity as you go from negative r to r = 0: you go from +infinity to -infinity.

What you can do is to say that the domain of the function is r >= 0, which causes the problems with negative r to vanish because the function would just not apply, just as it is not meaningful to "take the square root of a walnut". Then you are still stuck with what to do at r = 0 exactly. If you define T(0,z,t) = g(z) then you define a value for each (r,z,t), r >= 0, but unfortunately unless g(z) = -infinity always, then you still get a discontinuity because T(r,z,t) goes to -infinity as r approaches 0 from above and hence any finite value defined at r = 0 exactly would not be adjacent to the -infinity.

You need to rethink what should happen at r = 0 exactly. If r indicates a distance from the central axis of the cylinder, then g(z) is dependent only on the height and not on the time: is that truly reasonable? The value for every molecule touching those in the central core has a value that is dependent on height and time both; if any kind of solid is being modeled, one would expect those time-dependent temperatures to affect temperatures at the pole. Or is there a theorem that states that the inflow of heat from the adjacent molecules will always exactly balance the outflow of heat?