Analytical Solution and Numerical Simulation of Heat Transfer in Cylindrical and Spherical Shaped Bodies

1. Introduction

In nature and technologies, transport processes are essential because they drive an important number of phenomena [1]. One of them is diffusion, where energy or particles are transported in a specific way [2,3]. The phenomenon of spreading in the universe occur on a large scale from atoms to stars [4,5].

In the case of the models of two-dimensional diffusion, an important step has been realized by Machta and Zwanzig in Lorentz gas [6]. This work was followed by further simulations in Lorentz gas regarding the connection of diffusivity to certain phase space parameters [7] or to entropy [8]. Diffusion on surfaces is also a relevant topic. In the case when the surface may cause a chaotic dynamic, the phenomena has been studied in Refs. [9,10,11].

Recently, exciting results have been obtained regarding diffusion with large spatial extension, which covers important analytical and computational results [12,13,14,15,16]. However, the studied geometry of the system is either one-dimensional or if two-dimensional, then the Cartesian coordinate system is used. In reality, the geometry of the studied system often drives the researcher to use other coordinate systems, most importantly spherical and cylindrical. A typical example is the case when a piece of ice is melting in water [17]. In nuclear physics, neutron diffusion occurs in finite spherical reactors [18]. Reactive diffusion through triple layers in spherical geometry was studied by Erdélyi and Schmitz [19] experimentally as well as numerically by the finite volume method. Cylindrical geometry was used to study diffusion in nanowires and nanorods by Roussel et al. [20]. The classic textbook [21] and [22] covers the theory of heat conduction in solids and that of diffusion and reaction in permeable catalysts, respectively, including spherically symmetric cases.

Heat conduction and similar problems are routinely solved by well-established numerical methods. However, in our view, these are far from being the optimal ones since they have serious drawbacks. Conventional explicit finite difference schemes are unstable when the applied time step size exceeds the so-called Courant–Friedrichs–Lewy (CFL) limit, which is usually rather low. This is the main reason that implicit algorithms, which possess much better stability properties, are often used to tackle these equations [23,24,25,26]. These methods involve the solution of an algebraic equation system with the whole system matrix, which can be very large, especially if the number of spatial dimensions is larger than one. This is often time- and computer memory consuming and cannot be straightforwardly parallelized.

In the last years, our research team developed some explicit methods that are stable for the heat conduction equation for arbitrary time step sizes in any number of spatial dimensions [27,28]. The algorithms were examined in several systems with homogeneous and inhomogeneous material parameters. They are orders of magnitude faster than the conventional solvers, including the standard Runge-Kutta schemes and the common MATLAB 'ode' solvers. However, all of these tests were performed in Cartesian coordinate systems, and now it is high time to perform them in cylindrical and spherical systems.

In the next section, we derive the studied equation and present its discretization. It is followed by the derivation of the analytical solution, which is valid in both cylindrical and spherical cases. Section 4 is devoted to the description of the numerical algorithms, while Section 5 is about the verification of the methods using the analytical solution. Then we turn our attention to examining the performance of the methods in reproducing experimental results: first, the circumstances are described, and then the numerical results of the Ansys software and our code are presented. Finally, the last section summarizes our conclusions.

2. The Studied Problem

We derive the heat transfer equation (conduction, convection, and radiation) in cylindrical and spherical coordinates based on energy balance. Firstly, in cylindrical coordinates, consider a small 3D cylindrical element$\Delta V=\Delta \varphi \hspace{0.17em}(r+\raisebox{1ex}{$\Delta r$}\!\left/ \!\raisebox{-1ex}{$2$}\right.)\Delta r\times \Delta z$. In the case of full cylindrical symmetry, it is better to choose a full ring-shaped element, which yields $\Delta V=2\pi \hspace{0.17em}(r+\raisebox{1ex}{$\Delta r$}\!\left/ \!\raisebox{-1ex}{$2$}\right.)\Delta r\times \Delta z=\pi \hspace{0.17em}\left({(r+\Delta r)}^2-r^2\right)\Delta z$ .

Figure 1. Visualization of a cylindrical element.

Figure 1. Visualization of a cylindrical element.

We assume that physical quantities, including the temperature, are not changing in the ϕ -direction. It means we deal only with two-dimensional problems from the mathematical point of view, as shown in Fig. 1. The energy balance in this element during a time interval can be expressed as:

$$\begin{array}{l}\left(\begin{array}{l}\mathrm{Rate}\mathrm{of}\mathrm{heat}\mathrm{conduction}\\ \mathrm{at}r,z\end{array}\right)-\left(\begin{array}{l}\mathrm{Rate}\mathrm{of}\mathrm{heat}\mathrm{conduction}\mathrm{at}\\ r+\Delta r,z+\Delta z\end{array}\right)+\left(\begin{array}{l}\mathrm{Rate}\mathrm{of}\mathrm{heat}\mathrm{generation}\mathrm{inside}\\ \mathrm{and}\mathrm{on}\mathrm{the}\mathrm{surface}\mathrm{of}\mathrm{the}\mathrm{element}\end{array}\right)\\ \pm \left(\begin{array}{l}\mathrm{Rate}\mathrm{of}\mathrm{convection}\\ \mathrm{at}\mathrm{the}r,\hspace{0.17em}\hspace{0.17em}z\mathrm{element}\end{array}\right)\pm \left(\begin{array}{l}\mathrm{Rate}\mathrm{of}\mathrm{radiation}\\ \mathrm{at}\mathrm{the}r,\hspace{0.17em}\hspace{0.17em}z\mathrm{element}\end{array}\right)=\left(\begin{array}{l}\mathrm{Rate}\mathrm{of}\mathrm{change}\mathrm{of}\mathrm{energy}\\ \mathrm{content}\mathrm{of}\mathrm{the}\mathrm{element}\end{array}\right)\end{array}$$

or briefly,

$$Q_r+Q_{\varphi }+Q_z-Q_{r+\Delta r}-Q_{\varphi +\Delta \varphi }-Q_{z+\Delta z}+Q_{gen}\pm Q_{convection}\pm Q_{radiation}=\frac{\Delta E_{element}}{\Delta t}$$

(1)

To fill Eq. (1) with concrete formulas, the following three well-known laws are used.

Fourier’s law of heat conduction:

$$Q_r=-kS\frac{\Delta u}{\Delta r},\hspace{0.17em}\hspace{0.17em}{Q}_z=-kS\frac{\Delta u}{\Delta z},$$

(2)

In the given equation,$u=u(\overrightarrow{r},t)$ represents the temperature,$k=k(\overrightarrow{r})$ represents the thermal conductivity of the material, and S represents the surface area through which the heat flows.

Newton's law of heat convection can be stated as follows:

$$Q_{convection}=hS\Delta u=hS(u_a-u),$$

(3)

in the given equation, the symbol h denotes the convection heat transfer coefficient. The ambient temperature $u_a$ is independent of the temperature u of the body under examination. Hence, the term $hS{u}_a$is included in the heat generation term to account for this effect.

Stefan-Boltzmann law governs the incoming and outgoing heat radiation, and its expression is as follows:

$$Q_{radiation}={\sigma }^{\ast }S\left(u_a^4-u^4\right),$$

(4)

In the given context${\sigma }^{\ast }=SB\cdot \epsilon$, the universal Stefan-Boltzmann constant $SB=5.67\times {10}^{-8}{\mathrm{W}/\mathrm{m}}^2{\mathrm{K}}^4$needs to be multiplied by $\epsilon =0.85$the emissivity constant to account for the surface not being a black body. The incoming radiation ${\sigma }^{\ast }S{u}_a^4$is also considered in the heat source term q as the $hS{u}_a$absorptivity term. Consequently, the change in energy of an element over a specific time interval $\Delta t$ can be expressed as.

$$\Delta E_{element}=E_{t+\Delta t}-E_t=mc(u_{t+\Delta t}-u_t)=\rho c\Delta V(u_{t+\Delta t}-u_t)$$

(5)

where$\rho =\rho (\overrightarrow{r})$,$c=c(\overrightarrow{r})$ are the density and the specific heat, respectively.

From these equations one can derive the heat-transport equation in a 3D cylindrical coordinate system, which can be written as:

$$\frac{1}{r}\frac{\partial }{\partial r}\left(k\hspace{0.17em}r\frac{\partial u}{\partial r}\right)+\frac{1}{r^2}\frac{\partial }{\partial \varphi }\left(k\hspace{0.17em}r\frac{\partial u}{\partial \varphi }\right)+\frac{\partial }{\partial z}\left(k\frac{\partial u}{\partial z}\right)+\frac{Q_{gen}}{\Delta V}-\frac{hSu}{\Delta V}-\frac{{\sigma }^{*}S{u}^4}{\Delta V}=\rho c\frac{\partial u}{\partial t}$$

(6)

In the case of spherical coordinates, a small 3D spherical element can be seen in Fig. 2. The heat-transport equation for this case can be expressed as follows:

$$\frac{1}{r^2}\frac{\partial }{\partial r}\left(k\hspace{0.17em}{r}^2\frac{\partial u}{\partial r}\right)+\frac{1}{r^2{\sin }^2\theta }\frac{\partial }{\partial \varphi }\left(k\hspace{0.17em}r\frac{\partial u}{\partial \varphi }\right)+\frac{1}{r^2\sin \theta }\frac{\partial }{\partial \theta }\left(k\hspace{0.17em}\sin \theta \frac{\partial u}{\partial \theta }\right)+\frac{Q_{gen}}{\Delta V}-\frac{hSu}{\Delta V}-\frac{{\sigma }^{*}S{u}^4}{\Delta V}=\rho c\frac{\partial u}{\partial t}$$

Figure 2. Visualization of a spherical element.

Figure 2. Visualization of a spherical element.

If one does not consider the convection, radiation and source terms in Eq. (6) and assume that the material properties are homogeneous, one obtains the form of the heat conduction equation in cylindrical and spherical coordinate systems. We study only symmetrical systems, which means no relevant physical quantities depend on coordinate $\varphi$ in the cylindrical and on coordinates $\varphi$ and $\theta$ in the spherical case. If we – temporarily – also assume that nothing depends on the z coordinate in the cylindrical case, only the radius r remains as a spatial variable, which yields

$$\frac{\partial u}{\partial t}=\alpha \frac{1}{r^n}\frac{\partial }{\partial r}\left(\hspace{0.17em}{r}^n\frac{\partial u}{\partial r}\right),$$

(7)

where n=0, 1, and 2, which mean Cartesian, cylindrical and spherical coordinates respectively, while $\alpha =\frac{k}{c\rho }$ is the diffusivity. Eq. (7) is also used for particle-diffusion, where the diffusivity is usually denoted by D.

2.1. The spatial discretization of the problem

Instead of directly discretizing PDEs (6) or (7), we use an equivalent resistance-capacitance model. In the case of cylindrical geometry, we consider tube-shaped cells with height $\Delta z$ and thickness $\Delta r$. For spheres the cells have spherical-shell shapes with thickness $\Delta r$ again. The temperature is considered at the middle of the cell-layer, where the radial distance from the origin (the mean radius of the cells) is denoted by $r_i$, while the subsequent radius of the cell border is denoted by $r_i^{\ast }=r_i+\Delta r/2$.

The cells’ heat capacity in the cylindrical and in the spherical case is approximated as $C_i=c_i{\rho }_i\pi \left(r_{i+1}^{\ast 2}-r_i^{\ast 2}\right)\Delta z$, and $C_i=c_i{\rho }_i\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\pi \left(r_{i+1}^{\ast 3}-r_i^{\ast 3}\right)$, respectively.

Let us denote the area of the cylindrical cell-surface perpendicular to r with $S_r$, which can be given as $S_r=2\pi \hspace{0.17em}r\hspace{0.17em}\Delta z$. Now the thermal resistance in r-direction, the approximate formula

$$R_{i,i+1}\approx {\displaystyle {\int }_{r_i}^{r_{i+1}}\frac{dr}{k_{i,i+1}{S}_r}}={\displaystyle {\int }_{r_i}^{r_{i+1}}\frac{dr}{k_{i,i+1}2\pi r\Delta z}}=\frac{\ln (r_{i+1}-r_i)}{2\pi k_{i,i+1}\Delta z}$$

(8)

is used. For the thermal resistance in the z-direction, the approximate formula $R_{i,i+N_x}\approx \frac{\left(z_{i+N_r}-z_i\right)}{k_i\pi (r_{i+1}^2-r_i^2)}$ is used, where the cell i+N_r is below the cell i.

In the spherical case, $S_r$ can be given as $S_r=4\pi \hspace{0.17em}{r}^2$. Using this, the thermal resistance is calculated similarly as in the cylindrical case, but now the integration yields $R_{i,i+1}\approx \frac{1}{4\pi k_{i,i+1}}\frac{r_{i+1}-r_i}{r_i\hspace{0.17em}{r}_{i+1}}$. From Eq. (2) and (5) it is easy to obtain the ODE system

$$\frac{d{u}_i}{dt}={\displaystyle \sum _{j\ne i}\frac{u_j-u_i}{R_{i,j}{C}_i}}+\frac{Q_{gen}}{C_i}-\frac{hS{u}_i}{C_i}-\frac{{\sigma }^{*}S{u}_i^4}{C_i}$$

(9)

to determine the time-evolution of the cell temperatures. Here S is the area of the surface on which the convection and radiation occurs, which will be the outer surface of the cylinder in Section 7. In publications [29,30] the interested readers can find more details about this treatment of heat conduction. If one neglects the higher powers of $\Delta r$, one can easily derive that $C_i/S=c_i{\rho }_i\Delta r$ in both cases. Let us use the following notations:

$$K=\frac{h}{c\rho \Delta r}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\sigma =\frac{{\sigma }^{\ast }}{c\rho \Delta r}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}q=\frac{{\sigma }^{\ast }}{c\rho \Delta r}{u}_a^4+\frac{h}{c\rho \Delta r}\cdot u_{\mathrm{a}}$$

Inserting these into (9) we can write Eq. (9) into a simpler form:

$$\frac{d{u}_i}{dt}={\displaystyle \sum _{j\ne i}\frac{u_j-u_i}{R_{i,j}{C}_i}}+q_i-K{u}_i-\sigma u_i^4,$$

(10)

which will be solved numerically in Section 5 and Section 7.

3. The analytical solution

In this section we outline how the radial solutions of Eq. (7) look like for both cylindrical and spherical coordinate systems. The self-similar Ansatz we use:

$$u\left(r,\hspace{0.17em}t\right)=t^{-a}f\left(r/\sqrt{t}\right)=t^{-a}f\left(\eta \right)$$

Substituting the first and second derivative of the Ansatz into the original Eq. (2), we arrive at an ordinary differential equation (ODE) for$f\left(\eta \right)$

$$-af-\frac{\eta }{2}{f}^{\prime }=\alpha \left(\frac{n\hspace{0.17em}{f}^{\prime }}{\eta }+f^{″}\right),$$

(11)

with the constraint that a is an arbitrary real number. The software Maple12 gives the following solution to this ODE:

$$f\left(\eta \right)\hspace{0.17em}\hspace{0.17em}=\hspace{0.17em}\hspace{0.17em}{\mathrm{e}}^{\frac{-{\eta }^2}{4\alpha }}\left[c_{1}M\left(\frac{1}{2}+\frac{n}{2}+a,\hspace{0.17em}\hspace{0.17em}\frac{1}{2}+\frac{n}{2},\hspace{0.17em}\hspace{0.17em}\frac{{\eta }^2}{4\alpha }\right)\hspace{0.17em}\hspace{0.17em}+\hspace{0.17em}{c}_{2}U\left(\frac{1}{2}+\frac{n}{2}+a,\hspace{0.17em}\hspace{0.17em}\frac{1}{2}+\frac{n}{2},\hspace{0.17em}\hspace{0.17em}\frac{{\eta }^2}{4\alpha }\right)\hspace{0.17em}\right]\hspace{0.17em},$$

(12)

Unfortunately, for unclear reasons substituting the n = 0 we cannot get back the original Cartesian solution, some factor is missing. Luckily, Eq. (12) describes well the cylindrical and spherical symmetric solutions. Note that unlike the Cartesian case, here all the solutions have the even property. In our former theoretical publications [12,31] we analyzed the possible different kind of solutions. Depending on the numerical values of parameter a, four species of solution curves exist. It is generally true for Cartesian, cylindrical and spherical systems as well, the only difference is just the numerical values of a.

So for large negative a values, the solutions show an exploding behavior at infinite temporal and spatial values, which obviously violates energy conservation or thermodynamics. When a is close to zero, we get solutions which have a finite non-zero numerical value at large times and special points. For small positive a values, the solution has a maximum at zero and has a similar decay to zero as the fundamental Gaussian solution. Finally, for large positive a values, the solutions have a power law explosion in the origin (at zero time and zero spatial coordinate) with a drastic decay accompanied by additional oscillations. So, the solutions become negative at some $\eta$ values. These are the general new features of these kinds of solutions for Cartesian, cylindrical and spherical systems as well.

Figure 3. The Kummer's M shape functions f(η) from Eq. 9 for the cylindrical (left figure) and spherical (right) symmetry for three different self-similar exponents a. The black, red and blue lines are for a = 1/2, 1 and 3/2. The diffusion coefficient α was set to unity for all cases.

Figure 3. The Kummer's M shape functions f(η) from Eq. 9 for the cylindrical (left figure) and spherical (right) symmetry for three different self-similar exponents a. The black, red and blue lines are for a = 1/2, 1 and 3/2. The diffusion coefficient α was set to unity for all cases.

For our forthcoming numerical analysis, instead of the solution (12) of the ODE (11), we need the final form of the solution of the PDE (7), which reads as follows:

$$u^{\mathrm{exact}}(r,\hspace{0.17em}\hspace{0.17em}t)=t^{-a}{e}^{-\frac{r^2}{4\alpha t}}\left(c_{1\hspace{0.17em}}M\left(\frac{1}{2}+\frac{n}{2}+a,\hspace{0.17em}\hspace{0.17em}\frac{1}{2}+\frac{n}{2},\hspace{0.17em}\hspace{0.17em}\frac{r^2}{4\alpha t}\right)+c_{2\hspace{0.17em}}U\left(\frac{1}{2}+\frac{n}{2}+a,\hspace{0.17em}\hspace{0.17em}\frac{1}{2}+\frac{n}{2},\hspace{0.17em}\hspace{0.17em}\frac{r^2}{4\alpha t}\right)\right),$$

(13)

Figure 4 presents the $u\left(r,\hspace{0.17em}t\right)$ solution for Kummer's M function for cylindrical symmetry. One can observe certain similarities with the regular Gaussian.

4. The applied numerical methods

We present the numerical schemes adapted for Equation (10). We always use equidistant discretization of the time variable with time step size $\Delta t$, and the temperature of cell i at the time $n\Delta t$ is denoted by $u_i^n$. The following quantities will be extensively used:

$$w_i=\Delta t{\displaystyle \sum _{j\ne i}\frac{1}{C_i{R}_{ij}}} \mathrm{and}\hspace{0.17em}{A}_i=\Delta t{\displaystyle \sum _{j\ne i}\frac{u_j^n}{C_i{R}_{ij}}}+\Delta t\cdot q_i.$$

The quantity w is the generalization of the so-called mesh-ratio (which has the expression$\frac{\alpha \Delta t}{\Delta x^2}$), but it now depends on the cell index. The aggregated quantity A, on the other hand, responsible for the conductive heat exchange with the neighbors of cell i as well as the source term.

1. The UPFD (unconditionally positive finite difference) method was proposed in [32] for the diffusion-advection-reaction equation a decade ago. We recently adapted it [28] to Eq. (10) as follows:

$$u_i^{n+1}=\frac{u_i^n+A_i}{1+w_i+\Delta t\cdot K_i+\Delta t\cdot {\sigma }_i\cdot {(u_i^{\mathrm{n}})}^3}.$$

2. The pseudo-implicit (PI) method [28] has the following two stages:

$$\mathrm{Stage}1:u_i^{\mathrm{pred}}=\frac{u_i^n+A_i/2}{1+w_i+\Delta t\cdot K_i/2+\Delta t\cdot {\sigma }_i\cdot {(u_i^{\mathrm{n}})}^3/2},$$

$$\mathrm{Stage}2:u_i^{n+1}=\frac{\left(1-w_i\right)u_i^n+A_i^{\mathrm{new}}+\Delta t\cdot K_i\left(u_i^{\mathrm{pred}}-u_i^n\right)}{1+w_i+\Delta t\cdot K_i+\Delta t\cdot {\sigma }_i\cdot {(u_i^{\mathrm{n}})}^3},\mathrm{where}{A}_i^{\mathrm{new}}=\Delta t{\displaystyle \sum _{\mathrm{j}\ne \mathrm{i}}\frac{u_{\mathrm{j}}^{\mathrm{pred}}}{C_i{R}_{\mathrm{ij}}}}+\Delta t\hspace{0.17em}{q}_i.$$

3. A two-stage version of the Rational Runge-Kutta methods [33] will be applied as follows. First, a full step is taken by the standard FTCS (explicit Euler) scheme to calculate the predictor values:

$$u_i^{\mathrm{pred}}=u_i^n+g_i^1,$$

where the increment $g_i^1$ is calculated as

$$g_i^1=A_i-\Delta t\cdot K_i\cdot u_i^{\mathrm{n}}-\Delta t\cdot {\sigma }_i\cdot {(u_i^{\mathrm{n}})}^4.$$

Using these predictor values, the increment in the second stage is

$$g_i^2=A_i^{\mathrm{new}}-\Delta t\cdot K\cdot u_i^{\mathrm{pred}}-\Delta t\cdot \sigma \cdot {(u_i^{\mathrm{pred}})}^4.$$

Now the new values of the variable are as follows

$$u_i^{n+1}=u_i^n+\frac{2p_1{g}_i^1-2p_{12}{g}_i^1+p_1{g}_i^2}{4p_1-4p_{12}+p_2}.$$

where the following scalar products are used:

$$p_1=\left({\overrightarrow{g}}^1,{\overrightarrow{g}}^1\right)={\displaystyle \sum _{i=1}^N{g}_i^1{g}_i^1},\hspace{0.17em}\hspace{0.17em}{p}_{12}=\left({\overrightarrow{g}}^1,{\overrightarrow{g}}^2\right)={\displaystyle \sum _{i=1}^N{g}_i^1{g}_i^2},\hspace{0.17em}\hspace{0.17em}{p}_2=\left({\overrightarrow{g}}^2,{\overrightarrow{g}}^2\right)={\displaystyle \sum _{i=1}^N{g}_i^2{g}_i^2}.$$

A special, checkerboard-like spatial grid must be constructed if someone wants to use any version of the odd-even hopscotch methods [34]. The cells of this grid are labeled as odd and even, with the requirement that all the nearest neighbors of the even cells are odd and vice versa. In the case of the original version (denoted by OOEH here), the odd-even labels must be interchanged after each time step as it is displayed in Figure 5A. The standard FTCS formula is modified in the first stage to make it slightly more stable [30] by treating the convection and radiation terms in an “implicit” way:

$$u_i^{n+1}=\frac{\left(1-w_i\right)u_i^n+A_i}{1+\Delta t\cdot K_i+\Delta t\cdot {\sigma }_i\cdot {(u_i^{\mathrm{n}})}^3}$$

The implicit (Euler) scheme in the second stage is as follows:

$$u_i^{n+1}=\frac{u_i^n+A_i^{\mathrm{new}}}{1+w_i+\Delta t\cdot K_i+\Delta t\cdot {\sigma }_i\cdot {(u_i^{\mathrm{n}})}^3}$$

It is important that in the second stage, similarly to the other odd-even hopscotch scheme below, always the latest available values of the neighbors are used when the new values of u are calculated. Due to this, the $u_j^{n+1}$ values in the implicit formula have been just obtained at Stage 1, which makes it effectively explicit.

5. In the case of the shifted-hopscotch (SH) algorithm [29], five stages constitute a two-step long repeating block. Two of the stages are half- and the remaining three of them are full-length steps, see in Figure 5 (C). The first stage is for the odd cells only and symbolized by a dark green box with the number 1 in the figure. It uses the formula

(14)

Then, a full-time step with formula

$$u_i^{\mu +1}=\frac{\left(1-w_i/2\right)u_i^{\mu }+A_i^{\mu +\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{1+w_i/2+\Delta t\cdot K_i+\Delta t\cdot {\sigma }_i\cdot {\left(u_i^{\mu }\right)}^3}.$$

(15)

for the even, the odd, and the even cells come again (blue rectangles in the figure). The superscript µ equals n for the even and n+1 for the odd cells. Finally, a half-length time step (pink box with number 5 inside) for the odd cells closes the calculation with the formula

$$u_i^{n+2}=\frac{\left(1-w_i\right)u_i^{n+1}+A_i^{n+\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{1+\Delta t\cdot K+\Delta t\cdot \sigma \cdot {\left(u_i^{n+1}\right)}^3}.$$

(16)

6. The asymmetric hopscotch (ASH) algorithm is almost the same as the SH one, but the repeating block is only one-time step long as one can see in Fig 5D. It contains only three stages instead of five which use the same formulas (14), (15) and (16), respectively.

7. The procedure of the leapfrog-hopscotch (LH) algorithm [27] begins and ends by a half-length stage, but then all other stages have full time step length (blue boxes in Fig. 5B). At the first and intermediate stages, it uses formulas (14) and (15). However, at the last stage (pink rectangle) it uses (15) again, but $\Delta t$ must be divided by 2, including their appearances in the quantities w and A.

8. Although the Dufort–Frankel (DF) method [35] is a known explicit and unconditionally stable algorithm for the linear heat equation, it is rarely used. The formula adapted for the case of Eq. (10) is:

$$u_i^{\mathrm{n}+1}=\frac{\left(1-w_i\right)u_i^{n-1}+2A_i}{1+w_i+2\cdot \Delta t\cdot K_i-2\cdot \Delta t\cdot {\sigma }_i\cdot {(u_i^{\mathrm{n}})}^3}.$$

It is a two-step scheme which is not self-starter, thus we apply the UPFD formula of point 1 to start the DF method.

9. One of the most common algorithms to solve the heat conduction equation is the explicit-Euler based FTCS (forward time central space) algorithm. It can be adapted to our case in the standard way as follows:

$$u_i^{n+1}=\left(1-w_i\right)u_i^n+A_i-\Delta t\cdot K_i\cdot u_i^n-\Delta t\cdot {\sigma }_i\cdot {(u_i^n)}^4.$$

5. Verification using the analytical solution

In this section, we take the height of the cylinder as well as $\Delta z$ unity. It means that computationally there is one space dimension only in both the cylindrical and the spherical case. The solution parameters are:

$$\begin{array}{l}{N}_r=500\hspace{0.17em},\hspace{0.17em}Nz=1,\hspace{0.17em}\hspace{0.17em}N=Nr\times Nz=500,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{r}_0=0.0003,\hspace{0.17em}\hspace{0.17em}{r}_{\max }=0.999,\hspace{0.17em}\hspace{0.17em}\Delta r=0.002\hspace{0.17em},\hspace{0.17em}\\ \alpha =1\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}a\in \left\{1,\hspace{0.17em}1.2,\hspace{0.17em}\hspace{0.17em}2\right\}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{c}_1=1\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}{c}_2=0,\hspace{0.17em}\hspace{0.17em}{t}^0=0.1\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{t}^{\mathrm{fin}}=t^0+\mathrm{0.1.}\hspace{0.17em}\end{array}$$

Here N represents the total number of cells, while $\hspace{0.17em}{r}_0$ and $r_{\max }$ are the radial coordinates of the center of the first and last cells. The initial condition obtained by substituting the initial t and boundary r values into the analytical solution, respectively. The Dirichlet boundary conditions on the right side (the circumference of cylinder and sphere) obtained simply by substituting the radius $r_{\max }$ into the analytical solution and calculate the function value at each time step. On the left side (the cylinder and sphere center, $r=r_0$), zero-Neumann boundary applied, since no heat can disappear from the center of the cylinder or the sphere. This boundary actually exists only computationally but not physically. We remind the reader that the analytical solutions are constructed for Eq. (7), and here we accurately reproduce them by solving the apparently different Eq. (10) numerically.

We calculate the so-called maximum error:

$$MaxError=\underset{1\le i\le N}{\max }\left|u_{\mathrm{i}}^{\mathrm{analytic}}(t_{\mathrm{fin}})-u_{\mathrm{i}}^{\mathrm{num}}(t_{\mathrm{fin}})\right|,$$

(20)

The obtained maximum errors are displayed as a function of the time step size in Fig. 6 and 7 for two values of parameter a in cylindrical coordinates. Figure 8 presents the temperature value as a function of r. For the case of spherical coordinates, Figure 9 shows the maximum error as a function of the time step, and Figure 10 presents the temperature as a function of r. The fact that we obtain very small errors in all cases verifies not only the numerical algorithms, but the equivalence of the two mathematical treatments of the physical problem.

Figure 6. The maximum errors as a function of the time step size $\Delta t$ for the 9 numerical methods in case of cylindrical coordinates for a=1.

Figure 6. The maximum errors as a function of the time step size $\Delta t$ for the 9 numerical methods in case of cylindrical coordinates for a=1.

Figure 7. The maximum errors as a function of the time step size for the 9 numerical methods in case of cylindrical coordinate for a=2.

Figure 7. The maximum errors as a function of the time step size for the 9 numerical methods in case of cylindrical coordinate for a=2.

Figure 8. The values of temperature as a function of variable r in case of the initial function u⁰, the analytical solution Uexact, the DF method and the LH method in case of cylindrical coordinate for a=1. .

Figure 8. The values of temperature as a function of variable r in case of the initial function u⁰, the analytical solution Uexact, the DF method and the LH method in case of cylindrical coordinate for a=1. .

Figure 9. The maximum errors as a function of the time step size for the 9 numerical methods in case of spherical coordinates for a=1.2.

Figure 9. The maximum errors as a function of the time step size for the 9 numerical methods in case of spherical coordinates for a=1.2.

Figure 10. The values of temperature as a function of r variable in case of the initial function u⁰, the analytical solution Uexact, the DF method and the LH method in case of spherical coordinate for a=1.2.

Figure 10. The values of temperature as a function of r variable in case of the initial function u⁰, the analytical solution Uexact, the DF method and the LH method in case of spherical coordinate for a=1.2.

6. Setup of the reproduction of the experimental results

6.1. Material properties

We are going to reproduce the experimental results of Cabezas et al. [36], where heat transfer is studied in a steel C45 cylinder with properties shown in Table 1 below.

The volume of the cylinder segment is$V=1.4362\times {10}^{-4}\hspace{0.17em}{m}^3$,while $\left(r,z\right)\in \left[0,\hspace{0.17em}0.0165\hspace{0.17em}\mathrm{m}\right]\times \left[0,\hspace{0.17em}0.168\hspace{0.17em}\mathrm{m}\right]$. In our approximation, physical quantities are not changing in the ϕ-direction, thus that dimension is irrelevant and computationally we deal with a two-dimensional problem. The total area of the meshes is 0.0174 m² in the r direction, and $8.5487\times {10}^{-4}{m}^2$ in the z direction. The number of the cells along the r axis and z axis are set to N_r = 15 and N_z = 100, thus the total number of the cells in the system is $N=N_r{N}_z=1500$. We emphasize that we use the same code as in Section 5, only some parameters such as N_z are different.

6.2. The initial and the boundary conditions

We use a constant initial condition in all cases.

$$u\left(r,z,\hspace{0.17em}t=0\right)=30.7\hspace{0.17em}{C}^{\circ }$$

As in Section 4, we use different boundary conditions on different sides. On the left side, the center of the cylinder, we applied Neumann boundary condition in all cases, which don’t allow conductive heat transfer at the boundary,

$$u_r(r=0,z,t)=u_r(r=L_r,z,t)=u_z(r,z=L_z,t)=0$$

On the right and upper boundaries, we use two types of boundary conditions. The first one is zero-Neumann, when there is no heat exchange with the environment. The second one when there is a heat exchange with the environment by convection and radiation, considering the heat convection coefficient $h=3\left(\mathrm{W}\cdot {\mathrm{m}}^{-2}\cdot {\mathrm{K}}^{-1}\right)$, and the emissivity constant as 0.85 to obtain realistic values for ${\sigma }^{\ast }$. On the lower boundary, we applied the constant heat flux with P=1500 W for all cases to follow the experimental setup in the paper [36].

The heat generation contains a fraction of the ambient temperature radiation. The ambient temperature of the air is taken with three values depending on the experimental measurments (30.7,31.1, and 31.7 C). The term q contains also the convective heat gain due to the nonzero temperature $u_a$ of the air (in Kelvin). The convective and radiative energy transfer is perpendicular to the surface. The interior elements cannot gain or lose heat by the heat source, heat convection, or radiation.

7. Simulation results

In this section, the end of the examined time interval is defined as$t_{\mathrm{fin}}=1200,1440\hspace{0.17em}\mathrm{and}\hspace{0.17em}1800\hspace{0.17em}\mathrm{s}$. To ensure that the errors of the numerical algorithms are small, we used a reference solution of the spatially discretized Eq. (9) obtained by the ode15s solver in MATLAB with a small tolerance${10}^{-10}$. This reference solution is utilized in Eq. (20) to determine the maximum error of the time integration, which does not contain the error of the mathematical model and that of the space discretization.

7.1. Results of the numerical methods

For the simulation, we have chosen the top 5 algorithms, namely DF, OOEH, LH, SH, and ASH. Figure 11present the maximum error results obtained from these selected algorithms.

The simulation of a steel C45 cylinder was conducted using the selected algorithms, considering different boundary conditions as previously mentioned. Among these algorithms, the shifted-hopscotch method was chosen to visualize the temperature contour due to its high accuracy at small time step size. Figs 12-13 display the final temperature distribution obtained from this method.

Figure 12. The final temperature distribution contour for different time values (t=20, 24, and 30 min, respectively from left to right) presented by the SH method when there is no heat exchange with the environment.

Figure 12. The final temperature distribution contour for different time values (t=20, 24, and 30 min, respectively from left to right) presented by the SH method when there is no heat exchange with the environment.

Figure 13. The temperature distribution contour for different time values (t=20, 24, and 30 min) presented by the SH method when there is heat exchange with the environment by convection and radiation.

Figure 13. The temperature distribution contour for different time values (t=20, 24, and 30 min) presented by the SH method when there is heat exchange with the environment by convection and radiation.

7.2. Ansys Simulation Results

Ansys workbench steady state thermal analysis, with Mechanical APDL solver was used to simulate the steel C45 cylinder. The mesh size was $1\times {10}^{-3}$, the total number of elements was 197183, since it was a computationally 3D problem. In Figs. 14-16 we present the final temperature distribution at three measurement times.

Figure 14. The temperature contour at time (t=20 min) presented by Ansys workbench 19.2 when there is no heat exchange with the environment (left side of the figure) and when there is a heat exchange (right side of the figure).

Figure 14. The temperature contour at time (t=20 min) presented by Ansys workbench 19.2 when there is no heat exchange with the environment (left side of the figure) and when there is a heat exchange (right side of the figure).

Figure 15. The temperature contour at time (t=24 min) presented by Ansys workbench 19.2 when there is no heat exchange with the environment (left side of the figure) and when there is a heat exchange (right side of the figure).

Figure 15. The temperature contour at time (t=24 min) presented by Ansys workbench 19.2 when there is no heat exchange with the environment (left side of the figure) and when there is a heat exchange (right side of the figure).

Figure 16. The temperature contour at time (t=30 min) presented by Ansys workbench 19.2 when there is no heat exchange with the environment (left side of the figure) and when there is a heat exchange (right side of the figure).

Figure 16. The temperature contour at time (t=30 min) presented by Ansys workbench 19.2 when there is no heat exchange with the environment (left side of the figure) and when there is a heat exchange (right side of the figure).

7.3. Comparsion of the results

The results of the experimental measurements, the finite element method (FEM) using Ansys Workbench, and the explicit numerical methods (exemplified by the shifted hopscotch method) were compared. Both FEM and SH were subjected to two types of tests, one considering convection and radiation effects, and the other excluding them. In Table 1 and 2 the comparison was conducted at two specific spatial points, and the results were measured at three different time moments. The temperatures are compared at four space points via plots in Figure 17, Figure 18 and Figure 19.

Figure 17. The temperature at the selected 4 measurement points in z at time t=20min.

Figure 17. The temperature at the selected 4 measurement points in z at time t=20min.

Figure 18. The temperature at the selected 4 measurement points in z at time t=24min.

Figure 18. The temperature at the selected 4 measurement points in z at time t=24min.

Figure 19. The temperature at the selected 4 measurement points in z at time t=30min.

Figure 19. The temperature at the selected 4 measurement points in z at time t=30min.

8. Discussion and Summary

This work was devoted to solving heat transfer problems in cylindrical and spherical geometries. Using the self-similar Ansatz, novel analytical solutions of the heat-conduction PDE were constructed, which contain the Kummer functions. Nine numerical algorithms were presented, most of which are recently introduced unconditionally stable explicit methods. To perform the verification, the analytical solutions are reproduced by these methods.

After these, experimental work is considered from the literature where a cylinder is heated from below, and the results are attempted to be reproduced by the Ansys commercial software, but without considering convection and radiation on the surface of the cylinder. In contrast to that, we reproduced the experimental results by considering convection and radiation as well by not only the Ansys but the explicit methods as well. Since in reality, convection and radiation are present, taking them into account makes the results closer to the experimental ones. Moreover, the explicit and stable schemes were more accurate and effective than the finite element software in all cases.