PERISTALTIC MOTION OF A JOHNSON-SEGALMAN FLUID IN A PLANAR CHANNEL

This paper is devoted to the study of the two-dimensional flow of a Johnson-Segalman fluid in a planar channel having walls that are transversely displaced by an infinite, harmonic travelling wave of large wavelength. Both analytical and numerical solutions are presented. The analysis for the analytical solution is carried out for small Weissenberg numbers. (A Weissenberg number is the ratio of the relaxation time of the fluid to a characteristic time associated with the flow.) Analytical solutions have been obtained for the stream function from which the relations of the velocity and the longitudinal pressure gradient have been derived. The expression of the pressure rise over a wavelength has also been determined. Numerical computations are performed and compared to the perturbation analysis. Several limiting situations with their implications can be examined from the presented analysis.


Introduction
The dynamics of the fluid transport by peristaltic motion of the confining walls has received a careful study in the recent literature. The need for peristaltic pumping may arise in circumstances where it is desirable to avoid using any internal moving parts such as pistons in a pumping process. Moreover, the peristalsis is also well known to the physiologists to be one of the main mechanisms of fluid transport in a biological system. Specifically, peristaltic mechanism is involved in swallowing food through the oesophagus, urine transport from the kidney to the bladder through the ureter, movement of chyme in the gastrointestinal tract, and transport of spermatozoa in the ductus efferentes of the male reproductive tracts. Moreover, in the cervical canal, it is involved in the movement of ovum in the Fallopian tube, transport of lump in the lymphatic vessels, and vasomotion of small blood vessels such as arterioles, venules, and capillaries, as well as in the mechanical and neurological aspects of the peristaltic reflex. In plant physiology, such a mechanism is involved in phloem translocation by driving a sucrose solution along tubules by peristaltic contractions. In addition, peristaltic pumping occurs in many practical applications involving biomechanical systems such as roller and finger pumps. The of the solution describes the flow characteristics in the case of a Newtonian fluid. The perturbed parts of the solution are the contributions from the Johnson-Segalman fluid. Section 7 is concerned with the numerical solution of the nonlinear partial differential equation. In Section 8, numerical results and discussions are given. The analytical results are also compared to a numerical computation in order to determine the range of validity of the perturbation analysis.

Basic equations
The basic equations governing the flow of an incompressible fluid are the field equations div V = 0, div σ + ρf = ρ dV dt , (2.1) where V is the velocity, f the body force per unit mass, ρ the density, d/dt the material time derivative, and σ is the Cauchy stress. Johnson and Segalman [9] proposed an integral model which can also be written in the rate-type form. With an appropriate choice of kernel function and the time constants, the Cauchy stress σ in such a Johnson-Segalman fluid is related to the fluid motion through σ = −pI + T, (2.2)

3)
S + m dS dt + S(W − aD) + (W − aD) T S = 2ηD, (2.4) where D is the symmetric part of the velocity gradient and W the skew-symmetric part of the velocity gradient, that is, Also, −pI denotes the indeterminate part of the stress due to the constraint of incompressibility, µ and η are viscosities, m is the relaxation time, and a is called the slip parameter. When a = 1, the Johnson-Segalman model reduces to the Oldroyd-B model; when a = 1 and µ = 0, the Johnson-Segalman model reduces to the Maxwell fluid; and when m = 0, the model reduces to the classical Navier-Stokes fluid. Note that the bracketed term on the left-hand side of (2.4) is an objective time derivative.

Formulation of the problem and flow equations
Consider a two-dimensional infinite channel of uniform width 2n filled with an incompressible Johnson-Segalman fluid. We choose a rectangular coordinate system for the channel withX along the center line andȲ normal to it. LetŪ andV be the longitudinal and transverse velocity components of the fluid, respectively. We assume that an infinite train of sinusoidal waves progresses with velocity c along the walls in theX-direction. The 4 Peristaltic motion of a Johnson-Segalman fluid geometry of the wall surface is defined as where b is the amplitude and λ is the wavelength. We also assume that there is no motion of the wall in the longitudinal direction (extensible or elastic wall). For unsteady two-dimensional flows, In the fixed coordinate system (X,Ȳ ), the motion is unsteady because of the moving boundary. However, if observed in a coordinate system (x,ȳ) moving with the speed c, it can be treated as steady because the boundary shape appears to be stationary. The transformation between the two frames is given bȳ The velocities in the fixed and moving frames are related bȳ where (ū,v) are components of the velocity in the moving coordinate system.

Rate of volume flow and boundary conditions
The dimensional rate of fluid flow in the fixed frame is given by whereh is a function ofX and t. The rate of fluid flow in the moving frame is given by whereh is a function ofx alone. With the help of (3.9) and (3.10), one can show that these two rates are related through The time-averaged flow over a period T at a fixed positionX is given bȳ If we define the dimensionless time averaged flows Θ and F, respectively, in the fixed and moving frame as we find that (4.5) reduces to where If we choose the zero value of the streamline along the center line (y = 0) then the shape of the wave is given by the streamline of value The boundary conditions for the dimensionless stream function in the moving frame are Ψ = 0 (by convention), on the center line y = 0, and ∂Ψ ∂y = −1 (no slip condition), at the wall y = h. We also note that h represents the dimensionless form of the surface of the peristaltic wall: is the amplitude ratio or the occlusion and 0 < Φ < 1.
8 Peristaltic motion of a Johnson-Segalman fluid

Equations for large wavelength
A general solution of the dynamic equations (3.20) for arbitrary values of all parameters seems to be impossible to find. Even in the case of Newtonian fluids, all analytical solutions obtained so far by Shapiro et al. [25], and by L. M. Srivastava and V. P. Srivastava [29] are based on assumptions that one or some of the parameters are zero or small. Accordingly, we carry out our investigation on the basis that the dimensionless wave number in (3.18) is small, that is, which corresponds to the long-wavelength approximation [25]. Thus, to lowest order in δ, equations (3.20) give Substituting (5.4) and (5.6) into (5.5) yields and from (5.2), (5.3), and (5.7), we finally obtain

Perturbation solution
For small values of ᐃe 2 , (5.8) and (5.9) can be written using the binomial theorem as T. Hayat et al. 9 where the dimensionless parameters α 1 and α 2 are defined as Now, we seek the solution of (6.1) with boundary conditions (4.11) and (4.12) for a small Weissenberg number. We may expand flow quantities in a power series of ᐃe 2 . We write the stream function Ψ, the pressure field p, and the flow rate F in the following forms: If we substitute (6.3) into (4.11), (4.12), (5.3), and (6.1), and separate the terms of different orders in ᐃe 2 , we obtain the following systems of partial differential equations for the stream function and pressure gradients together with boundary conditions. 6.1. System of order ᐃe 0 . The following system of equations of zeroth order follows: with the boundary conditions This boundary value problem is that of linear viscous creeping flow in the long-wavelength approximation. No properties of the Johnson-Segalman fluid enter its formulation.
6.2. System of order ᐃe 2 . The first-order differential equations are 10 Peristaltic motion of a Johnson-Segalman fluid with the boundary conditions At this level, the material properties-manifested via a and η-now affect the flow inside the fluid layer.
6.3. System of order ᐃe 4 . The system of equations of the second order is composed of with the boundary conditions (6.9) In this system, further corrections due to the Johnson-Segalman constitutive equation enter. We now seek to solve the sequence of problems at each order and generate thereby the series solution.

Zeroth-order solution.
The solution to the zeroth-order problem (6.4) subject to the boundary conditions (6.5) is given by From the second and third equations in (6.4), it is clear that the transverse pressure gradient is zero and the longitudinal pressure gradient is given by The pressure rise per wavelength ( P λ ) in the longitudinal direction can be evaluated on the axis at y = 0. Thus, at the zeroth order, we have where The expressions (6.10), (6.11), and (6.12) are essentially the same as those of the Newtonian fluid given by Shapiro et al. [25]. This is obviously no surprise.

First-order solution (ᏻ(ᐃe 2 )
). Substituting the zeroth-order solution into (6.6), the system of order ᐃe 2 reduces the latter to (6.14) On solving (6.14) with boundary conditions (6.7), the expressions for the stream function Ψ 1 , the axial velocity component u 1 , the longitudinal pressure gradient dp 1 /dx, and the pressure rise per wavelength P λ1 turn out to be where I 4 = π 3Φ 2 + 2 (6.16) 6.6. Second-order solution (ᏻ(ᐃe 4 )). If we insert the zeroth-order and first-order solutions into (6.8), we find that the system of ᏻ(ᐃe 4 ) takes the form where I n is given by (6.16).

T. Hayat et al. 13
Now we summarize our results of the perturbation series through order ᐃe 4 . The expressions for Ψ, u, dp/dx, and P λ may, respectively, take the following forms: On substituting these expressions into (6.25) and retaining only terms up to order ᐃe 4 , we obtain (6.28)

Numerical method
We will solve the differential equation ( Because the differential equation (7.1) is nonlinear, we cannot solve this boundary value problem by the direct finite-difference method. In fact, in solving such nonlinear equations, iterative methods are usually used, and one of them is the asymptotic method.
A general stationary problem f Ψ(y), y = 0, y ∈ (a,b), accompanied with suitable boundary conditions, in which f may be a nonlinear higherorder differential function, can be treated as the steady asymptotic limit as t → ∞ of the nonstationary problem ∂φ ∂t + f φ(y,t), y = 0. In solving the stationary problem (7.3) by asymptotic methods, that is, in solving the nonstationary problem (7.4), we do not pay any attention to the transient behavior since it is of no interest whatsoever. Clearly, the nonstationary problem for φ, (7.4), can be solved by using finite-difference methods with respect to t. For example, where the counting index n indicates the time step and has nothing in common with the channel width n. As for the stationary problem (7.3), the solution is given by lim n→∞ φ n = Ψ. (7.8) In any case, from the point of view of solving the stationary problem, it is convenient to interpret the index n as denoting the iteration step rather than time, and in this case the real parameter τ indicates the iteration over-relaxation factor rather than the size of the time step. The parameter τ should be, on the one hand, sufficiently small to guarantee convergent iteration, and, on the other hand, as large as possible to minimize the number of iterations. If f is a differential function, we can obtain the discrete approximate solutions φ 1 (y 1 ), φ 2 (y 2 ),...,φ M (y M ) at discrete points a < y 1 < y 2 < ··· < y M < b by employing the difference equation φ n+1 j = φ n j − τg φ n 1 ,...,φ n M , y j , (7.9) in which the algebraic function g can be obtained by means of finite-difference approximations with respect to the algebra-differential function f in (7.7).
Here, we describe this method as applied to the boundary value problem (7.1) and (7.2). Instead of solving the stationary problem (7.1), we solve the nonstationary equation with the expression for c n j given by for the iteration steps n = 0,1,2,...,∞ and at space points j = 2,3,...,M − 1 (where y = 2∆y,3∆y,...,(M − 1)∆y) with the grid size ∆y = h/(M − 1). The boundary conditions (7.2) can be described in finite-difference forms as follows: We start with an initial trial solution φ 0 j ( j = 0,1,2,...,m + 1), which satisfies the boundary conditions (7.13). The iteration should be carried out until the relative differences of the computed φ n j and φ n+1 j between two successive iterative steps for all discrete points are smaller than a given error chosen to be 10 −8 . Therefore, the stationary solution, the solution of the stationary boundary value problem (7.1) and (7.2), is achieved.

Numerical results and discussions
A comparison of the direct numerical results as obtained for the boundary value problem (7.1) and (7.2) and its approximate perturbation solutions to order ᐃe 4 , (6.22) and (6.23), is shown in Figure 8.1. We display the dimensionless stream function Ψ (left panels) and the dimensionless velocity u obtained with three different values of the Weissenberg number ᐃe = 0.8,1.0,1.5 and a fixed total flux F = −2. It is surprising and pleasing that, for fairly large ᐃe, for example, ᐃe = 0.8, good agreement between the approximate solutions and the direct numerical solutions can be achieved. When ᐃe = 1, and more so when ᐃe > 1.0, their difference becomes conspicuously large. In fact, such a performance can be predicted. It is dependent on the choice of the total flux F. For F = −2 used in Figure 8.1, the distribution of the stream function is close to linear, as we can see in Figures 8.1(a), 8.1(c), and 8.1(e). In such a case, the term with ᐃe in (7.1), which is proportional to the second-order spatial derivative of the stream function, cannot become important. We can demonstrate this even more clearly if we choose F = −1. In this case, the solution of the boundary value problem yields a rigorously linear distribution of the stream function (i.e., constant velocity), and hence its second-order spatial derivative vanishes. Therefore, the results of the approximate perturbation solutions are in accordance with those of the direct numerical solutions which hold for any value of ᐃe. In fact, if we choose F = −4, for which a larger deviation of the stream function from the linear distribution is obtained, an obvious difference of the approximate perturbation solutions from the direct numerical solutions can be seen, although ᐃe takes a fairly small value ᐃe = 0.4, as displayed in Figure 8.2. Only when ᐃe = 0.1 the perturbation solution is an adequate approximation.
To quantitatively discriminate to what extent the perturbation solutions (6.22) and (6.23) can describe the boundary value problem (7.1) and (7.2), an error measure for a physical variable ϕ, for example, the stream function Ψ and the velocity u, is introduced: where ϕ num j denotes the direct numerical solution at the space position y j , while ϕ app j is the corresponding approximate value obtained by the perturbation solutions (6.22) and (6.23).
For various Weissenberg numbers ᐃe, the errors of the two solutions are listed in Table 8.1 with two different values of the total flux F. It is obvious that the errors increase with the increasing Weissenberg number ᐃe. For F = −2, if ᐃe > 0.8, the approximate solutions can no longer be used, while for a larger total flux value F = −4, the approximate perturbation solutions have to be abandoned for as small values as ᐃe > 0.2. In such cases, this boundary value problem must be solved by direct numerical methods.
In Figure 8.3 the distributions of the stream function Ψ(y) and the velocity u(y) are illustrated for various values of the total flux F. Results are obtained by the direct numerical method because here we choose a large value of ᐃe = 1.5 and hence the perturbation  solutions are no valid approximation. Obviously, the flow is dominantly influenced by the value of F. For a half-channel width of h = 1, if F > 1, the flow velocity near the channel center is larger than that at the top boundary due to an exerted pressure gradient in the direction opposite to the flow, while for F < 1, a pressure gradient in the flow direction is required to maintain such a flux value, hence the flow velocity near the channel center is smaller. For F = 1, a linear distribution of the stream function and a constant velocity are formed; in this case, no pressure gradient exists.
The distributions of the pressure gradient dp/dx within a wavelength x ∈ [0,2π] are exhibited in Figure 8.4 for various values of the dimensionless wave amplitude Φ ( Figure  8.4(a)) and the total flux F (Figure 8.4(b)). The dimensionless half-channel width is defined by (4.13), that is, h(x) = 1 + Φ sinx. The results are plotted for the approximate perturbation expression (6.24) because we here choose a small value of ᐃe = 0.2, and hence an adequate approximation by the perturbation method can still be ensured.  (Figure (a)) and the total flux F (Figure (b)). The other parameters are µ/η = 1, a = 0.8, and ᐃe = 0.2. larger ᐃe, one has to directly solve the boundary value problem (7.1) and (7.2) for various values of the channel width h within a wavelength and then substitute the results into (5.9) in order to obtain the pressure gradient at various points along the channel.
It can be clearly seen from Figure 8.4 that, on the one hand, in the wide part of the channel, x ∈ [0,π], the pressure gradient is relatively small, that is, the flow can easily pass without imposition of a large pressure gradient. On the other hand, in a narrow part of the channel, x ∈ [π,2π], a much larger pressure gradient is required to maintain the same flux to pass it, especially for the narrowest position near x = 3π/2 and when the flux F or the wave amplitude Φ is larger.
The corresponding pressure drops in the flow direction over a wavelength are listed in Table 8.2 for some values of F and Φ. These are still obtained by the perturbation method, that is, (6.25) due to the used small value of ᐃe = 0.2. For small values of F, a pressure rise occurs in the flow direction. With the increase of the flux F and the wave amplitude Φ, the pressure drops increase over a wavelength.

Conclusions
Plane steady peristaltic motions of a Johnson-Segalman fluid were analyzed for the conditions that the wall surface is sinusoidally deformed. This wall deformation has wavelength λ which is large in comparison to the undeformed channel width and propagates at a prescribed constant speed. The time-independent governing equations of this boundary value problem are expressed in terms of the stream function and pressure as basic unknowns with wall boundary conditions formulated in terms of the stream function alone. By prescribing the constant value of this stream function of the wall, the flux caused by the peristaltic wall deformation is prescribed. A scale analysis and associated nondimensionalization of the equations paired with a Galilee-transformation with speed c discloses a steady time-independent problem in which an aspect ratio parameter δ and the Weissenberg number ᐃe appear as the significant physical scales. The equations in the limit δ → 0 describe the large-wavelength approximation. These equations still contain the Weissenberg number ᐃe as a parameter. Perturbation solutions up to ᏻ(ᐃe 4 ) for the flow have been constructed for prescribed flux F and these solutions have been compared with numerical solutions that are valid for any value of the Weissenberg number. The following results are found.
(i) The deviation of the flow from a linear stream function and a constant velocity depends on both the flux F as well as the Weissenberg number ᐃe. For F = −1, the viscoelastic properties do not enter the flow no matter how large ᐃe is.
(ii) With growing |F| (decreasing F to large negative values) the influence of the viscoelastic properties of the fluid becomes more and more significant.
(iii) The profiles of the stream functions and the longitudinal velocities as computed with the perturbation expansion are only adequately predicted when the Weissenberg number ᐃe is small. The maximum values for which the perturbation solutions are valid approximations depend on the flux variable F. It is the smaller, the larger the deviation F from −1 is.
(iv) In general, the solutions ought to be numerically determined by using the full equations valid for all ᐃe.
(v) The pressure that develops for a certain flux depends very much on the flux parameter. Generally, rather localized maxima can arise (Figure 8.4). Such maxima can be controlled by adequately choosing the flux.
The next analysis will obviously be the corresponding peristaltic motion of a fluid in a circular flexible duct. This problem is presently under analysis.