Steady flow in a dividing pipe

J. Fluid Mech. (1999), vol. 401, pp. 339–364. c 1999 Cambridge University Press Printed in the United Kingdom 339 Steady flow in a dividing pipe By M. G. B L Y T H AND A. J. M E S T E L Mathematics Department, Imperial College, 180 Queen’s Gate, London SW7 2BZ, UK (Received 11 November 1998 and in revised form 17 May 1999) The high Reynolds number flow through a circular pipe divided along a diameter by a semi-infinite splitter plate is considered. Matched asymptotic expansions are used to analyse the developing flow, which is decomposed into four regions: a boundary layer of Blasius type growing along the plate, an inviscid core, a viscous layer close to the curved wall and a nonlinear corner region. The core solution is found numerically, initially in the long-distance down-pipe limit and thereafter the full problem is solved using down-pipe Fourier transforms. The accuracy in the corners of the semicircular cross-section is improved by subtracting out the singularity in the velocity perturbation. The linear viscous wall layer is solved analytically in terms of a displacement function determined from the core. A plausible structure for the corner region and equations governing the motion therein are presented although no solution is attempted. The presence of the plate has little effect ahead of the bifurcation, but wall shear on the curved wall is found to increase from its undisturbed value downstream. 1. Introduction We reconsider Smith’s (1977a) problem of the flow of a homogeneous incompressible Newtonian fluid along an infinite straight pipe of circular cross-section which is divided at some point by a semi-infinite flat plate across its diameter (see figure 1). Far upstream of the division we assume the motion to be steady and fully-developed. The main aim of this paper is to examine the entry flow problem when this uni-directional profile impinges on the divider. A long way downstream the flow will return to a Poiseuille motion, this time in two semicircular pipes. It is the intermediate development that is of concern here. The large number of industrial and biological applications has stimulated much study of pipe flows. In this paper we will concern ourselves chiefly with the problem’s physiological relevance. A review of the fluid mechanical aspects of the arterial system is given in Pedley (1980). In general, haemodynamic flows are rendered difficult to model by multifarious complicating factors. The problem of a bifurcating artery seems to have been originally addressed analytically by Zamir & Roach (1973). The effects of unsteadiness at pipe bifurcations have been considered numerically, for example in the computer simulations of Friedman, O’Brien & Ehrlich (1975) and Wille (1984), and experimentally, for example by Rieu & Pelissier (1991) who measured unsteady flow in a dividing pipe for different angles of bifurcation. The assumption of steady flow is perhaps the most significant simplification made in the present work, limiting its physiological applicability. However it may act as a first step in the model of a bifurcating artery on which future research could be based. Any reservations about its applicability notwithstanding, it is an interesting problem worthy of study in its own right. Nevertheless we adopt a biological standpoint 340 M. G. Blyth and A. J. Mestel y x z U0 Figure 1. Circular pipe with splitter-plate for x > 0. and in this spirit results are presented (see § § 3, 4) which would be of physiological interest but which we concede may be far removed from their likely counterparts in vivo. Caro, Fitz-Gerald & Schroter (1971) observed that the location on the artery walls of fatty streaks, which are believed to lead to the onset of atherosclerosis, may be influenced by certain haemodynamic factors since they are particularly found in regions associated with low wall shear. Preferential sites include the outer walls of dividing arteries and the inside walls of bends. Here we discuss a simple model of a bifurcating artery where the angle of bifurcation is taken to be zero. In reality, of course, branching arteries in the body will diverge with daughter tubes of generally different cross-sectional areas. Many nearly symmetric divisions occur in the airways in the lungs but the only close approximation to a symmetric branching in the systemic circulation occurs where the aorta splits into the iliac arteries (Pedley 1980). For simplicity we assume that the branching is symmetric, although the analysis is easily adaptable to allow for an asymmetric splitter plate. This is the most straightforward way of allowing for different pressures in the daughter tubes a long way downstream. A more difficult extension would be to admit different pressures with the geometrical symmetry maintained, in which case the uni-directionality of the leading order core flow would be lost. The assumption of a fully-developed upstream flow is also not necessarily valid. Given the many consecutive bifurcations in the arterial system, it may never be allowed the necessary room to evolve. Our steady motion occurs in a pipe which is taken to be perfectly straight and to suffer no bending or torsion. Curvature effects, and those of unsteadiness, have been considered analytically, notably by Lyne (1970) and Smith (1975), who studied pulsatile flow in a curved pipe driven by a sinusoidal pressure gradient. In reality more complicated three-dimensional geometrical effects may also be important and recently, for example, the effect of torsion has been taken into account by Zabielski & Mestel (1998a,b), who examined both steady and unsteady flow in a helically symmetric pipe. Zabielski & Mestel also included a realistic physiological pressure pulse. Gammack & Hydon (1999) considered pipes with curvature and torsion varying slowly with the downstream direction. Arteries are elastic in nature and are to some extent free to move. Some of the effects of the changing curvature of the coronary arteries during systole and diastole (the two-part cycle of a heart beat) have been given a preliminary treatment by Lynch, Waters & Pedley (1996) who investigated unsteady flow in pipes with time- Steady flow in a dividing pipe 341 dependent curvature. Distensibility is important in the propagation of the pressure pulse generated by the heart. Here we consider a rigid pipe of infinite length and to a first approximation believe that distensibility will be of little importance especially given the many other simplifications adopted. In this respect we follow the example of other similar studies of pipe flows that have been motivated at least in part by physiological considerations. This neglect may not be justified for joining veins, where blood pressure is in general much lower than in arteries, and which are subsequently more prone to collapse. Clearly in this eventuality elastic properties of the walls may not be ignored. The blood rheology is also to be considered. Blood effectively consists of a suspension of red blood cells in a plasma solution. On the lengthscale of a typical arterial diameter the assumption of homogeneity is justified. Experimental measurements (for example Whitmore 1968) suggest that blood behaves as a Newtonian fluid although this might be complicated once unsteadiness is incorporated into the model. The high Reynolds number limit adopted in the current work may also be supported by experimental observations. Pedley (1980) reports a peak value of 1000 for the Reynolds number in a canine femoral artery. In order to examine the entry flow in a bifurcation, the developing viscous layers and the three-dimensional dynamical effects of the branching, Smith (1977a) investigated the steady flow along an infinitely long straight circular pipe with a semi-infinite splitter plate sliced through the middle. He studied the flow with the down-pipe distance x scaled with Reynolds number which renders the equations of motion parabolic and of boundary layer type. Smith argued that the fluid flow may be divided up into eight distinct regions, with a Blasius-type boundary layer growing along the plate, a core region above it, and then an inertial and a viscous region above this on the pipe wall. The remaining regions effect the turning of the secondary flow in the corners. However, doubt was cast on the scalings Smith used at the upper edge of the core by his student Bennett (1987) in his PhD thesis. In fact, Bennett argues (as we do here) that there is a choice of two possible scalings at the top of the core, the alternative set requiring the core pressure to vanish on the curved wall (this would correspond to Smith’s P1 = 0). In this case, which we now believe to be correct, there is no requirement for an intermediate inertial layer and in fact the whole structure is simplified to just four main zones as detailed below. Our problem is thus essentially identical to Smith’s although we do not scale with Reynolds number in the down-pipe direction; instead we solve the flow in the regime x = O(1) and also consider the behaviour of the core flow when 1 ≪ x ≪ R. As usual with a boundary layer on a flat plate of infinitesimal thickness, around the leading edge a tiny region of width O(R −1 ), in this case broadening to O(R −2/3 ) in the corners, is needed in which the flow is described by the full Navier–Stokes equations. This zone is unimportant as far as the rest of the flow structure is concerned and a more detailed discussion is omitted. The paper is arranged as follows. In § 2 we describe the asymptotic structure of the flow and lay out the equations governing the separate regions of motion. The nonlinear equations in the corner region are given but no attempt is made to solve them. An analytic solution may be found for region III which is dependent on a function prescribed by the core. The core solution is obtained numerically and is outlined in § 3, where we consider both x = O(1) and the large x limit. Results and a discussion of the effect of the splitter plate on the undisturbed flow are presented in § 4. 342 M. G. Blyth and A. J. Mestel y III II r IV I h z Figure 2. The four distinct flow zones. 2. Flow structure We assume the steady flow to be Newtonian, incompressible and of constant viscosity. We therefore use the three-dimensional Navier–Stokes equations u · ∇u = −∇p + R −1 ∇2 u, ∇ · u = 0, (2.1a) (2.1b) where u = (u, v, w) in the cylindrical polar coordinates (x, θ, r) as in figure 2. Lengths have been non-dimensionalized with the pipe radius a∗ . For example x∗ = a∗ x, where an asterisk here denotes a dimensional quantity. Also u∗ = U ∗ u, v ∗ = U ∗ v, w ∗ = U ∗ w, p∗ = ρ∗ U ∗2 p. U ∗ is a characteristic velocity of the unperturbed flow and ρ∗ is the density of the fluid, assumed constant. Far upstream the flow is uni-directional and is given by u = U0 (r), v = w = 0, p = p0 − xR −1 , for constant p0 . The Reynolds number, R, (assumed R ≫ 1) is based on the pipe radius and is defined by R = U ∗ a∗ /ν ∗ , with ν ∗ the kinematic viscosity. We believe that the flow structure can be divided into four main regions (illustrated in figure 2). A short distance beyond the leading edge a boundary layer I starts to grow along the plate. Near the centre of the tube this will be of Blasius type, based on the local mainstream velocity. The boundary layer perturbs the core flow II above by a finite amount. The displacement thickness of the plate layer is of order R −1/2 and we find that the perturbation to the core is then of order R −1/2 , driven by the displacement effect of the plate layer. The equations governing the core flow are inviscid to leading order and consequently a slip velocity is left at the edge of II on r = 1. Thus another region, the viscous wall layer III, is needed to reduce the velocity to zero at the curved wall. This layer is linked with the core through the function A(θ) described below. A singularity arises in III as the corner of the semicircle is approached and this is smoothed out by a further viscous region IV which otherwise plays a passive role in the flow structure. 2.1. The core (II ) In the inviscid core the flow is perturbed as follows: u = U0 (r) + R −1/2 u1 (x, r, θ) + · · · , v = R −1/2 v1 (x, r, θ) + · · · , w = R −1/2 w1 (x, r, θ) + · · · , p = R −1/2 p1 (x, r, θ) + R −5/6 p2 (x, r, θ) + · · · , (2.2a) (2.2b) (2.2c) (2.2d) Steady flow in a dividing pipe 343 where r2 = y 2 + z 2 . The coordinate frame x, y, z is shown in figure 1, with θ = 0 on y = 0, z > 0. The origin of the additional pressure perturbation p2 is discussed in § 2.3. Far upstream the flow is given by U0 = 41 (1 − r2 ). (2.3) Sufficiently far downstream the core expansions (2.2) will become invalid. From continuity considerations (discussed in detail in § 2.1.1) we see that u1 ∼ x1/2 as x → ∞ and so (2.2a) will break down when x ∼ R. This corresponds to the boundary layer on the plate (I) growing large enough to fill the pipe. Therefore there exists a down-pipe lengthscale x = RX, where X is of order unity, for which X ≪ 1 will match to (2.2) as x → ∞, and for which the flow will develop into Poiseuille motion in a semicircular tube when X ≫ 1. It is the former lengthscale where 1 ≪ x ≪ R which was considered by Smith (1977a). Substituting (2.2) into (2.1a) and taking the limit R → ∞, we derive the following inviscid equations at leading order (1 − r2 )u1x − 2rw1 (1 − r2 )v1x (1 − r2 )w1x u1x + v1θ /r + (rw1 )r /r = −4p1x , = −4p1θ /r, = −4p1r , = 0. (2.4a) (2.4b) (2.4c) (2.4d) Following Smith (1976a) we can easily eliminate u1 , v1 , and w1 , and write (2.4) as a single equation for p1 4r p1r = 0, (2.5) 1 − r2 where ∇22 = ∂2 /∂r2 + (1/r)∂/∂r + (1/r2 )∂2 /∂θ2 , the two-dimensional Laplacian for the cross-section. This equation is elliptic and therefore the possibility of upstream influence is anticipated. The core flow is driven by the Blasius layer and so this provides the boundary condition for (2.5) on the base (θ = 0, π). We also need to consider what happens on the curved wall to determine the boundary condition at r = 1. As r → 1 we can approximate (2.5) by p1xx + ∇22 p1 + p1nn − 2p1n = 0, n (2.6) where n = 1 − r. Then, as n → 0, two different scalings are possible for p1 . Either p1 ∼ 1 or p1 ∼ n3 . Initially for the case x → ∞ Smith (1977a) proposed the former and solved the core problem with two boundary conditions applied on the base. However his student Bennett (1987), in his PhD thesis, suggested that in fact the latter choice is correct and it is these scalings which we adopt here, as argued below. It then seems reasonable that in the absence of some effect which is dominant for x ∼ 1 but which is asymptotically small when x is large, p1 = 0 on r = 1 for finite x also. In light of this we assume u1 ∝ 1, v1 ∝ n2 , w1 ∝ n, p1 ∝ n3 , (2.7) as n → 0. This then determines the second boundary condition required to solve (2.5), p1 = 0 on r = 1. (2.8) 344 M. G. Blyth and A. J. Mestel 2.1.1. The limit x → ∞ We use a continuity argument to show that the following system represents the correct limiting problem in region II for large x. We will see in § 2.2 that injection into the core is proportional to x−1/2 . If we integrate along the upper half of the pipe from a point upstream of the leading edge of the plate to some point downstream of it, we argue that as x → ∞ u1 will behave like x1/2 a long way downstream. From (2.4), either p1 ∼ x1/2 or p1 ∼ x−3/2 . To determine which is correct, we consider the first, writing u1 = u0 (r, θ)x1/2 , v1 = v0 (r, θ)x−1/2 , w1 = w0 (r, θ)x−1/2 , p1 = p0 (r, θ)x1/2 . Then as x → ∞, the core equations (2.4) become − r2 )u0 − 2rw0 = −2p0 , p0θ = p0r = 0, 1 u + v + (1/r)(rw 0θ 0 )r = 0. 2 0 1 (1 2 (2.9a) (2.9b) (2.9c) If we integrate over 0 < θ < π, 1 u 2 + 1 (1 − r2 )u − 2rw 2 [v0 ]π0 + (1/r)(rw)r = −2p0 π, = 0, (2.10a) (2.10b) where an overline indicates an averaged √ variable here and p0 , by (2.9b), is a constant. From (2.15), [v0 ]π0 = −β(1 − r2 )1/2 / 2. To achieve a balance in (2.10a) as r → 1 (or n → 0, writing n = 1 − r) we need either u ∼ 1/n or w ∼ 1. Eliminating w from (2.10) and integrating we find √ u = −2 2β(1 − r2 )1/2 + const. (2.11) and so it must be the case that w ∼ 1 as n → 0 unless p0 = 0. However this conclusion leads to an imbalance in mass transport inside the viscous wall layer. If w ∼ 1 as n → 0 then the fluid flux into III is of order unity. The boundary layer thickness is of order R −1/3 . By (2.11) there will be a mass flux downstream within the boundary layer of order R −1/3 at most which cannot cater for the larger influx. Therefore we conclude that neither w ∼ 1 nor u ∼ 1/n as n → 0 and the only option remaining to us is to set p0 = 0. This corresponds to the constant P1 in Smith (1977a) being zero. In this case Smith’s complicated asymptotic structure is unnecessary and we take the boundary condition for (2.5) on r = 1 to be that of vanishing pressure. 2.2. Base boundary layer (I ) The precise nature of the base condition is determined by matching to the plate layer I. As remarked above, the thickness of this layer is O(R −1/2 ) and so we set y = R −1/2 Y , where Y is of order unity. A Blasius-type similarity solution is possible for x > 0. We write u = (1 − z 2 )f ′ (η)/4 + O(R −1/2 ), v=R w= −1/2 −1/2 2 1/2 −1/2 (1 − z ) x √ (ηf − f)/(2 2) + · · · , R w1B (x, Y , z) + · · · , R −1/2 pB1 (x, Y , z) + · · · , ′ (2.12a) (2.12b) (2.12c) (2.12d) √ where η = U0 (z)1/2 Y x−1/2 / 2. Dhanak & Duck (1997) demonstrate that this solution is non-unique when the outer flow is non-zero on both walls, and that an O(1) crossflow is possible in this layer. Here, however U0 (1) = 0. At leading order z merely behaves as a parameter and the layer is governed according to (2.12) by the familiar p= 345 Steady flow in a dividing pipe Blasius equation f ′′′ + f f ′′ = 0, (2.13) with f(0) = 0, f ′ (0) = 0, f ′ → 1 as η → ∞. (2.14) Then, since as η → ∞, f ∼ η − β + exp, where ‘exp’ denotes exponentially small terms and β ≈ 1.21678 is the Blasius displacement, 1 1/2 v ∼ R −1/2 √ βU0 x−1/2 , 2 (2.15) as η → ∞, x > 0. By (2.4b), the boundary condition on the base of the core becomes d 1 p1θ = − √ βr(1 − r2 )3/2 (H(x)x−1/2 ) dx 8 2 on θ = 0, π, (2.16) valid ∀x. Here H(x) is the Heaviside function. As x → ∞, this condition may be written as 1 (2.17) p1θ = √ βr(1 − r2 )3/2 x−3/2 on θ = 0, π. 16 2 2.3. The viscous wall layer (III ) This region is necessarily viscous in order to bring the fluid at the top of the core to rest and so satisfy the no-slip condition on the curved wall. From (2.1a), the neglect of viscous terms in the core becomes invalid when inertial terms balance viscous ones. As r → 1, U0 ∼ 21 (1−r) and so viscous effects become important when (1−r) ∼ R −1/3 . By (2.2a), the core equations would break down only when (1 − r) ∼ R −1/2 , that is beyond where viscous effects must be taken into account. Region III is therefore taken to be of order thickness R −1/3 and we write (1 − r) = R −1/3 ζ, with ζ a new wall variable of order unity. Note that it is no longer necessary to include an extra layer here as in Smith (1977a). We then perturb the mainstream velocities in this viscous wall layer as follows u = 12 R −1/3 ζ + R −1/2 u1 (x, ζ, θ) + · · · , v = R −1/2 v 1 (x, ζ, θ) + · · · , w = R −5/6 w 1 (x, ζ, θ) + · · · , p = R −5/6 p1 (x, ζ, θ) + · · · , (2.18a) (2.18b) (2.18c) (2.18d) according to the scalings (2.7). The relevant equations of motion are then linear: 1 ζ 2 u1x − 21 w 1 = −p1x + u1ζζ , 1 ζ 2 v 1x = −p1θ + v 1ζζ , p1 = p1 (x, θ), u1x + v 1θ − w 1ζ = 0, (2.19a) (2.19b) (2.19c) (2.19d) with necessary matching conditions to the core u1 → 21 A(x, θ), and the no-slip condition v1 → 0 as ζ → ∞, u1 = v 1 = w 1 = 0 on ζ = 0. (2.20) (2.21) 346 M. G. Blyth and A. J. Mestel The function A(x, θ) comes from the solution of the core. In § 3 we find it using numerical methods. By (2.4) and (2.7), as r → 1 we find 1 A 2 xx = −p1rrr |r=1 . (2.22) Note that the azimuthal velocity v is larger than at the top of the core. Here it is of order R −1/2 , whereas in II, due to (2.7), v ∼ R −7/6 as r → 1, which is considerably smaller in magnitude. However (2.20) permits a match as we leave the wall layer. The pressure perturbation here, given by (2.18d ) and required to complete a consistent structure in the wall layer, is larger than that at the edge of the core and so provokes the extra perturbation p2 to II seen in (2.2d ). The equations and boundary conditions laid out above are all valid in the range −∞ < x < ∞ and it is convenient now to take a generalized Fourier transform in x. We write Z ∞ 1 ∗ u1 = √ u1 (x, ζ, θ)e−ikx dx, (2.23) 2π −∞ with the other variables similarly defined. Then, transforming (2.19) we obtain the solutions (cf. Smith 1976a, 1979; Cowley 1983) Z t Ai(q) dq − p∗1θθ M(t)/(2ω 5 ), u∗1 = 23 A∗ (k, θ) (2.24a) 0 v ∗1 = p∗1θ M(t)/ω 2 , Z ζ Z ωζ ′ ∗ ∗ 3 w 1 = 2 ikA Ai(q) dq dζ ′ , 0 (2.24b) (2.24c) 0 where Ai is the Airy function, and with Z q Z t 1 Ai(q ′ ) dq ′ dq. M(t) = Ai(t) 2 Ai (q) ∞ 0 (2.25) Here t = ωζ, ω = ( 21 ik)1/3 , and we have made the branch cut −3π/2 < arg(ω) 6 π/2 to ensure that Ai(t) → 0 as ζ → ∞. From (2.19), (2.20), and (2.21), p∗1 satisfies the equation p∗1θθ − k 2 p∗1 = 3ω 5 Ai′ (0)A∗ . (2.26) Assuming the flow has left/right symmetry we have the condition that p∗1θ = 0 at θ = π/2. (2.27) We solve (2.26) in § 4 but first we must determine the second boundary condition. This comes from a careful consideration of the flow structure in the corner. 2.4. The corner region (IV ) In order to detail the nature of the flow in the corner of the pipe, we need to consider what is happening in the corner of II. Part of our aim is to determine the second boundary condition needed to solve for the pressure in III. To this end (2.26) shows that we will need to consider how A∗ behaves as θ → 0. In the corner of II we are able to construct a local solution to (2.5) and we use this to examine the behaviour of A∗ there. As n → 0 and θ → 0 we may, after a Fourier transform, approximate (2.5) as ∗c p∗c 1nn + p1θθ = 2p∗c 1n , n (2.28) Steady flow in a dividing pipe 347 where the neglected −k 2 p∗c 1 term is assumed negligible in comparison with the others if we move sufficiently close to the corner. Here and throughout the rest of this paper the superscript or subscript c denotes a corner quantity. We can transform (2.28) to Laplace’s equation with the substitution p∗c 1n . (2.29) n Fixing local polar coordinates (ρ, γ) in the corner, with γ = 0 corresponding to n = 0, we have the following problem for the local solution: ψ= ∇2c ψ = 0, with (2.30) (2.31) ψ = 0 on γ = 0, ψγ = 3λc k 1/2 ρ1/2 on γ = π/2, √ 1 iπ/4 2 where λc = 16 βe 2, which comes from (2.16) (see also § 3.2.1), and ∇c = ∂2 /∂ρ2 + (1/ρ)∂/∂ρ + (1/ρ2 )∂2 /∂γ 2 . This has the solution √ (2.32) ψ = 6 2λc k 1/2 ρ1/2 sin(γ/2). Changing back to (n, θ) coordinates we can integrate once to give p∗c 1 = 12 λ k 1/2 (ξ 5 c − θ)1/2 (n2 − 31 θ(ξ − θ)), where ξ 2 = n2 + θ2 . Now, for θ ≪ n ≪ 1 p∗c 1 = 12 λ k 1/2 n5/2 5 c − 2λc k 1/2 n3/2 θ + O(θ2 ). (2.33) (2.34) This helps to determine the nature of the crossflow in the plate layer as |z| → 1. According to the scalings (2.12) the momentum equation in the transverse direction for the crossflow w1B in region I is B B B + vs w1Y = −pB1z + w1Y us w1x Y, (2.35) where us , vs denote the similarity solutions detailed in (2.12). As Y → ∞, the leading B ∼ −pB1z . Equation (2.34) implies the behaviour order balance suggests that us w1x B 3/2 B p1z ∼ (1 − z) as z → 1. Hence w1 ∼ (1 − z)1/2 as z → 1. A similarity solution for w1B can be found, valid as either |z| → 1 or x → ∞ (see Blyth 1999). Returning to the solution (2.33), we also find that for n ≪ θ ≪ 1 √ 2λc k 1/2 θ−1/2 n3 + O(n4 ), (2.36) p∗c 1 = in which case A∗ ∼ θ−1/2 as θ → 0 (see § 3). From (2.26) we must then have p∗1θθ ∼ θ−1/2 as θ → 0, requiring that p∗1 = α0 + α1 θ + α2 θ3/2 + O(θ2 ), (2.37) for small θ. As the corner must effect the turning of the fluid, pressure gradients should be expected in both the r- and θ-directions. Furthermore, if α1 6= 0 then by (2.24b) v ∗1 will be of order unity as the corner is approached from within the viscous wall layer. This is too large, as a consistent match with the base boundary layer requires the corner region to be at least as thick as the viscous wall layer. We therefore set α1 = 0, thereby providing the required second boundary condition for (2.26), p∗1θ = 0 on θ = 0. (2.38) We believe then that IV should form a square region of size O(R −1/3 ), with region I 348 M. G. Blyth and A. J. Mestel correspondingly thickening at each end. Accordingly we claim the following structure for the corner when x > 0 (cf. Walton & Smith 1997): u = R −1/3 uc1 (x, ζ, Θ) + R −1/2 uc2 (x, ζ, Θ) + · · · , (2.39a) w= (2.39c) R −2/3 v1c (x, ζ, Θ) + · · · , R −2/3 w1c (x, ζ, Θ) + · · · , R −5/6 pc1 (x) + R −4/3 pc2 (x, ζ, Θ) v= p= (2.39b) + ···, (2.39d) with (1 − r) = R −1/3 ζ, θ = R −1/3 Θ, and Θ ∼ 1. The nonlinear governing equations are uc1 uc1x + v1c uc1Θ − w1c uc1ζ c c c − w1c v1ζ + v1c v1Θ uc1 v1x c c c + v1c w1Θ − w1c w1ζ uc1 w1x c c − w1ζ uc1x + v1Θ = uc1ζζ + uc1ΘΘ , c c + v1ΘΘ , = −pc2Θ + v1ζζ c c c = p2ζ + w1ζζ + w1ΘΘ , = 0, (2.40a) (2.40b) (2.40c) (2.40d) with boundary conditions uc1 = v1c = w1c = 0 on ζ = 0, Θ = 0. Matching to region I we have, as ζ → ∞ keeping η = uc1 → 12 ζf ′ (η), With rc2 1 1/2 ζ Θx−1/2 2 (2.41) fixed, v1c → 12 ζ 1/2 x−1/2 (ηf ′ − f), (2.42a) (2.42b) w1c → x−1/2 ζ 1/2 h(η), pc2 → g(x)ζ 5/2 . = ζ + Θ we have, matching to region II as rc → ∞ for γ = tan−1 (ζ/Θ) 6= 0, 2 2 uc1 → 21 ζ + κ1 rc−1/2 , v1c → κ2 rc1/2 , Matching to region III we have, as Θ → ∞ uc1 → 12 ζ + κ5 Θ −1/2 , v1c → κ6 Θ 1/2 , w1c → κ3 rc1/2 , w1c → κ7 Θ −1/2 , pc2 → κ4 rc5/2 . pc2 → κ8 Θ 3/2 . (2.42c) (2.42d) The function h(η) comes from matching to the solution of (2.35) for w1B in the plate layer, which can be found with little difficulty; κi are functions of x and γ for i = 1, . . . , 4, and x and ζ for i = 5, . . . , 8, given in the Appendix along with g(x). The function pc1 is present in (2.39d ) because by (4.3) p∗1 is non-zero on θ = 0, π. It is purely a function of x since at next order in the x-momentum equation it is found that (2.43) pc1ζ = pc1Θ = 0. The next-order term in the expansion for u is required to be O(R −1/2 ) in order to balance −pc1x at second order in the x-momentum equation. In the corners, the tiny O(R −1 ) Navier–Stokes zone around the leading edge thickens to O(R −2/3 ) due to the smaller local downstream velocity there. Therefore, as in Rubin (1966), the solution in the corner would only be valid for 0 < x0 < x where x0 is unknown in advance but which is presumed to be small. Although it is nonlinear, and therefore interesting in itself, the corner region is essentially passive in nature with no controlling influence over the rest of the flow. It exists to ensure no-slip is satisfied in the corner and to effect the turning of the secondary motion. Justification of the structure proposed above could be obtained by a complete, presumably numerical, solution of the given equations which is not attempted here. Such a numerical solution might proceed along the lines of Rubin & Steady flow in a dividing pipe 349 Grossman (1971), who studied incompressible viscous flow along a right-angle corner. Guidance could also be obtained from the simpler, linear corner problem solved by Cowley (1983). 3. Numerical solution of the core region (II) Here we briefly describe the numerical procedure used to solve the core equation (2.5), and present some of our results. Two cases are considered. First we examine the solution of the core as x → ∞. Thereafter we consider the full x problem. The numerical procedure is similar in both cases. 3.1. Core solution for x → ∞ We showed in § 2.1.1 that p1 ∼ x−3/2 as x → ∞. We therefore write u1 v1 w1 p1 = x1/2e u1 + · · · , −1/2 e =x v1 + · · · , e1 + · · · , = x−1/2 w −3/2 e =x p1 + · · · , (3.1a) (3.1b) (3.1c) (3.1d) e in the large-x limit. Equations similar to (2.4) are then derived. A(θ), the core displacement function analogous to A(x, θ), is given by 1e A 2 = 4e p1rrr |r=1 . (3.2) Once more we can eliminate the other variables to obtain a single equation for e p1 (r, θ), ∇2 e p1 + with boundary conditions and 4r e p1r = 0, 1 − r2 (3.3) e p1 = 0 on r = 1, (3.4) e θ), e p1 = (1 − r2 )φ(r, (3.6) 1 √ βr(1 − r2 )3/2 on θ = 0, π, (3.5) 16 2 by (2.17). This is the same as Smith’s (1977a) equation (3.8a) with P1 = 0, although Smith applied two boundary conditions on the plate. We choose to solve it numerically using finite differences. It is expedient to make the substitution e p1θ = ± whereupon (3.3) becomes with boundary conditions e−4 e ≡ ∇2 φ Lφ eθ = ± φ 1 + r2 e φ = 0, (1 − r2 )2 1 √ βr(1 − r2 )1/2 16 2 on θ = 0, π, (3.7) (3.8) e = 0 on r = 1. In this way the equation is simplified somewhat in the and φ discretization. Also we anticipate that as r → 1, e p1 ∼ (1 − r)3 and consequently 2 e φ ∼ (1 − r) , and the finite difference approximations are more efficient at rendering this behaviour with the transformed equation. We use central differences and employ a second-order-accurate derivative approximation for the condition on the base, using 350 M. G. Blyth and A. J. Mestel ‘imaginary’ values one level below the base which are eliminated by application of the equation on θ = 0, π. Working in polar coordinates, we integrate over a small e there half-disc around the origin to obtain a condition for determining the value of φ (see for example Strikwerda 1989). The discretized equation was then solved using Gauss–Seidel point iteration with successive over-relaxation. Grids of various mesh lengths hr , hθ in the r-, θ-directions respectively were tried. Convergence of the grid was deemed to have been reached when the maximum difference over the grid between iterations was less than some prescribed tolerance multiplied by the maximum value over the grid at the present iteration. We used a tolerance of 10−7 . The scheme was found to be second-order accurate except in the corner of the semicircle. We know the local solution there exactly, given by (2.33), since the local analysis follows through in e a similar manner as x → ∞. From (2.36) we can see that close to the wall, as θ → 0, φ should behave like θ−1/2 . This square root singularity will clearly not be well resolved by finite difference approximations. We therefore attempt to subtract it out, writing where e−φ ec, e=φ ψ (3.9) e c = − 3β (ξ ′ − sin θ)1/2 (σ 2 − 1 sin θ(ξ ′ − sin θ)), (3.10) φ 3 40σ with 2σ = (1 − r2 ), and ξ ′ = σ 2 + sin2 θ. The derivative condition on the base is then e c has the correct behaviour in the corner and simplified (see (3.13b) below). Then φ the singularity should have been extracted to leading order. In addition, in order to e on r = 1, we need only approximate a second derivative on the wall since, estimate A by (3.6), e rr |r=1 . e (3.11) p1rrr |r=1 = −6φ The system we now have to solve is as follows: ec, Le ψ = −Lφ (3.12) e = 0 on r = 1, ψ e θ = 0 on θ = 0, π. ψ (3.13a) (3.13b) where L is the operator defined in (3.7). The boundary conditions are Finite differences are used as before and accuracy in the corner is somewhat improved e c behaves like although there remains a problem as now near the wall when θ → 0, Lφ 1/2 e θθ . However, we find nθ and so errors are still likely to arise in the discretization of ψ that taking hθ small enough results in good agreement in the corner with the known solution (2.33). Figure 3(a) makes a comparison of e p1rrr |r=1 used in the calculation of e according to (3.2). It shows plots of e A p1rrr |r=1 for the modified equation (3.7), the equation with the corner solution subtracted (3.12), and the corner solution (3.10) itself. We can see that the numerical results approach the local solution as θ → 0, the best approximation being given by the corner-subtracted solution. However, hθ has to be taken quite small to obtain this accuracy. Figure 3(b) shows how e p1 on θ = 0 compares with n5/2 as n → 0. We believe then that we have solved the core accurately over the whole semicircle except very far into the corners where, even having subtracted out the square root singularity in θ, finite differences are unable to reproduce the local solution for very small θ without taking an unfeasibly large number of grid points. 351 Steady flow in a dividing pipe 20 (a) 15 10 5 0 0.02 0.04 0.06 0.08 0.10 0.92 0.96 1.00 h 0 (b) – 0.001 – 0.002 – 0.003 0.80 0.84 0.88 r Figure 3. Numerical solution of II as x → ∞ in the corner. (a) n ≪ 1, local solution (——), solving (3.7) (- ⋄ - - ⋄ -), solving (3.12) (· · ·). (b) θ ≪ 1: e p1 (——) and −(3β/20)n5/2 (· · ·). e varies with θ. The singular θ−1/2 Figure 4 shows how the displacement function A behaviour is apparent near θ = 0, π. Figure 5(a) shows contours of pressure e p1 within the core. We can see that the pressure is greatest near the coordinate origin. In figure 5(b) the singular nature of the down-pipe velocity perturbation is reflected 352 M. G. Blyth and A. J. Mestel 3 2 1 2 Ã 1 0 1 2 3 h e varying with θ. Figure 4. 12 A(θ) (a) (b) Figure 5. Both figures refer to the limit x → ∞. Dashed lines indicate negative contours. (a) p1 |) = 0 on outer contour; Max(|e p1 |) = 0.019 at the origin. (b) Pressure contours of e p1 : Min(|e Contours of down-pipe velocity perturbation e u1 : Min(e u1 ) = 0.003 at the origin; Max(e u1 ) is at the corner singularities (see figure 4). in the way the contours bunch together in the corners. We note that despite the largest downstream velocities occurring near θ = 0, π, this does not correspond to the greatest fluid flux since the singularity disappears upon a single integration. Figure 6 displays contours of the secondary velocities in region II. Whilst the bulk of the fluid is being carried downstream, there is a secondary flux ‘upwards’ into the viscous wall layer III, in which region this additional fluid is flushed downstream. 3.2. The full x problem We now solve for the pressure in region II for all values of x. If we take the Fourier transform of (2.5) and make a similar substitution to that made in § 3.1, namely p∗1 = (1 − r2 )φ, we have ∇22 φ − 4 1 + r2 φ − k 2 φ = 0, (1 − r2 )2 (3.14) Steady flow in a dividing pipe (a) 353 (b) Figure 6. Secondary velocities in II as x → ∞ (dashed lines indicate negative contours). (a) Contours of v1 : Min(|v1 |) = 0 on r = 1, and |v1 | is largest around the origin. (b) Contours of w1 : Min(w1 ) = 0 on r = 1; w1 is largest around the origin. Note that the origin is purely a coordinate singularity. with φ = 0 on r = 1, and boundary condition on the symmetry line, φθ = ∓ 161 βeiπ/4 r(1 − r2 )1/2 k 1/2 (3.15) on θ = 0, π, (3.16) which comes from taking the transform of (2.16). For negative k the square root is taken according to the branch cut made in § 2.3. It is convenient to scale the k 1/2 from (3.16). This leaves us with the task of solving (3.14) according to φ = 0 on r = 1, (3.17) and (3.18) φθ = ∓ 161 βeiπ/4 r(1 − r2 )1/2 on θ = 0, π. Once again a numerical solution is sought using finite differences. The procedure is similar to the one employed previously, with convergence of the grid being reached using point iteration. There are similar problems to before in rendering the behaviour of the solution correctly in the corners. We now discuss the asymptotic behaviour of the solution for small and large k and show that, except far into the corners, our numerics approximate the known solutions well. 3.2.1. The case k ≪ 1 When k is small the local analysis has already been performed for p∗1 in § 2.4, the solution being given by (2.33). Therefore, since p∗1 ∼ 2nφ as r → 1, the local solution here is given by (3.19) φc = 56 λc (ξ − θ)1/2 (n2 − 13 θ(ξ − θ))/n, with λc = 161 βeiπ/4 21/2 arising from (3.18) as r → 1. A plot similar to that in figure 3(a) for x → ∞ reveals that if hθ is taken small enough, the numerical solution approaches the analytic solution for θ → 0. 3.2.2. The case k ≫ 1 From the form of (3.14) we expect φ to decay exponentially as |k| → ∞. Writing φ ∼ eb(r,θ)k for large positive k we derive the eikonal equation b2r + 1 2 b = 1, r2 θ (3.20) at O(k 2 ) and so, except near the corners, we anticipate that φ ∼ e−kr sin θ as k → ∞. (3.21) 354 M. G. Blyth and A. J. Mestel Near the curved wall we consider the approximation to (3.14) 2φ − k 2 φ = 0. (3.22) n2 Far enough into the corner the solution will be similar to (2.33). However when both n and θ are of order k −1 , the fourth term in (3.22) will become significant. In this case we write n = k −1 z and θ = k −1 y and expect the solution to be expanded as φnn + φθθ − φ = k −3/2 φ0 + · · · , (3.23) because of (3.18). Writing Ψ = φ0z + φ0 /z the problem reduces to 2 ∇ Ψ − Ψ = 0, with Ψ =0 Ψy = − 32 λc z −1/2 2 on z = 0, on y = 0, (3.24) (3.25a) (3.25b) and where ∇ = ∂2 /∂z 2 + ∂2 /∂y 2 . Taking a Fourier transform in z, and requiring that Ψ be an odd function of z, we identify the exact solution Z ∞ 2 1/2 sin(qz)e−y(q +1) 3λc dq. (3.26) Ψ=√ q 1/2 (q 2 + 1)1/2 2π 0 We are interested in comparing this with our numerical solution and so we use (3.26) to compute the local approximation to A∗ . Since we know that p∗1 ∼ n3 as n → 0, p∗1nnn |n=0 ≈ 4kΨz |z=0 , (3.27) where the k 1/2 factor has been scaled back into the right-hand side. The integral for Ψz |z=0 can be done analytically and thus, recalling (2.22), we find 3 = −k −1 27/4 Γ ( 34 )λc (kθ)−1/4 K1/4 (kθ), (3.28) π where K denotes the modified Bessel function of the second kind. From this we can see that as θ → ∞, A∗c ∼ θ−3/4 e−kθ , and when θ → 0, A∗c ∼ θ−1/2 as we would expect. Now A∗c is a purely local solution in the corner for large k and it is possible to construct a solution which is valid to leading order for large k over the entirety of region II. Writing Φ = φn + φ/n in (3.22) we approximate the solution to the resulting equation by solving (3.29) ∇22 Φ − k 2 Φ = 0 over the whole of the semicircle. Then as k → ∞ the most important terms around the walls of the semicircle will be encapsulated in (3.29). Its boundary conditions are Φ = 0 on r = 1, and on the base we take Φθ = g(z) where 1 ∗c A 2 3 g(z) = − √ λc (1 − z 2 )−1/2 sgn(z), 2 (3.30) which has been chosen so that Φθ approximates the base condition (3.25b) as z → 1. It is possible to solve this by the method of images and the use of Green’s theorem in the plane. Near r = 1 Z 1 1 (K0 (kr1 ) − K0 (kr2 ))g(z) dz, (3.31) Φ(r, θ) ≃ − π −1 355 Steady flow in a dividing pipe 0 (a) Large-k corner approxns –0.025 –0.050 (b) Numerical solution –0.075 –0.100 –0.125 0 0.05 0.10 0.15 0.20 h Figure 7. (a) Comparison in the corner of II at k = 10.0. Large-k corner approximations Re{ 21 A∗c }, Re{ 12 A∗G } (dotted lines) with the real part of the numerical solution (continuous line). (b) Contours of down-pipe velocity u1 in III at x = 5.0 (shown on a radially exaggerated scale). where r12 = z 2 + r2 − 2zr cos θ, r22 = z 2 + 1/r2 − 2z cos θ/r. (3.32a) (3.32b) Using Laplace’s method it is not hard to show that for large k, away from the corners, the exponential behaviour of Φ agrees with (3.21). Furthermore we demonstrate that (3.28) and (3.31) match up in the corner. As θ → 0, from (3.28) we have, 1 ∗c A 2 ∼ − 43 βeiπ/4 θ−1/2 k −3/2 . (3.33) ∗ Now for the Green’s function solution (3.31) we may approximate A as 1 ∗ A 2 G = 4k −2 Φr |r=1 (3.34) ∼ − 34 βeiπ/4 θ−1/2 k −3/2 , (3.35) and then as θ → 0 it is possible to show that 1 ∗ A 2 G and so the two solutions agree. Figure 7(a) shows how Re{ 21 A∗c } and Re{ 21 A∗G } compare with the numerical approximation of Re{ 12 A∗ } for small values of θ. The two analytic approximations coincide as predicted and both compare favourably with the numerics. The discrepancy between analytic and numerical solutions worsens as θ gets smaller. This is also to be expected since, as discussed above, the numerical solution becomes increasingly inaccurate in the corners and unlike the case x → ∞ we have made no attempt here to improve accuracy for small θ. 4. Discussion and results In this section we concentrate on the solution of the viscous wall layer III given by equations (2.24). We will primarily be interested in obtaining estimates for the skin frictions on the curved wall as it is these that, despite the gross simplifications 356 M. G. Blyth and A. J. Mestel (a) (b) Figure 8. Radial and azimuthal flow velocities in region III at x = 5.0 (on a radially exaggerated scale). (a) Contours of w 1 , which is positive everywhere and increases linearly with ζ as the core is approached. (b) Contours of v 1 . There is no swirl in the layer since v 1 > 0 everywhere, reaching a maximum near the corner. The fluid moves up along the curved wall and is carried downstream. that have been made in this model, are of most physiological interest. Now that the displacement function A∗ has been estimated, it is possible to proceed to obtain numerical approximations of the transforms of the flow variables and their derivatives in III. Once this has been done, a fast Fourier transform (FFT) routine is used to invert them numerically. From (2.24a) and (2.26) we have τ∗ax = u∗1ζ |ζ=0 = 21 µ0 ωA∗ + µ1 k 2 p∗1 /ω 4 , (4.1) for the axial skin friction, where µ0 = 3Ai(0) + Ai′ (0)/Ai(0), and µ1 = (6Ai(0))−1 . Also, from (2.24b) τ∗az = v ∗1ζ |ζ=0 = −2µ1 p∗1θ /ω, (4.2) for the azimuthal skin friction. Then (4.1) and (4.2), once inverted, will provide estimates of the down-pipe and azimuthal shear stresses respectively. In order to proceed further we need to solve (2.26) with (2.27) and (2.38) for the pressure p∗1 . We find that Z o n Z π/2 cosh(k(θ − π/2)) π/2 ∗ ∗ ∗ A (k, t) cosh(kt) dt , p1 = C A (k, t) sinh(k(θ − t)) dt + sinh(kπ/2) 0 θ (4.3) ′ 5 ∗ where C = 3Ai (0)ω /2k, which can then be calculated given A . We compute each transform for a range of θ values and then transform back with the FFT (see for example Brigham 1988). Contours of radial and azimuthal velocity perturbations in the viscous wall layer are shown in figure 8. The flow does not recirculate in this region since v 1 is always positive. The secondary flow in the core is sketched in figure 9(b). Down-pipe pressure gradients in the core and viscous wall layer are plotted in figure 10. Figure 11 shows plots of down-pipe and azimuthal shear stress perturbations for different θ. It is noticeable that around x = 0 there is a region where the axial shear stress perturbation has little azimuthal dependence. We now consider the asymptotic behaviour of the flow variables as x → ±∞. Since by (3.1d) as x → ∞, p1 ∼ x−3/2 we can deduce that A∗ ∼ k −3/2 as k → 0. Then it follows that for small k u∗1ζ |ζ=0 ∼ k −7/6 , v ∗1ζ |ζ=0 ∼ k −1/6 , p∗1 ∼ k −11/6 . (4.4a) (4.4b) 357 Steady flow in a dividing pipe 0.05 0 – 0.05 – 0.10 (b) p1h – 0.15 x increasing – 0.20 – 0.25 – 0.30 (a) 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 h Figure 9. (a) Azimuthal pressure gradients in III for x = 0.25, . . . , 1.5 in steps 0.25. (b) Schematic diagram of the secondary flow. Fluid is injected into the core from the plate layer I and swept around the curved wall. The radial velocity is positive throughout the core and the wall layer III, and mass conservation is maintained by the fluid being flushed downstream. Whereupon, as x → ∞ u1ζ |ζ=0 ∼ x1/6 , v 1ζ |ζ=0 ∼ x−5/6 , p1 ∼ x5/6 . (4.5a) (4.5b) How then do the variables behave as x → −∞? In the core if we assume the decay is algebraic in x, so p1 ∼ x−α , α > 0, we have (2.5) without the x-derivative term and boundary conditions p1 = 0 on r = 1 and p1θ = 0 on θ = 0, π by (2.16). In the absence of any eigensolutions to this problem we deduce that p1 = 0 and so it seems that for large negative x the decay in p1 will be exponential. Since regions II and III are linked through the displacement function A(x, θ) we would consequently expect everything there to decay in x at the same rate. We therefore propose in the core that p1 = eλx f(r, θ) + · · · as x → −∞, (4.6) and correspondingly in the viscous wall layer u1 = eλx u−∞ (r, θ) + · · · as x → −∞, (4.7) with similar expansions for the other variables in III. Then eigenvalues λ correspond to eigenfunctions f being solutions of ∇22 f + 4r fr + λ2 f = 0, 1 − r2 (4.8) with f = 0 on r = 1 and fθ = 0 on θ = 0, π. We Fourier analyse in θ, writing f= ∞ X m=0 fm (r) cos mθ, (4.9) 358 M. G. Blyth and A. J. Mestel 0 (a) – 0.1 – 0.2 p1x – 0.3 – 0.4 –2 0 2 4 6 10 8 x 0.0250 (b) 0.0125 0 p1x – 0.0125 – 0.0250 – 0.0375 –4 –2 0 2 4 6 8 x Figure 10. Down-pipe variation of pressure gradient perturbations. (a) In the viscous wall layer, p1x at θ = π/2. (b) In the core, p1x at (r, θ)=(0.5, π/2). which automatically satisfies the second condition on f. Then   4r3 2 ′′ r fm + r + fm′ − (m2 − λ2 r2 )fm = 0, 1 − r2 (4.10) 359 Steady flow in a dividing pipe 1.2 1.0 0.8 sax 0.6 h decreasing 0.4 0.2 (a) 0 –1 1 0 2 3 5 4 x 0 – 0.05 h decreasing – 0.10 saz – 0.15 – 0.20 – 0.25 – 0.30 (b) – 0.35 –4 –2 0 2 4 6 8 10 x Figure 11. Variation of wall stress perturbation. (a) Axial τax . (b) Azimuthal τaz . For both θ decreases in equal steps in the range π/4 6 θ 6 π/2. for which we impose that fm (0) be regular and fm (1) = 0. We solve this problem using the Fortran library Nag routine D02KEF. It is observed that when m = 0 the eigenvalues correspond to the zeros of J1 (r), the first-order Bessel function. This suggests an exact solution is possible in this case and indeed we find that for m = 0 360 M. G. Blyth and A. J. Mestel 4.0 (a) 3.5 m 3.0 = 0 2.5 fm 2.0 1.5 = 6 1.0 m 0.5 0.1 0 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 r –4 (b) –5 –6 –7 –8 –9 –10 –11 –12 –13 –2 – 1.8 – 1.6 – 1.4 – 1.2 – 1.0 – 0.8 – 0.6 – 0.4 – 0.2 0 x Figure 12. (a) Eigenfunctions fm for m = 0, . . . , 6. (b) Down-pipe comparison of log p1 at (r, θ) = (0.5, π/2) in the core (——) with a line of gradient 3.83 (· · ·). the eigenfunction is given by f0 = J0 (λr) + r2 J2 (λr), (4.11) where J1 (λ) = 0. Note that this solution appears in slightly different form in Smith Steady flow in a dividing pipe 361 (1976c). Despite extensive efforts, we have been unable to find an analytic solution to (4.10) for general m. The numerical solutions predict that the first non-zero eigenvalues λm0 for each mode increase more or less linearly with m. Figure 12(a) shows the eigenfunctions fm for the first few modes. We are interested in the smallest eigenvalue as this will then determine the slowest rate of decay of p1 . In fact we can prove that the smallest eigenvalue should occur when m = 0. If we write (4.10) in self-adjoint form !′ λ2 r2 − m2 rfm′ fm = 0, (4.12) + 2 2 (1 − r ) r(1 − r2 )2 then we can use the Rayleigh–Ritz method (see for example Duff & Naylor 1966) whereby the Rayleigh quotient Z 1 (p(wt′ )2 − qwt2 ) dr 0 (4.13) R= Z 1 2 swt dr 0 is minimized over a set of l trial functions wt which obey the boundary conditions on fm to obtain upper-bound estimates Λmj (j = 1, . . . , l) for the first l non-zero eigenvalues λmj corresponding to mode m with equivalence between the two occurring when wt is chosen to be the eigenfunction precisely. Here we have p = r/(1 − r2 )2 , q = −m2 /r(1 − r2 )2 , s = r/(1 − r2 )2 , (4.14a) (4.14b) (4.14c) λ00 < λ10 < · · · < λm0 < · · · , (4.15) and so p, s > 0, q 6 0 for 0 6 r 6 1, ∀m. We wish to estimate λ00 . In general we know = Λm+1 that λm0 = Λm0 if we choose wt = fm as a single trial function. Similarly λm+1 0 0 m m m if wt = fm+1 . Then if we select wt = fm+1 to estimate λ0 we have λ0 6 Λ0 6 λm+1 . 0 Therefore and so the first non-zero eigenvalue of the zeroth mode is the smallest. We conclude then that as x → −∞ 0 p1 = eλ0 x f(r, θ) + · · · as x → −∞, (4.16) with λ00 ≈ 3.8317. How does this compare with our numerical calculations? Figure 12(b) shows a logarithmic plot of p1 beside a straight line of gradient 3.83. We can see that even for relatively small x, p1 does seem to be decaying exponentially at approximately the predicted rate. 5. Concluding remarks We have extended the work of Smith (1977a) from the large down-pipe limit to solving over the entire range of the developing flow. The effect of larger constrictions and blockages on oncoming developed flow is considered by Smith (1977b) and those in a developing entry flow by Smith (1976b). In both of these cases, for pipe and channel flows under the correct conditions upstream influence may be significant and wide-ranging. However, with reference to the comment made in §2.1 that upstream 362 M. G. Blyth and A. J. Mestel influence may be anticipated, we have found in fact that it is minimal as the pressure in the core decays exponentially upstream. From a physiological point of view, bearing in mind the comments made in the introduction, it is interesting to note that on the curved wall the axial skin friction perturbation τax increases from its undisturbed Poiseuille value fairly sharply as the splitter plate is reached and thereafter increases algebraically downstream: τax = 1 2 + R −1/6 τax , (5.1) where τax ∼ x1/6 as x → ∞ by (4.5a). To leading order the axial shear stress near the middle of the splitter plate will approximately be that of the classical Blasius boundary layer, 1 (5.2) τBax = R 1/2 √ f ′′ (0)(1 − z 2 )3/2 x−1/2 + O(1) 8 2 with f ′′ (0) ≈ 0.4696. So the axial skin friction is lower on the curved wall than on the plate, with the two becoming of the same order of magnitude when x ∼ R. Of course in a dividing pipe with a non-zero angle of bifurcation, separation may occur on the upper wall which would affect the shear stress there dramatically. In this simplified scenario however, no such separation is encountered. There are many possible avenues of extension of the current study. Relaxation of some of the simplifications mentioned in the introduction to make it of more biological relevance could constitute much future work. The problem of merging flow over the trailing edge of a splitter plate might be considered. This has been examined in the two-dimensional case of a channel by Badr et al. (1985). They solved the flow over the trailing edge of the plate both numerically and analytically. Working in the large Reynolds number limit, they studied the flow on a lengthscale x = R 1/7 X, where X is O(1), surrounding the trailing edge as suggested by Smith (1977b). They then matched this to the following wake on an O(R) lengthscale, where a three-layer structure emerges, with a core flow sandwiched between two viscous layers which interact through a pressure-displacement relation in the core. Finally we note that the work presented here could fairly easily be adapted to pipes of different cross-section, for example an ellipse, which might better represent the shape of an artery. The assumption of up/down symmetry could be relaxed as the local solution to the core given in §2 can be simply altered for a corner of any angle. The analysis also extends to other oncoming uni-directional profiles. Appendix Here we detail the functions g, κi which are needed for the matching of region IV to regions I, II and III (see § 2.4). For i = 1, . . . , 4: √ (A 1a) κ1 = 3β sec(γ/2)x1/2 H(x)/ 2, √ 3 −1/2 κ2 = 2βcosec(γ) sin (γ/2)x H(x), (A 1b) √ −1/2 H(x)/ 2, (A 1c) κ3 = 3β sin(γ/2)x √ 3 −3/2 (A 1d) H(x)/(5 2), κ4 = −β sin (γ/2)(2 cos(γ) + 3)x where H(x) is the Heaviside function. Also, 3 g(x) = κ4 |(γ=π/2) = − 20 βx−3/2 H(x). (A 2) Steady flow in a dividing pipe 363 For i = 5, . . . , 8: κ∗5 = 32 Q nZ ωζ 0 o Ai(q) dq − Ai′ (0)M(ωζ) , κ∗6 = 6ω 3 QAi′ (0)M(ωζ), Z ζ Z ωζ ′ ∗ 3 Ai(q) dq dζ ′ , κ7 = 2 ikQ 0 where (A 3a) (A 3b) (A 3c) 0 √ κ8 = −βAi′ (0)Γ ( 16 )x−7/6 H(x)/(25/3 2π), √ Q = −12 2λc k −3/2 , ω = ( 21 ik)1/3 , (A 3d) (A 4) and where an asterisk denotes a Fourier transform according to (2.23). These should be inverted numerically. REFERENCES Badr, H., Dennis, S. C. R., Bates, S. & Smith, F. T. 1985 Numerical and asymptotic solutions for merging flow through a channel with an upstream splitter plate. J. Fluid Mech. 156, 63–81. Bennett, J. 1987 Theoretical properties of three-dimensional interactive boundary layers. PhD thesis, University College, London. Blyth, M. G. 1999 Steady flow in dividing and merging pipes. PhD thesis, Imperial College, London. Brigham, O. R. 1988 The Fast Fourier Transform and its Applications. Prentice-Hall. Caro, C. G., Fitz-Gerald, J. M. & Schroter, R. C. 1971 Atheroma and arterial wall shear: Observation, correlation and proposal of a shear dependent mass transfer mechanism of atherogenesis. Proc. R. Soc. Lond. B 177, 109–159. Cowley, S. J. 1983 Steady flow through distorted rectangular tubes of large aspect ratio. Proc. R. Soc. Lond. A 385, 107–127. Dhanak, M. R. & Duck, P. W. 1997 The effects of freestream pressure gradient on a corner boundary layer. Proc. R. Soc. Lond. A 453, 1793–1815. Duff, G. F. D. & Naylor, D. 1966 Differential Equations of Applied Mathematics. John Wiley & Sons. Friedman, M. H., O’Brien, V. & Ehrlich, L. W. 1975 Calculations of pulsatile flow through a branch. Implications for the haemodynamics of atherogenesis. Circ. Res. 36, 277–285. Gammack, D. & Hydon, P. E. 1999 Steady flow in pipes of non-uniform curvature and torsion. J. Fluid Mech. (submitted). Lynch, D. G., Waters, S. L. & Pedley, T. J. 1996 Flow in a tube with non-uniform, time-dependent curvature: governing equations and simple examples. J. Fluid Mech. 323, 237–265. Lyne, W. H. 1970 Unsteady viscous flow in a curved pipe. J. Fluid Mech. 45, 15–31. Pal, A. & Rubin, S. G. 1971 Asymptotic features of viscous flow along a corner. Q. Appl. Maths 29, 91–108. Pedley, T. J. 1980 The Fluid Mechanics of Large Blood Vessels. Cambridge University Press. Rieu, R. & Pelissier, R. 1991 In vitro study of a physiological type flow in a bifurcated vascular prosthesis. J. Biomechanics 24, No.10, 923–933. Rubin, S. G. 1966 Incompressible flow along a corner. J. Fluid Mech. 26, 97–110. Rubin, S. G. & Grossman, B. 1971 Viscous flow along a corner: numerical solution of the corner layer equations. Q. Appl. Maths 29, 169–186. Smith, F. T. 1975 Pulsatile flow in curved pipes. J. Fluid Mech. 71, 15–42. Smith, F. T. 1976a Pipeflows distorted by non-symmetric indentation or branching. Mathematika 23, 62–83. Smith, F. T. 1976b On entry-flow effects in bifurcating, blocked or constricted tubes. J. Fluid Mech. 78, 709–736. Smith, F. T. 1976c Flow through constricted or dilated pipes and channels: Part 2. Q. J. Mech. Appl. Maths 29, 365–376. 364 M. G. Blyth and A. J. Mestel Smith, F. T. 1977a Steady motion through a branching tube. Proc. R. Soc. Lond. A 355, 167–187. Smith, F. T. 1977b Upstream interactions in channel flows. J. Fluid Mech. 79, 631–655. Smith, F. T. 1979 Instability of flow through pipes of general cross-section. Part 1. Mathematika 26, 187–210. Strikwerda, J. C. 1989 Finite Difference Schemes and Partial Differential Equations. Wadswirth & Brooks/Cole. Walton, A. G. & Smith, F. T. 1997 Concerning three-dimensional flow past a tall building on flat ground. Q. J. Mech. Appl. Maths 50, 97–128. Whitmore, R. L. 1968 Rheology of the Circulation. Pergamon Press, Oxford. Wille, S. O. 1984 Numerical simulations of steady flow inside a three dimensional aortic bifurcation model. J. Biomed. Engng. 6, No. 1, 49–55. Zabielski, L. & Mestel, A. J. 1998a Steady flow in a helically symmetric pipe. J. Fluid Mech. 370, 297–320. Zabielski, L. & Mestel, A. J. 1998b Unsteady blood flow in a helically symmetric pipe. J. Fluid Mech. 370, 321–345. Zamir, M. & Roach, M. R. 1973 Blood flow downstream of a two-dimensional bifurcation. J. Theor. Biol. 42, 33–48.