2.1. Microfluidic device fabrication

Our rectangular channel devices (figures 2b and 2e) were made by moulding PDMS (Sylgard 184, Dow Chemical, MI, USA) on a patterned silicon wafer mould fabricated via standard lithography techniques. The PDMS was prepared by thoroughly mixing the base and curing agent in a 10 : 1 by weight ratio and was cured on the wafer in an oven for $1\ \mathrm {hr}$ at $65\,^\circ \mathrm {C}$. Inlet and outlet holes were created using a biopsy punch ($2\ \mathrm {mm}$, Integra LifeSciences, NJ, USA). The membrane was made by spin coating (WS-650-23, Laurell, PA, USA) PDMS on a clean wafer at $2500\ \mathrm {rpm}$ for $30$ seconds and was subsequently cured in the oven. The microfluidic device comprising two channels separated by the membrane was assembled by first removing the PDMS channel slabs from the moulds and punching the inlet and outlet holes. One slab was then bonded to the membrane (still attached to the wafer) via plasma activation (PDC-002, Harrick Plasma, NY, USA) on high radio frequency power for $30\ \mathrm {s}$. After bonding, the next step was to remove the membrane (now bound to a channel) from the wafer. To do this, we first filled the channel with water using a syringe. Then, we cut around the device's perimeter with a scalpel to release the membrane from the wafer, allowing us to peel off the membrane gently. By filling the channel with water, we mitigate the risk of the membrane collapsing into the channel when the membrane is peeled off the wafer. The other channel slab was bonded to the other side of the membrane via the same procedure. Note that we aligned the channels perpendicularly by eye, although a purpose-built alignment set-up (e.g. Li Reference Li2015) would be feasible in ensuring optimal alignment and centring.

To fabricate the rounded channel devices (figures 2c and 2f), we followed the procedure by Hongbin et al. (Reference Hongbin, Guangya, Siong, Shouhua and Feiwen2009) to make channel moulds. Briefly, the method consists of inflating a microchannel with a thin membrane lid and casting PDMS on top of the inflated membrane to yield a new channel with a height profile identical to the membrane deflection. However, instead of casting PDMS, we cast a fast-curing polymer (Elite Double 22, Zhermack, Italy) mixed 1 : 1 by weight. When curing was completed, we removed the polymer slab (with the channel imprint) and cast a liquid plastic resin (FormCast Burro, FormX, Netherlands) mixed 1 : 1 by weight to yield a rigid, rounded channel mould for our subsequent PDMS casting. This modification to Hongbin et al. (Reference Hongbin, Guangya, Siong, Shouhua and Feiwen2009) method allows for multiple devices to be cast from the same mould without repeating the inflation step. We used an inflation pressure of $\Delta p_I = 2\ \mathrm {kPa}$ to produce our rounded channel mould, which yielded a centre height of $h_0=250\pm 10\ \mathrm {\mu }\mathrm {m}$ (figure 2f). The advantage of this geometry is that it enables conforming contact between the deflected membrane and the bottom of the channel. Therefore, in this conforming geometry, the membrane can occlude the channel at a lower pressure relative to the rectangular channel geometry (Unger et al. Reference Unger, Chou, Thorsen, Scherer and Quake2000), thus allowing us to study the fluid–structure interactions in multiple connected junctions within our operating pressure range. To connect device channels, we used polyvinylchloride (PVC) tubing with internal and outer diameters of $1.0\ \mathrm {mm}$ and $2.0\ \mathrm {mm}$, respectively (GRA-GL0100005, Mikrolab Aarhus, Denmark). The PVC tubing's ends were cut with a razor to create a tapered tip, allowing easy insertion into the PDMS devices.

2.2. Experimental set-up

In characterizing our single or multiple junction devices, we measured the fluid flow rate $Q$ through the device as a function of applied pressure, $\Delta p \approx 0\unicode{x2013}40\ \mathrm {kPa}$. We used a pressure controller (LineUp Flow EZ^TM, Fluigent, France) to sweep the applied pressure. The inlet of the controller was connected to a source of pressurized air, and its outlet to a closed container containing water with viscosity $\eta = (0.9\pm 0.1)\times 10^{-3}\ \mathrm {Pa\ s}$. A straw tube inserted into the container directed the pressurized water into a pressure sensor (26PC Flow-through, HoneyWell, NC, USA), and then into our microfluidic device. The device's outlet was connected to another pressure sensor and finally into a flow meter (SLF3s-1300F, Sensirion, Switzerland) connected to a reservoir held at atmospheric pressure. Using two pressure sensors kept at a constant altitude, we could accurately measure the pressure drop across the device $\Delta p$, while the flow meter provided a reading for the fluid flow rate $Q$. The pressure sensors were amplified (HX711, SparkFun Electronics, CO, USA) and connected along with the flow meter to a microcontroller (Nano Every, Arduino, Italy), and data were acquired in MATLAB (V. 2022A, MathWorks, MA, USA) using custom-built software (available upon request).

2.3. Basic junction notation

Before we proceed with experimental observations of our fluidic devices, a basic junction labelling notation must be introduced to avoid confusion about the experimental configurations. To this end, we introduce a systematic labelling technique inspired by the Dowker–Thistlethwaite notation used in the mathematical knot literature (Dowker & Thistlethwaite Reference Dowker and Thistlethwaite1983; Adams Reference Adams1994). Briefly, our device encompasses two connected intersecting microchannels separated by a thin PDMS membrane. Fluid pressure causes one microchannel to deflect into the other. In our notation, we label the expanding channel by an odd integer and the contracted channel by an even integer. This study focuses on the hydraulic fingerprint of different sequential channel self-intersections. To this end, we will suppose that our microfluidic devices can be made with identical characteristics (i.e. identical internal channel dimensions and perfectly aligned and centred channels). For each identical device unit, we label the expanding and contracted channel, which, for three devices, yields unit one with channels [1] and [2], unit two with channels [3] and [4] and unit three with channels [5] and [6] (figure 3a). Later (§ 5), we will consider non-identical unit devices, but for now, we will label channels in a non-identical unit device by letters [A] and [B] (figure 3a).

Figure 3. Tabulation of junction connections. (a) Labelling of channels in identical and non-identical unit devices. Tabulation of (b) serial, (c) nested and (d) mixed sequences. (e) Commutative rule for identical devices and ( f) examples of non-commutative sequences.

Having labelled each channel, tabulating a connection between one or more unit devices is relatively straightforward, and the sequences we study can be roughly divided into three categories. The first category, serial sequences, is arguably the simplest and encompasses the connection of units where each unit's channels are individually connected. For instance, for unit one, the connection between channels [1] and [2] are made to yield the serial connection [12] (figure 3b). Connecting [12] to another individually connected unit (sequence [34]) yields another serial sequence [1234] (figure 3b), and similarly for the sequence [123456] by connecting an additional unit device. However, the first channel [1] could also be connected to [3] instead of looping directly to [2]. This yields the category of nested sequences, such as [1342] and [135642], where one or more intersections are nested inside each other (figure 3c). The third category, mixed sequences, encompasses serial and nested configurations that are connected to or within each other, such as [134562] and [123564] (figure 3d). For identically produced unit devices, the action of replacing or switching one unit with another does not change the sequence (figure 3e), nor does reversing the direction at which pressure is applied (e.g. by switching the sequence [12] into [21], see figure 3e). This commutative rule is not applicable for non-identical devices (sequence [AB]) where the flow rate depends on which direction the pressure is applied (figure 3f). For instance, a non-identical device can be made by stacking one rounded channel [A] and a rectangular channel [B] in a unit device (we will explore an application related to flow rectification in § 5 using a non-identical device). Moreover, changing a connection in a sequence breaks the commutative rule (e.g. [135642] differs from [135624], figure 3f).

It is worth pointing out that the number of possible sequences increases dramatically with the number of connected identical unit devices. Only one sequence, [12], is permitted for a single unit device, while for two unit devices, three sequences, [1234], [1342] and [1324], are permitted (starting with channel [1]). For three and four unit devices, $15$ and $105$, respectively, unique sequences can be made. To keep our experiments manageable, we limit ourselves to a sample of these possible sequences representing the three sequence categories. These are listed in table 1.

Table 1. Table of hydraulic knot configurations we test. The knot notation is introduced in § 2.3, and our experiments sample configurations from the three categories: serial, nested and mixed configurations.

2.4. Observations on elementary intertwined configurations

Having outlined the experimental methods and protocol, we now focus on measured flow–pressure relationships for elementary configurations of our fluidic devices. To explore elastohydrodynamic effects in self-intersecting channels, we will first examine the simplest case: the serial connection [12], which involves a single unit device. To understand the impact of membrane deformations, we first measured the characteristics of the [12] sequence with rounded channels when the membrane was relatively thick ($\tau \approx 3\ \mathrm {mm}$), inhibiting significant deformations. For this device, we found that the $Q-\Delta p$ relationship was approximately linear, in accordance with the Hagen–Poiseuille law (figure 4a). In contrast, when the membrane was thin, a deviation from the 1 : 1 relationship was observed (figure 4a). When the pressure exceeded $\Delta p\approx 10\ \mathrm {kPa}$, which was necessary for significant membrane deformation, the flow rate became approximately constant for the conforming device (flow limiting regime). A deviation from the 1 : 1 relationship was also observed for the device with rectangular channels. However, a constant flow rate was not reached within the pressure range available in the experiment. For the conforming channel devices, we denote the flow rate plateau value by $Q_{max}$. Below the actuation pressure, the $Q-\Delta p$ relationship was approximately linear, with a slope equal to that of the thick membrane experiment. With two unit devices, we access the serial configuration [1234] and two nested configurations [1324] and [1342]. For the serial configuration, [1234], we find that $Q_{max}$ is similar to that of the [12] sequence. However, the actuation pressure is $\Delta p\approx 20\ \mathrm {kPa}$, roughly double that of the [12] sequence (figure 4b). Interestingly, $Q_{max}$ values for the two nested sequences are approximately equal and lower than that of the serial sequence (figure 4b).

Figure 4. Flow rate vs pressure ($Q-\Delta p$) relationships for elementary hydraulic knots. The knot configurations are schematically drawn in each panel. (a) The serial sequence [12] using a single unit device with rectangular channel shapes and a relatively thick membrane (black data points), rectangular channel shapes and a relatively thin membrane (grey data points) and conforming channel shapes and a relatively thin membrane (purple data points). Both conforming (purple) and rectangular (grey) junction channels were tested. In (b–d), all channels have conforming cross-sections. (b–d) Show experimentally measured characteristics of two, three and four unit device sequences connected according to the diagrams in the legends.

Similar qualitative patterns are seen for the three and four unit device experiments (figure 4c,d). The largest $Q_{max}$ is found for the serial sequences, [123456] and [12345678], while the lowest are found for the nested sequences, [135642] and [13578642]. The mixed sequences, e.g. [123564], have a $Q_{max}$ that lies between that of the serial and nested sequences (figure 4c), and similarly for the four unit device mixed sequences. A pattern also emerges for $Q_{max}$ for nested sequences; the more nested devices, the lower the flow rate plateau, e.g. [1342] has a larger $Q_{max}$ than [135642], that again has a higher $Q_{max}$ than [13578642].

We briefly summarize our experimental findings with the following qualitative statements. Serial sequences yield approximately the same $Q_{max}$, and the pressure required for reaching $Q_{max}$ scales approximately linear with the number of unit devices. Nested sequences yield the lowest $Q_{max}$, decreasing with additional nested unit devices. Mixed sequences have a $Q_{max}$ higher than nested and lower than serial sequences. Finally, the initial slopes of the $Q-\Delta p$ diagrams are independent of sequence configurations and depend only on the number of unit devices.

3.1. Single junction

3.1.1. General considerations and approximations

We consider a single junction, comprising a top and a bottom channel of width $w$, crossing perpendicularly and separated by a flexible square membrane of thickness $\tau$. As a result of the transmural pressure, the membrane bends downwards, reducing the bottom channel cross-section and increasing its hydraulic resistance. We assume that the membrane is thin, in the sense that $\tau \ll w$, and that it is clamped at its boundaries. In practice, the transmural pressure drop changes as the pressure drops along the compressed channel and is thus a function of $z$ (figure 5a). Furthermore, the membrane is clamped on all four edges in the $w$-by-$w$ area in which the two channels intersect, which demands a complicated two-dimensional description of the deflection in the form of cumbersome series expansions (Timoshenko & Woinowsky-Krieger Reference Timoshenko and Woinowsky-Krieger1987). Even if we accounted for such solutions, solving the fluidic problem, which amounts to solving Reynolds equation $\boldsymbol {\nabla }\boldsymbol {\cdot }(H^3 \boldsymbol {\nabla } p) = 0$ (with $H$ the channel height, see figure 5) if the membrane slope is small (Bruus Reference Bruus2007), does not lead to analytical solutions because of the complicated form of $H$. To keep the simplest possible approach yet retain the constriction of the bottom channel as the essential ingredient, we ignore the dependence along the $z$-direction by assuming uniform transmural pressure $\Delta p_m$ and neglect end effects near the clamped region at $z=\pm w/2$. Moreover, we assume that the membrane deflection $\zeta (x)$ remains small, such that $|\zeta '(x)|\ll 1$. However, the deflection is often comparable to, or larger than, the membrane thickness. Hence, the bending of the membrane induces a tension force $T$, which must be retained in the analysis (Landau & Lifshitz Reference Landau and Lifshitz1986). In such a case, the deflection obeys the equation

(3.1)\begin{equation} B \frac{\mathrm{d}^4 \zeta}{\mathrm{d}\kern0.7pt x^4} - T \frac{\mathrm{d}^2 \zeta}{\mathrm{d}\kern0.7pt x^2} = \Delta p_m , \end{equation}

where $B=E\tau ^3 /[12(1-\nu ^2)]$ is the flexural rigidity (with $\nu$ the Poisson ratio), with clamped boundary conditions

(3.2)\begin{equation} \zeta = 0,\quad\frac{\mathrm{d}\zeta}{\mathrm{d}\kern0.7pt x} = 0 \quad \mathrm{at\ }x={\pm} \frac{w}{2}. \end{equation}

The solution to (3.1) with boundary conditions (3.2) is

(3.3)\begin{equation} \zeta(x) = \frac{\Delta p_m w^4}{B}\left( -\frac{x^2}{2w^2\beta^2}+c_1\cosh\frac{\beta x}{w} + c_2\right), \end{equation}

where $\beta = w\sqrt {T/B}$ is a parameter comparing stretching with bending, and

(3.4a,b)\begin{equation} c_1 = \frac{1}{2\beta^3\sinh(\beta/2)} \quad \mathrm{and} \quad c_2 = \frac{1}{8\beta^2}-\frac{1}{2\beta^3}\coth \frac{\beta}{2}, \end{equation}

are integration coefficients. When $\beta \ll 1$, bending dominates, and the deflection $\zeta (x) \simeq \Delta p_m (x-w/2)^2(x+w/2)^2/(96B)$ follows a quartic function. When $\beta \gg 1$, stretching dominates, and the deflection $\zeta (x) \simeq \Delta p_m (w^2/4 - x^2)/(2T)$ follows a parabolic function (except for a thin boundary layer close to the clamped ends). The tension force is related to the extension of the clamped membrane upon deflection. In the limit $|\zeta '(x)|\ll 1$, it equals (Landau & Lifshitz Reference Landau and Lifshitz1986)

(3.5)\begin{equation} T = \frac{E\tau}{2w} \int_{{-}w/2}^{w/2} \left( \frac{\mathrm{d}\zeta}{\mathrm{d}\kern0.7pt x} \right)^2 \mathrm{d}\kern0.7pt x , \end{equation}

which can be computed from (3.3).

Figure 5. (a) Top view sketch of the junction, highlighted in light shade. The top channel is in the $z$-direction, while the bottom channel is in the $x$-direction. (b) Side view sketch of the rectangular channel of width $w$ and height $h_0$. (c) Side view sketch of the conforming junction of width $w$ and centre height $h_0$.

We shall henceforth assume $\beta \gg 1$, an assumption which we will shortly justify. In this case, (3.5) yields $T \simeq E\tau \Delta p_m^2 w^6/(24B^2 k^4)$, whence

(3.6)\begin{equation} \beta = 3^{1/3}\sqrt{2(1-\nu^2)}\left(\frac{\Delta p_m}{E}\right)^{1/3}\left(\frac{w}{\tau}\right)^{4/3}. \end{equation}

For our geometric and material parameters (see § 2.1), and setting $\Delta p_m = \Delta p_I = 2.0\ \mathrm {kPa}$, the parameter is $\beta \approx 15$, indicating that stretching is the dominant mode of deflection in our range of pressures. Both bending and stretching could be important to describe the deflection at relatively low pressure accurately. However, according to our data (see, e.g. figure 4a), the flow–pressure relationships of our devices are mostly linear at relatively low pressure, indicating that the precise membrane deformation is not important at low pressure. To keep the approach simple, we therefore restrict the following analysis to the stretching-dominated mode of deflection, where the membrane deflection can be written as

(3.7)\begin{equation} \zeta(x) = \frac{4h_0}{w^2}\left(\frac{\Delta p_m}{\Delta p_I}\right)^{1/3}\left(\frac{w}{2} - x \right)\left(x+\frac{w}{2}\right), \end{equation}

where we have fitted the value of $\beta$ (3.6) to the rounded channel micrograph in figure 2( f) such that deflection equals the channel centre height, $\zeta (0) = h_0$, when the transmural pressure equals the inflation pressure used to generate the channel shape, $\Delta p_m=\Delta p_I=2.0\ \mathrm {kPa}$.

We denote $H(x)$ the height of the constricted bottom channel, taken independent of $z$, consistently with the aforementioned assumption; hence, the flow is unidirectional in the $z$-direction, and the pressure gradient $\boldsymbol {\nabla } P$ along $z$ is uniform. Once again, this is not true in our square junctions, but it is consistent with our very simple approximation of disregarding any $z$-dependence along the junction. We assume that the height slowly varies, $|H'(x)| \ll 1$, and that the Reynolds number is small. In this case, the relation between the flow rate $Q$ and the pressure gradient takes the following form, see e.g. Christov et al. (Reference Christov, Cognet, Shidhore and Stone2018):

(3.8)\begin{equation} Q ={-}\frac{\boldsymbol{\nabla} P}{12\eta}\int_{{-}w/2}^{w/2} H(x)^3\,\mathrm{d}\kern0.7pt x. \end{equation}

It should be noted that the channel height $H(x)$ can also vary with the flow-direction coordinate, $z$, when the transmural pressure drop is non-uniform. However, to keep the modelling approach simple, we have neglected such dependency. In our experiments, the effective Reynolds number (Christov et al. Reference Christov, Cognet, Shidhore and Stone2018) is $Re'=(h_0/L)(\rho Q/(\eta w))\approx 0.2$. If the membrane deflection and the channel height were truly independent of $z$, we could replace the pressure gradient $\boldsymbol {\nabla } P$ in (3.8) by $\Delta P_b/w$, with $\Delta P_b$ the pressure drop across the constricted bottom channel, and $w$ the length of the square junction. To account for the complicated three-dimensional structure of the membrane deflection and of the flow profile, we simply replace the relation $\boldsymbol {\nabla } P = \Delta P_b/w$ by

(3.9)\begin{equation} \boldsymbol{\nabla} P = \frac{\Delta P_b}{L_{eff}} , \end{equation}

where $L_{eff}$ is an effective junction length, which serves as a fitting parameter. Since the boundaries located at $z = \pm w/2$ tend to rigidify the membrane, the real deflection is overestimated by (3.7), and the bottom channel is less constricted. Hence, we expect $L_{eff}$ to be lower than $w$. Finally, inserting (3.9) in (3.8), the relation between the flow rate and the pressure drop across the constricted bottom channel then writes

(3.10)\begin{equation} Q ={-}\frac{\Delta P_b}{R_j}, \end{equation}

with the hydraulic resistance of the constricted bottom channel given by

(3.11)\begin{equation} R_j^{{-}1} = \frac{1}{12\eta L_{eff}} \int_{{-}w/2}^{w/2} H(x)^3 \,\mathrm{d}\kern0.7pt x . \end{equation}

The simplifications leading to the hydraulic resistance of the constricted bottom channel were also exploited by Ozsun, Yakhot & Ekinci (Reference Ozsun, Yakhot and Ekinci2013) to estimate the nonlinear resistance of deformed channels. We shall now calculate explicitly the flow rate–pressure drop relationship in the two cases encountered in our experiments: (i) a conforming junction, and (ii) a rectangular junction.

3.1.2. The conforming junction

The majority of our experiments are conducted on channels with a rounded cross-sectional shape (figure 2f), moulded from a membrane inflated by the pressure $\Delta p_I$, as discussed in § 2.1. The channel height follows approximately the parabolic shape $h(x)=(4h_0/w^2)(w/2-x)(w/2+x)$. As alluded to in the previous section, the compressed channel height then takes the parabolic form (figure 5b)

(3.12) \begin{equation} H(x) = \frac{4h_0}{w^2}\biggl[1-\left(\frac{\Delta p_m}{\Delta p_I}\right)^{1/3}\biggr]\left(\frac{w}{2}-x\right)\left(x+\frac{w}{2}\right). \end{equation}

Hence, the junction is conforming in the sense that, if contact were to be established between the membrane and the bottom wall of the channel, it would extend over all the channel, thereby closing it. Inserting (3.12) in (3.11) yields

(3.13) \begin{equation} R_j^{{-}1} = \frac{1}{R_1}\biggl[1-\left(\frac{\Delta p_m}{\Delta p_I}\right)^{1/3}\biggr]^3, \end{equation}

where

(3.14)\begin{equation} R_1 = \frac{105}{4} \frac{\eta L_{eff}}{wh_0^3}, \end{equation}

is the undeformed resistance of the compressed channel segment. The resistance of the straight segments of the rounded channel device, which are not influenced by membrane deflections, is given by

(3.15)\begin{equation} R_0 = \frac{105}{4} \frac{\eta L}{wh_0^3}. \end{equation}

3.1.3. The rectangular junction

Some of our devices present channels with standard rectangular cross-sections of uniform $h_0$ when the membrane is undeformed (figure 5c). When the membrane is deformed, two cases need to be analysed, according to whether the membrane makes or not contact with the bottom wall. This has been discussed by Gilet & van Loo (Reference Gilet and van Loo2022) but in the case of a membrane of constant curvature. The case without contact, respectively with contact, corresponds respectively to $\zeta (0) < h_0$ and $\zeta (0) \geq h_0$, where $\zeta (0)$ is computed from (3.7); hence, this corresponds respectively to $\Delta p_m < \Delta p_I$ and $\Delta p_m \geq \Delta p_I$.

In the absence of contact, the height profile is

(3.16) \begin{equation} H(x) = h_0 - \zeta(x) = h_0 \biggl[ 1 - 4 \left( \frac{\Delta p_m}{\Delta p_I} \right)^{1/3} \left( \frac{1}{2} - \frac{x}{w} \right) \left( \frac{1}{2} + \frac{x}{w} \right) \biggr]. \end{equation}

Inserting this expression in (3.11) yields:

(3.17) \begin{equation} R_j^{{-}1} = \frac{wh_0^3}{12\eta L_{eff}} \biggl[ 1 - 2 \left( \frac{\Delta p_m}{\Delta p_I} \right)^{1/3} + \frac{8}{5} \left( \frac{\Delta p_m}{\Delta p_I} \right)^{2/3} - \frac{16}{35} \frac{\Delta p_m}{\Delta p_I} \biggr]. \end{equation}

In the presence of contact, the expression (3.3) for the membrane deflection no longer holds. However, neglecting the bending term in (3.1) shows that the membrane profile remains parabolic away from the contact area. Moreover, the previous analysis showed that the curvature is given by $\mathrm {d}^2 \zeta /\mathrm {d}\kern0.7pt x^2 = -8h_0 (\Delta p_m/\Delta p_I)^{1/3}/w^2$, see (3.3). Hence, at given $\Delta p_m > \Delta p_I$, the portions of curved membrane outside of the contact area are parabolas of curvature $-8h_0 (\Delta p_m/\Delta p_I)^{1/3}/w^2$, tangent to the channel bottom wall at the (yet unknown) edge $x = \pm w_c/2$ of the contact area, and such that $\zeta = 0$ at $x = \pm w/2$. Straightforward algebra shows these conditions select both the contact edge location, $w_c = w[1 - (\Delta p_I/\Delta p_m)^{1/6}]$, and the following parabolic profile in the range $x\in [w_c/2,w/2]$:

(3.18)\begin{equation} \zeta(x) ={-}\frac{4h_0}{w^2}\left(\frac{\Delta p_m}{\Delta p_I}\right)^{1/3} \left( x^2 - w_c x + \frac{1}{2} w_c w - \frac{1}{4} w^2 \right) , \end{equation}

and a symmetric profile in the range $x\in [{-w/2},{-w_c/2}]$. Expressing the height channel $H(x) = h_0 - \zeta (x)$ and inserting in (3.11) then yields

(3.19)\begin{equation} R_j^{{-}1} = \frac{wh_0^3}{12\eta L_{eff}} \times \frac{1}{7} \left( \frac{\Delta p_I}{\Delta p_m} \right)^{1/6} . \end{equation}

The hitherto derived expressions for the hydraulic resistance of both the conforming and rectangular junctions share the same form

(3.20)\begin{equation} R_j = R_1 f\left(\frac{\Delta p_m}{\Delta p_I}\right), \end{equation}

where $R_1$ is the undeformed resistance. The crucial point is that they are functions of the transmural pressure (through the dimensionless function $f$), which depends on the rest of the hydraulic network. We now study the entire network to quantify the impact of junction peculiarities on flow rate–pressure drop correlation.

3.2. The simple overlap

3.2.1. General conditions

We now turn to the simple network involving a junction: the simple overlap (i.e. the single serial junction [12]). To connect channels, we fitted each straight device channel (with resistance $R_0$) with a cylindrical connecting tube. The connector has resistance $R_c = 8\eta L_c/({\rm \pi} a^4)$ (where $a$ and $L_c$ are the radius and length, respectively, with dimensions given in § 2.1). With our parameters $R_c \approx 0.5 R_0$, and thus the connecting tubes are not negligible in the total device resistance. When two channels are connected, e.g. by connecting [1] to [2] in the [12] knot (see figure 6) we are left with four segments of $R_0$ resistance, four segments of $R_c$ resistance and the junction resistance $R_j$. Following the standard additivity law of hydraulic resistors connected serially (Bruus Reference Bruus2007), the relation between applied pressure $\Delta p$ and the flow rate $Q$ is therefore

(3.21)\begin{equation} \Delta p = (4R_s + R_j)Q, \end{equation}

where $R_s = R_0 + R_c$ is the resistance of a straight channel segment independent of membrane deformations. Moreover, the transmural pressure across the membrane comes from the pressure drop along the path from the expanding top channel's junction outlet to the inlet of the compressed bottom junction channel, and hence it obeys $\Delta p_m = 2R_sQ$. Therefore, from (3.20), $\Delta p = [4R_s + R_1f(\bar {Q})]Q$ with a dimensionless flow rate $\bar {Q}=Q/Q_{max}$ with

(3.22)\begin{equation} Q_{max} = \frac{\Delta p_I}{2R_s}, \end{equation}

where the factor of $2$ arises from two resistive segments of $R_s = R_0 + R_c$ linking the bottom junction outlet to the top junction inlet. We can put the characteristics in dimensionless form by setting $\bar{\Delta p} = \Delta p/\Delta p_I$ and introducing the relative resistance $\alpha = R_1/(4R_s)$, and we get

(3.23)\begin{equation} \bar{\Delta p} = 2[1+\alpha f(\bar{Q})]\bar{Q}. \end{equation}

Notice that this relationship is nonlinear because of the coupling between the flow and the membrane deformation, quantified by the function $f$.

Figure 6. Schematic drawings of (a) the [12] and (b) 1342 junctions. Straight channel segments, connective tubes and junction resistances are termed $R_0$, $R_c$ and $R_j$, respectively. When multiple junctions are connected, e.g. in [1342] in b), the junction resistances are labelled with superscripts according to the device channels, see figure 3.

3.2.2. The conforming junction

In the case of the conforming junction, we get from (3.13) and (3.20) that

(3.24)\begin{equation} f(\bar{Q}) = (1-\bar{Q}^{1/3})^{{-}3}, \end{equation}

which when injected into (3.23) yields

(3.25)\begin{equation} \bar{\Delta p} = 2[1+\alpha(1-\bar{Q}^{1/3})^{{-}3}]\bar{Q}. \end{equation}

The natural outcome of this calculation is the prediction of a maximal possible flow rate, equal to $Q_{max}$. Here, flow limitation comes from the fact that at high pressure drop, any increase of the flow velocity induced by an increase of the pressure drop is exactly compensated by an increase of the hydraulic resistance of the occluding junction, leading to a saturation of the flow rate.

3.2.3. The rectangular junction

In the case of the rectangular junction, we get from (3.17), (3.19) and (3.20) that

(3.26) \begin{align} f(\bar{Q}) = \left\{ \begin{array}{@{}ll} \left( 1 - 2\bar{Q}^{1/3} + \dfrac{8}{5} \bar{Q}^{2/3} - \dfrac{16}{35} \bar{Q} \right)^{{-}1} & \mathrm{if} \ \bar{Q} \leq 1 \\ 7\bar{Q}^{1/6} & \mathrm{if} \ \bar{Q} \geq 1 \end{array} \right., \end{align}

which, injected into (3.23), yields the characteristics for rectangular channels. In marked contrast with the conforming junction, the flow rate can increase indefinitely if the pressure drop increases, according to the scaling law $\bar{\Delta p} \propto \bar {Q}^{7/6}$, which is a mild nonlinearity.

3.3. Multiple junctions

Multiple junction devices can be connected in various sequences. In our experiments, we found that the saturation flow rate $Q_{max}$ depends on the particular sequential connection of junctions, justifying the development of a mathematical model that can predict a given sequence's $Q-\Delta p$ relationship. What makes this exercise particularly challenging is that the transmural pressure drop, $\Delta p_m$, across one channel intersection can depend on the pressure-drop-dependent resistance of other intersections. To proceed, we will start by recognizing that the hydraulic knots we study are simple serial connections of resistances. Therefore, similarly to the simple junction, we can use the additivity law of hydraulic resistors connected in series (Bruus Reference Bruus2007) to write the relationship between applied pressure and output flow rate

(3.27) \begin{equation} \Delta p = Q\biggl[4NR_s + \sum_{{unit\ devices}} R_j^{{[ik]}}\biggl(\frac{\Delta p_m^{{[ik]}}}{\Delta p_I}\biggr)\biggr], \end{equation}

where $N$ is the number of unit devices in the knot. The first term stems from each device having four straight segments (independent of membrane deformations) of resistance $R_s$. The second term is the sum of junction resistances. Notice that we assigned the superscript ${[ik]}$ to junction resistance and transmural pressure, where, consistently with § 2.3, the odd integer $i$ corresponds to the expanded channel, and the even integer $k$ to the contracted channel, of a given unit device. For instance, for two unit devices [12] and [34], the second term in (3.27) yields $R_j^{{[12]}}+ R_j^{{[34]}}$, where the junction resistances depend on the transmural pressures $\Delta p_m^{{[12]}}$ and $\Delta p_m^{{[34]}}$, respectively (see figure 6b). As alluded to, evaluating (3.27) can be challenging because the junction resistances can couple. To help clarify, we shall first consider an example of a nested sequence, [1342].

In sequence [1342] the junction [34] is nested within [12] (figure 6b), and it is clear that the transmural pressure drop of [12], $\Delta p_m^{{[12]}}$ depends on the junction resistance of [34], $R_j^{{[34]}}$, which depends on the transmural pressure drop of [34], $\Delta p_m^{{[34]}}$. However, what enables us to proceed with developing the analytical model is that the transmural pressure across [34], $\Delta p_m^{{[34]}}$, does not depend on other junction resistances. In fact, for any sequence that follows our notation, there is always at least one junction whose $\Delta p_m$ does not depend on other junctions. For the example [1342], this means that $\Delta p_m^{{[34]}} = 2R_sQ$ (similarly to the simple overlap in § 3.2.1), which can be readily injected into (3.20) to find the junction resistance $R_j^{{[34]}}$. Knowing $R_j^{{[34]}}$ enables evaluating the pressure drop across the [12] junction as $\Delta p_m^{{[12]}} = Q(6R_s+R_j^{{[34]}})$, where the additional factors of $R_s$ come from the extra straight channel sections in the unit containing [34]. Again, $\Delta p_m^{{[12]}}$ can be injected into (3.20) to yield the junction resistance $R_j^{{[12]}}$ allowing evaluating the $\Delta p-Q$ relationship by summing all the resistive terms in (3.27).

Generally, the challenging step in evaluating (3.27) for a given sequence is identifying which junction to start with in terms of calculating the first junction-independent $\Delta p_m$ and $R_j$. However, using our knot notation, it is simply done by taking the first even integer number in a knot sequence. The next even integer junction may depend on the first or be trivial (which is true for all serial configurations). For example, in the sequence [135642], we recognize that [6] is the first even number, and we calculate $\Delta p_m^{{[56]}}=2R_sQ$ and $R_j^{{[56]}}$ (from (3.20)) which is then injected into $\Delta p_m$ and $R_j$ first for [34] and then [12], finally allowing evaluation of (3.27). In [135642] junction [34] depends on [56], and [12] depends on both [34] and [56]. In contrast, the serial connection [123456] contains only trivially solved junctions, as all of the transmural pressure drops can be written $\Delta p_m = 2R_sQ$.

We end this section by discussing the maximum attained flow rate in our hydraulic knot configurations. As was evident from our data (figure 3), nested sequences yield a lower $Q_{max}$ than serial sequences. In (3.25) we found that, for a single rounded junction, the maximum flow rate $Q_{max} = \Delta p_I/(2R_s)$ is set by the inflation pressure (for our experiments $\Delta p_I = 2.0\ \mathrm {kPa}$), as a result of the dimensionless flow–structure coupling function $f(\bar {Q})$. The implications for $Q_{max}$ in our knots are best appreciated by returning to the example of the [1342] sequence. The junction resistance $R_j^{{[34]}}$ is evaluated by inserting $\Delta p_m^{{[34]}} = 2R_sQ$ into (3.20), yielding

(3.28) \begin{equation} R_j^{{(34)}} = R_1 \biggl[1-\left(\frac{2R_sQ}{\Delta p_I}\right)^{1/3}\biggr]^{{-}3}, \end{equation}

which predicts the maximum flow rate $Q_{max} = \Delta p_I/(2R_s)$. To evaluate $R_j^{{[12]}}$, we insert $\Delta p_m^{{[12]}}=Q(6R_s+R_j^{{[34]}})$ into (3.20), which gives

(3.29) \begin{align} R_j^{{[12]}} &= R_1\biggl[1-\biggl(\frac{Q(6R_s+R_j^{{[34]}}}{\Delta p_I}\biggr)^{1/3}\biggr]^{{-}3} \nonumber\\ &= R_1 \biggl\{1-\left(\frac{Q}{\Delta p_I}\right)^{1/3}\biggl[6R_s+R_1\biggl(1-\biggl(\frac{2R_sQ}{\Delta p_I}\biggr)^{1/3}\biggr)^{{-}3} \biggr]^{1/3}\biggr\}^{{-}3}, \end{align}

where the reduced maximum flow rate $Q_{max}'$ can be found by equating the braced term to zero and solving for $Q$. For the geometric and material parameters used in our experiments, $Q_{max}' \approx 0.24 Q_{max}$, which is in reasonable agreement with our data comparing the maximum flow rates of [12] with [1342] in figure 4(a,b).

Before moving on with further comparing our mathematical model with our experimental data, we will briefly discuss modelling approaches for non-identical unit devices connected in a hydraulic knot.

3.4. Considerations on multiple junctions with non-identical geometries

In the previous section, we discussed modelling multiple connected identical junctions. In general, the junctions need not be identical and may differ in, e.g. channel height (linked to the inflation pressure for conforming channels), width or length. This can result in device-specific resistances, inflation pressures and different flow–structure coupling functions (i.e. for different channel geometries). For $N$ unit devices, the $\Delta p -Q$ relationship in (3.27) can be adapted to

(3.30) \begin{equation} \Delta p = Q\biggl[\sum_{j=1}^{2N} R_s^{{[j]}} + \sum_{{unit\ devices}} R_j^{{[ik]}} \biggl(\frac{\Delta p_m^{{[ik]}}}{\Delta p_I^{{[ik]}}}\biggr)\biggr], \end{equation}

where the first term sums the (different) straight channel resistances $R_s^{{[j]}}$ and the second term sums the junction resistances. We will not go further into details on modelling multiple non-identical junctions, but note that we have an experiment in § 5 where we explore the flow-rectification features of a device that has two different channel geometries.

5.1. Smoothing the flow rate output from a peristaltic pump

A significant drawback of peristaltic pumps is the ripples in flow rate output that arise due to the rotational compression of the pump tubing. This unsteady flow can be unfavourable for, e.g. flow cytometry (Piyasena & Graves Reference Piyasena and Graves2014) and infusion systems (Snijder et al. Reference Snijder, Konings, Lucas, Egberts and Timmerman2015), and considerable care must be taken in smoothing the flow rate. We will consider the basic application set-up where a peristaltic pump is connected to a device with resistance $R_D$. Introducing a compliant vessel between the pump and device can attenuate some, but not all, ripples in the peristaltic flow (Kang & Yang Reference Kang and Yang2012; Kang et al. Reference Kang, Yeom, Seo and Lee2014). Further improvements can be made by including a compliant vessel and a passive nonlinear resistance that only permits one specific flow rate within a range of pressure, akin to the flow limitation exhibited by our hydraulic knots (Doh & Cho Reference Doh and Cho2009; Zhang et al. Reference Zhang, Xiang, Tang, Huang, Wang, Yi and Ni2015; Biviano et al. Reference Biviano, Paludan, Christensen, Østergaard and Jensen2022). However, the introduced nonlinear resistances need to be carefully tuned such that the flow rate limitation level matches the application specifications. We hypothesize that we can overcome this obstacle by exploiting a feature of our hydraulic knots: combining multiple identical devices (each with the same flow limiting level) in different configurations can yield a spectrum of flow limiting levels (figure 4). To test our hypothesis, we conducted an experiment that comprised a standard peristaltic pump (WPX1 with four rollers and $0.50\ \mathrm {mm}$ internal diameter tubing, Welco, Japan). We connected one end of the pump tubing to a water-filled reservoir and the other end to a 20 cm long thin-walled silicone tubing (inner diameter $2.0\ \mathrm {mm}$, outer diameter $3.0\ \mathrm {mm}$, GRA-HS0200005, Mikrolab Aarhus, Denmark). We then connected the silicone tubing to one of our hydraulic knots and then a 10 cm long rigid resistor ($0.50\ \mathrm {mm}$ internal diameter polyetheretherketone tubing, Mikrolab Aarhus, Denmark), which served as a proxy for the application device (following the Hagen–Poiseuille law, $R_D \approx 5.9{\cdot } 10^{10}\ \mathrm {Pa\,s\ m}^{-3}$ in our experiment). Using a flow meter, we measured the flow rate output from the resistive tube (see § 2). The pump was controlled via a stepper drive (A4988, Allegro MicroSystems, NH, USA), and the pumping frequency was kept constant in our experiments at $f\approx 1.4\ \mathrm {Hz}$.

To establish a baseline measurement of the peristaltic pump output, we conducted the experiment without installing any of our knots, i.e. with the pump connected to the compliant tube connected to the resistive tube. Next, we tested the effect of our nonlinear resistive knots by installing the nested configurations [12], [1342] and [135642] between the compliant tube and the resistive tube. The reason for using our nested knots was the observation that the flow rate limiting level decreases with an increased number of nested devices. We measured the flow rate output as a function of time for each configuration. To quantify the attenuation efficacy, we measure the ratio between peak-to-peak amplitude in flow rate ($A$) and average flow rate ($\bar {Q}$). We define the smoothing efficacy $\varepsilon = 1-A/\bar {Q}$, where, for smooth flow, $\varepsilon \simeq 1$, and for noisy flow, $\varepsilon \simeq 0$.

For our baseline measurement, using only a compliant tube to attenuate the flow rate fluctuations, we measure a noisy flow and a low smoothing efficacy $\varepsilon \approx 0.05$ (figure 9a). The smoothing efficacy is improved by including the [12] knot, where we find $\varepsilon \approx 0.60$. Note that, with the inclusion of the [12] knot, the mean flow rate remains almost the same as in the baseline case ($\bar {Q}\approx 2.51\ \mathrm {ml\ min}^{-1}$ at baseline and $\bar {Q}\approx 2.22\ \mathrm {ml\ min}^{-1}$ with the [12] knot). The [1342] and [135642] knots permit lower flow rates than the [12] knots (see, e.g. figure 4). When we include the [1342] knot in the peristaltic set-up, the mean flow rate is reduced to $\bar {Q}\approx 1.17\ \mathrm {ml\ min}^{-1}$ and the smoothing efficacy is improved to $\varepsilon \approx 0.88$. The smoothing efficacy is further improved with the [135642] knot where, remarkably, $\varepsilon \approx 0.96$. Here, the measured flow rate oscillations are of the order of the flow rate sensor's noise level. Our experiment confirms that incorporating nested hydraulic knots into a standard peristaltic pump set-up can significantly reduce pulsatile flow (to the extent of $\varepsilon \approx 0.96$) and regulate the average flow rate output by adjusting the knot sequence. For the [12] and [1342] knots it is not unlikely that the smoothing can be improved in response to, e.g. increasing the pumping frequency, as demonstrated by Biviano et al. (Reference Biviano, Paludan, Christensen, Østergaard and Jensen2022) using similar nonlinear hydraulic resistances.

Figure 9. Soft hydraulic knots can enable microfluidics applications. (a) Output flow rate from a peristaltic pump connected to a compliant tube, our hydraulic knot and, finally, a resistive tube. The black data points are a baseline with no hydraulic knots. The purple, orange and blue data points are for the [12], [1342] and [135642] knots, respectively. The different series are arbitrarily shifted in time for a better rendering. (b) Flow–pressure characteristics of a direction-dependent resistor comprising an intersection between a rectangular and rounded channel. Schematic drawings show (not to scale) the approximate membrane position in the intersection. (c) Flow vs time characteristics of the direction-dependent resistor when an oscillating pressure source (inset) is applied. When the pressure amplitude exceeds $\approx 20\ \mathrm {kPa}$, a net significant net volume displacement is observed during each oscillation cycle.

5.2. Valveless flow rectification in a direction-dependent resistor

We now turn our attention to another feature of our device: the ability to intentionally mismatch the geometries of the overlapping channels. This can potentially be used to create a simple direction-dependent resistance such that the flow rate output will depend on the device's orientation (i.e. swapping the inlet and the outlet). Such diode-like fluidic resistors have been used to develop low-Reynolds-number microfluidic devices for flow rectification (Groisman & Quake Reference Groisman and Quake2004; Liu et al. Reference Liu, Chen, Taylor, Scherer and Kartalov2009a). Note that, at higher Reynolds numbers, fluid inertia can be exploited to accomplish flow rectification (Tesla Reference Tesla1920; Nguyen, Abouezzi & Ristroph Reference Nguyen, Abouezzi and Ristroph2021a; Nguyen et al. Reference Nguyen, Oza, Abouezzi, Sun, Childress, Frederick and Ristroph2021b).

We aim to test whether a single junction comprising one of our rounded and one of our rectangular channels can rectify the flow (see § 2 for channel dimensions). We connect the rounded and rectangular channel to form a [12] junction, and we then apply pressure to channel [1] first (positive $\Delta p$) or channel [2] first (negative $\Delta p$). We characterized the device similarly to our other hydraulic knots and found that the $Q-\Delta p$ characteristics are indeed asymmetric (figure 9b). At negative $\Delta p$, the flow is limited at $Q_{max}^{(-)}\approx -6.8 \ \mathrm {ml\ min}^{-1}$, while at positive $\Delta p$, the flow slowly converges to a limit at $Q_{max}^{(+)}\approx 11.0\ \mathrm {ml\ min}^{-1}$. The asymmetry is caused by the membrane occluding a larger ratio of the conforming channel than the rectangular channel at a similar (absolute) applied pressure. We then performed a time-dependent experiment, using two pressure controllers to control the pressure on either end of the asymmetric device. We coded the controllers to generate a rectangular pulse of applied pressure between positive $\Delta p=(5,10,20,40)\ \mathrm {kPa}$ and negative $-\Delta p=-(5,10,20,40)\ \mathrm {kPa}$, keeping the time-averaged applied pressure $\langle \Delta p\rangle =0\ \mathrm {kPa}$ (see inset in figure 9c). The period of each pulse was $\approx 25\ \mathrm {s}$. For small pressure pulses, $\Delta p =5\unicode{x2013}10\ \mathrm {kPa}$, the magnitude of flow rate output was approximately independent of direction. However, when we increased the pressure to $\Delta p = 20\ \mathrm {kPa}$ and $\Delta p=40\ \mathrm {kPa}$, we observed that the flow output depended on direction and that a net volume of fluid was displaced in the positive pressure direction, consistent with the measured $Q-\Delta p$ relationship in figure 9(b). We quantify the yield by the relative net flow rate output $Y=(Q^+-|Q^{-}|)/|Q^{-}|$, which, for our applied pressure pulses, yields $Y\approx 3\,\%,$ $6\,\%$, $21\,\%$ and $63\,\%$ in our experiment. We thus see that the yield increases with increasing applied pressure. It is not unlikely that the yield can be tuned by combining several asymmetric devices in optimized hydraulic knot configurations.