Bifurcation Analysis for Peristaltic Transport of a Power-Law Fluid

1 Introduction

Peristaltic transport is involved in many industrial and biological systems such as flow of lymph through lymphatic vessels, transport of urine in urethra, movement of blood in arterioles, venules and capillaries and flow of spermatozoa in male reproductive tract. This mechanism refers to the flow of the fluid in a tube or a channel by the movement of walls. It was observed in the living system for many centuries. Latham [1] was the first who discussed the peristaltic transport in context to mechanical pumping. Shapiro et al. [2] performed mathematical modelling and basic analytical study of peristaltic flow for Newtonian fluids in a wave frame of reference. Fung and Yih [3] adopted an alternative approach and discussed this mechanism in a laboratory frame of reference.

In most of the studies pertaining to peristaltic flows, the fluid rheological behaviour is characterised by Newtonian model. However, many industrial and biological fluids are non-Newtonian in nature, as they do not conform to the Newtonian behaviour due to their complex rheology. Typical effects associated with non-Newtonian fluids are shear-thinning, shear-thickening, Weissenberg effects, yield stress and memory effects. The simplest model used in the literature to capture shear-thinning and shear-thickening effects is the power-law model. Many researchers [4], [5], [6], [7], [8] have discussed the peristaltic flow problems for non-Newtonian fluids by considering power-law model under different geometric setups with low Reynolds number and long wavelength approximations. Recently, Wang et al. [9] highlighted the slip effects on the rotating electro-osmotic flow in a slowly varying microchannel by incorporating power-law model. The main objectives of the above-mentioned studies were to investigate flow and pumping features due to peristalsis and to characterise the conditions for which trapping and reflux occur. Trapping refers to the condition when a circulating bolus of fluid trapped inside closed streamline is transported with wave speed for certain values of the involved parameters. Siddiqui and Schwarz [10] investigated the peristaltic transport of non-Newtonian fluid in axisymmetric conduit and observed that for large flow rates, the trapping could break on the longitudinal axis. Two eddies are sustained below the wave crest and some fluid flow through the middle line. The same behaviour was observed by Pozrikidis [11] while investigating the Stokes flow through a channel. The above-mentioned studies reported the appearance and breakage of trapping for certain values of the involved parameters but did not present a rigorous bifurcation analysis using the methods of dynamical systems.

A bifurcation in peristaltic flow occurs when a change in the parameter causes a sudden topological change in its flow behaviour. The objective of this study is to reveal the nature and types of these bifurcations. In a previous study, Hartnack [12] discussed streamline patterns along with their bifurcation near a fixed wall for a viscous two-dimensional incompressible flow. He studied the change of streamline topologies near degenerate critical points. Brønse and Hartnack [13] investigated the streamline topologies near simple degenerate points of two-dimensional viscous flow away from the boundary. Jiménez and Sen [14] analysed the streamline patterns and their bifurcations for Newtonian fluid driven by peristaltic waves. Asghar and Ali [15], [16] extended the work of Jiménez and Sen by considering slip and heat convection effects. However, no such study is yet available for non-Newtonian fluids and thus our aim is to explore the streamline patterns and their bifurcations for peristaltic flow of shear-thinning and shear-thickening fluids characterised by the power-law model. The flow is considered in a plane symmetric channel. Local and global bifurcations are also analysed. Theoretical results are verified by the experimental data available in the literature.

2 Basic Equations

The basic equations for the flow of an incompressible fluid in the absence of body force are

where ρ is the density of fluid, V¯ is the velocity field of fluid, t is the time, d/dt is the material time derivative and T¯ is the Cauchy stress tensor.

For Ostwald-de Waele power-law model, the Cauchy stress tensor is defined as [17]

where p is the pressure, S¯ is the extra stress tensor, μ₀ is the consistency index, n is the fluid behaviour index, γ is the symmetric part of velocity gradient tensor, i.e.

The quantity μ0[12(𝜸¯:𝜸¯)]n−1 in (4) is called apparent viscosity. Equations (3) and (4) represent shear-thinning fluids for n < 1 and shear-thickening fluids for n > 1. When n = 1, (3) and (4) reduce to the Newtonian model.

3 Development of the Flow Field

We assume peristaltic flow of an incompressible power-law fluid in a two-dimensional symmetric channel with flexible walls. The geometry of the channel walls in Cartesian coordinate system (x^∗, y^∗) is shown in Figure 1 and given by the equation

Figure 1:

Geometry of the periodic wave.

where R₀, ω, d and c stand for the mean channel half-width, wavelength, wave amplitude and speed of the wave, respectively. Consider the motion in a wave frame of reference, (x¯,y¯), which is moving with velocity c relative to the fixed frame. The motion becomes time-independent in the wave frame of reference.

The coordinate transformation between fixed frame of reference, (x∗,y∗), and the wave frame of reference (x¯,y¯), is given by

where (u¯,v¯) and (u∗,v∗) are the velocity components in wave and fixed frames of reference, respectively. Equation (5), in wave frame of reference, becomes

For a channel flow, (2) in the wave frame of reference can be written as

Defining the dimensionless quantities as

Equations (8)–(11) and (7) become

where

In the above equations φ=d/R0 and Re=ρc2−nϵR0n/μ0 are amplitude ratio and Reynolds number, respectively, for power-law fluids. By assuming the long wavelength (ϵ ≈ 0) and low Reynolds number approximations (Re=O(1)), (13), (14) and (16) reduce to

The dimensionless no-slip boundary conditions in wave frame of reference are:

As the axial velocity u obtains the maximum value at the centreline, thus,

Further

In order to study the trapping phenomena, the stream-function ψ is introduced via the relations u=∂ψ/∂y and v=−∂ψ/∂x. Elimination of pressure term from (17) gives

Further, the conditions (18)–(21) become

Figure 2:

Streamlines for different flow rates q in wave frame of reference with n = 0.5 and φ = 0.8, (a) Backward flow, q = −4/5, (b) Trapping, q = −3/10, (c) Augmented flow, q = 1/5.

Figure 3:

Streamlines for different flow rates in wave frame with n = 1.5 and φ = 0.8, (a) Backward flow, q = −4/5, (b) Trapping, q = −3/10, (c) Augmented flow, q = 1/5.

The volume flow rate in wave frame of reference is given by

The symmetry at the centreline of the channel implies

First, the solution for the region above the centreline (y ≥ 0) shall be obtained. Integrating (22) twice with respect to y, using (26) to remove modulus sign and then applying the boundary condition (24), one finds

where A is the constant of integration. Now taking the nth root of (30) and ensuring that it must be taken only for positive quantities, we get

Integrating (31) and then utilising the boundary conditions (23), (28) and (29) yields the solution for the region above the centreline as

The above solution is not valid for y < 0 because it gives complex values for y < 0. Moreover, it does not satisfy the boundary condition at the lower wall. Repeating the procedure described above, the solution for the region y < 0 is given by

The non-uniqueness of the solution for peristaltic flow of power-law fluid in a channel and the derivation of corresponding non-unique solution is already discussed thoroughly by Subba Reddy et al. [18].

Taking n = 1, the solutions (32) and (33) reduce to the unique solution for Newtonian fluid, i.e. [14]

By analysing the solutions reported above one can identify three different types of flow situations: backward flow, trapping and augmented flow [2], [10]. If the entire fluid flows in a direction which is opposite to the wave direction, it is known as backward flow. Trapping exists when the streamlines split to contain bolus of fluid particles. Augmented flow refers to the condition when there is some fluid flowing forward and trapped bolus splits. These flows are shown in Figures 2 and 3 for shear-thinning and shear-thickening fluids, respectively.

4 Flow Field as a Non-linear Dynamical System

In this section, the qualitative theory of dynamical systems is used to study the change in the flow behaviour. At a particular instant (say t = t₀), the movement of individual particles moving in paths defined by x˙=v(x,t) is similar to instantaneous streamlines, i.e. x˙=v(x,t0). The present problem can be written as a system of non-linear autonomous differential equations by using the definition x˙=∂ψ∂y and y˙=−∂ψ∂x as:

For y ≥ 0,

For y ≤ 0,

where 𝜶=[φ,q] be the parameters, −∞<x<∞ and −h<y<h is our domain of interest and range of amplitude ratio is 0<φ<1. To obtain the critical points x~=(x~,y~), take

Now employ the Hartman–Grobman theorem [18] to investigate the qualitative nature of these critical points by evaluating the Jacobian J=∂(f±,g±)/∂(x,y) at these points (x~,y~). When the determinant of Jacobian becomes zero at (x~,y~), the critical point is classified as the degenerate point, which can be further classified into two subclasses [19]: simple degeneracy and non-simple degeneracy. The condition when eigenvalues of Jacobian matrix are zero refers to the simple degeneracy and if Jacobian matrix becomes null matrix, it corresponds to non-simple degeneracy. The Bakker’s notation of two-dimensional system [20] is used to classify the critical points. Trace: P~ab=λa+λb and Jacobian: d~ab=λaλb are used for the classification, where λ_a and λ_b are eigenvalues.

According to Seydel [21], solutions (x~,y~,𝜶c) of (39) and (40) are bifurcation points with respect to parameter 𝜶 where the number of quasi-periodic, periodic or equilibria solution changes when 𝜶 passes through 𝜶_c; 𝜶_c be the critical value.

Critical points satisfying (39) and (40) are:

1. {x~1,2,y~1,2}={mπ,±[(2n+1)q−2n(2n+1)(q+2)]nn+1},m∈ℤ

2. {x~3,4,y~3,4}={±arccos(2nφ−2n−(2n+1)q2nφ), 0},

3. {x~5,6,y~5,6}={(2m−1)π2,±(1−φ)[(2n+1)q+2n(1−φ)(2n+1)(q+2(1−φ))]nn+1},m∈ℤ

By taking n = 1, the above expressions reduce to the corresponding expressions for Newtonian fluid [14].

In the next section, the classification of critical points and their graphical representation of local and global bifurcations will be investigated. The critical values of the involved parameters will be also discussed.

4.1 Classification of the Critical Points {x~1,2,y~1,2} and their Bifurcation

The critical points, {x~1,2,y~1,2}={mπ,±[(2n+1)q−2n(2n+1)(q+2)]nn+1},m∈ℤ, are situated on a line which is vertically below the wave crest. At these points, Jacobian matrix is:

J|{x~1,2,y~1,2}=[0S1,2T1,20],

(41)whereS1,2=∓(2n+1)(q+2)2n((2n+1)q+2n(2n+1)(q+2))1n+1

(42)T1,2=±φ(2(n+1)+(2n+1)q)n+1((2n+1)q+2n(2n+1)(q+2))2n+1n+1∓(2n+1)qφn+1((2n+1)q+2n(2n+1)(q+2))1n+1,

and eigenvalues are

(43)λ1,2=±−2φ(2n+1)(2n+(2n+1)q)(q+2)(2n+1)(q+2),

(44)withP~12=0 and d~12=2φ(2n+(2n+1)q)(2n+1)(q+2).

It is noted that by changing the value of q, the stability and nature of critical points also change. The flow rate q is assumed to lie in the interval (−1, 1)

1. The critical point corresponds to a one-dimensional saddle in the range −1<q<−2n2n+1 as P~12=0 and d~12<0 in this range; see Figure 4(i).

2. An isolated critical point exists when q=qc1=−2n2n+1. It is known as non-hyperbolic degenerate point [22]. As, both the Jacobian matrix, J, and eigen values, λ1,2, are zero at q=qc1=−2n2n+1, therefore, the critical point at this value corresponds to a non-simple degeneracy [19]; see Figure 4(ii).

3. d~12>0 for q>−2n2n+1, so, each critical point in this range corresponds to a centre: see Figure 4(iii).

Figure 4:

Bifurcation for x~=mπ,m∈ℤ, and picture of topological changes: (i) q<−2n2n+1, (ii) q=−2n2n+1, and (iii) q>−2n2n+1.

A bifurcation diagram in the q−y~ plane for different fluid behaviour indices is presented in Figure 4. Bifurcation occurs on the vertical line at x~=mπ,m∈ℤ. This bifurcation is co-dimensional two as it depends upon flow rate q and fluid behaviour index n.

4.3 Classification of the Critical Points {x~5,6,y~5,6} and their Bifurcation

In the previous section, for −2n2n+1<q<−2n2n+1(1−φ), a fixed degenerate point bifurcates to saddle nodes on the longitudinal axis. When q approaches to qc2=−2n2n+1(1−φ), saddle nodes of the adjacent wave coincide below the wave troughs. Critical points combine at x~=(2m+1)π2 when q=qc2. This joint produces degenerate point having six heteroclinic connections. When q>qc2, the degenerate point bifurcates to saddle nodes on vertical line under wave trough.

The critical points, {x~5,6,y~5,6}={(2m−1)π2,±(1−φ)[(2n+1)q+2n(1−φ)(2n+1)(q+2(1−φ))]nn+1},m∈ℤ, lie on a line which is vertically below the wave trough. At these points, the Jacobian matrix is:

J|{x~5,6,y~5,6}=[0S5,6T5,60],

where

(49)S5,6=∓(2n+1)(q+2(1−φ))2n(φ−1)2((2n+1)q+2n(1−φ)(2n+1)(q+2(1−φ)))1n+1,

(50)T5,6=±4n(φ−1)(2n+1)(q+2(1−φ))((2n+1)q+2n(1−φ)(2n+1)(q+2(1−φ)))nn+1,

with eigenvalues

(51)λ5,6=±12n+12(2n+1)φ(1−φ)(q+2(1−φ))((2n+1)q+2n(1−φ))(4+q)φ−2φ2−2−q

(52)and P~56=0,d~56=−2φ((2n+1)q+2n(1−φ))(2n+1)(1−φ)(q+2(1−φ)).               

For q≥−2n2n+1(1−φ), the changes in nature of critical points at x~=(2m+1)π2 occur as follows:

1. As P~56=0 and d~56<0for q>−2n2n+1(1−φ), therefore, the critical points correspond to saddle nodes; see Figure 6(iii).

2. The critical point for q=qc2=−2n2n+1(1−φ) is degenerate because P~56=0=d~56. Also, the eigenvalues and Jacobian matrix at this critical point are zero, therefore, it corresponds to a non-simple degeneracy; see Figure 6(ii).

Figure 5:

Bifurcation diagram for y~=0 and topological changes: (i) q=−2n/(2n+1), (ii) q>−2n/(2n+1), and (iii) q=−2n(1−φ)/(2n+1).

Bifurcation diagrams for different values of fluid behaviour index n in the q−y~ plane at x~=(2m+1)π/2 are shown in Figure 6a–d.

5 Global Bifurcation and Streamline Patterns

For y = 0, the vector field becomes {x˙,y˙}={2n+12(n+1)(q+2h)(1h)−1,0}, from which

Critical points occur under wave crests at x~=mπ and wave troughs at x~=(2m−1)π/2; m∈ℤ. Bifurcation curves are obtained by setting

The global bifurcation diagram for different fluid behaviour indices n in the space φ−q has the following curves:

Figure 6:

Bifurcation for x~=(2m−1)π/2,m∈ℤ and topological changes: (i) q<−2n(1−φ)/(2n+1), (ii) q=−2n(1−φ)/(2n+1), and (iii) q=−2n(1−φ)/(2n+1).

Along the bifurcation curve M, there exist non-simple degenerate critical points under the wave crests. Whereas, along the bifurcation curve N, adjacent critical points coincide under the wave trough and heteroclinic connections containing non-simple degenerate points are formed. Bifurcation curves for different fluid behaviour indices n are shown in Figure 9. The flow region is divided as:

1. Backward flow: when all the fluid is flowing in a direction, which is opposite to the wave direction.

2. Trapping: when critical points, characterised as saddles, are connected by heteroclinic connections and there exist an interaction of vortices in the flow with opposite rotations.

3. Augmented flow: when the eddies join under the wave crests and produce heteroclinic connections along with their neighbours and some fluid flows through the centreline in the wave motion direction.

6 Results

Different types of streamline topologies and their bifurcation for shear-thinning and shear-thickening fluids are shown in Figures 4–8. Figure 4 explores the bifurcation, at q=qc1, corresponding to the critical points {x~1,2,y~1,2}. As fluid behaviour index n increases, the bifurcation point moves backward in the flow rate direction. That is, the bifurcation for shear-thickening fluid appears at a low flow rate as compared to shear-thinning fluid. This translation of bifurcation point results in narrowing down of the eddying region under wave crest, that is, the distance between two centres under same wave crest decreases with increasing the fluid behaviour index n. Figure 5 shows that, when flow rate further increases to qc1, saddle nodes on the longitudinal axis move towards neighbouring wave trough. This means the eddying region expands along vertical direction and spreads along the longitudinal direction with increasing q. For a small amplitude ratio φ, the eddying region spreads quickly. As q approaches to qc2, these saddle nodes merge under the wave trough to form a non-simple degenerate point with six heteroclinic connections as represented in Figure 6. Figure 6d represents the bifurcation at q=qc2 corresponding to the points {x~5,6,y~5,6} for different values of n.

Figure 7:

Global bifurcation (i)–(v) correspond to n = 0.5 and φ = 0.6, for different values of q: (i) q = −4/5, backward flow, (ii) q = −1/2, (iii) q = −1/3, trapping, (iv) q = −1/5 and (v) q = 1/5, augmented flow.

Figures 7 and 8 present the transition of streamline topologies for shear-thinning and shear-thickening fluids, respectively, by taking different flow rate values q and fixed amplitude ratio φ. In each figure, the streamline patterns for degenerate cases are shown in panels (ii) and (iv). Two possible bifurcations occur as volume flow rate, q, and amplitude ratio, φ, are varied. Panels (i)–(iii) of both Figures 7 and 8 represent the formation of symmetric eddying zone on the axis. These eddying zones are enclosed by heteroclinic connections between saddle nodes. Panels (iii)–(v) of both figures describe the merging of adjacent eddying zone, where the critical points, situated on the longitudinal axis, coincide and lift off under wave trough to produce heteroclinic connections joining saddles. The global bifurcation diagram is shown in Figure 9. It is noted that the trapping region increases by increasing the fluid behaviour index n. That is, the range of volume flow rate q for which trapping occur is the greatest for shear-thickening fluids and least for shear-thinning fluids. Further, it is observed that the streamline patterns shown in Figures 7 and 8 are in accordance with the predictions of Figure 9.

Figure 8:

Global bifurcation diagram (i)–(v) correspond to n = 1.5 and φ = 0.6, for different values of q: (i) q = −4/5, backward flow, (ii) q = −3/4, (iii) q = −1/2, trapping, (iv) q = −3/10 and (v) q = 1/5, augmented flow.

7 Comparison with the Experimental Results

In this section, the obtained theoretical results are verified by comparing them with the experimental data of Weinberg et al. [23]. Weinberg et al. performed experimental analysis of two-dimensional peristaltic transport induced by sinusoidal waves with Reynolds number ranging from the inertia-free limit to value in which the inertial effects are significant. In order to visualise the reflux and trapping phenomena, dyed fluid was used and observations were carried out in the moving frame by employing a camera attached with the experimental setup. They measured mean flow, mean-pressure rise, pressure pulses and trapping phenomena for φ = 0.7 and 0.9. Based on these measurements, they calculated the experimental values of Q∗/Q0 for which trapping occurs, where Q^∗ and Q₀ be the time-mean flow rate and dimensionless time-mean flow rate at zero pressure rise, respectively. For establishing the theoretical trapping limits of Q∗/Q0, Weinberg et al. [23] used the dimensionless wall form given by

Figure 9:

Global bifurcation diagram for different fluid behavior indices n, (i) backward flow region, (ii) trapping region, (iii) augmented flow region.

Using the analysis reported in this article, the critical values for trapping at the moving frame of reference are given by

The non-dimensional time-mean flow rate is

Replacing (59) and (60) in (61) gives

These trapping limits agree with those given by Subba Reddy et al. [18] for symmetric channel.

The dimensionless time-mean flow rate at zero pressure rises Q₀, can be obtained by integrating the pressure gradient ∂p/∂x, over x per wavelength. Under long wavelength and low Reynolds number approximations, the pressure gradient is given by

The pressure rise per wavelength is

In order to calculate the dimensionless time-mean flow rate at zero pressure rises, i.e. Q0=Q∗(△pω=0), the integral in (65) is evaluated numerically by using MATHEMATICA. For trapping to occur, we must have

The above inequality gives the range in which trapping appears. Qc1∗/Q0 and Qc2∗/Q0 are the critical conditions corresponding to the branches of global bifurcation. These relations are plotted in (Q∗/Q0)−φ plane and compared them with the experimental results of Weinberg et al. for viscous fluid as shown in Figure 10. Clearly, there is an excellent agreement between the experimental and theoretical predictions for n = 1. It is important to mention that Weinberg et al. predicted experimental results only for viscous Newtonian fluids (n = 1). No study is yet available where such measurements are reported for shear-thinning and/or shear-thickening fluids. Moreover, at present, we are not involved (our work is theoretical and computational) in the experimental validations due to lack of sophisticated laboratories required for such experiments. In this scenario, we are not in a position to comment on the validity of our theoretical results. However, the authors would welcome experimentalists to compliment the present theoretical predictions with the laboratory investigations of this phenomenon.

Figure 10:

Comparison of observed trapping limits and experimental results of Weinberg et al. [23] for n = 1: • trapping; ∘ no trapping. Region I: backward flow, region II: trapping and region III: augmented flow.

8 Conclusion

The peristaltic flow of shear-thinning and shear-thickening fluids was analysed by considering the power-law model in order to discuss the streamline patterns and their bifurcations. The flow is taken in a two-dimensional channel. The system of non-linear autonomous differential equations was developed for application of the theory of dynamical system. Eigenvalues of the Jacobian matrix were used to classify the critical points and it was found that the critical points were either center or saddle. The bifurcation analysis of critical points suggested three different types of flow situations: backward flow, trapping and augmented flow. The conversions of backward flow to trapping and then trapping to augmented flow correspond to bifurcations where non-simple degenerate points change their stability to produce heteroclinic connections between saddle nodes. It was observed that increase in fluid behaviour index caused the appearance of trapping at low flow rates. The trapping region expands and backward region shrinks by increasing the fluid behaviour index. The trapping limits for different n were well defined in global bifurcation diagram. By this approach, one can easily explore the ranges of parameter in which each kind of flow occur and characterise the stagnation points by their position and stability. This analysis is necessary to disclose the complete information about the flow, i.e. the fluid transports directly from entrance to exit or some fluid traps within the wall and also where it traps. The theoretical results of this article are compared with the experimental results available in the literature.