# A New Approach for the Solution of Three-Dimensional Magnetohydrodynamic Rotating Flow over a Shrinking Sheet

Authors: S. S. Motsa; S. Shateyi
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

The numerical solution of magnetohydrodynamic (MHD) and rotating flow over a porous shrinking sheet is obtained by the new approach known as spectral homotopy analysis method (SHAM). Using a similarity transformation, the governing equations for the momentum are reduced to a set of ordinary differential equations and are solved by the SHAM approach to determine velocity distributions and shear stress variations for different governing parameters. The SHAM results are analysed and validated against numerical results obtained using MATLAB's built-inbvp4c routine, and good agreement is observed.


---

## Body

## 1. Introduction

The study of flow, heat, and mass transfer problems due to stretching boundary/surface has many applications in technological processes, particularly in polymer systems involving drawing of fibres and films or thin sheets, production of paper, linoleum, roofing shingles, insulting material, and many other applications. In most cases, the polymer sheet is stretched while it is extruded from the dye. The sheet is pulled through viscous liquid with a cooling system to obtain the final product with prescribed characteristics. The moving sheet may introduce a motion in the neighbouring fluid, or alternatively, the fluid may have an independent forced-convection motion which is parallel to that of the sheet. Sakiadis [1] was the first to investigate the flow due to a sheet issuing with constant speed from a slit into a fluid at rest. Since then, many investigators have considered various aspects of this problem and have obtained similarity solutions, and a good amount of references can be found in papers by Crane [2], Magyari and Keller [3–5], Liao and Pop [6], Sparrow, and Abraham [7] and Abraham and Sparrow [8], among others.On the other hand, the laminar incompressible boundary layer flow caused by the stretching of a flat surface in rotating fluid has been studied by Wang [9], Rajeswari and Nathi [10] and Nazar et al. [11]. Ariel [12] presented a noniterative numerical scheme which computes the steady, three-dimensional flow of a viscous incompressible fluid past a stretching sheet in single integration. Rashidi and Dinarvand [13] found a totally analytic solution for the problem of condensation or spraying on an inclined rotating disk.In recent years, problems involving magnetic field have become important. Many metallurgical processes such as drawing, annealing and thinning of copper wire involve the cooling of continuous strips or filaments by drawing them through an ambient fluid. By drawing these filaments in an electrically conducting fluid under the influence of an applied magnetic field, controls the rate of cooling. Kumari and Nath [14] studied, using the homotopy analysis method, the unsteady magnetohydrodynamic viscous fluid and heat transfer of Newtonian fluids induced by an impulsively stretched plane surface in two lateral directions. Fang et al. [15] analytically solved the MHD flow under slip condition over a permeable stretching surface.The boundary layer flow and heat transfer problem over a moving surface differs from that over a stationary surface caused by the free stream velocity due to the entrainment of the fluid. The moving surface prevents or delays the separation of the boundary layer from the wall by injecting momentum in the existing boundary layer. The magnetic field and the rotation of the fluid increase the surface shear stress for primary flow, but reduce the surface heat transfer. Takhar et al. [16] studied the nonsimilar boundary layer flow of a viscous incompressible electrically conducting fluid over a moving surface in a rotating fluid, in the presence of a magnetic field, Hall currents and the free stream velocity. Vajravelu and Kumar [17] analyzed hydromagnetic flow between two horizontal plates in a rotating system, where the lower is a stretching sheet and the upper is a porous solid sphere.Literature survey indicates that very little attention has been given to the shrinking flow. Wang [18] developed unsteady shrinking sheet for a specific value of the suction parameter. The rotating flow of an electrically fluid occurs in cosmical and geophysical fluid dynamics. It is also important in the solar cycle and the structure of rotating magnetic stars. Hayat et al. [19] obtained series solution of magnetohydrodynamic and the rotating flow over a porous shrinking sheet using a homotopy analysis method. Sajid and Hayat [20] considered the MHD viscous flow due to a shrinking sheet. The study obtained series solution valid for both two dimensional and axisymmetric shrinking sheet by using homotopy analysis method. Yao and Chen [21] applied the homotopy analysis method to investigate analytically the laminar incompressible viscous flow for a moving semi-infinite flat, or a flat plate continuously shrinking into a slot in a stationary fluid with mass transfer governed by the Blasius equation.Noor et al. [22] examined analytically the magnetohydrodynamic (MHD) viscous flow due to a shrinking sheet using the Adomian decomposition method (ADM) coupled with Padé approximants to handle the condition at infinity. Muhaimina et al. [23] studied the effect of the thermophoresis particle deposition on nonlinear MHD mixed convective heat and mass transfer over a porous shrinking sheet in the presence of suction.The main objective of the present study is to find the solution for the problem of three-dimensional rotating flow induced by shrinking sheet with suction using the recently developed SHAM approach [24]. The problem was previously considered in [19] using the standard homotopy analysis method and in this work we use the new SHAM approach to solve the same problem. The SHAM method was shown to produce more improved results than the traditional HAM, see Motsa et al. [24]. In this work we compare the SHAM results for velocity distributions and shear stresses at the bounding walls against results generated using the very efficient MATLAB bvp4c in-built routine. The comparison indicates that there is excellent agreement between the two results proving that the SHAM is at least as good as the bvp4c and can be used in place of traditional numerical approaches such as Runge-Kutta methods, finite differences, Keller-Box method, for solving nonlinear boundary value problems.


## 2. Mathematical Formulation

We consider the steady, incompressible, three-dimensional flow of an electrically conducting viscous fluid between two horizontal parallel plates aty=h. Both the fluid and the plates rotate in unison with a constant angular velocity Ω=Ωj, where j is a unit vector in the y-direction. The plate y=+h is rigid and stationary. The flow in the fluid system is caused due to shrinking of a porous plate at y=-h. The equations governing the rotating flow are Vajravelu and Kumar [17]:(2.1)∂u∂x+∂v∂y=0,(2.2)u∂u∂x+v∂u∂y+2Ωw=-1ρ∂p*∂x+ν[∂2u∂x2+∂2u∂y2]-σBo2ρu,(2.3)u∂u∂y=-1ρ∂p*∂y+ν[∂2v∂x2+∂2v∂y2],(2.4)u∂w∂x+v∂w∂y-2Ωw=ν[∂2w∂x2+∂2w∂y2]-σBo2ρw.
The boundary conditions for the problem considered here are:(2.5)u=-ax,v=-V,w=0,aty=-h,u=0,v=0,w=0,aty=+h,
where u,v, and w are the velocity components in x-, y-, and z-directions, respectively, ρ is the density, ν is the kinematic viscosity, σ is the electrical conductivity, Bo is the magnetic induction, p* is the modified pressure, a>0 is the shrinking constant and V>0 is the suction velocity. In order to reduce (2.1)–(2.4) into a set of convenient ordinary differential equations, we introduce the similarity variable η and the dimensionless variable f and g as follows:(2.6)η=yh,u=-axf′(η),v=ahf(η),w=axg(η).
The continuity equation (2.1) is automatically satisfied and (2.2)–(2.4), after eliminating the modified pressure, are characterized by the following model equations:(2.7)fiv-M2f′′-2Kg′-R(f′f′′-ff′′′)=0,(2.8)g′′-M2g+2K2f′-R(f′g-fg′)=0,
where primes indicates differentiation with respect to η. In view of equation (2.6), the boundary equations (2.5), transform into:(2.9)f=λ,f′=-1,g=0,atη=-1,(2.10)f=0,f′=0,g=0,atη=1,
in which the suction parameter λ, the viscosity parameter R, the Hartman number M, and the rotating parameter K2 are:(2.11)λ=-Vah,R=ah2ν,M2=σBo2h2ρν,K2=Ωh2ν.In the next section we will solve the nonlinear ordinary equation (2.7) by using the spectral homotopy analysis method (SHAM).


## 3. Spectral Homotopy Analysis Method Solution

In this section, we apply the SHAM approach to solve the governing equations (2.7)–(2.10). We begin by introducing the following transformation(3.1)f(η)=f0(η)+F(η),g(η)=g0(η)+G(η),
where(3.2)f0(η)=(λ2-14)+(1-3λ4)η+η24+(λ-14)η3,g0(η)=1-η2
are the initial approximations which are chosen to satisfy the boundary conditions (2.9)–(2.10). Equation (3.1) is substituted into the governing equations (2.7)–(2.10) with the resulting equations written as a sum of their linear and nonlinear components as(3.3)L1[F(η),G(η)]+N1[F(η),G(η)]=Φ(η),(3.4)L2[F(η),G(η)]+N2[F(η),G(η)]=Ψ(η),
subject to the boundary conditions(3.5)F=0,F′=0,G=0,atη=-1,(3.6)F=0,F′=0,G=0,atη=1,
where(3.7)L1[F(η),G(η)]=Fiv+a1F′′′+a2F′′+a3F′+a4F-2K2G′,(3.8)L2[F(η),G(η)]=G′′+a1G′+a2G+b1F′+b2F,(3.9)N1[N(η),G(η)]=-R(F′F′′-FF′′′),(3.10)N2[N(η),G(η)]=-R(F′G-FG′),(3.11)Φ(η)=-[f0iv-M2f0′′-2K2g0′-R(f0′f0′′-f0f0′′′)],(3.12)Ψ(η)=-[g0′′-M2g0+2K2f0′-R(f0′g0-f0g0′)].
In the above definitions, the coefficient parameters are defined as(3.13)a1=R0f0(η),a2=-M2-Rf0′(η),a3=-Rf0′′,a4=Rf0′′′,b1=2K2-Rg0,b2=Rg0′.
The SHAM approach builds on the basic ideas of the homotopy analysis method (HAM). However, for brevity, details of the HAM are omitted in this paper. For a detailed exposition of the HAM approach interested readers can refer to [25, 26] for a general description on the method and to [27–46] for the application of the HAM in boundary value problems over bounded domains. Thus, importing the ideas of the HAM approach, we construct the so-called zero-order deformation equations as(3.14)(1-q)L1{[F̃(η),G̃(η)]-[F0(η),G0(η)]}=qℏ{L1[F̃(η),G̃(η)]+N1[F̃(η),G̃(η)]-Φ(η)},(1-q)L2{[F̃(η),G̃(η)]-[F0(η),G0(η)]}=qℏ{L2[F̃(η),G̃(η)]+N2[F̃(η),G̃(η)]-Ψ(η)},
where ℏ is the convergence controlling parameter, q∈[0,1] is the embedding parameter, F̃(η) and G̃(η) are unknown functions and F0(η) and G0(η) are initial approximations which are obtained as solutions of the linear part of equations (3.3)–(3.6) given as(3.15)L1[F0(η),G0(η)]=Φ(η),(3.16)L2[F0(η),G0(η)]=Ψ(η),
subject to the boundary conditions(3.17)F0=0,F0′=0,G0=0,atη=-1,(3.18)F0=0,F0′=0,G0=0,atη=1.
Following the HAM approach, the zero-order deformation deformation equations are differentiated m times with respect to the embedding parameter q then divided by m! with q=0 being set to the resulting equations to obtain the so-called higher-order deformation equations given by(3.19)L1{[Fm(η),Gm(η)]-χm[Fm-1(η),Gm-1(η)]}=ℏRm,1,L2{[Fm(η),Gm(η)]-χm[Fm-1(η),Gm-1(η)]}=ℏRm,2,
subject to the boundary conditions(3.20)Fm(-1)=Fm′(-1)=Fm(1)=F′(1)=Gm(-1)=Gm(1)=0,
where(3.21)Rm,1=Fm-1iv+a1Fm-1′′′+a2Fm-1′′+a3Fm-1′+a4Fm-2K2Gm-1′-R∑n=0m-1(Fn′Fm-n-1′′-FnFm-n-1′′′)-(1-χm)Φ(η),Rm,2=Gm-1′′+a1Gm-1′+a2Gm-1+b1Fm-1′+b2Fm-R∑n=0m-1(Fn′Gm-n-1-FnGm-n-1′)-(1-χm)Ψ(η),χm={0,m≤1,1,m>1.
We remark that, unlike in the standard HAM approach, the higher-order deformation equations (3.19)–(3.20) form a set of coupled ordinary differential equations (ODEs) instead of the decoupled set of ODEs that are generated in the HAM, that is in the SHAM approach the linear operators depend on both F and G (see (3.7) and (3.8)) whereas in the case of the HAM the linear operator would depend on one variable at a time. The SHAM technique also doe not depend on the rule of solution expression and the rule of ergodicity unlike the standard HAM. We use the Chebyshev pseudospectral method (see, e.g., [47, 48]) to solve equations (3.19)–(3.20). The unknown functions Fm(η) and Gm(η) are approximated as truncated series of Chebyshev polynomials of the forms(3.22)Fm(η)≈∑k=0NF̂m(ηk)Tk(ηj),Gm(η)≈∑k=0NĜm(ηk)Tk(ηj),
where j=0,1,…,N, Tk is the kth Chebyshev polynomial, and η0,η1,…,ηN are Gauss-Lobatto collocation points (see [47]) defined by(3.23)ηj=cos⁡πjN,j=0,1,…,N.
Derivatives of the functions Fm(η) and Gm(η) at the collocation points are represented as(3.24)drFmdηr=∑k=0NDkjrFm(ηj),drGmdηr=∑k=0NDkjrGm(ηj)
where r is the order of differentiation and 𝒟 is the Chebyshev spectral differentiation matrix ([47, 48]).Substituting equations (3.22)–(3.24) in (3.19)–(3.20) yields(3.25)[A11A12A21A22][FmGm]=(χm+ℏ)[A11A12A21A22][Fm-1Gm-1](3.26)-ℏ(1-χm)[ΦΨ]+[Pm-1Qm-1],
subject to the boundary conditions (3.27)Fm(ηN)=0,Gm(ηN)=0,Fm(η0)=0,Gm(η0)=0,(3.28)∑k=0NDNkFm(ηk)=0,∑k=0ND0kFm(ηk)=0,
where(3.29)A11=D4+a1D3+a2D2+a3D+a4,A12=-2K2D,a21=b1D+b2,A22=D2+a1D+a2,Fm=[Fm(η0),Fm(η1),…,Fm(ηN)]T,Gm=[Gm(η0),Gm(η1),…,Gm(ηN)]T,Φ=[Φ(η0),Φ(η1),…,Φ(ηN)]T,Ψ=[Ψ(η0),Ψ(η1),…,Ψ(ηN)]T,Pm-1=-R∑n=0m-1(DFnD2Fm-1-n-FnD3Fm-1-n),Qm-1=-R∑n=0m-1(DFnGm-1-n-FnDGm-1-n).
In the above definitions the superscript T denotes transpose, and ai,bi denotes diagonal matrices, I is an identity matrix of size (N+1)×(N+1). The boundary conditions (3.25) and (3.27) are imposed on equation (3.24) are the resulting equation is solved for Fm and Gm iteratively using F0(η) and G0(η), which are obtained as solutions of (3.15)–(3.18), as a starting point.


## 4. Results and Discussion

In this section we give the SHAM results for the four main parameters affecting the flow. We remark that, all the SHAM results presented in this work were obtained usingN=50 collocation points. Tables 1–4 give a comparison of the SHAM results for f′′(-1), f′′(1), g′(-1) and g′(1) at different orders of approximation against the numerical results. The numerical results are obtained using the MATLAB routine bvp4c. Table 1 shows that full convergence of the SHAM is achieved by as early as the second-order, substantiating the claim that SHAM is a very powerful technique. We observe that convergence is achieved at second-order of approximation for all parameter values or combinations of these parameters as depicted in all these tables. We observe in Table 1 that the suction parameter λ significantly affects the shear stress exerted by the shrinking sheet at η=-1. Increasing the values of λ causes much reductions in the shear stress at η=-1 as shown by both values of f′′(-1) and g′(-1). This is because blowing gives rise to a thicker velocity boundary layer, thereby causing a decrease in the velocity gradient at the surface.Table 1
Comparison of the values of wall shear stressesf′′(-1), g′(-1) with the numerical solution for different orders of the SHAM approximation when λ is varied with ℏ=-1, M=0.5, K=0.5, and R=0.2.
λf′′(-1)g′(-1)1st-order2nd orderNumerical1st-order2nd orderNumerical0.02.0367552.0367552.036755-0.086956-0.086956-0.0869560.51.2907681.2907681.290768-0.204070-0.204070-0.2040701.00.5173180.5173180.517318-0.337831-0.337831-0.3378311.5-0.285247-0.285247-0.285247-0.490317-0.490317-0.4903172.0-1.118646-1.118646-1.118646-0.663878-0.663879-0.6638792.5-1.984669-1.984669-1.984669-0.861178-0.861182-0.8611823.0-2.885172-2.885172-2.885172-1.085248-1.085259-1.0852593.5-3.822080-3.822080-3.822080-1.339543-1.339572-1.3395724.0-4.797383-4.797383-4.797383-1.628018-1.628084-1.628084Table 2
Comparison of the values of wall shear stressesf′′(-1), g′(-1) with the numerical solution for different orders of the SHAM approximation when M is varied with ℏ=-1, λ=0.5, K=0.5, and R=0.5.
Mf′′(-1)g′(-1)1st-order2nd orderNumerical1st-order2nd orderNumerical0.01.2145391.2145391.214539-0.226272-0.226273-0.2262730.51.2711251.2711251.271125-0.213930-0.213930-0.2139301.01.4312381.4312371.431237-0.186192-0.186200-0.1862001.51.6716891.6716871.671687-0.156767-0.156789-0.1567892.01.9666841.9666781.966678-0.131647-0.131681-0.1316812.52.2952722.2952622.295262-0.111685-0.111723-0.1117233.02.6433252.6433102.643310-0.096121-0.096157-0.0961573.53.0023403.0023223.002322-0.083955-0.083986-0.0839864.03.3675163.3674963.367496-0.074329-0.074355-0.074355Table 3
Comparison of the values of wall shear stressesf′′(-1), g′(-1) with the numerical solution for different orders of the SHAM approximation when K is varied with ℏ=-1, λ=0.5, R=0.2, and M=0.5.
Kf′′(-1)g′(-1)2nd order4th orderNumerical2nd order4th orderNumerical0.01.2907781.2907781.290778-0.000000-0.000000-0.0000000.51.2907681.2907681.290768-0.204070-0.204070-0.2040701.01.2896291.2896291.289629-0.808721-0.808721-0.8087211.51.2678761.2678761.267876-1.758288-1.758288-1.7582882.01.1381351.1381371.138137-2.939868-2.939865-2.9398652.50.7838600.7838690.783869-4.302754-4.302744-4.3027443.00.1767550.1767710.176771-5.885655-5.885639-5.8856393.5-0.664164-0.664146-0.664146-7.721068-7.721051-7.7210514.0-1.740066-1.740049-1.740049-9.812337-9.812322-9.812322Table 4
Comparison of the values of wall shear stressesf′′(-1), g′(-1) with the numerical solution for different orders of the SHAM approximation when R is varied with ℏ=-0.98, λ=0.5, K=0.2, and M=0.5.
Rf′′(-1)g′(-1)2nd order4th orderNumerical2nd order4th orderNumerical0.01.3033321.3033321.303332-0.197931-0.197931-0.1979310.21.2907681.2907681.290768-0.204070-0.204070-0.2040700.41.2777821.2777821.277782-0.210552-0.210552-0.2105520.61.2643541.2643541.264354-0.217406-0.217406-0.2174060.81.2504631.2504631.250463-0.224665-0.224665-0.2246651.01.2360851.2360851.236085-0.232365-0.232365-0.2323652.01.1559581.1559581.155958-0.279118-0.279118-0.2791183.01.0590781.0590781.059078-0.346195-0.346202-0.346202From Table2, it is observed that the Hartman number tends to greatly increase the local skin friction at the shrinking sheet (η=-1). This is because the increase in the magnetic field strength leads to a thinner velocity boundary layer, thereby causing an increase in the velocity gradient at the wall. In Table 3 we observe the influence of the rotation parameter K on the shear stress f′′(-1) and g′(-1). We observe that both f′′(-1) and g′(-1) decrease as the values of K increase. In Table 4 we observe that f′′(-1) decreases by increasing R and g′(-1) increases as R increases.Figures1–4 have been plotted to depict the influence of suction parameter λ, the Hartman number M, rotation parameter K and viscosity parameter R. On these figures, we also give comparisons between the numerical results and the second-order SHAM solutions and excellent agreement between the two sets of results was always achieved. In Figure 1, we have the effects of varying the values of suction parameter λ on f, f′ and g. From Figure 1 it is found that f increases as λ increases and f has maximum values at the lower end of the plate (shrinking sheet). It is clearly depicted in Figure 1 that f′ decreases when values of the suction parameter increase. It is also observed that for small values of λ, f′  has large values near the center of the channel. Figure 1 also elucidates the effects of λ on g. We observe that g decreases as λ increases and the decrease is more pronounced at the center of the channel as compared to near the plates.Influence of the suction parameterλ on f(η), f′(η) and g(η), when ℏ=-1. Numerical solution (solid line) is compared against the SHAM 1st-order approximation (diamonds).
(a)(b)(c)Influence ofK on f(η), f′(η) and g(η), when ℏ=-1. Numerical solution (lines) is compared against the SHAM 1st-order approximation (open circles).
(a)(b)(c)Influence ofM on f(η), f′(η) and g(η), when ℏ=-1. Numerical solution (lines) is compared against the SHAM 1st-order approximation (open circles).
(a)(b)(c)Influence ofR on f(η),  f′(η) and g(η), when ℏ=-1. Numerical solution (lines) is compared against the SHAM 1st-order approximation (open circles).
(a)(b)(c)In Figure2 we depict the effects of the rotation parameter K on f, f′ and g. We observe in this figure that near the shrinking plate, the rotation parameter has no effect on f. However, as we move towards the center of the channel f increases as K increases. Figure 2 indicates that f′ increases near the shrinking sheet and also that the boundary layer thickness decreases near this sheet. As we approach the nonpermeable plate, we observe that f′ is now a decreasing function of K. In this figure we observe that f′ is not a monotonous function of K. We also have the effects of K on g depicted in Figure 2. We clearly see that g decreases as K increases.Figure3 depicts the effects of M on f, f′, and g. We observe in this figure that f is an increasing function of the Hartman number M. It is observed in Figure 3 that f′ initially increases but then decreases after the center of the channel as values of M increase. The Hartman number M significantly reduces the values of g. We observe in this Figure 3 that g has quite opposite behaviour when compared with the suction parameter λ.Lastly, in Figure4 we show the effects of viscosity parameter R on f, f′ and g. We observe that viscosity reduces the velocity f and increases the boundary layer thickness. The minimum values of f are observed near the center of the channel. It is noted in Figure 4 that increasing the values of R initially decreases f′ but increases it after the channel center. It can also be observed in Figure 4 that R significantly affects g. As R increases, g values are greatly reduced attaining their minimum values near the shrinking sheet.


## 5. Conclusion

The three-dimensional rotating flow in a channel generated by a shrinking sheet is studied. The spectral-homotopy analysis method is used to solve the nonlinear system of ordinary differential equations. The variations of the four main parameters on the velocity(f,f′,g) and wall shear stress (f′′(-1),-g′(-1)) are discussed through graphs and tables, respectively. The following observations have been made.  (i)
The SHAM rapidly converges to the numerical results generated by MATLABbvp4c routine.(ii)
The velocityf increases for λ, M, and K but decreases for R.(iii)
The velocityf′ decreases for increasing values of λ but is not a monotonous function of K, M, and R.(iv)
The velocityg decreases for increasing values of λ, K and R but increases for increasing values of M.



---

*Source: 586340-2010-12-08.xml*