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. Reference Smink, Venner, Visser and Hagmeijer2023), 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 Reference Weisbach1845)

(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 Reference Garcia and Steffe1986)

(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 Reference Moody1944)). 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 Reference Hanks and Ricks1974)

(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 Reference Everts, Robbertse and Spitholt2022). 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 Reference Swamee and Aggarwal2011) 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. Reference Smink, Venner, Visser and Hagmeijer2023)

(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 Reference Liu, Li and Smits2019). 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 Reference Herschel and Bulkley1926; Chilton & Stainsby Reference Chilton and Stainsby1998; Smink et al. Reference Smink, Venner, Visser and Hagmeijer2023)

(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 Reference Reiner1926; Smink et al. Reference Smink, Venner, Visser and Hagmeijer2023)

(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. Reference Smink, Venner, Visser and Hagmeijer2023)

(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 Reference Colebrook1939) 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 Reference Moody1944) 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 Reference Dodge and Metzner1959). 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 Reference Garcia and Steffe1986) 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 Reference Uylings1977; Sherman Reference Sherman1981; Woldenberg & Horsfield Reference Woldenberg and Horsfield1986; Williams et al. Reference Williams, Trask, Weaver and Bond2008; Kou et al. Reference Kou, Chen, Zhou, Lu, Wu and Fan2014). 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 Reference Murray1926b ; Smink et al. Reference Smink, Venner, Visser and Hagmeijer2023) 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 Reference Uylings1977). (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. Reference Kou, Chen, Zhou, Lu, Wu and Fan2014)) 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 Reference Sherman1981; Stephenson & Lockerby Reference Stephenson and Lockerby2016) ( $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 Reference Uylings1977; Bejan et al. Reference Bejan, Rocha L.A. and Lorente2000; Kou et al. Reference Kou, Chen, Zhou, Lu, Wu and Fan2014; Stephenson & Lockerby Reference Stephenson and Lockerby2016; Zhou et al. Reference Zhou2024). 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.