Analytical Solution of Flow and Heat Transfer over a Permeable Stretching Wall in a Porous Medium

Abstract

Boundary layer flow through a porous medium over a stretching porous wall has seen solved with analytical solution. It has been considered two wall boundary conditions which are power-law distribution of either wall temperature or heat flux. These are general enough to cover the isothermal and isoflux cases. In addition to momentum, both first and second laws of thermodynamics analyses of the problem are investigated. The governing equations are transformed into a system of ordinary differential equations. The transformed ordinary equations are solved analytically using homotopy analysis method. A comprehensive parametric study is presented, and it is shown that the rate of heat transfer increases with Reynolds number, Prandtl number, and suction to the surface.

1. Introduction

Boundary layer flows over a stretching surface have application in engineering processes such as liquid composite molding, extrusion of plastic sheets, paper production, glass blowing, metal spinning, wire drawing, and hot rolling [1–3]. More importantly, the quality of the products, in the abovementioned processes, depends on the kinematics of stretching and the simultaneous heat and mass transfer rates during the fabrication process. Sakiadis [4, 5] and Crane [6] were the pioneers in the investigations of boundary layer flow over continuously moving surfaces that are quite different from the free stream flow over stationary flat plates. Pop and Na [7] studied free convection heat transfer of non-Newtonian fluids along a vertical wavy surface in a porous medium.

The flow field of a stretching surface with a power-law velocity variation was discussed by Banks [8]. Elbashbeshy [9] investigated heat transfer over a stretching surface with variable surface heat flux. Elbashbeshy and Bazid [10] analyzed heat transfer in a porous medium over a stretching surface with internal heat generation and suction or injection. This work was extended by Cortell [11] to include power-law temperature distribution. Steady flow and heat transfer of a viscous incompressible fluid flow through porous medium over a stretching sheet were studied by Sriramalu et al. [12]. Ali [13] investigated thermal boundary layer on a power-law stretched surface with suction or injection for uniform and variable surface temperatures. Recently, Pantokratoras [14] published analytical solution for velocity distribution inside a Darcy-Brinkman porous medium with a stretching boundary. Elbashbeshy [15] included thermal radiation effects in a stretching surface problem. Tamayol and Bahrami [16] understood that porous materials can be used to enhance the heat transfer rate from stretching surfaces to improve processes such as hot rolling and composite fabrication.

In view of the above, an analytical solution is developed in the present study to solve fluid flow, heat transfer in a porous medium over a porous plate with linear velocity, subjected to different power-law thermal boundary conditions. The analytical solution is successfully validated in comparison with numerical analysis. There are many effective methods for obtaining the solutions of nonlinear equation such as variational iteration method [17], Adomian method [18], and homotopy perturbation method [19, 20], and we use one of a powerful technique, namely, the homotopy analysis method (HAM) which was expressed by Liao [21–24]. This method has been successfully applied to solve many types of nonlinear problems [25–29] and provides us with great freedom to express solutions of a given nonlinear problem by means of different base functions. Secondly, unlike all previous analytic techniques, the homotopy analysis method always provides us with a family of solution expressions in the auxiliary parameter ℏ, even if a nonlinear problem has a unique solution. Thirdly, unlike perturbation techniques, the homotopy analysis method is independent of any small or large quantities. So, the homotopy analysis method can be applied no matter if governing equations and boundary/initial conditions of a given nonlinear problem contain small or large quantities or not.

Above all, there are no rigorous theories to direct us to choose the initial approximations, auxiliary linear operators, auxiliary functions, and auxiliary parameter ℏ. From the practical viewpoints, there are some fundamental rules such as the rule of solution expression, the rule of coefficient ergodicity, and the rule of solution existence, which play important roles within the homotopy analysis method. Unfortunately, the rule of solution expression implies such an assumption that we should have, more or less, some knowledge about a given nonlinear problem a prior. So, theoretically, this assumption impairs the homotopy analysis method, although we can always attempt some base functions even if a given nonlinear problem is completely new for us.

2. Governing Equation

Consider a steady, constant property, two-dimensional flow through a homogenous porous medium of permeability 𝐾, over a stretching surface with linear velocity distribution; that is, 𝑢_𝑤=𝑢₀𝑥/𝐿 (Figure 1). The transport properties of the medium can be assumed independent of temperature when the difference between wall and ambient temperatures is not significant [4]. The origin is kept fixed, while the wall is stretching, and the 𝑦-axis is perpendicular to the surface. The governing equations are [4, 5] 𝜕𝑢+𝜕𝑥𝜕𝜐𝜌𝑢𝜕𝑦=0,𝜕𝑢𝜕𝑥+𝜐𝜕𝑢𝜕𝑦=𝜇_eff𝜕²𝑢𝜕𝑦²−𝜇𝐾^′𝑢𝑢,𝜕𝑇𝜕𝑥+𝜐𝜕𝑇𝜕𝑦=𝛼_eff𝜕²𝑇𝜕𝑦², (1) where 𝑢and 𝜐 are velocity components in the 𝑥 and 𝑦 directions, respectively. 𝜇_eff is the effective viscosity which for simplicity in the present study is considered to be identical to the dynamic viscosity, 𝜇. The transport properties of the porous medium such as permeability depend on their microstructure and can be calculated either using existing correlations in the literature or through experimental measurements. 𝛼_eff is the effective thermal diffusivity of the medium. The hydrodynamic boundary conditions are 𝑢𝑥^∗,0=𝑢₀𝑥^∗,𝜐𝑥^∗,0=𝜐_𝑤,𝑢𝑥^∗,∞=0, (2) where 𝑥^∗=𝑥/𝐿 is the nondimensional 𝑥-coordinate and 𝐿 is the length of the porous plate. The following thermal boundary conditions are considered: 𝑇𝑥^∗,0=𝑇_∞+𝑇₀𝑥^∗𝑛,𝑇𝑥^∗,∞=𝑇_∞,−𝜅𝜕𝑇𝜕𝑦_(𝑥∗,0)=𝑞₀𝑥^∗𝑛,𝑇𝑥^∗,∞=𝑇_∞, (3) where 𝜅 is the effective thermal conductivity of the medium and is a function of thermal conductivities of the fluid and solid phases and the porous medium microstructure.

Figure 1 

Schematic diagram of problem.

Using dimensionless parameters, 𝑦𝜂=√𝐾,𝑢=𝑢₀𝑥^∗𝑓^′𝜂,𝑢𝜐=−⁰𝐿√𝐾𝑓𝜂,𝜓=𝑢₀𝑥^∗√𝐾𝑓𝜂, (4) where 𝑓^′ is 𝑑𝑓/𝑑𝜂.

The transformed nonlinear ordinary differential equations are 𝑓^′′′+Re𝑓𝑓^′′−𝑓^′2−𝑓^′𝜃=0,_′′+Re⋅Pr𝑓𝜃^′−𝑛⋅Re⋅Pr𝑓^′𝜃=0, (5) where Re=𝜌𝑢₀𝐾/𝐿𝜇 is the Reynolds number. Equation (5) should be solved subject to the following boundary conditions: 𝑓^′(0)(0)𝜐=1,𝑓=−^𝑤𝐿𝑢₀√𝐾=𝑓_𝑤,𝑓^′(∞)=0, (6)𝜃(0)(∞)𝜃=1,𝜃=0Power-lawtemperature,_′(0)=−1,𝜃^′(∞)=0Power-lawheatflux, (7a)(7b)where 𝑓_𝑤 is the injection parameter. Positive/negative values of 𝑓_𝑤 show suction/injection into/from the porous surface, respectively. The wall shear stress term can be calculated as 𝜏_𝐿=−𝜇𝑢₀𝐿𝑓^′′(0)2√𝐾. (8)

For power-law fluid wit constant temperature and heat flux boundary conditions, respectively. Employing the definition of convective heat transfer coefficient, the local Nusselt numbers become Nu_𝑥=ℎ𝑥𝑘=−𝜃^′(0)𝑥√𝐾𝑞,Power-lawtemperature,_𝑤𝑥𝑘𝑇_𝑤−𝑇_∞=𝑥𝜃(0)√𝐾,Power-lawheatflux. (9)

3. Solution of Problem by Homotopy Analysis Method (HAM)

As mentioned by Liao, a solution may be expressed with different base functions, among which some converge to the exact solution of the problem faster than others. Such base functions are obviously better suited for the final solution to be expressed in terms of. Noting these facts, we have decided to express 𝑔(𝜂) by a set of base functions of the following form: 𝑓𝜂=^∞∞_𝑚=0^∞_𝑛=0𝑘=0𝑏^𝑘_𝑚,𝑛𝜂^𝑘,𝜃𝜂=exp−𝑛𝜂_{∞∞𝑚=0∞𝑛=0𝑘=0}𝑏^𝑘_𝑚,𝑛𝜂^𝑘exp−𝑛𝜂. (10) The rule of solution expression provides us with a starting point. It is under the rule of solution expression that initial approximations, auxiliary linear operators, and the auxiliary functions are determined. So, according to the rule of solution expression, we choose the initial guess and auxiliary linear operator in the following form: 𝑓₀𝜂=𝑓_𝑤+1−exp−𝜂,𝜃₀𝜂=exp−𝜂. (11) We choose linear operator as follows: ℒ₁𝑓=𝑓^′′′+𝑓^′′,ℒ₂(𝜃)=𝜃^′′+𝜃^′,ℒ₁𝑐₁+𝑐₂𝜂+𝑐₃ℒexp−𝜂=0,₂𝑐₄+𝑐₅exp−𝜂=0, (12) where 𝑐_𝑖(𝑖=1--5) are constants. Let 𝑃∈[0,1] denote the embedding parameter, and let ℏ indicate nonzero auxiliary parameters. Then, we construct the following equations.

3.1. Consider Zeroth–Order Deformation Equations

(1−𝑃)ℒ₁𝑓𝜂;𝑝−𝑓₀𝜂=𝑝ℏ₁𝐻𝜂𝑁₁𝑓,𝑓𝜂;𝑝0;𝑝=𝑓_𝑤,𝑓^′0;𝑝=1,𝑓^′()ℒ∞;𝑝=0,1−𝑃₂𝜃𝜂;𝑝−𝜃₀𝜂=𝑝ℏ₂𝐻𝜂𝑁₂𝜃,𝜃𝜃𝜂;𝑝0;𝑝=1,𝜃∞;𝑝=0Power-lawtemperature,_′𝑁0;𝑝=−1,𝜃∞;𝑝=0Power-lawheatflux,₁𝑓=𝜕𝜂;𝑝,𝜃𝜂;𝑝₃𝑓𝜂;𝑝𝜕𝜂³𝑓𝜕+Re𝜂;𝑝²𝑓𝜂;𝑝𝜕𝜂²−𝜕𝑓𝜂;𝑝𝜕𝜂²−𝜕𝑓𝜂;𝑝,𝑁𝜕𝜂₂𝑓=𝜕𝜂;𝑝,𝜃𝜂;𝑝₂𝜃𝜂;𝑝𝜕𝜂²+Re⋅Pr𝑓𝜂;𝑝𝜕𝜃𝜂;𝑝𝜃𝜕𝜂−𝑛⋅Re⋅Pr𝜂;𝑝𝜕𝑓𝜂;𝑝.𝜕𝜂 (13) For 𝑝=0 and 𝑝=1, we have 𝑓𝜂;0=𝑓₀𝜂𝜂,𝜃,𝑓𝜂;1=𝑓𝜂;0=𝜃₀𝜂,𝜃𝜂;1=𝜃𝜂. (14) When 𝑝 increases from 0 to 1, then 𝑓(𝜂;𝑝) and 𝜃(𝜂;𝑝) vary from 𝑓₀(𝜂) and 𝜃₀(𝜂) to 𝑓(𝜂) and 𝜃(𝜂). By Taylor’s theorem and using (14), we can write the following: 𝑓𝜂;𝑝=𝑓₀𝜂+^∞_𝑚−1𝑓_𝑚𝜂𝑝^𝑚,𝑓_𝑚𝜂=1𝜕𝑚!^𝑚𝑓𝜂;𝑝𝜕𝑝^𝑚,𝜃𝜂;𝑝=𝜃₀𝜂+^∞_𝑚−1𝜃_𝑚𝜂𝑝^𝑚,𝜃_𝑚𝜂=1𝜕𝑚!^𝑚𝜃𝜂;𝑝𝜕𝑝^𝑚. (15) In which ℏ₁ and ℏ₂ are chosen in such a way that these two series are convergent at 𝑝=1; therefore we have the following through (15): 𝑓𝜂=𝑓₀𝜂+^∞_𝑚=1𝑓_𝑚𝜂,𝜃𝜂=𝜃₀𝜂+^∞_𝑚=1𝜃_𝑚𝜂. (16)

3.2. Consider Mth-Order Deformation Equations

ℒ₁𝑓_𝑚𝜂−𝜒_𝑚𝑓_𝑚−1𝜂=ℏ₁𝐻𝜂𝑅^𝑓_𝑚𝜂,𝑓_𝑚(0)=0,𝑓^′_𝑚(0)=0,𝑓^′_𝑚(∞)𝑅=0,_𝑓𝑚𝜂=𝑓^′′′_𝑚−1+^𝑚−1_𝑘=0𝑓_𝑘𝑓^′′_{𝑚−1−𝑘}−^𝑚−1_𝑘=0𝑓^′_𝑘𝑓^′_{𝑚−1−𝑘}−𝐾𝑓^′_𝑚−1,ℒ₂𝜃_𝑚𝜂−𝜒_𝑚𝜃_𝑚−1𝜂=ℏ₂𝐻𝜂𝑅^𝜃_𝑚𝜂,𝜃_𝑚(0)=0,𝜃_𝑚(∞)𝜃=0,_′𝑚(0)=0,𝜃_𝑚(∞)𝑅=0,_𝜃𝑚𝜂=𝜃^′′_𝑚−1+Re⋅Pr^𝑚−1_𝑘=0𝑓_𝑘𝜃^′_{𝑚−1−𝑘}−𝑛⋅Re⋅Pr^𝑚−1_𝑘=0𝜃_𝑘𝑓^′_{𝑚−1−𝑘},𝜒_𝑚=0,𝑚≤1,1,𝑚>1. (17)(18)(19)(20)(21)(22)(23)(24) The general solutions of (17)–(23) are 𝑓_𝑚𝜂−𝜒_𝑚𝑓_𝑚−1𝜂=𝑓^∗_𝑚𝜂+𝐶^𝑚₁+𝐶^𝑚₂𝜂+𝐶^𝑚₃,𝜃exp−𝜂_𝑚𝜂−𝜒_𝑚𝜃_𝑚−1𝜂=𝜃^∗_𝑚𝜂+𝐶^𝑚₄+𝐶^𝑚₅exp−𝜂, (25) where 𝐶^𝑚₁ to 𝐶^𝑚₅ are constants that can be obtained by applying the boundary condition in (18), (21), and (22).

As discussed by Liao, the rule of coefficient ergodicity and the rule of solution existence play important roles in determining the auxiliary function and ensuring that the high-order deformation equations are closed and have solutions. In many cases, by means of the rule of solution expression and the rule of coefficient ergodicity, auxiliary functions can be uniquely determined. So we define the auxiliary function 𝐻(𝜂) in the following form: 𝐻𝜂=exp−𝜂. (26)

4. Convergence of the HAM Solution

As was mentioned in introduction, the convergence and the rate of approximation of the HAM solution strongly depend on the values of auxiliary parameter ℏ. By means of the so-called ℏ curves, it is easy to find out the so-called valid regions of ℏ to gain a convergent solution series. According to Figures 2 and 3, the acceptable range of auxiliary parameter for Pr=1,𝑛=0,Re=1,and𝑓_𝑤=0 is −0.3<ℏ₁<−0.1 and −1.8<ℏ₂<−0.3.

Figure 2 

The ℏ1 validity for Pr=1, 𝑛=0, Re=1, and 𝑓𝑤=0.

Figure 3 

The ℏ2 validity for Pr=1, 𝑛=0, Re=1, and 𝑓𝑤=0.

Figures 4 and 5 show how auxiliary parameters ℏ₁ and ℏ₂ vary with changing 𝑓_𝑤. If 𝑓_𝑤 increases, the range of convergency of solution is restricted and then decreased.

Figure 4 

The ℏ1 validity for various 𝑓𝑤, Pr=1, 𝑛=0, and Re=1.

Figure 5 

The ℏ2 validity for various 𝑓𝑤, Pr=1, 𝑛=0, Re=1, and 𝑓𝑤=0.

5. Results and Discussion

In the present study, the Homotopy analysis method has been used to solve the nonlinear equations of heat transfer and fluid flow over a permeable stretching wall in a porous medium. The nondimensional numbers introduced in the present analysis are Reynolds number (Re), Prandtl number (Pr), and the injection number (𝑓_𝑤). Another important parameter is the power of the surface temperature/heat flux distribution (𝑛), which is considered in the present study.

In order to ensure the convergence of the solution series, the variation of 𝑓 has been plotted at different orders of approximation in Figure 6. The comparison between results of HAM and numerical solution (Runge-Kutta) has been displayed in Table 1. It can be observed that there is a good agreement between HAM method and numerical solution. Figures 7, 8, and 9 illustrate effect of variation of wall injection/suction parameter (𝑓_𝑤) on velocity and temperature distribution. It can be observed that all boundary layer thickness decreased by increasing 𝑓_𝑤 from negative to positive (from injection to suction). Figures 10, 11, and 12 show the effect of Reynolds number on the velocity profile when Pr=1, 𝑛=0, and 𝑓_𝑤=0. This figure shows that the boundary layer thickness and thermal boundary layer thickness are quite opposite to that of Reynolds number. The influence of 𝑛 on temperature field for both types of the thermal boundary conditions considered has been studied in Figures 13 and 14. We notice that increasing 𝑛 reduces the thermal boundary layer thickness regardless of the boundary condition type leading to a heat transfer augmentation. Effect of Prandtl number (isothermal) on the temperature field is plotted in Figure 15. This figure shows that thermal boundary layer thickness directly depends on Prandtl number.

Table 1 

The results of HAM and NS for 𝑓(𝜂), 𝑓^′(𝜂), and 𝜃(𝜂) when 𝑓_𝑤=0, Pr=1, 𝑛=0, and Re=1.

Figure 6 

The variation of 𝑓 at different orders of approximations, Pr=1, 𝑛=0, Re=1, and 𝑓𝑤=0.

Figure 7 

Velocity profile 𝑓′ for various 𝑓𝑤 when Pr=1, 𝑛=0, and Re=1.

Figure 8 

Temperature profile (isothermal) 𝜃 for various 𝑓𝑤 when Pr=1, 𝑛=0, and Re=1.

Figure 9 

Temperature profile (isoflux) 𝜃 for various 𝑓𝑤 when Pr=1, 𝑛=0, and Re=1.

Figure 10 

Velocity profile 𝑓′ for various Re when Pr=1, 𝑛=0, and 𝑓𝑤=0.

Figure 11 

Temperature profile (isothermal) 𝜃 for various Re when Pr=1, 𝑛=0, and 𝑓𝑤=0.

Figure 12 

Temperature profile (isoflux) 𝜃 for various Re when Pr=1, 𝑛=0, and 𝑓𝑤=0.

Figure 13 

Temperature profile (isothermal) 𝜃 for various 𝑛 when Pr=1, 𝑛=0, and 𝑓𝑤=0.

Figure 14 

Temperature profile (isoflux) 𝜃 for various 𝑛 when Pr=1, 𝑛=0, and 𝑓𝑤=0.

Figure 15 

Temperature profile 𝜃 (isothermal) for various Pr when Pr=1, 𝑛=0, and 𝑓𝑤=0.

Tables 2 and 3 compare the results of HAM and numerical solution when the Reynolds number varies for 𝑓_𝑤=0, Pr=1, and 𝑛=0 for isothermal and heat flux boundary conditions, respectively. It can be seen that for isothermal boundary condition with increasing Reynolds number, the wall shear stress for both boundary conditions considered here (isothermal and heat flux) and Nusselt number increase with Reynolds number for isothermal state and independent of Reynolds number for heat flux boundary condition.

Table 2 

The results of HAM and NS for 𝑓^′′(0), 𝜃^′(0) when 𝑓_𝑤=0, Pr=1, and 𝑛=0 for temperature.

Table 3 

The results of HAM and NS for 𝑓^′′(0), 𝜃^′(0) when 𝑓_𝑤=0, Pr=1, and 𝑛=0 for heat flux.

Comparison of the results of HAM and numerical solution has been shown in Tables 4 and 5 for various Prandtl numbers when Re=1, 𝑛=0 and 𝑓_𝑤=−0.3,𝑓_𝑤=0.3, respectively. It can be observed that the Prandtl number has no effect on the wall shear stress for both boundary conditions and Nusselt number for heat flux boundary condition.

Table 4 

The results of HAM and NS for 𝑓^′′(0), 𝜃^′(0) for various Pr when, Re=1, 𝑛=0, and 𝑓_𝑤=−0.3.

Table 5 

The results of HAM and NS for 𝑓^′′(0), 𝜃^′(0) for various Pr when, Re=1, 𝑛=0, and 𝑓_𝑤=0.3.

6. Summary and Conclusions

Homotopy analysis method (HAM) is applied to compute wall driven flow through a porous medium over a stretching permeable surface subjected to power-law temperature and heat flux boundary conditions. The proper range of the auxiliary parameter ℏ to ensure the convergency of the solution series was obtained through the so-called ℏ curves. When compared with other analytical methods, it is clear that HAM provides highly accurate analytical solutions for nonlinear problems. Moreover, second-law (of thermodynamics) aspects of the problem are investigated. The highlights of this study are the following. (i)The nondimensional viscous boundary layer thickness has a direct relationship with Reynolds number; thus Nusselt number rate increases with Re. (ii)Nusselt number, wall shear stress have a reverse relationship with and mass transfer from the wall, 𝑓_𝑤.(iii)Increasing the Prandtl number results in reduction of thermal boundary layer thickness.

Consequently, Nusselt number increase with Pr.

Nomenclature

𝑓:	Similarity function for velocity
𝑓𝑤:	Injection parameter, 𝑓𝑤=−𝜐𝑤𝐿/𝑢0√𝐾
ℏ:	Auxiliary parameter
HAM:	Homotopy analysis method
𝐻:	Auxiliary function
ℒ:	Linear operator of HAM
𝐾:	Permeability of the porous medium, m2
𝑁:	Non-linear operator
𝑛:	Power of temperature/heat flux distribution
Nu:	Local Nusselt number
Nu𝐿:	Averaged Nusselt number
Pr:	Prandtl number, Pr=𝑣/𝛼eff
𝑞0:	Wall heat flux coefficient, W/m2
Re:	Reynolds number, Re=𝜌𝑢0𝐾/𝐿
𝑇:	Temperature
𝑇0:	Wall temperature coefficient, 𝐾
𝑢:	Velocity in 𝑥 direction
𝑢0:	Wall velocity coefficient, m/s
𝜐:	Velocity in 𝑦 direction
𝜐𝑤:	Injection velocity, m/s
𝐵:	Positive constant
𝜃:	Similarity function for temperature
𝑥:	Coordinate system, m
𝑦:	Coordinate system, m
𝑄:	Volumetric rate of heat generation
𝜆:	Heat generation
𝑁:	Radiation parameter
𝜌:	Density of the fluid
Ψ:	Stream function
𝜇:	Dynamic viscosity
𝜎1:	Effective viscosity
𝜅1:	Absorption coefficient
𝜂:	Dimensionless similarity variable
𝑣:	Kinematic viscosity.