# Heat Transfer in a Porous Radial Fin: Analysis of Numerically Obtained Solutions

Authors: R. Jooma; C. Harley
Journal: Advances in Mathematical Physics
 (2017)

---

## 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.


---

## Body

## 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 (M) in terms of its proportionality to the length of the fin (L) was observed. They found that the heat transfer in the fin seemed to be unstable for small values of (M) due to the fact that M∝L. 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 Section2. 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 radiusrb, tip radius rt, and thickness s as shown in Figure 1. The fin comprises an effective thermal conductivity porous material keff and permeability K. 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 Tb. The internal energy per unit volume with the absolute temperature T is denoted by ρCvT. 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 ρf, a specific heat of the porous fin Cp,f, the kinematic viscosity of the ambient fluid νf, the thermal conductivity of the ambient fluid kf, and the volumetric thermal expansion coefficient of the ambient fluid βf. The top and bottom surfaces are presumed to have a constant surface emissivity and emit radiation to the ambient fluid at temperature Tb. This also serves as the radiation heat sink.Figure 1
Porous radial fin geometry incorporating the energy balance of a cylindrical cross section.Constructing the energy balance of a fin element (Figure1) of thickness s, circumference 2πr, and radial width dr at position r, we obtain(1)qr-qr+dr-qconv-qrad=I,where(2)I=ρCv2πrsdrdTdt,(3)qr-qr+dr=ddr2πrktdTdrdr,(4)qrad=2ϵσFf-a2πrdrT4-Ta4,(5)qconv=ρfvr2πrdrcp,fT-Ta,where (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 r, the velocity of the buoyancy driven flow vr is obtained by the implementation of Darcy’s as follows:(6)vr=gβfT-Tavf.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: (7)ρCvkeff∂T∂t=1r∂∂rr∂T∂r-2ϵσFf-akeffsT4-Ta4-ρfCp,fKgβfkeffsVfT-Ta2,(8)rb≤r≤rt,t≥0.The boundary conditions are given by(9)∂T∂rrb=0,Trt,t=Tb.The fin is initially kept at the temperature of the fluid (ambient temperature) which is given by(10)T0,r=Ta.We introduce the following nondimensional quantities:(11)θ=T-TaTb-Ta,θa=TaTb-Ta,r^=r-rtrb-rt,r^∗=rbrt,τ=kefftρCvrt2,(12)Nc=ρfCp,fKgβfrt2Tb-TaVfkeffs,Nr=2ϵσFf-art2Tb-Ta3keffs.Upon appropriate substitution and rearrangement (8) becomes(13)∂θ∂τ=1r^r^∗-1+1r^∗-1∂∂r^r^r^∗-1+1r^∗-1∂θ∂r^-Ncθ2-Nrθ+θa4-θa4,where r^∗ is defined as rb/rt and the parameters Nc and Nr are stated as per (12).Employing the nondimensional transformations given by (12) we find that at r=rb we obtain r^=0 and at rt we obtain r^=1. Thus, we now work on the interval 0≤r^≤1 with τ≥0. As such the nondimensional boundary conditions are given by(14)∂θ∂r^r^=0=0,θ1,τ=1,with the initial condition(15)θr^,0=0.Equation (13) is a nonlinear partial differential equation comprising two nonlinear terms. The first nonlinear term Nc refers to the buoyancy or natural convection transport of energy by the infiltrate whereas the second nonlinear term Nr refers to the surface radiative heat transfer from the fin to the ambient fluid. The parameter Nc is a combination of Rayleigh number Ra, Darcy’s number Da, the thermal conductivity ratio Kr, and the ratio of the fin tip radius to fin thickness. The parameter Nr indicates the role of surface radiation to conduction in the fin. The parameter θa is the ratio of the ambient fluid temperature and the difference between the base temperature and the fluid temperature Tb-Ta.For the rest of this research we drop the hat on the independent variabler 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 r∗, Nc, Nr, and θa. The value of r∗=rb/rt is however defined as the inverse of r∗ in the investigations of Darvishi et al. [2] (in their work they define this parameter as R∗=rt/rb>1). To allow for a means of comparison we prescribe values for r∗ as 1/R∗.Remark 1.
The nondimensional ambient temperatureθa is defined as Ta/Tb 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 θa as Ta/Tb-Ta, 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 θa. Instead of a linear relationship between Ta and θa we have a nonlinear one given by Ta=θaTb/1+θa, 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 functionfx in the domain D can be expanded via the Taylor series about a point x=0 as(16)fx=∑k=0∞xkk!dkfxdxkx=0.The differential transform offx can thus be defined as(17)Fk=1k!dkfxdxkx=0,with the inverse differential transform as(18)fx=∑k=0∞xkFk.The fundamental theorems of the one-dimensional DTM, some of which will be implemented in this paper, are given as follows.Theorem 2.
Iffx=gx±hx, then Fk=Gk±Hk.Theorem 3.
Iffx=cgx, then Fk=cGk, where c is a constant.Theorem 4.
Iffx=dngx/dxn, then Fk=k+n!/kGk+n.Theorem 5.
Iffx=gxhx, then Fk=∑k1=0k‍Gk1Hk-k1.Theorem 6.
Iffx=xn, then Fk=δk-n, where (19)δk-n=1,k=n,0,k≠n.Theorem 7.
Iffx=g1xg2x⋯gn-1xgnx, then (20)Fk=∑kn-1=0k∑kn-2=0kn-1⋯∑k2=0k3∑k1=0k2G2k1G1k2-k1×⋯×Gn-1kn-1-kn-2Gnkn-kn-1.Through the use of these properties of the one-dimensional DTM, we obtain the following iterative procedure for the steady state case of (13):(21)Θk+2=-r∗-12k+1k+2k+1Θk+1r∗-1r∗-1+1δk-1-Nc∗∑k1=0kΘk1Θk-k1-Nr∗∑k4=0k∑k3=0k4∑k2=0k3∑k1=0k2Θk1Θk2-k1Θk3-k2Θk4-k3Θk-k4+4θa∑k3=0k∑k2=0k3∑k1=0k2Θk1Θk2-k1Θk3-k2Θk-k3Θk-k4+6θa2∑k2=0k∑k1=0k2Θk1Θk2-k1Θk-k2+4θa3Θk-θa4.This method requires initial conditions; that is, we require the values ofθ0,τ and ∂θ0,τ/∂r 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 r=0. Therefore we specify our initial conditions for iterative procedure equation (4) as follows:(22)Θ0=Θ0,Θ1=0,where the initial condition for θ is given an assumed value while the first derivative condition becomes Θ1=0. Through the use of (21) and the conditions given by (22) we are able to obtain (23)Θ0=-12-1+r∗2-Θ02Nc-NrΘ05+4Θ04θa+6Θ03θa2+4Θ0θa3-θa4,Θ1=-16-1+r∗2Nrθa4--1+r∗-Θ02Nc-NrΘ05+4a4θa+6Θ03θa2+4Θ0θa3-θa4,Θ2=-112-1+r∗2Θ0Nc-1+r∗2-Θ02Nc-NrΘ05+4Θ04θa+6Θ03θa2+4Θ0θa3-θa4-12r∗-1+r∗Nrθa4--1+r∗-Θ02Nc-NrΘ05+4Θ04θa+6Θ03θa2+4Θ0θa3-θa4-Nr-θa4-52Θ04-1+r∗2-Θ02Nc-NrΘ05+4Θ04θa+6Θ03θa2+4Θ0θa3-θa4-8Θ03-1+r∗2θa-Θ02Nc-NrΘ05+4Θ04θa+6Θ03θa2+4Θ0θa3-θa4-9Θ02-1+r∗2θa2-Θ02Nc-NrΘ05+4Θ04θa+6Θ03θa2+4Θ0θa3-θa4-2-1+r∗2θa3-Θ02Nc-NrΘ05+4Θ04θa+6Θ03θa2+4Θ0θa3-θa4,which allows us upon substitution into (18) to obtain(24)θr=∑k=0nr-r0kΘk,where r0=0. However, in order to determine the solution given by (24) we need to first find the value of Θ0. In order to determine a value for Θ0, we could use the false position method. The boundary condition which has not yet been employed is θ1=1. Thus we substitute the values obtained as per (23) into (24) for r=1. This provides us with the equality WΘ0=1, where WΘ0 can be written as vΘ0/uΘ0 and hence we need to solve the equation ZΘ0=vΘ0-uΘ0=0. By plotting vΘ0 and uΘ0 we are able to obtain an interval Θ0a,Θ0b, where Θ0a and Θ0b are two guesses to Θ0, within which the two functions are equal. Using these values we are able to obtain an improved value which satisfies ZΘ0=0 as follows:(25)Θ0new=Θ0bZΘ0b-Θ0aZΘ0aZΘ0b-ZΘ0a.This procedure iterates until a chosen error tolerance is reached such thatZΘ0new≪ϵ. 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 Nc=10, Nr=1, r∗=1/1.75, θa=0.2, and n=3 to obtain Θ0=0.564821 resulting in(26)θr=0.564821+0.311221r2+0.0444111r3+0.0685314r4+0.0110162r5.In this work however, we have employed a polynomial of ordern=10 as our steady state solution; see Section 5.2 for discussions on solutions obtained. We are able to obtain a polynomial of degree n=16 via this methodology; however for any degree larger we find that limited memory becomes a stumbling block for Mathematica. Furthermore, for n=10 we find that our tip temperature θ0=0.55725 whereas for n=16, θ0=0.55678; 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 τm which gives(27)∂θ∂τ≈θnm+1-θnmΔτ,where Δτ is the step size in the time domain τ≥0. Furthermore we approximate the spatial terms ∂θ/∂r and ∂2θ/∂r2 using the central difference approximations at rn giving(28)∂θ∂r≈θn+1m-θn-1m2Δr,∂2θ∂r2≈θn+1m-2θnm+θn-1mΔr2,respectively, where Δr is the step size in the spatial domain r∈0,1.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(29)∂θ∂r≈θn+1m+1-θn-1m+14Δr+θn+1m-θn-1m4Δr,∂2θ∂r2≈θn+1m+1-2θnm+1+θn-1m+12Δr2+θn+1m-2θnm+θn-1m2Δr2.Given our nondimensional boundary conditions (14)-(15) and the finite difference approximation given by (27) we find that(30)∂θ∂rr0=0⟹θ1m=θ-1m∀m,θ1,τ=1⟹θNm=1∀m,whererN=1.The initial condition is given as(31)θn0=0∀n.Having employed the Crank-Nicolson scheme on spatial terms and a forward difference approximation for the time derivative on (13), we have(32)θnm+1-θnmΔτ=1rnr∗-1+1r∗-1θn+1m+1-θn-1m+14Δr+θn+1m-θn-1m4Δr+1r∗-12θn+1m+1-2θnm+1+θn-1m+12Δr2+θn+1m-2θnm+θn-1m2Δr2-Ncθnm2-Nrθnm+θa4-θa4.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 obtain(33)-αnθn+1m+1+1+δθnm+1-βnθn-1m+1=αnθn+1m+1-δθnm+βnθn-1m-ΔτNcθnm2-ΔτNrθnm+θa4-θa4,where we have defined(34)γn=1rnr∗-1+1r∗-1,such that(35)αn=Δτγn4Δr+Δτ2Δr2r∗-12,δ=ΔτΔr2r∗-12,βn=-Δτγn4Δr+Δτ2Δr2r∗-12,for the sake of simplicity. This equation can be given in the following matrix notation:(36)Aθ_m+1=Bθ_m+Gθ_m,where A and B are the relevant coefficient matrices and Gθ_m 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 them+1th 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:(37)-αnθn+1m+1+1+δθnm+1-βnθn-1m+1+ΔτNc2θnm+12+ΔτNr2θnm+1+θa4-θa4=αnθn+1m+1-δθnm+βnθn-1m-ΔτNc2θnm2-ΔτNr2θnm+θa4-θa4,where we again use the definitions given by (34)-(35). In matrix notation,(38)Aθ_m+1-12Gθ_m+1=Bθ_m+12Gθ_m,where, once again, A and B are the relevant coefficient matrices and Gθ_m and Gθ_m+1 are defined as the vectors of nonlinear terms only.We obtain an approximation for the nonlinear vector at the(m+1)th time level in an iterative fashion as done by Britz et al. [44]. The system(39)Jθ_δθ_=-Fθ_is solved where F is the set of difference equations created as per the right hand side of (37) and J is the Jacobian of the left hand side of (37). The starting vector at τ=0 is chosen as per initial condition (15) such that θ_0=0. We solve for δθ_k using (39) and then iterate as follows:(40)θ_k+1=θ_k+δθ_k,until convergence is reached which provides us with an approximation to θ_m+1. 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 JacobianJ 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.


## 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 τm which gives(27)∂θ∂τ≈θnm+1-θnmΔτ,where Δτ is the step size in the time domain τ≥0. Furthermore we approximate the spatial terms ∂θ/∂r and ∂2θ/∂r2 using the central difference approximations at rn giving(28)∂θ∂r≈θn+1m-θn-1m2Δr,∂2θ∂r2≈θn+1m-2θnm+θn-1mΔr2,respectively, where Δr is the step size in the spatial domain r∈0,1.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(29)∂θ∂r≈θn+1m+1-θn-1m+14Δr+θn+1m-θn-1m4Δr,∂2θ∂r2≈θn+1m+1-2θnm+1+θn-1m+12Δr2+θn+1m-2θnm+θn-1m2Δr2.Given our nondimensional boundary conditions (14)-(15) and the finite difference approximation given by (27) we find that(30)∂θ∂rr0=0⟹θ1m=θ-1m∀m,θ1,τ=1⟹θNm=1∀m,whererN=1.The initial condition is given as(31)θn0=0∀n.Having employed the Crank-Nicolson scheme on spatial terms and a forward difference approximation for the time derivative on (13), we have(32)θnm+1-θnmΔτ=1rnr∗-1+1r∗-1θn+1m+1-θn-1m+14Δr+θn+1m-θn-1m4Δr+1r∗-12θn+1m+1-2θnm+1+θn-1m+12Δr2+θn+1m-2θnm+θn-1m2Δr2-Ncθnm2-Nrθnm+θa4-θa4.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 obtain(33)-αnθn+1m+1+1+δθnm+1-βnθn-1m+1=αnθn+1m+1-δθnm+βnθn-1m-ΔτNcθnm2-ΔτNrθnm+θa4-θa4,where we have defined(34)γn=1rnr∗-1+1r∗-1,such that(35)αn=Δτγn4Δr+Δτ2Δr2r∗-12,δ=ΔτΔr2r∗-12,βn=-Δτγn4Δr+Δτ2Δr2r∗-12,for the sake of simplicity. This equation can be given in the following matrix notation:(36)Aθ_m+1=Bθ_m+Gθ_m,where A and B are the relevant coefficient matrices and Gθ_m 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 them+1th 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:(37)-αnθn+1m+1+1+δθnm+1-βnθn-1m+1+ΔτNc2θnm+12+ΔτNr2θnm+1+θa4-θa4=αnθn+1m+1-δθnm+βnθn-1m-ΔτNc2θnm2-ΔτNr2θnm+θa4-θa4,where we again use the definitions given by (34)-(35). In matrix notation,(38)Aθ_m+1-12Gθ_m+1=Bθ_m+12Gθ_m,where, once again, A and B are the relevant coefficient matrices and Gθ_m and Gθ_m+1 are defined as the vectors of nonlinear terms only.We obtain an approximation for the nonlinear vector at the(m+1)th time level in an iterative fashion as done by Britz et al. [44]. The system(39)Jθ_δθ_=-Fθ_is solved where F is the set of difference equations created as per the right hand side of (37) and J is the Jacobian of the left hand side of (37). The starting vector at τ=0 is chosen as per initial condition (15) such that θ_0=0. We solve for δθ_k using (39) and then iterate as follows:(40)θ_k+1=θ_k+δθ_k,until convergence is reached which provides us with an approximation to θ_m+1. 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 JacobianJ 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 (CNNR). 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 T, where τ∈0,T. 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 CNNR 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, termedT, 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 T=1; 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 Δr=0.1 and Δτ=0.0001. 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,θ0, 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 θ0 for set parameter values and furthermore do not converge at the same rate. At T=100 the difference in the values of θ0 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].Table 1
Temperature atr=0 for varied values of θa.
θa
θ
0
C
N
θ
0
N
R
θ0DTM
0.2
0.54482
0.67327
0.55725
0.4
0.51953
0.64956
0.54037
0.6
0.48191
0.61510
0.51752
0.8
0.43321
0.57044
0.49348
Parameter values:Nc=10,  Nr=1,  r∗=1/1.75, and  T=1.Table 2
Maximum error computed between the respective numerical solutions obtained, whereE=max⁡CNsol-CNsolNR, and their respective temperatures at the tip of the fin for varied values of Nc.
Nc
E
θ
0
C
N
θ
0
N
R
1
0.07624
0.79980
0.87604
10
0.12845
0.54482
0.67327
50
0.11464
0.27244
0.38697
Parameter values:θa=0.2,  Nr=1,  r∗=1/1.75, and  T=1.Table 3
Maximum error computed between the respective numerical solutions obtained, whereE=max⁡CNsol-CNsolNR, and their respective temperatures at the tip of the fin for varied values of Nc.
N
r
E
θ
0
C
N
θ
0
N
R
5
0.10318
0.62343
0.72661
10
0.10698
0.52686
0.63384
50
0.09496
0.30490
0.39982
Parameter values:θa=0.2,  r∗=1/1.75,  Nc=1, and  T=100.Table 4
Maximum error computed between the respective numerical solutions obtained, whereE=max⁡CNsolCN-CNsolNR, and their respective temperatures at the tip of the fin, given as θ0, for varied values of r∗.
r
∗
E
θ
0
C
N
θ
0
N
R
0.2
0.10849
0.27903
0.38752
0.4
0.12751
0.40088
0.52840
0.6
0.12615
0.57382
0.69998
0.8
0.07434
0.81944
0.89379
Parameter values:θa=0.2,  Nr=1,  Nc=10, and  T=1.Figure 2
Plot of the temperature withNr=1, Nc=10, and r∗=1/1.75 for T=1 and varied values of θa.While there are discrepancies we can see behaviour patterns emerging from the different solutions. In Figure2 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 θa is defined as Ta/Tb in the work of Darvishi et al. [2]. In this work we altered our nondimensionalisation to define θa as Ta/Tb-Ta, as above. This change in ratio impacts on the behaviour observed by the solution with regard to changes in the value of θa. Instead of a linear relationship between Ta and θa we have a nonlinear one given by Ta=θaTb/1+θa, 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 Figure3 and Table 2 shows the effect of the natural convection heat loss (Nc) on the temperature distribution in the fin for set parameter values. As the buoyancy effect increases, that is, Nc, 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 (Nc) 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, Nr, leads to a decrease in the local fin temperature. Furthermore, the impact of radiation (Nr) influences the temperature field by increasing the heat transfer rates from the surface. We find that when the convection parameter (Nc) and the radiation parameter (Nr) 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].Figure 3
Plot of the temperature withθa=0.2, r∗=1/1.75, and Nr=1 for T=1 and varied values of Nc.Figure 4
Plot of the temperature withθa=0.2, r∗=1/1.75, and Nc=1 for T=100 and varied values of Nr.


### 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 degreen=10 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 Table5, for increasing θa, that both the value of θ0 for the DTM and the maximum error produced between the DTM and the CN scheme increase, whereas the error between the DTM and the CNNR 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 Nc=1; the CN method performs well across the various values of this parameter. Furthermore, this table shows that for increasing values of Nc the error between the DTM and the CN method decreases when T=5; in turn Table 7 indicates that for increasing values of Nr at T=100 we obtain an increase in the error. These relationships are inverted when considering the comparison between the semianalytical solution and the CNNR scheme: the error increases for increasing values of Nc at T=5 and decreases for increasing values of Nr at T=100. These occurrences could be indicative of the impact of the parameters Nr and Nc on the CN and CNNR methods, respectively. For the CN scheme an increase in the values of θa also leads to an increase in the error, while for the CNNR scheme the inverse is true. We thus observe that an increase in the values of Nr and θa (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 Nc, which is the coefficient of the second sink term, leads to an increase in the error between the DTM and the CNNR method.Table 5
Numerical solutions obtained for DTM and the maximum error computed between the DTM solution and the respective numerical solutions obtained, whereECN=max⁡CNsol-DTMsol and ENR=max⁡CNsolNR-DTMsol, for varied values of Nc.
θ
a
θ
0
D
T
M
E
C
N
E
N
R
0.2
0.55725
0.01284
0.11602
0.4
0.54037
0.02087
0.10919
0.6
0.51752
0.03561
0.09758
0.8
0.49348
0.06028
0.07696
Parameter values:Nc=10,  Nr=1,  r∗=1/1.75, and  T=5.Table 6
Numerical solutions obtained for DTM and the maximum error computed between the respective numerical solutions obtained, whereECN=max⁡CNsol-DTMsol and ENR=max⁡CNsolNR-DTMsol, and their respective temperatures at the tip of the fin for varied values of Nc.
N
c
θ
0
D
T
M
E
C
N
E
N
R
1
0.81615
0.01635
0.05989
10
0.55725
0.01284
0.11602
50
0.27747
0.01177
0.10950
Parameter values:θa=0.2,  Nr=1,  r∗=1/1.75, and  T=5.Table 7
Numerical solutions obtained for DTM and the maximum error computed between the respective numerical solutions obtained, whereECN=max⁡CNsol-DTMsol and ENR=max⁡CNsolNR-DTMsol, and their respective temperatures at the tip of the fin for varied values of Nc.
N
r
θ
0
D
T
M
E
C
N
E
N
R
1
0.81615
0.01672
0.04911
5
0.66993
0.047281
0.03450
10
0.59043
0.06465
0.01524
Parameter values:θa=0.2,  r∗=1/1.75,  Nc=1, and   T=100.The impact of a change in the value ofr∗ 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.Table 8
Numerical solutions obtained for DTM and the maximum error computed between the respective numerical solutions obtained, whereECN=max⁡CNsol-DTMsol and ENR=max⁡CNsolNR-DTMsol, and their respective temperatures at the tip of the fin, given as θ0, for varied values of r∗.
r
∗
θ
0
D
T
M
E
C
N
E
N
R
0.2
0.30038
0.05459
0.07799
0.4
0.41843
0.02504
0.09325
0.6
0.58528
0.01230
0.08572
0.8
0.82286
0.00361
0.02814
Parameter values:θa=0.2,  Nr=1,  Nc=10, and   T=100.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; n=10 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 T one can see that the scheme does tend to the values obtained by the CN method; however in terms of computational time the CNNR scheme has shown itself to be inefficient.


## 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 (CNNR). 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 T, where τ∈0,T. 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 CNNR 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, termedT, 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 T=1; 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 Δr=0.1 and Δτ=0.0001. 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,θ0, 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 θ0 for set parameter values and furthermore do not converge at the same rate. At T=100 the difference in the values of θ0 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].Table 1
Temperature atr=0 for varied values of θa.
θa
θ
0
C
N
θ
0
N
R
θ0DTM
0.2
0.54482
0.67327
0.55725
0.4
0.51953
0.64956
0.54037
0.6
0.48191
0.61510
0.51752
0.8
0.43321
0.57044
0.49348
Parameter values:Nc=10,  Nr=1,  r∗=1/1.75, and  T=1.Table 2
Maximum error computed between the respective numerical solutions obtained, whereE=max⁡CNsol-CNsolNR, and their respective temperatures at the tip of the fin for varied values of Nc.
Nc
E
θ
0
C
N
θ
0
N
R
1
0.07624
0.79980
0.87604
10
0.12845
0.54482
0.67327
50
0.11464
0.27244
0.38697
Parameter values:θa=0.2,  Nr=1,  r∗=1/1.75, and  T=1.Table 3
Maximum error computed between the respective numerical solutions obtained, whereE=max⁡CNsol-CNsolNR, and their respective temperatures at the tip of the fin for varied values of Nc.
N
r
E
θ
0
C
N
θ
0
N
R
5
0.10318
0.62343
0.72661
10
0.10698
0.52686
0.63384
50
0.09496
0.30490
0.39982
Parameter values:θa=0.2,  r∗=1/1.75,  Nc=1, and  T=100.Table 4
Maximum error computed between the respective numerical solutions obtained, whereE=max⁡CNsolCN-CNsolNR, and their respective temperatures at the tip of the fin, given as θ0, for varied values of r∗.
r
∗
E
θ
0
C
N
θ
0
N
R
0.2
0.10849
0.27903
0.38752
0.4
0.12751
0.40088
0.52840
0.6
0.12615
0.57382
0.69998
0.8
0.07434
0.81944
0.89379
Parameter values:θa=0.2,  Nr=1,  Nc=10, and  T=1.Figure 2
Plot of the temperature withNr=1, Nc=10, and r∗=1/1.75 for T=1 and varied values of θa.While there are discrepancies we can see behaviour patterns emerging from the different solutions. In Figure2 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 θa is defined as Ta/Tb in the work of Darvishi et al. [2]. In this work we altered our nondimensionalisation to define θa as Ta/Tb-Ta, as above. This change in ratio impacts on the behaviour observed by the solution with regard to changes in the value of θa. Instead of a linear relationship between Ta and θa we have a nonlinear one given by Ta=θaTb/1+θa, 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 Figure3 and Table 2 shows the effect of the natural convection heat loss (Nc) on the temperature distribution in the fin for set parameter values. As the buoyancy effect increases, that is, Nc, 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 (Nc) 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, Nr, leads to a decrease in the local fin temperature. Furthermore, the impact of radiation (Nr) influences the temperature field by increasing the heat transfer rates from the surface. We find that when the convection parameter (Nc) and the radiation parameter (Nr) 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].Figure 3
Plot of the temperature withθa=0.2, r∗=1/1.75, and Nr=1 for T=1 and varied values of Nc.Figure 4
Plot of the temperature withθa=0.2, r∗=1/1.75, and Nc=1 for T=100 and varied values of Nr.


## 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 degreen=10 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 Table5, for increasing θa, that both the value of θ0 for the DTM and the maximum error produced between the DTM and the CN scheme increase, whereas the error between the DTM and the CNNR 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 Nc=1; the CN method performs well across the various values of this parameter. Furthermore, this table shows that for increasing values of Nc the error between the DTM and the CN method decreases when T=5; in turn Table 7 indicates that for increasing values of Nr at T=100 we obtain an increase in the error. These relationships are inverted when considering the comparison between the semianalytical solution and the CNNR scheme: the error increases for increasing values of Nc at T=5 and decreases for increasing values of Nr at T=100. These occurrences could be indicative of the impact of the parameters Nr and Nc on the CN and CNNR methods, respectively. For the CN scheme an increase in the values of θa also leads to an increase in the error, while for the CNNR scheme the inverse is true. We thus observe that an increase in the values of Nr and θa (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 Nc, which is the coefficient of the second sink term, leads to an increase in the error between the DTM and the CNNR method.Table 5
Numerical solutions obtained for DTM and the maximum error computed between the DTM solution and the respective numerical solutions obtained, whereECN=max⁡CNsol-DTMsol and ENR=max⁡CNsolNR-DTMsol, for varied values of Nc.
θ
a
θ
0
D
T
M
E
C
N
E
N
R
0.2
0.55725
0.01284
0.11602
0.4
0.54037
0.02087
0.10919
0.6
0.51752
0.03561
0.09758
0.8
0.49348
0.06028
0.07696
Parameter values:Nc=10,  Nr=1,  r∗=1/1.75, and  T=5.Table 6
Numerical solutions obtained for DTM and the maximum error computed between the respective numerical solutions obtained, whereECN=max⁡CNsol-DTMsol and ENR=max⁡CNsolNR-DTMsol, and their respective temperatures at the tip of the fin for varied values of Nc.
N
c
θ
0
D
T
M
E
C
N
E
N
R
1
0.81615
0.01635
0.05989
10
0.55725
0.01284
0.11602
50
0.27747
0.01177
0.10950
Parameter values:θa=0.2,  Nr=1,  r∗=1/1.75, and  T=5.Table 7
Numerical solutions obtained for DTM and the maximum error computed between the respective numerical solutions obtained, whereECN=max⁡CNsol-DTMsol and ENR=max⁡CNsolNR-DTMsol, and their respective temperatures at the tip of the fin for varied values of Nc.
N
r
θ
0
D
T
M
E
C
N
E
N
R
1
0.81615
0.01672
0.04911
5
0.66993
0.047281
0.03450
10
0.59043
0.06465
0.01524
Parameter values:θa=0.2,  r∗=1/1.75,  Nc=1, and   T=100.The impact of a change in the value ofr∗ 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.Table 8
Numerical solutions obtained for DTM and the maximum error computed between the respective numerical solutions obtained, whereECN=max⁡CNsol-DTMsol and ENR=max⁡CNsolNR-DTMsol, and their respective temperatures at the tip of the fin, given as θ0, for varied values of r∗.
r
∗
θ
0
D
T
M
E
C
N
E
N
R
0.2
0.30038
0.05459
0.07799
0.4
0.41843
0.02504
0.09325
0.6
0.58528
0.01230
0.08572
0.8
0.82286
0.00361
0.02814
Parameter values:θa=0.2,  Nr=1,  Nc=10, and   T=100.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; n=10 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 T one can see that the scheme does tend to the values obtained by the CN method; however in terms of computational time the CNNR 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 parametersNr and Nc 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:(41)-αnθn+1m+1+1+δθnm+1-βnθn-1m+1=αnθn+1m+1-δθnm+βnθn-1m,where we define(42)αn=Δτγn4Δr+Δτ2Δr2r∗-12,βn=-Δτγn4Δr+Δτ2Δr2r∗-12,with γn and δ are defined as before (34) and (35), respectively. Upon substitution of θnm=ξmeinΔrω, where i2=-1, into (41) we obtain (43)-αnξm+1ein+1Δrω+1+δξm+1einΔrω-βnξm+1ein-1Δrω=αnξmein+1Δrω+1-δξmeinΔrω+βnξmein-1Δrω.Introducing Euler’s identities into (43) and simplifying provide (44)ξΔτγn4Δr-2isin⁡ωΔr+Δτ2Δr2r∗-12-2cos⁡ωΔr+1+δ=Δτγn4Δr2isin⁡ωΔr+Δτ2Δr2r∗-122cos⁡ωΔr+1-δ.Further simplification allows us to find the amplification factor(45)ξ=1-δ+Δτγn/2Δrisin⁡ωΔr+Δτ/Δr2r∗-12cos⁡ωΔr1+δ-Δτγn/2Δrisin⁡ωΔr-Δτ/Δr2r∗-12cos⁡ωΔr.We note at this point that the von Neumann stability analysis is applicable to problems with constant coefficients, and yetγn is dependent on rn. However, we may still deem the matrix to be a matrix of constant coefficients since we know the domain of rn (i.e., rn∈0,1). Hence, upon substitution of rn for n=0,1,…,N we obtain constant coefficients as required. We also note that 0<δ<1.Given thatr∗=rb/rt<1 for rn∈0,1 and having defined(46)γn=1rnr∗-1+1r∗-1,we find that since r∗-1<0 and rnr∗-1+1>0, γn<0. We also note that at r=0 we obtain γn=1/r∗-1<0 and when r=1 similarly γn=1/r∗r∗-1<0. Now we let γ¯n=-γn>0, providing the amplification factor(47)ξ=1-δ-Δτγ¯n/2Δrisin⁡ωΔx+Δτ/Δr2r∗-12cos⁡ωΔx1+δ+Δτγ¯n/2Δrisin⁡ωΔx-Δτ/Δr2r∗-12cos⁡ωΔx.We may conclude that sinceξ<1, 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 wherer=0; 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:(48)θ_˙=Fθ_,r,so that (49)θ˙1=θ2,(50)θ˙2=1-r∗θ2rr∗-1+1+Nrr∗-12θ1+θa4-θa4+Ncr∗-12θ12.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,(51)Fθ_,r=0_.Theorems for the stability of fixed points of the system(52)θ_˙=Aθ_+fθ_,rare presented by Coddington and Levinson [50], where A is the autonomous part of the Jacobian Dθ_Fθ_,r and a constant matrix (upon substitution of the relevant equilibrium points and parameter values) and the vector function f is continuous in θ_ and r. In order to employ these theorems we write our system of equations (50) in the form (52):(53)θ_˙=014Nrr∗-12θ1+θa3+2Ncr∗-12θ10θ_+01-r∗rr∗-1+1θ2.As per (51) we obtain the equilibrium points θ1,θ2:(54)0,0,(55)18NcNr2θa-10Nr3θa3+5.19615Nc3Nr3+14Nc2Nr4θa2-12NcNr5θa4+4Nr6θa633Nr-3NcNr+2Nr2θa23Nr18NcNr2θa-10Nr3θa3+5.19615Nc3Nr3+14Nc2Nr4θa2-12NcNr5θa4+4Nr6θa63-4θa3,0.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 A from (53) which we can write with general parameter values as follows: (56)λ1=-2r∗-1Ncθ1+2Nrθ13+6Nrθ12θa+6Nrθ1θa2+2Nrθa3,λ2=2r∗-1Ncθ1+2Nrθ13+6Nrθ12θa+6Nrθ1θa2+2Nrθa3,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 ofA at equilibrium point (54) are (57)λ1=-2Nrr∗-1θa3/2,λ2=2Nrr∗-1θa3/2.We can see that this equilibrium point (given the ranges for the relevant parameter values) represents a saddle point which is unstable; see Figure5(b).Figure 5
Phase plots indicating the trajectories around the equilibrium points (54) and (55) for θa=0.2, r∗=1/1.75, and Nr=1. (a) has Nc=0.1 whereas (b) has Nc=0.01. (a) is a closer view of the second equilibrium point.
(a)
(b)Coddington and Levinson [50] have shown that if at least one eigenvalue of A has a real part which is positive, then upon fθ_,r satisfying the requirement that, given any ϵ>0, there exist σ and R such that(58)fθ_,r≦ϵθ_θ_≦σ,r≧R,the solution θ_=0_ of (52) is not stable. In our case, we find that for the equilibrium point 0,0 at least one eigenvalue of A is positive and furthermore(59)fθ_,r=0001-r∗rr∗-1+1θ1θ2.We need to show that1-r∗/rr∗-1+1θ2≦ϵθ2. In an attempt to do so, we maximise the term on the left of the inequality; that is, we consider r=0 such that (60)1-r∗rr∗-1+1θ2≦1-r∗rr∗-1+1θ2≦1-r∗θ2=ϵθ2,where, given that 0<1-r∗<1, we define 0<ϵ<1. This proves our requirement (58), confirming that equilibrium point (54) is not stable.


#### 6.2.2. Case  2: Equilibrium Point (55)

The eigenvalues ofA at this equilibrium point are lengthy for arbitrary values of Nc, Nr, r∗, and θa. However, as pointed out above, the second case needs to be for Ncθ1+2Nrθ13+6Nrθ12θa+6Nrθ1θa2+2Nrθa3<0. A variety of values were considered for the parameters Nc, Nr, r∗, and θa, all leading to eigenvalues of the form expected, which are (61)λ1=μi,λ2=-μi,wherei2=-1,providing us with a stable centre; see Figure 5(a). At this stage we note that this case requires θ1=θ<0 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 ofF• near the equilibrium point θ_≈Θ_e [53]. This may be constructed as follows:(62)θ˙tθ_≈Θ_e=Fθ_,rθ_≈Θ_e=FΘ_e,r+∂Fθ_,r∂θ_θ_≈Θ_eθ_-Θ_e+higherorderterms,where we need to ensure that the higher order terms (h.o.t.) are small enough for θ_ close to Θ_e, the equilibrium point, for all r. A commonly used criterion based on the L2-norm is(63)sup⁡Fh.o.t.θ_,rθ_θ_→Θ_e≈0_∀r,which 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Θ_e (54) we recall that(64)Fix=Fix0+∂Fi∂xjx=x0xj-x0+12∂2Fi∂xj∂xkx=x0xj-x0xk-x0+⋯,which for our system (50) provides us with the following higher order terms 1/2!θ_THθ_+⋯, where H represents the Hermitian matrices: (65)HF1=∂2F1∂θ12∂2F1∂θ1∂θ2∂2F1∂θ2∂θ1∂2F1∂θ22=0000,HF2=∂2F2∂θ12∂2F2∂θ1∂θ2∂2F2∂θ2∂θ1∂2F2∂θ22=12Nrr∗-12θ1+θa2+2Ncr∗-12000.We are interested in the equilibrium pointΘ_e=0_ (54) due to its physical significance. Thus we now consider θ_→Θ_e=0_, providing(66)HF2=12Nrr∗-12θa2+2Ncr∗-12000,which implies that 12Nrr∗-12θa2+2Ncr∗-12=r∗-1212Nrθa2+2Nc 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 r∗≠1 we need 12NrUa2+2Nc≈0 such that(67)Nr≈-Nc6θa2,which is unrealistic. The only manner in which (67) could hold is if both parameters were to approximate zero; that is, Nc,Nr≈0. 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 Nc and Nr. It is for this reason that the phase portraits presented are for Nc=0.01,0.1 and Nr=1.


### 6.4. Remarks

We find that forθ1=θ2=0 we are in fact referring to the point where θ=θ′=0, which is a point that would be located at the origin r=0 (see the boundary conditions). We find as either Nr→∞ or Nc→∞ that θ→0 at the point where θ′0=0. Thus as Nr and Nc 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 nearr=1. For the CN scheme we note that this behaviour occurs at lower values of Nc than Nr, indicating once again the more prominent impact the former parameter has on the behaviour of the scheme. In a similar fashion the CNNR scheme also produces negative temperature values for large values of Nr and Nc, particularly near r=1. Increasing θa we find that the value of Nr 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 θa with the same value now assigned to Nc 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 θa and Nr when employing the CNNR 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 parameterNr 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 Nc, which is the convective heat loss, would also cause a severe drop in the temperature. Thus, the instability of the equilibrium point 0,0 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 Nr and Nc are small, will the linear system approximate the nonlinear system.


## 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 parametersNr and Nc 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:(41)-αnθn+1m+1+1+δθnm+1-βnθn-1m+1=αnθn+1m+1-δθnm+βnθn-1m,where we define(42)αn=Δτγn4Δr+Δτ2Δr2r∗-12,βn=-Δτγn4Δr+Δτ2Δr2r∗-12,with γn and δ are defined as before (34) and (35), respectively. Upon substitution of θnm=ξmeinΔrω, where i2=-1, into (41) we obtain (43)-αnξm+1ein+1Δrω+1+δξm+1einΔrω-βnξm+1ein-1Δrω=αnξmein+1Δrω+1-δξmeinΔrω+βnξmein-1Δrω.Introducing Euler’s identities into (43) and simplifying provide (44)ξΔτγn4Δr-2isin⁡ωΔr+Δτ2Δr2r∗-12-2cos⁡ωΔr+1+δ=Δτγn4Δr2isin⁡ωΔr+Δτ2Δr2r∗-122cos⁡ωΔr+1-δ.Further simplification allows us to find the amplification factor(45)ξ=1-δ+Δτγn/2Δrisin⁡ωΔr+Δτ/Δr2r∗-12cos⁡ωΔr1+δ-Δτγn/2Δrisin⁡ωΔr-Δτ/Δr2r∗-12cos⁡ωΔr.We note at this point that the von Neumann stability analysis is applicable to problems with constant coefficients, and yetγn is dependent on rn. However, we may still deem the matrix to be a matrix of constant coefficients since we know the domain of rn (i.e., rn∈0,1). Hence, upon substitution of rn for n=0,1,…,N we obtain constant coefficients as required. We also note that 0<δ<1.Given thatr∗=rb/rt<1 for rn∈0,1 and having defined(46)γn=1rnr∗-1+1r∗-1,we find that since r∗-1<0 and rnr∗-1+1>0, γn<0. We also note that at r=0 we obtain γn=1/r∗-1<0 and when r=1 similarly γn=1/r∗r∗-1<0. Now we let γ¯n=-γn>0, providing the amplification factor(47)ξ=1-δ-Δτγ¯n/2Δrisin⁡ωΔx+Δτ/Δr2r∗-12cos⁡ωΔx1+δ+Δτγ¯n/2Δrisin⁡ωΔx-Δτ/Δr2r∗-12cos⁡ωΔx.We may conclude that sinceξ<1, 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 wherer=0; 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:(48)θ_˙=Fθ_,r,so that (49)θ˙1=θ2,(50)θ˙2=1-r∗θ2rr∗-1+1+Nrr∗-12θ1+θa4-θa4+Ncr∗-12θ12.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,(51)Fθ_,r=0_.Theorems for the stability of fixed points of the system(52)θ_˙=Aθ_+fθ_,rare presented by Coddington and Levinson [50], where A is the autonomous part of the Jacobian Dθ_Fθ_,r and a constant matrix (upon substitution of the relevant equilibrium points and parameter values) and the vector function f is continuous in θ_ and r. In order to employ these theorems we write our system of equations (50) in the form (52):(53)θ_˙=014Nrr∗-12θ1+θa3+2Ncr∗-12θ10θ_+01-r∗rr∗-1+1θ2.As per (51) we obtain the equilibrium points θ1,θ2:(54)0,0,(55)18NcNr2θa-10Nr3θa3+5.19615Nc3Nr3+14Nc2Nr4θa2-12NcNr5θa4+4Nr6θa633Nr-3NcNr+2Nr2θa23Nr18NcNr2θa-10Nr3θa3+5.19615Nc3Nr3+14Nc2Nr4θa2-12NcNr5θa4+4Nr6θa63-4θa3,0.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 A from (53) which we can write with general parameter values as follows: (56)λ1=-2r∗-1Ncθ1+2Nrθ13+6Nrθ12θa+6Nrθ1θa2+2Nrθa3,λ2=2r∗-1Ncθ1+2Nrθ13+6Nrθ12θa+6Nrθ1θa2+2Nrθa3,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 ofA at equilibrium point (54) are (57)λ1=-2Nrr∗-1θa3/2,λ2=2Nrr∗-1θa3/2.We can see that this equilibrium point (given the ranges for the relevant parameter values) represents a saddle point which is unstable; see Figure5(b).Figure 5
Phase plots indicating the trajectories around the equilibrium points (54) and (55) for θa=0.2, r∗=1/1.75, and Nr=1. (a) has Nc=0.1 whereas (b) has Nc=0.01. (a) is a closer view of the second equilibrium point.
(a)
(b)Coddington and Levinson [50] have shown that if at least one eigenvalue of A has a real part which is positive, then upon fθ_,r satisfying the requirement that, given any ϵ>0, there exist σ and R such that(58)fθ_,r≦ϵθ_θ_≦σ,r≧R,the solution θ_=0_ of (52) is not stable. In our case, we find that for the equilibrium point 0,0 at least one eigenvalue of A is positive and furthermore(59)fθ_,r=0001-r∗rr∗-1+1θ1θ2.We need to show that1-r∗/rr∗-1+1θ2≦ϵθ2. In an attempt to do so, we maximise the term on the left of the inequality; that is, we consider r=0 such that (60)1-r∗rr∗-1+1θ2≦1-r∗rr∗-1+1θ2≦1-r∗θ2=ϵθ2,where, given that 0<1-r∗<1, we define 0<ϵ<1. This proves our requirement (58), confirming that equilibrium point (54) is not stable.


### 6.2.2. Case  2: Equilibrium Point (55)

The eigenvalues ofA at this equilibrium point are lengthy for arbitrary values of Nc, Nr, r∗, and θa. However, as pointed out above, the second case needs to be for Ncθ1+2Nrθ13+6Nrθ12θa+6Nrθ1θa2+2Nrθa3<0. A variety of values were considered for the parameters Nc, Nr, r∗, and θa, all leading to eigenvalues of the form expected, which are (61)λ1=μi,λ2=-μi,wherei2=-1,providing us with a stable centre; see Figure 5(a). At this stage we note that this case requires θ1=θ<0 which is not physically relevant.


## 6.2.1. Case  1: Equilibrium Point (54)

The eigenvalues ofA at equilibrium point (54) are (57)λ1=-2Nrr∗-1θa3/2,λ2=2Nrr∗-1θa3/2.We can see that this equilibrium point (given the ranges for the relevant parameter values) represents a saddle point which is unstable; see Figure5(b).Figure 5
Phase plots indicating the trajectories around the equilibrium points (54) and (55) for θa=0.2, r∗=1/1.75, and Nr=1. (a) has Nc=0.1 whereas (b) has Nc=0.01. (a) is a closer view of the second equilibrium point.
(a)
(b)Coddington and Levinson [50] have shown that if at least one eigenvalue of A has a real part which is positive, then upon fθ_,r satisfying the requirement that, given any ϵ>0, there exist σ and R such that(58)fθ_,r≦ϵθ_θ_≦σ,r≧R,the solution θ_=0_ of (52) is not stable. In our case, we find that for the equilibrium point 0,0 at least one eigenvalue of A is positive and furthermore(59)fθ_,r=0001-r∗rr∗-1+1θ1θ2.We need to show that1-r∗/rr∗-1+1θ2≦ϵθ2. In an attempt to do so, we maximise the term on the left of the inequality; that is, we consider r=0 such that (60)1-r∗rr∗-1+1θ2≦1-r∗rr∗-1+1θ2≦1-r∗θ2=ϵθ2,where, given that 0<1-r∗<1, we define 0<ϵ<1. This proves our requirement (58), confirming that equilibrium point (54) is not stable.


## 6.2.2. Case  2: Equilibrium Point (55)

The eigenvalues ofA at this equilibrium point are lengthy for arbitrary values of Nc, Nr, r∗, and θa. However, as pointed out above, the second case needs to be for Ncθ1+2Nrθ13+6Nrθ12θa+6Nrθ1θa2+2Nrθa3<0. A variety of values were considered for the parameters Nc, Nr, r∗, and θa, all leading to eigenvalues of the form expected, which are (61)λ1=μi,λ2=-μi,wherei2=-1,providing us with a stable centre; see Figure 5(a). At this stage we note that this case requires θ1=θ<0 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 ofF• near the equilibrium point θ_≈Θ_e [53]. This may be constructed as follows:(62)θ˙tθ_≈Θ_e=Fθ_,rθ_≈Θ_e=FΘ_e,r+∂Fθ_,r∂θ_θ_≈Θ_eθ_-Θ_e+higherorderterms,where we need to ensure that the higher order terms (h.o.t.) are small enough for θ_ close to Θ_e, the equilibrium point, for all r. A commonly used criterion based on the L2-norm is(63)sup⁡Fh.o.t.θ_,rθ_θ_→Θ_e≈0_∀r,which 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Θ_e (54) we recall that(64)Fix=Fix0+∂Fi∂xjx=x0xj-x0+12∂2Fi∂xj∂xkx=x0xj-x0xk-x0+⋯,which for our system (50) provides us with the following higher order terms 1/2!θ_THθ_+⋯, where H represents the Hermitian matrices: (65)HF1=∂2F1∂θ12∂2F1∂θ1∂θ2∂2F1∂θ2∂θ1∂2F1∂θ22=0000,HF2=∂2F2∂θ12∂2F2∂θ1∂θ2∂2F2∂θ2∂θ1∂2F2∂θ22=12Nrr∗-12θ1+θa2+2Ncr∗-12000.We are interested in the equilibrium pointΘ_e=0_ (54) due to its physical significance. Thus we now consider θ_→Θ_e=0_, providing(66)HF2=12Nrr∗-12θa2+2Ncr∗-12000,which implies that 12Nrr∗-12θa2+2Ncr∗-12=r∗-1212Nrθa2+2Nc 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 r∗≠1 we need 12NrUa2+2Nc≈0 such that(67)Nr≈-Nc6θa2,which is unrealistic. The only manner in which (67) could hold is if both parameters were to approximate zero; that is, Nc,Nr≈0. 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 Nc and Nr. It is for this reason that the phase portraits presented are for Nc=0.01,0.1 and Nr=1.


## 6.4. Remarks

We find that forθ1=θ2=0 we are in fact referring to the point where θ=θ′=0, which is a point that would be located at the origin r=0 (see the boundary conditions). We find as either Nr→∞ or Nc→∞ that θ→0 at the point where θ′0=0. Thus as Nr and Nc 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 nearr=1. For the CN scheme we note that this behaviour occurs at lower values of Nc than Nr, indicating once again the more prominent impact the former parameter has on the behaviour of the scheme. In a similar fashion the CNNR scheme also produces negative temperature values for large values of Nr and Nc, particularly near r=1. Increasing θa we find that the value of Nr 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 θa with the same value now assigned to Nc 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 θa and Nr when employing the CNNR 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 parameterNr 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 Nc, which is the convective heat loss, would also cause a severe drop in the temperature. Thus, the instability of the equilibrium point 0,0 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 Nr and Nc 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:(68)θ=TTb,θa=TaTb,r^=rrb,r^∗=rtrb,Nc=ρfCp,fKgβfrb2TbVfkeffs,Nr=2ϵσFf-arb2Tb3keffs.Given that we are considering the time dependent case we also incorporate(69)τ=kefftρCvrb2,providing the following equation for consideration:(70)∂θ∂τ=1r∂∂rr∂θ∂r-Ncθ-θa2-Nrθ4-θa4,where we have once again dropped the hat on r for convenience, r∗>1 is defined as rt/rb, and the parameters Nc and Nr are stated as per (68).Employing the nondimensional transformations given by (68) we obtain the following boundary conditions:(71)dθdrr=r∗=0,θ1,τ=1,such that r∈1,r∗. We use the step sizes Δr=0.01 and Δτ=0.00001 due to the domain of r. The numerical solutions obtained mirror those provided above. The only difference is the behaviour of the solutions for θa; 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:(72)θ˙1=θ2,θ˙2=-θ2r+Ncθ1-θa2+Nrθ14-θa4.It must be noted that we consider the case wherer=r∗; 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]:(73)θ_˙=Aθ_+fθ_,r,so that we may employ theorems for the stability of fixed points of the system by considering the eigenvalues of A. We obtain the equilibrium points θ1,θ2:(74)θa,0,(75)5.19615Nc3Nr3+14.Nc2Nr4θa2-12.NcNr5θa4+4.Nr6θa6+18.NcNr2θa-10.Nr3θa32/3-3.NcNr+2.Nr2θa23Nr5.19615Nc3Nr3+14.Nc2Nr4θa2-12.NcNr5θa4+4.Nr6θa6+18.NcNr2θa-10.Nr3θa31/3-θa3,0.In considering the eigenvalues ofA from (73), we can write them with general parameter values as follows: (76)λ1=-2Ncθ1+2Nrθ13-Ncθa,λ2=2Ncθ1+2Nrθ13-Ncθa,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 ofA at equilibrium point (74) are (77)λ1=-2Nrθa3/2,λ2=2Nrθa3/2,which represent an unstable saddle point; see Figures 6(b) and 7(b) corresponding to the saddle point obtained for the equilibrium point 0,0 for the previous nondimensionalisation.Figure 6
Phase plots indicating the trajectories around the equilibrium points (75) for θa=0.6, r∗=1.75, Nc=0.001, and Nr=0.01. (a) is a closer view of the second equilibrium point for the alternate nondimensionalisation.
(a)
(b)Figure 7
Phase plots indicating the trajectories around the equilibrium points (75) for θa=0.6, r∗=1.75, Nc=100, and Nr=100. (a) is a closer view of the second equilibrium point for the alternate nondimensionalisation.
(a)
(b)This comparison to the equilibrium point0,0 obtained via the previous nondimensionalisation is important given the numerical solutions we obtain for certain values of θa. 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 θa, we were unable to do so in this instance. For θa→1 we found that we were unable to obtain numerical solutions for certain values of Nc and Nr. 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: (78)HF1=∂2F1∂θ12∂2F1∂θ1∂θ2∂2F1∂θ2∂θ1∂2F1∂θ22=0000,HF2=∂2F2∂θ12∂2F2∂θ1∂θ2∂2F2∂θ2∂θ1∂2F2∂θ22=12Nrθ12+2Nc000.Now asθ_→θa,0, we have(79)HF2=12Nrθa2+2Nc000,which means that for the linear system to approximate the nonlinear dynamics of the system we now require Nc≈0 and either Nr or θa 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 CNNR schemes.In considering the CN scheme, we observe that for the caseNc⪅5 and θa<1 we were able to obtain solutions as Nr→∞. Once θa→1 we found that there were values for Nc and Nr (not deemed large in magnitude) for which we could not obtain solutions. In particular we noticed that as Nr→∞ with θa→1 the temperature profiles oscillated, often attaining negative temperature values. The moment we diminished the value of the ambient temperature θa we were able to obtain realistic solutions once more. Furthermore, for reasonably small values of Nr and θa we found that for values of Nc which could not be deemed excessively high (Nc⪆40) we were unable to obtain solutions to our problem. As Nc→∞ we needed Nr and θa to decrease in value in a comparative fashion: for Nc=100 with Nr=0.5 we could not obtain a solution to θa=0.8 but we could do so for θa=0.1. As such, the impact of the convection parameter Nc on the ability of the CN scheme to obtain coherent solutions is once again confirmed.With regard to theCNNR scheme we find that as Nc→∞ for small values of Nr and θa the scheme is unable to obtain solutions. Furthermore, as Nr→∞ and θa→1 we are unable to obtain solutions for small values of Nc; in particular we found oscillatory behaviour at r=1. 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 Nc, Nr, and θa 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θa for which limitations are placed upon Nc, Nr 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 ofA at this equilibrium point are lengthy for arbitrary values of Nc, Nr, r∗, and θa. However, as before we realise that the eigenvalues would have the form (80)λ1=μi,λ2=-μi,wherei2=-1,providing us with a stable centre; see Figures 6(a) and 7(a). For small values of Nc and Nr 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 θ1=θ>0 we indeed have a physically relevant equilibrium point. We can in this fashion confirm that the boundary dynamics are at the fin tip (where θ2=θ′=0), while the temperature is such that θ1=θ<1; 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 1 (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 θa provide stable dynamics (Case 2: 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 θ2=θ′ (Case 1: saddle point).


## 7.1. Case  1: Equilibrium Point (74)

The eigenvalues ofA at equilibrium point (74) are (77)λ1=-2Nrθa3/2,λ2=2Nrθa3/2,which represent an unstable saddle point; see Figures 6(b) and 7(b) corresponding to the saddle point obtained for the equilibrium point 0,0 for the previous nondimensionalisation.Figure 6
Phase plots indicating the trajectories around the equilibrium points (75) for θa=0.6, r∗=1.75, Nc=0.001, and Nr=0.01. (a) is a closer view of the second equilibrium point for the alternate nondimensionalisation.
(a)
(b)Figure 7
Phase plots indicating the trajectories around the equilibrium points (75) for θa=0.6, r∗=1.75, Nc=100, and Nr=100. (a) is a closer view of the second equilibrium point for the alternate nondimensionalisation.
(a)
(b)This comparison to the equilibrium point0,0 obtained via the previous nondimensionalisation is important given the numerical solutions we obtain for certain values of θa. 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 θa, we were unable to do so in this instance. For θa→1 we found that we were unable to obtain numerical solutions for certain values of Nc and Nr. 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: (78)HF1=∂2F1∂θ12∂2F1∂θ1∂θ2∂2F1∂θ2∂θ1∂2F1∂θ22=0000,HF2=∂2F2∂θ12∂2F2∂θ1∂θ2∂2F2∂θ2∂θ1∂2F2∂θ22=12Nrθ12+2Nc000.Now asθ_→θa,0, we have(79)HF2=12Nrθa2+2Nc000,which means that for the linear system to approximate the nonlinear dynamics of the system we now require Nc≈0 and either Nr or θa 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 CNNR schemes.In considering the CN scheme, we observe that for the caseNc⪅5 and θa<1 we were able to obtain solutions as Nr→∞. Once θa→1 we found that there were values for Nc and Nr (not deemed large in magnitude) for which we could not obtain solutions. In particular we noticed that as Nr→∞ with θa→1 the temperature profiles oscillated, often attaining negative temperature values. The moment we diminished the value of the ambient temperature θa we were able to obtain realistic solutions once more. Furthermore, for reasonably small values of Nr and θa we found that for values of Nc which could not be deemed excessively high (Nc⪆40) we were unable to obtain solutions to our problem. As Nc→∞ we needed Nr and θa to decrease in value in a comparative fashion: for Nc=100 with Nr=0.5 we could not obtain a solution to θa=0.8 but we could do so for θa=0.1. As such, the impact of the convection parameter Nc on the ability of the CN scheme to obtain coherent solutions is once again confirmed.With regard to theCNNR scheme we find that as Nc→∞ for small values of Nr and θa the scheme is unable to obtain solutions. Furthermore, as Nr→∞ and θa→1 we are unable to obtain solutions for small values of Nc; in particular we found oscillatory behaviour at r=1. 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 Nc, Nr, and θa 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θa for which limitations are placed upon Nc, Nr 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 ofA at this equilibrium point are lengthy for arbitrary values of Nc, Nr, r∗, and θa. However, as before we realise that the eigenvalues would have the form (80)λ1=μi,λ2=-μi,wherei2=-1,providing us with a stable centre; see Figures 6(a) and 7(a). For small values of Nc and Nr 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 θ1=θ>0 we indeed have a physically relevant equilibrium point. We can in this fashion confirm that the boundary dynamics are at the fin tip (where θ2=θ′=0), while the temperature is such that θ1=θ<1; 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 1 (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 θa provide stable dynamics (Case 2: 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 θ2=θ′ (Case 1: 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 ifNr and Nc 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 ofNr and Nc, showing that as they tend to infinity erratic and unstable behaviour may be observed at the equilibrium point which represents the fin tip (θ→0, θ′0=0). 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 Nr and the ambient temperature impact more severely on the stability of the CNNR 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 θ>0 and θ′0=0 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 θa,0 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 ofNr and Nc we find that θ1=θ>0, 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 0<θ<1 and θ′=0, are indeed stable. Given the other equilibrium point, it would seem however that for cases where θ=θa we would revert to unstable behaviour at the fin tip. What is interesting to note is that since in the latter nondimensionalised equation θa 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.



---

*Source: 1658305-2017-06-27.xml*