Abstract
A time dependent nonlinear partial differential equation modelling heat transfer in a porous radial fin is considered. The Differential Transformation Method is employed in order to account for the steady state case. These solutions are then used as a means of assessing the validity of the numerical solutions obtained via the Crank-Nicolson finite difference method. In order to engage in the stability of this scheme we conduct a stability and dynamical systems analysis. These provide us with an assessment of the impact of the nonlinear sink terms on the stability of the numerical scheme employed and on the dynamics of the solutions.
1. Introduction
Circular annular fins are extensively used to increase the rate of heat transfer from a heat source for a given temperature difference in heat exchange devices or to reduce the temperature difference between the heat source and the given heat flow rate of heat sink [1]. The use of fins is widely ranging in engineering applications where it is necessary to enhance the heat transfer from a surface to an adjacent coolant so that the fin can perform within the acceptable temperature limits. These engineering applications range from considerably large systems such as industrial heat exchangers to smaller systems such as transistors. Radial fins have been conventionally used as a coolant for internal combustion engines, heat exchanges, compressors, and so forth. Due to these extensive practical applications this is considered a vibrant field of research; this is true even more so since the problems that arise are nonlinear and hence not always solvable via analytical techniques.
Many research articles have investigated the use of porous fins [2]. Even though a porous material fin has low thermal conductivity, a vast area of the material comes into contact with the cooling agent enabling the porous fin to give superior performance [2]. Over the past decades, numerous studies have been conducted on the performance of annular fins [3–7]. Aziz and Rahman [8] examined a fin comprising functionally graded material and analysed the performance on the radial fin with a continuously increasing thermal conductivity in the radial direction. They discovered that the heat transfer as well as the fin efficiency and effectiveness are at their highest maximum values when the thermal conductivity of the fin varies inversely with the square of the radius. Furthermore, they found that the use of a spatially averaged thermal conductivity model is not recommended due to large errors occurring in some cases. Kiwan [9] conducted a thermal analysis of natural convection porous fins by introducing Darcy’s model to construct the energy equation governing the distribution of temperature. He further discovered that, by choosing a precise value for the thermal conductivity ratio and the fin length to thickness ratio, the performance of the porous fin exceeded the performance of the solid fin. A study was conducted by Kiwan and Zeitoun [10] to test the performance of rectangular porous fins mounted around the inner cylinder of a cylindrical annulus by performing a finite volume type numerical study. It was concluded that, in comparison to solid fins, porous fins provided higher transfer rates for similar configurations and that the heat transfer rate from the cylinder equipped with porous fins decreased as the fin inclination increased. Gorla and Bakier [11] investigated natural convection and radiation in porous fins. They found that the radiation transfers more heat in comparison to a similar model without radiation. Abu-Hijleh [12] analysed the effects of using permeable fins on the forced convection heat transfer from a horizontal cylinder. The results obtained were similar to results obtained as per Kiwan and Zeitoun [10] in terms of the permeable fins providing much higher heat transfer rates. To study radial fins, a combination of the Taylor transformation and finite difference approximation was implemented by Yu and Chen [13, 14]. They further performed a study on the optimization of a circular fin with variable thermal parameters. Naidu et al. [15] set forth a numerical study of natural convection from a cylindrical fin placed in a cylindrical porous enclosure. Hence, they conducted a conjugate conduction-convection analysis by solving the heat conduction equation. Moitsheki and Harley [16] studied the transient heat transfer through a longitudinal fin of various profiles by employing classical Lie point symmetry methods. They observed that for long periods of time the temperature profile becomes unusual for the heat transfer in longitudinal triangular and concave parabolic fins. This, however, was corrected by increasing the thermogeometric fin parameter. In recent studies, Darvishi et al. [17] accounted for the effects of radiation and convection heat transfer in a rectangular radial porous fin. This allowed for the heat flow to infiltrate the porous fin enabling a solid-fluid interaction to occur. They concluded that in a model containing radiation more heat is present than in a similar model without radiation. In a similar context, our model takes into account the time rate of change of internal energy and the heat flow due to conduction as well as the heat due to radiation and convection.
Extensive analytical studies have been done via the Differential Transformation Method (DTM) for the solution of problems such as the one under discussion. The DTM, which was first proposed by Zhou [18], is a seminumerical-analytical method applied to linear and nonlinear systems of ordinary differential equations. The method captures the exact solution in terms of a Taylor series expansion. This method has been successfully implemented in engineering applications [19–25]. Ndlovu and Moitsheki [26] derived approximate analytical solutions for the temperature distribution in a longitudinal rectangular and convex parabolic fin with temperature dependent thermal conductivity and heat transfer coefficients. These authors, for the first time, used a two-dimensional DTM for the transient heat conduction problem. Ertürk [27] constructed seminumerical-analytical solutions for a linear sixth-order boundary value problem using the DTM. It was observed that the method served as an effective and reliable tool for such problems. Recently, Torabi et al. [28] analysed the radiative radial fin with temperature dependent thermal conductivity by implementing the DTM as well as the Boubaker polynomials expansion scheme (BPES). Similar to the results obtained by Ertürk [27], suitable results were obtained in predicting the solution for both BPES and DTM. A study of a radial fin in terms of the fin’s thickness with convection heating at the base and the convection-radiative cooling at the tip was conducted by Aziz et al. [29]. Furthermore, they conducted an analysis using DTM and verified the results by comparing it to an exact analytical solution. The preceding literature clearly shows that the Differential Transformation Method has been applied to problems relating to many different fins, but no attempt has been made to apply it when investigating the heat transfer in a porous radial fin.
In terms of numerical investigations, an efficient, accurate, extensively validated, and unconditionally stable method was developed in the mid-20th century by Crank and Nicolson [30] in order to evaluate numerical solutions for nonlinear partial differential equations. Rani et al. [31] obtained a solution for the time dependent nonlinear coupled governing equations with the help of an unconditionally stable Crank-Nicolson scheme for the transient couple stress fluid flowing over a vertical cylinder. They observed that the time taken for the flow to reach steady state increases as the Schmidt and Prandtl values increase and decreases with respect to the buoyancy ratio parameter. Furthermore, Ahmed et al. [32] employed the Crank-Nicolson finite difference scheme to the conservative equations in modelling porous media transport for magnetohydrodynamic unsteady flow. They found that the flow velocity and temperature decrease with an increase in the Darcian drag force. The concept of the Crank-Nicolson scheme combined with the Newton-Raphson method was used by Qin et al. [33] to model the heat flux and to estimate the evaporation in applied hydrology and meteorology. The Crank-Nicolson method was used to expand the differential equations whereas the iterative Newton-Raphson method was used to approximate latent heat flux and surface temperatures. Both these methods proved to be successful. In a similar context, Janssen et al. [34] implemented the Crank-Nicolson scheme to transform a system of differential equations into algebraic equations. The Newton-Raphson method is used to implicitly enhance the model’s efficiency by improving the poor convergence rate. Once again, it can be seen that the Crank-Nicolson scheme with the Newton-Raphson method has not been implemented for heat transfer in a porous radial fin.
As far as we know, there has been no or very little work that has been done on obtaining asymptotic solutions or employing a dynamical systems analysis to the problem presented in this research. The purpose of the asymptotic solution is to reveal the dominant physical mechanisms of the model. It can be seen in Moitsheki and Harley [16, 35] that the impact of the thermogeometric parameter () in terms of its proportionality to the length of the fin () was observed. They found that the heat transfer in the fin seemed to be unstable for small values of () due to the fact that . By investigating the asymptotic solution to the steady state heat transfer in a rectangular longitudinal fin, they were able to validate the above relationship and establish the importance of the fin length. Furthermore, the same authors [16, 35] conducted a small scale dynamical analysis. In order to expand the analysis in [35], Harley [36] employed an in-depth dynamical analysis to monitor the behaviour especially at the fin tip. This dynamical analysis also served as a means of investigating the role and effect of the thermogeometric parameter. In this work we do not derive an asymptotic solution; such a solution was derived for the problem but we did not deem it useful in terms of providing deeper insight into the dynamics observed. We will however employ a dynamical systems analysis, which we find to be of immense use in classifying the dynamics of specific points of the solution and the engagement between parameters.
A vast amount of work has been conducted on the steady state case of problems in this field since analytical methods lend themselves more easily to the solution of these equations. In this research, a time dependent partial differential equation modelling the heat transfer in a porous radial fin will be considered. The equation will be derived and nondimensionalised appropriately in Section 2. As a means of comparison to the computational methods employed we structure a semianalytical solution via the Differential Transformation Method in Section 3. In the sections to follow, Sections 4 and 5, the partial differential equation is solved numerically using the Crank-Nicolson scheme combined with the Newton-Raphson method as a predictor-corrector. In order to assess the effectiveness of this scheme and its limitations we employ a dynamical systems analysis in Section 6; in this manner we are able to engage with the limitations placed on parameter values with regard to obtaining solutions. Section 7 provides insight into the comparative dynamics and stability of the equation obtained via an alternate means of nondimensionalisation. Concluding remarks are made in Section 8.
The importance of this work relates to the detailed analysis of the dynamics of the temperature at the fin tip. Furthermore, we are able to investigate the impact of the nonlinear source terms on the stability of the schemes employed and the solutions that can be obtained. This highlights the care than needs to be taken when choosing to solve equations of this nature, particularly when solving via numerical tools. Parameters of physical importance are shown to have value limitations due to stability requirements. Consequently, while numerical solutions may have been obtained here and by other authors noting that these do not indicate an ability to obtain solutions for all relevant parameter values or that it may be assumed that solutions that seemed to have converged are indeed dynamically stable and physically accurate is needed.
2. Model Derivation
Consider a cylindrical porous radial fin with base radius , tip radius , and thickness as shown in Figure 1. The fin comprises an effective thermal conductivity porous material and permeability . It is assumed that the tip of the fin is adiabatic (i.e., a process that occurs without the transfer of heat/matter between the system and its surroundings) and the base of the fin is maintained at a constant temperature . The internal energy per unit volume with the absolute temperature is denoted by . In accordance with Darcy’s law, the fin makes contact with an ambient fluid which infiltrates the fin. The ambient fluid comprises an effective density of the porous fin , a specific heat of the porous fin , the kinematic viscosity of the ambient fluid , the thermal conductivity of the ambient fluid , and the volumetric thermal expansion coefficient of the ambient fluid . The top and bottom surfaces are presumed to have a constant surface emissivity and emit radiation to the ambient fluid at temperature . This also serves as the radiation heat sink.
Constructing the energy balance of a fin element (Figure 1) of thickness , circumference , and radial width at position , we obtainwherewhere (2) is the time rate of change of internal energy in the volume element, (3) is a basis of the application of Fourier’s law of heat conduction, (4) denotes the heat loss due to radiation through the top and bottom surfaces of the fin, and (5) denotes the heat loss due to convection due to Darcy’s flow through a porous fin. At any radiation location , the velocity of the buoyancy driven flow is obtained by the implementation of Darcy’s as follows:
The Darcy model stimulates the fluid-solid interaction in the porous medium. By substituting (2)–(6) into (1), the following nonlinear partial differential equation governing the temperature distribution in the fin is obtained:
The boundary conditions are given by
The fin is initially kept at the temperature of the fluid (ambient temperature) which is given by
We introduce the following nondimensional quantities:
Upon appropriate substitution and rearrangement (8) becomeswhere is defined as and the parameters and are stated as per (12).
Employing the nondimensional transformations given by (12) we find that at we obtain and at we obtain . Thus, we now work on the interval with . As such the nondimensional boundary conditions are given bywith the initial condition
Equation (13) is a nonlinear partial differential equation comprising two nonlinear terms. The first nonlinear term refers to the buoyancy or natural convection transport of energy by the infiltrate whereas the second nonlinear term refers to the surface radiative heat transfer from the fin to the ambient fluid. The parameter is a combination of Rayleigh number , Darcy’s number , the thermal conductivity ratio , and the ratio of the fin tip radius to fin thickness. The parameter indicates the role of surface radiation to conduction in the fin. The parameter is the ratio of the ambient fluid temperature and the difference between the base temperature and the fluid temperature .
For the rest of this research we drop the hat on the independent variable for simplicity. Furthermore, the model employed closely follows the work by Darvishi et al. [2]; in order to be able to compare the results obtained by those authors we use similar parameter values for , , , and . The value of is however defined as the inverse of in the investigations of Darvishi et al. [2] (in their work they define this parameter as ). To allow for a means of comparison we prescribe values for as .
Remark 1. The nondimensional ambient temperature is defined as in the work of Darvishi et al. [2]. Their nondimensionalisation has also been implemented in this research (see Section 7); however a dynamical systems analysis led to the conclusion that there are limitations on certain parameters when obtaining solutions. As such, we altered our nondimensionalisation to define as , as above, as a means of comparison. This change in ratio impacts on the behaviour observed by the solution with regard to changes in the value of . Instead of a linear relationship between and we have a nonlinear one given by , such that the behaviour inversely corresponds to that of Darvishi et al. [2]. As a means of checking the method and analysis conducted in this research, we compared solutions and the dynamics of these solutions to those obtained via the alternate nondimensionalisation employed in [2]; in Section 7 we will briefly comment on the results.
3. Differential Transformation Method
The DTM is an approximation method based on the Taylor series expansion that can easily be applied to several linear and nonlinear problems. It constructs an analytical solution in terms of a polynomial and is capable of reducing the size of the computational work required for solution of the equation. The predominant advantage of this method is that it can be applied directly to linear and nonlinear ordinary differential equations without requiring linearisation, discretisation, or perturbation [37]. This method has been used in a similar manner by Ghasemi et al. [38] to solve the nonlinear temperature distribution equation with temperature dependent internal heat generation in porous and solid longitudinal fins. Furthermore, Fidanoglu et al. [39] used the method to construct a general expression for the heat distribution profile of a fin with spine geometry. Based on the fin effectiveness, entropy values, and efficiency, it was found that the cylindrical fin had the highest rate of heat transfer.
An arbitrary analytical function in the domain can be expanded via the Taylor series about a point as
The differential transform of can thus be defined aswith the inverse differential transform as
The fundamental theorems of the one-dimensional DTM, some of which will be implemented in this paper, are given as follows.
Theorem 2. If , then
Theorem 3. If , then , where is a constant.
Theorem 4. If , then
Theorem 5. If , then
Theorem 6. If , then , where
Theorem 7. If , then
Through the use of these properties of the one-dimensional DTM, we obtain the following iterative procedure for the steady state case of (13):
This method requires initial conditions; that is, we require the values of and since the equation is of second order. We notice that even though we have an initial value for the first derivative, we do not have a value for at . Therefore we specify our initial conditions for iterative procedure equation (4) as follows:where the initial condition for is given an assumed value while the first derivative condition becomes . Through the use of (21) and the conditions given by (22) we are able to obtain which allows us upon substitution into (18) to obtainwhere . However, in order to determine the solution given by (24) we need to first find the value of . In order to determine a value for , we could use the false position method. The boundary condition which has not yet been employed is . Thus we substitute the values obtained as per (23) into (24) for . This provides us with the equality , where can be written as and hence we need to solve the equation . By plotting and we are able to obtain an interval , where and are two guesses to , within which the two functions are equal. Using these values we are able to obtain an improved value which satisfies as follows:
This procedure iterates until a chosen error tolerance is reached such that . We employed this method as well as NSolve in Mathematica as a means of verifying the accuracy of the value. As a brief example, we consider a case where , , , , and to obtain resulting in
In this work however, we have employed a polynomial of order as our steady state solution; see Section 5.2 for discussions on solutions obtained. We are able to obtain a polynomial of degree via this methodology; however for any degree larger we find that limited memory becomes a stumbling block for Mathematica. Furthermore, for we find that our tip temperature whereas for , ; thus increasing the degree of the polynomial does not add much in terms of accuracy. This discrepancy is also not of extreme concern since this semianalytical solution is not our only means of comparison for the numerical solutions obtained below (we also consider the work by Darvishi et al. [2]). We can remark here that many terms are likely to be needed in order for the series solution to match solutions obtained in [2]; however it may still be used as a reasonable measure for the numerical solutions we obtain.
4. Numerical Investigation
4.1. Crank-Nicolson Scheme
The Crank-Nicolson method [40] is a finite difference method which is implicit in time, often provides unconditional stability, and has a high order of accuracy. Given the method’s advantages we will implement the scheme for the equation under consideration. The complexity of our problem lies in the nonlinear sink terms which are likely to influence not only the stability of the scheme, but also the order of accuracy it is able to obtain. We will at first employ the scheme with an explicit consideration of the nonlinear terms. In the subsection to follow we will develop the method further as a means of comparison.
In terms of the finite difference schemes, we approximate using the forward time difference approximation at which giveswhere is the step size in the time domain . Furthermore we approximate the spatial terms and using the central difference approximations at givingrespectively, where is the step size in the spatial domain .
The Crank-Nicolson scheme makes use of approximations to values halfway between nodes by taking averages over two time points. Thus the finite difference approximations on the spatial terms are given by
Given our nondimensional boundary conditions (14)-(15) and the finite difference approximation given by (27) we find that
The initial condition is given as
Having employed the Crank-Nicolson scheme on spatial terms and a forward difference approximation for the time derivative on (13), we have
It should be noticed that we have not applied the Crank-Nicolson method to the nonlinear terms. This is for simplicity at present; however this will be revisited in the next subsection. From (32) we thus obtainwhere we have definedsuch thatfor the sake of simplicity. This equation can be given in the following matrix notation:where and are the relevant coefficient matrices and is defined as the vector of nonlinear terms only.
4.2. Crank-Nicolson Scheme: Newton-Raphson Implementation
In numerical computation, a predictor-corrector method is an algorithm that implements two steps. First, the prediction step calculates a rough approximation of the desired quantity (in this case the nonlinear term on the th time level) whereas the second step refines the initial approximation using other means. In this research we will require such a technique, known as the Newton-Raphson method, given the presence of the nonlinear sink terms. The Crank-Nicolson method incorporating the Newton-Raphson method has been used to determine the solution of a nonlinear diffusion equation (see Kouhia [41] and Habibi et al. [42]) to model the electron heat transfer process in laser inertial fusion for the energy transport mechanism from the region of energy decomposition into the ablation surface. The method proved to be efficient as per the numerical results obtained and hence is a natural extension of the method employed in the previous subsection. We implement the Newton-Raphson method [43] as done by Britz et al. [44] and discussed in [45]. This method is iterative and has been known for its fast convergence to the root as well as its dependence on the initial condition. Furthermore, in order to find the explicit inverse of a general tridiagonal matrix, we shall implement Usmani’s method [46]; the reasoning behind this decision will be discussed at a later stage.
In this subsection we will now employ the Crank-Nicolson method for the nonlinear sink terms as well. As such we rewrite our scheme (33) as follows:where we again use the definitions given by (34)-(35). In matrix notation,where, once again, and are the relevant coefficient matrices and and are defined as the vectors of nonlinear terms only.
We obtain an approximation for the nonlinear vector at the th time level in an iterative fashion as done by Britz et al. [44]. The systemis solved where is the set of difference equations created as per the right hand side of (37) and is the Jacobian of the left hand side of (37). The starting vector at is chosen as per initial condition (15) such that . We solve for using (39) and then iterate as follows:until convergence is reached which provides us with an approximation to . The Newton-Raphson iteration converges within 2-3 steps given that changes are usually relatively small.
The convergence rate of the Crank-Nicolson method incorporating the Newton-Raphson method was found to be extremely slow, making the method infeasible. We discovered that this was due to the fact that the Jacobian J is either close to being singular or badly scaled. As such, we decided to implement Usmani’s method, which is a simple algorithm for finding the explicit inverse of a general Jacobi tridiagonal matrix [46]. This implementation provided the expected results in terms of obtaining the inverse of the Jacobian. Nonetheless, the computational time taken by the method to reach convergence proved inefficient.
5. Remarks
5.1. Comparison of Numerical Schemes
In this subsection we compare the results we obtained via our two numerical methods, that is, the Crank-Nicolson scheme (CN) and the Crank-Nicolson scheme with Newton-Raphson predictor-corrector and Usmani method (). It is important to note that only the first scheme matches the results and behaviour obtained and observed in Darvishi et al. [2]. While the second scheme provides appropriate behaviour, the values at the origin indicate clear discrepancies with the work conducted in [2] even for large values of , where . Furthermore, the length of computational time required for near convergence exceeds that required when the CN scheme is implemented which converges swiftly. As such, aside from the discrepancies in terms of the solutions obtained, we conclude that the method is not efficient in terms of running time.
In order to make a comparison between the two numerical schemes we have employed, we allow our numerical schemes to iterate for large time, termed , so that we allow for convergence to a steady state. In this instance we find that convergence has occurred for the Crank-Nicolson scheme at ; however this is not the case for the Crank-Nicolson method with the predictor-corrector. To maintain our CFL number as less than one for stability purposes, we choose and . As can be surmised, the consequence of this is that the Crank-Nicolson scheme which employs the Newton-Raphson method as a means of updating the nonlinear terms requires a lot of computational time.
We find that while the two numerical schemes produce solutions similar in shape and behaviour, there are discrepancies when we consider the errors between the methods. The temperature at the origin, that is, , is not the same across the two methods; see Tables 1–4. From a consideration of the temperature at the fin tip we find that the solutions do not converge to the same value of for set parameter values and furthermore do not converge at the same rate. At the difference in the values of across the schemes and the maximum error (obtained via comparisons between the numerical schemes) diminish, however not to the extent that the difference is numerically inconsequential. The convergence rate we suspect is heavily influenced by the fact that the Jacobian required for the predictor-corrector method is close to singular or badly scaled. In Figures 2 to 4 we employ the Crank-Nicolson method without the predictor-corrector method due to the speed with which it is able to reach convergence and its accuracy upon comparison with the work by Darvishi et al. [2].
While there are discrepancies we can see behaviour patterns emerging from the different solutions. In Figure 2 as well as Table 1 we observe that an increase in the dimensionless ambient temperature leads to a decrease in the tip temperature; however this can be misleading. We note, as before, that the dimensionless ambient temperature is defined as in the work of Darvishi et al. [2]. In this work we altered our nondimensionalisation to define as , as above. This change in ratio impacts on the behaviour observed by the solution with regard to changes in the value of . Instead of a linear relationship between and we have a nonlinear one given by , such that the behaviour inversely corresponds to that of Darvishi et al. [2]. Hence, our results are consistent with the physical dynamics of the problem.
A consideration of Figure 3 and Table 2 shows the effect of the natural convection heat loss () on the temperature distribution in the fin for set parameter values. As the buoyancy effect increases, that is, , the local temperature in the fin decreases. This is indicative of the effect of natural convective heat loss on the temperature distribution in the fin when the radiation heat loss, the ambient temperature, and the ratio of inner to outer radius are kept fixed. Since buoyancy is principally a macroscale effect we see that the buoyancy force () influences the velocity and temperature fields. Similarly, for Figure 4 and Table 3 the increasing surface radiative heat transfer from the fin to the ambient fluid, that is, , leads to a decrease in the local fin temperature. Furthermore, the impact of radiation () influences the temperature field by increasing the heat transfer rates from the surface. We find that when the convection parameter () and the radiation parameter () are allowed to vary, the local fin temperature decreases because of the increasing strength of radiative heat exchange between the exposed surface of the fin and the ambient [2].
5.2. Comparison of Numerical Results to Semianalytical Solution
We will now compare the results we obtained via our two numerical methods against the series solution obtained via our steady state DTM. We calculate a polynomial of degree via the DTM and also consider the work in [2] to be a benchmark. We find that all three solutions perform well; however only the Crank-Nicolson scheme performs as accurately and efficiently as required.
It can be seen in Table 5, for increasing , that both the value of for the DTM and the maximum error produced between the DTM and the CN scheme increase, whereas the error between the DTM and the scheme decreases. Table 6 indicates that the maximum difference in the temperature between the DTM and the CN scheme with and without the Newton-Raphson method, respectively, agrees to the first decimal place only for ; the CN method performs well across the various values of this parameter. Furthermore, this table shows that for increasing values of the error between the DTM and the CN method decreases when ; in turn Table 7 indicates that for increasing values of at we obtain an increase in the error. These relationships are inverted when considering the comparison between the semianalytical solution and the scheme: the error increases for increasing values of at and decreases for increasing values of at . These occurrences could be indicative of the impact of the parameters and on the CN and methods, respectively. For the CN scheme an increase in the values of also leads to an increase in the error, while for the scheme the inverse is true. We thus observe that an increase in the values of and (which are both contained in only one of the sink terms) leads to an increase in the error between the DTM and the CN scheme. In turn, an increase in the value of , which is the coefficient of the second sink term, leads to an increase in the error between the DTM and the method.
The impact of a change in the value of provides behaviour which is quite stable; see Table 8. The error decreases in both instances; however the error for the CN scheme is much less than that obtained with the introduction of the Newton-Raphson predictor-corrector.
The solutions obtained in Darvishi et al. [2] indicate that the CN scheme provides excellent numerical solutions for the problem under consideration. As such, the errors we obtained when comparing our numerical solutions to the DTM are in part due to a poor approximation of this series solution. Due to computational limitations we are able to obtain a limited number of terms before a lack of memory becomes an issue; are very few terms if one considers the work of other authors such as Darvishi et al. [2] and Aziz et al. [29], for instance. A comparison with the results obtained by Darvishi et al. [2] confirms that the CN scheme produces accurate solutions quickly and efficiently, without the requirement of an excessive amount of terms to be calculated and evaluated as per the semianalytical methods which have been employed in the literature. Unfortunately the addition of a Newton-Raphson predictor-corrector slows the scheme down excessively. By increasing one can see that the scheme does tend to the values obtained by the CN method; however in terms of computational time the scheme has shown itself to be inefficient.
6. Dynamical Systems Analysis
6.1. Stability Analysis on the Linear System
A stability analysis of a system allows us to determine whether or not a system is stable or will be stable if perturbed. In this subsection, we implement a von Neumann stability analysis which is a procedure used to analyse the stability of finite difference schemes by Fourier methods as applied to linear differential equations [47]. It can be noted that the method requires the solution to be periodic and a constant spacing . The method does not account for the boundary conditions.
The structure of the equation incorporates nonlinear sink terms; these do not allow for the use of a von Neumann stability analysis as a means of investigating the numerical scheme. As such we consider the equation without these terms; that is, we consider the linear scheme excluding sink terms. The purpose of this is to determine the extent to which the sink terms may impact on the stability of the schemes we have implemented. We have already noted the impact of the parameters and on the relevant schemes with regard to the maximum error obtained when comparing to the DTM. We will now continue to investigate this phenomenon.
Consider the linearised structure of (33) for the Crank-Nicolson scheme:where we definewith and are defined as before (34) and (35), respectively. Upon substitution of , where , into (41) we obtain
Introducing Euler’s identities into (43) and simplifying provide
Further simplification allows us to find the amplification factor
We note at this point that the von Neumann stability analysis is applicable to problems with constant coefficients, and yet is dependent on . However, we may still deem the matrix to be a matrix of constant coefficients since we know the domain of (i.e., ). Hence, upon substitution of for we obtain constant coefficients as required. We also note that .
Given that for and having definedwe find that since and , We also note that at we obtain and when similarly . Now we let , providing the amplification factor
We may conclude that since , without any required conditions on the relevant step sizes or parameters, the linear system is unconditionally stable.
This conclusion indicates that the structures of the coefficient matrices employed should not impact on the stability of the numerical solutions obtained once the nonlinear terms have been incorporated; rather we need to engage with the impact of the sink terms on the stability of the schemes considered. Thus, to fully investigate the stability of the dynamics of the equation, we now turn to a nonlinear dynamical systems analysis. In this manner we aim to gain insight into the impact of the nonlinear terms on the stability of the relevant schemes.
6.2. Nonlinear Dynamical Systems Analysis
Before we are able to ascertain whether a simple linearisation of the steady state equation may be used as a means of analysing the dynamics of the nonlinear system, we need to determine what this linearised system and the relevant equilibrium points are. It must be noted that we consider the case where ; this is due to the physical relevance placed on the dynamics at the fin tip.
Firstly, as per [48], we restructure our steady state second-order ordinary differential equation into a system of first-order ordinary differential equations as follows:so that
We now turn to fixed point solutions or equilibrium solutions with the intention of investigating the dynamics of the system [49]. In the case of this nonlinear, nonautonomous system the fixed points (or equilibrium points) are defined by the vanishing of the vector field; that is,
Theorems for the stability of fixed points of the systemare presented by Coddington and Levinson [50], where is the autonomous part of the Jacobian and a constant matrix (upon substitution of the relevant equilibrium points and parameter values) and the vector function is continuous in and . In order to employ these theorems we write our system of equations (50) in the form (52):
As per (51) we obtain the equilibrium points :
We follow the work of Coddington and Levinson [50] as discussed in Nayfeh and Balachandran [48]. As such, we need only consider the eigenvalues of the matrix from (53) which we can write with general parameter values as follows: to be evaluated at the equilibrium points in order to classify them. From the form they take we can see two possible cases (depending on whether the term inside the square root is positive or negative): a saddle point or a centre.
6.2.1. Case 1: Equilibrium Point (54)
The eigenvalues of at equilibrium point (54) are
We can see that this equilibrium point (given the ranges for the relevant parameter values) represents a saddle point which is unstable; see Figure 5(b).
(a)
(b)
Coddington and Levinson [50] have shown that if at least one eigenvalue of has a real part which is positive, then upon satisfying the requirement that, given any , there exist and such thatthe solution of (52) is not stable. In our case, we find that for the equilibrium point at least one eigenvalue of is positive and furthermore
We need to show that . In an attempt to do so, we maximise the term on the left of the inequality; that is, we consider such that where, given that , we define . This proves our requirement (58), confirming that equilibrium point (54) is not stable.
6.2.2. Case 2: Equilibrium Point (55)
The eigenvalues of at this equilibrium point are lengthy for arbitrary values of , , , and . However, as pointed out above, the second case needs to be for . A variety of values were considered for the parameters , , , and , all leading to eigenvalues of the form expected, which are providing us with a stable centre; see Figure 5(a). At this stage we note that this case requires which is not physically relevant.
6.3. Consideration of Viability of Linearised System
In the previous subsection we were able to determine the equilibrium points and classify them. We will now consider the feasibility of the dynamics of this linearised system describing the nonlinear dynamics of the system under discussion. Furthermore, we expand our stability analysis to consider the physical significance of the nonlinear terms and the degree, if at all, to which they impact on the scheme’s stability. This analysis has been used by Charrier-Mojtabi et al. [51] in the study of Soret-driven convection in a horizontal porous layer. Furthermore, Patel and Meher [52] performed a stability analysis via the fixed point theorem method for the Adomian decomposition Sumudu transformation method, when solving for the heat transfer in a solid and porous fin with temperature dependent internal heat generation.
In order to employ a linearisation approach to study the dynamics of this nonlinear nonautonomous system we consider a Taylor series expansion of near the equilibrium point [53]. This may be constructed as follows:where we need to ensure that the higher order terms (h.o.t.) are small enough for close to , the equilibrium point, for all . A commonly used criterion based on the -norm iswhich specifies that the higher order terms (h.o.t.) of the Taylor series expansion are approximated by zero at the relevant equilibrium point; this allows us to approximate the nonlinear nonautonomous system by a linearised system.
In order to consider the higher order terms for the equilibrium point (54) we recall thatwhich for our system (50) provides us with the following higher order terms , where represents the Hermitian matrices:
We are interested in the equilibrium point (54) due to its physical significance. Thus we now consider , providingwhich implies that should tend to 0 or be approximately 0 in order for us to consider the dynamics of the linear system as a guide for those of the nonlinear nonautonomous system. Thus, given that we need such thatwhich is unrealistic. The only manner in which (67) could hold is if both parameters were to approximate zero; that is, . Thus, we conclude that the higher order terms have an influence on the stability of the system and cannot be ignored; if we are to engage with the linear system as a means of analysis then we can only do so for small values of and . It is for this reason that the phase portraits presented are for and .
6.4. Remarks
We find that for we are in fact referring to the point where , which is a point that would be located at the origin (see the boundary conditions). We find as either or that at the point where . Thus as and increase we would have a situation which reflects this equilibrium point, albeit a rare case.
Under these circumstances the solutions obtained via the CN scheme provide physically meaningless solutions. We obtain negative temperature values, oscillations in the solutions, and particularly exaggerated negative temperature values near . For the CN scheme we note that this behaviour occurs at lower values of than , indicating once again the more prominent impact the former parameter has on the behaviour of the scheme. In a similar fashion the scheme also produces negative temperature values for large values of and , particularly near . Increasing we find that the value of for which we obtain unrealistic temperature profiles diminishes in magnitude, and we observe severe oscillations and negative temperature values. In turn, implementing large values of with the same value now assigned to does not lead to such physically unreasonable behaviour; we are in fact able to obtain realistic solutions. Hence an increase in the parameter values of and when employing the scheme has a noticeable impact on the stability of the scheme, whereas the buoyancy force’s influence is diminished in comparison.
The behaviour of the solutions obtained for large values of the radiation and convection parameters for either scheme make physical sense. The parameter serves as the radiative sink and as such large values thereof would mean that the surface radiation heat transfer would decrease the temperature unrealistically. Similarly , which is the convective heat loss, would also cause a severe drop in the temperature. Thus, the instability of the equilibrium point can be verified via these dynamics and match the results obtained pertaining to the higher order terms of the Taylor series expansion: only if we assume that and are small, will the linear system approximate the nonlinear system.
7. Alternate Nondimensionalisation
The nondimensionalisation employed in Darvishi et al. [2] is given as follows:
Given that we are considering the time dependent case we also incorporate providing the following equation for consideration:where we have once again dropped the hat on for convenience, is defined as , and the parameters and are stated as per (68).
Employing the nondimensional transformations given by (68) we obtain the following boundary conditions:such that . We use the step sizes and due to the domain of . The numerical solutions obtained mirror those provided above. The only difference is the behaviour of the solutions for ; we now find that they are consistent with [2] as expected. Once more, the addition of the predictor-corrector increases the computational time required to reach the steady state; as such we can again conclude that the Crank-Nicolson scheme performs best out of the two numerical methods considered.
In order to consider the dynamics of the equation we structure the system of first-order ordinary differential equations as before:
It must be noted that we consider the case where ; this is due to the physical relevance placed on the dynamics at the fin tip. We once again structure this system as per the work of Coddington and Levinson [50]:so that we may employ theorems for the stability of fixed points of the system by considering the eigenvalues of . We obtain the equilibrium points :
In considering the eigenvalues of from (73), we can write them with general parameter values as follows: to be evaluated at the equilibrium points we want to classify. Again as per the work conducted above we have two cases: a saddle point or a centre.
7.1. Case 1: Equilibrium Point (74)
The eigenvalues of at equilibrium point (74) are which represent an unstable saddle point; see Figures 6(b) and 7(b) corresponding to the saddle point obtained for the equilibrium point for the previous nondimensionalisation.
(a)
(b)
(a)
(b)
This comparison to the equilibrium point obtained via the previous nondimensionalisation is important given the numerical solutions we obtain for certain values of . While for the initial nondimensionalised equation considered we were able to obtain solutions (whose behaviour provided sensible insight into the problem) for a range of values of , we were unable to do so in this instance. For we found that we were unable to obtain numerical solutions for certain values of and . In order to understand why this is the case, we consider the impact of the higher order terms as done before. We obtain the following:
Now as , we havewhich means that for the linear system to approximate the nonlinear dynamics of the system we now require and either or to ≈0. We find that once we have defined these parameters to be relatively small we were able to obtain solutions for the equation under consideration (70) without any difficulty for both the CN and schemes.
In considering the CN scheme, we observe that for the case and we were able to obtain solutions as . Once we found that there were values for and (not deemed large in magnitude) for which we could not obtain solutions. In particular we noticed that as with the temperature profiles oscillated, often attaining negative temperature values. The moment we diminished the value of the ambient temperature we were able to obtain realistic solutions once more. Furthermore, for reasonably small values of and we found that for values of which could not be deemed excessively high () we were unable to obtain solutions to our problem. As we needed and to decrease in value in a comparative fashion: for with we could not obtain a solution to but we could do so for . As such, the impact of the convection parameter on the ability of the CN scheme to obtain coherent solutions is once again confirmed.
With regard to the scheme we find that as for small values of and the scheme is unable to obtain solutions. Furthermore, as and we are unable to obtain solutions for small values of ; in particular we found oscillatory behaviour at . For this scheme isolating the parameter which played the dominant role in determining the stability of the scheme was very difficult. The structure of (70) is such that there is an intricate interaction between the three parameters. This makes it difficult to isolate, under all circumstances, whether solutions can be obtained or not, unless the parameters , , and are chosen to be small.
Given the fact that we obtained an unstable saddle point at this equilibrium value, we can now via the analysis conducted here confirm that for this problem structure there are certain values of for which limitations are placed upon , in order to obtain solutions. We note here that the behaviour of this equilibrium point matches that of (54) considered for the original nondimensionalisation. Once more, it comments on unstable behaviour at the fin tip under very specific values of and .
7.2. Case 2: Equilibrium Point (75)
The eigenvalues of at this equilibrium point are lengthy for arbitrary values of , , , and . However, as before we realise that the eigenvalues would have the form providing us with a stable centre; see Figures 6(a) and 7(a). For small values of and this equilibrium point provides a case which is not of any physical relevance. When we increase the values of these parameters we observe that since we indeed have a physically relevant equilibrium point. We can in this fashion confirm that the boundary dynamics are at the fin tip (where ), while the temperature is such that ; we indeed have stable behaviour as would be expected. We notice however that the value of at which this occurs is less than the value at which the equilibrium point discussed under Case (where the equilibrium point is dependent on the ambient temperature) is observed. This provides us with the understanding that at the fin tip temperature values less than provide stable dynamics (Case : centre), whereas an increase in the temperature at that point starts to lead to unstable dynamics where stable behaviour becomes dependent on the value of the gradient (Case : saddle point).
8. Conclusion
We employed the Crank-Nicolson method with and without the Newton-Raphson predictor-corrector (employing Usmani’s method for the former) to solve a model describing the heat transfer in a porous radial fin. We show that the results obtained via the Crank-Nicolson scheme match those in the work by Darvishi et al. [2]; the Crank-Nicolson method, incorporating the Newton-Raphson and Usmani method, and the DTM do not match these results even though the behaviour observed is appropriate and comparable.
We obtain semianalytical steady state solutions via the DTM which provide a guide for the numerical solutions obtained. We note that in order to obtain the relevant accuracy for this method a large number of terms would need to be employed (as per the work by Darvishi et al. [2] and Aziz et al. [29], e.g.); however the required memory to do so becomes a limiting factor. The accuracy of the Crank-Nicolson scheme when compared to the work by Darvishi et al. [2] meant that it was the most appropriate one to consider in terms of an analysis of the scheme’s stability and the dynamics of the steady state system.
In conducting a von Neumann stability analysis for the linear Crank-Nicolson scheme we obtain unconditional stability. However, the presence of nonlinear sink terms necessitated an investigation into the impact they may have on the stability of the scheme and hence the validity of the solutions obtained. We find that the only manner in which the linearisation of the system of first-order ordinary differential equations (obtained from the steady state equation) is able to describe the dynamics of the nonlinear system is if and approximate zero, that is, are very small. Thus, we conclude that the higher order terms (and hence nonlinear sink terms) have an influence on the stability of the system and cannot be ignored.
Furthermore, via our analysis, we could ascertain the impact of and , showing that as they tend to infinity erratic and unstable behaviour may be observed at the equilibrium point which represents the fin tip (, ). More interestingly, we were able to show that the convection parameter has a more prominent impact on the stability of the solutions obtained via the CN scheme, whereas the radiative parameter and the ambient temperature impact more severely on the stability of the scheme. This behaviour is sensible given the nature of these sink terms. Hence, the dynamics at the fin tip can be said to be stable in general since an instance where the parameters would obtain such values leading to a case where the equilibrium point is representative of the values of the necessary temperature and temperature gradient is rare. Under the circumstances at the fin tip where and we observe that we do have stable dynamics (see Figure 5(b)).
We also compared the dynamics between the nondimensionalised equation introduced here and that considered in [2]. Considering the second form of nondimensionalisation, the equilibrium point is found to be unstable as well. This case would correspond to an instance where the tip temperature is equated to the ambient temperature (to compare to the alternate case, we previously had an ambient temperature which corresponded to zero so this equilibrium point is equivalent to the one mentioned in the paragraph above). As such the dynamics at the fin tip can again be considered to be unstable under certain circumstances; our numerical solutions confirmed this and displayed the intricate interactions between the relevant parameters. We found that the buoyancy force had a more dominant role to play in the schemes’ ability to obtain solutions. This seems reasonable given the nature of the higher order terms obtained.
The interesting additional insight with regard to the latter nondimensionalisation was that for larger values of and we find that , which provided a second physically relevant equilibrium point. The dynamical analysis conducted assisted us in confirming that the boundary dynamics at the fin tip, that is, for and , are indeed stable. Given the other equilibrium point, it would seem however that for cases where we would revert to unstable behaviour at the fin tip. What is interesting to note is that since in the latter nondimensionalised equation engages with both the buoyancy force/convection parameter and the radiation parameter in the equation, the dynamics are more intricate. As the convection increases, small values of the radiation parameter lead to no solution. Then as the latter increases the solution will achieve the ambient temperature at the fin tip.
For both nondimensionalised equations we discovered that under certain circumstances there could be unstable behaviour at the fin tip; the convection and radiation parameters as well as the ambient temperature have a role to play regarding the stability and effectiveness of the numerical scheme employed. Future research will aim to establish the distinct manner in which these parameters interact, so that we may more fully understand the impact of the parameters on solution accuracy and stability.
Nomenclature
| : | Specific heat of the ambient fluid |
| : | Shape factor for radiation heat transfer |
| : | Darcy number |
| : | Acceleration due to gravity |
| : | Effective thermal conductivity of porous fin |
| : | Thermal conductivity of ambient fluid |
| : | Thermal conductivity ratio |
| : | Buoyancy or natural convection parameter |
| : | Auxiliary linear operator |
| : | Radiation parameter |
| : | Base heat flow |
| : | Length of the fin |
| : | Nondimensional base heat flow |
| : | Radial coordinate |
| : | Base radius |
| : | Tip radius |
| : | Nondimensional radius |
| : | Nondimensional ratio of the tip radius to the base radius |
| : | Fin thickness |
| : | Local fin temperature |
| : | Ambient temperature |
| : | Base temperature. |
| : | Thermal diffusivity of ambient fluid |
| : | Volumetric thermal expansion coefficient of the ambient fluid |
| : | Nondimensional temperature |
| : | Surface emissivity of fin |
| : | Stefan-Boltzmann constant |
| : | Kinematic viscosity of the ambient fluid |
| : | Density of the ambient fluid |
| : | Nondimensional ambient temperature. |
Disclosure
Unless otherwise indicated the views expressed in this paper are the authors’ views only and not the views of the DST/NRF Centre of Excellence in the Mathematical and Statistical Sciences.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
C. Harley acknowledges support from the National Research Foundation, South Africa, under Grant no. 91114.