Natural Convection Flow over a Vertical Permeable Circular Cone with Uniform Surface Heat Flux in Temperature-Dependent Viscosity with Three-Fold Solutions within the Boundary Layer

Abstract

The aim of this study is to investigate the effects of temperature-dependent viscosity on the natural convection flow from a vertical permeable circular cone with uniform heat flux. As part of numerical computation, the governing boundary layer equations are transformed into a non-dimensional form. The resulting nonlinear system of partial differential equations is then reduced to local non-similarity equations which are solved computationally by three different solution methodologies, namely, (i) perturbation solution for small transpiration parameter (ξ), (ii) asymptotic solution for large ξ, and (iii) the implicit finite difference method together with a Keller box scheme for all ξ. The numerical results of the velocity and viscosity profiles of the fluid are displayed graphically with heat transfer characteristics. The shearing stress in terms of the local skin-friction coefficient and the rate of heat transfer in terms of the local Nusselt number (Nu) are given in tabular form for the viscosity parameter (ε) and the Prandtl number (Pr). The viscosity is a linear function of temperature which is valid for small Prandtl numbers (Pr). The three-fold solutions were compared as part of the validations with various ranges of Pr numbers. Overall, good agreements were established. The major finding of the research provides a better demonstration of how temperature-dependent viscosity affects the natural convective flow. It was found that increasing Pr, ξ, and ε decrease the local skin-friction coefficient, but ξ has more influence on increasing the rate of heat transfer, as the effect of ε was erratic at small and large ξ. Furthermore, at the variable Pr, a large ξ increased the local maxima of viscosity at large extents, particularly at low Pr, but the effect on temperature distribution was found to be less significant under the same condition. However, at variable ε and fixed Pr, the temperature distribution was observed to be more influenced by ε at small ξ, whereas large ξ dominated this scheme significantly regardless of the variation in ε. The validations through three-fold solutions act as evidence of the accuracy and versatility of the current approach.

1. Introduction

Natural convection occurs in the environment. Different industrial applications, closed containers, or any type of chamber, such as a greenhouse, are also good examples of locations where convective heat transfer occurs. In order to characterize different heat transfer applications, heat flux is one of the major indices to be considered. Uniform heat flux is often considered in geophysics as the soil is a complex porous material. The ability of a soil specimen or soil surface to conduct heat, i.e., thermal conductivity, is significant to agricultural systems, as the growth of plant canopies is directly dependent on soil temperature. Soil has moisture contents, and moisture (water) exhibits temperature-dependent viscosity. The viscosity of water hence directly influences the surface energy balance. Therefore, heat flux through soil or any porous medium is significant in environmental or agricultural applications. Due to different destructive geophysical applications, such as irrigations, compactions, and sampling, the soil often loses its texture. Due to the complexities in the geometry, the soil is often assigned a representative element to discuss different properties, such as electrical or thermal conduction [1,2,3]. Therefore, heat transfer applications in different geometries representing the soil, or any porous material, have attracted great interest. Researchers have been studying natural convection for decades to solve different heat transfer applications in different industries, such as renewable energy, agriculture, or environment, to name a few [4,5,6,7,8,9]. In computational fluid dynamics (CFD) research, the study of natural convective flow comprises different geometries and representative elements. The physical models or geometries are often considered to be square or rectangular cavities with heated or cooled walls satisfying certain boundary conditions [10,11,12,13]. On the other hand, the representation of elements is particularly defined by different dimensionless parameters which are often varied to conduct sensitivity tests to understand the behavior of the fluid flow under various conditions in order to perform the computation.

The study of heat transfer with different cross-sections has gained strong interest. It is realistic to observe the patterns of fluid flow under various circumstances around surfaces of different cross-sections, such as cylinders or cones, with possible inclusions of permeability [14,15,16,17,18]. Merk and Prins [19,20] are often credited to be the pioneers in developing solutions for the fluid flow past a vertical cone considering an axisymmetric form. Later, Braun et al. [21] and Hering and Grosh [22] studied and developed the numerical model further by keeping the Prandtl number (Pr) between 0.7 and 0.72. In any case, all the studies mentioned above agreed on the existence of similar solutions for natural convective flow from the vertical cone. Later, in the 2000s, the series of works done by Hossain and Paul [23,24], as well as Hossain et al. [25], highlighted the investigation of the natural convective flow from a heated vertical permeable circular cone by considering both uniform [25] and non-uniform [23,24] surface heat flux and temperature. The three-fold solutions were obtained by the finite difference method, series solution method, and asymptotic solution method.

Most of the relevant published research works focused on sensitivity analyses on constant viscosity. It is a well-known fact that the viscosity of the fluid does not remain constant as a function of temperature. For example, Cebeci and Bradshaw [26] provided a detailed chart outlining the changes in the water viscosity as a function of temperature, where the value of viscosity decreased by approximately 240% as the temperature of the fluid was elevated from 10 °C to 50 °C. A similar conclusion was drawn for other fluids as well, such as air [26]. Therefore, a proper numerical establishment and computational model in state-of-the-art research should take variable viscosity into account along with other physical attributes. Some of the works have been found in the literature where viscosity was varied, such as [27,28]. Ling and Dybbs [29] mentioned varying viscosity inversely for large Pr numbers, but the possible implications of including small Pr numbers were not properly investigated. In that case, many important fluids, such as liquid sodium (Pr ≈ 0.004) and mercury (Pr ≈ 0.03), are being left out of consideration. Fluids with low Pr numbers are thought to have small kinematic viscosity or a greater heat diffusivity. Therefore, fluids with Pr << 1 will have heat being diffused quicker than the velocity. As a result, the numerical models that are valid for high Pr numbers have significantly less versatility. In terms of variable viscosity with temperature, Rahman et al. [18] investigated the natural convective flow along the vertical wavy cone where viscosity was considered an exponential function of temperature. In a similar geometry, Thohura et al. [30] emphasized the effects of the temperature-dependent thermal conductivity. Recently, Khan et al. [31] studied the entropy generation of fluids within the incompressible boundary layer to understand the changes in velocity and concentration profiles, which could only be explored due to the consideration of temperature-dependent viscosity. Gladys and Reddy [32] also accentuated the role of the temperature-dependent viscosity of non-Newtonian nanofluids through the accelerating vertical plate to analyze the nonlinear buoyancy impacts.

In the literature, most of the works consider the constant viscosity of the fluid within the boundary layer and failed to prove the accuracy of the model in various ranges of Pr numbers. The current study considers both small (for example, 0.05) and large (for example, 0.7) Pr numbers. Therefore, it is essential to validate the current approach with different numerical parametric tests to gain more confidence in the accuracy and versatility of the approach. Therefore, more than one numerical solution/validation should be conducted to gain a better understanding of the accuracy. Among different numerical techniques, the Keller box method has been a popular and widely accepted numerical technique for nearly five decades [33], with more applications being added to this scheme recently. For example, Kamran et al. [34] have considered the Keller box approach to describe the Jefferey–Hamel flow by considering different non-dimensional parameters to obtain solutions. Reddy et al. [35] also obtained implicit finite difference results by the Keller box technique to study the Joule heating and associated chemical reactions on the magneto Casson nanofluid. Therefore, another evidence of expanding applications of the Keller box method was duly noted. On the other hand, perturbation and asymptotic solutions were also found to be highly efficient for small and large transpiration parameters (ξ), respectively, for the last few years [36]. Therefore, the validations through a combination of conventional and modern approaches should be adequate to prove the accuracy of the current approach and will provide a new dimension to future study.

The present study aims to investigate natural convection flow in a vertical permeable cone with uniform heat surface flux as viscosity varies in the computational model. The present model is also valid for fluids with a low Prandtl number (Pr), thus resolving the shortcomings of most of the literature in the past decades. The non-dimensional viscosity-variation parameter (ε) has been included in the model, along with the pseudo-similarity variable (η), to tune the fluid characteristics in the model. The suction parameter (ξ) has been included at various ranges to observe the effect of changes in the rate of heat transfer and shear stress. The sensitivity analyses have been conducted with both constant and variable Pr numbers to observe the behaviors of the temperature and velocity distribution of the fluid within the boundary layer. The model has been validated with both fixed and variable Pr numbers to showcase the uniformity as well as the accuracy of the approach. The three-fold numerical solutions have been conducted by implementing the aforementioned parameters. The implicit finite difference method together with the Keller box scheme for all ξ has been conducted, followed by perturbation for small ξ and, finally, the asymptotic solution for large ξ. To the authors’ knowledge, there has not been any published work that considers all the characteristics and scopes mentioned above.

2. Formulation of the Problem

A steady two-dimensional free convective flow in laminar form is considered over a vertical permeable circular cone of radius r with uniform heat flux. The cone is immersed in a viscous and incompressible fluid with temperature-dependent viscosity, for the numerical investigation. It is assumed that the surface heat flux of the cone is q_w. Here, T_∞ is the ambient temperature of the fluid, and the considered geometry is illustrated in Figure 1.

The equations governing the flow are [23,24,25,30,37]:

(1)$\mathrm{Continuity} \mathrm{equation}: \frac{\partial \left(ru\right)}{\partial x}+\frac{\partial \left(rv\right)}{\partial y}=0$

(2) $\mathrm{Momentum} \mathrm{equation}: \rho \left(u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=\frac{\partial }{\partial y}\left(\mu \frac{\partial u}{\partial y}\right)+\rho \mathit{g}\beta \left(T-T_{\infty }\right)\sin \gamma$

(3) $\mathrm{Energy} \mathrm{equation}: u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\frac{k}{\rho C_p}\frac{{\partial }^{2}T}{\partial y^2}$

The boundary conditions of Equations (1)–(3) are provided below [23]:

(4)$\begin{array}{c}u=0,\hspace{0.17em}v=-V_w,\hspace{1em}{q}_w=-k\left(\frac{\partial T}{\partial y}\right),\hspace{0.17em}\hspace{0.17em}\mathrm{at}\hspace{0.17em}y=0\\ u\to 0,\hspace{1em}T\to T_{\infty }\hspace{1em}\mathrm{as}\hspace{1em}y\to \infty \end{array}\}$

where $\left(u,v\right)$ are velocity components along the $\left(x,y\right)$ axes, g is the acceleration due to gravity, ρ is the density, γ is the cone apex half-angle, k is the thermal conductivity, β is the coefficient of thermal expansion, μ(T) is the temperature-dependent viscosity of the fluid, where T denotes the temperature, and V is the transpiration velocity, which is positive for suction and negative for the injection of fluid through the cone. In the current research, only the suction case has been considered, and therefore V_w has been considered to be positive throughout.

Viscosity variations have multifarious forms. However, the one proposed by Charraudeau [38] is one of the most accepted ones and has been taken into account in this study. The expression could be presented as the following:

(5)$\mu ={\mu }_{\infty }\left[1+\frac{1}{{\mu }_{fil}}{\left(\frac{\partial \mu }{\partial T}\right)}_{fil}\left(T-T_{\infty }\right)\right]$

where the suffix fil represents the film temperature of the fluid.

If suction is non-existent in the model, the boundary layer acts as the free convection boundary layer near the leading edge, although improved downstream suction will be able to dominate the flow to a greater extent. As a result, the following transformations are required [23]:

(6)$$\begin{gathered} \psi ={\nu }_{\infty }rG{r}_x^{1/5}\left[f\left(\xi ,\eta \right)+\frac{1}{2}\xi \right],\hspace{1em}T-T\infty =\frac{q_{w}x}{k}G{r}_x^{-1/5}\theta \left(\xi ,\eta \right) \\ \eta =\frac{y}{x}G{r}_x^{1/5},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\xi =\frac{V_{w}x}{{\nu }_{\infty }}G{r}_x^{-1/5},\hspace{0.17em}\hspace{1em}G{r}_x=\frac{\mathit{g}\beta \cos \gamma x^4}{k{\nu }_{\infty }^2},\hspace{0.17em}r=x\sin \gamma \end{gathered}$$

where ν_∞(=μ_∞/ρ) is the reference kinematic viscosity, Gr_x is the local Grashof number, ξ is the dimensionless transpiration parameter, η is the pseudo-similarity variable, f(ξ,η) and θ(ξ,η) are the non-dimensional stream and temperature function, respectively, and ψ is the stream function defined by the following:

(7)$u=\frac{1}{r}\frac{\partial \psi }{\partial y},\hspace{1em}v=-\frac{1}{r}\frac{\partial \psi }{\partial x}$

Substituting (6) into Equations (1)–(5) and after some algebraic calculation, the following transformed equations are obtained:

(8)$\left(1+\epsilon \xi \theta \right)f^{‴}+\frac{9}{5}f{f}^{″}-\frac{3}{5}{f^{\prime }}^2+\epsilon \xi {\theta }^{\prime }{f}^{″}+\theta +\xi f^{″}=\frac{1}{5}\xi \left(f^{\prime }\frac{\partial f^{\prime }}{\partial \xi }-f^{″}\frac{\partial f}{\partial \xi }\right)$

(9)$\frac{1}{\mathit{Pr}}{\theta }^{″}+\frac{9}{5}f{\theta }^{\prime }-\frac{1}{5}{f}^{\prime }\theta +\xi {\theta }^{\prime }=\frac{1}{5}\xi \left(f^{\prime }\frac{\partial \theta }{\partial \xi }-{\theta }^{\prime }\frac{\partial f}{\partial \xi }\right)$

And the boundary conditions are presented hereby:

(10)$\begin{array}{c}f=f^{\prime }=0,\hspace{1em}{\theta }^{\prime }=-1,\hspace{1em}\mathrm{at}\hspace{1em}\eta =0\\ f^{\prime }\to 0,\hspace{1em}\theta \to 0,\hspace{1em}\mathrm{as}\hspace{1em}\eta \to \infty \hspace{1em}\end{array}\}$

where ε is the viscosity-variation parameter, and Pr is the Prandtl number defined as the following:

(11)$\epsilon =\frac{1}{{\mu }_{fil}}{\left(\frac{\partial \mu }{\partial T}\right)}_{fil}\frac{q_w{\nu }_{\infty }}{kV}\hspace{1em}\mathrm{and}\hspace{0.17em}\mathit{Pr}=\frac{{\mu }_{\infty }{C}_p}{k}$

The local non-similar partial differential equations presented in Equations (8) and (9) need to be solved after finding the boundary conditions (Equation (10)), where the latter could be obtained by the Keller box method, as verified in the literature [25]. The results are computed by considering uniform grids in the ξ-direction with 1001 grid-points considered in the 0 ≤ ξ ≤ 20. Further information on iterations and simulation strategies has been included in Section 3.1.

After calculating the respective values of f, θ, and their derivatives, the fundamental approach is to calculate the local skin-friction coefficient and the local Nu number from the following expressions:

(12)$C_{fx}G{r}_x^{1/5}=f^{″}\left(\xi ,0\right)$

(13)$N{u}_{x}G{r}_x^{-1/5}=\frac{1}{\theta \left(\xi ,0\right)}$

The results obtained by this method are presented in tabular form in Table 1 for different values of the viscosity variation parameter ε (=0.0, 0.5, 2.0) and Prandtl number Pr (=0.05, 0.1, 0.7). In the following sections, the solutions for small and large ξ have been presented, and are numerically valid and accurate in both the neighboring leading edge (small ξ) and downstream region (ξ).

2.1. Solution for Small ξ (ξ << 1)

The value ξ for x is small near the leading edge or small V or both, and the series solution of Equations (8) and (9) can be found by applying the perturbation method considering x as a perturbation parameter. In order to establish such an objective, the functions f(ξ,η) and θ(ξ,η) in powers of ξ are expanded, and therefore:

(14)$f\left(\xi ,\eta \right)={\displaystyle \sum _{i=0}^{\infty }{\xi }^i{f}_i}\left(\eta \right) \mathrm{and} \theta \left(\xi ,\eta \right)={\displaystyle \sum _{i=0}^{\infty }{\xi }^i{\theta }_i}\left(\eta \right)$

Substituting the expansions from Equation (14) into Equations (8) and (9) and equating the various powers of ξ up to 0(ξ²), the following batches of equations are obtained:

(15)$f_0^{‴}+\frac{9}{5}{f}_0{f}_0^{″}-\frac{3}{5}{f_0^{\prime }}^2+{\theta }_0=0$

(16)$\frac{1}{\mathit{Pr}}{\theta }_0^{″}+\frac{9}{5}{f}_0{\theta }_0^{\prime }-\frac{1}{5}{f}_0^{\prime }{\theta }_0=0$

(17)$\begin{array}{l}{f}_0\left(0\right)=f_0^{\prime }\left(0\right)=0,\hspace{1em}{\theta }_0^{\prime }\left(0\right)=-1\\ f_0^{\prime }(\infty )={\theta }_0(\infty )=0\end{array}$

(18)$f_1^{‴}+\epsilon \left({\theta }_0^{\prime }{f}_0^{‴}+{\theta }_0^{\prime }{f}_0^{″}\right)+\frac{9}{5}{f}_0{f}_1^{″}-\frac{7}{5}{f}_0^{\prime }{f}_1^{\prime }+2f_0^{″}{f}_1+{\theta }_1+f_0^{″}=0$

(19)$\frac{1}{\mathit{Pr}}{\theta }_1^{″}+\frac{9}{5}{f}_0{\theta }_1^{\prime }-\frac{2}{5}{f}_0^{\prime }{\theta }_1-\frac{1}{5}{\theta }_0{f}_1^{\prime }+2{\theta }_0^{\prime }{f}_1+{\theta }_0^{\prime }=0$

(20)$\begin{array}{l}{f}_1\left(0\right)=f_1^{\prime }\left(0\right)=0,\hspace{1em}{\theta }_1^{\prime }\left(0\right)=0\\ f_1^{\prime }(\infty )={\theta }_1(\infty )=0\end{array}$

And

(21)$\begin{array}{l}{f}_2^{‴}+\epsilon \left({\theta }_0{f}_1^{‴}+{\theta }_1{f}_0^{‴}+{\theta }_0^{\prime }{f}_1^{″}+{\theta }_1^{\prime }{f}_0^{″}\right)+\frac{9}{5}{f}_0{f}_2^{″}+\frac{11}{5}{f}_0^{″}{f}_2\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\frac{4}{5}{f_1^{\prime }}^2-\frac{8}{5}{f}_0^{\prime }{f}_2^{\prime }+2f_1^{″}{f}_1+{\theta }_2+f_1^{″}=0\end{array}$

(22)$\frac{1}{\mathit{Pr}}{\theta }_{21}^{″}+\frac{9}{5}{f}_0{\theta }_2^{\prime }+\frac{11}{5}{\theta }^{\prime }{f}_2{}_0-\frac{2}{5}{f}_1^{\prime }{\theta }_1-\frac{1}{5}{\theta }_0{f}_2^{\prime }-\frac{3}{5}{f}_0^{\prime }{\theta }_2+2{\theta }_1^{\prime }{f}_1+{\theta }_1^{\prime }=0$

(23)$\begin{array}{l}{f}_2\left(0\right)=f_2^{\prime }\left(0\right)=0,\hspace{1em}{\theta }_2^{\prime }\left(0\right)=0\\ f_2^{\prime }(\infty )={\theta }_2(\infty )=0\end{array}$

Equations (15)–(23) are solved pair-wise one after another. The solutions are obtained by combining the famous Runge–Kutta–Butcher [39] initial value solver with the Nachts–eim–Swigert iteration scheme [40]. Thus, solutions are found for f_i and θ_i (i = 0, 1, 2) and their respective derivatives.

Once the values of f_i and θ_i for i = 0, 1, 2 and their derivatives are obtained, the local skin-friction coefficient and the local Nusselt number are calculated from the following expressions:

(24)$C_{fx}G{r}_x^{1/5}=f^{″}\left(\xi ,0\right)=f_0^{″}\left(0\right)+\xi f_1^{\prime }\left(0\right)+{\xi }^2{f}_2^{″}\left(0\right)$

(25)$N{u}_{x}G{r}_x^{-1/5}=\frac{1}{\theta \left(\xi ,0\right)}=1/\left[{\theta }_0\left(0\right)+\xi {\theta }_1\left(0\right)+{\xi }^2{\theta }_2\left(0\right)\right]$

The comprehensive values calculated from Equations (24) and (25) are presented in Table 1. The comparison also served as part of the validations of the present approach.

2.2. Solution for Large ξ

Attention has been given in this part to the behavior of the solution to the Equations (8) and (9) when ξ is considerably large. By the order of magnitude analysis of the various terms, ξθ′f″ in (8) and ξθ′ in (9) were found to be the largest. However, in both of their equations, both numerical terms need to be balanced mathematically. The balancing part is performed by assuming η to be a small parameter, which will eventually make η-derivatives larger. It is also important to determine the standard scaling approach considering θ = O(ξ − 1) as ξ→ ∞. After balancing the f‴, θ, and ξθ′f″ terms in (9), η = O(ξ − 1) and f = O(ξ − 4) as ξ→ ∞ would be found, which would also serve as a confirmation of the accuracy in the balancing. As a result, the following expressions are substituted:

(26)$f={\xi }^{-4}\tilde{f}\left(\xi ,\tilde{\eta }\right),\hspace{1em}\theta ={\xi }^{-1}\tilde{\theta }\left(\xi ,\tilde{\eta }\right),\hspace{1em}\tilde{\eta }=\xi \eta$

Substituting this transformation into Equations (8) and (9), we get the following equations:

(27)$\left(1+\epsilon \tilde{\theta }\right){\tilde{f}}^{‴}+{\xi }^{-5}\tilde{f}{\tilde{f}}^{″}+\tilde{\theta }+{\tilde{f}}^{″}+\epsilon {\tilde{\theta }}^{\prime }{\tilde{f}}^{″}=\frac{1}{5}{\xi }^{-4}\left({\tilde{f}}^{\prime }\frac{\partial {\tilde{f}}^{\prime }}{\partial \xi }-{\tilde{f}}^{″}\frac{\partial \tilde{f}}{\partial \xi }\right)$

(28)$\frac{1}{\mathit{Pr}}{\tilde{\theta }}^{″}+{\xi }^{-5}\tilde{f}{\tilde{\theta }}^{\prime }+{\tilde{\theta }}^{\prime }=\frac{1}{5}{\xi }^{-4}\left({\tilde{f}}^{\prime }\frac{\partial \tilde{\theta }}{\partial \xi }-{\tilde{\theta }}^{\prime }\frac{\partial \tilde{f}}{\partial \xi }\right)$

The corresponding boundary conditions are

(29)$\begin{array}{l}\tilde{f}\left(\xi ,0\right)={\tilde{f}}^{\prime }\left(\xi ,0\right)=0,\hspace{1em}{\tilde{\theta }}^{\prime }\left(\xi ,0\right)=-1\\ {\tilde{f}}^{\prime }\left(\xi ,\infty \right)=0,\hspace{1em}\tilde{\theta }\left(\xi ,\infty \right)=0\end{array}$

where primes describe differentiation with respect to $\tilde{\eta }$. Equations (27) and (28) are solved in terms of an inverse power of series in ξ.

The functions $\tilde{f}\left(\xi ,\tilde{\eta }\right)$ and $\tilde{\theta }\left(\xi ,\tilde{\eta }\right)$ are expanded in the power series in the negative powers of ξ, considering that ξ is large, which yields the following:

(30)$\tilde{f}\left(\xi ,\tilde{\eta }\right)={\displaystyle \sum _{i=0}^{\infty }{\xi }^{-5i}{\tilde{f}}_i\left(\tilde{\eta }\right)} \mathrm{and} \tilde{\theta }\left(\xi ,\tilde{\eta }\right)={\displaystyle \sum _{i=0}^{\infty }{\xi }^{-5i}{\tilde{\theta }}_i\left(\tilde{\eta }\right)}$

Now substituting the above expansions into Equations (27)–(29) and equating the coefficient of power of ξ⁰ and ξ¹, we get the following equations

(31)$\left(1+\epsilon {\tilde{\theta }}_0\right){\tilde{f}}_0^{‴}+{\tilde{f}}_0^{″}+{\tilde{\theta }}_0+\epsilon {\tilde{\theta }}_0^{\prime }{\tilde{f}}_0^{″}=0$

(32)$\frac{1}{\mathit{Pr}}{\tilde{\theta }}_0^{″}+{\tilde{\theta }}_0^{\prime }=0$

(33)$\begin{array}{l}{\tilde{f}}_0\left(0\right)={\tilde{f}}_0^{\prime }\left(0\right)=0,\hspace{1em}{\tilde{\theta }}_0^{\prime }\left(0\right)=-1\\ {\tilde{f}}_0^{\prime }(\infty )=0,\hspace{1em}{\tilde{\theta }}_0(\infty )=0\end{array}$

(34)$\left(1+\epsilon {\tilde{\theta }}_0\right){{\tilde{f}}^{‴}}_1+{\tilde{f}}_0{\tilde{f}}_0^{″}+{\tilde{\theta }}_1+\epsilon \left({\tilde{\theta }}_1{\tilde{f}}_0^{″}+{\tilde{\theta }}_1^{\prime }{\tilde{f}}_0^{″}+{\tilde{\theta }}_0^{\prime }{\tilde{f}}_1^{″}\right)=0$

(35)$\frac{1}{\mathit{Pr}}{\tilde{\theta }}_1^{″}+2{\tilde{\theta }}_0^{\prime }{\tilde{\theta }}_1^{\prime }+{\tilde{f}}_0^{\prime }{\tilde{\theta }}_0^{\prime }+{\tilde{\theta }}_1^{\prime }=0$

(36)$\begin{array}{l}{\tilde{f}}_1\left(0\right)={\tilde{f}}_1^{\prime }\left(0\right)=0,\hspace{1em}{\tilde{\theta }}_1^{\prime }\left(0\right)=0\\ {\tilde{f}}_1^{\prime }(\infty )=0,\hspace{1em}{\tilde{\theta }}_1(\infty )=0\end{array}$

The local skin-friction and the local Nusselt number are as follows:

(37)$C_{fx}G{r}_x^{1/5}=f^{″}\left(\xi ,0\right)={\xi }^{-2}\left[{\tilde{f}}_0^{″}\left(0\right)+{\xi }^{-5}{\tilde{f}}_1^{″}\left(0\right)\right]$

(38)$N{u}_{x}G{r}_x^{-1/5}=\frac{1}{\theta \left(\xi ,0\right)}=\frac{\xi }{{\tilde{\theta }}_0\left(0\right)+{\xi }^{-5}{\tilde{\theta }}_1\left(0\right)}$

The asymptotic solutions obtained from (31)–(36) in terms of local skin-friction and local Nusselt number and compared with the solution of the finite difference method in Table 1.

3. Results and Discussion

3.1. Computational Methods

The analyses have been conducted by tuning the non-dimensional empirical parameters at different stages with both fixed and variable Pr numbers. The computation was done by Fortran 90 [41]. The finite difference solutions were obtained by the Keller box scheme [31,32,33,34,35,36,37,38,39,40,41,42], which is considered to be one of the most accurate implicit finite difference techniques in computational mathematics. Refer to [37] for a detailed explanation of the Keller box method. In this method, the non-linear system of the partial differential equations governing the fluid flow is solved, and an assumption is made in terms of the functions and the derivatives to express the first-order equations. The derivatives in terms of the central differences are approximated in both co-ordinate directions. This part is particularly obtained by denoting the mesh points in the (ξ, η) plane by x_i and η_I, where i = 1, 2, 3, …, M and j = 1, 2, 3, …, N [37], followed by central difference approximations where the equations containing x are centred explicitly at (ξ_i−1/2, η_j−1/2) and the rest at (ξ_i, η_j−1/2), where η_j−1/2 = (η_j + η_j−1)/2, for instance. It yields a batch of non-linear differential equations for the unknowns at xi in as a function of ξ_i−1. Then, those equations are subject to linearization by Newton’s quasi-linearization method with a view to solving them by taking a block-tridiagonal algorithm into account. The initial iteration of the converged solution is considered to be at ξ = ξ_i−1. At the commencement, i.e., ξ = 0, the guess profiles for all the considered variables are provided, and then the Keller box method comes into the computation to solve the governing ordinary differential equations. By obtaining the lower stagnation point solution, the steps could be performed along with the boundary layer of the geometry [37]. The iteration is terminated once the target difference in velocity and temperature computation is reached, and in this study it was |δf ⁱ| ≤ 10⁻⁶, where i represents the number of iterations. It should be mentioned here that computations were not conducted by a uniform grid in the y-direction. However, a non-uniform grid was considered, which was then defined by η_j = sinh(( j − 1)/p), with j = 1, 2, …, 301 and p = 100.

3.2. Overview of Numerical Analyses

At first, the validations were conducted by comparing the three different types of solutions considered in this study. A different type of 3D analyses was conducted to assess the correlation coefficient, and good R² values were obtained overall. As part of the discussion, the streamlines and isothermal behavior have been presented. After that, the influence of parameters has been investigated on the viscosity, velocity, and temperature distributions. At the end of each analysis, a conclusion has been drawn outlining the significant and non-significant impact of specified parameters under different circumstances.

3.3. Comparing Finite Difference Solutions with Perturbation and Asymptotic Solutions

3.3.1. Validation and Comparison at Fixed Pr

To assess the accuracy of the approach, the typical finite difference solutions have been compared against the present approach with small and large suction parameters (ξ). To achieve such an objective, two different values of ε were assigned (0 and 2), and Pr = 0.7 was considered. In other words, the comparisons have been made in both the absence and presence of ε for the same type of fluid. The cross-validations were conducted by observing the values of two different parameters namely, the local skin-friction coefficient (C_fxGr_x^1/5) and the local Nusselt number (Nu_xGr_x^−1/5). Table 1 contains the comprehensive comparisons of the aforementioned conditions. In general, there is hardly any big difference between the finite difference solutions and the asymptotic solutions. The slight deviation in the values could be attributed to the percentage error in the approach.

According to Table 1, as ξ increased, C_fxGr_x^1/5 decreased, regardless of the ε values. However, C_fxGr_x^1/5 values were found to be lower in the presence of ε, thus confirming the role of the viscosity-variation parameter on the reduction of the local skin-friction coefficient. This behavior could be explained in terms of the viscosity and temperature of the fluid. As ε increased, the fluid temperature within the boundary layer increased, which led to increased viscosity. Therefore, C_fxGr_x^1/5 decreased further. Similar patterns were observed for both finite difference and asymptotic solutions.

Meanwhile, Nu_xGr_x^−1/5 exhibited erratic trends. In general, Nu_xGr_x^−1/5 increased as ξ increased when ε = 0 was considered. However, as ε = 2 was assigned, Nu_xGr_x^−1/5 values had a different pattern. Initially, Nu_xGr_x^−1/5 kept on increasing as ξ increased when ε = 2, and the values of local Nu number were lower than those achieved at ε = 0. However, at ξ ≥ 2, the values of Nu_xGr_x^−1/5 exhibited values higher than Nu_xGr_x^−1/5 at ε = 0.

Furthermore, surface analyses were conducted to understand the rational behavior as well as the accuracy of the approach. The purpose was to add one type of extra validation in the current approach to check the model predictions for a different type of fluid. Therefore, the local skin-friction coefficients and local Nu number were calculated at Pr = 0.7 and the numerical accuracy was assessed by observing the coefficient of correlation (R²) at ε = 0 and ε = 2. The surface analyses were performed by the data analysis software Origin Pro, developed by OriginLab Corporation [42]. In the 3D analyses, solutions from the finite difference method and perturbation-asymptotic solutions were illustrated in the same frame to observe the analytical behavior as a function of the suction parameter (ξ). In general, an R² between 0.97 and 0.98 was obtained, which provides more confidence in the calculative approach. Figure 2 depicts the predicted surfaces for C_fxGr_x^1/5 at different ξ in the absence of ε (Figure 2a) and with an assigned value of ε (Figure 2b). Most of the obtained points were found on the surface. In short, the marginal calculation error was noticed. On the other hand, a similar analysis was conducted in terms of Nu_xGr_x^−1/5 in Figure 3 under the same ε conditions. It is evident that Nu_xGr_x^−1/5 is increasing concurrently as ξ was increasing in Figure 3a,b, in line with the argument presented in Table 1, where Pr = 0.2. Furthermore, the calculated values from the asymptotic and perturbation solutions for small and large ξ were found to be closer than those from the finite difference approach. The missing values have not been included in the illustrations, and therefore the input sample size is not consistent in all figures.

3.3.2. Comparison at Variable Pr

Figure 4 depicts the comparison of solutions at two different Pr values to confirm the accuracy of the approach for two different types of fluids. It could be observed that the agreements were excellent in terms of both C_fxGr_x^1/5 (Figure 4a) and Nu_xGr_x^−1/5 at ξ ∈ [0,20] (Figure 4b)

It was observed that as Pr increased C_fxGr_x^1/5 values decreased, while Nu_xGr_x^−1/5 increased concurrently. As ξ increased, the temperature of the fluid increased as mentioned earlier, and hence at ξ = 20 the values of C_fxGr_x^1/5 were close to 0 due to the dominance of the suction parameter, whereas at the same suction parameter value, Nu_xGr_x^−1/5 was found to be the highest regardless of Pr values. However, it should be mentioned that as ξ kept increasing, the difference between the Nu_xGr_x^−1/5 of the fluid corresponding to Pr = 0.1 became wider than the former, with Pr = 0.05.

Since the temperature of the fluid kept increasing at increasing ξ values, the shear stress, i.e., the local skin-friction coefficient, decreased, and to balance the temperature difference, the rate of the heat transfer, i.e., the local Nu number, increased.

3.4. Development of Streamlines

This part of the discussion refers to fluids corresponding to two different types of Pr numbers (Pr = 0.05, 0.1) to depict the changes in the streamlines as a function of the suction parameter ξ. The streamlines illustrate the flow pattern within the momentum boundary layer. The changes were analyzed by observing the variations in the values of the flow rate (ψ). Two different viscosity-variation parameters (ε = 0.0, 0.5) were considered for each type of fluid to include the effect of ε on the streamlines. It was anticipated that the differences in the shapes of streamlines in both ε values would be quite marginal at low ξ(ξ < 4), due to the lower effect of the suction parameter. However, as ξ increases to a considerably larger value (ξ > 4), the differences between the shapes of streamlines of ε = 0.0 and ε = 0.5 increase and do not overlap. Furthermore, since ε is the viscosity-variation parameter, its absence will lead to the maximum flow rate, regardless of the Pr values. This behavior was expected due to the effect of viscosity on the fluid movement, which has been explained in the later sections, where the distribution of velocity profiles has been comprehensively discussed.

According to Figure 5, at lower ξ the streamlines corresponding to both ε values overlap, indicating the non-significant impact of the suction parameter on the fluid characteristics. A similar behavior was found for both Pr = 0.05 (Figure 5a) and Pr = 0.1 (Figure 5b). However, as ξ kept on increasing, the differences in the maximum values of the fluid flow rate of each segment corresponding to ε = 0.0 and ε = 0.5 were increasing gradually, which was mentioned in the first paragraph. This behavior could be also explained in terms of the effect of ε on the streamline’s development. In the presence of the viscosity-variation parameter, the streamlines of the fluid seem to occupy less surface area due to the effect of ε, which eventually restricts the local maximum of the velocity distribution of the fluid. As per Figure 5a, while the maximum flow rate was found to be 42.90 in the absence of ε, the value plummeted to 13.11 (Figure 5b). On the other hand, the fluid corresponding to Pr = 0.1 exhibited a different trend in the pattern of streamlines. While the streamlines corresponding to the maximum flow rate were more aligned to the right side of the boundary layer, the maximum flow rate curve at Pr = 0.1 shifted to the center of the boundary, indicating the effect of Pr on the fluid characteristic. If Pr is increased, the momentum and thickness of the thermal boundary layer decreased, and hence the flow rate got reduced.

3.5. Isothermal Behavior within the Boundary Layer

After investigating the behavior of the streamlines, a similar parametric effect was studied in terms of isotherms, illustrated in Figure 6. The isotherms depict the temperature distribution within the boundary layer. It was expected that the isothermal lines closer to the wall would be equal to the boundary temperature, which would lead to lines without any significant peak values. The isothermal line considerably far from the wall boundary would have peak values before plummeting gradually as the lines tended to come near the boundary layer.

The isotherms appended in Figure 6 confirm the aforementioned statement. The isothermal lines far from the boundary layer could be seen with a peak before gradually starting to reduce. There is a slight deviation in the values between ε = 0.0 and ε = 0.5, where the presence of ε reasonably increased the isothermal lines due to the influence of the viscosity-variation. However, this was only visible when the lines were at an optimal distance from the thermal boundary layer, where the influence of the latter would be less significant. A similar behavior was obtained in terms of streamlines as well where the presence of ε shifted the streamlines up compared to the streamlines corresponding to ε = 0.0. However, as the fluid motion decreases with an increase from Pr = 0.05 (Figure 6a) to Pr = 0.1 (Figure 6b), the peak values of isothermal lines decrease (visible ones) due to the reduction of the flow rate. The significant reduction in the maximum flow rate was also recorded in Figure 6b, and hence the consistency of the streamlines and isothermal analyses could be confirmed before heading towards the detailed parametric analyses in the following sections.

3.6. Impact of ε, η, and ξ Fluid Characteristics at Fixed Pr

3.6.1. Influence on Viscosity Distribution

Figure 7 depicts the effects of both the viscosity-variation parameter (ε) and the suction parameter (ξ) on the dimensionless viscosity. The pseudo-similarity variable (η) was varied from 0 to a minimum of 6 (for ξ = 10) and a maximum of 10.0 (ξ = 1), depending on where the viscosity started to remain indifferent regardless of the η values. In Figure 7a,b, Pr = 0.1 was considered to maintain the consistency of the fluid characteristics.

In general, it was found that as η increased, viscosity exhibited the local maximum values at the maximum ε, in this case ε = 2, as well as at η = 0. However, the local maximum viscosity was obtained at approximately 20.5 at ξ = 10, whereas it was found to be around 7.6 at ξ = 1. This behavior could be attributed to the fact that as the suction parameter increases, the viscosity exhibits the maximum values due to the absence of a pseudo-similarity variable which acts as a counterincentive to viscosity. On the other hand, ξ works as a catalyst to the viscosity at an increased suction parameter, and at η = 0, the temperature of the fluid within the respective boundary layer decreases, leading to the increased value of viscosity. Meanwhile, the curves corresponding to ε ∈ [0,2] in Figure 7a,b showed similar decreasing trends in viscosity. At ε = 0, there is no variation in the viscosity, and therefore η will always be 1, regardless of the values of other influential parameters. As ε began to increase, peak values started to increase from η = 0, before merging towards the unity, indicating an insignificant impact of η on the viscosity. Due to the dominance of ξ = 10, the viscosity reached the unity at a lower η than that of ξ = 1, which was also anticipated.

3.6.2. Effect on Velocity Distribution (f′)

In this part of the study, the changes in the non-dimensional velocity profiles were investigated by varying η and ε ∈ [0,2], at Pr = 0.1.

It was anticipated that, in the absence of ε, the velocity distribution would keep increasing as η increased and quickly reached the peak value, before η started to dominate, which would lead to a rapid decrease in the velocity (f′). As a consequence, at one specific value of η, the fluid would not exhibit any rapid mobility within the boundary. However, the local maximum of velocity would occur at a later stage of η if ε kept increasing. Furthermore, the static condition of fluid was also expected to be quicker for a higher ξ (in this case, ξ = 10), as a higher suction parameter demobilizes the fluid motion to a greater extent. Figure 8 illustrates the influence of the aforementioned parameters on the velocity distributions. As per Figure 8, as ε keeps increasing, the curves corresponding to the velocity exhibited peak values and later started to decrease significantly, heading to the static state as η increased.

At ε = 0 and ξ = 1 (Figure 8a), the velocity reached its maximum value of 0.955 at η = 1, whereas it was 0.85 at ε = 0.5. The lowest peak was recorded at the maximum ε values considered in this section (ε = 2), and due to the dominance of η and ε, the fluid had the lowest peak. In any case, after reaching the peak, velocity started to decrease significantly before reaching the completely static state at η = 10. On the other hand, at ξ = 10 (Figure 8b), similar trends were observed for the velocity at different ε and η, but due to the 10-time-increasing high suction parametric value, the peak velocity values were significantly lower, which indicates the superiority of the high ξ values. It is also expected that further increases in ξ will lead to further reduced peak values, and, at one stage, the fluid will remain static.

3.6.3. Changes in Temperature Distribution

Studying the changes in the temperature was crucial to the current research, since it provided more information on the effect of ε and η. To maintain consistency in the presentation, the changes in the non-dimensional temperature were studied by varying η and ε ∈ [0,2], at Pr = 0.1.

Figure 9 reports the consistent decreasing trend of temperature as η increased. However, the peak θ was recorded at the highest ε (in this case ε = 2) and lowest pseudo-similarity variable (η = 0), where ξ was unity. This behavior could be well-explained in the light of a conceptual understanding of the effect of the viscosity-variation parameter on the temperature. In the absence of η and ε, the temperature exhibited the lowest peak value of 2.9, whereas at ε = 2 the value of dimensionless temperature was found to be 3.35. The differences in the θ outline the impact of the viscosity-variation parameter on the temperature. However, as η started to have non-zero input, the temperature kept on decreasing rapidly and it was close to 0 as η reached maximum input values corresponding to the case studies (η = 10 for Figure 9a, η = 6 for Figure 9b). The lowest temperature values could also be linked with the velocity profiles obtained in Figure 8, which showed that in similar η values, the velocity reached 0, thus almost demobilizing the fluid within the boundary layer. Since there was no velocity or marginal velocity recorded at η = 10 in Figure 8a and η = 6 in Figure 8b, this backs up the concept that temperature will be close to 0 or equal to 0 as the fluid remains static. In addition, the influential impact of the suction parameter ξ largely reduced the temperature of the fluid, and the impact of ε was quite marginal, indicating the dominance of ξ (Figure 9b). While the curves seemed to overlap in Figure 9b, the magnified frame demonstrated the marginal differences at the narrowed axes.

3.7. Impact of ε, η, and ξ Fluid Characteristics at Variable Pr

3.7.1. Impact on Viscosity Distribution

In general, the Pr number distinguishes different types of fluids. For example, Pr = 0.01 refers to sodium, Pr between 0.02 to 0.03 refers to mercury. The purpose of different Pr tests in this study is to explore the uniformity of the numerical solution.

Figure 10 shows the changes in the non-dimensional viscosity values as a function of Pr. The value of ε was kept constant at 1, whereas two different suction parameters were appended in the numerical simulations. In addition, η ∈ [0,10] was also considered. It was observed that the fluid corresponding to Pr = 0.01 exhibited the highest peak value when η = 0, indicating the absence of η or the fact that the low viscosity-variation parameter increased the viscosity value. As Pr got increased, the peak values got reduced as well. While the local maximum viscosity of the fluid of Pr = 0.01 was found to be 9.45 at η = 0, the fluid corresponding to Pr = 0.1 exhibited a local maximum viscosity of 4.2. However, as η kept on increasing, the viscosity kept decreasing due to the dominance of η on the viscosity values. However, these phenomena were obtained at a low suction parametric value ξ = 1, as presented in Figure 10a. As ξ increased to 5, the peak values of each fluid at different Pr exhibited a rapid increase in the local maxima. The fluid corresponding to Pr = 0.01 demonstrated a peak viscosity value of approximately 40, which was 9.45 when ξ = 1 was considered. This 323.28% increase in the value could be attributed to the increase in ξ, which eventually increased the viscosity value in the absence of η. The rest of the fluids with further increased Pr also demonstrated similar characteristics, as depicted in Figure 10b.

3.7.2. Changes in the Velocity Distribution

Velocity distribution corresponding to different Pr implies the influence of Pr on different types of fluids, as other parameters, such as η and ξ, are varied. The viscosity-variation parameter ε was fixed at 1. In general, the patterns of velocity distribution had similarities with those obtained in Figure 11, which also indicates consistency in the current approach. As η started to increase from 0, the velocity started to increase rapidly before hitting the peak, followed by a sharp decrease in a fluid motion. However, within the η ∈ [0,10] boundary, only fluids corresponding to Pr = 0.01, 0.02 remained mobile, although the values of f′ were significantly lower at η = 10, as shown in Figure 11a. On the other hand, the rest of the fluids headed towards the static condition within the boundary layer at η = 10, thus establishing the preeminence of η over the fluid mobility. A similar behavior was observed at an increased ξ = 5 (Figure 8b), but the influence on f′ was not robust. As ξ increased to 5 from 1 in Figure 11b, the local maxima decreased. While at ξ = 1, Pr = 0.02, the local maximum was recorded to be approximately 1.4, and the value was found to be around 0.85 when ξ = 5 was considered. It should be mentioned here that while ξ attenuated the peak values of the velocity, it also influenced the local maxima of fluids with lower Pr values to occur at higher η. For example, fluid corresponding to Pr = 0.02 had a peak at η ≈ 2 at ξ = 1 (Figure 11a), whereas it was found to be at η ≈ 3 when ξ = 5 (Figure 11b).

3.7.3. Effect on Temperature Distribution

In the final part of sensitivity analyses, the effect of dimensionless parameters η and ξ has been investigated on the temperature distribution as Pr varied from 0.01 to 0.1 under no viscosity variation, i.e., ε was constant 1 throughout this analysis. In general, it was expected that η would have more influence on the temperature distribution with varied Pr numbers, as fluids’ temperature would likely to be more influenced by Pr values, as it would characterize the type of fluid, and η would be able to dominate, since the viscosity-variation parameter was kept constant. The suction parameter ξ was expected to have a marginal influence on the temperature distribution of fluids with low Pr (in this case 0.01 and 0.02) and a pronounced effect on the rest. In any case, the influence of ξ will not be significantly large on temperature distribution.

Figure 12 supports the aforementioned claims. For instance, the fluid corresponding to Pr = 0.01 had the highest peak temperature value, regardless of the ξ values. However, the peak was only visible at η = 0 in both cases (ξ = 1 or ξ = 5). As η increased gradually, the temperature value θ kept decreasing, indicating a reduction in fluid temperature within the boundary layer. However, fluids with Pr ≥ 0.05 exhibited a θ close to 0 as η = 10 was considered in the input. This behavior was also observed in both ξ values. At Pr = 0.01, θ marginally decreased from 7.9 (Figure 12a) to 7.6 (Figure 12b) as ξ increased from 1 to 5 when η was non-existent in the computation. A similar attribute was noticed at Pr = 0.02 as well. However, at Pr = 0.1, the plummeting rate was significantly larger. While θ = 3.6 was observed for P r= 0.1 under no η condition (Figure 12a), the value was found to be decreased to 1.8 when ξ = 5 was considered. This decrease could be attributed to the thermal properties of the fluid. Increased Pr values imply the fluids which have low or decreasing thermal conductivity. Therefore, as the Pr increased, the fluid was not able to conduct heat significantly, which eventually decreases the thickness of the boundary layer.

4. Conclusions

The present study discusses the effect of viscosity-variation in the natural convective flow inside a permeable vertical cone with a uniform surface heat flux. The numerical analyses performed in this study will provide elaborate explanations on the novelty of considering temperature-dependent viscosity to aid in different environmental and agricultural applications targeting uniform heat flux of a porous medium such as soil. Soil contains different types of fluids and may be represented in different geometries to explore all possible scenarios that occurred by either natural or human-made phenomena. Different types of fluids have been assigned by varying Pr numbers. Later, further investigations have been conducted in different parametric studies, which could be subcategorized into two major categories: (i) fixed Prandtl number (Pr) and (ii) variable Prandtl number (Pr). The numerical model was developed by transforming the governing boundary layer equations into a dimensionless form. In addition, the local non-similarity equations, obtained as a result of a non-linear system of partial differential equations’ reduction, were solved by three different solution techniques. The validations were conducted at both fixed Pr and variable Pr to assess the accuracy of the model. The model assessment was conducted by comparing the finite difference solution with the asymptotic solutions for large and small ξ values, and good agreements were established. The summary of the present study is mentioned in the following:

• Increasing the suction parameter (ξ) leads to decreasing shear stress (local skin-friction coefficient) and increasing the rate of heat transfer (local Nusselt number). The increasing and decreasing characteristics could be attributed to the temperature difference of the fluid within the boundary layer, which requires balancing the physical difference.

• As the viscosity-variation parameter (ε) increases, the local skin-friction coefficient decreases concurrently due to the effect of viscosity. On the other hand, increasing the rate of heat transfer as a function of ε does not remain consistent due to the effect of ξ. A small ξ rate of heat transfer decreases as ε increases, but the opposite behavior is observed with a large ξ, which indicates the superiority of the suction parameter over viscosity.

• If Pr increases, the local skin-friction coefficient decreases, and the rate of heat transfer increases rapidly due to the changes in the physical characteristics of the fluid and its viscosity.

• At ε ≠ 0, as η increases, viscosity decreases rapidly and heads towards 0 due to the dominance of η. However, at fixed ε and variable Pr, the curves corresponding to viscosity values exhibit the lowest value at a much later stage, indicating the high viscous characteristics of fluids with low Pr number.

• At any value of ε, at η ∈ [0,1], velocity increases the maximum at one point and then sharply plummets towards static condition at η > 6. However, a large suction parameter (ξ = 10) significantly lowers the peak values regardless of ε or types of fluids.

• An increased value of ε leads to the highest local maximum of temperature distribution in the absence of η and at a fixed Pr number. However, at η ≠ 0, the temperature starts to decrease gradually. The large suction parameter (ξ = 10) also suggestively lowers the peak values regardless of ε or types of fluids. Meanwhile, at the variable Pr number, the local maxima of temperature distribution get marginally affected at a low Pr number (<0.05).

Author Contributions

Conceptualization, M.F.H. and M.M.M.; methodology, M.F.H. and M.M.M.; software, M.M.M.; validation, M.F.H. and M.M.M.; formal analysis, M.F.H.; investigation, M.F.H. and M.M.M.; resources, M.M.M., M.K., S.S.; data curation, M.M.M.; writing—original draft preparation, M.F.H. and M.M.M.; writing—review and editing, M.F.H., M.M.M., M.K. and S.S.; visualization, M.F.H. and M.M.M.; supervision, M.M.M.; project administration, M.F.H., M.M.M., M.K. and S.S.; funding acquisition, M.M.M. All authors have read and agreed to the published version of the manuscript.

Funding

The second author thanks the North South University (NSU) for the financial support as a Faculty Research Grant (CTRG-21-SEPS-12). The second author also acknowledges the Ministry of Science and Technology (MOST) of the government of Bangladesh for providing the financial support for buying a desktop computer (Grant No.:474-EAS).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data availability on request.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Nomenclature

C	Specific heat
Cf	Skin-friction
f	Dimensionless stream function
Gr	Grashof number
k	Thermal conductivity
Nu	Nusselt number
Pr	Prandtl number
q	Heat flux
r	Radius of the cone
u	Velocity component in the x-direction
v	Velocity component in the y-direction
x	Coordinate along a cone ray
y	Coordinate normal to cone surface
g	Acceleration due to gravity
T	Temperature
V	Transpiration velocity
Greek Symbols
β	Thermal expansion coefficient
$\mathrm{ϵ}$	Viscosity-variation parameter
γ	Cone apex half-angle
$\mathsf{\Sigma }$	Summation
$\mathsf{\theta }$	Dimensionless temperature function
η	Pseudo-similarity variable
ξ	Dimensionless transpiration/suction parameter
$\mathrm{\rho }$	Fluid density
µ	Viscosity of the fluid
$\mathsf{\psi }$	Stream function
Subscripts
fil	Film temperature
i	Sequence of term
p	Constant pressure or isobaric
x	Differentiation with respect to x
w	Surface
∞	Ambient state
Superscript
‘	Differentiation with respect to $\tilde{\eta }$
i	Number of iterations

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. Physical geometry and coordinate system of the vertical permeable circular cone considered in this study.

Figure 2. Comparison of the finite difference solution with asymptotic solutions for the small and large suction parameter ξ in terms of the CfxGrx1/5 while Pr = 0.7 for two different values of ε, (a) ε = 0, and (b) ε = 2.0.

Figure 3. Comparison of the finite difference solution with asymptotic solutions for the small and large suction parameter ξ in terms of the NuxGrx−1/5 while Pr = 0.7 for two different values of ε namely, (a) ε = 0.0, and (b) ε = 2.0.

Figure 4. Comparison of the finite difference solution with asymptotic solutions for small and large suction parameters ξ in terms of (a) the local skin-friction coefficient (CfxGrx1/5) and (b) the local Nusselt number (NuxGrx−1/5) while ε = 0.5 and for two different values of Pr.

Figure 5. Streamlines for the two different viscosity-variation parameters ε at (a) Pr = 0.05 (b) Pr = 0.1. Here, the solid lines stand for ε = 0.0 and the dashed lines for ε = 0.5. The values on each curve represent the flow rate (ψ).

Figure 6. Isotherms for the two different viscosity-variation parameters ε at (a) Pr = 0.05 (b) Pr = 0.1. Here, the solid lines stand for ε = 0.0 and the dashed lines for ε = 0.5.

Figure 7. Viscosity distribution profiles for different viscosity-variation parameters ε at (a) suction parameter ξ = 1 and (b) suction parameter ξ = 10 while Pr = 0.1.

Figure 8. Changes in velocity distribution profiles for different viscosity-variation parameters ε at (a) suction parameter ξ = 1 and (b) suction parameter ξ = 10 while Pr = 0.1.

Figure 9. Temperature distribution for different viscosity-variation parameters ε at (a) suction parameter ξ = 1 and (b) suction parameter ξ = 10 while Pr = 0.1.

Figure 10. Viscosity distribution for the different Prandtl number (Pr): at (a) suction parameter ξ = 1 and (b) suction parameter ξ = 5 while ε = 1.0.

Figure 11. Velocity distribution for the different Prandtl number (Pr): at (a) suction parameter ξ = 1 and (b) suction parameter ξ = 5 while ε = 1.0.

Figure 12. Temperature distribution for the different Prandtl number (Pr): at (a) suction parameter ξ = 1 and (b) suction parameter ξ = 5 while ε = 1.0.

Table 1

Comparison of the finite difference solution with asymptotic solutions for the small and large suction parameter ξ in terms of the C_fxGr_x^1/5 and Nu_x Gr_x^−1/5 while Pr = 0.2 for two different values of ε.

ε = 0.0			ε = 1.0		ε = 0.0		ε = 1.0
Cf x Gr x 1/5					Nu x Gr x −1/5
ξ	Finite Diff. For All ξ	Small & Large ξ	Finite Diff. For All ξ	Small & Large ξ	Finite Diff. for All ξ	Small & Large ξ	Finite Diff. for All ξ	Small & Large ξ
0.0	1.9905	1.9929	1.9905	1.9929	0.3821	0.3821	0.3821	0.3821
0.02	1.9919	1.9958	1.9319	1.9354	0.3838	0.3836	0.3823	0.3822
0.04	1.9939	1.9987	1.8766	1.8823	0.3855	0.3852	0.3826	0.3824
0.06	1.9969	2.0014	1.8256	1.8337	0.3871	0.3867	0.3828	0.3826
0.08	1.9999	2.0042	1.7778	1.7895	0.3887	0.3883	0.3829	0.3827
0.10	2.0028	2.0069	1.7330	1.7499	0.3903	0.3898	0.3833	0.3831
0.14	2.0028	2.0122	1.6909	1.6842	0.3934	0.3929	0.3836	0.3833
3.0	1.6938		0.3584	-	0.6893		0.6591	-
5.0	0.9349	0.9388	0.1657	0.1620	1.0185	1.0169	1.0119	1.0095
6.0	0.6741	0.6774	0.1189	0.1145	1.2099	1.2081	1.2077	1.2046
7.0	0.5020	0.5044	0.0898	0.0846	1.4064	1.4044	1.4069	1.4025
8.0	0.3064	0.3284	0.0709	0.0649	1.8044	1.6026	1.6082	1.6014
10.0	0.2488	0.2495	0.0496	0.0416	2.0042	2.0011	2.0163	2.0006

References

1. Mojid, M.A.; Cho, H. Estimating the fully developed diffuse double layer thickness from the bulk electrical conductivity in clay. Appl. Clay Sci.; 2006; 33, pp. 278-286. [DOI: https://dx.doi.org/10.1016/j.clay.2006.06.002]

2. Zhang, N.; Wang, Z. Review of soil thermal conductivity and predictive models. Int. J. Therm. Sci.; 2017; 117, pp. 172-183. [DOI: https://dx.doi.org/10.1016/j.ijthermalsci.2017.03.013]

3. Hasan, M.F.; Abuel-Naga, H.; Leong, E.-C. A modified series-parallel electrical resistivity model of saturated sand/clay mixture. Eng. Geol.; 2021; 290, 106193. [DOI: https://dx.doi.org/10.1016/j.enggeo.2021.106193]

4. Himika, T.A.; Hassan, S.; Hasan, M.F.; Molla, M.M. Lattice Boltzmann Simulation of MHD Rayleigh–Bénard Convection in Porous Media. Arab. J. Sci. Eng.; 2020; 45, pp. 9527-9547. [DOI: https://dx.doi.org/10.1007/s13369-020-04812-z]

5. Simo-Tagne, M.; Ndukwu, M.C.; Zoulalian, A.; Bennamoun, L.; Kifani-Sahban, F.; Rogaume, Y. Numerical analysis and validation of a natural convection mix-mode solar dryer for drying red chilli under variable conditions. Ren. Energy; 2020; 151, pp. 659-673. [DOI: https://dx.doi.org/10.1016/j.renene.2019.11.055]

6. Pangavhane, D.R.; Sawhney, R.L.; Sarsavadia, P.N. Design, development and performance testing of a new natural convection solar dryer. Energy; 2002; 27, pp. 579-590. [DOI: https://dx.doi.org/10.1016/S0360-5442(02)00005-1]

7. Kumar, D.; Mahanta, P.; Kalita, P. Energy and exergy analysis of a natural convection dryer with and without sensible heat storage medium. J. Energy Stor.; 2020; 29, 101481. [DOI: https://dx.doi.org/10.1016/j.est.2020.101481]

8. Calautit, J.K.; O’Connor, D.; Sofotasiou, P.; Hughes, B.R. CFD simulation and optimisation of a low energy ventilation and cooling system. Computation; 2014; 3, pp. 128-149. [DOI: https://dx.doi.org/10.3390/computation3020128]

9. Hireche, Z.; Himrane, N.; Nasseri, L.; Hamrioui, Y.; Ameziani, D.E. Analysis of Thermal Performances in a Ventilated Room Using LBM-MRT: Effect of a Porous Separation. Computation; 2022; 10, 4. [DOI: https://dx.doi.org/10.3390/computation10010004]

10. Himika, T.A.; Hasan, M.F.; Molla, M.M. Lattice Boltzmann simulation of airflow and mixed convection in a general ward of hospital. J. Comp. Eng.; 2016; 2016, 5405939. [DOI: https://dx.doi.org/10.1155/2016/5405939]

11. Hasan, M.F.; Himika, T.A.; Molla, M.M. Lattice Boltzmann simulation of airflow and heat transfer in a model ward of a hospital. J. Therm. Sci. Eng. Appl.; 2017; 9, 011011. [DOI: https://dx.doi.org/10.1115/1.4034817]

12. Chen, L.; Zhao, J.; Yuan, Y.; Lee, J. Numerical Study on the Thermal Field and Heat Transfer Characteristics of a Hexagonal-Close-Packed Pebble Bed. Computation; 2022; 10, 1. [DOI: https://dx.doi.org/10.3390/computation10010001]

13. Molla, M.M.; Nag, P.; Thohura, S.; Khan, A. A graphics process unit-based multiple-relaxation-time lattice Boltzmann simulation of non-Newtonian fluid flows in a backward facing step. Computation; 2020; 8, 83. [DOI: https://dx.doi.org/10.3390/computation8030083]

14. Ahmed, N.S.; Badruddin, I.A.; Zainal, Z.A.; Khaleed, H.M.T.; Kanesan, J. Heat transfer in a conical cylinder with porous medium. Int. J. Heat Mass Transf.; 2009; 52, pp. 3070-3078. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2008.12.030]

15. Khan, T.Y.; Badruddin, I.A.; Quadir, G.A. Heat transfer in a conical porous cylinder with partial heating. IOP Conference Series: Materials Science and Engineering; IOP Publishing: Bristol, UK, 2016; 012211.

16. Hanif, H.; Khan, I.; Shafie, S. Heat transfer exaggeration and entropy analysis in magneto-hybrid nanofluid flow over a vertical cone: A numerical study. J. Therm. Analys. Cal.; 2020; 141, pp. 2001-2017. [DOI: https://dx.doi.org/10.1007/s10973-020-09256-z]

17. Rao, M.; Gangadhar, K.; Chamkha, A.J.; Surekha, P. Bioconvection in a convectional nanofluid flow containing gyrotactic microorganisms over an isothermal vertical cone embedded in a porous surface with chemical reactive species. Arab. J. Sci. Eng.; 2021; 46, pp. 2493-2503. [DOI: https://dx.doi.org/10.1007/s13369-020-05132-y]

18. Rahman, A.; Molla, M.M.; Sarker, M.M.A. Natural convection flow along the vertical wavy cone in case of uniform surface heat flux where viscosity is an exponential function of temperature. Int. Comm. Heat Mass Transf.; 2011; 38, pp. 774-780. [DOI: https://dx.doi.org/10.1016/j.icheatmasstransfer.2011.03.021]

19. Merk, H.J.; Prins, J.A. Thermal convection in laminary boundary layers I. Appl. Sci. Res. Sect. A; 1953; 4, pp. 11-24. [DOI: https://dx.doi.org/10.1007/BF03184660]

20. Merk, H.J.; Prins, J.A. Thermal convection in laminar boundary layers II. Appl. Sci. Res. Sect. A; 1954; 4, pp. 195-206. [DOI: https://dx.doi.org/10.1007/BF03184951]

21. Braun, W.H.; Ostrach, S.; Heighway, J.E. Free-convection similarity flows about two dimensional and axisymmetric bodies with closed lower ends. Int. J. Heat Mass Transf.; 1961; 2, pp. 121-135. [DOI: https://dx.doi.org/10.1016/0017-9310(61)90020-5]

22. Hering, R.G.; Grosh, R.J. Laminar free convection from a non-isothermal cone. Int. J. Heat Mass Transf.; 1962; 5, pp. 1059-1068. [DOI: https://dx.doi.org/10.1016/0017-9310(62)90059-5]

23. Hossain, M.A.; Paul, S.C. Free convection from a vertical permeable circular cone with non-uniform surface heat flux. Heat Mass Transf.; 2001; 37, pp. 167-173. [DOI: https://dx.doi.org/10.1007/s002310000129]

24. Hossain, M.A.; Paul, S.C. Free convection from a vertical permeable circular cone with non-uniform surface temperature. Acta Mech.; 2001; 151, pp. 103-114. [DOI: https://dx.doi.org/10.1007/BF01272528]

25. Hossain, M.A.; Paul, S.C.; Mandal, A.C. Natural convection flow along a vertical circular cone with uniform surface temperature and surface heat flux in a thermally stratified medium. Int. J. Num. Methods Heat Fluid Flow; 2002; 12, pp. 290-305. [DOI: https://dx.doi.org/10.1108/09615530210422965]

26. Cebeci, T.; Bradshaw, P. Physical and Computational Aspects of Convective Heat Transfer; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012.

27. Gary, J.; Kassoy, D.R.; Tadjeran, H.; Zebib, A. The effects of significant viscosity variation on convective heat transport in water-saturated porous media. J. Fluid Mech.; 1982; 117, pp. 233-249. [DOI: https://dx.doi.org/10.1017/S0022112082001608]

28. Mehta, K.N.; Sood, S. Transient free convection flow with temperature dependent viscosity in a fluid saturated porous medium. Int. J. Eng. Sci.; 1992; 30, pp. 1083-1087. [DOI: https://dx.doi.org/10.1016/0020-7225(92)90032-C]

29. Ling, J.X.; Dybbs, A. Forced Convection over a Flat Plate Submersed in a Porous Medium: Variable Viscosity Case; American Society of Mechanical Engineers: New York, NY, USA, 1987.

30. Thohura, S.; Rahman, A.; Molla, M.M.; Sarker, M.M. Effects of temperature dependent thermal conductivity on natural convection flow along a vertical wavy cone with heat flux. Proc. Eng.; 2014; 90, pp. 497-503. [DOI: https://dx.doi.org/10.1016/j.proeng.2014.11.763]

31. Khan, M.; Salahuddin, T.; Altanji, M. A viscously dissipated Blasisus boundary layer flow with variable thermo-physical properties: An entropy generation study. Int. Comm. Heat Mass Transf.; 2022; 131, 105873. [DOI: https://dx.doi.org/10.1016/j.icheatmasstransfer.2021.105873]

32. Gladys, T.; Reddy, G.V.R. Contributions of variable viscosity and thermal conductivity on the dynamics of non-Newtonian nanofluids flow past an accelerating vertical plate. Partial Diff. Eq. Appl. Math.; 2022; 5, 100264. [DOI: https://dx.doi.org/10.1016/j.padiff.2022.100264]

33. Keller, H.B. Numerical methods in boundary-layer theory. Ann. Rev. Fluid Mech.; 1978; 10, pp. 417-433. [DOI: https://dx.doi.org/10.1146/annurev.fl.10.010178.002221]

34. Kamran, A.; Azhar, E.; Akmal, N.; Mehmood, Z.; Iqbal, Z. Finite Difference Approach for Critical Value Analysis to Describe Jeffery–Hamel Flow toward an Inclined Channel with Microrotations. Arab. J. Sci. Eng.; 2022; pp. 1-8. [DOI: https://dx.doi.org/10.1007/s13369-021-06532-4]

35. Reddy, Y.D.; Goud, B.S.; Chamkha, A.J.; Kumar, M.A. Influence of radiation and viscous dissipation on MHD heat transfer Casson nanofluid flow along a nonlinear stretching surface with chemical reaction. Heat Transfer.; 2022; pp. 1-17. [DOI: https://dx.doi.org/10.1002/htj.22460]

36. Vajravelu, K.; Kerehalli, V.P. Keller Box Method and Its Application; Walter de Gruyter GmbH & Co KG: De Gruyter, Germany, 2014; Volume 8.

37. Molla, M.M.; Hossain, M.A. Radiation effect on mixed convection laminar flow along a vertical wavy surface. Int. J. Therm. Sci.; 2007; 46, pp. 926-935. [DOI: https://dx.doi.org/10.1016/j.ijthermalsci.2006.10.010]

38. Charraudeau, J. Influence de gradients de proprietes physiques en convection forcee—Application au cas du tube. Int. J. Heat Mass Transf.; 1975; 18, pp. 87-95. [DOI: https://dx.doi.org/10.1016/0017-9310(75)90011-3]

39. Butcher, J.C. Implicit Runge-Kutta processes. Math. Comput.; 1964; 18, pp. 50-64. [DOI: https://dx.doi.org/10.1090/S0025-5718-1964-0159424-9]

40. Nachtsheim, P.R.; Swigert, P. Statisfaction of Asymptotic Boundary Conditions in Numerical Solution of Systems of Nonlinear Equations of Boundary-Layer Type. No. NASA-TN-D-3004. 1965; Available online: https://ntrs.nasa.gov/api/citations/19650026350/downloads/19650026350.pdf (accessed on 18 March 2022).

41. Akin, E.D. Object-Oriented Programming via Fortran 90/95; Cambridge University Press: Cambridge, UK, 2003; Volume 1.

42. OriginPro, V. OriginLab Corporation; Northampton, MA, USA. 2016; Available online: https://www.originlab.com (accessed on 18 March 2022).