A generalized optimization principle for asymmetric branching in fluidic networks.

When applied to a branching network, Murray's law states that the optimal branching of vascular networks is achieved when the cube of the parent channel radius is equal to the sum of the cubes of the daughter channel radii. It is considered integral to understanding biological networks and for the biomimetic design of artificial fluidic systems. However, despite its ubiquity, we demonstrate that Murray's law is only optimal (i.e. maximizes flow conductance per unit volume) for symmetric branching, where the local optimization of each individual channel corresponds to the global optimum of the network as a whole. In this paper, we present a generalized law that is valid for asymmetric branching, for any cross-sectional shape, and for a range of fluidic models. We verify our analytical solutions with the numerical optimization of a bifurcating fluidic network for the examples of laminar, turbulent and non-Newtonian fluid flows.

Introduction

The optimal branching of fluidic networks has been the subject of numerous studies owing to its importance in understanding the behaviour of biological vessels and for the biomimetic design of artificial systems. Much of the research stems from Murray’s law [1], who posited that there were two main energy requirements for blood flow through a cylindrical vessel of radius R: (i) the energy required to overcome viscous drag and drive the flow and (ii) the energy required to metabolically maintain the fluid and vessel. Assuming the flow to be laminar, Newtonian, steady, incompressible and fully developed, Murray used the Hagen–Poiseuille equation to model the driving power requirement (i.e. the driving power is proportional to R−4), and assumed that the maintenance power was proportional to the volume of the channel (i.e. proportional to R2). Applying the principle of minimum work to the total power requirement, Murray surmised that in an optimized cylindrical channel, the volumetric flow rate Q is proportional to R3. By applying mass conservation over a branching network, and assuming that local pressure losses through the junction (owing to bends and channel contractions) are negligible compared with the pressure losses over the channel lengths, this principle is most commonly expressed as a power law between a parent channel and N daughter channel branches where the subscripts p and d denote the parent and ith daughter, respectively. Although originally targeting blood transport through the cardiovascular system [2-9], experimental data have shown Murray’s law to be a decent approximation for a number of other biological networks, e.g. in the bronchial trees of humans and dogs [10-12]; in the chick embryo [13], and in the leaf veins of plants [14-16].

Whole blood (plasma and cells) is a non-Newtonian fluid that exhibits shear-thinning behaviour, i.e. its viscosity decreases with increased shear-strain rate. To more accurately model vascular networks, Murray’s law has been applied to non-Newtonian fluids [17,18] using the popular power-law fluid model [19,20]. In both studies, it was found that the optimal radius relation is unaffected by shear-thinning or shear-thickening fluid behaviour, and equation (1.1) is maintained for the whole range of non-Newtonian fluids. Murray’s law has also been applied to turbulent flows [21-23], which can be found in the upper airways of the lung [24] and, under some circumstances, in blood flow through the aorta [25], as well as in a number of hydraulic and pneumatic civil engineering applications. For fully rough-wall turbulent flow, the flow rate was found to be proportional to R7/3, leading to the relation for branching networks.

While Murray’s law has been most often applied to circular or elliptical [18] cross sections (owing to the shape of biological networks), optimized branching is also useful for the design of artificial systems [5], such as for fuel cells [26] or heat exchangers [27,28], which are often constrained to certain shapes by manufacturing procedures. To this end, Murray’s law has been adapted for networks of rectangular and trapezoidal cross sections [29].

Although originally derived from the principle of minimum work, it has been noted that the application of other optimization principles results in the same relationship between flow rate and channel radius: minimizing the total mass of the channel [22], minimizing volume for a constant pressure drop and flow rate [30], minimizing pumping power [31], maintaining a constant shear stress in all channels [32], and minimizing flow resistance for a constant volume [5,23,33].

However, despite many developments to Murray’s law, we submit that it is, in fact, suboptimal for asymmetric branching. In this paper, we derive a generalized law that is applicable to symmetric and asymmetric branching, for any cross-sectional shape, and for a range of fluidic models (e.g. Newtonian and non-Newtonian, laminar and turbulent).

Analytical solutions

The conditions for optimal branching can be generalized as a maximization of flow conductance per unit volume through each branch of the network, for a variety of constraint combinations. As with Murray’s law, we assume the flow to be steady and fully developed. For a two-level network (consisting of a single parent channel branching into multiple daughter channels), this can be expressed as where Q is the volumetric flow rate through the parent channel, is the set of pressure drops ΔP over each of the N network branches—from the inlet of the parent to the outlet of the corresponding daughter—V is the network volume, and is the jth daughter–parent cross-sectional area ratio. Note, for each constraint option, two parameters are fixed and the third parameter is optimized (volume minimization, flow-rate maximization and pressure-drop minimization, respectively). Therefore, all three constraint options lead to an identical optimal relation because regardless of the constraints chosen. Note, noting that the last term means that d(ΔP)/dΓ=0 for all i. For our optimization, we assume that the channel lengths L are large compared with the size of the parent–daughters junction, so that (i) the localized pressure losses at the junction are negligible compared with the pressure drops over individual channels (as in Murray’s law), and (ii) the volume of the network can be considered to be the sum of the channel volumes If we consider the optimization of the jth daughter channel, then inserting equation (2.4) into the fitness function of equation (2.1), and noting equation (2.3), gives The pressure drop over the parent and each of the i daughter channels can be expressed in terms of the volumetric flow rate and where Ψ=Qd/Q is the fraction of the total flow rate taken by the ith daughter channel and k is flow resistance per unit length, e.g. kp=ΔPp/(QLp). The pressure drop over each network branch ΔP=ΔPp+ΔPd is then Differentiating equation (2.8) with respect to Γ, and noting equation (2.3), gives Substituting equation (2.9) into equation (2.5), via the chain rule, gives our generalized optimal area ratio which, for brevity, will be referred to as the generalized law to distinguish it from Murray’s law. It should be noted that the subscript j is not present in equation (2.10), so this relationship is not specific to a particular daughter channel; it relates the properties of all daughter channels to that of the parent.

This generalized law is valid for any cross-sectional shape, for any fluid (e.g. non-Newtonian) and for any Reynolds number (e.g. for turbulent flow). It is also valid for flows through nanoscale networks where the fluid is dominated by velocity slip at the walls; however, in the paper, we restrict our attention to the continuum-flow limit. We now consider some important cases where A can be expressed easily as an analytical function of k.

Laminar flow

The steady-state incompressible Navier–Stokes momentum equation describes laminar flow through a long channel with an arbitrary cross-sectional shape, i.e. where μ is the dynamic viscosity (constant for a Newtonian fluid) and u is the streamwise channel velocity. This can be non-dimensionalized using ΔP/L, μ and cross-sectional area A, such that where and tilde denotes a dimensionless quantity or operator. The axes of the cross-sectional plane are defined as y,z and Provided the boundary conditions are fixed (which is the case for the continuum-flow limit, where the no-slip boundary condition applies), the solution of equation (2.12), , is independent of A, ΔP, L and μ, and is thus a property of the cross-sectional shape alone. Similarly, so is An expression for the volumetric flow rate is obtained by integrating the fluid momentum over the cross-sectional area Substitution of equations (2.13)–(2.15) into (2.16) gives the volumetric flow rate through a channel with an arbitrary cross-sectional shape and flow resistance per unit length It is assumed that the pressure drop over the network is small such that the viscosity and density are constants for a Newtonian fluid. For cylindrical channels S=1/(8π), equation (2.17) becomes the Hagen–Poiseuille flow rate. If the cross-sectional shape is constant throughout the network, i.e. S=const., substituting (2.18) into the generalized law (2.10) and cancelling the constant terms gives From equation (2.7), it can be seen that for all daughter channels. Combining equations (2.18) and (2.20) produces the cross-sectional area relationship between the ith and jth daughter channels where is the pressure-gradient ratio between the jth and ith daughter channels. Note that, because the shape property S cancels in equation (2.21), the generalized law will be independent of the cross-sectional shape of the channels. Substituting equation (2.21) into equation (2.19) and rearranging for Γ as defined by equation (2.2) gives Equation (2.23) relates the area of the parent channel to the area of the jth daughter channel in an optimized two-level network of laminar flow. It is valid for any cross-sectional shape, provided the shape is constant through the network. Equation (2.23) is only equivalent to Murray’s law (which is when posed in terms of an area ratio) if the daughter channels branch symmetrically, i.e. Ψ=Ψ=1/N and Φ=1. By inserting these constraints into equation (2.23), the symmetric generalized law for laminar flow is This means that for symmetric branching, Murray’s law is valid for any cross-sectional shape, not just circles. The reason Murray’s law produces a suboptimal result for asymmetric branching is that it was derived to optimize a single channel in isolation. However, as shown above (and verified later), for asymmetric branching, the global optimum is not the same as the optimum for each channel considered separately. One important reason for this is that the result of Murray’s single-channel optimization (Q∝R3) is independent of the pressure drop; so when applied to a branching network, the relative pressure drops over the daughter channels are not considered, and the optimization is under-constrained. Murray’s original principle, over time, has been misinterpreted as a general branching law (for symmetric and asymmetric configurations), leading to the prevalence of the incorrect form shown in equation (1.1). This misinterpretation has endured in subsequent literature regarding turbulent flow [21] and non-Newtonian fluids [17,18], as we shall now demonstrate.

Turbulent flow

Turbulent flow is described by the phenomenological Darcy–Weisbach equation, which relates the pressure drop to the mean velocity for an incompressible fluid in a channel of arbitrary cross-sectional shape where f is the Darcy friction factor, is the mean streamwise velocity, is the hydraulic diameter, and is the wetted perimeter. Note that equation (2.25) is applicable to gravity-driven open channels, e.g. rivers, as well as closed pipes. In a river, the pressure drop is a function of the channel slope [23]. Making the substitutions (which is a property of the cross-sectional shape, like S) and , equation (2.25) can be rewritten as and the flow resistance per unit length k is For turbulent flows, the pressure drop is proportional to the square of the volumetric flow rate, so k is a function of Q. However, for all constraint options of the optimization described by equation (2.1), dQ/dΓ=0. For the sake of deriving an analytical expression comparable to Murray’s law, we restrict our interest to fully rough-wall turbulent flow, where the friction factor is also approximately constant[1] —i.e. it is independent of the Reynolds number and the volumetric flow rate. In this regime, the main applications are civil engineering hydraulic and pneumatic systems. So, when the shape is constant through the network, substituting equation (2.27) into the generalized law (equation (2.10)) gives Combining equations (2.20) and (2.27) produces the cross-sectional area relationship between the ith and jth daughter channels Substituting equation (2.29) into equation (2.28) and rearranging for the daughter–parent area ratio gives Equation (2.30) relates the area of the parent channel to the area of the jth daughter channel in an optimized two-level network for turbulent flow. It is valid for channels of any cross-sectional shape, provided the shape is constant through the network, and is only equivalent to the turbulent Murray’s law [21] (equation (1.2)) for symmetric branching, i.e. Ψ=Ψ=1/N and Φ=1, where (2.30) reduces to This also agrees with the turbulent flow symmetric branching results from previous studies [22,23]. Comparing equations (2.24) and (2.31) shows that, for symmetric branching, the optimal daughter–parent area ratio is smaller for turbulent flow than it is for laminar flow.

Non-Newtonian fluid flow

Non-Newtonian fluids are typically characterized by a nonlinear relationship between shear stress and shear-strain rate. The power-law constitutive model [19,20] is one of the most popular, enabling a wide range of engineering problems to be solved analytically. For fluid flow through a circular channel, this is where m is the flow consistency index, dd/dr is the shear strain rate and n is the flow behaviour index. This relationship leads to an effective viscosity of Power-law fluids can be divided into three classes based on their flow behaviour index: (i) pseudo-plastic (shear thinning) fluids (n<1) exhibit a decrease in viscosity with increased shear strain rate; (ii) Newtonian fluids (n=1) exhibit a constant viscosity; and (iii) dilatant (shear thickening) fluids (n>1) exhibit an increase in viscosity with increased shear strain rate. Owing to the difficulty of applying the power-law model in two dimensions, the optimal branching of non-Newtonian fluids are considered here only for circular cross sections.

In cylindrical coordinates, the steady-state Navier–Stokes equation for incompressible laminar flow through a circular cross section is where u is the streamwise velocity (the radial and swirl velocity components are assumed to be zero). Substituting equation (2.33) into (2.34) and integrating with respect to r gives where C1 is a constant. At the midpoint of the cross section, when r=0, the velocity is at a maximum and thus dd/dr=0; therefore, C1=0. Integrating with respect to r once more produces where C2 is another constant. Inserting the no-slip condition at the wall into equation (2.36) produces the non-Newtonian velocity profile The volumetric flow rate is obtained by integrating the fluid momentum over the cross-sectional area, which, in cylindrical coordinates, is where θ is the azimuth. Substituting equation (2.37) into (2.38) gives By setting n=1 in equation (2.39) the Hagen–Poiseuille equation is recovered. Noting that A=πR2, the flow resistance per unit length is For non-Newtonian fluids, there is a nonlinear relationship between the pressure gradient and the volumetric flow rate, so k is again a function of Q. As explained previously, dQ/dΓ=0 for all the constraint options in the optimization described by equation (2.1), so substituting equation (2.40) into the asymmetric generalized law (equation (2.10)) gives Combining equations (2.20) and (2.40) produces the cross-sectional area relationship between the ith and jth daughter channels Substituting equation (2.42) into equation (2.41), and rearranging for the daughter–parent area ratio gives Equation (2.43) relates the area of the parent channel to the area of the jth daughter channel in an optimized two-level network of circular channels transporting a non-Newtonian fluid. By setting n=1, the asymmetric generalized law for Newtonian fluid flows (equation (2.23)) is retrieved. Equation (2.43) shows that Γ is dependent on the flow behaviour index n, contrary to results from previous studies on non-Newtonian branching flows that used Murray’s law [17,18]. The optimal daughter–parent area ratio is independent of n only for symmetric branching, i.e. Ψ=Ψ=1/N and Φ=1: This is exactly the same as equation (2.24) for symmetric Newtonian flows and agrees with previous studies [17,18] for symmetric non-Newtonian flows. Equation (2.43) can also be used to determine the optimal area ratio for the theoretical limits of the flow behaviour index n. For the shear-thickening-fluid limit, when , equation (2.43) reduces to Interestingly, equation (2.45) is independent of the daughter–daughter pressure-gradient ratio Φ and is exactly the same as Murray’s law for asymmetric branching (1.1) when posed in terms of areas. For the shear-thinning-fluid limit, when , equation (2.43) reduces to For this limit, the optimal area ratio is independent of the daughter flow-rate fraction Ψ, so when the daughter channel pressure gradients are equal, i.e. Φ=1, the daughter channel areas are also equal and Γ=N−2/3.

Numerical verification and discussion

In this section, we construct a numerical model of a two-level branching network which adopts the same fluid-physics assumptions used in Murray’s original paper, its subsequent extensions for turbulent flows and non-Newtonian fluid flows, and our own generalized law. These are (i) the flow through each channel is steady-state, incompressible and fully developed; (ii) the pressure is continuous throughout the network; and (iii) the pressure linearly varies over the entire length of each channel from inlet/outlet to a common branching point. The purpose of using the same physical model for the fluidic network is solely to demonstrate that its optimization does not lead to Murray’s law, but to the generalized law we derived above; verification of Murray’s fluid-physics assumptions is beyond the scope of this paper.

For clarity, we present the numerical model in a form specific to laminar flow through cylindrical channels (as per Murray’s original case), and refer the reader to appendix A for a more general description. The flow through each channel of the fluidic network is determined by momentum conservation, and is treated as being positive if it flows towards the point of branching—i.e. flow through the parent channel will be positive and flow through daughters will be negative. In the specific case of laminar flow through a cylindrical channel, and given the previously stated fluid-physics assumptions, this is the Hagen–Poiseuille law where q is the volumetric flow rate through the ith channel in the network, a is the cross-sectional area of the ith channel, p is the pressure at the end of the ith channel (i.e. the inlet of the parent channel or the outlet of a daughter channel), l is the length of the ith channel and pB is the pressure at the point of branching (which is common to all channels). Note, here, unlike in our analytical derivation, the subscript i could denote either the parent channel (i=1) or one of the daughter channels (i=2,3…,M, where M is the total number channels that comprise the network). The model is completed by mass continuity at the branching point, i.e. calculation of the total network volume v, i.e. and the definition of the cross-sectional area ratio between the (i+1)th and ith channels

The system of equations (3.1)–(3.4) can be solved for the mass flow rates q(1: if the following are fixed: pressure at boundaries p(1:, channel lengths l(1:, fluid viscosity μ, network volume v, and the channel cross-sectional area ratios γ(1:. In this paper, the solution is obtained using the trust-region dogleg algorithm [34] in MATLAB.

To enable a comparison with Murray’s law and our generalized law, we now optimize the network model (i.e. equations (3.1)–(3.4)) using a brute-force approach. For otherwise fixed properties (e.g. fixed volume, boundary pressures, etc.), the cross-sectional area ratio between the first daughter channel and the parent channel area γ1 is varied, through all physically viable values, to locate that which maximizes the volumetric flow rate through the network. This result can then be compared directly with Murray’s law and the generalized law, as the definition of γ1 is equivalent to that of Γd.

The first set of optimization results demonstrate that Murray’s law is suboptimal for asymmetrically branching networks of any cross-sectional shape, even for laminar flow. To demonstrate that our generalized law is valid for any cross-sectional shape (as long as the shape remains constant through the network), the numerical verification is performed for three different cross sections: circular, square, and rectangular with an aspect ratio α=5. For circles, S=1/(8π), and for rectangles, an accurate approximation of S is calculated from simulations, using a standard central-difference solution of the laminar Navier–Stokes equations (2.12).[2] For Murray’s law, mass conservation provides the closure , which leads to

In figure 1, to induce asymmetry, the daughter flow-rate fraction Ψ is varied whereas the daughter–daughter pressure-gradient ratio is kept constant at Φ=1. All solutions that the greater the fraction of flow through the daughter channel, the greater the optimum daughter’s area (relative to the parent); as is intuitive. The results confirm the finding that Murray’s law is only optimal for symmetric bifurcations (Ψ=0.5); for a flow-rate percentage of 10% (Ψ=0.1), Murray’s law under-predicts the optimum daughter area by as much as 26%. In contrast, the generalized law is accurate for all values of Ψ for all cross-sectional shapes tested. This confirms the analytical finding that Murray’s law has been mistakenly applied to asymmetrically branching networks, where the optimized result for each individual channel is not optimal for the network as a whole.

Figure 1.

Optimal daughter–parent area ratio (Γ) against daughter flow-rate fraction (Ψ) for laminar flow through a two-level bifurcating network with equal daughter channel pressure gradients (Φ=1). Plotted for Murray’s law (equation (3.5)), the generalized law (equation (2.23)), and the results from the numerical optimization for circles, squares, and rectangles of aspect ratio α=5. (Online version in colour.)

This can be further demonstrated by inducing asymmetry by varying the daughter–daughter pressure-gradient ratio Φ, and maintaining a constant daughter flow-rate fraction Ψ=0.5, as shown in figure 2. Murray’s law does not consider Φ to be a variable that affects the optimal daughter–parent area ratio Γ and shows a notable departure from the numerical optimization results; for a pressure-gradient ratio of Φ=2, Murray’s law over predicts the optimum daughter area by 24%. In contrast, the generalized law is accurate for all values of Φ, for all cross-sectional shapes and, as expected, shows that the optimal area of the jth daughter channel is smaller when it has a larger pressure gradient relative to the other daughter channel, as the flow rate flowing through each daughter is equal (Ψ=Ψ=0.5). This result is the same whether the pressure gradient is altered by varying the relative daughter channel lengths or the pressure drops. Figures 1 and 2 both show that as the extent of asymmetry increases, Murray’s law provides a poorer estimate of the optimal area ratio.

Figure 2.

Optimal daughter–parent area ratio (Γ) against daughter–daughter pressure-gradient ratio Φ=(ΔPd/Ld)/(ΔPd/Ld) for laminar flow through a two-level bifurcating network with equal flow through each daughter channel, i.e. Ψ=Ψ=0.5. Plotted for Murray’s law (equation (3.5)), the generalized law (equation (2.23)), and the results from the numerical optimization for circles, squares, and rectangles of aspect ratio α=5. (Online version in colour.)

The next set of optimization results are for fully rough-wall turbulent flow through an asymmetrically bifurcating network of channels with arbitrary, but constant, cross-sectional shape. For Murray’s law, the closure is provided by mass conservation, which leads to

The optimization results shown in figures 3 and 4 have asymmetry induced by varying the daughter flow-rate fraction Ψ and daughter–daughter pressure-gradient ratio Φ, respectively. Again, there is excellent agreement between the numerical optimization and the turbulent generalized law for all values of Ψ and Φ. Both figures 3 and 4 show that, except in the case of large asymmetries, the optimal daughter–parent area ratio for laminar flow is larger than the optimal area ratio for turbulent flow. This trend can broadly be explained by considering the equations for volumetric flow rate for laminar and turbulent flow (equations (2.17) and (2.26), respectively). For laminar flow Q∝A2, whereas for turbulent flow Q∝A5/4. Considering these relationships for the parent and the jth daughter channel, then and for laminar and turbulent flows, respectively (this is confirmed by the generalized laws for laminar flow (2.23) and turbulent flow (2.30)). As Ψ is always less than one, and thus the optimal area ratio for laminar flow will generally be larger than optimal area ratio for turbulent flow for the same fixed parameters.

Figure 3.

Optimal daughter–parent area ratio (Γ) against daughter flow-rate fraction (Ψ) for fully rough turbulent flow through a two-level bifurcating network of arbitrary, but constant, cross-sectional shape with equal daughter channel pressure gradients (Φ=1). Plotted for Murray’s law (equation (3.6)), the laminar generalized law (equation (2.23)), the turbulent generalized law (equation (2.30)), and results from the numerical optimization. (Online version in colour.)

Figure 4.

Optimal daughter–parent area ratio (Γ) against daughter–daughter pressure-gradient ratio (Φ) for fully rough-wall turbulent flow through a two-level bifurcating network of arbitrary, but constant, cross-sectional shape. The volumetric flow rate through each daughter channel is equal, i.e. Ψ= 0.5. Plotted for Murray’s law (equation (3.6)), the laminar generalized law (equation (2.23)), the turbulent generalized law (equation (2.30)), and results from the numerical optimization. (Online version in colour.)

Murray’s law proves to be a more accurate approximation for asymmetrically branching turbulent flows (compared with laminar flows), but still errs by 10% when Ψ=0.1 (and Φ=1), and by 19% when Φ=2 (and Ψ=0.5). For symmetric branching (Ψ=0.5 and Φ=1), the numerical and analytical solutions both agree with Murray’s law [21] and the results by [22,23].

The final set of optimization results are for non-Newtonian fluid flow through an asymmetrically bifurcating network of circular channels. For Murray’s law, equation (3.5) is used. In figure 5, asymmetry is induced by varying Ψ, whereas Φ=1 is constant. The results demonstrate that the optimal daughter–parent area ratio Γ is dependent on the flow behaviour index n, contrary to the results of previous studies based on Murray’s law [17,18].

Figure 5.

Optimal daughter–parent area ratio (Γ) against daughter flow-rate fraction (Ψ) in a two-level bifurcating network of circular channels with equal daughter channel pressure gradients (Φ=1). Comparison of Murray’s law (equation (3.5)), the generalized law (equation (2.43)), and results from the numerical optimization. Plotted for n=10−4, 0.1, 0.74 (blood), 1 (Newtonian fluid), 2 and 100. (Online version in colour.)

The Newtonian fluid case (n=1) is highlighted with a filled marker, and it is noted that this solution is the same as that shown in figure 1. It is observed that, for all n, the gradient (dΓ/dΨ) increases monotonically with increasing n. As , the fluid approaches the shear-thickening-fluid limit (equation (2.45)) where Murray’s law is correct for all Ψ. For smaller values of n, Murray’s law is correct only for symmetric bifurcations (Ψ=0.5) and becomes increasingly inaccurate as n decreases; for a flow-rate percentage of 10% (Ψ=0.1), Murray’s law under predicts Γ by 66% for n=10−4. The increasing error in the Murray’s law solution as n decreases is also shown in figure 6, where asymmetry is induced by varying Φ and Ψ=0.5 is fixed. Here, for a pressure-gradient ratio of Φ=2, Murray’s law over predicts the optimum daughter area by 172% when n=10−4. In contrast, the generalized law is accurate for all values of Ψ and n. The plot for n=0.74 is an approximation of the optimal area ratio for the cardiovascular system, based on the measurements of a falling-ball viscometer [35]. As , the fluid approaches the shear-thinning-fluid limit (equation (2.46)) and Γ becomes independent of Ψ.

Figure 6.

Optimal daughter–parent area ratio (Γ) against daughter–daughter pressure-gradient ratio (Φ) in a two-level bifurcating network of circular channels with Ψ=0.5. Comparison of Murray’s law (equation (3.5)), the generalized law (equation (2.43)), and results from the numerical optimization. Plotted for n=10−4, 0.1, 0.74 (blood), 1 (Newtonian fluid), 2, and 100. (Online version in colour.)

The reason for the shear-thickening- and shear-thinning-fluid limits can be found by examining the volumetric flow rate of a non-Newtonian fluid. When , by raising all terms to the power of n, equation (2.39) becomes 1=ΔPR/(2Lm) and the area ratio is only a function of the pressure gradient; hence, Γ does not vary with Ψ. When the daughter–daughter pressure-gradient ratio Φ decreases, the area must increase (and vice versa), as shown in figure 6. Conversely, when , equation (2.39) becomes Q=πR3/3 and the optimal area ratio is only a function of the volumetric flow rate; hence, Γ does not vary with Φ. This expression, with Q∝R3, is equivalent to Murray’s law (equation (1.1)), where the local optimization is the same as the global optimization, as shown in figures 5 and 6.

Conclusion

We have derived a generalized optimization principle that leads to analytical expressions for the optimum daughter–parent area ratio Γ for asymmetrically branching networks of any cross-sectional shape and for a range of fluidic models. This new optimal relation will enable deeper understanding of biological network behaviour and provide a generalized biomimetic design principle that can be applied to a variety of artificial branching systems to maximize their efficiency.

We have verified analytical solutions using a numerical optimization procedure and shown that, for symmetric branching of laminar and Newtonian fluids, our generalized law is equivalent to Murray’s law. However, when applied to an asymmetrically branching network, Murray’s law is suboptimal, as the global optimization of the entire network is not equal to the local optimization of each individual channel, which Murray’s law presumes. We further demonstrate that this mistake in the application of Murray’s law to asymmetric branching networks has endured for non-Newtonian fluids (e.g. in the cardiovascular system) and turbulent flows (e.g. in hydraulic or pneumatic civil engineering applications).

In non-Newtonian fluidic networks, Γ is dependent on the flow behaviour index n for asymmetric branching, contrary to what previous studies based on Murray’s law have stated. Murray’s law is only retrieved for non-Newtonian fluid networks at the shear-thickening limit, when and Γ is no longer dependent on the relative pressure gradients over the daughter channels. At the shear-thinning limit, when , Γ becomes independent of the relative flow rates through each daughter channel.

Table 1.

Volumetric flow rate constants for laminar flow, fully rough-wall turbulent flow, and non-Newtonian flow.

flow model	b1	b2	b3
laminar	Sμ_A_A	2	1
turbulent	8RfA_AA^A^A	54	12
non-Newtonian	n(3n+1)12mπ^(n+1)/2nA_A_AA^A^A	(3n+1)2n	1n