$$\newcommand{\D}{\displaystyle} \renewcommand{\eqref}[1]{Eq.~(\ref{#1})}$$

## 7.7.2 Example: Cooling of a sphere

Figure 102: Problem setup for example on the cooling of a sphere.

Figure 102 shows our problem setup, which consists of a sphere immersed in water and exchanging heat with it. The sphere radius is $$b=5$$ and its temperature before it is immersed in water is $$T_k$$. Moreover, we assume that the water keeps a constant temperature $$T_v$$ throughout the entire process and we neglect heat exchange with the environment.

• Heat conductivity: $$k=0.1W/(cm\cdot^\circ C)$$
• Heat transfer coefficient: $$\bar h =0.2W/(cm\cdot^\circ C)$$
• Thermal diffusivity: $$\alpha=0.04 cm^2/s$$
We have chosen the above reported coefficients so that $$\dfrac{\bar h \cdot b}{k}=1$$, since this choice leads to a simple analytical solution. Moreover, such values correspond to characteristic values found for nickel alloys.

We must now solve the following problem with $$T=T(r,t)$$: $$$$\tag{7.122} \frac{\partial T}{\partial t}=\alpha \left(\frac{\partial^2T}{\partial r^2}+\frac{2}{r}\frac{\partial T}{\partial r} \right),$$$$

with boundary conditions: $$$$\tag{7.123} \frac{\partial T}{\partial r}(0,t)=0,\ \text{(symmetribetingelse)}$$$$

For $$r=b:$$ $$$$\tag{7.124} k\frac{\partial T}{\partial r}= \bar h \cdot (T_v-T_b)$$$$

Initial conditions are: $$$$\tag{7.125} T(r,0)=T_k$$$$

In this example we use the FTCS scheme, but the derivation can be easily extended to the $$\theta$$-scheme as shown in Example 7.7.1 Example: Start-up flow in a tube. With $$r_j=\Delta r\cdot j,\ \Delta r= \frac{1}{N+1},\ j=0,1,\dots,N+1$$ and $$D=\alpha \frac{\Delta t}{(\Delta r)^2}$$ from (7.92) we get that when $$u(r,t)\to T(r,t),\ \theta=0$$ and $$\lambda = 2$$: $$$$\tag{7.126} T^{n+1}_j=(1-2D)\cdot T^n_j+D[(1-1/j)\cdot T^n_{j-1}+(1+1/j)\cdot T^n_{j+1}],\ j=1,2,\dots$$$$

As in Example 7.7.1 Example: Start-up flow in a tube, we employ two versions of the symmetry condition for $$r=0$$. 1)

From (7.97) with $$\lambda = 2$$: $$$$\tag{7.127} T^{n+1}_0=(1-6D)\cdot T^n_0+6D\cdot T^n_1$$$$

2) From (7.99): $$$$\tag{7.128} T^n_0=\frac{1}{3}(4T^n_1-T_2^n),\ \text{alle } n$$$$

Boundary conditions for $$r=b$$.

Discretizing (7.123) using 2nd order backward differences: $$\begin{equation*} k\cdot \left( \frac{3T^n_{N+1}-4T^n_N+T^n_{N-1}}{2\cdot\Delta r} \right) = \bar h\cdot (T_v-T_b) \end{equation*}$$

solving for $$T^n_{N+1}$$ yields: $$$$\tag{7.129} T^n_{N+1}=\frac{4T^n_N-T^n_{N-1}+2\delta\cdot T_v}{3+2\delta}$$$$ with $$$$\tag{7.130} \delta = \frac{\Delta r\cdot \bar h}{k}$$$$

We have previously shown that we must have $$D < 1/3$$ when we use boundary condition (7.127). In Example 7.7.1 Example: Start-up flow in a tube for a cylindrical domain, we found that the stability limit of $$D$$ increased when we reduced $$\Delta r$$ using the FTCS scheme. On the other hand, for the spherical case it appears that the condition $$D < 1/3$$ is independent of $$\Delta r$$. The reason for this can be understood by writing (7.126) for $$j=1$$: $$\begin{equation*} T^{n+1}_1=(1-2D)\cdot T^n_1+2D\cdot T^n_2 \end{equation*}$$

The term $$(1-1/j)\cdot T^n_{j-1}=(1-1)\cdot T^n_0$$ vanishes for $$j=1$$, so that the temperature in the center of the ball does not affect any other point. This explains why the stability limit is independent of $$\Delta r$$. We can actually solve for $$j=1,2,\dots,N$$ without bothering about boundary conditions (7.127) and (7.128) (see [14] for further details). For boundary condition (7.128) we only need to find the temperature $$T^n_0$$ at the center of the sphere.

von Neumann stability analysis without considering the boundary conditions showed that (7.126) is stable for $$D\leq 1/2$$. In addition, the analysis showed that the influence of the variable coefficient disappeared for both, the cylindrical and the spherical cases (see discussion in connection with (7.93)). We have not shown that $$D\leq 1/2$$ is an adequate condition for the entire system, since that requires to include the boundary condition (7.129).

We can summarize the usage of the FTCS scheme for the spherical case as follows:

The scheme in (7.129) is stable for $$D < 1/2$$ for $$j=1,2,\dots$$

If the temperature at the center of the sphere is also required we must ensure that $$D < 1/3$$ when (7.127) is used.

Using (7.128), the temperature at the center of the sphere can be computed for $$D < 1/2$$.

This is in agreement with Eisen's analysis, mentioned in connection to (7.99).

Below we show the printout resulting from running program kule for a computation where we used $$D=0.4$$ and $\Delta r=0.1$cm. $$T_k=300^\circ C$$, $$T_v=20^\circ C$$ and a simulation time of 10 minutes, with time steps of 1 second. The results for the analytical solution are computed using the program function kanalyt. We see that there is a good agreement between analytical and numerical values (the analytical solution is derived in Appendix G.9 of the Numeriske Beregninger.)

 r(cm) T ($$^\circ C$$) $$T_a$$ r(cm) T ($$^\circ C$$) $$T_a$$ 0.00 53.38 53.37 2.60 49.79 49.78 0.10 53.37 53.36 2.70 49.52 49.51 0.20 53.36 53.35 2.80 49.24 49.23 0.30 53.33 53.32 2.90 48.95 48.94 0.40 53.29 53.28 3.00 48.65 48.64 0.50 53.24 53.23 3.10 48.35 48.34 0.60 53.18 53.17 3.20 48.03 48.03 0.70 53.11 53.10 3.30 47.71 47.71 0.80 53.03 53.02 3.40 47.38 47.38 0.90 52.93 52.93 3.50 47.05 47.04 1.00 52.83 52.82 3.60 46.70 46.70 1.10 52.72 52.71 3.70 46.35 46.35 1.20 52.59 52.58 3.80 46.00 45.99 1.30 52.46 52.45 3.90 45.63 45.63 1.40 52.31 52.30 4.00 45.26 45.26 1.50 52.16 52.15 4.10 44.89 44.88 1.60 51.99 51.98 4.20 44.50 44.50 1.70 51.81 51.81 4.30 44.11 44.11 1.80 51.63 51.62 4.40 43.72 43.71 1.90 51.43 51.42 4.50 43.32 43.31 2.00 51.22 51.22 4.60 42.92 42.91 2.10 51.01 51.00 4.70 42.51 42.50 2.20 50.78 50.78 4.80 42.09 42.09 2.30 50.55 50.54 4.90 41.67 41.67 2.40 50.30 50.30 5.00 41.25 41.24 2.50 50.05 50.04

Animation for example on the cooling of a sphere using the FTCS scheme, as well as analytical solution.