A series solution for melting heat transfer characteristics of hybrid Casson fluid under thermal radiation

In the present paper, we focus on the melting heat transfer characteristics of Casson fluid involving thermal radiation and viscous dissipation. To this end, the governing partial differential equations (PDEs) are transformed into the ordinary differential equations (ODEs) via the similarity variables. Besides establishing a homotopy-based methodology and its optimization performed in MATHEMATICA package BVPh2.0, the present findings are compared and validated by those available results in the literature. It can be shown that regardless of the variable fluid properties, this methodology predicts the heat transfer rate with and without melting effect at any Prandtl number. Furthermore, it is seen that the velocity distribution is significantly affected by the melting parameter.


Introduction
T he main aim of fluid mechanics, as a multiform phenomenon, is to perform modelling of the liquids and gases in their common temperature and pressure ranges. Conservation of mass, momentum and energy can be taken as three basic laws for any engineering problem involving in the boundary layer theory. Recently, the complexity of such problems can be greatly simplified by using the approximate solution methodologies. In this way, Raees et al. [1] proposed a multiple solution for homogeneous-heterogeneous reactions on the Oberbeck-Boussinesq buoyancy-driven flow of dilute nanofluids. They derived the governing nonlinear equations based on the Buongiorno's mathematical model, and showed that the dilute nanofluid gives a higher stability of chemical reactions than a typical one. Ghiasi and Saleh [2] analyzed the Casson type fluid over an unsteady shrinking surface in the presence of Joule heating and magnetic inclination. They showed that the effect of viscous force in their model will be usually negligible for large values of the Hartmann number. Gaffar et al. [3] employed the Keller box finite difference method (FDM) to simulate the thermomechanical behaviour of viscoelastic Jeffery's fluid inside a non-parallel vertical channel so that an excellent consistency with those of Hossain and Paul [4] was captured. They also proved that the predicted behaviour is quite elastic with large Deborah numbers. Sayyed et al. [5] took the influence of permeability and slip boundary condition into account to study the geothermal aspects of flow over a wedge-shaped configuration. They combined the differential transformation method (DTM) with the Padé approximation, and showed that the suction in this case has the same effect as the injection on the permeability. Besthapu et al. [6] extended the double stratification (i.e., the stratified medium due to the temperature and nanoparticle concentration fields) to the laminar boundary layer flow along an exponentially stretching surface, which previously was reported by Khan and Pop [7]. Furthermore, they showed that their results are in good agreement with those obtained by Ishak [8].
Melting heat transfer can be briefly defined as follows: The melt liquid has not been diffused into the ambient liquid [9,10]. It has found many applications in engineering and physics, for example, to magma solidification, to semiconductor materials, thawing permafrost, etc. However, as we will now see, there exist only a few research studies concerning this issue. Hayat et al. [11] performed a melting heat transfer analysis of Powell-Eyring fluid with convective boundary conditions. Sheikholeslami and Rokni [12] used a similar idea applied to the forced convection flows considering the Buongiorno's mathematical model. They showed that the drag force distribution generated from the magnetic field affects the momentum boundary layer thickness.
Valipour et al. [13] investigated the heat and mass transfer characteristics of two-phase nanofluids induced by nonlinearly rotating plates carried out numerically in the standard MAPLE software. They found that the rotational velocity and temperature distributions enhance with an increase in the Reynolds number. Hayat et al. [14] employed the perturbation technique to analyze the melting heat transfer of peristaltic flows in the presence of viscous dissipation and Joule heating. They also analyzed the oscillatory heat transfer caused by the sinusoidal traveling wave in a channel.
It is noteworthy to mention that although some sort of solution methodologies has been employed to investigate the melting heat transfer characteristics to date [15][16][17][18][19][20][21][22][23][24][25][26], the homotopy-based methodology, due to its optimization performed in MATHEMATICA package BVPh2.0 [27], can be developed as an efficient procedure compared to any other one [11,28]. It is found that this methodology does not suffer from long run-time and can provide high accuracy estimates. The rest of the paper is organized as follows: Section 2 deals exclusively with the mathematical formulation. Section 3 outlines the homotopy-based methodology and its optimization. Section 4 focuses on the results and discussion. The concluding remarks are summarized in Section 5.

Governing equations
Let us denote the constitutive equation for an isotropic and incompressible flow of Casson fluid as follows [29]: where τ ij is the Cauchy stress tensor, µ B is the plastic dynamic viscosity, p y = µ B √ 2π λ is the fluid yield stress, e ij is the (i, j)th component(s) of the deformation rate, π(= e ij e ij ) is the product of the deformation rate component(s), π c is the critical value of π and λ is termed as the Casson fluid parameter. In the case of π > π c , Equation (1) reduces to [30], where µ n f is the dynamic viscosity. According to the above formula, the kinematic viscosity can be defined as where ρ n f is the density. The velocity and temperature fields for a two-dimensional laminar flow in the Cartesian coordinates take the form where u and v are the velocity components along the x and y directions, respectively, and T is the temperature. By eliminating the pressure terms, one can obtain a set of governing PDEs in the following form with the boundary conditions and [10], where c p is the specific heat at constant pressure, k n f is the thermal conductivity, q r is the radiation heat flux, U w is the velocity at the wall, a > 0 is an arbitrary constant, T m is the temperature of the melting surface, T ∞ is the ambient temperature, δ is the latent heat of the fluid, c s is the heat capacity of the solid surface and T 0 is the temperature of the solid surface. According to the basic hypothesis of Rosseland approximation [31], the radiation heat flux presented in Equation (5) becomes where σ * and k * are the Stefan-Boltzmann constant and the mean absorption coefficient, respectively. Assuming that the fluid-phase temperature discrepancy within the flow is sufficiently small [32], T 4 can be expanded in the Taylor series as follows By substituting Equation (9) into Equation (8), one would get The similarity variables, which convert Equation (5) into the ODEs are given by where η is the similarity parameter, ψ is the stream function which satisfies the continuity equation (i.e. u = ∂ψ ∂y and v = − ∂ψ ∂x ), f is the similarity function and θ is the dimensionless temperature.
By substituting Equations (10) and (11) into Equation (5), the following set of governing ODEs can be expressed as with the boundary conditions where Pr = Here, the physical quantities of interest i.e., the skin friction coefficient and local Nusselt number can be defined as, where By substituting Equation (11) into Equations (14) and (15) we see that where Re x = xU w v n f is the local Reynolds number.

Solution Methodology
As it was mentioned earlier, finding an exact closed-form solution for the third-order nonlinear ODEs given in Equation (12) is not generally possible. According to the basic concept of homotopy-based methodology, the governing ODEs should be discretized to an infinite sequence of linear functions employing a typical series expansion such as Taylor, Maclaurin, etc. [33,34]. It is worth mentioning that choosing an appropriate auxiliary parameter provides a convenient way to adjust the convergence of series solution [33,35]. Furthermore, unlike perturbation techniques, this methodology can be applied to any type of ODEs and associated boundary conditions consisting of small or large parameters [36]. In this way, we choose the initial guesses and auxiliary linear operators as follows, with the properties where C 1 − C 5 are the integration constants. By introducing q ∈ [0, 1] as an embedding parameter, the zeroth order deformation equations are constructed as [33], where f and θ are the so-called auxiliary parameters, and N f and N θ are the nonlinear differential operators which can be expressed as follows, with the boundary conditions, As q increases from 0 to 1,f (η; q) andθ(η; q) vary from the initial guesses f 0 (η) and θ 0 (η) to the exact solutions f (η) and θ(η), respectively. Expandingf (η; q) andθ(η; q) in the Taylor series with respect to q gives, where f n and θ n are the nth-order deformation derivatives, that is, It is to be noted here that the solution of Equation (20) depends only on the auxiliary linear operators, initial guesses and auxiliary parameters [33]. Accordingly, if these quantities are properly selected so the series Equation (23) converges at q = 1 in the following form By Differentiating Equation (20) n times with respect to q, then setting q = 0 and dividing them by n!, one can obtain the nth-order deformation equations as follows where χ n = 0, n ≤ 1, To complete the solution procedure, it is sufficient to set the given boundary conditions given in Equation (22) equal to zero, obtain the arbitrary constants C 1 − C 5 , and thereby derive the general solutions of Equation (26) as where f * n (η) and θ * n (η) are the particular solutions. Hence, after solving the nth-order deformation Equation (26), the nth−order approximate solutions are given by Here, in the context of optimization, the following squared residual errors are defined as [37], where j = 20 and δη = 0.5 [38,39].

Results and Dicussion
To find the optimal values of auxiliary parameters and number of iterations, a convergence study performed in MATHEMATICA package BVPh2.0 [27] is provided here. To this end, the governing physical parameters involved in Equations (12) and (13) are given in Table 1. According to the results of Table 2, the optimal values of auxiliary parameters decrease when the order of approximation is increased. Hence, in what follows, only f = −0.64492 and θ = −0.85516 are employed to ensure the high accuracy estimates. Furthermore, Table 3 represents the values of squared residual errors at different orders. It is observed form this table that the series solution will thoroughly be converged if the number of iterations is set equal to p = 20. As it can be seen from Table 4, increasing the value of melting parameter decreases the skin friction coefficient. This is due to the fact that by increasing the melting parameter, a temperature discrepancy between the solid and melting surfaces generates which is led to the molecular motion of the particles. Furthermore, Table 4 provides a comparison between the present findings and those reported by Krishnamurthy et al. [15].  2.322 × 10 −9 14 6.881 × 10 −9 9.170 × 10 −10 16 3.230 × 10 −9 6.989 × 10 −10 18 1.014 × 10 −9 4.721 × 10 −10 20 9.175 × 10 −10 1.990 × 10 −10 Based on the results shown in this table, the present findings agree well with those of Krishnamurthy et al. [15] in all cases; because the relative error between them does not exceed 0.08%. Table 5 investigates the variation of heat transfer rate versus the Prandtl number for some different melting parameters. According to this table, the heat transfer rate increases monotonically as the Prandtl number is increased from 0.07 to 7. This is because of the correspondence between the thermal conductivity and diffused heat energy [40]. In addition, as it can be observed from Table 5, the heat transfer rate decreases as the melting parameter increases which is only due to the existing energy dissipation. However, at very small melting parameters there is a low-temperature enhancement [19].
It is to be mentioned here that the accuracy of the present findings without its melting effect is verified by comparison with the numerical observations reported by Khan and Pop [7]. Hence, in view of the results listed in Tables 4 and 5, one can emphasize that employing the homotopy-based methodology and its optimization is sufficient to accelerate the convergence on any large interval and thereby ensure the accuracy and efficiency of the series solution. Here, for the sake of brevity, only the velocity and temperature distributions with and without the melting effect are provided in Table 6.
As it was mentioned in Section 2, the fluid yielding behavior is greatly affected by the shear and thermal history. Furthermore, the Casson fluid parameter plays a crucial role in the fluid yielding threshold of such rheological systems [41]. The reason is that the ratio of viscous force to yield stress increases with an increase in the Casson fluid parameter. It is worth noting that this increase of the Casson fluid parameter leads to the free flow movement with a specific pressure gradient. This issue is clearly illustrated in Figure 1.
According to Figure 1, one can observe that the skin friction coefficient decreases with an increase in the Casson fluid parameter; however, Mabood et al. [22] have provided some evidence to illustrate that this coefficient in some cases is relatively insensitive to the Casson fluid parameter. The discrepancy between the present paper and that reported by Mabood et al. [22] is due to the presence of temperature-dependent dynamic viscosity and its roughly exponential function. Table 4. Comparison between the present findings and those reported by Krishnamurthy et al. [15] for the values of -Re 1/2 x C f with λ → ∞ and N = 0.   Table 6. Velocity and temperature distributions with parameters presented in Table 1. Based on the results presented in Figure 2, neglecting the melting effect decreases slightly the local Nusselt number. Furthermore, increasing the value of radiation parameter decreases the effect of melting on increasing the local Nusselt number. It is due to a reduction in the mean absorption coefficient and consequently divergence of the radiation heat flux. Thus, it is desirable to consider this observation in some cases as the basis of controlling the thermal radiation from solid and melting surfaces.

Concluding remarks
The present paper aimed to provide an investigation of how the melting effect can be applied to a set of governing PDEs which govern a two-dimensional laminar flow of Casson fluid. In this way, the conservation of mass, momentum and energy, combined with the Rosseland approximation, were employed and reduced to a system of ODEs through the similarity variables. The homotopy-based methodology performed in MATHEMATICA package BVPh2.0 was developed because the series solution to be optimized depends significantly on choosing the appropriate auxiliary parameters. The accuracy of the present findings was also verified by comparison with the numerical results in the literature. Since the author's perspective was to introduce the idea of melting effect, the main inferences in this case can be summarized as follows: • In general, the velocity and temperature distributions are affected by the melting parameter more than that of the other pertinent parameters. • Accounting for the effect of Casson fluid as well as increasing the melting parameter decreases the skin friction coefficient. • Accounting for the effect of thermal conductivity as well as increasing the melting parameter decreases the heat transfer rate. • Accounting for the effect of thermal radiation as well as increasing the melting parameter increases the local Nusselt number.