# Spherical Heat Equation and Convection Boundary Conditions

I'm trying to solve transient heat conduction equation in spherical domain in Mathematica. I made the simplifying hypothesis that temperature varies only with time and radial coordinate. The code is:

rmax = 10^(-3);(*sphere external radius(m)*)
rmin = 0;(*center of sphere(m)*)
h = 10^4;(*heat transfer coefficient - W/m^2.K*)
\[Rho] = 8000;(*sphere density - kg/m^3*)
c = 502;(*sphere specific heat - J/kg.K*)
tinf = 25;(*external temperature - \[Degree]C*)
ti = 1200;(*Initial temperature - \[Degree]C*)
k = 50(*thermal conductivity - W/m.K*);
eq1 = D[temp[r, t], t] == k*(1/r^2)*D[r^2 D[temp[r, t], r], r];
cc1 = h*(tinf - temp[rmax, t]) (1 - E^(-10^4 t)) == -k *
Evaluate[D[temp[r, t], r] /. r -> rmax];
cc2 = Evaluate[D[temp[r, t], r] /. r -> rmin] == 0;
sol = NDSolve[{eq1, temp[r, 0] == ti, cc1, cc2},
temp, {t, 0, 10}, {r, rmin, rmax}];
Plot[Evaluate[temp[0.1, t] /. sol], {t, 0, 10}, PlotRange -> All];


where eq1 is heat equation and cc1 and cc2 are boundary conditions. The initial condition of uniform temperature (Ti=1200) is in NDSolve.

With this code I had the following errors:

Power::infy: Infinite expression 1/0.^2 encountered.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

So I changed rmin variable to not be zero, but a number small enough so that it represents the center of sphere without having indetermination.

rmax = 10^(-3);(*sphere external radius(m)*)
rmin = 10^(-6);(*center of sphere(m)*)
h = 10^4;(*W/m^2.K*)
\[Rho] = 8000;(*kg/m^3*)
c = 502;(*J/kg.K*)
tinf = 25;(*\[Degree]C*)
ti = 1200;(*\[Degree]C*)
k = 50;
eq1 = D[temp[r, t], t] == k*(1/r^2)*D[r^2 D[temp[r, t], r], r];
cc1 = h*(tinf - temp[rmax, t]) (1 - E^(-10^4 t)) == -k *
Evaluate[D[temp[r, t], r] /. r -> rmax];
cc2 = Evaluate[D[temp[r, t], r] /. r -> rmin] == 0;
sol = NDSolve[{eq1, temp[r, 0] == ti, cc1, cc2},
temp, {t, 0, 10}, {r, rmin, rmax}];
Plot[Evaluate[temp[0.1, t] /. sol], {t, 0, 10}, PlotRange -> All];


But now I have the following errors:

General::munfl: 1/4.51156*10^307 is too small to represent as a normalized machine number; precision may be lost.

General::munfl: 1/4.53138*10^307 is too small to represent as a normalized machine number; precision may be lost.

General::munfl: 1/4.56494*10^307 is too small to represent as a normalized machine number; precision may be lost.

General::stop: Further output of General::munfl will be suppressed during this calculation.

NDSolve::mxst: Maximum number of 70421 steps reached at the point t == 0.0007513588858766584.

Does anybody know how to set up this problem properly in Mathematica?

• What is the physical meaning of the boundary condition cc1 ? Left side seems to be convection, right side is unclear to me. Commented Aug 9, 2023 at 8:46
• cc1 represents convection boundary condition. Letf side is convenction heat transfer which must be equal to conduction heat transfer in the external boundary. Commented Aug 9, 2023 at 12:09
• In my view Very unusual! Look for NeumannValue! Commented Aug 9, 2023 at 16:20
• The general boundary conditon at r=1 is a T(t,1)+ b Derivative[0,1[]T]t,1]=0. This in effect generates a the appropriate linear combinations a_n sin(n \pi r) + b_n cos(n \pi r) for the particular solutions. But I did not want the wait if Mathematica can DSolve for that case, too. I was even astonished that 13.3.0, DSolve finds the exact Fourier series for Dirichlet condition. Commented Aug 10, 2023 at 8:14

Here a numerical solution approach (3D).

It's not necessary to set boundary condition cc2! Instead of boundary condition cc1 we'll define NeumannValue.

rmax = 10^(-3) ;(*sphere external radius(m)*)
rmin = 0;(*center of sphere(m)*)
h = 10^4;(*heat transfer coefficient-W/m^2.K*)
\[Rho] = 8000;(*sphere density-kg/m^3*)
c = 502;(*sphere specific heat-J/kg.K*)
tinf = 25;(*external temperature-\[Degree]C*)
ti = 1200;
k = 50(*thermal conductivity-W/m.K*);

tsim = 1  ;
reg = ImplicitRegion[x^2 + y^2 + z^2 <= rmax^2, {x, y, z}]
U = NDSolveValue[{\[Rho] c Derivative[1, 0, 0, 0][u][t, x, y, z] ==k Laplacian[u [t, x, y, z], {x, y, z}] + h(*h/(4/3 Pi rmax^2)*)
NeumannValue[tinf - u [t, x, y, z] ,True],
u[0, x, y, z] == ti}, u, {t, 0, tsim}, Element[{x, y, z}, reg]]

Plot3D[U[t, x, 0, 0], {t, 0, tsim }, {x, 0, rmax},AxesLabel -> {t, x, u},MeshFunctions -> {#1 &}, PlotRange -> {0, ti}]


The equation can be solved in general. Particular solutions are products of $$\cos (n r) /r$$ and $$e^{- p t}$$ factors, that yield zero derivative at $$r=0$$

 (D[#, t] - k (D[#, r, r] + 2/r D[#, r]) &)[  E^( -k t p^2) Cos[  n r]/r] // Simplify


All conditions together and the $$1/r$$ factor explicitely introduced to cancel the first derivative $$2/r\partial_r$$

     DSolve[{  D[T[t, r]/r, t] == k (  D[T[t, r]/r, r, r] + 2/r  D[T[t, r]/r, r]),
Derivative[0, 1][T][t, 0] == 0, T[t, 1] == T1, T[0, r] == T0},
T [t, r], {t, r}][[1]] /. {K[1] -> n}


$$T(t,r)\to \text{T1} + \frac{4 (\text{T0}-\text{T1})}{\pi \ r}{ \underset{n=0}{\overset{\infty }{\sum }}\frac{(-1)^n}{2n+1} \ e^{\frac{1}{2} (-k) \left(n+\frac{1}{2}\right)^2 \pi ^2 t} \ \cos \left(\left(n+\frac{1}{2}\right) \ \pi \ r\right)}$$

The same works with NDSolve and numerical factors. But dont use the physical SI units of diffusion constants, their absolute values produce underflow in the exponents.

• I'm missing boundarycondition cc1` in your solution. Commented Aug 9, 2023 at 8:50
• I think boundary condition T[t,1]==T1 should be something like h*(T1- T[t,1]) (1 - E^(-10^4 t)) == -k * Evaluate[D[T[r, t], r] /. r ->1] Commented Aug 9, 2023 at 12:14
• I was happy that DSolve can now yields Fourier sum solutions for a Dirichlet problem even i it needs some help to get rid of first oder derivatives. The particular solutions for mixed byoundary conditions are $$0==\left(\frac{\partial \text{\#1}}{\partial t}-\frac{\frac{\partial }{\partial r}\left(r^2 \frac{\partial \text{\#1}}{\partial r}\right)}{r^2}\&\right)\left(\frac{(\alpha -\beta ) \sin (k-k r)+\beta k \cos (k-k r)}{e^{k^2 t} (k r ((\beta -\alpha ) \cos (k)+\beta k \sin (k)))}\right)$$ Commented Aug 11, 2023 at 10:05