Abstract

The governing nonlinear equation for unidirectional flow of a Sisko fluid in a cylindrical tube due to translation of the tube wall is modelled in cylindrical polar coordinates. The exact steady-state solution for the nonlinear problem is obtained. The reduction of the nonlinear initial value problem is carried out by using a similarity transformation. The partial differential equation is transformed into an ordinary differential equation, which is integrated numerically taking into account the influence of the exponent n and the material parameter b of the Sisko fluid. The initial approximation for the fluid velocity on the axis of the cylinder is obtained by matching inner and outer expansions for the fluid velocity. A comparison of the velocity, vorticity, and shear stress of Newtonian and Sisko fluids is presented.

1. Introduction

Over the past few decades, use of the Newtonian fluid model to analyze and predict the behaviour of many real fluids has been extensively adopted in industry. However, the flow characteristics of many real fluids have been found to be quite different from those of the Newtonian fluid and hence researchers have proposed many non-Newtonian fluid models to explain the deviation in the behaviour of real fluids from that of the Newtonian fluid. Several rheological models of non-Newtonian fluids have been proposed to represent the viscosity function of these fluids. Amongst these models is the Sisko model [1], which is the most suitable for the flow of greases. The appropriateness of the Sisko model has been successfully extended to the shear thinning rheological behaviour of concentrated non-Newtonian slurries [2]. Some polymeric suspensions such as waterbone coating are known to be non-Newtonian in nature and follow the Sisko model [3]. The viscosity of such coatings depends on the shear rate and the strain history. Many rheological fluids such as drilling fluids and cement slurries without yield stress also obey this model.

The properties of Sisko fluids have been investigated by the study of a range of problems. In [4] the problem of a Sisko fluid in Taylor’s scraping problem has been considered and magnetohydrodynamic [MHD] peristaltic motion of a Sisko fluid in symmetric and asymmetric channels has been considered in [5]. In [6] a Lie group analysis of the boundary layer equations for a Sisko fluid has been performed, and thin film flow of non-Newtonian and second grade fluids on a moving belt has been analysed in [7, 8]. Flow of a Sisko fluid in a porous medium has also been investigated. Solutions for MHD flow have been obtained in [9] and an analysis of heat transfer of an MHD flow in a porous medium has been performed in [10]. Stokes’ first problem for a rotating Sisko fluid with porous space has been studied in [11]. The Rayleigh problem has been investigated for a rotating Sisko fluid in [12] and for an MHD Sisko fluid in [13], while Stokes’ first problem for a Sisko fluid over a porous wall has been considered in [14].

There have been several investigations of non-Newtonian and Sisko fluids in cylindrical geometry. Steady flow and heat transfer of a Sisko fluid in an annular pipe have been investigated in [15]. Exact solutions for rotating flows of a generalised Burgers’ fluid have been derived in cylindrical geometry in [16] and exact solutions for the helical flow of a generalised Oldroyd-B fluid in a circular cylinder have also been obtained in [17]. The unsteady flow of an MHD Sisko fluid between two concentric tubes due to a prescribed pressure gradient along the tube, with the tube walls at rest, has been investigated in [18]. The problem was solved numerically by using the fourth-order Runge-Kutta method. More recently the steady flow and heat transfer of an MHD fluid in the porous space between concentric tubes have been considered in [19]. The flow was due to the motion of the outer cylinder and a constant pressure gradient along the tube. An analytical solution was derived using homotopy analysis and a numerical solution was obtained by an iterative method.

In this paper we consider the flow of a Sisko fluid in a cylindrical tube due to the translation of the tube wall parallel to the axis of the tube. The velocity of the wall is not prescribed but is determined from a similarity solution of the partial differential equation. There is no inner tube on which a no-slip boundary condition could be imposed as in [16, 17]. The flow is unsteady and the fluid velocity on the axis of the tube is estimated by matching inner and outer expansions for the fluid velocity. Using this initial estimate a numerical solution for the fluid velocity is obtained by a shooting method. This numerical method for a Sisko fluid is a new feature of the paper.

In the reduction of a partial differential equation subject to initial and boundary conditions to an ordinary differential equation with boundary conditions by a similarity transformation, the initial and boundary conditions cannot be arbitrary because they must be expressible in terms of the similarity variable. The value of such a transformation is the great simplification achieved by the reduction of a partial differential equation to an ordinary differential equation. The present investigation is an extension to cylindrical geometry of the two-dimensional problem of the flow induced by an infinite sliding solid plate on a half-space of viscous fluid [20]. When the plate is impulsively set in motion with constant speed , the flow is referred to as Stokes’ first problem or the Rayleigh problem. The analytical solution was used by Rayleigh as a model to study the diffusion of vorticity in a boundary layer on a flat plate. In the two-dimensional flow the plate velocity cannot be arbitrary for a similarity solution to exist but must be a power law of time . If the velocity of the plate is proportional to then the applied stress on the plate which induces the flow is constant, while if it is proportional to , the acceleration of the plate is constant. Other power laws can be considered leading to numerical solutions. We will find that, for a similarity solution to exist for a Sisko fluid in a cylindrical tube undergoing translation, the velocity of the tube wall must depend on time in a determined way. The initial velocity of the Sisko fluid across the tube cannot be arbitrary but must have a -shaped profile. Although these conditions would be difficult to realise in practice, the physical relevance of the similarity transformation is that it does yield a model to investigate the evolution with time of a Sisko fluid in a tube undergoing translation and to study the diffusion of its vorticity and shear stress from the translating wall to the axis of the tube.

There are several ways to derive similarity solutions of partial differential equations. We will derive the similarity solution by first obtaining the Lie point symmetries of the partial differential equation. This is a powerful systematic method which does not assume a form for the solution. Only one Lie point symmetry of the partial differential equation will be used which presents the possibility of other forms of solution for different boundary conditions. Other methods could be applied which do not require a knowledge of Lie group analysis of differential equations, for example, the approach of Dresner [21].

Although it is difficult to obtain exact solutions of the equations of motion of a non-Newtonian fluid, travelling wave and similarity solutions of nonlinear equations are desirable as such solutions play a very important role in the study of nonlinear wave and fluid flow phenomena. The analytical solutions, if available, facilitate the verification of numerical solvers and are also helpful in the stability analysis of solutions. In the literature, there are very few analytical solutions for non-Newtonian fluids. It is due to the fact that the governing equations of such fluids are much more complicated and of higher order than the Navier-Stokes equations. Unlike the Navier-Stokes equations which have nonlinear terms only in the inertia term, the equations for a non-Newtonian fluid have higher order nonlinear terms in the viscous term. Although an analytical solution may not be derived, it may be possible to reduce the partial differential equations to ordinary differential equations. With this background, the investigation of unsteady flow of a Sisko fluid in a cylindrical tube subject to initial and boundary conditions is carried out in the present study.

In this paper we concentrate on the reduction and numerical solution of the partial differential equation for the unsteady flow of a Sisko fluid in a cylindrical tube. This partial differential equation is transformed into an ordinary differential equation by using one of the Lie point symmetries of the partial differential equation. Numerical solutions of the ordinary differential equation are derived for values of the exponent corresponding to a shear thinning, Newtonian, and shear thickening fluid and for a fixed value of the material parameter .

The underlying physical process that the problem seeks to clarify is diffusion in a Sisko fluid. It is an ideal problem for investigation of this process. The diffusion of velocity, vorticity, and shear stress from the wall to the axis of the tube due to the translation of the wall will be studied. The process will be illustrated by computer generated graphs.

A gap in the literature for Sisko fluids which this investigation attempts to fill is the extension of Stokes’ first problem (Rayleigh problem), from a flat plate set in motion to the wall of a cylindrical tube set in motion. For a Newtonian fluid Stokes’ first problem yielded important insights into boundary layers.

There have been investigations of steady and unsteady flow of a Sisko fluid between concentric tubes [15, 18] with the no-slip boundary condition on the inner tube. A novel feature of the present study is the absence of the inner tube for unsteady flow driven by the translation of the tube wall. A matching procedure used to obtain an estimate for the boundary condition on the axis of the cylinder and the shooting method take the place of the no-slip boundary condition. The solution is possible because a similarity transformation is found that not only reduces the partial differential equation to an ordinary differential equation, but also determines the initial condition in the form of a -shaped velocity profile that has to be imposed.

The remainder of the paper is organized as follows. Section 2 deals with the formulation of the nonlinear initial boundary value problem. In Section 3 steady-state solutions are investigated, while in Section 4 the reduction of the partial differential equation and the formulation of the problem in terms of similarity variables are given. In Section 5 the problem is reformulated as a boundary value problem suitable for numerical computation. In Section 6 the numerical results are presented and discussed. Finally concluding remarks are made in Section 7.

2. Problem Formulation

Consider the unsteady unidirectional flow of an incompressible Sisko fluid in a circular cylinder parallel to the axis of the cylinder (see Figure 1). Cylindrical polar coordinates are chosen with the -axis along the axis of the cylinder. We assume the velocity, the stress fields, and pressure are of the form The fluid flow is generated by the translation of the wall of the cylindrical tube and not by a pressure gradient along the tube. The incompressibility condition is identically satisfied. The fluid flow is illustrated in Figure 1.

The Cauchy stress tensor for a Sisko fluid [7, 9] is In 3 and 4, is the pressure, is the identity tensor, is the trace-free nonisotropic part of the stress tensor, is the first Rivlin-Ericksen tensor, and , , and are positive material constants defined differently for different fluids.

The substitution of 1 into 3 yields the nonzero component of stress The body force due to gravity is neglected. The -component of the momentum balance equation in the absence of body forces gives The -component of the momentum balance equation shows that is independent of and the -component is identically satisfied. From 5 and 6, we obtain The boundary conditions are no-slip of the viscous fluid at the wall of the cylinder, , and symmetry of the velocity at : where is the velocity of the wall of the cylinder; is the reference velocity; and is dimensionless. The initial condition is where is the radius of the cylinder. The functions and are as yet arbitrary functions. We choose the reference velocity to be the velocity of the wall of the cylinder at , so that

We now introduce dimensionless variables and define the dimensionless quantities , , , , and by The parameter depends on . Expressed in nondimensional form and suppressing the asterisks, the problem is to solve the nonlinear diffusion equation subject to the boundary conditions where and to the initial condition where We will consider the solution of the problem for When 18 is satisfied, 12 becomes in expanded form where subscripts denote partial differentiation.

3. Steady-State Solution

For the steady-state solution . We set in 12 assuming 18, which gives Integration of 20 with respect to and imposition of the boundary condition 14 result in Since and imposing the boundary condition 16 that , gives which is the constant velocity of the wall of the cylinder.

4. Lie Symmetry Analysis and Similarity Solution

The Lie approach to derive similarity solutions of partial differential equations is widely used and has been applied to problems in fluid mechanics. We refer the reader to [22, 23] in which the prolongation formulae and method are given in detail. In this paper we present the main results. The Lie point symmetry generator of 19 is of the form There are four cases to consider.

Case (i) .

The coefficient functions in 23 satisfy the determining equations The Lie point symmetries are Case (ii) .

The determining equations are The Lie point symmetries 25 are again derived.

Case (iii) .

This describes a Newtonian fluid and the underlying linear parabolic equation (19 with ) has the Lie point symmetries where satisfies the linear partial differential equation itself. In this case there are infinitely many Lie point symmetries.

Case (iv) .

For   19 again has infinitely many Lie point symmetries; namely, where satisfies the linear partial differential equation 19 for .

For all four cases one of the Lie point symmetries is the scaling symmetry . Now is a group invariant solution generated by the Lie point symmetry provided that is, provided satisfies the first order linear partial differential equation The differential equations of the characteristic curves of 30 are The first pair of terms in 31 gives where is a constant. The second and third terms in 31 give where is a constant. The general solution of 30 is where is an arbitrary function. Since the similarity solution for generated by is of the form Substituting 35 into 19 reduces 19 to the ordinary differential equation where dash denotes differentiation with respect to .

The similarity form 35 can also be found by using the similarity ansatz and determining the constants , , and by substituting 37 into the partial differential equation 19. The Lie group analysis is of value, however, because it has shown that there may be other useful solutions of the partial differential equation generated by a linear combination of its Lie point symmetries.

By considering the similarity form 35, we take as initial condition which satisfies the boundary conditions 16 and 17, that , and also assumption 18. Expressed in terms of the similarity variables 35, 38 becomes

Next consider the boundary conditions. Expressed in terms of the similarity variables, 14 becomes The boundary condition, 13 and 15, becomes The condition 15 that is satisfied because of 39. Equation 41 gives the velocity of the wall of the cylinder after the solution for has been derived. The velocity of the wall of the cylinder cannot be prescribed arbitrarily in this similarity solution.

The problem is to solve the ordinary differential equation 36 for subject to the boundary condition 40 at and 39 at . The boundary condition 40 is not in a form suitable for numerical computation. In Section 5 the boundary value problem will be reformulated in a way suitable for numerical computation.

5. Numerical Method

Consider first the boundary conditions at ,  , which corresponds to . From 35, the -component of the velocity is Since the fluid velocity must be finite on the axis of the cylinder, , it follows that where is a constant which will be nonzero because the viscous fluid is set in motion by the translation of the wall of the cylinder. We therefore introduce the new function defined by Then and satisfies the boundary condition The constant is still to be determined. Also from 39, and the boundary condition 14 becomes

Consider next the boundary condition at which corresponds to and . From 39, and therefore Thus the boundary conditions at and can be expressed in terms of in a form suitable for numerical computation. The problem will therefore be formulated in terms of instead of and the ordinary differential equation 36 will be expressed in terms of .

The problem is to solve the differential equation for subject to the boundary conditions where the constant has still to be determined.

A shooting method will be used which will require an initial estimate for . To obtain an approximate expression for we first consider the asymptotic expansion for as . From 49, we consider the expansion as . Substituting 52 into 50 and equating the coefficients of like powers of give and therefore, for , as . From 11, depends on through . It first depends explicitly on at the term in the expansion 54. To obtain an estimate for it is sufficient to consider the first two terms in 53 which give a simple result. It follows from 45 that as . In Figure 2 the two-term expansion for is plotted against for and at . We see that decreases steadily from as decreases from , reaches a minimum value at , and then steadily increases and tends to infinity as . Using only the first two terms in 55, we have The minimum value of occurs at where We match 58 with an expansion for small . Since from 48, , a Taylor expansion of at gives the two-term expansion as and therefore, from 45, since , as . We take as a first approximation The matching condition is and therefore, from 58 and 61,

The problem, 50 and 51, was treated as an initial value problem with initial conditions and . The initial estimate for was 63. A shooting method was utilized and the value of was adjusted until the boundary condition was attained to sufficient accuracy. A comparison of the initial estimate for   63 and the final value for obtained with the shooting method is given in Section 6.

Once has been calculated, for is given by From the no-slip boundary condition, the velocity of the wall of the cylinder, , is obtained by setting in 64: For small values of time, the asymptotic expansion for large values of applies and from 54, as . For large values of time, which corresponds to small values of , it follows from 60 that as , where , the final value of derived in the shooting method, depends on . The wall velocity cannot be prescribed. It increases approximately linearly for small , while for large it increases as .

In cylindrical polar coordinates the vorticity is given by where is the unit base vector parallel to the -coordinate line in the direction of increase of . From 47 and therefore, assuming 18, the magnitude of the vorticity, , is given by Since , it follows that and since , we have Consider the vorticity at the wall of the cylinder . From expansion 54 as it follows that as . We see that the magnitude of the vorticity at the wall decreases approximately linearly for small values of time.

Finally, consider the shear stress given by 5. The characteristic quantities are defined by 11. The characteristic stress is defined by Expressed in dimensionless form and suppressing the asterisks, the shear stress is Using 47, 75 becomes Since and , it follows that The wall shear stress is given by The expansion of for small values of time is obtained from the expansion 54 as . We find that as . The wall shear stress, like the vorticity at the wall, decreases approximately linearly with time for small values of time.

6. Results and Discussion

When comparing the evolution of shear thinning, Newtonian, and shear thickening fluids, the same characteristic time and characteristic stress must be used and the dimensionless parameter , where needs to be specified. We will make a preliminary comparison by assuming that , , and are the same for the fluids under consideration. The characteristic time and characteristic stress are therefore the same for all fluids considered. Furthermore we will assume that the ratio in the SI units used so that is also the same for all fluids considered. To carry out a more detailed comparison, the parameters , , , , and would need to be specified.

Using  Matlab ode23s, the shooting method outlined in Section 5 gave physically acceptable results for shear thickening fluids, Newtonian fluids, and for shear thinning fluids with values down to . For we did not obtain physically acceptable results with for . To obtain feasible results we resorted to  Matlab bvp4c for . We are of the opinion that there may be a type change in the solution at some small value of . The  Matlab ordinary differential equation solvers will generally be better than anything one would program oneself. They are able to estimate the error in the solution at each time step and decide whether or not the time step is too large (error too high) or too small (inefficient).  Matlab routine  bvp4c is an adaptive Lobatto iterative scheme. Boundary value problems arise in most diverse forms. Just about any boundary value problem can be formulated for solution with  bvp4c. The first step is to write the ordinary differential equations as a system of first order ordinary differential equations. Details of the solution method can be found in [24] and the references therein.

The initial estimate for , 63, was used for all values of . In Figure 3 the final value obtained for using the shooting method is plotted against for a range of values of . The graphs for , , and show that starts to increase rapidly for below a critical value. This critical value of increases slowly with . It can be expected that the graphs for and would also show this property if the solution had been taken to smaller values of . This rapid increase in demonstrates the difficulty in obtaining a numerical solution for very small values of .

In Figure 2, which illustrates the matching principle, and Figures 4 to 7 which investigate the physical properties of the solution, the value of used to perform the numerical calculations is . In Table 1 the initial estimate for given by 63, the final value obtained for using the shooting method and the relative error percentage for a range of values of are listed for . We see that 63 underestimates for very small values of , but otherwise it overestimates . The relative error percentage in the initial estimate for increases rapidly as decreases below and increases steadily with for . The shooting method with the initial estimate 63 still converges when the relative error in the initial estimate for exceeds .

The final value of has important physical significance. From 45 and 46 the fluid velocity on the axis of the cylinder, , is and from 67 we see that for large values of time the wall velocity is given by For large values of time the velocity profile becomes almost independent of and the wall velocity tends to the value of the fluid velocity on the axis.

Consider next the fluid velocity for . In Figure 4   is plotted against for shear thinning fluids with and , for a Newtonian fluid with , and for a shear thickening fluid with . For all values of considered, the initial -shaped velocity profile evolves into a smooth -shaped profile for which gradually flattens out as time increases. Because of the no-slip boundary condition, fluid velocity is generated at the wall of the cylinder and diffuses radially from the wall towards the axis of the cylinder which produces the flattening of the velocity profile.

In Figures 5(a), 5(b), and 5(c) the scaled velocity is plotted against for , and . Figure 5 clearly shows the diffusion of velocity in the radial direction from the cylindrical wall to the axis. For the shear thinning fluid the velocity profile in the neighbourhood of the axis of the tube flattens out at a much earlier time than that for the Newtonian fluid and shear thickening fluid. We see clearly from the graphs for that the fluid velocity on the axis is less for the shear thickening fluid with than for the Newtonian and shear thinning fluids with and , respectively. In Figure 5(d) the velocity of the cylindrical wall, , is plotted against for , and . Because of the no-slip boundary condition at the wall, . The wall velocity is determined by the similarity solution and cannot be prescribed. We see that increases approximately linearly for small values of in accordance with expansion 65. For large values of time, behaves approximately like . From Table 1 for , , and and for large time, respectively, in agreement with Figure 5(d). For large time the wall velocity decreases as increases from values for a shear thinning fluid to a shear thickening fluid, except for a small range of which may be due to a small computational error.

Consider next the magnitude of the vorticity, , given by 70. In Figures 6(a), 6(b), and 6(c),   is plotted against for a range of values of time for , and . The vorticity generated at the wall due to the no-slip boundary condition diffuses inwards towards the axis of the tube. The initial vorticity profile, given by 72, has the form of a right angle. The vorticity decreases steadily with time and vanishes as when the fluid velocity becomes uniform across the tube. In Figure 6(d) the magnitude of the vorticity at the wall, , is plotted against for , and . At , and it decreases approximately linearly with time for small values of time in agreement with expansion 73. It vanishes as . Since the viscosity of the shear thickening fluid is increased by the shear at the wall, the vorticity generated at the wall is greater than that of the Newtonian fluid. The viscosity of the shear thinning fluid is decreased by the wall shear stress and the vorticity generated at the wall is less than that of the Newtonian fluid.

Finally, consider the shear stress given by 76. In Figures 7(a), 7(b), and 7(c), the shear stress, scaled by its initial value, is plotted against for a range of values of and for , and . The initial profile of the shear stress is given by 77 and has the form of a right angle. For the shear stress is greatest at the wall of the tube where the velocity gradient is greatest and decreases steadily to zero on the axis of the tube where the velocity gradient vanishes. In Figure 7(d), the scaled wall shear stress is plotted against for , , and . The decrease of the wall shear stress with time is approximately linear for small values of time in agreement with expansion 79. From 78, as . Figure 7(d) and 84 show that the wall shear stress decreases faster with time for a shear thickening fluid than for a shear thinning fluid. A significant wall shear stress needs to be maintained longer for a shear thinning fluid than for a shear thickening fluid.

The results in the literature for a Sisko fluid generally do not consider vorticity or shear stress. However, the axial velocity profiles in the annular flow of a Sisko fluid between concentric tubes driven by a negative pressure gradient have been obtained in [18] and it is of interest to compare their results with our results for the axial velocity profiles of a Sisko fluid in a cylindrical tube undergoing translation. The two problems are very different and are almost the opposite of each other. For the problem considered in [18] the initial axial velocity profile is flat with no-slip boundary conditions on each tube wall. This flat velocity profile evolves into the familiar parabolic-like profile for flow between concentric tubes. The maximum axial velocity decreases as the material parameter increases. The axial velocity profile stays flat in the central region for a longer time for a shear thinning fluid. This compares with the evolution of our initial -shaped velocity profile into a flat profile. From Figure 3, the velocity on the axis of the tube, , increases as the material parameter increases. The velocity profile in the central region becomes flat at an earlier time for a shear thinning fluid. The reason for the different behaviour in the two problems is the small rate-of-shear in the nearly flat velocity profile compared with the large rate-of-shear in - and -shaped velocity profiles.

7. Conclusions

A preliminary investigation of the properties of the numerical solution was undertaken. Since we took , , and to be the same for all the fluids considered, the investigation mainly demonstrated the effect of the power law exponent, , on the numerical solution. As we expected, the velocity of the tube wall and the initial velocity profile of the Sisko fluid across the tube are determined by the similarity solution and cannot be assigned arbitrarily. The initial -shaped velocity profile 38 may be difficult to realise in practice. However the similarity solution provided an instructive model to investigate the time evolution of a Sisko fluid and the diffusion of its vorticity and shear stress from a translating tube wall to the axis of the tube.

The shooting method described in Section 5 was applicable for a large range of values of . It gave physically acceptable results for shear thickening and Newtonian fluids and for shear thinning fluids for values of as low as . This is a strong indication that the numerical solution is correct. If the solution for values of close to each other is being considered then instead of using the initial estimate 63 for the final value for from the previous solution could be used.

The  Matlab ODE suite is a collection of five user-friendly finite difference codes for solving initial value problems given by first-order systems of ordinary differential equations [25]. They are able to estimate the error in the solution at each time step and decide whether or not the time step is too large (error too high) or too small (inefficient). The code  ode23s is a triple of modified implicit Rosenbrock methods of orders and with error control for stiff systems. It advances from time step to time step with the second-order method (i.e., without local extrapolation) and controls the local error by taking the difference between the third- and second-order numerical solutions. Because it is a one-step solver, it may be more efficient than  ode15s at crude tolerances. The code  bvp4c is a finite difference code that implements the three-stage Lobatto formula. This is a collocation formula and the collocation polynomial provides a -continuous solution that is fourth-order accurate uniformly in interval . Mesh selection and error control are based on the residual of the continuous solution [24].

The physical mechanism in the problem is diffusion. Fluid velocity generated by the no-slip boundary condition when the wall of the cylinder is impulsively set in motion diffuses in the radial direction towards the axis of the cylinder. This causes the velocity profiles to flatten out and the vorticity and shear stress across the tube to steadily decrease and vanish as . The problem has some features in common with Stokes’ first problem [20] for flow induced in a half-space of viscous fluid when a wall is impulsively set in motion.

The graphs of the fluid velocity, vorticity, and wall shear stress within the fluid showed no significant differences in form for shear thickening, Newtonian, and shear thinning fluids, although the velocity profile near the axis of the tube for a shear thinning fluid flattened out at a much earlier time. Quantitative differences were observed which were most apparent on the axis of the tube and on the wall of the tube.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

S. Abelman, D. P. Mason, and F. M. Mahomed acknowledge support from the National Research Foundation of South Africa and the University of the Witwatersrand, Johannesburg. S. Abelman acknowledges Dr. Herven Abelman for drawing Figure 1 in the paper. Taha Aziz is thanked for his valuable comments on the original paper. The authors thank anonymous referees for constructive comments which have improved the paper. F. M. Mahomed is Visiting Professorial Fellow at UNSW for 2014.