2
$\begingroup$

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?

$\endgroup$
4
  • $\begingroup$ What is the physical meaning of the boundary condition cc1 ? Left side seems to be convection, right side is unclear to me. $\endgroup$ Aug 9, 2023 at 8:46
  • $\begingroup$ cc1 represents convection boundary condition. Letf side is convenction heat transfer which must be equal to conduction heat transfer in the external boundary. $\endgroup$ Aug 9, 2023 at 12:09
  • $\begingroup$ In my view Very unusual! Look for NeumannValue! $\endgroup$ Aug 9, 2023 at 16:20
  • $\begingroup$ 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. $\endgroup$
    – Roland F
    Aug 10, 2023 at 8:14

2 Answers 2

3
$\begingroup$

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}]

enter image description here

$\endgroup$
5
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ I'm missing boundarycondition cc1 in your solution. $\endgroup$ Aug 9, 2023 at 8:50
  • $\begingroup$ 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] $\endgroup$ Aug 9, 2023 at 12:14
  • $\begingroup$ 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)$$ $\endgroup$
    – Roland F
    Aug 11, 2023 at 10:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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