## Qualitative and Quantitative Techniques for Differential Equations Arising in Mathematical Physics

View this Special IssueResearch Article | Open Access

R. Jooma, C. Harley, "Heat Transfer in a Porous Radial Fin: Analysis of Numerically Obtained Solutions", *Advances in Mathematical Physics*, vol. 2017, Article ID 1658305, 20 pages, 2017. https://doi.org/10.1155/2017/1658305

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

**Academic Editor:**Igor L. Freire

#### Abstract

A time dependent nonlinear partial differential equation modelling heat transfer in a porous radial fin is considered. The Differential Transformation Method is employed in order to account for the steady state case. These solutions are then used as a means of assessing the validity of the numerical solutions obtained via the Crank-Nicolson finite difference method. In order to engage in the stability of this scheme we conduct a stability and dynamical systems analysis. These provide us with an assessment of the impact of the nonlinear sink terms on the stability of the numerical scheme employed and on the dynamics of the solutions.

#### 1. Introduction

Circular annular fins are extensively used to increase the rate of heat transfer from a heat source for a given temperature difference in heat exchange devices or to reduce the temperature difference between the heat source and the given heat flow rate of heat sink [1]. The use of fins is widely ranging in engineering applications where it is necessary to enhance the heat transfer from a surface to an adjacent coolant so that the fin can perform within the acceptable temperature limits. These engineering applications range from considerably large systems such as industrial heat exchangers to smaller systems such as transistors. Radial fins have been conventionally used as a coolant for internal combustion engines, heat exchanges, compressors, and so forth. Due to these extensive practical applications this is considered a vibrant field of research; this is true even more so since the problems that arise are nonlinear and hence not always solvable via analytical techniques.

Many research articles have investigated the use of porous fins [2]. Even though a porous material fin has low thermal conductivity, a vast area of the material comes into contact with the cooling agent enabling the porous fin to give superior performance [2]. Over the past decades, numerous studies have been conducted on the performance of annular fins [3–7]. Aziz and Rahman [8] examined a fin comprising functionally graded material and analysed the performance on the radial fin with a continuously increasing thermal conductivity in the radial direction. They discovered that the heat transfer as well as the fin efficiency and effectiveness are at their highest maximum values when the thermal conductivity of the fin varies inversely with the square of the radius. Furthermore, they found that the use of a spatially averaged thermal conductivity model is not recommended due to large errors occurring in some cases. Kiwan [9] conducted a thermal analysis of natural convection porous fins by introducing Darcy’s model to construct the energy equation governing the distribution of temperature. He further discovered that, by choosing a precise value for the thermal conductivity ratio and the fin length to thickness ratio, the performance of the porous fin exceeded the performance of the solid fin. A study was conducted by Kiwan and Zeitoun [10] to test the performance of rectangular porous fins mounted around the inner cylinder of a cylindrical annulus by performing a finite volume type numerical study. It was concluded that, in comparison to solid fins, porous fins provided higher transfer rates for similar configurations and that the heat transfer rate from the cylinder equipped with porous fins decreased as the fin inclination increased. Gorla and Bakier [11] investigated natural convection and radiation in porous fins. They found that the radiation transfers more heat in comparison to a similar model without radiation. Abu-Hijleh [12] analysed the effects of using permeable fins on the forced convection heat transfer from a horizontal cylinder. The results obtained were similar to results obtained as per Kiwan and Zeitoun [10] in terms of the permeable fins providing much higher heat transfer rates. To study radial fins, a combination of the Taylor transformation and finite difference approximation was implemented by Yu and Chen [13, 14]. They further performed a study on the optimization of a circular fin with variable thermal parameters. Naidu et al. [15] set forth a numerical study of natural convection from a cylindrical fin placed in a cylindrical porous enclosure. Hence, they conducted a conjugate conduction-convection analysis by solving the heat conduction equation. Moitsheki and Harley [16] studied the transient heat transfer through a longitudinal fin of various profiles by employing classical Lie point symmetry methods. They observed that for long periods of time the temperature profile becomes unusual for the heat transfer in longitudinal triangular and concave parabolic fins. This, however, was corrected by increasing the thermogeometric fin parameter. In recent studies, Darvishi et al. [17] accounted for the effects of radiation and convection heat transfer in a rectangular radial porous fin. This allowed for the heat flow to infiltrate the porous fin enabling a solid-fluid interaction to occur. They concluded that in a model containing radiation more heat is present than in a similar model without radiation. In a similar context, our model takes into account the time rate of change of internal energy and the heat flow due to conduction as well as the heat due to radiation and convection.

Extensive analytical studies have been done via the Differential Transformation Method (DTM) for the solution of problems such as the one under discussion. The DTM, which was first proposed by Zhou [18], is a seminumerical-analytical method applied to linear and nonlinear systems of ordinary differential equations. The method captures the exact solution in terms of a Taylor series expansion. This method has been successfully implemented in engineering applications [19–25]. Ndlovu and Moitsheki [26] derived approximate analytical solutions for the temperature distribution in a longitudinal rectangular and convex parabolic fin with temperature dependent thermal conductivity and heat transfer coefficients. These authors, for the first time, used a two-dimensional DTM for the transient heat conduction problem. Ertürk [27] constructed seminumerical-analytical solutions for a linear sixth-order boundary value problem using the DTM. It was observed that the method served as an effective and reliable tool for such problems. Recently, Torabi et al. [28] analysed the radiative radial fin with temperature dependent thermal conductivity by implementing the DTM as well as the Boubaker polynomials expansion scheme (BPES). Similar to the results obtained by Ertürk [27], suitable results were obtained in predicting the solution for both BPES and DTM. A study of a radial fin in terms of the fin’s thickness with convection heating at the base and the convection-radiative cooling at the tip was conducted by Aziz et al. [29]. Furthermore, they conducted an analysis using DTM and verified the results by comparing it to an exact analytical solution. The preceding literature clearly shows that the Differential Transformation Method has been applied to problems relating to many different fins, but no attempt has been made to apply it when investigating the heat transfer in a porous radial fin.

In terms of numerical investigations, an efficient, accurate, extensively validated, and unconditionally stable method was developed in the mid-20th century by Crank and Nicolson [30] in order to evaluate numerical solutions for nonlinear partial differential equations. Rani et al. [31] obtained a solution for the time dependent nonlinear coupled governing equations with the help of an unconditionally stable Crank-Nicolson scheme for the transient couple stress fluid flowing over a vertical cylinder. They observed that the time taken for the flow to reach steady state increases as the Schmidt and Prandtl values increase and decreases with respect to the buoyancy ratio parameter. Furthermore, Ahmed et al. [32] employed the Crank-Nicolson finite difference scheme to the conservative equations in modelling porous media transport for magnetohydrodynamic unsteady flow. They found that the flow velocity and temperature decrease with an increase in the Darcian drag force. The concept of the Crank-Nicolson scheme combined with the Newton-Raphson method was used by Qin et al. [33] to model the heat flux and to estimate the evaporation in applied hydrology and meteorology. The Crank-Nicolson method was used to expand the differential equations whereas the iterative Newton-Raphson method was used to approximate latent heat flux and surface temperatures. Both these methods proved to be successful. In a similar context, Janssen et al. [34] implemented the Crank-Nicolson scheme to transform a system of differential equations into algebraic equations. The Newton-Raphson method is used to implicitly enhance the model’s efficiency by improving the poor convergence rate. Once again, it can be seen that the Crank-Nicolson scheme with the Newton-Raphson method has not been implemented for heat transfer in a porous radial fin.

As far as we know, there has been no or very little work that has been done on obtaining asymptotic solutions or employing a dynamical systems analysis to the problem presented in this research. The purpose of the asymptotic solution is to reveal the dominant physical mechanisms of the model. It can be seen in Moitsheki and Harley [16, 35] that the impact of the thermogeometric parameter () in terms of its proportionality to the length of the fin () was observed. They found that the heat transfer in the fin seemed to be unstable for small values of () due to the fact that . By investigating the asymptotic solution to the steady state heat transfer in a rectangular longitudinal fin, they were able to validate the above relationship and establish the importance of the fin length. Furthermore, the same authors [16, 35] conducted a small scale dynamical analysis. In order to expand the analysis in [35], Harley [36] employed an in-depth dynamical analysis to monitor the behaviour especially at the fin tip. This dynamical analysis also served as a means of investigating the role and effect of the thermogeometric parameter. In this work we do not derive an asymptotic solution; such a solution was derived for the problem but we did not deem it useful in terms of providing deeper insight into the dynamics observed. We will however employ a dynamical systems analysis, which we find to be of immense use in classifying the dynamics of specific points of the solution and the engagement between parameters.

A vast amount of work has been conducted on the steady state case of problems in this field since analytical methods lend themselves more easily to the solution of these equations. In this research, a time dependent partial differential equation modelling the heat transfer in a porous radial fin will be considered. The equation will be derived and nondimensionalised appropriately in Section 2. As a means of comparison to the computational methods employed we structure a semianalytical solution via the Differential Transformation Method in Section 3. In the sections to follow, Sections 4 and 5, the partial differential equation is solved numerically using the Crank-Nicolson scheme combined with the Newton-Raphson method as a predictor-corrector. In order to assess the effectiveness of this scheme and its limitations we employ a dynamical systems analysis in Section 6; in this manner we are able to engage with the limitations placed on parameter values with regard to obtaining solutions. Section 7 provides insight into the comparative dynamics and stability of the equation obtained via an alternate means of nondimensionalisation. Concluding remarks are made in Section 8.

The importance of this work relates to the detailed analysis of the dynamics of the temperature at the fin tip. Furthermore, we are able to investigate the impact of the nonlinear source terms on the stability of the schemes employed and the solutions that can be obtained. This highlights the care than needs to be taken when choosing to solve equations of this nature, particularly when solving via numerical tools. Parameters of physical importance are shown to have value limitations due to stability requirements. Consequently, while numerical solutions may have been obtained here and by other authors noting that these do not indicate an ability to obtain solutions for all relevant parameter values or that it may be assumed that solutions that seemed to have converged are indeed dynamically stable and physically accurate is needed.

#### 2. Model Derivation

Consider a cylindrical porous radial fin with base radius , tip radius , and thickness as shown in Figure 1. The fin comprises an effective thermal conductivity porous material and permeability . It is assumed that the tip of the fin is adiabatic (i.e., a process that occurs without the transfer of heat/matter between the system and its surroundings) and the base of the fin is maintained at a constant temperature . The internal energy per unit volume with the absolute temperature is denoted by . In accordance with Darcy’s law, the fin makes contact with an ambient fluid which infiltrates the fin. The ambient fluid comprises an effective density of the porous fin , a specific heat of the porous fin , the kinematic viscosity of the ambient fluid , the thermal conductivity of the ambient fluid , and the volumetric thermal expansion coefficient of the ambient fluid . The top and bottom surfaces are presumed to have a constant surface emissivity and emit radiation to the ambient fluid at temperature . This also serves as the radiation heat sink.

Constructing the energy balance of a fin element (Figure 1) of thickness , circumference , and radial width at position , we obtainwherewhere (2) is the time rate of change of internal energy in the volume element, (3) is a basis of the application of Fourier’s law of heat conduction, (4) denotes the heat loss due to radiation through the top and bottom surfaces of the fin, and (5) denotes the heat loss due to convection due to Darcy’s flow through a porous fin. At any radiation location , the velocity of the buoyancy driven flow is obtained by the implementation of Darcy’s as follows:

The Darcy model stimulates the fluid-solid interaction in the porous medium. By substituting (2)–(6) into (1), the following nonlinear partial differential equation governing the temperature distribution in the fin is obtained:

The boundary conditions are given by

The fin is initially kept at the temperature of the fluid (ambient temperature) which is given by

We introduce the following nondimensional quantities:

Upon appropriate substitution and rearrangement (8) becomeswhere is defined as and the parameters and are stated as per (12).

Employing the nondimensional transformations given by (12) we find that at we obtain and at we obtain . Thus, we now work on the interval with . As such the nondimensional boundary conditions are given bywith the initial condition

Equation (13) is a nonlinear partial differential equation comprising two nonlinear terms. The first nonlinear term refers to the buoyancy or natural convection transport of energy by the infiltrate whereas the second nonlinear term refers to the surface radiative heat transfer from the fin to the ambient fluid. The parameter is a combination of Rayleigh number , Darcy’s number , the thermal conductivity ratio , and the ratio of the fin tip radius to fin thickness. The parameter indicates the role of surface radiation to conduction in the fin. The parameter is the ratio of the ambient fluid temperature and the difference between the base temperature and the fluid temperature .

For the rest of this research we drop the hat on the independent variable for simplicity. Furthermore, the model employed closely follows the work by Darvishi et al. [2]; in order to be able to compare the results obtained by those authors we use similar parameter values for , , , and . The value of is however defined as the inverse of in the investigations of Darvishi et al. [2] (in their work they define this parameter as ). To allow for a means of comparison we prescribe values for as .

*Remark 1. *The nondimensional ambient temperature is defined as in the work of Darvishi et al. [2]. Their nondimensionalisation has also been implemented in this research (see Section 7); however a dynamical systems analysis led to the conclusion that there are limitations on certain parameters when obtaining solutions. As such, we altered our nondimensionalisation to define as , as above, as a means of comparison. This change in ratio impacts on the behaviour observed by the solution with regard to changes in the value of . Instead of a linear relationship between and we have a nonlinear one given by , such that the behaviour inversely corresponds to that of Darvishi et al. [2]. As a means of checking the method and analysis conducted in this research, we compared solutions and the dynamics of these solutions to those obtained via the alternate nondimensionalisation employed in [2]; in Section 7 we will briefly comment on the results.

#### 3. Differential Transformation Method

The DTM is an approximation method based on the Taylor series expansion that can easily be applied to several linear and nonlinear problems. It constructs an analytical solution in terms of a polynomial and is capable of reducing the size of the computational work required for solution of the equation. The predominant advantage of this method is that it can be applied directly to linear and nonlinear ordinary differential equations without requiring linearisation, discretisation, or perturbation [37]. This method has been used in a similar manner by Ghasemi et al. [38] to solve the nonlinear temperature distribution equation with temperature dependent internal heat generation in porous and solid longitudinal fins. Furthermore, Fidanoglu et al. [39] used the method to construct a general expression for the heat distribution profile of a fin with spine geometry. Based on the fin effectiveness, entropy values, and efficiency, it was found that the cylindrical fin had the highest rate of heat transfer.

An arbitrary analytical function in the domain can be expanded via the Taylor series about a point as

The differential transform of can thus be defined aswith the inverse differential transform as

The fundamental theorems of the one-dimensional DTM, some of which will be implemented in this paper, are given as follows.

Theorem 2. *If , then *

Theorem 3. *If , then , where is a constant.*

Theorem 4. *If , then *

Theorem 5. *If , then *

Theorem 6. *If , then , where *

Theorem 7. *If , then *

Through the use of these properties of the one-dimensional DTM, we obtain the following iterative procedure for the steady state case of (13):

This method requires initial conditions; that is, we require the values of and since the equation is of second order. We notice that even though we have an initial value for the first derivative, we do not have a value for at . Therefore we specify our initial conditions for iterative procedure equation (4) as follows:where the initial condition for is given an assumed value while the first derivative condition becomes . Through the use of (21) and the conditions given by (22) we are able to obtain which allows us upon substitution into (18) to obtainwhere . However, in order to determine the solution given by (24) we need to first find the value of . In order to determine a value for , we could use the false position method. The boundary condition which has not yet been employed is . Thus we substitute the values obtained as per (23) into (24) for . This provides us with the equality , where can be written as and hence we need to solve the equation . By plotting and we are able to obtain an interval , where and are two guesses to , within which the two functions are equal. Using these values we are able to obtain an improved value which satisfies as follows:

This procedure iterates until a chosen error tolerance is reached such that . We employed this method as well as NSolve in Mathematica as a means of verifying the accuracy of the value. As a brief example, we consider a case where , , , , and to obtain resulting in

In this work however, we have employed a polynomial of order as our steady state solution; see Section 5.2 for discussions on solutions obtained. We are able to obtain a polynomial of degree via this methodology; however for any degree larger we find that limited memory becomes a stumbling block for Mathematica. Furthermore, for we find that our tip temperature whereas for , ; thus increasing the degree of the polynomial does not add much in terms of accuracy. This discrepancy is also not of extreme concern since this semianalytical solution is not our only means of comparison for the numerical solutions obtained below (we also consider the work by Darvishi et al. [2]). We can remark here that many terms are likely to be needed in order for the series solution to match solutions obtained in [2]; however it may still be used as a reasonable measure for the numerical solutions we obtain.

#### 4. Numerical Investigation

##### 4.1. Crank-Nicolson Scheme

The Crank-Nicolson method [40] is a finite difference method which is implicit in time, often provides unconditional stability, and has a high order of accuracy. Given the method’s advantages we will implement the scheme for the equation under consideration. The complexity of our problem lies in the nonlinear sink terms which are likely to influence not only the stability of the scheme, but also the order of accuracy it is able to obtain. We will at first employ the scheme with an explicit consideration of the nonlinear terms. In the subsection to follow we will develop the method further as a means of comparison.

In terms of the finite difference schemes, we approximate using the forward time difference approximation at which giveswhere is the step size in the time domain . Furthermore we approximate the spatial terms and using the central difference approximations at givingrespectively, where is the step size in the spatial domain .

The Crank-Nicolson scheme makes use of approximations to values halfway between nodes by taking averages over two time points. Thus the finite difference approximations on the spatial terms are given by

Given our nondimensional boundary conditions (14)-(15) and the finite difference approximation given by (27) we find that

The initial condition is given as

Having employed the Crank-Nicolson scheme on spatial terms and a forward difference approximation for the time derivative on (13), we have

It should be noticed that we have not applied the Crank-Nicolson method to the nonlinear terms. This is for simplicity at present; however this will be revisited in the next subsection. From (32) we thus obtainwhere we have definedsuch thatfor the sake of simplicity. This equation can be given in the following matrix notation:where and are the relevant coefficient matrices and is defined as the vector of nonlinear terms only.

##### 4.2. Crank-Nicolson Scheme: Newton-Raphson Implementation

In numerical computation, a predictor-corrector method is an algorithm that implements two steps. First, the prediction step calculates a rough approximation of the desired quantity (in this case the nonlinear term on the th time level) whereas the second step refines the initial approximation using other means. In this research we will require such a technique, known as the Newton-Raphson method, given the presence of the nonlinear sink terms. The Crank-Nicolson method incorporating the Newton-Raphson method has been used to determine the solution of a nonlinear diffusion equation (see Kouhia [41] and Habibi et al. [42]) to model the electron heat transfer process in laser inertial fusion for the energy transport mechanism from the region of energy decomposition into the ablation surface. The method proved to be efficient as per the numerical results obtained and hence is a natural extension of the method employed in the previous subsection. We implement the Newton-Raphson method [43] as done by Britz et al. [44] and discussed in [45]. This method is iterative and has been known for its fast convergence to the root as well as its dependence on the initial condition. Furthermore, in order to find the explicit inverse of a general tridiagonal matrix, we shall implement Usmani’s method [46]; the reasoning behind this decision will be discussed at a later stage.

In this subsection we will now employ the Crank-Nicolson method for the nonlinear sink terms as well. As such we rewrite our scheme (33) as follows:where we again use the definitions given by (34)-(35). In matrix notation,where, once again, and are the relevant coefficient matrices and and are defined as the vectors of nonlinear terms only.

We obtain an approximation for the nonlinear vector at the th time level in an iterative fashion as done by Britz et al. [44]. The systemis solved where is the set of difference equations created as per the right hand side of (37) and is the Jacobian of the left hand side of (37). The starting vector at is chosen as per initial condition (15) such that . We solve for using (39) and then iterate as follows:until convergence is reached which provides us with an approximation to . The Newton-Raphson iteration converges within 2-3 steps given that changes are usually relatively small.

The convergence rate of the Crank-Nicolson method incorporating the Newton-Raphson method was found to be extremely slow, making the method infeasible. We discovered that this was due to the fact that the Jacobian** J** is either close to being singular or badly scaled. As such, we decided to implement Usmani’s method, which is a simple algorithm for finding the explicit inverse of a general Jacobi tridiagonal matrix [46]. This implementation provided the expected results in terms of obtaining the inverse of the Jacobian. Nonetheless, the computational time taken by the method to reach convergence proved inefficient.

#### 5. Remarks

##### 5.1. Comparison of Numerical Schemes

In this subsection we compare the results we obtained via our two numerical methods, that is, the Crank-Nicolson scheme (CN) and the Crank-Nicolson scheme with Newton-Raphson predictor-corrector and Usmani method (). It is important to note that only the first scheme matches the results and behaviour obtained and observed in Darvishi et al. [2]. While the second scheme provides appropriate behaviour, the values at the origin indicate clear discrepancies with the work conducted in [2] even for large values of , where . Furthermore, the length of computational time required for near convergence exceeds that required when the CN scheme is implemented which converges swiftly. As such, aside from the discrepancies in terms of the solutions obtained, we conclude that the method is not efficient in terms of running time.

In order to make a comparison between the two numerical schemes we have employed, we allow our numerical schemes to iterate for large time, termed , so that we allow for convergence to a steady state. In this instance we find that convergence has occurred for the Crank-Nicolson scheme at ; however this is not the case for the Crank-Nicolson method with the predictor-corrector. To maintain our CFL number as less than one for stability purposes, we choose and . As can be surmised, the consequence of this is that the Crank-Nicolson scheme which employs the Newton-Raphson method as a means of updating the nonlinear terms requires a lot of computational time.

We find that while the two numerical schemes produce solutions similar in shape and behaviour, there are discrepancies when we consider the errors between the methods. The temperature at the origin, that is, , is not the same across the two methods; see Tables 1–4. From a consideration of the temperature at the fin tip we find that the solutions do not converge to the same value of for set parameter values and furthermore do not converge at the same rate. At the difference in the values of across the schemes and the maximum error (obtained via comparisons between the numerical schemes) diminish, however not to the extent that the difference is numerically inconsequential. The convergence rate we suspect is heavily influenced by the fact that the Jacobian required for the predictor-corrector method is close to singular or badly scaled. In Figures 2 to 4 we employ the Crank-Nicolson method without the predictor-corrector method due to the speed with which it is able to reach convergence and its accuracy upon comparison with the work by Darvishi et al. [2].

| ||||||||||||||||||||||||||||||||

Parameter values: , , , and . |

| ||||||||||||||||||||||||||||

Parameter values: , , , and . |

| ||||||||||||||||||||||||||||

Parameter values: , , , and . |

| ||||||||||||||||||||||||||||||||

Parameter values: , , , and . |

While there are discrepancies we can see behaviour patterns emerging from the different solutions. In Figure 2 as well as Table 1 we observe that an increase in the dimensionless ambient temperature leads to a decrease in the tip temperature; however this can be misleading. We note, as before, that the dimensionless ambient temperature is defined as in the work of Darvishi et al. [2]. In this work we altered our nondimensionalisation to define as , as above. This change in ratio impacts on the behaviour observed by the solution with regard to changes in the value of . Instead of a linear relationship between and we have a nonlinear one given by , such that the behaviour inversely corresponds to that of Darvishi et al. [2]. Hence, our results are consistent with the physical dynamics of the problem.

A consideration of Figure 3 and Table 2 shows the effect of the natural convection heat loss () on the temperature distribution in the fin for set parameter values. As the buoyancy effect increases, that is, , the local temperature in the fin decreases. This is indicative of the effect of natural convective heat loss on the temperature distribution in the fin when the radiation heat loss, the ambient temperature, and the ratio of inner to outer radius are kept fixed. Since buoyancy is principally a macroscale effect we see that the buoyancy force () influences the velocity and temperature fields. Similarly, for Figure 4 and Table 3 the increasing surface radiative heat transfer from the fin to the ambient fluid, that is, , leads to a decrease in the local fin temperature. Furthermore, the impact of radiation () influences the temperature field by increasing the heat transfer rates from the surface. We find that when the convection parameter () and the radiation parameter () are allowed to vary, the local fin temperature decreases because of the increasing strength of radiative heat exchange between the exposed surface of the fin and the ambient [2].

##### 5.2. Comparison of Numerical Results to Semianalytical Solution

We will now compare the results we obtained via our two numerical methods against the series solution obtained via our steady state DTM. We calculate a polynomial of degree via the DTM and also consider the work in [2] to be a benchmark. We find that all three solutions perform well; however only the Crank-Nicolson scheme performs as accurately and efficiently as required.

It can be seen in Table 5, for increasing , that both the value of for the DTM and the maximum error produced between the DTM and the CN scheme increase, whereas the error between the DTM and the scheme decreases. Table 6 indicates that the maximum difference in the temperature between the DTM and the CN scheme with and without the Newton-Raphson method, respectively, agrees to the first decimal place only for ; the CN method performs well across the various values of this parameter. Furthermore, this table shows that for increasing values of the error between the DTM and the CN method decreases when ; in turn Table 7 indicates that for increasing values of at we obtain an increase in the error. These relationships are inverted when considering the comparison between the semianalytical solution and the scheme: the error increases for increasing values of at and decreases for increasing values of at . These occurrences could be indicative of the impact of the parameters and on the CN and methods, respectively. For the CN scheme an increase in the values of also leads to an increase in the error, while for the scheme the inverse is true. We thus observe that an increase in the values of and (which are both contained in only one of the sink terms) leads to an increase in the error between the DTM and the CN scheme. In turn, an increase in the value of , which is the coefficient of the second sink term, leads to an increase in the error between the DTM and the method.

| ||||||||||||||||||||||||||||||||

Parameter values: , , , and . |

| ||||||||||||||||||||||||||||

Parameter values: , , , and . |

| ||||||||||||||||||||||||||||

Parameter values: , , , and . |

The impact of a change in the value of provides behaviour which is quite stable; see Table 8. The error decreases in both instances; however the error for the CN scheme is much less than that obtained with the introduction of the Newton-Raphson predictor-corrector.

| ||||||||||||||||||||||||||||||||

Parameter values: , , , and . |

The solutions obtained in Darvishi et al. [2] indicate that the CN scheme provides excellent numerical solutions for the problem under consideration. As such, the errors we obtained when comparing our numerical solutions to the DTM are in part due to a poor approximation of this series solution. Due to computational limitations we are able to obtain a limited number of terms before a lack of memory becomes an issue; are very few terms if one considers the work of other authors such as Darvishi et al. [2] and Aziz et al. [29], for instance. A comparison with the results obtained by Darvishi et al. [2] confirms that the CN scheme produces accurate solutions quickly and efficiently, without the requirement of an excessive amount of terms to be calculated and evaluated as per the semianalytical methods which have been employed in the literature. Unfortunately the addition of a Newton-Raphson predictor-corrector slows the scheme down excessively. By increasing one can see that the scheme does tend to the values obtained by the CN method; however in terms of computational time the scheme has shown itself to be inefficient.

#### 6. Dynamical Systems Analysis

##### 6.1. Stability Analysis on the Linear System

A stability analysis of a system allows us to determine whether or not a system is stable or will be stable if perturbed. In this subsection, we implement a von Neumann stability analysis which is a procedure used to analyse the stability of finite difference schemes by Fourier methods as applied to linear differential equations [47]. It can be noted that the method requires the solution to be periodic and a constant spacing . The method does not account for the boundary conditions.

The structure of the equation incorporates nonlinear sink terms; these do not allow for the use of a von Neumann stability analysis as a means of investigating the numerical scheme. As such we consider the equation without these terms; that is, we consider the linear scheme excluding sink terms. The purpose of this is to determine the extent to which the sink terms may impact on the stability of the schemes we have implemented. We have already noted the impact of the parameters and on the relevant schemes with regard to the maximum error obtained when comparing to the DTM. We will now continue to investigate this phenomenon.

Consider the linearised structure of (33) for the Crank-Nicolson scheme:where we definewith and are defined as before (34) and (35), respectively. Upon substitution of , where , into (41) we obtain

Introducing Euler’s identities into (43) and simplifying provide

Further simplification allows us to find the amplification factor

We note at this point that the von Neumann stability analysis is applicable to problems with constant coefficients, and yet is dependent on . However, we may still deem the matrix to be a matrix of constant coefficients since we know the domain of (i.e., ). Hence, upon substitution of for we obtain constant coefficients as required. We also note that .

Given that for and having definedwe find that since and , We also note that at we obtain and when similarly . Now we let , providing the amplification factor

We may conclude that since , without any required conditions on the relevant step sizes or parameters, the linear system is unconditionally stable.

This conclusion indicates that the structures of the coefficient matrices employed should not impact on the stability of the numerical solutions obtained once the nonlinear terms have been incorporated; rather we need to engage with the impact of the sink terms on the stability of the schemes considered. Thus, to fully investigate the stability of the dynamics of the equation, we now turn to a nonlinear dynamical systems analysis. In this manner we aim to gain insight into the impact of the nonlinear terms on the stability of the relevant schemes.

##### 6.2. Nonlinear Dynamical Systems Analysis

Before we are able to ascertain whether a simple linearisation of the steady state equation may be used as a means of analysing the dynamics of the nonlinear system, we need to determine what this linearised system and the relevant equilibrium points are. It must be noted that we consider the case where ; this is due to the physical relevance placed on the dynamics at the fin tip.

Firstly, as per [48], we restructure our steady state second-order ordinary differential equation into a system of first-order ordinary differential equations as follows:so that

We now turn to fixed point solutions or equilibrium solutions with the intention of investigating the dynamics of the system [49]. In the case of this nonlinear, nonautonomous system the fixed points (or equilibrium points) are defined by the vanishing of the vector field; that is,

Theorems for the stability of fixed points of the systemare presented by Coddington and Levinson [50], where is the autonomous part of the Jacobian and a constant matrix (upon substitution of the relevant equilibrium points and parameter values) and the vector function is continuous in and . In order to employ these theorems we write our system of equations (50) in the form (52):

As per (51) we obtain the equilibrium points :

We follow the work of Coddington and Levinson [50] as discussed in Nayfeh and Balachandran [48]. As such, we need only consider the eigenvalues of the matrix from (53) which we can write with general parameter values as follows: to be evaluated at the equilibrium points in order to classify them. From the form they take we can see two possible cases (depending on whether the term inside the square root is positive or negative): a saddle point or a centre.

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

The eigenvalues of at equilibrium point (54) are

We can see that this equilibrium point (given the ranges for the relevant parameter values) represents a saddle point which is unstable; see Figure 5(b).

**(a)**

**(b)**

Coddington and Levinson [50] have shown that if at least one eigenvalue of has a real part which is positive, then upon satisfying the requirement that, given any , there exist and such thatthe solution of (52) is not stable. In our case, we find that for the equilibrium point at least one eigenvalue of is positive and furthermore

We need to show that . In an attempt to do so, we maximise the term on the left of the inequality; that is, we consider such that where, given that , we define . This proves our requirement (58), confirming that equilibrium point (54) is not stable.

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

The eigenvalues of at this equilibrium point are lengthy for arbitrary values of , , , and . However, as pointed out above, the second case needs to be for . A variety of values were considered for the parameters , , , and , all leading to eigenvalues of the form expected, which are providing us with a stable centre; see Figure 5(a). At this stage we note that this case requires which is not physically relevant.

##### 6.3. Consideration of Viability of Linearised System

In the previous subsection we were able to determine the equilibrium points and classify them. We will now consider the feasibility of the dynamics of this linearised system describing the nonlinear dynamics of the system under discussion. Furthermore, we expand our stability analysis to consider the physical significance of the nonlinear terms and the degree, if at all, to which they impact on the scheme’s stability. This analysis has been used by Charrier-Mojtabi et al. [51] in the study of Soret-driven convection in a horizontal porous layer. Furthermore, Patel and Meher [52] performed a stability analysis via the fixed point theorem method for the Adomian decomposition Sumudu transformation method, when solving for the heat transfer in a solid and porous fin with temperature dependent internal heat generation.

In order to employ a linearisation approach to study the dynamics of this nonlinear nonautonomous system we consider a Taylor series expansion of near the equilibrium point [53]. This may be constructed as follows: