## Abstract

Multiply branching fluid flows are modelled in two contexts. The first (type I) is for one-to-many branching. Computations are described for flow through a channel, with fully developed motion upstream, which branches abruptly into a number of subchannels downstream. The differences in pressure between the upstream end of the channel and the downstream ends of the subchannels are substantial. Comparisons with recent analytical predictions show fair agreement for Reynolds numbers in the low tens and above. The second context (type II) has successive generations of bifurcation in a network. Modelling, computations and analysis include the effects of many bifurcations.

## 1. Introduction

There are two main kinds of multi-branching: (I) a single one-to-many branching and (II) a succession of generations of bifurcations (i.e. of one-to-two branchings). Type II is well known because of its relevance to comparatively normal networks such as in the lung, cardiovascular or cranial systems. Type I is also of practical significance in relation to arteriovenous malformations (cerebral abnormalities in which a blood vessel subdivides over a relatively short distance into many subvessels; Hademenos *et al*. 1996; Smith & Jones 2003) and to emerging issues, including increased mass flux and an apparent absence of flow reversal. Several theoretical conjectures remain largely untested, while direct computation has only given solutions either accurately over a few generations or with limited detail over more generations (see Hademenos *et al*. 1996; Pries *et al*. 1998; Goldman & Popel 2000; Smith & Jones 2003). Related branching or network studies have been conducted with some success by Handa *et al*. (1993), Miyasaka *et al*. (1993), Pedley *et al*. (1994), Gatlin *et al*. (1995), Young *et al*. (1996), Gao *et al*. (1997), Hademenos & Massoud (1997), Wilquem & Degrez (1997), Zhai *et al*. (1997), Pries *et al*. (1998), Brada & Kitchen (2000), Kassab *et al*. (2000), Lorthois *et al*. (2000), McEvoy *et al*. (2000), Cassidy *et al*. (2001) and Comer *et al*. (2001). The presence of an arteriovenous malformation in a cranial system eventually requires study of the application of type I within the wider framework of an application similar to type II. The significant Reynolds numbers in such physiology vary from quite small to quite large values, although no more than a few hundred or thousand, whereas much larger values are often encountered in industrial processes, such as in oil pipeline design and water dispersal networks.

To test conjectures concerning the one-to-many branching (type I), we apply direct numerical simulation and make comparisons with theoretical predictions. This was designed to gain insight for theory and computation, to understand flow physics, and to quantify the variation in flux from one daughter to the next. We are also interested in the consequences for total through flow, individual fluxes and wall shears, given that multi-branching can produce significant alterations in pressure and flow direction (as in impacting motions). Moreover, it is useful to discover the range of Reynolds numbers over which the approximations in Smith & Jones (2003) apply satisfactorily. Simulations and theoretical work by Tadjfar & Smith (2004) for three-dimensional motion show fair agreement, even at relatively low Reynolds numbers and high divergence angle. The response at a branching, especially of type I, can be a vital component in a cranial or cardiovascular network (see references cited above), and this emphasizes the need to assess the accuracy of local branching analysis. Concerning networks (type II), we also apply a combination of local inviscid and global viscous flow balancing similar to that in Smith *et al*. (2003*b*). The local part is nonlinear, connected with the Smith & Jones study, whereas the global part is linear, akin to Hademenos *et al*. (1996). Thus, a nonlinear system of network equations holds.

The two contexts (type I and type II) are investigated here as largely separate issues, even though they have a number of common features. Thus, multi-branching junctions are treated by direct simulations and by inviscid theory, whereas longer-scale networks bring in the additional influences of viscous dominated vessel flows. Section 2 describes the configuration in the type I context for steady motion (given that the cranial blood flow of interest here is generally much steadier than elsewhere in the arterial system for instance) along with the computational method and results for various multi-branching geometries. This is followed by comparisons with theoretical predictions. These are for two or more daughters, which we assume to be parallel channels. We also note that the Smith & Jones (2003) work implies that even a few daughters yield properties close to those for many daughters. Section 3 then describes our treatment for the type II context, in which there is an abundance of parameters and flow phenomena. Section 4 provides further comments.

## 2. The model for one-to-many branching

Near a multi-branching junction, the flow response over streamwise lengths comparable to the tube width is expected to matter most when the Reynolds number (*Re*) is medium to large. To test this expectation, we first describe a numerical investigation of the full viscous problem for such a junction. The velocities, lengths and pressures non-dimensionalized with respect to *u*_{0}, *ℓ*, are (*u*,*v*), (*x*,*y*), *p*, respectively, and *Re* is *u*_{0}*ℓ*/*ν*, where *ρ*, *ν* denote the fluid density and kinematic viscosity in turn, *u*_{0} is a representative measure of flow speed and *ℓ* is typical of a tube half width. As the characteristic pressure difference Δ*p* can be taken to be prescribed, we set *u*_{0} as (|Δ*p*|/*ρ*)^{1/2}, which leaves *Re* as (*ℓ*^{2}|Δ*p*|/*ν*^{2}*ρ*)^{1/2}. The Navier–Stokes equations then govern the internal two-dimensional steady laminar motion:(2.1)(2.2)(2.3)where ∇^{2} denotes the Laplacian operator (∂^{2}/∂*x*^{2} + ∂^{2}/∂*y*^{2}). The system is subject to the boundary conditions of no slip at all the fixed solid surfaces and of unidirectional flow sufficiently far upstream. Thus,(2.4)(2.5)(2.6)

Solutions of the full viscous problem in equations (2.1)–(2.6) at finite *Re* values are presented in figure 1*a*–*e*. Here, *x*_{−∞}, *x*_{∞} are the end station values, suitably far upstream and downstream, respectively, while the zeros of ∂*u*/∂*x*, *v* in equations (2.4) and (2.5) correspond strictly to those stations being infinitely far upstream and downstream. This allows necessary flexibility regarding the inflow and outflow, although near Poiseuille flow is generally found to emerge anyway in the results at medium *Re* values. The values *p*_{−∞}, *p*_{∞} are the prescribed end pressures upstream and downstream, *N* is the number of daughter channels, and the walls of the mother channel upstream are given by *y*=±1, for *x*<0. The branching itself occurs at *x*=0. In each daughter, where *x*>0, the coordinate *y* may be defined as perpendicular to the daughter wall. The mass flux in total through the mother and through the daughters is unknown. The branching geometry is to allow for expansions and contractions but by means of straight sections of channel. The branching flow was treated by a direct finite difference method based on Dennis & Hudson (1995) and Smith *et al*. (2003*a*). The end stations were usually taken at −1, 1 as fundamental cases. Reduced downstream areas were obtained (i.e. *A*<1) by closing off some of the downstream area nearest the outer wall, effectively shutting off a daughter tube.

A typical case required around 500 iterations to converge to within 10^{−8}. As *Re* increased, substantial under-relaxation was necessary to achieve a converged solution. The grid size was generally taken to be 0.05 in both directions, but refinement was used to ensure accuracy of the solution. To help resolve the thin wall layers present, the results below *Re*=200 have a grid size refined to 0.02, although they are virtually identical to those from the coarser mesh in terms of velocity profiles and through flow.

The results will be compared with the limit analysis in Smith & Jones (2003), which investigates a nonlinear inviscid model for high *Re*. The model holds formally over the short length-scale where |*x*| is *O*(1). It seems reasonable to neglect viscous effects to a first approximation near a junction, as *Re* increases, and indeed this is exploited in the work of the next section. Longer length-scales, with distances |*x*| of 0(*Re*) or more, lead to viscous filling of the daughter tube flows. This filling yields the fully developed velocity profile, which helps to form the starting condition for the short-scale problem. Over the short length-scale, recurrence relations between the unknown flow downstream in the daughter tubes and the incoming rotational flow in the mother are derived by balancing the mass flux and the pressure head since they remain constant along each streamline. The relations are then analysed in detail to determine the flows in the individual daughters as well as the total flow rate. The analytical predictions complement those of the current study. Thus, figure 1*a* compares the through flow generated over a range of downstream areas *A* as predicted by analysis (Smith & Jones 2003) and by the current finite *Re* method. It shows the case *N*=11 for *Re* ranging from 10 to 60. The comparisons indicate fair agreement as *Re* increases. Further comparisons are presented in Smith *et al*. (2003*a*), including cases where flow is forced down a daughter tube against the pressure gradient; that is, with one daughter pressure higher than the upstream pressure while the other daughters have a low pressure. Forward flow is obtained in all daughters in some cases. This phenomenon is also found in the inviscid theory. Figure 1*b* shows velocity profiles obtained at *Re*=70 for a case with *N*=5. Figure 1*c–e* presents results for a mother splitting into just two daughter branches at *Re*=200. Figure 1*c* shows the *u* velocity profiles for various divider positions with downstream pressures being equal in each daughter. At these medium Reynolds numbers near Poiseuille flow emerges at the inlet and, after some deformation as the branch splits, the solution appears to settle back into near Poiseuille flow very quickly after the bifurcation. We note here that such a fully developed motion upstream and downstream is not assumed a priori; it emerges at these *Re* values. At large *Re*, the downstream influence length (or entry length) has the well-known long scale 0(*Re*), whereas the upstream influence length is usually of order *Re*^{1/7} (Smith 1977), consistent with the present numerical results. Figure 1*d* provides the *u* profiles with equal downstream sizes but varying the lower daughter pressure from *p*=0 (the symmetric case) to *p*=−30. As might be expected, when *p* is decreased the flow through that daughter, and the entire system, is increased. Finally, figure 1*e* displays the through flow in each case. For the varying area ratios, the nearer to the symmetric case the less through flow is generated for a given pressure distribution. In the varying pressure case, the through flow increases seemingly linearly with the magnitude of the overall pressure drop.

## 3. Modelling of successive bifurcation networks

For the network (type II; figure 2*a*), we presume the pressure drop between the original entrance and the end vessels furthest downstream is known, as well as the details of all vessels (resistance and diameter, etc.), and we seek the total flux. The flux and pressure drop through individual vessels can then be calculated. The generic model makes some central assumptions. First, the flow is planar. Second, a long–short scale split is exploited. The dynamics in individual vessels are over such a long scale that fluid inertia may be neglected, so that a fully developed Poiseuille profile holds. Hence, the pressure drop over a vessel is proportional to the flux through it, or to the mean velocity at its entrance. In a sense, that is the far-field view. In contrast, the flow locally at the junctions (i.e. in the near-field) is so (spatially) rapid as to be governed by inviscid dynamics as in Smith & Jones (2003). Here, vorticity and pressure head are conserved along a streamline, which, with mass conservation, allows calculation of the downstream velocity profiles (Smith & Jones 2003), along with nonlinear expressions for the pressure drops suffered on passing into the downstream vessels. These profiles then develop into the Poiseuille profile, owing to viscous action, on an intermediate scale that is short compared with the present vessel length (as noted in §2). The short–long dynamics are combined through the global condition of given total pressure drop from the original mother to any end vessel, giving a nonlinear system for the flux in every vessel. In principle, this system is solvable but here, in an attempt to understand the entire network, we make a third assumption of replacing the incoming Poiseuille profile at a given branch by a uniform plug flow with the same mass flux. This assumption approximates the second above as Smith & Jones (2003) shows that it emerges naturally for increasingly large numbers of downstream branches for *N* above approximately three. Additional comparisons by N. C. Ovenden (personal communication) show that this simplification gives results similar to those from the full equations in cases where the flow remains forward and, equally important, is asymptotically correct in the limit of an equal division of the oncoming mass flux into the downstream vessels. In addition, we restrict consideration to junctions with net cross-sectional area decrease to prevent separation (see figure 1). Nevertheless, area increase can be accommodated here by means of the longer-scale (far-field) geometry, where the model allows slow variation in an individual vessel width over its length where the flow remains attached, thereby accommodating networks where the net area increases.

The network thus consists of a concatenation of units, with a single parent vessel upstream. A unit has the downstream end of that parent vessel, splitting into two over a short scale, and the full extent of the two downstream vessels. A vessel identified by subscript *n* (with the original mother having *n*=1) has upstream width *s*_{n} and downstream width *σ*_{n}. The pressure drop along its length (effectively in the far-field) is *U*_{n}*r*_{n} say, where *U*_{n} is the cross-sectionally averaged velocity at the upstream end and(3.1)is the vessel resistance. The integration is along the vessel length *l*_{n}, measured by *x*, and *s*_{n}(*x*) is the local cross-section. Lengths are normalized on the network's mother vessel cross-section *d**, so that *s*_{1}=1. We require *l*_{n}≫1, and now *Re*=*d***U**/*ν*, so that *l*_{n}/*Re* is presumed to be 0(1) or more. The scale *U** is implicitly defined through an insistence that *r*_{1}=1. The far-field result in equation (3.1) then stems from the Reynolds lubrication equation. In contrast, across a junction (the near-field) the pressure head *p*+*U*^{2}/2 is constant, where *U* is the averaged velocity. The pressure drop between the downstream ends of the parent, with subscript *n*, and of a daughter, with subscript *m*, is(3.2)where *e*_{n}=*σ*_{n}/*s*_{n}, which is greater than unity if a vessel widens along its length. The similar drop for the other daughter vessel, with subscript *l*, and mass conservation respectively give(3.3)

These equations determine *U*_{n}, given a knowledge (or guess) of Π_{m}, Π_{l}. Knowing the total pressure drop across the network allows the simultaneous calculation of these local pressure drops and the velocities in each vessel. The resulting system is quadratically nonlinear and may be solved by Newton iteration.

Analytical progress may be made if the imposed pressure drops across the network are identical and the vessels at each generation *i* are similar. The limit *N*→∞, with *N* now the total number of generations and *i*=*qN*, 0≤*q*<1, , , , yields the ordinary differential equations , , with *U*_{i}=*u*(*q*) and with *P*(*q*) denoting the pressure drop between the mouth and the current generation, as illustrated by equations (3.1)–(3.3). The change in subscript is to distinguish between generations rather than between vessels. The velocity *U*_{0} into the parent is therefore related to the pressure drop Δ*P* across the network through the quadratic relation(3.4)in this simple example. In less simple cases, computation is necessary.

In two calculations here, based on equations (3.1)–(3.3), a single input vessel splits into 32 downstream (see figure 2*a*). In figure 2*b*, the viscous resistance *r* for each vessel is half that at the previous generation, and at each bifurcation, the total cross-section of the network is reduced by a factor of 0.9. We should recall that the individual vessel resistance depends on both width and length. Further, the pressure drop across the left-hand side of the tree is fixed at 5, while that across the right, Δ*P*_{3}, is varied. The fluxes through each half, *Q*_{2}, *Q*_{3}, and through the whole network, *Q*_{1}, are plotted together with the velocities into each vessel *U*_{1,2,3}. The minimum Δ*P*_{3} is approximately 1.25, corresponding to no flow through the right side (zero flow at the first bifurcation). The upper value is approximately 6.65. Above this value, the flow decelerates into the left half and we can expect flow separation. In figure 2*c*, the pressure drop is 5 for all vessels, but in the right half each vessel has its cross-section increased by a factor *F*. The vessel lengths are varied so that the viscous resistance remains unaltered. The velocities at corresponding generations through the two halves of the network remain equal, although the fluxes do not. As *F*→0, the flux into the right branch vanishes as expected, although the velocity does not. For *F* >1.22, the flow decelerates into both halves at the first junction.

## 4. Further comments

Although this article has treated the cases of a multi-branching junction (type I) and a bifurcating network (type II) as almost separate issues, they have common features as well as distinctions. First, the comparisons of type I in §2 are a potential boost to the theoretical and numerical aspects of the modelling. Next, forward flow can occur within a daughter despite a pressure rise, according to both theory and the present simulations. This may seem counter-intuitive; however, it arises because of combined pressure gradient effects on the complete motion. Under different conditions, those effects can provoke back flow according to the simulations, in which case theory (Smith & Jones 2003) requires modification owing to the vorticity associated with the back flow. Again, on broader scales, the mother and daughter end pressures depend on longer-scale pressures whenever the total vessel lengths exceed the representative viscous length-scale, a feature in common with the type II context. Then, local pressure differences (as discussed in §2) can readily be set up by the total differences imposed over the longer scale. In the type II context, on the other hand, the model cannot allow the above counter flow, and issues remain on such junction properties. Nevertheless, simulations as in figure 1 provide some guidance. Moreover, the network response altogether is perhaps not unlike that in type I. Naturally, the real physiological and other systems, such as those with arteriovenous malformations and cranial networks, are extremely complicated. Although, as remarked earlier, cranial blood flow may be relatively steady, three-dimensional effects are important, as is wall flexibility. The specific flow solutions given in this study, coupled with earlier analyses and extensions (e.g. inserting context type I in a type II network), are aimed at promoting wider understanding of these complicated systems.

## Acknowledgments

Thanks are due to NSERC of Canada and EPSRC of the UK for financial support.

## Footnotes

One contribution of 19 to a Theme ‘New developments and applications in rapid fluid flows’.

- © 2005 The Royal Society