Steady-state heat transfer analysis in a spherical domain revisited

The paper discusses the numerical solution of the onedimensional radially axi-symmetric non-linear second-order differential equation to model the conduction and radiation transfer through a spherical domain as a result of an exothermic heat source. The equation is transformed to a non-dimensional form. The dimensionless numbers emanating from the transformation represent the effect of the reaction rate, reaction type, activation energy, radiation and the convection on the temperature. The non-dimensional differential equation for the temperature distribution was previously solved using the Runge-Kutta-Fehlberg method coupled with a Shooting technique. In this paper the solution of the nondimensional differential equation using an iterative Galerkin finite element method approach employing the Picard method is described. The commercial finite element code Comsol is also employed to solve the nondimensional differential equation. The current study was motivated by inconsistencies that were observed in the previous results that were presented. Although the assumed underlying physics is used to evaluate the results, the study focuses purely on the numerical solution of the nondimensional differential equation. The results obtained by the Galerkin finite element code and Comsol were found to be in exact agreement and also exhibit no inconsistencies.


Introduction
Lebelo et al. [1] discussed a one-dimensional radially axi-symmetric non-linear secondorder differential equation to model the conduction and radiation transfer through a spherical domain as a result of an exothermic heat source. A zero heat flux was applied at the centre of the domain and at the boundary a convection heat flux. The equation was transformed to a non-dimensional form. The dimensionless numbers emanating from the transformation represent the effect of the reaction rate, reaction type, activation energy, radiation and the convection on the temperature. The non-dimensional differential equation for the temperature distribution was solved using the Runge-Kutta-Fehlberg method coupled with a Shooting technique. The results for solutions obtained for selected combinations of the dimensionless numbers were presented by Lebelo et al. [1].
In this paper the solution of the non-dimensional differential equation using an iterative Galerkin finite element method approach employing the Picard method is described [2][3]. The commercial finite element code Comsol [4] was also employed to solve the nondimensional differential equation. Solutions for the same selected combinations of the dimensionless numbers as used by Lebelo et al. [1] were obtained.
The current study was motivated by inconsistencies that were observed in the results presented by Lebelo et al. [1]. Although the assumed underlying physics is used to evaluate the results, no attempt was made to determine the validity of the mathematical formulation of the physical problem. The study focused purely on the numerical solution of the nondimensional differential equation.
The results obtained by the Galerkin finite element code and Comsol were found to be in exact agreement and also exhibit no inconsistencies. The paper will provide a full comparison and evaluation of the results presented by Lebelo et al. [1] and the results obtained by the finite element codes.

Governing equations
The conduction and radiation heat transfer through a spherical domain due to an exothermic heat source can be described by the one-dimensional radially axi-symmetric non-linear second-order differential equation [1]: Where k is the thermal conductivity, r the radial distance, T the temperature, Q the heat of reaction, A the rate constant,  is the Stefan-Boltzmann constant,  the vibration frequency, l the Planck number, E the activation energy, R universal gas constant,  a dimensionless activation energy parameter, and w T the ambient temperature. The parameter m represents the chemical reaction kinetics type. The equation is subject to the boundary conditions [1]: Where a is the radius of the sphere and h the heat transfer coefficient.  Where Bi is the Biot number and 0 1 r   . Lebelo et al. [1] solved eq. (3) numerically subject to the boundary conditions given in eq. (4) employing the Runge-Kutta-Fehlberg method coupled with the Shooting technique [5].

Finite element formulation
The weak formulation of eq. (3) is formed by multiplying the differential equation with 2 r and an arbitrary but suitable weighting function v and integrating over the domain [2]. This then gives: Integrating eq. (5) by parts, rearranging the result, expanding the boundary terms and substituting the boundary conditions then leads to : The equation can be solved using an iterative Galerkin finite element (GFE) approach employing the Picard method [2]. For the (i+1)-th iteration with the solution for the i-th iteration available we get : The last equation accounts for the boundary condition on the right-hand side. The solution for the (i+1)-th iteration is then given by : The solution is considered to be converged when : For all the simulations 25 quadratic elements were employed with the relaxation parameter taken as 0.5
The dimensionless non-linear second-order differential equation was also solved using the finite element method as implemented in the Comsol Multiphysics software (version 3.5a) [4]. The domain 0 1 r   was divided into 60 uniformly sized elements on which quadratic Lagrange interpolation was employed resulting in a total is 121 nodal points. Numerical integration was performed using Guassian quadrature with an integration order of four. The problem was solved iteratively by making use of a Damped Newton method with a relative error tolerance of 6 10  .

Results
Lebelo et al. [1] studied the effect on the dimensionless temperature distribution of selected values for  , m ,  , Ra and Bi . The following arbitrary reference values were selected for the parameters, namely 0.1 Lebelo et al. [1] stated that all the results are due to exothermic reactions that are taking place in the domain and where heat is thus generated inside the domain. In steady-state conditions conservation of energy requires that heat be released to the surrounding environment at the boundary of the domain. In order for the generated heat to be transferred to the boundary and released to the surrounding environment, the temperature gradient in the domain must be 0 d dr Five sets of simulations were performed by Lebelo et al. [1] and repeated in the current study. In each set the value of one of the parameters was varied whilst the values of the other parameters were kept constant.
In the first set of simulations the values of m ,  , Ra and Bi are kept constant at the reference values, whilst the values 0.1, 0.2, 0.3 and 0.4 were selected for  . According to Lebelo et al. [1] increasing the value of  represents an acceleration in the rate of the reaction that supports the exothermic chemical reaction.
The results obtained by Lebelo et al. [1] are depicted in Fig. 1, whilst the results obtained in the current study are depicted in Fig. 2.  It can be seen in both Figs. 1 and 2 that the gradients of the temperature profiles are negative as expected for an exothermic chemical reaction. Also in both Figs. 1 and 2 the level of the temperature increases as the rate of reaction increases and more heat is released. In Fig. 2 it can be seen that the GFE and Comsol results are in exact agreement. However, the temperature levels obtained by Lebelo et al. [1] are considerably higher than the corresponding temperature levels obtained by the GFE and Comsol simulations.  Fig. 4 shows the corresponding results that were obtained in the current study. Again it can be observed in Fig. 4 that the GFE and Comsol results are in exact agreement. However, the GFE and Comsol results differ markedly from the corresponding results shown in Fig. 3. The results obtained by Lebelo et al. [1] depicted in Fig. 3 are indicative of an endothermic chemical reaction where heat is absorbed from the surrounding environment. In contrast the GFE and Comsol results shown in Fig. 4 correspond to the assumed exothermic chemical reaction. It can be observed in Fig. 4 that the reaction kinetics type has an almost negligible effect on the temperature level. However, according to the results shown in Fig. 3 the reaction kinetics type has a marked effect on the temperature level.
In the third set of simulations the values of  , m , Ra and Bi are kept constant at the reference values, whilst the values 0.1, 0.2, 0.3 and 0.4 were selected for  . According to Lebelo et al. [1]  is the activation energy parameter and a reduction in the temperature level due to  implies a reduction in the activation energy to stimulate the exothermic chemical reaction.
The results obtained by Lebelo et al. [1] are shown in Fig. 5, whilst the results obtained from the GFE and Comsol simulations are shown in Fig. 6. It can again be observed in Fig.  6 that the GFE and Comsol results are in exact agreement.  [1] Ra deals with the rate at which heat is released from the surface of the sphere to the surrounding environment due to radiation. When Ra increases the rate at which heat is released also increases and results in the reduction of the temperature level.
The results obtained by Lebelo et al. [1] are shown in Fig. 7, whilst the results obtained from the GFE and Comsol simulations are shown in Fig. 8. It can again be observed in Fig.  8 that the GFE and Comsol results are in exact agreement.  It can be seen in both Figs. 7 and 8 that the temperature level decreases as the value of Ra increases. However, in Fig. 7 the results for 1.0 Ra  and 2.0 imply that the chemical reaction is exothermic, whilst the results for 3.0 Ra  and 4.0 imply that the chemical reaction is endothermic. In Fig. 8 [1] Bi deals with the rate at which heat is released for the surface of the sphere to the surrounding environment due to convection. When Bi increases the rate at which heat is released also increases and results in the reduction of the temperature level. The results obtained by Lebelo et al. [1] are shown in Fig. 9, whilst the results obtained from the GFE and Comsol simulations are shown in Fig. 10.
It can be seen in both Figs. 9 and 10 that the gradients of the temperature profiles are negative as expected for an exothermic chemical reaction. Also in both Fig. 9 and 10 the level of the temperature decreases as the rate of convection heat transfer to the surrounding environment increases. In Fig. 10 it can be seen that the GFE and Comsol results are in exact agreement. However, the temperature levels obtained by Lebelo et al. [1] are considerably higher than the corresponding temperature levels obtained by the GFE and Comsol simulations.
In an attempt to resolve the difference between the temperature levels of the solutions obtained by Lebelo et al. [1] and the GFE and Comsol solutions, additional GFE simulations were performed where the simulations were initialized by taking    

Conclusions
In this paper the numerical solution of the one-dimensional radially axi-symmetric nonlinear second-order differential equation to model the conduction and radiation transfer through a spherical domain as a result of an exothermic heat source presented by Lebelo et al. [1] was evaluated. The evaluation was performed by generating corresponding solutions using a Galerkin Finite Element approach and the Comsol Multiphysics code [4].
In all the cases for which simulations were performed the GFE and Comsol solutions were in exact agreement and consistent with the described underlying physics of the problem. In particular the results indicated that the chemical reaction is exothermic with heat being released to the surrounding environment.
The solutions presented by Lebelo et al. [1] contain a number of inconsistencies. The set of simulations demonstrating the effect of the reaction kinetics on the temperature profile lead to results that are inconsistent with the assumed exothermic chemical reactions and are also inconsistent with the boundary condition prescribed at the surface of the sphere. The sets of results associated with the effects of the activation energy and the radiation parameter on the temperature profile also contain inconsistencies. Half of the results in each set are indicative of exothermic chemical reactions, whilst the other half of the results are indicative of endothermic chemical reactions. In the latter cases the results are also inconsistent with the prescribed boundary condition at the surface of the sphere.
It is recommended that the implementation of the Runge-Kutta-Fehlberg method coupled with the Shooting technique [5] by Lebelo et al. [1] be revisited and clarified.