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.
Additional data are
- 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 must now solve the following problem with \( T=T(r,t) \): $$ \begin{equation} \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), \end{equation} $$
with boundary conditions: $$ \begin{equation} \tag{7.123} \frac{\partial T}{\partial r}(0,t)=0,\ \text{(symmetribetingelse)} \end{equation} $$
For \( r=b: \) $$ \begin{equation} \tag{7.124} k\frac{\partial T}{\partial r}= \bar h \cdot (T_v-T_b) \end{equation} $$
Initial conditions are: $$ \begin{equation} \tag{7.125} T(r,0)=T_k \end{equation} $$
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 \): $$ \begin{equation} \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 \end{equation} $$
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 \): $$ \begin{equation} \tag{7.127} T^{n+1}_0=(1-6D)\cdot T^n_0+6D\cdot T^n_1 \end{equation} $$
2) From (7.99): $$ \begin{equation} \tag{7.128} T^n_0=\frac{1}{3}(4T^n_1-T_2^n),\ \text{alle } n \end{equation} $$
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: $$ \begin{equation} \tag{7.129} T^n_{N+1}=\frac{4T^n_N-T^n_{N-1}+2\delta\cdot T_v}{3+2\delta} \end{equation} $$ with $$ \begin{equation} \tag{7.130} \delta = \frac{\Delta r\cdot \bar h}{k} \end{equation} $$
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.