Engineering of branched fluidic networks that minimise energy dissipation

Abstract

Power minimisation of fluid transport in branched fluidic networks has become of paramount importance for microfluidics, additive manufacturing and hierarchical functional materials. For fully developed laminar flow of Newtonian fluids, Murray's theory provides a solution for the channel and network dimensions that minimise power consumption. However, design and optimisation of networks that transport complex fluids is still challenging. Here, we generalise Murray's theory towards fluid rheologies, including non-Newtonian (power-law) and yield-stress fluids (Bingham, Herschel–Bulkley, Casson). A straightforward graphical approach is presented that provides the optimal radii in a branching network, and the angles between these branches. The wall shear stress is found to be uniform over the entire network, and the velocity profile is self-similar. Furthermore, the effect of non-optimal channel radii on the power consumption of the network is investigated. Finally, examples illustrate how this approach applies to a wide variety of systems.

1. Introduction

Fluidic networks are ubiquitous in biological systems and engineering applications, as shown in figure 1(a). In biology, these networks are found in e.g. vascular networks (Murray 1926b; Hutchins, Miner & Boitnott 1976; Lee & Lee 2010), bronchial trees of the lungs (Hooper 1977; Xu et al. 2016a), and leaf veins of plants (McCulloh, Sperry & Adler 2003; Gleason et al. 2018). Here, evolution has resulted in fluidic networks that efficiently transport large amounts of heat or mass, with low dissipation during fluid transport and limited volumes (Bejan & Lorente 2013). These characteristics are also desirable for three recent fluidic-based platforms: microfluidics (Whitesides 2006), including microreactors (Dong et al. 2021), additive manufacturing (Visser et al. 2019), and functional ‘Murray’ materials in which heat or mass is diffused through the channel walls to the surrounding material (Sciubba 2016; Zheng et al. 2017). However, the per-channel throughput of these platforms is intrinsically low, for example because long reaction times are required in microreactors (Su et al. 2016; Madhawan et al. 2018). The per-channel flow rate of microfluidic systems can be increased by 1–3 orders of magnitude by entering the inertial regime (Di Carlo 2009) or by in-air microfluidics (Visser et al. 2018), but even this is insufficient for large-scale use. Paralellisation of these technologies is therefore required to harvest their full benefit (Holtze 2013; Skylar-Scott et al. 2019; Dong et al. 2022). For example, in microfluidics, Jeong et al. (2015) achieved scaling of channels up to 1000$\times$ to realise commercial throughput of particle-generating flows. In additive manufacturing of complex soft-robotic walkers, Skylar-Scott et al. (2019) achieved direct-write via 128 parallel channels by repeated bifurcation from the main channel for the ‘ink’ supply. In materials science, reactors for heat or mass exchange were parallelised both top-down by Dong et al. (2021) (resembling microfluidic chips) and bottom-up by self-assembly by Zheng et al. (2017) (resembling biological materials). Upscaling by paralellisation is typically achieved by distributing flow to individual channels via branched fluidic networks. Figure 1(b) shows that several levels of fluidic branching were achieved for the aforementioned applications, but also that biological systems still exhibit a far larger number of branchings. Therefore, further optimisation of engineered networks may help in realising the enormous application potential of microfluidics, additive manufacturing, and Murray materials. Similarly, optimised fluidic networks would benefit upscaling of emerging methods for catalysis (Yi et al. 2022), carbon capture (Moore et al. 2019), droplet and particle production (Holtze 2013; Yadavali, Jeong & Lee 2018), three-dimensional (3-D) printing of polymer foams (Visser et al. 2019), artificial arterial networks (Kinstlinger & Miller 2016), and microfluidic blood oxygenation (Dabaghi et al. 2020). The purpose of this paper is to provide optimisation tools that are suitable for a broad range of fluid rheologies and readily applicable by scientists and engineers who design fluidic networks.

Figure 1. (a) Examples of branched fluidic networks in engineering applications and biological systems. (b) Radius of the channel $R_i$ in the network as a function of the relative flow rate $Q_i/Q_0$ at level $i$ for the networks corresponding to the examples in (a). Every marker indicates a level in the branching. The error bars indicate the deviation from the marker (Singhal, Henderson & Horsfield 1973; Gan et al. 1993; Kassab et al. 1993; Nordsletten et al. 2006; Su et al. 2016; Carvalho et al. 2017; Zheng et al. 2017; Zhao et al. 2018; Skylar-Scott et al. 2019). Adapted with permission from Skylar-Scott et al. (2019), copyright (2019) The Authors, published by Springer Nature; and adapted under terms of the CC-BY licence from Zheng et al. (2017), copyright (2017) The Authors, published by Springer Nature.

Theoretical analyses of fluidic networks, including the present work, are usually constructed in two steps (Murray 1926b).

First, an expression is derived for the radius of a channel such that power consumption is minimised for a given flow rate. The key idea is to consider not only the power needed to maintain the flow, i.e. the product of pressure drop and flow rate, but also to minimise the channel's volume by considering the costs introduced by the volume in the channels. For that, a cost factor $\alpha$, representing the power per volume needed to maintain the fluid, is multiplied with the volume and included in the power function. For example, Murray (1926b) proposed the use of the metabolic energy to maintain blood in the body as a cost factor. For more details on the cost factor, see § 2.4.

Assuming a laminar flow of an incompressible Newtonian fluid in an artery of circular cross-section, the cube of the radius $R$ is proportional to the flow rate $Q$ for a power-optimised channel:

(1.1)\begin{equation} \frac{R^3}{Q}= \text{const.}\end{equation}

Second, the network is optimised by considering a bifurcation consisting of a parent channel (index ‘0’) and two daughter channels (indices ‘1’ and ‘2’) (Murray 1926a), as shown schematically in figures 2(a)–2(c). Employing mass conservation and assuming incompressibility ($Q_0=Q_1+Q_2$), the radii of the tubes must satisfy

(1.2)\begin{equation} R_0^3 = R_1^3 + R_2^3, \end{equation}

if (1.1) is satisfied in all channels.

Figure 2. (a) Branching of parent channel into $N$ daughter channels. The location of the branching point $\boldsymbol {x}$ follows from the analysis and determines the lengths $L_i$ of the channels. The grey channels indicate that it is possible to have many channels that originate from the branching point. (b) Schematic of a single branch. (c) Schematic of velocity profile of fully developed laminar flow of non-yield-stress fluid in pipe. (d) Schematic of velocity profile of fully developed laminar flow of yield-stress fluid in pipe. (e) The fluid models as analysed in this work. The colours are consistent across (c–f). ( f) Different fluid models characterised by $n$ and $\phi$. For the characterisation of Casson-like fluids, see Appendix A.

Murray's theory for this flow regime has been analysed further by e.g. Murray (1926a,b), Zamir (1977), Sherman (1981), Kamiya, Togawa & Yamamota (1974), Stephenson et al. (2015), Oka & Nakai (1987) and Miguel & Rocha (2018), and verified experimentally by e.g. Rossitti & Löfgren (1993) and Hutchins et al. (1976). In addition, it has been extended to high turbulent flow of Newtonian fluids by Uylings (1977), Kou et al. (2014) and Stephenson & Lockerby (2016), and to laminar flow of non-Newtonian (power-law) fluids (Mayrovitz 1987; Revellin et al. 2009; Miguel & Rocha 2018; Stephenson & Lockerby 2016). Flow regimes of more complex fluid models – Bingham (Reiner 1926), Herschel–Bulkley (Herschel & Bulkley 1926) and Casson (Venkatesan et al. 2013) – have been investigated for single pipe flow (Chhabra & Richardson 2008; Chilton & Stainsby 1998) (figure 2d). Ponalagusamy (2012) explored optimal branching of fluidic networks containing yield-stress fluids (Bingham and Herschel–Bulkley). Next to channels with circular cross-sections, channels with elliptical (Tesch 2010) and rectangular (Emerson & Barber 2012) cross-sections have been introduced, from which it appeared that the cross-section does not affect Murray's proportionality, in contrast to the flow regime (i.e. laminar or turbulent), which can change Murray's cube rule. Other effects, such as the effect of curved channels (Miguel 2018), pulsating flow (Painter, Edén & Bengtsson 2006), asymmetric branching (Zamir 1978) and the efficiency difference between bifurcation and trifurcations (Rosenberg 2020), have been investigated. The state of the art is discussed in more detail in Revellin et al. (2009), Bejan & Lorente (2006), Miguel (2016) and Xu et al. (2016b).

Next to the radii of the channels, the location of the branching point $\boldsymbol {x}$ is an important variable in the optimisation of the network, as it determines the lengths of the channels for a network with given end points $\boldsymbol {x}_i$, $i=0,1,\ldots,N$ (see figure 2a). Mostly, when theoretically optimising a fluidic network, this is done by Murray's proposal to minimise both the power of driving the fluid (minimisation of resistance) and the power to maintain the fluid (minimisation of volume) (Murray 1926a,b; Uylings 1977). Consequently, $\boldsymbol {x}$ has to be chosen such that the volume of the channels is minimal. The angles formed between the parent and daughter channels have gained special attention, and were investigated for the case of a bifurcation or a symmetric trifurcation, theoretically (Murray 1926a; Horsfield & Cumming 1967; Kamiya & Togawa 1972; Uylings 1977; Zamir 1976, 1977) and experimentally (Hutchins et al. 1976; Horsfield & Cumming 1967). These angles depend on the channel radii and the type of minimisation (Zamir 1976), and on other constraints such as fixed pressures at end points (Kamiya & Togawa 1972). Although the theoretical optimal radii of the channels match well with the observed results in biology, the correlation between theoretical optimal and observed angles between the channels is weaker (Hutchins et al. 1976), possibly related to additional geometrical constraints beyond network optimisation (Horsfield & Cumming 1967; Hutchins et al. 1976).

Despite all of the previous results, engineering branched fluidic networks remains challenging. First, many sources aim to describe biological fluidic networks, but biological networks are not only designed for energy minimisation (Koçillari 2021) and current methods are less suitable for engineering purposes. Second, complex mathematical calculations are often required to optimise networks that transport non-Newtonian fluids (see e.g. Ponalagusamy 2012). Third, alternative fluid models such as Casson have not been applied in branched fluidic networks before, despite their applicability to e.g. food processing and blood.

Therefore, here we introduce a procedure for calculating the geometry of an optimal branching for fully developed laminar flow of simple and complex fluid models (figure 2e), which is described in § 2. The losses of non-optimised networks, the wall shear stress, the self-similar velocity profile, and approaches to determine the cost factor are discussed. Section 3 shows example cases of the optimisation method, after which the validity of the theory is discussed in § 4, and the conclusions are presented in § 5.

2. Optimisation of a branching using Murray's theory

We consider a branching consisting of a parent channel connected to $N$ daughter channels in a branching point $\boldsymbol {x}$; see figure 2(a). The channels are numbered from $0$ to $N$, with $0$ indicating the parent channel. In this study, we choose to fix the flow rates $Q_i$, $i=0,1,\ldots,N$, and the position of the termination points of the channels $\boldsymbol {x}_i$, $i=0,1,\ldots,N$. These parameters being fixed, we aim at minimising the power required to perform this flow, with respect to the effective radii $R_i$, $i=0,1,\ldots,N$, of the channels and to the location of the branching point $\boldsymbol {x}$. Furthermore, $Q_0$ is taken positive towards the branching point, whereas the flow rates in the daughter channels $Q_i$, $i=1,\ldots,N$, are taken positive away from the branching point. To satisfy mass conservation (assuming incompressibility), the flow rates satisfy

(2.1)\begin{equation} Q_0=\sum_{i=1}^N Q_i. \end{equation}

The effective radius is defined as $R\equiv D_H/2$, where $D_H \equiv 4A/p$ is the hydraulic diameter, with $A$ being the cross-sectional area of the channel, and $p$ being the wetted perimeter of the cross-section. This can be used for deriving the flow rate through elliptical and rectangular cross-sections, and comparing channel characteristics, but will result in different expressions for the flow rate $Q$ than for circular cross-sections. In the present study, the focus is on only circular cross-sections, and here, the effective radius reduces to the channel radius. For channel flow with elliptical and rectangular cross-sections, the reader is referred to e.g. Tesch (2010), Emerson & Barber (2012), Lekner (2007) and Cornish (1928).

Furthermore, the lengths of the channels $L_i$ are functions of the branching location:

(2.2)\begin{equation} L_i\equiv|\boldsymbol{x}_i-\boldsymbol{x}|, \quad i=0,1,\ldots,N. \end{equation}

The power depends on the radii and lengths of the channels, and is the sum of the individual channel contributions:

(2.3)\begin{equation} P(\boldsymbol{R},\boldsymbol{x}) \equiv \sum_{i=0}^N \left\{\left|\frac{{\rm d}p}{{\rm d}z}\right| QL + \alpha V\right\}_i. \end{equation}

Here, $V$ is the channel volume, with

(2.4)\begin{equation} V_i = {\rm \pi}R_i^2 L_i, \end{equation}

the pressure gradient $|{\rm d}p/{\rm d}z|$ is a fluid-type dependent function derived in § 2.1, and $z$ is the coordinate in the flow direction of a channel. The pressure at the nodes is not defined or constrained. The cost factor is $\alpha$, representing the power per unit volume required to maintain the fluid ($\alpha \approx 1\ {\rm kW}\ {\rm m}^{-3}$; Murray 1926b).

Optimisation of the branched network is carried out by differentiation of (2.3) to $\boldsymbol {R}$ and $\boldsymbol {x}$, and equating these expressions to 0. From this, one obtains optimisation problems for: (1) the channel radii $\boldsymbol {R}$ (§ 2.2), and (2) the location of the branching point $\boldsymbol {x}$ (§ 2.3). The flow profile, dissipation characteristics and wall shear stress of the branched network are discussed in § 2.4.

2.1. Fully developed laminar channel flow

The pressure difference needed to drive a fluid through a channel is governed by the ‘fluid model’ that describes the rheological behaviour of the fluid. The (effective) dynamic viscosity $\mu$ describes the resistance of the fluid against shear, which can be shear-rate $\dot {\gamma }$ dependent according to $\mu = \mu '\,|\dot {\gamma }|^{n-1}$. Here, $\mu '$ is the flow consistency index, and $n>0$ is the flow index, with $n>1$ representing shear-thickening, and $0< n<1$ representing shear-thinning behaviour of the fluid. Furthermore, a fluid may have a yield stress $\tau _0 \geqslant 0$, which means that the fluid shears only if the local shear stress $\tau _{rz}$ exceeds the yield stress. The shear rate of the fluid as a function of the applied local shear stress for a Herschel–Bulkley fluid is described as

(2.5)\begin{equation} \dot{\gamma}(\tau_{rz}) = \begin{cases} \text{sign}(\tau_{rz})\left(\dfrac{|\tau_{rz}|-\tau_0}{\mu'} \right)^{1/n} & \text{if $|\tau_{rz}|\geqslant\tau_0$},\\[8pt] 0 & \text{if $|\tau_{rz}|<\tau_0$}. \end{cases} \end{equation}

Here, $\dot {\gamma }\equiv {\partial u}/{\partial r}\leqslant 0$, with $u$ the axial velocity. The velocity in the radial and azimuthal directions is assumed to be 0. In the following, we assume that fluid properties $\mu '$, $\tau _0$ and $n$ are constant, and that there is fully developed laminar flow in cylindrical channels.

The local shear stress (with $R\ll L$) is then described as a function of the axial pressure gradient ${\rm d}p/{\rm d}z$:

(2.6)\begin{equation} \tau_{rz} ={-}\frac{1}{2}\,\frac{{\rm d}p}{{\rm d}z}\,r, \end{equation}

with $r$ the radial coordinate. Let $R$ be the radius of the channel, and let $R_p$ be the plug radius in case of yield-stress fluids:

(2.7)\begin{equation} R_p \equiv \frac{2\tau_0}{\left|\dfrac{{\rm d}{p}}{{\rm d}z}\right|}. \end{equation}

A schematic image of the plug radius $R_p$ is presented in figure 2(d). Scaling this plug radius with the radius of the tube results in the dimensionless plug radius $\phi$:

(2.8)\begin{equation} \phi \equiv \frac{R_p}{R} = \frac{2\tau_0}{\left|\dfrac{{\rm d}{p}}{{\rm d}z}\right|R}. \end{equation}

Using these definitions, (2.5) can be rewritten in terms of $\phi$ and $n$:

(2.9)\begin{equation} \dot{\gamma}(r) = \begin{cases} \text{sign}(\tau_{rz})\left( \left(\dfrac{R}{2\mu'}\left|\dfrac{{\rm d}{p}}{{\rm d}z}\right|\right) \Bigg(\dfrac{r}{R}-\phi\Bigg)\right)^{1/n} & \text{if $r\geqslant R_p$},\\[12pt] 0 & \text{if $r< R_p$}. \end{cases} \end{equation}

Figure 2( f) gives an overview of the different values assigned to $\phi$ and $n$ for each fluid model. For Newtonian fluids, (2.9) collapses to the well-known equation $\dot {\gamma }(r) = \text {sign}(\tau_{rz})\,({r}/{2\mu })\,|{\rm d}p/{\rm d}z|$. For a derivation of Casson-like fluids (such as blood), see Appendix A.

Integration of $\dot {\gamma }$ and applying a no-slip boundary condition at the wall ($u(r=R)=0$) results in the velocity profile

(2.10) \begin{align} u(r) = \begin{cases} \dfrac{nR}{n+1}\left(\dfrac{R}{2\mu'}\left|\dfrac{{\rm d}p}{{\rm d}z}\right|\right)^{1/n} \Biggl((1-\phi)^{{(1+n)}/{n}} - \Biggl(\dfrac{r}{R}-\phi\Biggr)^{({1+n})/{n}} \Biggr) & \text{if $r\geqslant R_p$},\\[12pt] \dfrac{nR}{n+1}\left(\dfrac{R}{2\mu'}\left|\dfrac{{\rm d}p}{{\rm d}z}\right|\right)^{1/n} (1-\phi)^{({1+n})/{n}} & \text{if $r< R_p$}. \end{cases} \end{align}

Here, $u$ is restricted to a velocity profile directed positively in the $z$ direction. The negative direction is simply $-u(r)$. The flow field is integrated over the cross-section to obtain the flow rate $Q$:

(2.11)\begin{equation} Q \equiv \iint_A u(r) \,\text{d}A = 2{\rm \pi} \int^R_0 u(r)\,r \,\text{d}r, \end{equation}

so $Q$ is a function of $|{\rm d}p/{\rm d}z|$, $R$, $\mu '$, $\tau _0$ and $n$. Dimension analysis shows that $Q$ can be written in the form

(2.12)\begin{equation} Q = {\rm \pi}R^3 \left(\frac{R}{2\mu'}\left|\frac{{\rm d}p}{{\rm d}z}\right|\right)^{1/n} \left(\frac{n}{3n+1}\right)\times \psi(\phi,n), \end{equation}

or, in terms of $\phi$ as

(2.13)\begin{equation} Q = {\rm \pi}R^3\left(\frac{\tau_0}{\mu'\phi}\right)^{1/n} \left(\frac{n}{3n+1}\right)\times \psi(\phi,n), \end{equation}

where we have included the factor ${{\rm \pi} n}/({3n+1})$ for convenience. Here, $\psi (\phi,n)$ is the dimensionless flow rate – the ratio of the flow rate $Q$ and the flow rate of a non-yield fluid $Q_{NY}$ – which is dependent solely on $\phi$ and $n$, and has a value between 0 and 1. The pressure gradient in relation to the flow rate $Q$ is obtained by rewriting (2.12):

(2.14)\begin{equation} \left|\frac{{\rm d}p}{{\rm d}z}\right| = \frac{2\mu'}{R}\left(\frac{Q}{{\rm \pi} R^3}\right)^n \left(\frac{3n+1}{n}\right)^n\frac{1}{\psi(\phi,n)^n}. \end{equation}

The general expression for $\psi (\phi,n)$ can be found by computing the integral in (2.11):

(2.15)\begin{equation} \psi = \frac{(1-\phi)^{(n+1)/n}}{(3n+1)^{{-}1}} \times\left( \frac{(1-\phi)^2}{3n+1}+\frac{2\phi(1-\phi)}{2n+1}+\frac{\phi^2}{n+1}\right). \end{equation}

For Bingham fluids ($n=1$), the expression for $\psi$ reduces to

(2.16)\begin{equation} \psi = 1-\tfrac{4}{3}\phi+\tfrac{1}{3}\phi^4 , \end{equation}

and for the case of non-yield-stress fluids ($\phi =0$), $\psi$ always reduces to 1.

Equation (2.13) shows that $Q/R^3$ is solely a function of the dimensionless plug radius $\phi$, and vice versa. The same holds for the term $|{\rm d}p/{\rm d}z|\,R$, which is obtained if $\phi$ is known. This property will be used in the characterisation of the optimised networks.

The Darcy–Weisbach equation is commonly used for calculation of the pressure drop for pipe flow, which is given as

(2.17)\begin{equation} \Delta p = f\,\frac{\rho}{4{\rm \pi}^2 }\,\frac{Q^2}{R^5}\,L, \end{equation}

in which the Darcy friction factor $f$ for the laminar flows of all treated fluid models is then provided by

(2.18)\begin{equation} f = \frac{64}{Re'\,\psi^n}, \end{equation}

where (Garcia & Steffe 1987)

(2.19)\begin{equation} Re' = 8\left(\frac{n}{3n+1}\right)^n\frac{\rho}{\mu'}\left(\frac{Q}{{\rm \pi} R^3}\right)^{2-n} R^2. \end{equation}

For optimised fluidic networks with laminar flow in the first branch ($R^3/Q = \text {const}.$, as discussed in the next subsection), $Re'$ scales with $R^2$, a generalisation of what Cohn (1955) derived for symmetric bifurcations. Therefore, all daughter channels will have a lower Reynolds number than the parent channel, ensuring laminar flow throughout the network if $Re'_0< Re'_{crit}$. The critical Reynolds number $Re'_{crit}$ is a function of $\phi$ and $n$ that was derived to be (Hanks & Ricks 1974):

(2.20)\begin{equation} Re'_{crit} =\frac{6464n}{(1+3n)^2}\times (2+n)^{({2+n})/({1+n})}\times\frac{\psi^{2-n}}{(1-\phi)^{({n+2})/{n}}}, \end{equation}

which reduces to $Re'_{crit} \approx 2100$ for Newtonian fluids.

2.2. Optimisation of the channel radii

Every channel within a branching contributes to the total power. When differentiating the total power to the individual radii $R_i$ for given flow rate $Q_i$, it is found that the optimisation condition for that radius $R_i$ depends only on the power contribution of the corresponding channel. Furthermore, the optimal radii $R_i$ do not depend on the lengths of the channels $L_i$. Therefore, calculation of the optimal channel radius is a decoupled problem, in which the power defined by (2.3) attains a global minimum with respect to the radius of each individual channel if

(2.21)\begin{equation} \frac{\partial P_i}{\partial R_i}=\frac{\partial}{\partial R_i}\left(\left|\frac{{\rm d}p}{{\rm d}z}\right|_i\right)Q_iL_i+2{\rm \pi}\alpha R_iL_i=0, \quad i=0,1,\ldots,N. \end{equation}

The optimisation condition for $R_i$ is obtained by differentiating (2.14) to $R_i$ (see § B.1), resulting in

(2.22)\begin{equation} \frac{R_{i,*}^3}{Q_i} = \frac{1}{\rm \pi}\left(\frac{\mu'}{\alpha}\left(\frac{3n+1}{n}\right)^n\frac{J_*}{\psi_*^n}\right)^{1/(n+1)}, \quad i=0,1,\ldots,N, \end{equation}

where functions evaluated at optimal radius are denoted by a subscript $*$, and $J$ is defined as

(2.23)\begin{equation} J\equiv 1+\frac{3n}{1-\dfrac{n\phi}{\psi}\,\dfrac{\partial \psi}{\partial\phi}}. \end{equation}

For a Herschel–Bulkley fluid, substitution of $\psi$ from (2.15) into (2.23) results in the expression

(2.24)\begin{equation} J = \frac{3n+1}{\dfrac{6n^3}{(2n+1)(n+1)}\,\phi^3+\dfrac{6n^2}{(2n+1)(n+1)}\,\phi^2+\dfrac{3n}{2n+1}\,\phi +1}. \end{equation}

For Bingham fluids ($n=1$), the expression for $J$ reduces further to

(2.25)\begin{equation} J = \frac{4}{\phi^3+\phi^2+\phi +1}, \end{equation}

and for non-yield fluids ($\phi =0$), it reduces to $J=3n+1$. Finally, for $\phi =1$, one obtains $J=1$.

For the case of Newtonian and power-law fluids ($\phi =0$), the optimisation condition for $R_*^3/Q$ is calculated relatively easily (with $J_*=3n+1$ and $\psi _*=1$). However, for yield-stress fluids, determining the optimal value for $R_*^3/Q$ is more difficult. Therefore, the following procedure is proposed.

First, the following dimensionless numbers are introduced:

(2.26)$$\begin{gather} \frac{\tilde{R}^3}{\tilde{Q}}\equiv \frac{R^3}{Q}\left(\frac{\alpha}{\mu'}\right)^{{1}/({n+1})}, \end{gather}$$

(2.27)$$\begin{gather}\tilde{\tau}_0 \equiv \frac{\tau_0}{(\mu'\alpha^n)^{{1}/({n+1})}}. \end{gather}$$

Using these dimensionless numbers, (2.22) is rewritten as

(2.28)\begin{equation} \frac{\tilde{R}_{i,*}^3}{\tilde{Q_i}} = \frac{1}{\rm \pi}\left(\left(\frac{3n+1}{n}\right)^n\frac{J_*}{\psi_*^n}\right)^{1/(n+1)}, \quad i=0,1,\ldots,N. \end{equation}

Applying the same non-dimensionalisation to (2.13), and rewriting it as a function of $\psi$, then substitution into (2.28), results in the alternative expression

(2.29)\begin{equation} \frac{\tilde{R}_{i,*}^3}{\tilde{Q_i}} = \frac{J_*\tilde{\tau}_{0}}{{\rm \pi}\phi_*}, \quad i=0,1,\ldots,N. \end{equation}

When eliminating $\tilde {R}_*^3/\tilde {Q}$ from (2.28) and (2.29), we obtain the expression

(2.30)\begin{equation} \tilde{\tau}_{0} = \phi_* \left(\frac{3n+1}{n\,J_*(\phi_*,n)\,\psi_*(\phi_*,n)}\right)^{{n}/({n+1})}. \end{equation}

Here, $\tilde {\tau }_0$ is composed of fluid and system properties, and is therefore considered as known for an optimisation process. As both $\psi$ and $J$ are solely functions of $\phi$ and $n$, the value of $\phi _*$ corresponding to $\tilde {\tau }_{0}$ can be found implicitly using (2.15) and (2.24). Figure 3 shows the calculated optimal values of $\phi _*$ as functions of $\tilde {\tau }_0$ and $n$. Using $\phi _*$, the corresponding optimisation condition for $\tilde {R}_{i,*}^3/\tilde {Q}_i$ is calculated readily using (2.28) or (2.29), and plotted in figure 4. The optimal channel diameter $R_{i,*}$ at desired flow rate $Q_i$ then follows directly from (2.26).

Figure 3. Contour plot of the dimensionless plug radius for the optimised network $\phi _*$ as a function of $\tilde {\tau }_{0}$ and $n$.

Figure 4. Contour plot of $\tilde {R}_*^3/\tilde {Q}$ as a function of $\tilde {\tau }_0$ and $n$. The colours represent the different fluid models and correspond to figure 2. The Herschel–Bulkley model covers the entire parameter space of $\tilde {\tau }_0$ and $n$.

For the non-yield limit ($\tilde {\tau }_0\rightarrow 0$), the optimisation condition becomes $\tilde {R}_{i,*}^3/\tilde{Q_i}\approx ({3n+1}){{\rm \pi} }^{-1}\, n^{(-n/(n+1))}$, which has a maximum $\tilde {R}_{i,*}^3/\tilde{Q_i}\approx 1.41$. This simplified expression is valid for $\tilde {\tau }_0<10^{-2}$, because then the fluid yield-stress effect in the optimised network becomes negligible. For larger values of $\tilde {\tau }_0$, plug formation becomes more relevant, and $\tilde {R}_{i,*}^3/\tilde{Q_i}$ becomes larger than in the case of a non-yield fluid. For the yield limit ($\tilde {\tau }_0\rightarrow \infty$), the optimisation condition tends to go to $\tilde {R}_{i,*}^3/\tilde{Q_i}\rightarrow \tilde {\tau }_0/{\rm \pi}$. Figure 4 also shows the parameter space of the different fluid models. The graphical approach of figure 4 for the optimal value of $R_{i,*}^3/Q_i$ for complex fluids prevents extensive calculations, enabling straightforward optimisation of the geometry of fluidic networks.

2.3. Optimisation of the network topology

When the optimal channel radii are calculated according to § 2.2, the second optimisation step is made to provide the location of the branching point $\boldsymbol {x}$. As shown in figure 2(a), the location of the branching point determines the length of all channels. Therefore, the optimised branching point $\boldsymbol {x}_*$ minimises the total (volume- and dissipation-induced) power of the network. Now the total power consumption of the network in (2.3) is differentiated to $\boldsymbol {x}$ and set to 0, resulting in

(2.31)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{x}}{P}=\sum_{i=0}^N\left(\left|\frac{{\rm d}p}{{\rm d}z}\right|_i Q_i + \alpha{\rm \pi} R_i^2\right) \boldsymbol{\nabla} L_i=0. \end{equation}

Substituting the optimal radii (2.22) results in the optimisation condition of the branching point $\boldsymbol {x}=\boldsymbol {x}_*$ (for a full derivation, see § B.2):

(2.32a,b)\begin{equation} \sum_{i=0}^N R_{i,*}^2 \boldsymbol{e}_{i,*}=\textbf{0}, \quad \boldsymbol{e}_{i,*}\equiv\left(\boldsymbol{\nabla} L_i\right)_*=\frac{\boldsymbol{x}_*-\boldsymbol{x}_i}{|\boldsymbol{x}_*-\boldsymbol{x}_i|}. \end{equation}

From this implicit equation, the coordinates of $\boldsymbol {x}_*$ are solved by simple numerical methods, providing $V_i$. Note that $\boldsymbol {x}_*$ can be located in a 3-D space, depending on the location of the end points $\boldsymbol {x}_i$. The resulting $\boldsymbol {x}_*$ determines a network that is optimised with respect to both the channel radii (2.22) and the channel lengths (2.32a,b).

The global minimum power of the branched network $P$ is obtained by inserting the result of (2.22) into (2.3):

(2.33)\begin{equation} P_*=\left(\frac{J_*+2}{J_*}\right)\alpha \sum_{i=0}^N V_{i,*}. \end{equation}

The angles between the channels are calculated by taking the inner product of (2.32a,b) with unit vectors that originate in the branching point and point towards the nodes. The resulting cosines of the corresponding angles between the two unit vectors can be calculated using $\cos (\theta _{ij})=\boldsymbol {e}_i\boldsymbol {\cdot }\boldsymbol {e}_j$:

(2.34)\begin{equation} \sum_{i=0}^N R_{i,*}^2 \boldsymbol{e}_{i,*}\boldsymbol{\cdot} \boldsymbol{e}_{j,*}=0. \end{equation}

This equation results in a linear system of $N+1$ equations for $\boldsymbol {e}_{i}\boldsymbol {\cdot } \boldsymbol {e}_{j}$ ($i,j=0,1,\ldots,N$, with $i\neq j$). For a bifurcation and a symmetric trifurcation, the angles between the channels for the optimal branching are independent of the coordinates of the nodes $\boldsymbol {x}_i$ (see figure 5). For both symmetric and asymmetric bifurcations, i.e. $N=2$, the optimal branching point $\boldsymbol {x}_*$ lies in a plane in the 3-D space spanned by $\boldsymbol {x}_0,\boldsymbol {x}_1,\boldsymbol {x}_2$, as shown in figure 5(a). The cosines of the smallest angles between each pair of channels involved are given by

(2.35)\begin{equation} \left.\begin{gathered} \boldsymbol{e}_0\boldsymbol{\cdot}\boldsymbol{e}_1=\cos(\theta_{01})={-}\frac{R_0^4+R_1^4-R_2^4}{2R_0^2R_1^2}, \\ \boldsymbol{e}_0\boldsymbol{\cdot}\boldsymbol{e}_2=\cos(\theta_{02})={-}\frac{R_0^4-R_1^4+R_2^4}{2R_0^2R_2^2}, \\ \boldsymbol{e}_1\boldsymbol{\cdot}\boldsymbol{e}_2=\cos(\theta_{12})=\frac{R_0^4-R_1^4-R_2^4}{2R_1^2R_2^2}. \end{gathered}\right\} \end{equation}

Figure 5. (a) Angles $\theta _{ij}$ in a bifurcation. (b) Angles $\theta _{ij}$ in a trifurcation

For a symmetric trifurcation, i.e. $N=3$, the optimal branching point $\boldsymbol {x}_*$ lies in a plane in the 3-D space spanned by $\boldsymbol {x}_0$, $\boldsymbol {x}_1$ or $\boldsymbol {x}_3$, and $\boldsymbol {x}_2$, where symmetry is assumed, i.e. channels 1 and 3 are mirrored in the line from $\boldsymbol {x}_0$ to $\boldsymbol {x}_2$ (figure 5b). This means that $R_1=R_3$, and the angles that they make with channels 0 and 2 are equal. In addition, $\boldsymbol {x}_0$, $\boldsymbol {x}_*$ and $\boldsymbol {x}_2$ are on the same line, reducing the unknown angles even further. The cosines of the smallest angles between each pair of channels involved are given by

(2.36)\begin{equation} \left.\begin{gathered} \boldsymbol{e}_0\boldsymbol{\cdot}\boldsymbol{e}_1=\frac{R_2^2-R_0^2}{2R_1^2}, \\ \boldsymbol{e}_1\boldsymbol{\cdot}\boldsymbol{e}_2=\frac{R_0^2-R_2^2}{2R_1^2}, \\ \boldsymbol{e}_1\boldsymbol{\cdot}\boldsymbol{e}_3=\frac{(R_0^2-R_2^2)^2-2R_1^4}{2R_1^4}. \end{gathered}\right\}\end{equation}

The results in (2.35) and (2.36) are in line with the findings of e.g. Zamir (1976).

For other cases, the coordinates of $\boldsymbol {x}_i$ are needed for calculating the angles between the channels because the system of equations from (2.34) is underdetermined. For a derivation of this general case as well as the above equations for bifurcations and trifurcations, see Appendix C.

2.4. The effect of non-optimised branching, the velocity profile, the wall shear stress, and the cost factor

Existing fluidic networks are not always optimised, as their dimensions may depend on physical constraints or industry standards. In this subsection, we compare the energy dissipation of such non-ideal networks to optimised fluidic networks as described in §§ 2.2 and 2.3.

For single channels within a network, the minimum power for a channel optimised with respect to the channel radius is obtained by taking the $i$th element of (2.3) and substituting (2.14) and the optimisation condition (2.22):

(2.37)\begin{equation} P_{i,*} = \frac{J_*+2}{J_*}\,\alpha {\rm \pi}R_{i,*}^2 L_i, \quad i=0,1,\ldots,N. \end{equation}

Next, the ratio of the non-optimised power $P_i$ and its minimum value $P_{i,*}$ per channel is found by dividing the $i$th element of (2.3) by (2.37) (see also § B.1):

(2.38)\begin{equation} \frac{P_i}{P_{i,*}} = \frac{2}{J_{*}+2}\left(\frac{\psi_{*}}{\psi_i}\right)^n \left(\frac{R_i}{R_{i,*}}\right)^{-(3n+1)} + \frac{J_{*}}{J_{*}+2}\left(\frac{R_i}{R_{i,*}}\right)^{2}, \quad i=0,1,\ldots,N. \end{equation}

Note that the optimal parameters $J_*$, $\phi _*$ and $\psi _*$ are without index $i$, because these are constant over the whole optimised branching (see also below). Substitution of the optimal values for a certain fluid and system ($J_*$ and $R_*$) results in an expression dependent on $\phi$ and $R$. As $\phi$ and $R$ are related by (2.13), the dimensionless power consumption is calculated by knowing the actual radii.

Substitution of $\psi _*=\psi =1$ and $J_*=4$ for laminar flow of a Newtonian fluid in (2.38) recovers the relations for power dissipation in a network as obtained by Murray (1926b) and Uylings (1977). The proof of this result is provided in § B.1.

The ratio of the actual power for the whole network $P$ and the optimal power $P_*$ is given by

(2.39)\begin{equation} \frac{P}{P_*} = \sum_{i=0}^N \frac{P_i}{P_*} = \sum_{i=0}^N \frac{P_i}{P_{i,*}}\,\frac{P_{i,*}}{P_*}, \end{equation}

where

(2.40)\begin{equation} \frac{P_{i,*}}{P_*} = \frac{R_{i,*}^2L_{i,*}}{\sum_{i=0}^N R_{i,*}^2L_{i,*}}, \quad i=0,1,\ldots,N, \end{equation}

and $P_{i}/P_{i,*}$ given by (2.38) is represented graphically in figure 6 as a function of $R_i/R_{i,*}$ and $\tilde {\tau }_0$, for several values of $n$. Choosing radii that are only 0.3 times smaller than optimal can result in orders of magnitude higher power consumption in comparison to optimised networks (especially for $n=1.5$ and fluids with a low yield stress). Section 3 provides an example of how figure 6 can be used to determine the power consumption of a fluidic network.

Figure 6. Contour plots of the power consumption of a channel, relative to the power of an optimal channel $P_i/P_{i,*}$ (2.38) as functions of $R_i/R_{i,*}$ and $\tilde {\tau }_0$ for $n=0.5$, $n=1.0$ and $n=1.5$.

When the channel radii are non-ideal, this affects the location of the optimal branching point. In that case, one solves the following equation instead of (2.32a,b):

(2.41)\begin{equation} \sum_{i=0}^N \Bigg( \frac{2}{J_{*}} \left(\frac{R_i}{R_{i,*}}\right)^{-(3n+1)} \left(\frac{\psi_{*}}{\psi_i}\right)^n+ \left(\frac{R_i}{R_{i,*}}\right)^{2}\Bigg)R_{i,*}^2\boldsymbol{e}_{i,*}=\textbf{0}. \end{equation}

When $R_i=R_{i,*}$, (2.41) reduces to (2.32a,b). For a derivation, see § B.2.

Equation (2.10) suggests a self-similar velocity profile, as it scales only with $R$ and has the same shape in every optimised channel due to the constant $|\textrm {d}{p}/{\textrm {d}z}|\,R$ and $\phi$. To assess this property, we combine (2.22) into (2.13) to show that $R_{i,*}^{3}/Q_i = \textrm {const.}$ holds for every channel within an optimised branching for all treated fluid rheologies, in the following manner. Substitution of (2.22) into (2.13) gives an equation containing only constant fluid properties and $\phi$. Consequently, $\phi$ is constant for an optimised branching, and this is also a proof that $R^3/Q$ is constant in every channel of the branching. As $\phi$ is constant, $\psi$ and $|\textrm {d}p/{\textrm {d}z}|\,R$ are also constant in every channel. As a result of these invariants, the velocity profiles in the channels of an optimised branching are self-similar.

Furthermore, the wall shear stress averaged over the perimeter $\langle \tau _w \rangle$ is a function of the fluid properties and the terms $R^3/Q$ and $\psi (\phi )$. As a result, the average wall shear stress is constant over the whole branching for laminar flow of all treated fluid models when the branching is optimised by (2.22):

(2.42)\begin{equation} \langle \tau_w \rangle \sim \left(\frac{Q}{R^{3}}\right)^n\frac{\mu'}{\psi^n}. \end{equation}

For the derivation, see Appendix D.

An estimate of the cost factor based on governing costs should be made per situation. One approach is similar to Murray's metabolic cost factor, based on the energy needed to maintain a fluid (such as maintenance of the blood or maintenance of a temperature by heating). Also, as elaborated in the second example in § 3, one could think of a print nozzle network filled with expensive ink, where costs are based on the electricity and material costs. Other cost factors have been described for plants, where the conduit wall volume determines the costs, as the wall should be strong enough to withstand the negative internal pressure (McCulloh et al. 2003), diffusive systems (Zheng et al. 2017), and systems that require drag minimisation (Woldenberg & Horsfield 1986).

Even if a cost factor $\alpha$ is not known explicitly, it is possible to design an optimised network based on other constraints – for example, in 3-D printing, where a flow rate and a nozzle radius at an outlet are specified and are assumed to be optimal. A second approach is making a channel to require a certain average flow velocity or dimensionless number such as the Weber number, which serves as a relation to obtain the corresponding optimal $R$ for given $Q$. For both approaches, the ratio $R^3/Q$ is then known, and based on (2.22), the whole optimal network can be calculated. The cost factor can be made explicit by calculating the optimal $\phi$ from (2.13) and calculating $\alpha$ from (2.22).

3. Examples: optimisation of a branched fluidic network

This first example shows the convenience of the graphical approach in determination of $R^3/Q$. Here, a branched fluidic network inspired on a lubrication system is optimised, which has a laminar flow of grease described by the Bingham model ($\rho =1000\ \textrm {kg}\ \textrm {m}^{-3}$, $\mu '=1.85\ \textrm {Pa}\ \textrm {s}$, $n=1$, $\tau _0=1.0$ Pa), inspired by Westerberg et al. (2010). The temperature of the lubricant is maintained by adding heat, resulting in a cost factor $\alpha = 100\ \textrm {W}\ \textrm {m}^{-3}$. The network has one parent channel (channel 0) and two daughter channels (channels 1 and 2, $N=2$). The nodes are given by $\boldsymbol {x}_0 = (0,0)$ m, $\boldsymbol {x}_1 = (0.0625,0.025)$ m, and $\boldsymbol {x}_2 = (0.0625,-0.025)$ m. The flow rate in the parent channel is $Q_0=10\ \textrm {ml}\ \min ^{-1}$, and the flow rates in the daughter channels are determined to be $Q_1 =Q_2= \frac {1}{2}Q_0$, so the bifurcation is symmetric.

The dimensionless yield stress $\tilde {\tau }_0$ is calculated as $\tilde {\tau }_0=\tau _0/(\mu '\alpha ^n)^{1/(n+1)}$ = 0.074, after which the corresponding optimisation condition $R^3/Q \approx 1.29$ is easily read from figure 4. Alternatively, the optimal value of $\phi$ is calculated from (2.16), (2.25) and (2.30), resulting in $\phi _*=0.068$. With $\phi _*$, the corresponding optimal value for $\psi$ can be calculated using (2.15): $\psi _*=0.910$. With these parameters, (2.28) provides an analytic value for $\tilde {R}^3/\tilde {Q}$, being equal to 1.2889, which is readily rewritten to real units using (2.26), where $(\mu '/\alpha )^{1/(n+1)} = 0.136$ s: we obtain $R^3/Q=0.175$ s. Based on this condition, the main channel radius is calculated to be $R_0 = 3.08$ mm. Knowing that $R^3/Q$ must be constant, based on (2.22), the radii of the daughter channels are determined:

(3.1)\begin{equation} R_1=R_2= R_0 \left(\frac{Q_1}{Q_0}\right)^{{1}/{3}}=0.794 R_0 = 2.45\ \text{mm}. \end{equation}

As a check on whether the flow is still laminar, the largest Reynolds number and critical Reynolds number in the network are calculated from (2.19) and (2.20): $Re' = 0.0186$ and $Re'_{crit}=2357$, respectively. This shows that the flow in the whole network will stay laminar.

From (2.32a,b), it is calculated that the branching point is located at $\boldsymbol {x}_*= (0.030,0.000)$ m. The resulting topology is visualised in figure 7(a), where we make use of the velocity profile given in (2.10). The angles between the channels for this optimised network are calculated using (2.34) and also shown in the figure.

Figure 7. Topology of branching with lubrication flow with $N=2$ daughter channels. (a) Symmetric bifurcation. (b) Asymmetric bifurcation with $Q_1 = \frac {1}{3}Q_0$ and $Q_2 = \frac {2}{3}Q_0$.

To show the effect of an asymmetric flow rate division, the same example is taken, but now with flow rates $Q_1 = \frac {1}{3}Q_0$ and $Q_2 = \frac {2}{3}Q_0$. In addition, the end points are chosen to be $\boldsymbol {x}_0 = (0,0)$ m, $\boldsymbol {x}_1 = (0.053,0.017)$ m and $\boldsymbol {x}_2 = (0.0625,-0.025)$ m. Then the optimisation condition $R^3/Q$ is the same as before, but the daughter channel radii are different. These are $R_1 = R_0(Q_1/Q_0)^{1/3} = 2.14$ mm and $R_2 = R_0(Q_2/Q_0)^{1/3} = 2.69$ mm. As a consequence, the location of the branching point and the corresponding angles between the channels become different. The angles are calculated using (2.34), and the topology is shown in figure 7(b). Here, one sees that the angle between the parent channel and the largest daughter channel has become larger, tending more to the limit of $180^\circ$.

In a second example, we optimise the channel radius of a hypothetical branching nozzle, as for example used in parallel direct-write (Skylar-Scott et al. 2019; see figure 1). A model fluid is described by the Herschel–Bulkley model ($\rho = 1200\ \textrm {kg}\ \textrm {m}^{-3}$, $\mu '= 15\ \textrm {Pa}\ \textrm {s}^n$, $n=0.4$, $\tau _0=2000$ Pa). The cost factor is estimated to be the ratio of the material costs per volume and the electricity costs to pump the fluid for a certain operation time: $\alpha = {c\rho }/{E\,\Delta t}$. The fluid costs , electricity costs , and a print session takes $\Delta t = 18$ min. Then the cost factor is estimated at $\alpha = 1\ \textrm {GW}\ \textrm {m}^{-3}$. The flow rate of the parent channel is $Q_0 = 10\ \textrm {ml}\ \min ^{-1}$, and the channel splits up into $K=4$ levels of symmetric bifurcations. Calculating $\tilde {\tau }_0=\tau _0/(\mu '\alpha ^n)^{1/(n+1)}=0.78$ and using figure 4, we obtain a value for the optimisation condition as $\tilde {R}_{i,*}^3/\tilde {Q}_i \approx 1.14$. Therefore, the optimisation condition in real units becomes $R_{i,*}^3/Q_i = \tilde {R}_{i,*}^3/\tilde {Q}_i (\mu '/\alpha )^{1/(n+1)} = 2.9\times 10^{-6}$ s. Consequently, the optimal channel radii for the four different levels (where due to symmetry, every channel at the same level has the same size) are $R_0 = 79\ \mathrm {\mu } \textrm {m}$, $R_1 = 62\ \mathrm {\mu } \textrm {m}$, $R_2 = 50\ \mathrm {\mu } \textrm {m}$ and $R_3 = 39\ \mathrm {\mu } \textrm {m}$. These channels are unrealistically small, because their manufacture is challenging and the material cost within the channels may be small compared to e.g. depreciation of equipment. However, this example shows clearly how expensive fluids could drive miniaturisation of the system.

In a third example, we investigate how a microfluidic branched network with equally sized channels performs in comparison to an optimised network, for symmetric bifurcations and $K=5$ levels. The fluid in the network is described by the Herschel–Bulkley model ($n=0.5$, $\tau _0=1000$ Pa, $\mu '=200\ \textrm {Pa}\ \textrm {s}^n$; Skylar-Scott et al. 2019). The radius and flow rate of the channel were considered to be optimal for the channels at outlet and have values $R_{opt} = 250\ \mathrm {\mu } \textrm {m}$ and $Q_{opt} = 2\ \textrm {ml}\ \min ^{-1}$. Consequently, all channels have radius $250 \ \mathrm {\mu } \textrm {m}$, and the main parent channel has flow rate $Q_0=32\ \textrm {ml}\ \min ^{-1}$. The flow rate in the other channels is calculated according to equal distribution of $Q_0$. Furthermore, the channel lengths are assumed to be a function of the flow rate: $L_i/L_0 = (Q_i/Q_0)^{2/3}$. The optimisation condition based on the optimal channel becomes $R_{i,*}^3/Q_i = 4.69$ s. From that, the optimal radii are calculated, based on the known flow rate in each channel. Together with the optimisation condition $R_{i,*}^3/Q_i$, one can calculate the properties of the optimal network. Based on the optimal channel, one can calculate $\phi _*=0.14$ and $\psi _*=0.68$ using (2.13), which also yields $J_*=2.24$ using (2.24) and $\tilde {\tau }_{0,*}=0.210$ using (2.30). Calculation of the radius ratio $R_i/R_{i,*}$ for each channel provides the needed input for the power ratio per channel via (2.38) or figure 6. Multiplying each power ratio with the corresponding term $P_{i,*}/P_*$ from (2.40), and summing everything, results in the overall network power ratio $P/P_*= 2.38$. This shows that application of the optimal radii instead of equally sized channels could reduce 58 % of the power consumption in the network.

4. Discussion

The optimal geometry derived above does not incorporate the constraint that the pressure at all end points must be equal. Therefore, flow rates may deviate from the calculated flows for open-ended systems such as nozzles that exit into the air. This issue can be mitigated in general by a pressure control system at the end of every node. Fortunately, however, bifurcating channels that are symmetric with respect to both the flow rate and the geometry always exhibit the same pressure drop, and will therefore provide optimal flow through the network without active control. In addition, this problem is prevented for radius-wise optimised networks by requiring $\sum _i L_i/R_i = \textrm {const.}$ for both branches that originate from a branching point.

Furthermore, Murray's theory (§§ 2.2 and 2.3) is sometimes simplified to the so-called Murray law stating that the cubed parent channel radius is equal to the sum of the cubed daughter channel radii (see e.g. Sherman 1981; Stephenson et al. 2015):

(4.1)\begin{equation} R_{0}^3=\sum_{i=1}^N R_{i}^3. \end{equation}

This relation is obtained by substituting $Q_i$ from the optimisation condition (2.22) into the mass conservation equation (2.1), where the right-hand side of (2.22) is constant for an optimised network. It is important to note that a branching optimised by (2.22) automatically fulfils (4.1), but satisfying (4.1) does not necessarily give the optimal branching. This was already pointed out by Kamiya et al. (1974). Rosenberg (2020) analysed this point critically, and invalidated the conclusion by Sherman (1981) that conservation of the cube of the radii in a branching point is the determinant condition for network optimisation for both symmetric and asymmetric branches, for whatever number of daughter channels. Therefore, (4.1) is not sufficient for determining the optimal radii, except for branchings in which the flow rate in all daughter channels is the same (e.g. symmetric bifurcations). In that case, (4.1) reduces to

(4.2)\begin{equation} \frac{R_i}{R_0}=N^{-{1}/{3}}, \quad i=1,2,\ldots,N. \end{equation}

Finally, flow effects such as entrance effects, non-axisymmetric velocity profiles at junctions, and channel curvature effects have not been taken into account due to the assumption of the channel flows being fully developed. However, if the slenderness condition in the channels is satisfied ($R\ll L$), then the major part of the flow is still governed by the fully developed flow equations. Therefore, the theorem is suitable for overall network design or as an analysis tool, but for design details near junctions, simulations or experiments are recommended.

5. Conclusions and outlook

A generalised approach to minimise power consumption of branched fluidic networks was developed for Newtonian, power-law, Bingham, Herschel–Bulkley and Casson-like fluids in the steady laminar regime.

First, the flow rate in a single channel was analysed, providing the pressure drop over a channel as a function of the fluid model and channel dimensions. Subsequently, the optimal network geometry was obtained by minimising power with respect to the channel radii and the location of the branching point. For all treated fluid rheologies, Murray's optimisation condition for all channels within a network was recovered:

(5.1)\begin{equation} \frac{\tilde{R}_{i,{*}}^{3}}{\tilde{Q}_i}=\mbox{const.,} \quad i=0,1,2,\ldots,N. \end{equation}

The value of this constant is represented graphically in figure 4, enabling optimisation of fluid networks for non-Newtonian systems with a minimum of mathematical analysis for the first time. Detailed analysis of the network provided the following additional results.

1. (i) Insight into the increase of power consumption in case of non-optimal channel radii, revealing potentially large energy gains.

2. (ii) A self-similar velocity profile across the network, resulting in a constant wall shear stress.

3. (iii) Multiple methods to estimate the cost factor to maintain the liquid flow.

4. (iv) Conditions for which optimal flow is maintained without active control of the pressure at every junction.

5. (v) Proof that if the parent channel is laminar, the entire optimised network exhibits laminar flow. The Darcy friction factor is provided.

Future work may further address situation-specific cost factors for engineering systems (Woldenberg & Horsfield 1986), as well as systems with a transition from turbulent to laminar flow regimes.