Optimising branched fluidic networks: a unifying approach including laminar and turbulent flows, rough walls and non-Newtonian fluids

Abstract

Power minimisation in branched fluidic networks has gained significant attention in biology and engineering. The optimal network is defined by channel radii that minimise the sum of viscous dissipation and the volumetric energetic cost of the fluid. For limit cases including laminar flows, high-Reynolds-number turbulence or smooth-channel approximations, optimal solutions are known. However, current methods do not allow optimisation for a large intermediate part of the parameter space which is typically encountered in realistic fluidic networks that exhibit turbulent flow. Here, we present a unifying optimisation approach based on the Darcy friction factor, which has been determined for a wide range of flow regimes and fluid models and is applicable to the entire parameter space: (i) laminar and turbulent flows, including networks that exhibit both flow types, (ii) non-Newtonian fluids (powerlaw, Bingham and Herschel–Bulkley) and (iii) networks with arbitrary wall roughness, including non-uniform relative roughness. The optimal channel radii are presented analytically and graphically. All existing limit cases are recovered, and a concise framework is presented for systematic optimisation of fluidic networks. Finally, the parameter $x$ in the optimisation relationship $Q\propto R^{x}$, with $Q$ the flow rate and $R$ the channel radius, was approximated as a function of the Reynolds number, revealing in which case the entire network can be optimised based on one optimal channel radius, and in which case all radii must be optimised individually. Our approach can be extended to a wide range of fluidic networks for which the friction factor is known, such as different channel curvatures, bubbly flows or specific wall slip conditions.

1. Introduction

Branched networks for fluid transport are omnipresent in biology and engineering. Biological systems, such as vascular networks (Hutchins, Miner & Boitnott 1976; Reichold et al. 2009; Kelch et al. 2015) and the bronchial trees of the lungs (Hooper 1977; Xu et al. 2016; Sznitman 2022), are capable of efficiently transporting heat or mass with low dissipation within limited volumes (Bejan & Lorente 2013). Similarly, fluidic networks provide efficient transport in emerging engineering applications, such as microfluidics (Whitesides 2006), additive manufacturing (Skylar-Scott et al. 2019; Shaqfeh et al. 2023), creation of synthetic vasculatures (Sexton et al. 2023), microreactors (Dong et al. 2021) and multifunctional materials (Zheng et al. 2017). These networks were comprehensively analysed and optimised for Newtonian (Murray 1926a , b ; Kamiya, Togawa & Yamamota 1974; Zamir 1977; Oka & Nakai 1987), power-law (Mayrovitz 1987; Revellin et al. 2009; Stephenson & Lockerby 2016; Miguel 2018) and yield-stress fluids (Ponalagusamy 2012; Smink et al. 2023) in the laminar flow regime. However, a transition to turbulence commonly occurs in natural flows (Olson, Dart & Filley 1970; Stein & Sabbah 1976; Ku 1997; Berger & Jou 2000; Calmet et al. 2016; Ha et al. 2018, 2019) and industrial flows such as district heating (Gumpert, Wieland & Spliethoff 2019; Steinegger et al. 2023), water distribution networks, heat exchangers (Siddiqui & Zubair 2017), the paper making process (Lundell, Söderberg & Alfredsson 2011) and inertial microfluidics (Wang, Yang & Zhao 2014). Similarly, fluidic networks that exhibit both turbulent flow and non-Newtonian behaviour may benefit upscaling of nozzle-based technologies via parallelisation, such as in-air microfluidics (Visser et al. 2018), spray drying (van Deventer, Houben & Koldeweij 2013), prilling (Kamis et al. 2023) or 3D printing of polymer foams (Visser et al. 2019). Finally, analytical optimisation of fluidic networks in various turbulent regimes would provide key validation data for computational fluid dynamics, which is increasingly used to describe (Morris et al. 2016) or design (Sexton et al. 2023) networks with e.g. complex geometries.

The parameter space spanned by the Reynolds number, the roughness and different fluid models is shown in figure 1(a,b). Applications are plotted in figure 1(a), showing the broad range of relevant Reynolds numbers, $Re$ , and relative channel wall roughnesses, $\delta$ , especially in the intermediate- $Re$ turbulent regime (the Reynolds number is defined in (2.4)). Figure 1(b) shows that the energy consumption by the network has been minimised (Murray 1926b ; Smink et al. 2023) only for specific limits of this parameter space. This total energy consists of the sum of the power needed to maintain the flow (product of the pressure drop $\Delta p$ and the volume flow rate $Q$ ) and a volumetric cost function $\alpha V$ , with $\alpha$ a cost factor and $V$ the volume of the channel, for networks as schematically shown in figure 1(c, d). For example, in the case of a vascular network, the cost factor expresses the metabolic energy needed to maintain blood ( $\alpha \approx$ 1 kW m $^{-3}$ ) (Murray 1926b ); for other situations, see Smink et al. (2023). Note that optimising the relative channel radii for a given volume is a related problem, which has been pursued for many cases with e.g. a constructional law (Bejan, Rocha L.A. & Lorente 2000; Miguel & Rocha 2018; Miguel 2018). Instead, we focus on optimisation against a cost factor, which is generalisable for different boundary conditions such as a cost function of the wall surface area (Woldenberg & Horsfield 1986) and readily provides access to the optimal channel radii.

Figure 1. (a) Parameter space of the Reynolds number $Re$ and the relative channel wall roughness $\delta$ for fluidic networks in coronary arteries (Singhal, Henderson & Horsfield 1973; Kassab et al. 1993; Burton & Espino 2019), paper making (Forgacs, Robertson & Mason 1957; Grossman & Carpenter 1968; Moller 1976), water distribution (van der Schans et al. 2015; NEN 2018), inertial microfluidics (Di Carlo 2009; Zhang et al. 2016; Lu et al. 2017), heat exchangers (Towler & Sinnott 2013) and district heating (Gumpert et al. 2019; Steinegger et al. 2023). The colour indicates a Newtonian (blue), power-law (green) or yield-stress (red) fluid. Here, $Re_{crit}$ is assumed to be independent of $\delta$ for this parameter space. (b) Previously optimised parts of the parameter space include laminar flow and low- $Re$ turbulent flow in a smooth channel and a specific condition for complete rough-channel turbulence. (c) Schematic of a single branch with fully developed laminar flow profiles within the channels. (d) Schematic of a branched fluidic network, where a parent channel splits up 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. (e, f) Examples of optimised branched fluidic networks. The colour indicates the Reynolds numbers in the channel. The position coordinates of the begin and endnodes are given as an input; the coordinates of intermediate nodes follows from the optimisation. Also the flow rates in all channels are given as an input. The following quantities are kept constant: $Q_0 = 20$ l min $^{-1}$ , $\rho =1000$ kg m $^{-3}$ , $\alpha =10^{3}$ W m $^{-3}$ and $\tau _0=0$ Pa. Fluid and system parameters are defined in § 2. (e) Newtonian fluid in rough channel ( $\mu '=10^{-3}$ Pa s, $n=1.0$ , $\varepsilon = 10^{-5}$ m, 5 levels, symmetric branching). (f) Power-law fluid in smooth channel ( $\mu '=5\times 10^{-5}$ Pa s $^{1.5}$ , $n=1.5$ , 4 levels, asymmetric branching 1:2).

For a circular pipe with radius $R$ , the power is generally minimised if

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

with $x$ a flow-regime and fluid-model-dependent power; $x = 3$ for laminar flow of an incompressible Newtonian fluid, as further analysed theoretically (e.g. revealing the shear forces, the optimal angles between channels or special cases such as curved pipes or porous channels) (Kamiya et al. 1974; Zamir 1977; Sherman 1981; Miguel 2018) and verified experimentally for e.g. human coronary and cerebral arteries (Hutchins et al. 1976; Rossitti & Löfgren 1993). The value $x = 3$ also holds for laminar flows of different fluid models, including power-law (Mayrovitz 1987; Revellin et al. 2009; Stephenson & Lockerby 2016; Miguel 2018) and yield-stress fluids, such as the Bingham, Herschel–Bulkley and Casson models (Smink et al. 2023), as well as for channels with elliptical (Tesch 2010) and rectangular (Emerson & Barber 2012), or even arbitrary cross-sections (Emerson et al. 2006; Stephenson et al. 2015; Zhou et al. 2024).

For turbulent flows, fluidic networks were optimised only for limit cases and approximations (Uylings 1977; Woldenberg & Horsfield 1986; Bejan et al. 2000; Stephenson & Lockerby 2016). For turbulent flow of Newtonian fluids at intermediate Reynolds numbers $2\times 10^{4} \lt Re\lt 10^{6}$ in a hydraulically smooth channel (shaded for relative wall roughness $\delta \rightarrow 0$ in figure 1 b), $x = 17 / 7$ was derived from an empirical relation for the friction factor (Bejan 2013; Kou et al. 2014). For complete turbulence (the light-shaded area in figure 1(b), which is defined as $Re\gt 3500/\delta$ Moody 1944), it was derived that $x = 7 / 3$ (Uylings 1977; Bejan et al. 2000; Williams et al. 2008; Stephenson & Lockerby 2016). Several previous studies (e.g. Revellin et al. 2009 and Zhou et al. 2024) have generalised the optimisation procedure for friction factors that are a power-law function of $R$ , and show that all these limit cases can be described by a single approach. However, the intermediate part of the parameter space is described by more complex friction factors that cannot be implemented in these existing approaches, such as the Colebrook–White equation. Stephenson et al. (2015) and Stephenson & Lockerby (2016) use a different approach allowing friction factors in arbitrary functional form, but they do not apply this to the intermediate part of the parameter space. This terra incognita represents a major knowledge gap, both fundamentally, as it covers several orders of magnitude of both $Re$ and $\delta$ , and from an application perspective, as we will show that turbulent limit cases are rarely reached for realistic flows.

In this study, we introduce a universal approach for optimising fluidic networks in the entire parameter space. The optimisation procedure is applicable to any flow regime, wall roughness or fluid model for which the Darcy friction factor is known, even if its formulation is implicit. First, § 2 describes non-dimensionalisation of pipe flows, enabling calculation of the optimal Reynolds number within a channel and the corresponding channel radius. Subsequently, the optimal channel radius is derived for laminar and turbulent flow at arbitrary $Re$ , arbitrary channel roughness $\delta$ and for both Newtonian and non-Newtonian fluids in § 3. Section 4 synthesises the results from §§ 2 and 3, resulting in approximations of $x$ as a function of $Re$ for all flow regimes and fluid models. The design procedure is presented in § 5. The conclusions are presented in § 6. The design of the optimal location of branching points within the network is described in Appendix A.1, resulting in fully optimised networks for e.g. rough-wall channels (figure 1 e) or for non-Newtonian fluids (figure 1 f). Finally, an application example is provided in Appendix F.

2. Optimisation of fluidic branching

Consider a branching configuration comprising a parent channel connected to $N$ daughter channels at a branching point denoted as $\boldsymbol{x}$ , as depicted in figure 1(d). The channels are labelled $0$ to $N$ , with index $0$ indicating the parent channel. The effective radii of the channels are represented as $\boldsymbol{R}\equiv (R_0, R_1,\ldots, R_N)$ , while the termination points of the channels are fixed and denoted as $\boldsymbol{x}_i$ for $i=0,1,\ldots, N$ . The flow rates in the daughter channels are specified as $Q_i$ for $i=1,\ldots, N$ , with $Q_0$ corresponding to the flow rate in the parent channel; $Q_0$ is considered 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 ensure mass conservation, assuming incompressibility, the flow rates must satisfy $Q_0=\sum _{i=1}^{N} Q_i.$ The lengths of the channels, $L_i$ , are functions of the branching location, given by $L_i\equiv |\boldsymbol{x}_i-\boldsymbol{x}|$ for $i=0,1,\ldots, N$ . In this work, it is assumed that there is only an axial velocity $u$ . The channels are assumed to be cylindrical with $R\ll L$ . The fluid properties are taken constant and the flow is fully developed.

For optimisation of a fluidic network, the power to maintain the flow and to maintain a fluid is minimised. The power is represented by a cost function depending on the radii and lengths of the channels, and is the sum of the individual channel contributions

(2.1) \begin{equation} P(\boldsymbol{R},\boldsymbol{x}) \equiv \sum _{i=0}^{N} \left (\Delta p_i\, Q_i + \alpha V_i\right ). \end{equation}

Here, $V_i \equiv \pi R_i^{2} L_i$ is the channel volume and the pressure gradient ${\Delta p}/{L}$ is a fluid-type and flow-regime-dependent function. The pressure at the nodes is not defined or constrained. Differentiation of (2.1) to $\boldsymbol{R}$ and $\boldsymbol{x}$ and equating these expressions to 0 provides the optimisation problems for: (i) the channel radii $\boldsymbol{R}$ , and (ii) the location of the branching point $\boldsymbol{x}$ (optimised in Appendix A.1). Optimisation of the channel radius is independent of $\boldsymbol{x}$ , because it is a decoupled problem (Smink et al. 2023), in which the power defined by (2.1) attains a global minimum with respect to the radius if all channels satisfy ${\partial P}/{\partial \boldsymbol{R}}=\boldsymbol{0}$ . Therefore, the optimisation condition for any radius $R_i$ only depends on the power contribution of the corresponding channel and $R_i$ does not depend on the lengths of the channels $L_i$ .

For minimising the power consumption in the network, every single channel has to be optimised individually according to

(2.2) \begin{equation} \frac {\partial {P}}{\partial {R_i}}=\frac {\partial {\Delta p_i}}{\partial {R_i}}Q_i+2\pi \alpha R_iL_i=0, \quad i=0,1,\ldots, N. \end{equation}

The pressure drop $\Delta p$ is represented as a function of a flow rate $Q$ and a radius $R$ using the Darcy–Weisbach equation (Weisbach 1845)

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

Here, $\rho$ is the fluid density, $L$ is the length of the channel and $f$ the Darcy friction factor. The friction factor $f$ as function of the Reynolds number for the different fluid models is presented in figures 2(b)–2(d). The specific equations of $f$ for different flow regimes and fluid models are provided in § 3. In this study, a generalised Reynolds number is defined as (Garcia & Steffe 1986)

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

Figure 2. (a) Fluid models as analysed in this work. (b–d) Friction factor as function of the Reynolds number for different fluid models. The dashed black line indicates the transition from laminar to turbulent flow. (b) Newtonian fluid in rough channels (Moody diagram (Moody 1944)). The blue lines represent constant relative roughness $\delta$ from $10^{-7}, 10^{-6}, \ldots 10^{-1}$ . (c) Power-law fluid in smooth channels. The green-scale lines represent values of the flow index $n$ from $0.2, 0.4, \ldots 1.8$ . (d) Herschel–Bulkley fluid with $n=1.0$ in smooth channels. the red-scale lines represent constant values of the Hedström number from $10^{1}, 10^{2}, \ldots, 10^{10}$ .

The value of $f$ also depends on the Hedström number $He$ (commonly used for description of yield-stress fluids), the relative wall roughness $\delta \equiv \varepsilon /2R$ (where $\varepsilon$ is the absolute channel wall roughness) and the flow index $n$ .

The effective dynamic viscosity $\mu$ describes the resistance of a fluid against shear, and is defined as $\mu = \mu '|\dot {\gamma }|^{n-1}$ , where $\dot {\gamma }$ is the shear rate, $\mu '$ is the flow consistency index and $n$ the flow index, where $0\lt n\lt 1$ represents a shear-thinning and $n\gt 1$ represents a shear-thickening fluid, as shown in figure 2(a). In addition, a fluid may show yield-stress behaviour if $\tau _0\geqslant 0$ (e.g. Bingham and Herschel–Bulkley fluids, see figure 2 a), where the fluid only shears once an effective local shear stress exceeds the yield stress $\tau _0$ . The rheological behaviour for the shear rate of a Herschel–Bulkley fluid as function of the effective local shear stress is described as

(2.5) \begin{equation} \dot {\gamma }(\tau ) = \begin{cases} \text{sign}(\tau _{rz})\left ( \frac {|\tau _{rz}|-\tau _0}{\mu '} \right )^{1/n} & \text{if $|\tau _{rz}|\geqslant \tau _0$},\\ 0 & \text{if $|\tau _{rz}|\lt \tau _0$},\\ \end{cases} \end{equation}

where $\tau _{rz}$ is the local shear stress induced by the flow. Figures 2(b)–2(d) also show the transition from laminar to turbulent flows at the critical Reynolds number (Hanks & Ricks 1974)

(2.6) \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}

where $\psi (\phi, n)$ and $\phi$ are defined in § 3.1. Although an increased relative roughness results in a decrease in $Re_{crit}$ , this effect becomes only dominant around $\delta =0.1$ and larger (Everts, Robbertse & Spitholt 2022). In order to ensure the validity of the used equations, the present study will limit the discussed optimisation results to $\delta \leqslant 0.1$ .

The optimisation condition for the radius of every channel is obtained by substituting (2.3) into (2.2), providing

(2.7) \begin{equation} R^{7}=\frac {\rho Q^{3}}{8\pi ^{3}\alpha } B, \end{equation}

where $B$ is defined as

(2.8) \begin{equation} B \equiv 5f -R\frac {\partial {f}}{\partial {R}}. \end{equation}

Here, $B$ is a function of $f$ and $R$ and thus indirectly also a function of $Re$ , $He$ , $\delta$ and $n$ .

The goal is to obtain an expression for the optimal channel radius as a function of parameters that are independent of $R$ . Therefore, dimensional analysis is carried out to find a dimensionless form of the channel radius as function of radius-independent dimensionless groups (for more details, see Appendix A.2). Scaling of the optimisation problem results in the following 5 dimensionless numbers:

(2.9) \begin{align} \tilde {R} &\equiv \frac {R}{\rho ^{-\frac {1}{2}} \mu ^{\prime \frac {3}{2(n+1)}} \alpha ^{\frac {n-2}{2(n+1)}}},\quad \tilde {Q} \equiv \frac {Q}{\rho ^{-\frac {3}{2}} \mu ^{\prime \frac {7}{2(n+1)}} \alpha ^{\frac {3n-4}{2(n+1)}}},\quad \tilde {\tau }_0 \equiv \frac {\tau _0}{\mu ^{\prime \frac {1}{n+1}} \alpha ^{\frac {n}{n+1}}}, \nonumber\\\tilde {\varepsilon } & \equiv \frac {\varepsilon }{\rho ^{-\frac {1}{2}} \mu ^{\prime \frac {3}{2(n+1)}} \alpha ^{\frac {n-2}{2(n+1)}}} \text{ and } n .\end{align}

To align with the literature, the Reynolds number will be used instead of $\tilde {R}$ . This choice will simplify the later computations in § 3 and enables direct characterisation of the flow in the channel. The Reynolds number is rewritten in dimensionless groups as follows:

(2.10) \begin{equation} Re=a(n)\frac {\tilde {Q}^{2-n}}{\tilde {R}^{4-3n}}, \end{equation}

where $a(n)$ is defined as

(2.11) \begin{equation} a(n) \equiv \frac {8}{\pi ^{2-n}}\left (\frac {n}{3n+1}\right )^n. \end{equation}

The network is described by the dimensionless numbers $Re$ , $\tilde {Q}$ , $\tilde {\tau }_0$ , $\tilde {\varepsilon }$ and $n$ . As $\tilde {\tau }_0$ , $\tilde {\varepsilon }$ and $n$ are fluid or system parameters that are constant in a network, $Re$ and $\tilde {Q}$ are the only flow-rate-dependent dimensionless numbers, which is convenient for network optimisation. Described with the new dimensionless numbers, $f$ will now be a function of 5 dimensionless numbers instead of 4. Therefore, in the new domain, there will be a curve along which $f$ is constant, removing one degree of freedom (see Appendix A.3).

For given fluid, flow and system properties, an optimal channel radius can be calculated via the Reynolds number using the non-dimensionalised optimisation condition (2.7)

(2.12) \begin{equation} (8\pi ^{3})^{4-3n}\, a^{7} \, \tilde {Q}^{2(n+1)}= \tilde {B}^{4-3n}Re^{7}, \end{equation}

where $\tilde {B}$ is a function of $\tilde {f}$ , which is a function of $Re$ , $\tilde {\tau }_0$ , $\tilde {Q}$ , $\tilde {\varepsilon }$ and $n$ . Hence, (2.8) can be rewritten as

(2.13) \begin{equation} \tilde {B} = 5\tilde {f}-R\frac {\partial {\tilde {f}}}{\partial {R}} = 5\tilde {f}-R\,\frac {\partial {\tilde {f}}}{\partial {Re}}\frac {\partial {Re}}{\partial {R}} = \tilde {f}\left (5 +(4-3n)\frac {Re}{\tilde {f}} \,\frac {\partial {\tilde {f}}}{\partial {Re}}\right ), \end{equation}

where differentiation of (2.4) provides ${\partial {Re}}/{\partial {R}} = -(4-3n) {Re}/{R}$ . Therefore, knowing $\tilde {\tau }_0$ , $\tilde {Q}$ , $\tilde {\varepsilon }$ and $n$ , one can determine the optimal $Re$ for a channel.

We conclude this section by introducing dimensionless groups that will be used in § 3. The Hedström number $He$ is defined (Swamee & Aggarwal 2011) and rewritten as follows:

(2.14) \begin{equation} He \equiv \frac {4R^{2}\rho }{\mu '}\left (\frac {\tau _0}{\mu '}\right )^{\frac {2-n}{n}} =4\left (a(n)\frac {\tilde {Q}^{2-n}}{Re}\right )^{\frac {2}{4-3n}}\tilde {\tau }_0^{\frac {2-n}{n}}. \end{equation}

The relative roughness $\delta$ becomes in the new dimensionless groups

(2.15) \begin{equation} \delta \equiv \frac {\varepsilon }{2R} =\frac {1}{2}\, \tilde {\varepsilon } \left ( \frac {Re}{a(n) \tilde {Q}^{2-n}}\right )^{\frac {1}{4-3n}}. \end{equation}

Characteristic for a yield-stress fluid is the formation of a plug in the centre of the channel, where the local shear stress $\tau _{rz}$ does not exceed the yield stress $\tau _0$ . The associated plug radius $R_p$ is characterised by

(2.16) \begin{equation} R_p \equiv \frac {2\tau _0}{\dfrac {\Delta p}{L}}. \end{equation}

Non-dimensionalisation of the plug radius yields the dimensionless plug radius $\phi$ , which often appears in descriptions for $f$ in the case of yield-stress fluids

(2.17) \begin{equation} \phi \equiv \frac {R_p}{R} = \frac {2\tau _0}{\dfrac {\Delta p}{L}R}. \end{equation}

By definition, $\phi$ can also be expressed in terms of the Reynolds and Hedström number and in the new dimensionless groups

(2.18) \begin{equation} \phi = \frac {64}{f\, Re}\left (\frac {2He}{Re}\left (\frac {n}{3n+1}\right )^{2}\right )^{\frac {n}{2-n}} = 8\pi ^{2} a(n)^{\frac {4}{4-3n}}\frac {\tilde {\tau }_0}{\tilde {f}\, Re}\left (\frac {\tilde {Q}^{2}}{Re^{3}}\right )^{\frac {n}{4-3n}}. \end{equation}

Up to here, no assumptions are made with respect to the flow type.

3. Network optimisation for different flow regimes and fluid models

The following flow regimes and fluid models are discussed in subsections:

• § 3.1 Laminar flow of a Newtonian and non-Newtonian fluids.

• § 3.2 Turbulent flow of Newtonian fluids, rough and smooth channels.

• § 3.3 Turbulent flow of non-Newtonian fluids, smooth channel, for power-law (§ 3.3.1) and Herschel–Bulkley fluids (§ 3.3.2).

For each fluid model, we will obtain expressions for $\tilde {f}$ , and thereafter for $\tilde {B}$ , which via (2.12) provides the optimal Reynolds number and therefore the optimal radius $R$ for each channel via (2.4) and (2.9) . As many expressions for $\tilde {f}$ and $\tilde {B}$ are implicit, finding an optimal channel radius analytically is often impossible. Therefore, we also provide these results in the form of contour plots for the optimisation condition for $Re$ as a function of $\tilde {Q}$ and $\tilde {\tau }_0$ , $\tilde {\varepsilon }$ or $n$ .

3.1. Laminar flow of Newtonian and non-Newtonian fluids

For laminar flow through a circular channel results in the friction factor (Smink et al. 2023)

(3.1) \begin{equation} \tilde {f} =\frac {64}{Re \,\psi (\phi, n)^n}, \end{equation}

where $\psi$ is a dimensionless flow rate that ranges from 0 to 1. For laminar flow, the friction factor is often considered to be independent of the channel wall roughness, but this is only true for $\delta \lt 0.01$ . When the roughness becomes of the order of the channel size, then the resulting constricted diameter amplifies the friction factor (Liu, Li & Smits 2019). In this section, it is assumed that the relative channel wall roughness is sufficiently small ( $\delta \lt 0.01$ ), such that the influence of the wall roughness on the friction factor is negligible. For a Herschel–Bulkley fluid, $\psi$ is computed by integration of the velocity field over the channel’s cross-section as (Herschel & Bulkley 1926; Chilton & Stainsby 1998; Smink et al. 2023)

(3.2) \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$ ), $\psi$ reduces to (Reiner 1926; Smink et al. 2023)

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

and for the case of $\tau _0=0$ (Newtonian and power-law fluids, $\phi =0$ ), $\psi$ further reduces to 1.

Calculation of the optimisation condition results in the following expression for $\tilde {B}$ :

(3.4) \begin{equation} \tilde {B} = \tilde {f}J(\phi, n), \end{equation}

where the function $J(\phi, n)$ is defined as

(3.5) \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 (3.2) into (3.5) results in the following expression for $J$ :

(3.6) \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 to

(3.7) \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 the yield limit ( $\phi =1$ ), one obtains $J=1$ .

Substitution of $\tilde {B}$ into (2.12) reduces the optimisation condition to

(3.8) \begin{equation} Re^{3(n+1)}=\left (\frac {\pi ^3}{8}\right )^{4-3n}\, a^7\, \tilde {Q}^{2(n+1)} \left (\frac {\psi (\phi, n)^n}{J(\phi, n)}\right )^{4-3n}. \end{equation}

Together with (2.18) and (3.1), $Re$ is then explicitly solved for given $\tilde {Q}$ , $\tilde {\tau }_0$ and $n$ , providing the optimal radius via (2.4) and (2.9). Figures 8(a)–8(c) in Appendix E present a contour plot for the optimal $Re$ as function of $\tilde {Q}$ and $\tilde {\tau }_0$ for different values of $n$ . Note that the number of parameters can be reduced even further, by combining $\tilde {Q}$ and $Re$ , to a parameter which is only a function of given $\tilde {\tau }_0$ and $n$ . Then the optimisation condition becomes (Smink et al. 2023)

(3.9) \begin{equation} \frac {\tilde {R}^3}{\tilde {Q}} =\left ( a^3\frac {\tilde {Q}^2}{Re^3}\right )^{\frac {1}{4-3n}}= \frac {1}{\pi }\left (\left (\frac {3n+1}{n}\right )^n \frac {J(\phi, n)}{\psi (\phi, n)^n}\right )^{1/(n+1)}, \end{equation}

where the right-hand side is constant for an optimised network. Plotting ${\tilde {R}^3}/{\tilde {Q}}$ as function of $\tilde {\tau }_0$ and $n$ gives the contour plot in figure 8(d) in Appendix E. This figure can be made for laminar flows because $\tilde {R}^3/\tilde {Q}$ is constant over the entire network, which is highly useful for design of optimised branched fluidic networks, as only one optimisation condition has to be calculated.

3.2. Turbulent flow of Newtonian fluid, rough and smooth channel

We now demonstrate that the proposed optimisation method is universal if the friction factor is known, by optimising networks for diverse turbulent flow regimes. For turbulent flow of a Newtonian fluid, the Colebrook-White equation (Colebrook 1939) describes the friction factor for both low- and high-turbulent flow ( $Re = [Re_{crit},\infty \rangle$ ) in hydraulically smooth and rough pipes using the following implicit equation:

(3.10) \begin{equation} f = \left \{-2.0\log _{10}\left (\frac {\delta }{3.7}+\frac {2.51}{Re\sqrt {f}}\right )\right \}^{-2}. \end{equation}

Equation (3.10) underlies the well-known Moody diagram (Moody 1944) as shown in figure 2(b). Rewriting (3.10) in terms of the dimensionless numbers gives

(3.11) \begin{equation} \tilde {f} = \left \{-2.0\log _{10}\left (\frac {1}{3.7}\,\frac {\pi }{4}\frac {\tilde {\varepsilon }Re}{\tilde {Q}}+\frac {2.51}{Re\sqrt {\tilde {f}}}\right )\right \}^{-2}. \end{equation}

Differentiation to $Re$ gives

(3.12) \begin{equation} \frac {\partial {\tilde {f}}}{\partial {Re}} = -\frac {\tilde {f}}{Re}\frac {2\left (\dfrac {1}{3.7}\,\dfrac {\pi }{4}\dfrac {\tilde {\varepsilon }Re}{\tilde {Q}}-\dfrac {2.51}{Re\sqrt {\tilde {f}}}\right )}{\left (\dfrac {1}{3.7}\,\dfrac {\pi }{4}\dfrac {\tilde {\varepsilon }Re}{\tilde {Q}}+\dfrac {2.51}{Re\sqrt {\tilde {f}}}\right )\ln \left (\dfrac {1}{3.7}\,\dfrac {\pi }{4}\dfrac {\tilde {\varepsilon }Re}{\tilde {Q}}+\dfrac {2.51}{Re\sqrt {\tilde {f}}}\right )-\dfrac {2.51}{Re\sqrt {\tilde {f}}}}, \end{equation}

resulting in

(3.13) \begin{equation} \tilde {B} = \tilde {f}\left (5-\frac {2\left (\dfrac {1}{3.7}\,\dfrac {\pi }{4}\dfrac {\tilde {\varepsilon }Re}{\tilde {Q}}-\dfrac {2.51}{Re\sqrt {\tilde {f}}}\right )}{\left (\dfrac {1}{3.7}\,\dfrac {\pi }{4}\dfrac {\tilde {\varepsilon }Re}{\tilde {Q}}+\dfrac {2.51}{Re\sqrt {\tilde {f}}}\right )\ln \left (\dfrac {1}{3.7}\,\dfrac {\pi }{4}\dfrac {\tilde {\varepsilon }Re}{\tilde {Q}}+\dfrac {2.51}{Re\sqrt {\tilde {f}}}\right )-\dfrac {2.51}{Re\sqrt {\tilde {f}}}}\right ). \end{equation}

Substituting $\tilde {B}$ into (2.12) enables implicit solving $Re$ for given $\tilde {Q}$ . Contour plots of the optimal $Re$ as function $\tilde {Q}$ and $\tilde {\varepsilon }$ are presented in figure 3(a). Limit cases of the Colebrook–White equation include complete turbulence (for which Von Kármán’s formula applies), turbulence in smooth channels and low-Reynolds-number turbulence in smooth channels (e.g. Blasius’ formula). Expressions of the optimal channel radius for these cases are derived in Appendix B.

Figure 3. Optimal $Re$ for different fluid models. The colour-scale contour lines represent the Reynolds numbers $2\times 10^{x}$ , $3\times 10^{x}$ , … $9\times 10^{x}$ with decreasing brightness. (a) Contour plot of the optimal $Re$ as function of $\tilde {Q}$ and $\tilde {\varepsilon }$ for a Newtonian fluid in a rough-wall channel. (b) Contour plot of the optimal $Re$ as function of $\tilde {Q}$ and $n$ for a power-law fluid in a smooth-wall channel. (c) Contour plot of the optimal $Re$ as function of $\tilde {Q}$ and $\tilde {\varepsilon }$ for a Herschel–Bulkley fluid ( $n=1$ ) in a smooth-wall channel.

3.3. Turbulent flow of non-Newtonian fluids, smooth channel

3.3.1. Power-law fluid

Turbulent pipe flow of non-Newtonian power-law fluids at low Reynolds numbers in a smooth circular channel ( $\varepsilon \rightarrow 0$ ) was analysed by, amongst others, Dodge and Metzner (Dodge & Metzner 1959). They modified the Von Kármán equation for turbulent Newtonian pipe flow, resulting in an implicit relation for the Darcy friction factor $f$

(3.14) \begin{equation} \frac {2}{\sqrt {f}}=\frac {4}{n^{0.75}}\log _{10}\left (Re\left (\frac {f}{4}\right )^{1-n/2}\right )-\frac {0.4}{n^{1.2}} .\end{equation}

The value of $f$ is shown as a function of $Re$ and $n$ in figure 2(c).

Implicit differentiation of $\tilde {f}$ to $Re$ gives

(3.15) \begin{equation} \frac {\partial {\tilde {f}}}{\partial {Re}} = -\frac {\tilde {f}}{Re}\left ( \frac {n^{0.75}\ln 10}{4\sqrt {\tilde {f}}}+1-\frac {n}{2} \right )^{-1}, \end{equation}

resulting in

(3.16) \begin{equation} \tilde {B} = \tilde {f}\left (5+(3n-4)\left ( \frac {n^{0.75}\ln 10}{4\sqrt {\tilde {f}}}+1-\frac {n}{2} \right )^{-1}\right ). \end{equation}

Here, $\tilde {B}$ is still a function of $\tilde {f}$ , which is in turn an implicit function of $Re$ governed by (3.14). Therefore, $\tilde {B}$ is governed by $Re$ and the flow index $n$ . Consequently, with this expression for $\tilde {B}$ , the optimal value of $Re$ in a channel can be calculated as function of $\tilde {Q}$ and $n$ . A contour plot for the optimal Reynolds number in a channel is given in figure 3(b).

3.3.2. Herschel–Bulkley fluid

For a turbulent flow of a Herschel–Bulkley fluid in a smooth circular channel ( $\varepsilon \rightarrow 0$ ), Torrance (Garcia & Steffe 1986) developed a relationship for the friction factor. This relation for the Darcy friction factor is given by

(3.17) \begin{equation} \frac {2}{\sqrt {\tilde {f}}}=0.45-\frac {2.75}{n}+\frac {1.97}{n}\ln (1-\phi ) +\frac {1.97}{n}\ln \left (Re\left (\frac {3n+1}{4n}\right )^n \left (\frac {\tilde {f}}{4}\right )^{1-\frac {n}{2}}\right ). \end{equation}

Here, $\phi$ as given in (2.18) is used in the calculations. Figure 2(d) shows $f$ as function of $Re$ for $n=1.0$ and a range of values of $\tilde {\tau }_0$ .The $f$ - $Re$ plots for other values of $n$ are presented in figure 7 in Appendix E.

Implicit differentiation of $\tilde {f}$ to $Re$ leads to

(3.18) \begin{equation} \frac {\partial \tilde {f}}{\partial Re} = -\frac {\tilde {f}}{Re}\frac {1+\dfrac {\phi }{1-\phi }\dfrac {4}{4-3n}}{\dfrac {n}{1.97}\dfrac {1}{\sqrt {\tilde {f}}}+\dfrac {\phi }{1-\phi }+1-\dfrac {n}{2}}, \end{equation}

giving an expression for $\tilde {B}$

(3.19) \begin{equation} \tilde {B} = \tilde {f}\left (5-\frac {4-3n+\dfrac {4\phi }{1-\phi }}{\dfrac {n}{1.97}\dfrac {1}{\sqrt {\tilde {f}}}+\dfrac {\phi }{1-\phi }+1-\dfrac {n}{2}}\right ). \end{equation}

Contour plots for the optimal Reynolds number in a channel as function of $\tilde {\tau }_0$ and $\tilde {Q}$ are given for $n=1.0$ in figure 3(c) and the results for $n=0.5$ and $n=1.5$ are provided in figure 9 in Appendix E. An additional fluid model covering laminar and turbulent flow of a Bingham fluid is presented in Appendix C.

4. Scaling of the channel radii

In the literature, many attempts have been made to find the proportionality between the flow rate $Q$ and the channel radius $R$ for optimised networks, in the form of (1.1) (Uylings 1977; Sherman 1981; Woldenberg & Horsfield 1986; Williams et al. 2008; Kou et al. 2014). In the optimisation method of the present study, only $Re$ and $\tilde {Q}$ contain parameters involving $R$ and $Q$ , while $\tilde {\tau }_0$ , $\tilde {\varepsilon }$ and $n$ are independent of $R$ and $Q$ . Therefore, when knowing the proportionality between $\tilde {Q}$ and $Re$ , $x$ can be calculated analytically. However, this only holds if $\tilde {B} \propto Re^{c_1}\tilde {Q}^{c_2}$ , which is only true for relatively simple descriptions of $f$ applicable to laminar flows or limit cases of the turbulent regime (e.g. Blasius’ formula, Appendix B.3). Therefore, in the following, $x$ will be calculated by locally approximating $\tilde {B}$ by a power function via

(4.1) \begin{equation} x = \displaystyle\frac {R}{Q}\displaystyle\frac {\partial Q}{\partial R}= \displaystyle\frac {(4-3n)\dfrac {Re}{\tilde {Q}}\dfrac {\partial {\tilde {Q}}}{\partial {Re}}}{(2-n)\dfrac {Re}{\tilde {Q}}\dfrac {\partial {\tilde {Q}}}{\partial {Re}}-1}. \end{equation}

By calculating the fluid-model-specific $({Re}/{\tilde {Q}})\, {\partial {\tilde {Q}}}/{\partial {Re}}$ , one obtains $x$ in the relevant ( $Re,\tilde {\tau }_0,\tilde {\varepsilon },n$ )-space, as shown in figure 4 for the different fluid models (the underlying equations are shown in Appendix D). Known limit cases are shown in black; all coloured results are new. For laminar flow of all treated fluid models, $x$ becomes 3, as expected (Murray 1926b ; Smink et al. 2023) For networks that span ranges of $Re$ for which $x$ is (almost) constant, the value of $x$ is readily used to scale all channels in the network if one optimal diameter is known, using $Q\propto R^{x}$ .

Figure 4. Plot of $x$ (4.1) as function of $Re$ for the different fluid models as discussed in § 3. The dashed lines for $x=3$ and $x=7/3$ show the expected limit cases for laminar flow and high-turbulent flow, respectively (Uylings 1977). (a) Turbulent flow of Newtonian fluids, described by Colebrook–White ( $\tilde {\varepsilon } \in [10^{-5},10^{1},10^{2},10^{3},10^{4},10^{5}]$ ), Blasius’ formula (B12) and an empirical relation (B14). The blue $\times$ -symbols correspond with the optimised channels in figure 1(e), showing that the values for $x$ change significantly for that network. (b) Plot of $x$ as function of $Re$ for turbulent flow of a Newtonian fluid in a hydraulically smooth channel, described by (3.10) for the limit of $\tilde {\varepsilon }=0$ . For optimised networks, high- $Re$ turbulent flow in smooth channels will for realistic $Re$ never reach $x=7/3$ . (c) Plot of $x$ as function of $\delta = \varepsilon /2R$ for high- $Re$ turbulent flow of a Newtonian fluid in a rough channel, described by (B1). For optimised networks, high- $Re$ turbulent flow in channels with finite roughness will never reach $x=7/3$ . (d) Turbulent flow of a power-law fluid (Metzner & Dodge) ( $n\in [0.2,0.4,1.0,1.4,1.8]$ ). (e) Turbulent flow of a Herschel–Bulkley fluid (Torrance) ( $\tilde {\tau }_0 \in [0.1,0.25,1,2.5,10,25]$ with $n= 1.0$ ). Plots for $n=0.5$ and $n=1.5$ are presented in figure 10 in Appendix E.

For a Newtonian fluid (figure 4 a), at the transition from laminar to turbulent flow, $x$ drops quickly from 3 to around 2.5 for hydraulically smooth channels. We find that the low- $Re$ approximations (Blasius, $x=27/11$ , and an empirical relation, $x=17/7$ (Kou et al. 2014)) correspond to the average value of $x$ over their applicable $Re$ -range of the Colebrook–White equation for a smooth channel. Although these approximations and the resulting values of $x$ are valuable for networks that span a limited range of the Reynolds number, the network shown in figure 1(e) represents a case where $x$ is not applicable for determination of the optimal channel radii. Here, the optimal network partly corresponds to laminar flow and partly to turbulent flow. As these regimes result in strongly different values of $x$ , as shown by the markers ( $\times$ ) in figure 4(a), the so-called Murray’s law (Sherman 1981; Stephenson & Lockerby 2016) ( $R_0^{x} = \sum _i R_i^{x}$ , for one branching point) will not hold. Therefore, the optimal radius of each channel must be calculated individually.

Figure 4(a) also shows that $x$ increases for increasing wall roughness. The lines are clipped for $\delta = 0.1$ , corresponding to $x \approx 2.71$ . Analysing the Colebrook–White equation for the smooth-channel limit ( $\tilde {\varepsilon }=0$ ) shows that, for extremely large $Re$ , $x$ will indeed approach $x=7/3$ , as shown in figure 4(b). The high- $Re$ limit for non-zero roughness, however (figure 4(c), $\tilde {\varepsilon }= \text{const.}$ ), shows that $x=7/3$ will never be reached if a finite value of $\delta$ is included in the optimisation process.

When using the Colebrook–White equation, one has to make an assumption about the wall roughness. Figure 4(a) shows the results for a fixed dimensionless absolute roughness $\tilde {\varepsilon }$ . One can also choose to fix the relative wall roughness $\delta$ , which has often been an implicit assumption in previous works (Uylings 1977; Bejan et al. 2000; Kou et al. 2014; Stephenson & Lockerby 2016; Zhou et al. 2024). Then $\delta$ is excluded from the optimisation procedure. Only in this case, it holds for complete turbulence in rough-wall channels that $x=7/3$ (figure 4 c, $\delta = \text{const.}$ ). The assumption of a fixed $\delta$ requires that the absolute roughness $\varepsilon$ is proportional to the optimised channel radius for every channel in the network. However, many networks consist of a given wall material with uniform absolute roughness $\varepsilon$ . Even when there is a changing absolute roughness in the network, the choice for the use of $\tilde {\varepsilon }$ is preferred and then the graphical approach (figure 3 a) can be used for optimisation of the channel radii.

Figure 4(d) shows $x$ for turbulent flow of a power-law fluid in smooth channels. The behaviour of $x$ is similar to that of a Newtonian fluid in a smooth channel, but increases slightly with increasing $n$ . Figure 4(e) shows the results for a Herschel–Bulkley fluid ( $n=1.0$ ), which appears to have a complete transition from $x=3$ to $x\rightarrow 7/3$ , with a remarkable bump in the transition, which is related to the sharp corner around $10^1\lt \tilde {\tau }_0\lt 10^2$ as observed in the optimisation plot (figure 3 c). The transition occurs at larger $Re$ for increasing $\tilde {\tau }_0$ , and for large $Re$ , all lines converge. Results for $n=0.5$ and $n=1.5$ are presented in figure 10 in Appendix E.

The transitional behaviour of $x$ as presented in figure 4 stresses its importance for network optimisation when considering turbulent flow. For domains of $Re$ where $x$ is approximately constant, figure 4 enables making a quick estimation of the optimal channel radii within a network. If one channel within a network has been optimised, then the other channels can be optimised using $Q\propto R^{x}$ . For domains of $Re$ where $x$ clearly non-constant, the full computations of the equations or the use of the graphical approach of figure 3 is required.

5. Synthesis: design of an optimised fluidic network

The framework described in §§ 2, 3 and 4 leads to the following approach for optimisation of branched fluidic networks (such as shown in figure 1 e,f):

1. (i) Determine the fluid model and the cost factor $\alpha$ , and calculate $\tilde {\tau }_0$ , $\tilde {\varepsilon }$ and $n$ . Calculate $\tilde {Q}$ for all the to be optimised channels (§ 2).

2. (ii) Determine the optimal radius $R$ in one channel.

   1. (a) First, determine the optimal $Re$ by one of the following options:

1. – Read $Re$ from an optimisation contour plot (e.g. Figure 3 a) for given $\tilde {Q}$ .

2. – Calculate $\tilde {B}(\tilde {f},Re,\tilde {Q},\tilde {\tau }_0,\tilde {\varepsilon },n)$ via (2.13) and then solve (2.12) for $Re$ for given $\tilde {Q}$ .

1. (b) Subsequently, calculate $R$ from the calculated $Re$ by using (2.4).

1. (iii) Obtain all channel radii:

   1. – By repeating step 2 for each channel, or

   2. – by obtaining $x$ from figure 4 and calculating all channel radii using $Q\propto R^{x}$ . This option only applies if $x$ is (almost) constant over the applicable range of $Re$ .

2. (iv) Calculate the optimal location of the branching points (Appendix A.1).

This approach is illustrated by an example of the optimisation of a fluidic network with turbulent flow of a Newtonian fluid in rough-wall channels, presented in Appendix F.

6. Conclusion

A unifying design procedure for optimisation of branched fluidic networks was presented. The optimal channel radius is obtained by minimising the power consumption for maintaining the flow and minimising the volume simultaneously, which is carried out for all channel segments individually. The approach is applicable to every flow configuration for which the Darcy friction factor is available, even for non-power-law or implicit formulations.

Network optimisation was carried out for the parameter space spanned by the Reynolds number (laminar and turbulent flows), the wall roughness and for different rheologies. This parameter space is widely applicable to e.g. heat exchangers or (inertial) microfluidics, but was hardly covered in the literature, as shown in figure 1(a,b). First, the entire optimisation problem is non-dimensionalised in terms of $Re$ , $\tilde {Q}$ , $\tilde {\tau }_0$ , $\tilde {\varepsilon }$ and $n$ . For every channel, the optimal Reynolds number and the optimal channel radius are obtained as a function of the dimensionless flow rate $\tilde {Q}$ . To minimise the required mathematical analysis, graphs are presented that readily provide the optimal Reynolds number for every fluid model.

Subsequently, the value of $x$ in the proportionality $Q\propto R^{x}$ was evaluated for the full parameter space of the Reynolds number and channel roughness for different fluid rheologies. The well-known limit values $x=3$ for laminar flow and $x=7/3$ for high- $Re$ turbulence were recovered (Uylings 1977; Bejan et al. 2000; Kou et al. 2014; Stephenson & Lockerby 2016; Zhou et al. 2024). However, as the latter value is hardly reached in practice, the intermediate values of $x$ were calculated and plotted. This revealed in which case the entire network can be optimised based on one optimal channel radius, and in which case all radii must be optimised individually. Finally, a design procedure of optimised fluidic networks that incorporates all aforementioned elements was presented.

Our unifying approach for optimisation of fluidic networks can be extended to other configurations with known friction factors, including non-circular cross-sections, bended pipes, bubbly flows, porous walls (Miguel 2018), porous media (MacDonald et al. 1991) or slip flow (Morini & Spiga 1998). As such, it enables optimisation of many realistically encountered configurations in nature and engineering.