On a linearisation method for Reiner-Rivlin swirling flow

On a linearisation method for Reiner-Rivlin swirling flow

Zodwa G. Makukula; Precious Sibanda^*; Sandile S. Motsa

School of Mathematical Sciences, University of KwaZulu-Natal, Private Bag X01, Scottsville 3209, Pietermaritzburg, South Africa E-mails: [e-mail] / [e-mail] / [e-mail]

ABSTRACT

The steady flow of a Reiner-Rivlin fluid with Joule heating and viscous dissipationis studied. We present a novel technique for accelerating the convergence of the spectral-homotopy analysis method. Solutions of the nonlinear momentum and energy equations are obtained using the improved spectral homotopy analysis method. Solutions were also generated using the spectral-homotopy analysis method and benchmarked against results in the literature.

Mathematical subject classification: Primary: 76A05, 76N05; Secondary: 76M25.

Key words: Reiner-Rivlin fluid, Chebyshev spectral method, spectral-homotopy analysis method, improved spectral-homotopy analysis method.

1 Introduction

The boundary layer induced by a rotating disk arises in many engineering applications, for example, in computer storage devices, viscometry, turbo-machinery and in crystal growth processes (Attia [4]). Since the pioneering study by von Kármán [39], research on swirling flows has been carried out by, among others, Cochran [15] who proposed an improved solution to the von Kármán formulation based on a mixture of analytical and numerical techniques. Benton [12] studied the impulsive rotation from rest of a disk in an infinite viscous fluid. He improved Cochran's solutions by first expanding the variables in a power series and solving for the first two orders analytically, and then numerically computing the next two orders.

The shooting method was used to solve the von Kármán equations and to investigate heat transfer in porous medium in [19, 29, 30, 33, 35]. Numericalschemes involving Runge-Kutta methods, finite element and finite difference approximations were used in [14, 18, 24] to study the effects of a rough disk surface on the flow. The Crank-Nicholson implicit scheme was used in studies involving non-Newtonian characteristic of the fluid [4, 5, 6, 8, 10, 33] and in Newtonian fluids [7, 9, 11]. Perturbation techniques including the differential transform method and Padé approximations (DTM-Padé), the variational iteration method (VIM) and the homotopy perturbation method (HPM) were used in [1, 2, 31, 32, 34]. Analytical methods such as the DTM and the homotopy analysis method (HAM) were applied in, among other studies, [40, 3, 17, 37, 38]. These methods may result in secular terms in the solutions, and may converge very slowly or may even fail to converge for problems with strong non-linearity and/or with very large parameter values. Some approaches may not be applicable at all for certain problems, for example, the DTM for unbounded domainproblems. Motsa et al. [22, 23, 25, 26, 27] proposed and applied a modification of the homotopy analysis method that improves the performance of this method and removes some restrictions associated with it. In general, the convergence of many numerical methods depends on how good the initial approximation is to the true solution. In this paper we present an algorithm that first seeks to improve the initial "guess" and then uses the spectral-homotopy analysis method to find solutions to systems of nonlinear equations that governthe Reiner-Rivlin swirling flow. This procedure considerably accelerates the convergence rate of the spectral homotopy analysis method. Here we apply the improved spectral-homotopy analysis method (ISHAM) to solve the nonlinear equations that govern the flow of an electrically conducting Reiner-Rivlin fluid in the presence of Joule heating and viscous dissipation. The governing equations were solved earlier by Sahoo [33] using a numerical scheme that blends the finite difference scheme and the shooting method. Solutions obtained are compared with those of Sahoo [33] and against the 'standard' spectral homotopyanalysis method.

2 Equations

We consider an infinite rotating disk coinciding with the plane z = 0 with the space z > 0 occupied by a viscous, incompressible Reiner-Rivlin fluid. Thefluid motion and heat transfer are governed by the equations (see [4, 5, 33]);

with the following no-slip boundary conditions

where the disk is rotating with a constant angular velocity Ω about the line r = 0 and an external uniform magnetic field is applied perpendicular to theplane of the disk with a constant magnetic flux density B₀. The velocity components in the directions of increasing r, ϕ, z are u, v, w respectively. ρ is the density of the fluid, σ is the electrical conductivity of the fluid, µ is the coefficient of viscosity, κ is the thermal conductivity, c_p is the specific heat at constant pressure of the fluid. The temperature of the fluid T, equals T_w at the surface of the disk. At large distances from the disk, T tends to T_∞ where T_∞ is the temperature of the ambient fluid. The second term on the right hand side of equation (5) represents the viscous dissipation while the last term represents the Joule heating. The constitutive equation for the Reiner-Rivlin fluidis given by

where p represents the pressure, is the stress tensor, is the rate of straintensor and µ_c is the coefficient of cross viscosity. The Reiner-Rivlin model is a simple model which can provide some insight into predicting the flow characteristics and heat transfer performance for viscoelastic fluid above a rotating disk [6]. The first term on the right hand side of (8) represents the viscous property of the fluid and the third term, the elastic property of the fluid. We introduce the non-dimensional distance η = z measured along the axis of rotation and the von Kármán transformations [39];

where F, G, H, P and Θ are non-dimensional functions of η, ν = µ/ρ is the kinematic viscosity. With these transformations equations (1)-(5) take the form

with

where K = µ_c Ω/µ is the parameter that describes the non-Newtonian characteristic of the fluid, M = σ/ρΩ is the magnetic interaction number, Pr is the Prandtl number and Ec is the Eckert number. The system (10)-(12) with the prescribed boundary conditions (15) are sufficient to solve for the three velocity components. Equation (13) can be used to find the pressure distribution at any point if required. Simplifying the equation system by substituting equation (10) into (11), (12) and (14) yields

subject to the boundary conditions

In the following section we solve the nonlinear coupled system (17)-(19) with boundary conditions (20) by the ISHAM.

3 Method of solution

The main thrust of the method of solution [21, 28], is the improvement of the initial approximation used in the higher order deformation equations of the spectral homotopy analysis method. A systematic approach is used to find optimal initial "guesses" which are then used in the SHAM algorithm to accelerate convergence. In the first instance we assume that solutions for H(η), G(η) and Θ(η) in equations (17)-(19) can be found in the form

where h_i, g_i and θ_i are unknown functions whose solutions are obtained usingthe SHAM approach at the ith iteration and h_m, g_m and θ_m (m > 1) are known from previous iterations. For m = 0, suitable initial guesses satisfying the boundary conditions (20) are

The initial guesses (22) are improved upon as follows. Substituting (21) intothe governing equations (17)-(19) gives

subject to the boundary conditions

The coefficient parameters a_{k, i - 1}, b_{k, i - 1}, c_{k, i - 1} (k = 0, ..., 6), r_{1, i - 1}, r_{2,i - 1} and r_{3,i - 1} are defined as

Starting from the initial guesses (22), the subsequent solutions h_i, g_i and θ_i (i > 1) are obtained by recursively solving equations (23)-(25). To solve equations (23)-(25), we start by defining the following linear operators

where q ∈ [0, 1] is the embedding parameter, and

are unknown functions. The zeroth order deformation equations are given by

where ħ is the non-zero convergence controlling auxiliary parameter and _h, _g and _θ are nonlinear operators given by

Differentiating (43)-(45) m times with respect to q and then setting q = 0 and finally dividing the resulting equations by m! yields the mth order deformation equations

subject to the boundary conditions

where

and

The initial approximations h_i,0, g_i,0 and θ_i,0 that are used in the higher order equations (46)-(48) are obtained by solving the linear part of equations (23)-(25) given by

with the boundary conditions

It is worthwhile to note at this stage that the initial approximate solutions are no longer just h₀, g₀ and θ₀ but h_i,0, g_i,0 and θ_i,0 at the ith iteration. Essentially, this procedure allows for the improvement of the initial guesses at each iteration. Since the right hand side of equations (54)-(56) for i = 1, 2, 3, ..., areknown from previous iterations, the equations may be solved using any numerical method. In this work, we apply the Chebyshev spectral collocation method to integrate equations (54)-(56). The method is based on the Chebyshev polynomials defined on the interval [-1, 1] by

We first transform the physical region [0, ∞) into the region [-1, 1] using the domain truncation technique. The problem is solved in the interval [0, L] instead of [0, ∞). This leads to the following algebraic mapping

where L is the scaling parameter used to invoke the boundary condition at infinity. The Chebyshev nodes in [-1, 1] are defined by the Gauss-Lobatto collocation points [13, 36] given by

where N is the number of collocation points. The variables h_i,0(η), g_i,0(η) and θ_i,0(η) are approximated as truncated series of Chebyshev polynomials of the form

where T_1,k, T_2,k and T_3,k are the k^th Chebyshev polynomials. Derivatives of the variables at the collocation points are represented as

where r is the order of differentiation, D = D and D is the Chebyshev spectral differentiation matrix. Substituting equations (61)-(64) in (53)-(56) yields

subject to the boundary conditions

where

In the above definitions T stands for transpose, I is an (N + 1) × (N + 1) identity matrix and a_{k,i - 1}, b_{k,i - 1} and c_{s,i - 1} (k = 0,...,5, s = 0,...,6) are diagonal matrices of size (N + 1) × (N + 1). After modifying the matrix system (65) to incorporate the boundary conditions, the solution is obtained as

Similarly, applying the Chebyshev spectral transformation on the higher order deformation equations (46)-(48) gives

subject to the boundary conditions

where B_{i - 1} and Q_{i - 1}, are as defined in (69) and

The boundary conditions (72)-(74) are implemented in matrix B_{i - 1} on the left hand side of equation (71) in rows 1, N, N + 1, N + 2, 2(N + 1) 2N + 3 and 3(N + 1) respectively as before with the initial solution above. The corresponding rows, all columns, of B_{i - 1} on the right hand side of (71), Q_{i - 1}and P_{m - 1} are all set to be zero. This results in the following recursive formula for m > 1.

where is the modified matrix B_{i - 1} on the right hand side of (71) afterincorporating the boundary conditions (72)-(74). Thus starting from the initial approximation, which is obtained from (70), higher order approximations X_i,m(ξ) for m > 1, can be obtained through the recursive formula (77). The solutions for h_i, g_i and θ_i are then generated using the solutions for h_i,m, g_i,m and θ_i,m as follows

The [i, m] approximate solutions for h(η), g(η) and θ(η) are then obtained by substituting h_i, g_i and θ_i which are obtained from (78), (79) and (80) into equation (21), where i represents the ith iteration of the initial approximation and m represents the mth iteration of the spectral homotopy analysis method.

4 Convergence theorem

The approximate solutions of the nonlinear equations are generated using the higher order deformation equations (46)-(48). The right hand sides of theseequations are governed by the unknown functions

_i(η; q), _i(η; q) and Θ(η; q). As the embedding parameter q gradually increases from 0 to 1, the solutions vary from the initial approximations to the exact solutions, i.e.

Expanding

_i(η; q), _i(η; q) and Θ_i(η; q) using the Taylor series expansion about q yields

We note that at q = 1 the series becomes the exact solutions

For validity of the solutions generated by these equations, it is important to show that these series converge at q = 1. As stated earlier, the SHAM is a hybrid method founded on the HAM. We kindly refer readers to Liao's proof [20, ch. 3] since the higher order deformation equations are similar in the two methods.

5 Results and discussion

In this section we present and discuss results computed using the improvedspectral homotopy analysis method, the original spectral homotopy analysis method and the numerical bvp4c routine which is based on Runge-Kuttaschemes. Comparison is also made between the current results and those in the literature. For our simulations we used ħ = -1, L = 30 and N = 150. The CPU run times (RT) in seconds are shown for the ISHAM and SHAM for comparison of computational efficiency.

Tables 1 and 2 present approximate solutions of the shear stresses in the radial F'(0) and tangential -G'(0) directions respectively. The results are computed for a Newtonian fluid (K = 0) and for different values of the magnetic parameter M. We note that the ISHAM approximate solutions for both F'(0) and -G'(0) converge to the numerical solutions at 2^nd order approximations for up to 8 decimal places. Comparison with Sahoo [33] shows a good agreement. The effect of the magnetic parameter on the Newtonian fluid shows that F'(0) decreases while -G'(0) increases as M is increased.

The radial shear stress when M = 2 and for different values of K is presented in Table 3. For validation of the current method, and to determine the effect of improving the initial guesses, we compare the results against the numerical solution. For convergence of the method the results are compared with the 'standard' spectral homotopy analysis method for the same values of N, L and ħ. For 0 < K < 2 the ISHAM converges at 2^nd order while the SHAM would not have converged even at the 8^th order for some values of K. Comparatively therefore the ISHAM converges much faster than the SHAM. This is clearly seen in Table 4 where the absolute errors in the solution are given.

In Table 5, the tangential stress results obtained using the ISHAM and the SHAM are compared with the numerical results when M = 2 and for different values of K. Convergence of the ISHAM to the numerical solutions was achieved at the 2^nd order of approximation. When using the SHAM convergence to the numerical solution was achieved at the 8^th order for values of K up to 2. It is also clear that the ISHAM is computationally much more efficient compared to the SHAM. A comparison of the absolute errors is made in Table 6 where the fast convergence of the ISHAM when compared with the SHAM is confirmed. It is worth noting that the SHAM converges faster for -G'(0) compared to F'(0). This is due to the differences in the level of nonlinearity of the equations of F(η) and G(η). For the ISHAM however, convergence has not been affected by this difference in the nonlinearity of functions. This there appears to be an added advantage of the ISHAM over the SHAM.

Table 7 gives a comparison of the convergence rate of the ISHAM and the SHAM versus the numerical solutions for F'(0) for different values M when K = 1. Convergence to the numerical solution is achieved at 2^nd order of the ISHAM and at the 8^th order of the SHAM. The absolute errors are shown in Table 8. In Table 9 the ISHAM and the SHAM solutions are compared with the numerical solutions for different values M when K = 1. An increase in the tangential shear stress is observed with an increase in M. The absolute errors are shown in Table 10. The nonlinearity of the equations has no effect on the convergence of the ISHAM.

Table 11 shows the heat transfer coefficient -Θ'(0) at different orders of the ISHAM compared against numerical results for different values of M, Pr, and Ec when K is fixed. For all cases convergence of the ISHAM approximate solutions to the numerical solution is achieved at 2^nd orders of approximations. For a fixed value of K, increasing M, Pr, and Ec decreases -Θ'(0). For all results with CPU run times, the ISHAM has shown great computer efficiency as well. In all simulations the run times are significantly lower with the ISHAM than with the SHAM.

Figures 1-2 show the effect of M and K on the radial and axial velocity profiles respectively. Increasing M reduces the radial component of the velocity while increasing K enhances F. The axial velocity H(η) increases with the magnetic parameter but decreases when K is increased. There is excellent agreement between the second order ISHAM solutions for F(η) and H(η) and thenumerical result.

The tangential velocity component and the temperature profiles are show in Figures 3-4 for different values of M and K respectively. An increase in M reduces G(η) while enhancing the Θ(η). When K values are increased both G(η) and Θ(η) decrease.

In Figure 5 temperature profiles are presented for varying values of Pr and Ec. The temperature decreases with increasing Prandtl numbers while an increase in the Eckert number enhances the temperature.

6 Conclusion

A novel approach for accelerating the convergence of the spectral homotopyanalysis method that is used to solve nonlinear equations in science and engineering has been proposed and applied successfully to the nonlinear system of equations governing the Reiner-Rivlin fluid in with Joule heating and viscousdissipation. The primary objective of the algorithm is to improve the initial approximate solution. The improved approximations are then used in the algorithm of the spectral-homotopy analysis method to reduce the number of iterations required to achieve convergence and better accuracy. The shear stresses in the radial and azimuthal directions were computed and the corresponding absolute errors determined. Convergence to the numerical solutions of the ISHAM approximate solutions was achieved at the 2^nd orders for all flow parameters while the SHAM converged at the 8^th order for some of the flow parameters.

The effects of flow parameters was also investigated for the radial and tangential shear stresses for both the Newtonian (K = 0) and non-Newtonian (K ≠ 0) cases. For the Newtonian case, increasing M reduces F'(0) and enhances -G'(0) while in the non-Newtonian case increasing M enhances both F'(0) and -G'(0).

The effect of K and M was determined and it was observed that an increase in K results in an increase in F(η), and a decrease in H(η), G(η) and Θ(η) while increasing M increases H(η) and Θ(η) while both F(η) and G(η) decreases. Θ(η) decreased with an increase in the Ec and decreased with an increase in Pr. The success of the ISHAM in solving the non-linear equations governing the von Kármán flow of an electrically conducting non-Newtonian Reiner-Rivlin fluid in the presence of viscous dissipation, Joule heating and heat transfer proves that the ISHAM fits as a newly improved method of solution that can be used to solve non-linear problems arising in science and engineering.

Acknowledgements. The authors wish to acknowledge financial support from the University of KwaZulu-Natal and the National Research Foundation (NRF).

Received: 04/I/11.

Accepted: 10/XI/11.

#CAM-311/11.

• [1] M.A. Abdou, New Analytical Solution of Von Karman swirling viscous flow Acta. Appl. Math., 111 (2010), 7-13.
• [2] P.D. Ariel, The homotopy perturbation method and analytical solution of the problem of flow past a rotating disk Computers and Mathematics with Applications, 58 (2009), 2504-2513.
• [3] A. Arikoglu, I. Ozkol and G. Komurgoz, Effect of slip on entropy generation in a single rotating disk in MHD flow Applied Energy, 85 (2008), 1225-1236.
• [4] H.A. Attia, The effect of ion slip on the flow of Reiner-Rivlin fluid due to a rotating disk with heat transfer J. of Mechanical Science and Technology, 21 (2007), 174-183.
• [5] H.A. Attia, Rotating disk flow and heat transfer through a porous medium of a non-Newtonian fluid with suction and injection Communications in Nonlinear Science and Numerical Simulation, 13 (2008), 1571-1580.
• [6] H.A. Attia and M.E.S. Ahmed, Non-Newtonian conducting fluid flow and heat transfer due to a ratating disk ANZIAM Journal, 46 (2004), 237-248.
• [7] H.A. Attia, Ion-slip effect on the flow due to a rotating disk The Arabian J. for Sc. and Eng., 29 (2004), 165-172.
• [8] H.A. Attia, Rotating Disk Flow and Heat Transfer of a Conducting Non-Newtonian Fluid with Suction-Injection and Ohmic Heating J. of the Braz. Soc. of Mech. Sci. & Eng., (XXIX) (2007).
• [9] H.A. Attia, Steady Flow over a Rotating Disk in Porous Medium with Heat Transfer Nonlinear Analysis: Modelling and Control, 14 (2009), 21-26.
• [10] H.A. Attia, Rotating disk flow with heat transfer of a non-Newtonian fluid inporous medium Turk. J. Phys., 30 (2006), 103-108.
• [11] H.A. Attia, On the effectiveness of uniform suction and injection on steady rotating disk flow in porous medium with heat transfer Turkish J. Eng. Env. Sci., 30 (2006), 231-236.
• [12] E.R. Benton, On the flow due to a rotating disk J. Fluid Mech., 24 (1966), 781-800.
• [13] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral Methods in Fluid Dynamics, Springer-Verlag, Berlin (1988).
• [14] C.S. Chien and Y.T. Shih, A cubic Hermite finite element-continuation method for numerical solutions of the von Kármán equations Appl. Math. and Comp., 209 (2009), 356-368.
• [15] W.G. Cochran, The flow due to a rotating disk Math. Proc. Cambridge Philos. Soc., 30 (1934), 365-375.
• [16] W.S. Don and A. Solomonoff, Accuracy and speed in computing the Chebyshev Collocation Derivative SIAM J. Sci. Comput., 16 (1995), 1253-1268.
• [17] A. El-Nahhas, Analytic approximations for von Kármán swirling flow Proc.Pakistan Acad. Sci., 44(3) (2007), 181-187.
• [18] P. Hatzikonstantinou, Magnetic and viscous effects on a liquid metal flow due to a rotating disk Astrophysics and Space Science, 161 (1989), 17-25.
• [19] N. Kelson and A. Desseaux, Note on porous rotating disk flow ANZIAM Journal, 42(E) (2000), C837-C855.
• [20] S.J. Liao, Beyond perturbation: Introduction to homotopy analysis methodChapman & Hall/CRC Press (2003).
• [21] Z.G Makukula, P. Sibanda and S.S. Motsa, On new solutions for heat transfer in a visco-elastic fluid between parallel plates, International journal of mathematical models and methods in applied sciences, 4(4) (2010), 221-230.
• [22] Z.G Makukula, P. Sibanda and S.S. Motsa, A Novel Numerical Technique for Two-Dimensional Laminar Flow between Two Moving Porous Walls, Mathematical Problems in Engineering, 2010 (2010), Article ID 528956, 15 pages doi: 10.1155/2010/528956.
• [23] Z.G Makukula, P. Sibanda and S.S. Motsa, A Note on the Solution of the Von Kármán Equations Using Series and Chebyshev Spectral Methods, Boundary Value Problems, 2010 (2010), Article ID 471793, 17 pages doi: 10.1155/2010/471793.
• [24] M. Miklavcic and C.Y. Wang, The flow due to a rough rotating disk Z. Angew. Math. Phys., 54 (2004), 1-12.
• [25] S.S. Motsa, P. Sibanda and S. Shateyi, A new spectral-homotopy analysis method for solving a nonlinear second order BVP, Communications in Nonlinear Science and Numerical Simulation, 15 (2010), 2293-2302.
• [26] S.S. Motsa, P. Sibanda, F.G. Awad and S. Shateyi, A new spectral-homotopy analysis method for the MHD Jeffery-Hamel problem Computers & Fluids, 39 (2010), 1219-1225.
• [27] S.S. Motsa and P. Sibanda, On the solution of MHD flow over a nonlinear stretching sheet by an efficient semi-analytical technique Int. J. Numer. Meth. Fluids, (2011), doi: 10.1002/fld.
• [28] S.S. Motsa, G.T. Marewo, P. Sibanda and S. Shateyi, An improved spectral homotopy analysis method for solving boundary layer problems Boundary Value Problems, 3 (2011), doi: 10.1186/1687-2770-2011-3, 11 pages.
• [29] E. Osalusi, Effects of thermal radiation on MHD and slip flow over a porous rotating disk with variable properties Romanian Journal Physics, 52 (2007),217-229.
• [30] E. Osalusi and P. Sibanda, On variable laminar convective flow properties due to a porous rotating disk in a magnetic field Romanian Journal Physics, 51 (2006), 937-950.
• [31] M.M. Rashidi and S.A.M. Pour, A novel analytical solution of steady flow over a rotating disk in porous medium with heat transfer by DTM-padé African J. of Math. and Comp. Sc. Research, 3 (2010), 93-100.
• [32] M.M. Rashidi and H. Shahmohamadi, Analytical solution of three-dimensional Navier-Stokes equations for the flow near an infinite rotating disk Communications in Nonlinear Science and Numerical Simulation, 14 (2009), 2999-3006.
• [33] B. Sahoo, Effects of partial slip, viscous dissipation and Joule heating on von Kármán flow and heat transfer of an electrically conducting non-Newtonian fluid Communications in Nonlinear Science and Numerical Simulation, 14 (2009), 2982-2998.
• [34] B.K. Sharma, A.K. Jha and R.C. Chaudhary, MHD forced flow of a conducting viscous fluid through a porous medium induced by an impervious rotating disk Rom. Journ. Phys., 52 (2007), 73-84.
• [35] P. Sibanda and O.D. Makinde, On steady MHD flow and heat transfer past a rotating disk in a porous medium with ohmic heating and viscous dissipation Int. J. of Num. Methods for Heat & Fluid Flow, 20 (2010), 269-285.
• [36] L.N. Trefethen, Spectral Methods in MATLAB SIAM (2000).
• [37] M. Turkyilmazoglu, Purely analytic solutions of magnetohydrodynamic swirling boundary layer flow over a porous rotating disk Computers & Fluids, 39 (2010), 793-799.
• [38] M. Turkyilmazoglu, Purely analytic solutions of the compressible boundary layer flow due to a porous rotating disk with heat transfer Physics of fluids, 21 (2009), 106104-12.
• [39] T. von Kármán, U berlaminare und turbulence reibung ZAMM, 1 (1921), 233-252.
• [40] C. Yang and S.J. Liao, On the explicit, purely analytic solution of von Kármán swirling viscous flow Comm. in Nonlinear Science and Num. Simulation, 11 (2006), 83-93.

*

Corresponding author

Publication Dates

• Publication in this collection
  26 Apr 2012
• Date of issue
  2012

History

• Received
  04 Jan 2011
• Accepted
  10 Nov 2011