# Self-Similar Unsteady Flow of a Sisko Fluid in a Cylindrical Tube Undergoing Translation

Authors: M. Khan; S. Abelman; D. P. Mason; F. M. Mahomed
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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 exponentn 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.


---

## Body

## 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 U
0, 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 t. If the velocity of the plate is proportional to t
1
/
2 then the applied stress on the plate which induces the flow is constant, while if it is proportional to t, 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 V-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 exponentn corresponding to a shear thinning, Newtonian, and shear thickening fluid and for a fixed value of the material parameter b.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 V-shaped velocity profile that has to be imposed.The remainder of the paper is organized as follows. Section2 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 Figure1). Cylindrical polar coordinates (
r
,
θ
,
z
) are chosen with the z-axis along the axis of the cylinder. We assume the velocity, the stress fields, and pressure are of the form
(1)
V
=
0,0
,
w
r
,
t
,
S
=
S
r
,
t
,
p
=
p
(
r
,
t
)
.
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
(2)
div
⁡
V
=
0
is identically satisfied. The fluid flow is illustrated in Figure 1.Figure 1
Unidirectional flow induced in a cylindrical tube of incompressible Sisko fluid when the wall of the cylinder is impulsively set in motion.The Cauchy stress tensor for a Sisko fluid [7, 9] is
(3)
T
=
-
p
I
+
S
,
S
=
a
+
b
1
2
tr
⁡

A
1
2
n
-
1
A
1
,
(4)
A
1
=
L
+
L
T
,
L
=
grad
V
.
In 3 and 4, p is the pressure, I is the identity tensor, S is the trace-free nonisotropic part of the stress tensor, A
1 is the first Rivlin-Ericksen tensor, and n, a, and b are positive material constants defined differently for different fluids.The substitution of1 into 3 yields the nonzero component of stress
(5)
S
r
z
=
a
+
b
∂
w
∂
r
n
-
1
∂
w
∂
r
.
The body force due to gravity is neglected. The z-component of the momentum balance equation in the absence of body forces gives
(6)
ρ
∂
w
∂
t
=
1
r
∂
∂
r
r
S
r
z
.
The r-component of the momentum balance equation shows that p is independent of r and the θ-component is identically satisfied. From 5 and 6, we obtain
(7)
ρ
∂
w
∂
t
=
1
r
∂
∂
r
r
a
+
b
∂
w
∂
r
n
-
1
∂
w
∂
r
.
The boundary conditions are no-slip of the viscous fluid at the wall of the cylinder, r
=
R, and symmetry of the velocity at r
=
0:
(8)
w
R
,
t
=
U
0
V
t
,
t
≥
0
,
∂
w
∂
r
0
,
t
=
0
,
t
≥
0
,
where U
0
V
t is the velocity of the wall of the cylinder; U
0 is the reference velocity; and V
t is dimensionless. The initial condition is
(9)
w
r
,
0
=
W
r
,
0
≤
r
≤
R
,
where R is the radius of the cylinder. The functions V
t and W
r are as yet arbitrary functions. We choose the reference velocity U
0 to be the velocity of the wall of the cylinder at t
=
0, so that
(10)
V
0
=
1
,
W
1
=
U
0
.We now introduce dimensionless variables and define the dimensionless quantitiesw
*, t
*, r
*, b
*, and W
* by
(11)
w
=
U
0
w
*
,
t
=
ρ
R
2
a
t
*
,
r
=
R
r
*
,
b
=
a
R
U
0
n
-
1
b
*
,
W
r
=
U
0
W
*
r
*
.
The parameter b
* depends on n. Expressed in nondimensional form and suppressing the asterisks, the problem is to solve the nonlinear diffusion equation
(12)
∂
w
∂
t
=
1
r
∂
∂
r
r
1
+
b
∂
w
∂
r
n
-
1
∂
w
∂
r
,
subject to the boundary conditions
(13)
w
1
,
t
=
V
t
,
t
≥
0
,
(14)
∂
w
∂
r
0
,
t
=
0
,
t
≥
0
,
where
(15)
V
0
=
1
,
and to the initial condition
(16)
w
r
,
0
=
W
r
,
0
≤
r
≤
1
,
where
(17)
W
1
=
1
.
We will consider the solution of the problem for
(18)
∂
w
∂
r
=
∂
w
∂
r
,
0
≤
r
≤
1
.
When 18 is satisfied, 12 becomes in expanded form
(19)
w
t
=
w
r
r
+
b
n
w
r
r
w
r
n
-
1
+
1
r
w
r
+
b
w
r
n
,
where subscripts denote partial differentiation.


## 3. Steady-State Solution

For the steady-state solutionw
=
w
r. We set w
t
=
0 in 12 assuming 18, which gives
(20)
d
d
r
r
w
r
+
b
w
r
n
=
0
.
Integration of 20 with respect to r and imposition of the boundary condition 14 result in
(21)
w
r
+
b
w
r
n
=
0
,
n
>
0
.
Since b
>
0 and imposing the boundary condition 16 that w
1
=
1, gives
(22)
w
r
=
constant
=
1
,
n
>
0
,
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
(23)
X
=
τ
t
,
r
,
w
∂
∂
t
+
ξ
t
,
r
,
w
∂
∂
r
+
η
t
,
r
,
w
∂
∂
w
.
There are four cases to consider.Case (i)n
≠
0,1
,
2.The coefficient functions in23 satisfy the determining equations
(24)
τ
r
=
τ
w
=
η
t
=
η
r
=
ξ
w
=
ξ
t
=
0
,
η
w
=
ξ
r
,
2
ξ
r
-
τ
t
=
0
,
η
w
=
τ
t
-
1
r
ξ
.
The Lie point symmetries are
(25)
X
1
=
∂
∂
t
,
X
2
=
∂
∂
w
,
X
3
=
2
t
∂
∂
t
+
r
∂
∂
r
+
w
∂
∂
w
.
Case (ii) n
=
2.The determining equations are(26)
τ
r
=
τ
w
=
ξ
w
=
0
,
η
w
=
3
ξ
r
-
τ
t
,
2
b
η
r
+
τ
t
=
2
ξ
r
,
-
2
ξ
r
r
+
1
r
η
w
-
2
r
ξ
r
-
1
r
2
ξ
+
1
r
τ
t
=
0
,
η
t
=
η
r
r
+
1
r
η
r
,
2
b
η
r
r
-
1
r
ξ
r
+
2
b
r
η
r
-
1
r
2
ξ
+
1
r
τ
t
+
ξ
t
=
0
.
The Lie point symmetries 25 are again derived.Case (iii)n
=
1.This describes a Newtonian fluid and the underlying linear parabolic equation (19 with n
=
1) has the Lie point symmetries
(27)
X
1
=
∂
∂
t
,
X
2
=
w
∂
∂
w
,
X
3
=
2
t
∂
∂
t
+
r
∂
∂
r
+
w
∂
∂
w
,
X
4
=
w
∂
∂
w
,
X
d
=
d
(
t
,
r
)
∂
∂
w
,
where d
(
t
,
r
) satisfies the linear partial differential equation itself. In this case there are infinitely many Lie point symmetries.Case (iv)n
=
0.Forn
=
0  19 again has infinitely many Lie point symmetries; namely,
(28)
X
1
=
∂
∂
t
,
X
2
=
w
∂
∂
w
,
X
3
=
2
t
∂
∂
t
+
r
∂
∂
r
+
w
∂
∂
w
,
X
4
=
w
+
b
r
∂
∂
w
,
X
d
=
d
(
t
,
r
)
∂
∂
w
,
where d
(
t
,
r
) satisfies the linear partial differential equation 19 for n
=
0.For all four cases one of the Lie point symmetries is the scaling symmetryX
3. Now w
=
W
(
r
,
t
) is a group invariant solution generated by the Lie point symmetry X
3 provided
(29)
X
3
w
-
W
r
,
t
|
w
=
W
=
0
,
that is, provided W
(
r
,
t
) satisfies the first order linear partial differential equation
(30)
2
t
∂
W
∂
t
+
r
∂
W
∂
r
=
W
.
The differential equations of the characteristic curves of 30 are
(31)
d
t
2
t
=
d
r
r
=
d
W
W
.
The first pair of terms in 31 gives
(32)
r
t
=
a
1
,
where a
1 is a constant. The second and third terms in 31 give
(33)
W
r
=
a
2
,
where a
2 is a constant. The general solution of 30 is
(34)
a
2
=
f
a
1
,
where f is an arbitrary function. Since w
=
W
r
,
t the similarity solution for w generated by X
3 is of the form
(35)
w
(
r
,
t
)
=
r
f
(
ζ
)
,
ζ
=
r
t
.
Substituting 35 into 19 reduces 19 to the ordinary differential equation
(36)
ζ
f
′
′
+
2
+
1
2
ζ
2
f
′
+
b
n
2
f
′
+
ζ
f
′
′
f
+
ζ
f
′
n
-
1
+
1
ζ
f
+
ζ
f
′
+
1
ζ
b
f
+
ζ
f
′
n
=
0
,
where dash denotes differentiation with respect to ζ.The similarity form35 can also be found by using the similarity ansatz
(37)
w
(
r
,
t
)
=
r
c
1
f
(
ζ
)
,
ζ
=
r
c
2
t
c
3
and determining the constants c
1, c
2, and c
3 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 form35, we take as initial condition
(38)
w
r
,
0
=
r
,
0
≤
r
≤
1
,
which satisfies the boundary conditions 16 and 17, that w
1,0
=
1, and also assumption 18. Expressed in terms of the similarity variables 35, 38 becomes
(39)
f
∞
=
1
.Next consider the boundary conditions. Expressed in terms of the similarity variables,14 becomes
(40)
∂
w
∂
r
0
,
t
=
lim
⁡
ζ
→
0
f
ζ
+
ζ
f
′
ζ
=
0
.
The boundary condition, 13 and 15, becomes
(41)
V
t
=
f
1
t
,
V
0
=
1
.
The condition 15 that V
(
0
)
=
1 is satisfied because of 39. Equation 41 gives the velocity of the wall of the cylinder after the solution for f
(
ζ
) 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 equation36 for f
(
ζ
) subject to the boundary condition 40 at ζ
=
0 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 atr
=
0,  t
>
0, which corresponds to ζ
=
0. From 35, the z-component of the velocity is
(42)
w
r
,
t
=
t
ζ
f
ζ
,
ζ
=
r
t
.
Since the fluid velocity must be finite on the axis of the cylinder, r
=
0, it follows that
(43)
lim
⁡
ζ
→
0
ζ
f
ζ
=
B
0
,
where B
0 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 g
(
ζ
) defined by
(44)
g
ζ
=
ζ
f
ζ
.
Then
(45)
w
r
,
t
=
t
g
ζ
,
ζ
=
r
t
,
and g
ζ satisfies the boundary condition
(46)
g
0
=
B
0
.
The constant B
0 is still to be determined. Also from 39,
(47)
∂
w
∂
r
r
,
t
=
d
g
d
ζ
,
and the boundary condition 14 becomes
(48)
d
g
d
ζ
0
=
0
.Consider next the boundary condition atζ
=
∞ which corresponds to 0
<
r
≤
1 and t
=
0. From 39, f
∞
=
1 and therefore
(49)
d
g
d
ζ
∞
=
1
.
Thus the boundary conditions at ζ
=
0 and ζ
=
∞ can be expressed in terms of g
(
ζ
) in a form suitable for numerical computation. The problem will therefore be formulated in terms of g
(
ζ
) instead of f
(
ζ
) and the ordinary differential equation 36 will be expressed in terms of g
(
ζ
).The problem is to solve the differential equation(50)
1
+
n
b
d
g
d
ζ
n
-
1
d
2
g
d
ζ
2
+
b
ζ
d
g
d
ζ
n
+
ζ
2
+
1
ζ
d
g
d
ζ
-
1
2
g
=
0
,
for g
ζ subject to the boundary conditions
(51)
g
0
=
B
0
,
d
g
d
ζ
0
=
0
,
d
g
d
ζ
∞
=
1
,
where the constant B
0 has still to be determined.A shooting method will be used which will require an initial estimate forB
0. To obtain an approximate expression for B
0 we first consider the asymptotic expansion for g
ζ as ζ
→
∞. From 49, we consider the expansion
(52)
g
ζ
=
ζ
+
A
0
+
A
1
ζ
+
A
2
ζ
2
+
A
3
ζ
3
+
A
4
ζ
4
+
O
1
ζ
5
as ζ
→
∞. Substituting 52 into 50 and equating the coefficients of like powers of ζ give
(53)
A
0
=
0
,
A
1
=
1
+
b
,
A
2
=
0
,
A
3
=
1
2
1
+
b
1
+
n
b
,
A
4
=
0
,
and therefore, for 0
<
n
<
∞,
(54)
g
ζ
=
ζ
+
1
+
b
ζ
+
1
+
b
1
+
n
b
2
ζ
3
+
O
1
ζ
5
,
as ζ
→
∞. From 11, g
ζ depends on n through b. It first depends explicitly on n at the O
ζ
-
3 term in the expansion 54. To obtain an estimate for B
0 it is sufficient to consider the first two terms in 53 which give a simple result. It follows from 45 that
(55)
w
O
r
,
t
=
r
+
1
+
b
t
r
+
O
t
2
r
3
,
as t
→
0. In Figure 2 the two-term expansion for w
O
(
r
,
t
) is plotted against r for 0
<
r
≤
1 and b
=
2 at t
=
0.05. We see that w
O
(
r
,
t
) decreases steadily from w
O
(
1
,
t
) as r decreases from r
=
1, reaches a minimum value at r
=
r
min
⁡, and then steadily increases and tends to infinity as r
→
0. Using only the first two terms in 55, we have
(56)
∂
w
O
∂
r
r
,
t
=
1
-
1
+
b
t
r
2
.
The minimum value of w
O
r
,
t occurs at r
=
r
min
⁡ where
(57)
r
min
⁡
=
1
+
b
1
/
2
t
,
(58)
w
O
r
min
⁡
,
t
=
2
1
+
b
1
/
2
t
.
We match 58 with an expansion for small r. Since from 48, g
′
0
=
0, a Taylor expansion of g
ζ at ζ
=
0 gives the two-term expansion
(59)
g
ζ
=
g
0
+
O
ζ
2
,
as ζ
→
0 and therefore, from 45, since g
0
=
B
0,
(60)
w
I
r
,
t
=
B
0
t
+
O
r
2
t
,
as r
→
0. We take as a first approximation
(61)
w
I
r
,
t
=
B
0
t
.
The matching condition is
(62)
w
I
r
min
⁡
,
t
=
w
O
r
min
⁡
,
t
,
and therefore, from 58 and 61,
(63)
B
0
=
2
1
+
b
1
/
2
.Figure 2
Matching principle. The two-term outer expansionw
O
(
r
,
t
) plotted against r for b
=
2 and t
=
0.05. The two-term inner expansion is w
I
(
r
,
t
)
=
B
0
t
+
O
(
r
2
).The problem,50 and 51, was treated as an initial value problem with initial conditions g
0
=
B
0 and g
′
0
=
0. The initial estimate for B
0 was 63. A shooting method was utilized and the value of B
0 was adjusted until the boundary condition g
′
∞
=
1 was attained to sufficient accuracy. A comparison of the initial estimate for B
0  63 and the final value for B
0 obtained with the shooting method is given in Section 6.Onceg
ζ has been calculated, w
r
,
t for t
>
0 is given by
(64)
w
r
,
t
=
t
g
r
t
,
0
≤
r
≤
1
.
From the no-slip boundary condition, the velocity of the wall of the cylinder, V
t, is obtained by setting r
=
1 in 64:
(65)
V
t
=
t
g
1
t
.
For small values of time, the asymptotic expansion for large values of ζ applies and from 54,
(66)
V
t
=
1
+
1
+
b
t
+
1
2
1
+
b
1
+
n
b
t
2
+
O
t
3
,
as t
→
0. For large values of time, which corresponds to small values of ζ, it follows from 60 that
(67)
V
t
=
B
0
t
+
O
1
t
,
as t
→
∞, where B
0, the final value of B
0 derived in the shooting method, depends on n. The wall velocity cannot be prescribed. It increases approximately linearly for small t, while for large t it increases as t.In cylindrical polar coordinates the vorticityω is given by
(68)
ω
=
curl
⁡
v
=
-
∂
w
∂
r
r
,
t
e
θ
,
where e
θ is the unit base vector parallel to the θ-coordinate line in the direction of increase of θ. From 47
(69)
ω
=
-
d
g
d
ζ
e
θ
,
and therefore, assuming 18, the magnitude of the vorticity, ω, is given by
(70)
ω
r
,
t
=
d
g
d
ζ
,
ζ
=
r
t
.
Since g
′
0
=
0, it follows that
(71)
ω
r
,
∞
=
0
,
0
≤
r
≤
1
,
and since g
′
∞
=
1, we have
(72)
ω
r
,
0
=
1
,
0
<
r
≤
1
.
Consider the vorticity at the wall of the cylinder r
=
1. From expansion 54 as ζ
→
∞ it follows that
(73)
ω
1
,
t
=
1
-
1
+
b
t
-
3
2
1
+
b
1
+
n
b
t
2
+
O
t
3
,
as t
→
0. We see that the magnitude of the vorticity at the wall decreases approximately linearly for small values of time.Finally, consider the shear stressS
r
z given by 5. The characteristic quantities are defined by 11. The characteristic stress is defined by
(74)
S
r
z
=
a
U
0
R
S
r
z
*
.
Expressed in dimensionless form and suppressing the asterisks, the shear stress is
(75)
S
r
z
r
,
t
=
1
+
b
∂
w
∂
r
n
-
1
∂
w
∂
r
.
Using 47, 75 becomes
(76)
S
r
z
r
,
t
=
1
+
b
d
g
d
ζ
n
-
1
d
g
d
ζ
,
ζ
=
r
t
.
Since g
′
∞
=
1 and g
′
0
=
0, it follows that
(77)
S
r
z
r
,
0
=
1
+
b
,
0
<
r
≤
1
,
S
r
z
r
,
∞
=
0
,
0
≤
r
≤
1
.
The wall shear stress is given by
(78)
S
r
z
1
,
t
=
1
+
b
d
g
d
ζ
1
t
n
-
1
d
g
d
ζ
1
t
.
The expansion of S
r
z
1
,
t for small values of time is obtained from the expansion 54 as ζ
→
∞. We find that
(79)
S
r
z
1
,
t
=
1
+
b
1
-
1
+
n
b
t
+
O
t
2
,
as t
→
0. 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 parameterb
*, where
(80)
b
*
=
b
a
U
0
R
n
-
1
,
needs to be specified. We will make a preliminary comparison by assuming that ρ, a, and b are the same for the fluids under consideration. The characteristic time a
/
ρ
R
2 and characteristic stress R
/
a
U
0 are therefore the same for all fluids considered. Furthermore we will assume that the ratio U
0
/
R
=
1 in the SI units used so that b
* is also the same for all fluids considered. To carry out a more detailed comparison, the parameters ρ, a, b, R, and U
0 would need to be specified.UsingMatlab ode23s, the shooting method outlined in Section 5 gave physically acceptable results for shear thickening fluids, Newtonian fluids, and for shear thinning fluids with n values down to 0.3. For 0
<
n
<
0.3 we did not obtain physically acceptable results with (
∂
w
/
∂
r
)
r
,
t
>
0 for 0
<
r
≤
1. To obtain feasible results we resorted to Matlab bvp4c for n
<
0.3. We are of the opinion that there may be a type change in the solution at some small value of n. 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 forB
0, 63, was used for all values of n. In Figure 3 the final value obtained for B
0 using the shooting method is plotted against n for a range of values of b. The graphs for b
=
2, b
=
5, and b
=
10 show that B
0 starts to increase rapidly for n below a critical value. This critical value of n increases slowly with b. It can be expected that the graphs for b
=
0.2 and b
=
1 would also show this property if the solution had been taken to smaller values of n. This rapid increase in B
0 demonstrates the difficulty in obtaining a numerical solution for very small values of n.Figure 3
The final value obtained forB
0 using the shooting method plotted against n for different values of b. The initial estimate for B
0 is B
0
=
2
1
+
b.In Figure2, which illustrates the matching principle, and Figures 4 to 7 which investigate the physical properties of the solution, the value of b used to perform the numerical calculations is b
=
2. In Table 1 the initial estimate for B
0 given by 63, the final value obtained for B
0 using the shooting method and the relative error percentage for a range of values of n are listed for b
=
2. We see that 63 underestimates B
0 for very small values of n, but otherwise it overestimates B
0. The relative error percentage in the initial estimate for B
0 increases rapidly as n decreases below n
=
0.3 and increases steadily with n for n
>
1. The shooting method with the initial estimate 63 still converges when the relative error in the initial estimate for B
0 exceeds 20
%.Table 1
Shooting method withb
=
2: n values, initial B
0
=
2
(
1
+
b
)
1
/
2 values, final B
0 values, and relative error percentage.
n
B
0 initial
B
0 final
Relative error%
0.25
3.4641
4.4237
21.69
0.3
3.4641
3.8554
10.15
0.3333
3.4641
3.2734
5.83
0.40
3.4641
3.2183
7.64
0.5
3.4641
3.1758
9.08
0.75
3.4641
3.1258
10.82
0.8
3.4641
3.1248
10.86
0.9
3.4641
3.1351
10.49
1
3.4641
3.1598
9.63
1.05
3.4641
3.1412
10.28
1.1
3.4641
3.1202
11.02
1.2
3.4641
3.0811
12.43
1.5
3.4641
2.9796
16.26
1.55
3.4641
2.9648
16.84
1.6
3.4641
2.9504
17.41
1.75
3.4641
2.9101
19.04
2
3.4641
2.8514
21.49
2.25
3.4641
2.8010
23.67
2.5
3.4641
2.7564
25.67
2.75
3.4641
2.7191
27.40
3
3.4641
2.6864
28.95
5
3.4641
2.6863
28.95Figure 4
Fluid velocityw
(
r
,
t
) derived numerically and plotted against r for a range of t values with b
=
2 and n
=
0.3333,0.5,1 and 3.
(a)
(b)
(c)
(d)The final value ofB
0 has important physical significance. From 45 and 46 the fluid velocity on the axis of the cylinder, r
=
0, is
(81)
w
0
,
t
=
B
0
t
,
and from 67 we see that for large values of time the wall velocity is given by
(82)
V
t
=
B
0
t
.
For large values of time the velocity profile becomes almost independent of r and the wall velocity tends to the value of the fluid velocity on the axis.Consider next the fluid velocityw
r
,
t for 0
≤
r
≤
1. In Figure 4  w
r
,
t is plotted against r for shear thinning fluids with n
=
0.3333 and 0.5, for a Newtonian fluid with n
=
1, and for a shear thickening fluid with n
=
3. For all values of n considered, the initial V-shaped velocity profile evolves into a smooth U-shaped profile for t
>
0 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 Figures5(a), 5(b), and 5(c) the scaled velocity w
r
,
t
/
w
1
,
t is plotted against r for n
=
0.3333,1, and 3. 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 t
=
0.10 that the fluid velocity on the axis r
=
0 is less for the shear thickening fluid with n
=
3 than for the Newtonian and shear thinning fluids with n
=
1 and n
=
0.3333, respectively. In Figure 5(d) the velocity of the cylindrical wall, V
t, is plotted against t for n
=
0.3333,1, and 3. Because of the no-slip boundary condition at the wall, V
t
=
w
1
,
t. The wall velocity is determined by the similarity solution and cannot be prescribed. We see that V
t increases approximately linearly for small values of t in accordance with expansion 65. For large values of time, V
t behaves approximately like B
0
t. From Table 1 for n
=
0.3333, n
=
1, and n
=
3 and for large time,
(83)
V
(
t
)
=
3.2734
t
,
V
(
t
)
=
3.1598
t
,
V
(
t
)
=
2.6864
t
,
respectively, in agreement with Figure 5(d). For large time the wall velocity decreases as n increases from values for a shear thinning fluid to a shear thickening fluid, except for a small range of n
(
0.5
<
n
<
1
) which may be due to a small computational error.Figure 5
In (a), (b), and (c) the numerical solution for the scaled fluid velocityw
(
r
,
t
)
/
w
(
1
,
t
) is plotted against r for b
=
2, n
=
0.3333, 1, and 3 and for a range of t values. In (d) the wall velocity V
(
t
)
=
w
(
1
,
t
) is plotted against t for b
=
2, n
=
0.3333, 1, and 3.
(a)
(b)
(c)
(d)Consider next the magnitude of the vorticity,ω
r
,
t, given by 70. In Figures 6(a), 6(b), and 6(c),  ω
r
,
t is plotted against r for a range of values of time for n
=
0.5,1, and 3. 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 t
→
∞ when the fluid velocity becomes uniform across the tube. In Figure 6(d) the magnitude of the vorticity at the wall, ω
1
,
t, is plotted against t for n
=
0.5,1, and 3. At t
=
0, ω
1
,
t
=
1 and it decreases approximately linearly with time for small values of time in agreement with expansion 73. It vanishes as t
→
∞. 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.Figure 6
In (a), (b), and (c) the magnitude of the vorticityω
(
r
,
t
) is plotted against r for t
=
0.0,0.01,0.10,1
,
2,3 and ∞ with b
=
2, n
=
0.5, 1, and 3. In (d) the vorticity at the wall ω
(
1
,
t
) is plotted against t for b
=
2 and n
=
0.5,1, and 3.
(a)
(b)
(c)
(d)Figure 7
In (a), (b), and (c) the scaled shear stressS
r
z
(
r
,
t
)
/
S
r
z
(
r
,
0
) is plotted against r for t
=
0.0,0.01,0.10,1
,
2,3 and ∞ with b
=
2, n
=
0.5,1, and 3. In (d) the scaled wall shear stress S
r
z
(
1
,
t
)
/
S
r
z
(
1,0
) is plotted against t for b
=
2 and n
=
0.5,1, and 3.
(a)
(b)
(c)
(d)Finally, consider the shear stressS
r
z
r
,
t given by 76. In Figures 7(a), 7(b), and 7(c), the shear stress, scaled by its initial value, is plotted against r for a range of values of t and for n
=
0.5,1, and 3. The initial profile of the shear stress is given by 77 and has the form of a right angle. For t
>
0 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 S
r
z
1
,
t
/
S
r
z
1,0 is plotted against t for n
=
0.5, 1, and 3. The decrease of the wall shear stress with time is approximately linear for small values of time in agreement with expansion 79. From 78,
(84)
d
d
t
S
r
z
1
,
t
=
-
1
+
b
1
+
n
b
+
O
t
,
as t
→
0. 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 b 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 V-shaped velocity profile into a flat profile. From Figure 3, the velocity on the axis of the tube, B
0
t, increases as the material parameter b 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 V- and U-shaped velocity profiles.


## 7. Conclusions

A preliminary investigation of the properties of the numerical solution was undertaken. Since we tookρ, a, and b to be the same for all the fluids considered, the investigation mainly demonstrated the effect of the power law exponent, n, 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 V-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 Section5 was applicable for a large range of values of n. It gave physically acceptable results for shear thickening and Newtonian fluids and for shear thinning fluids for values of n as low as n
=
0.25. This is a strong indication that the numerical solution is correct. If the solution for values of n close to each other is being considered then instead of using the initial estimate 63 for B
0 the final value for B
0 from the previous solution could be used.TheMatlab 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 3 and 2 with error control for stiff systems. It advances from time step k to time step k
+
1 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 I
I
I
a formula. This is a collocation formula and the collocation polynomial provides a C
1-continuous solution that is fourth-order accurate uniformly in interval [
a
,
b
]. 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 ast
→
∞. 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.



---

*Source: 135424-2015-09-17.xml*