1. Introduction
In a study by Boucher & Friot (Reference Boucher and Friot2017) for the International Union for Conservation of Nature, it was found that between
$15\,\%$
and
$31\,\%$
of all plastics in the ocean originate from household and industrial appliances, termed primary microplastics, and up to
$35\,\%$
of these particles are composed of microfibres from our clothes. The term ‘microplastic’ was introduced by Thompson et al. (Reference Thompson, Olsen, Mitchell, Davis, Rowland, John, McGonigle and Russell2004), who identified the danger of the build up of these micron-sized plastic particles in the oceans, causing a risk to marine life. Napper & Thompson (Reference Napper and Thompson2016) and De Falco et al. (Reference De Falco, Di Pace, Cocca and Avella2019) found that the microplastics removed from textiles take the form of microfibres, with diameters of
$11.9$
–
$17.7$
$\unicode{x03BC}$
m and lengths of
$360$
–
$660$
$\unicode{x03BC}$
m. Recent work on quantifying the number of microfibres that end up in the ocean after each wash has been reviewed by several authors, with a consensus of at least
$700\,000$
(Dris et al. Reference Dris, Imhof, Sanchez and Gasperi2015; Napper & Thompson Reference Napper and Thompson2016; Acharya et al. Reference Acharya, Shaida, Yang and Noureddine2021; Hazlehurst et al. Reference Hazlehurst, Tiffin, Sumner and Taylor2023).
It is becoming vital to remove these fibres at the source of the discharge, i.e. from washing machine wastewater. Some conventional washing machines incorporate dead-end filters which are effective at removing microfibres but clog relatively quickly (Enten et al. Reference Enten, Leipner, Bellavia, King and Sulchek2020; Akarsu, Kumbur & Kideys Reference Akarsu, Kumbur and Kideys2021). Enten et al. (Reference Enten, Leipner, Bellavia, King and Sulchek2020) also examine periodic back-flushing through the dead-end filter to wash fibres off the surface and reduce the need for filter cleaning. However, consumers often postpone cleaning such filters for as long as possible, which results in the fluid bypassing the filter through the emergency overflow mechanism, thus negating the purpose of the filter. An alternative method for increasing the lifespan of a washing machine dead-end filter, being explored by Beko PLC, involves diverting as much water away from the dead-end filter as possible whilst retaining as many microfibre particles as possible flowing into the dead-end filter (see figure 1).
Layout schematic of a branched channel filter preceding a dead-end filter. Microfibre particles, trajectories and foulant are indicated in red and water flow is indicated in blue. The operating directions are indicated by black arrows.
Beko’s prototype designs are motivated by a fluid–particle separation technique identified in manta ray fish, a type of suspension feeder, coining the term ricochet separation (Divi, Strother & Paig-Tran Reference Divi, Strother and Paig-Tran2018). Within the manta ray’s mouth, there is a gill-like pore structure resembling branched channels. Plankton-rich water flows over this structure, with clean water flowing through the pores while plankton collide with the rigid structure and ricochet back into the main free-stream flow above the pores, thus providing an efficient filtration mechanism. Divi et al. (Reference Divi, Strother and Paig-Tran2018) conclude that filtration efficiency increases for high-Reynolds-number flows, and that particles much smaller than the size of the pores may be filtered, with increasing efficiency for increasing size. They simulate a quasi-steady flow and model the particle dynamics using Newton’s second law of motion, and compare the results with experiments.
Beko believe that the ricochet method might be an excellent way for removal of a substantial amount of water whilst retaining microplastic particles, removing the challenge of constructing a device with pores sufficiently small to remove the particles via size exclusion. Their aim is to use ricochet separation prior to the dead-end filtration, to reduce the amount of water – but not microfibres – that flows into a dead-end filter – a trade-off between maximising the flow and minimising the number of particles through the branched channel structure.
The method of ricochet separation in manta rays, and various other mobula rays, has been further considered by Mao, Bischofberger & Hosoi (Reference Mao, Bischofberger and Hosoi2024). They consider a framework to describe the trade-off between fluid leakage rate and the critical particle size through the branched channels. For low Reynolds numbers, they derive an analytic solution for the leakage rate through each branched channel, whereas for larger Reynolds numbers, they rely on numerical simulations. In the latter regime, numerical simulations indicate that altering the branch angle while keeping the branch width constant does not substantially affect the leakage rate. In their model, as in a mobula ray’s mouth, they assume the pressure at all outlets is atmospheric. This differs from a scenario in which a dead-end filter is placed at the main outlet, as this will experience a pressure build-up. They provide a relationship between the flow leakage rate and particle cutoff size in both the low- and high-Reynolds-number cases, concluding that the particle cutoff size is smaller for smaller leakage rates. In a different approach, Hamann (Reference Hamann2023) designed and analysed ‘semi’-cross-flow filters, also motivated by suspension feeders, to explore whether a single filter is possible for microplastic capture without the need for a further dead-end filter – later proposing such a filter in Hamann et al. (Reference Hamann, Reuß, Herzog, Schreiber, Geitner and Blanke2025).
In this paper, motivated by the industrial idea of ricochet separation for filtration, we build and solve a mathematical model to understand the trade-off between leakage rate and particle size that avoids the need to perform computationally challenging numerical simulations in the high-Reynolds-number laminar regime. We will exploit the fact that there are many branched channels in the device and so the proportion of fluid that passes through each one is low, which will give rise to multi-scale behaviour.
To facilitate a mathematical approach, we consider a simpler geometry than Divi et al. (Reference Divi, Strother and Paig-Tran2018) and Mao et al. (Reference Mao, Bischofberger and Hosoi2024), removing their rounded-lobe structure and replacing this with a series of T-junction branched channels equidistantly spaced along a channel. We exploit the fact that there are many branched channels in the device and so the proportion of fluid that passes through each one is low, which gives rise to multi-scale behaviour. We use a multiple-scale analysis to derive an effective boundary condition that encapsulates the key behaviour in the device.
Dalwadi et al. (Reference Dalwadi, King, Dyson and Arkill2020) study the normal, but not tangential, flow in a single T-junction and they reduce the problem to that of a point sink. We use this idea to transform the branched channels into a series of point sinks of appropriate strengths. The idea of using homogenisation theory to derive effective boundary conditions is elegantly presented in Chapman, Hewett & Trefethen (Reference Chapman, Hewett and Trefethen2015), and further developed in Hewett & Hewitt (Reference Hewett and Hewitt2016), where effective boundary conditions for the electrical potential in a discrete metal cage are derived, which smooths out the effects caused by individual point sources. We use the ideas in these papers to solve our resulting cell problem using complex variable theory to conformally map and solve for the velocity potential, which we then match into the outer flow to give us our effective leakage boundary condition.
The concept of asymptotic homogenisation was formally introduced by Bensoussan, Lions & Papanicolaou (Reference Bensoussan, Lions and Papanicolaou1978). When using asymptotic expansions in a microscale cell structure, they were able to find effective equations to describe the macroscale behaviour, which incorporated microscale behaviour into the material parameters. This approach drastically simplifies the numerical complexity of problems. We take a similar approach in this paper; however, rather than deriving effective equations, we find effective boundary conditions, and so our method is instead a multi-scale and boundary layer analysis (Van Dyke Reference Van Dyke1975; Hinch Reference Hinch1991).
Such an approach has also been considered by Bruna, Chapman & Ramon (Reference Bruna, Chapman and Ramon2015) to derive an effective diffusive flux condition for the concentration through a thin-film composite membrane. Similarly, Bottaro & Naqvi (Reference Bottaro and Naqvi2020) have used homogenisation techniques to derive an effective boundary condition for purely tangential flow past a wall with periodic imperfections, with no fluid removal. Their result recovers well-understood phenomena such as Navier slip in their leading-order results. This is extended in the work by Naqvi (Reference Naqvi2021), where they further consider homogenisation of the flow over and through porous micro-structures and elastic surfaces. However, our set-up is different, since we have pressure gradients down our branched channels.
We begin in § 2 by deriving a high-Reynolds-number laminar flow model in a branched channel filter. We simplify this model in § 3 to an inviscid flow everywhere between two parallel fixed walls, with point sinks along the bottom wall; the sink strength is related to the flux found via Poiseuille flow in each branched channel. We solve for the outer flow far away from the wall in § 4 and the inner flow close to the wall in § 5; in § 6 we return to find the full explicit outer solution and further find composite solutions and an effective boundary condition in § 7. We compare these asymptotic results with the numerical simulations in § 8. In §§ 9 and 10, we consider a particle model, exploring how collisions with the wall take place under the influence of our asymptotic prediction for the flow. We model the motion of individual particles via a force balance model with a bouncing/removal wall condition to mimic ricochet or branched channel capture, respectively. For the application of our flow model in particle filtration, we explore the trade-off between the fraction of particles and the fluid flux fraction that passes through the branched channels. Finally, in § 11 we discuss our findings, draw conclusions and provide extensions.
2. Flow model derivation
We consider a two-dimensional domain of length
$L$
and inlet height
$h_1$
, with a series of
$N$
repeated T-junctions of length
$L_1$
, width
$h_2$
and spacing
$L_2$
, where
$N \boldsymbol{\cdot }L_2 = L$
, protruding from the bottom wall (see figure 2). Throughout this work, we use
$\hat {\boldsymbol{\cdot }}$
to indicate dimensional quantities. We use
$\hat {x}$
to denote distance along the domain from the left-hand side, and
$\hat {y}$
to denote distance above the bottom surface of the main channel. The centre of the top of each T-junction is located at
$(\hat {x}_i, 0)$
, where
\begin{equation} \hat {x}_i = \frac {(2 i - 1) L_2}{2}, \end{equation}
for
$i = 1, 2, {\cdots} , N$
, and we assume that the T-junctions are laterally symmetric about this centre point. We denote the full domain by
$\hat {\varOmega }$
, as the union of the main channel
$ ( \hat {x}, \hat {y} ) \in [0,L] \times [0,h_1]$
and branched channels
$ ( \hat {x}, \hat {y} ) \in [\hat {x}_i-({h_2}/{2}), \hat {x}_i +( {h_2}/{2})] \times [-L_1,0]$
. The domain,
$\hat {\varOmega }$
, has fixed boundary walls given by
$\partial \hat {\varOmega }_w$
, along with inlet at
$\hat {x}=0$
, main outlet at
$\hat {x}=L$
and outlets at the bottom of the
$N$
branched channels,
$\hat {y}=-L_1$
(illustrated in figure 2).
Two-dimensional repeatable T-junction domain,
$\hat {\varOmega }$
, given by a main channel compartment with
$N$
perpendicular branched channels on the bottom wall. Inlet and outlets are indicated by dashed black lines, the T-junction spacing is indicated by dashed red lines and boundary walls are denoted by
$\partial \hat {\varOmega }_w$
, in solid black lines. The domain design parameters are indicated as
$h_1$
,
$h_2$
,
$L$
,
$L_1$
,
$L_2$
and
$N$
.
We consider a high-Reynolds-number, laminar, steady-state flow, in which the fluid velocity
$\hat {\boldsymbol{u}} = (\hat {u}, \hat {v})$
and pressure
$\hat {p}$
satisfy the steady, incompressible Navier–Stokes equations
\begin{align} \hat {\boldsymbol{\nabla }} \boldsymbol{\cdot }\hat {\boldsymbol{u}} &= 0, \\[-12pt]\nonumber \end{align}
\begin{align} \rho _{\!f} (\hat {\boldsymbol{u}} \boldsymbol{\cdot }\hat {\boldsymbol{\nabla }} ) \hat {\boldsymbol{u}} &= - \hat {\boldsymbol{\nabla }} \hat {p} + \mu {\hat {{\nabla} }^2} \hat {\boldsymbol{u}}, \end{align}
where
$\rho _{\!f}$
is the fluid density and
$\mu$
is the viscosity. We impose a uniform plug inlet flow
\begin{align} \hat {u} = \frac {\mathcal{Q}}{h_1}, \quad &\text{on} \quad \hat {x} = 0, \quad \hat {y} \in [0, h_1], \\[-12pt]\nonumber \end{align}
\begin{align} \hat {v} = 0, \quad &\text{on} \quad \hat {x} = 0, \quad \hat {y} \in [0, h_1], \end{align}
such that
$\mathcal{Q}$
is the constant inlet flux at
$\hat {x}=0$
. We prescribe a no-slip and no-flux condition on the domain boundary,
$\partial \hat {\varOmega }_w$
, given by
\begin{equation} \hat {\boldsymbol{u}} = \boldsymbol{0}. \end{equation}
We apply an outlet pressure, given by
\begin{equation} \hat {p} = \hat {\mathcal{P}}_{\textit{out}} + p_{\textit{atm}}, \quad \text{on} \quad \hat {x} = L, \quad \hat {y}\in [0, h_1], \end{equation}
where
$\hat {\mathcal{P}}_{\textit{out}}$
is a given constant and
$p_{\textit{atm}}$
denotes constant atmospheric pressure. Finally, we impose that the pressure at the branch outlets is atmospheric and thus we write
\begin{equation} \hat {p} = p_{\textit{atm}}, \quad \text{on} \quad \hat {y} = - L_1, \quad \hat {x} \in \bigcup _{i=1}^N \left [ \hat {x}_i - \frac {h_2}{2}, \hat {x}_i + \frac {h_2}{2} \right ]\!. \end{equation}
2.1. Non-dimensionalisation
We non-dimensionalise the model, (2.2)–(2.8), using the high-Reynolds-number scalings
\begin{equation} \hat {\boldsymbol{x}} = L \boldsymbol{x}, \quad \hat {\boldsymbol{u}} = \frac {\mathcal{Q}}{L} \boldsymbol{u}, \quad \hat {p} = p_{\textit{atm}} + \frac {\rho _{\!f} \mathcal{Q}^2}{\epsilon ^{2} L^2} p, \end{equation}
where we have picked the pressure scaling so that the pressure at the end of the main outlet is
$O(1)$
, and where
\begin{equation} \epsilon = \frac {L_2}{L} \end{equation}
is defined as the dimensionless width of each T-junction, such that
$N \boldsymbol{\cdot }\epsilon = 1$
. We assume that the hole-width-to-T-junction-width aspect ratio,
\begin{equation} \delta = \frac {h_2}{L_2} = \frac {h_2 N}{L}, \end{equation}
is a fixed constant. In this case, if we increase the number of branches,
$N$
, then the width of each T-junction decreases. Hence each individual branched channel of width
$\delta \epsilon$
will also decrease with
$\epsilon$
and so a pressure increase of
$1/\epsilon ^2$
is required to force the same fluid flux through each hole. We define the Reynolds number, the dimensionless length of the branched channels and the dimensionless height of the main channel as
\begin{equation} \textit{Re} = \frac {\rho _{\!f} \mathcal{Q}}{\mu }, \qquad \lambda = \frac {L_1}{L}, \qquad \gamma = \frac {h_1}{L}, \end{equation}
respectively. We assume that in all cases,
$\textit{Re} \gg 1$
,
$\epsilon \ll 1$
,
$\gamma = {O}(1)$
,
$\lambda = {O}(1)$
and
$\delta \epsilon \ll \lambda$
.
The governing equations (2.2) and (2.3) become
\begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} &= 0, \\[-12pt]\nonumber \end{align}
\begin{align} (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla } ) \boldsymbol{u} &= - \frac {1}{\epsilon ^2} \boldsymbol{\nabla }\! p + \frac {1}{\textit{Re}} {\nabla} ^2 \boldsymbol{u} , \end{align}
and the boundary conditions (2.4)–(2.8) become
\begin{align} u = \frac {1}{\gamma }, \quad &\text{on} \quad x = 0, \quad y \in [0, \gamma ], \\[-12pt]\nonumber \end{align}
\begin{align} v = 0, \quad &\text{on} \quad x = 0, \quad y \in [0, \gamma ], \\[-12pt]\nonumber \end{align}
\begin{align} u = v = 0, \quad &\text{on} \quad \partial \varOmega _w, \\[-12pt]\nonumber \end{align}
\begin{align} p = \mathcal{P}_{\textit{out}}, \quad &\text{on} \quad {x} = 1, \quad y \in [0, \gamma ], \\[-12pt]\nonumber \end{align}
\begin{align} p = 0, \quad &\text{on} \quad y = - \lambda , \quad x \in \bigcup _{i=1}^N \left [ x_i - \frac {\delta \epsilon }{2}, x_i + \frac {\delta \epsilon }{2} \right ]\!, \end{align}
where
\begin{equation} x_i = \frac {(2 i - 1) \epsilon }{2}, \end{equation}
for
$i = 1, 2, {\cdots} , N$
and
\begin{equation} \mathcal{P}_{\textit{out}} = \frac {\epsilon ^2 L^2}{\rho _{\!f} \mathcal{Q}^2} \hat {\mathcal{P}}_{\textit{out}}, \end{equation}
which we assume is
${O}(1)$
.
3. Flow solution structure
We begin the solution procedure by finding the flux through each branched channel via the well-known problem of unidirectional fully developed flow.
3.1. Branched channel flow
The dimensionless Navier–Stokes equations given by (2.13) and (2.14) hold in each individual branched channel. Since we assume that
$\delta \epsilon \ll \lambda$
, each branched channel is long and thin, driven by a constant pressure drop with no slip on the walls due to viscous effects. The inlet, outlet and no-slip boundary conditions on an isolated branched channel are given by
\begin{align} p = p(x_i, 0), \quad &\text{on} \quad y = 0, \\[-12pt]\nonumber \end{align}
\begin{align} p = 0, \quad &\text{on} \quad y = - \lambda , \\[-12pt]\nonumber \end{align}
\begin{align} \boldsymbol{u} = \boldsymbol{0}, \quad &\text{on} \quad x = x_i \pm \frac {\delta \epsilon }{2}, \end{align}
respectively, where we assume that
$p(x_i, 0)$
is the constant pressure at the top and centre of each branched channel via Taylor expansion. We note that there will be a small region (see Appendix A) near the top of the T-junction where the flow develops into unidirectional flow, but we neglect this here. The leading-order solution for channel flow with no slip on the walls is
\begin{align} p &= p(x_i, 0) \left ( 1 + \frac {y}{\lambda } \right )\!, \\[-12pt]\nonumber \end{align}
\begin{align} u &= 0, \\[-12pt]\nonumber \end{align}
\begin{align} v &= \frac {\textit{Re}}{2 \epsilon ^2} \frac {\mathrm{d} p}{\mathrm{d} y} \left ( (x - x_i)^2 - \frac {\left ( \delta \epsilon \right )^2}{4} \right ) = \frac {\textit{Re}}{2 \epsilon ^2 \lambda } p(x_i,0) \left ( (x - x_i)^2 - \frac {\left ( \delta \epsilon \right )^2}{4} \right )\!. \end{align}
Therefore, the dimensionless flux through each
$i$
-th branched channel is given by
\begin{align} Q_i^{\textit{branch}} &= - \int _{x_i - \frac {\delta \epsilon }{2}}^{x_i + \frac {\delta \epsilon }{2}} v (x,0) \, \mathrm{d}x \\[-12pt]\nonumber \end{align}
\begin{align} &= - \left [ \frac {\textit{Re}}{2 \epsilon ^2 \lambda }p(x_i,0) \int ^{x_i + \frac {\delta \epsilon }{2}}_{x_i -\frac {\delta \epsilon }{2}}\left ( (x - x_i)^2 - \frac {(\delta \epsilon )^2}{4} \right ) \mathrm{d}x \right ] \\[-12pt]\nonumber \end{align}
\begin{align} &= \frac {\epsilon \delta ^3 \textit{Re}}{12 \lambda } p(x_i, 0). \end{align}
We suppose that the total fluid flux through all the branched channels over
$x \in [0,1]$
is
${O}(1)$
, so the fluid flux through each branched channel,
$Q_i^{\textit{branch}}$
, is
${O}(\epsilon )$
. Therefore,
\begin{equation} Q_i^{\textit{branch}} = \epsilon \kappa p(x_i,0), \end{equation}
where we define
\begin{equation} \kappa = \frac {\delta ^3 \textit{Re}}{12 \lambda } = {O}(1). \end{equation}
We now take the limit as
$\delta \rightarrow 0$
, with
$\kappa = {O}(1)$
, so that the width of each branched channel vanishes. This allows us to consider the remaining problem on a simplified domain
$(x,y) \in [0,1] \times [0, \gamma ]$
, where the branched channels are approximated by two-dimensional point sinks in the bulk flow located along the bottom wall. Since each sink only draws in fluid from
$y\gt 0$
, each has strength given by twice
$Q_i^{\textit{branch}}$
(Batchelor Reference Batchelor2000). Thus we replace (2.13) with
\begin{equation} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} =- 2\epsilon \kappa \sum _{i=1}^N p(x_i,0) \Delta \left (x-x_i, \, y \right )\!, \end{equation}
where
$\varDelta (\boldsymbol{\cdot },\boldsymbol{\cdot })$
is the two-dimensional delta function, and we henceforth focus on the geometry shown in figure 3.
Reduced dimensionless geometry, with point sinks replacing each branched channel. Each point sink has coordinates
$(x_i, 0)$
, where
$x_i = (i-1/2) \epsilon$
for
$i = 1, 2, {\cdots} , N$
, and has strength
$2 Q_i^{\textit{branch}}$
, where
$Q_i^{\textit{branch}}$
is the flux through a single channel, as in (3.10), (Batchelor Reference Batchelor2000). The outer problem views the point sinks as an effective boundary condition, capturing the overall average behaviour. Both boundary layers are indicated in the regime
$\epsilon \gg 1/\sqrt {\textit{Re}}$
.
3.2. Inviscid approximation
Since we assume that the Reynolds number is large, there exist viscous boundary layers, of thickness
$1/\sqrt {\textit{Re}}$
, on the top and bottom walls of the main channel. In addition, there is a further layer, of thickness
$\epsilon$
, over which the effects due to fluid loss down the branched channels will be smoothed out.
Outside of the viscous boundary layers and away from the bottom wall, we simplify the governing equations (2.13) and (2.14) by supposing that the flow is inviscid and irrotational. This leaves us with the system of equations
\begin{align} \frac {\partial {u}}{\partial {x}} + \frac {\partial {v}}{\partial {y}} &= 0, \\[-12pt]\nonumber \end{align}
\begin{align} \frac {\partial {v}}{\partial {x}} - \frac {\partial {u}}{\partial {y}} &= 0, \\[-12pt]\nonumber \end{align}
\begin{align} \frac {1}{\epsilon ^2} p + \frac {1}{2} |{\boldsymbol{u}}|^2 &= \text{constant}, \\[10pt] \nonumber \end{align}
outside of the viscous boundary layers, where we have not included the right-hand side of (3.12) since we are far away from the point sinks. We note that, at leading order in
$\epsilon$
, the pressure,
$p$
, given by (3.15), is constant everywhere.
When considering the combination of the two boundary layers on the bottom wall, there are three possible regimes: (i) the viscous boundary layer is larger such that
$\epsilon \ll 1/\sqrt {\textit{Re}}$
; (ii) the branched-channel-inducing boundary layer is larger, such that
$\epsilon \gg 1/\sqrt {\textit{Re}}$
; (iii) the layers are of the same order, i.e.
$\epsilon = O(1/\sqrt {\textit{Re}})$
. We focus our attention on case (ii), as indicated in figure 3, so that the flow inside the
$\epsilon$
-boundary layer is also governed by (3.13)–(3.15), and we neglect the viscous boundary layer for simplicity, since the focus of our interest is in the effect of the branched channels. The inviscid approximation of the flow everywhere reduces the complexity of the problem, and so we require fewer boundary conditions (see §§ 4 and 5).
Across the
$\epsilon$
-boundary layer, the flow generated by the point sinks will be smoothed out and our aim is to determine an effective boundary condition to be imposed on the outer flow. We adopt the following solution procedure. We first partially solve for the outer flow. Scaling into the
$\epsilon$
-boundary layer, we solve for the flow close to the point sinks. We then match between the two regions to find the effective boundary condition and a composite flow solution.
4. Outer flow problem
In the outer flow,
${\boldsymbol{u}}^{(o)} = ( {u}^{(o)},{v}^{(o)} )$
and
${p}^{(o)}$
denote the outer velocity and pressure, respectively, which satisfy (3.13)–(3.15). The boundary conditions for this outer problem (see figure 4) become
\begin{align} u^{(o)} = \frac {1}{\gamma }, \quad &\text{on} \quad x = 0, \\[-12pt]\nonumber \end{align}
\begin{align} v^{(o)} = 0, \quad &\text{on} \quad {y} = \gamma , \\[-12pt]\nonumber \end{align}
\begin{align} p^{(o)} = \mathcal{P}_{\textit{out}}, \quad &\text{on} \quad {x} = 1, \\[-12pt]\nonumber \end{align}
\begin{align} v^{(o)}_0 = - v^*(x), \quad &\text{on} \quad y = 0, \end{align}
Outer flow domain with boundary conditions, including the effective boundary condition,
$v^{(o)} (x, 0) = - v^* (x)$
.
where
$v^*(x)$
is an unknown function to be determined via matching to the inner problem.
4.1. Outer solution
To solve for the leading-order outer solution, we consider an asymptotic expansion in
$\epsilon$
of the form
\begin{equation} f^{(o)} (x,y) = f^{(o)}_0 (x,y) + \epsilon f^{(o)}_1 (x,y) + {\cdots} , \end{equation}
for the dependent variables
$u^{(o)}$
,
$v^{(o)}$
and
$p^{(o)}$
. Considering (3.15) with boundary condition (4.3), we find that the leading-order pressure is given by
\begin{equation} p^{(o)}_0 (x,y) = \mathcal{P}_{\textit{out}}, \end{equation}
everywhere. Since the outer flow is irrotational, we introduce a velocity potential,
$\phi ^{(o)}$
, such that
\begin{align} u^{(o)}_0 = \frac {\partial \phi ^{(o)}_0}{\partial x}, \quad v^{(o)}_0 = \frac {\partial \phi ^{(o)}_0}{\partial y}, \quad \end{align}
which will be useful for matching with the inner problem. To solve for the leading-order outer velocities,
$u^{(o)}_0$
and
$v^{(o)}_0$
, we first need to determine
$v^*(x)$
. This involves solving the inner problem and so we address this next, before returning to complete the solution in the outer flow in § 6.
5. Inner flow problem
Having partially solved for the outer problem, we now scale into the inner region via the scaling
\begin{equation} y = \epsilon Y. \end{equation}
In this region, the discrete nature of the point sinks becomes apparent, as illustrated in figure 3. We denote the inner velocity by
$\boldsymbol{u}^{(i)} = ( u^{(i)},v^{(i)})$
and pressure
$p^{(i)}$
, and the scaled versions of (3.12), (3.14), and (3.15) in the inner region are
\begin{align} \frac {\partial {u}^{(i)}}{\partial x} + \frac {1}{\epsilon } \frac {\partial v^{(i)}}{\partial Y} &= - 2\epsilon \kappa \sum _{i=1}^N p(x_i,0) \Delta \left (x-x_i, \, \epsilon Y \right ) \\[-12pt]\nonumber \end{align}
\begin{align} &= - 2 \kappa \sum _{i=1}^N p(x_i,0) \Delta \left (x-x_i, \, Y \right )\!, \\[-12pt]\nonumber \end{align}
\begin{align} \frac {\partial v^{(i)}}{\partial {x}} - \frac {1}{\epsilon } \frac {\partial u^{(i)}}{\partial Y} &= 0, \\[-12pt]\nonumber \end{align}
\begin{align} \frac {1}{\epsilon ^2} p^{(i)} + \frac {1}{2} \left ( \big (u^{(i)} \big )^2 + \big ( v^{(i)} \big )^2 \right) &= \text{constant} , \end{align}
where we have used the property that
$\Delta (\alpha x, \beta y) = \Delta (x, y) / \alpha \beta$
. The inner region boundary conditions are
\begin{align} u^{(i)} = \frac {1}{\gamma }, \quad &\text{on} \quad {x} = 0, \\[-12pt]\nonumber \end{align}
\begin{align} v^{(i)} = 0, \quad &\text{on} \quad Y = 0, \\[-12pt]\nonumber \end{align}
\begin{align} p^{(i)} = {\mathcal{P}}_{\textit{out}}, \quad &\text{on} \quad x = 1, \\[10pt] \nonumber \end{align}
with matching conditions to the outer region,
\begin{align} u^{(i)} (x, Y) \rightarrow u^{(o)}(x, 0), \quad &\text{as} \quad Y \rightarrow + \infty , \\[-12pt]\nonumber \end{align}
\begin{align} v^{(i)} (x, Y) \rightarrow v^{(o)} (x, 0), \quad &\text{as} \quad Y \rightarrow + \infty . \end{align}
As in the outer problem, since the problem is irrotational, we choose to express the inner problem in terms of a velocity potential,
$\phi ^{(i)}$
, defined by the ansatz
\begin{align} u^{(i)} = u^{(o)}_0 (x, \epsilon Y) + \epsilon \frac {\partial \phi ^{(i)}}{\partial x}, \quad v^{(i)} = v^{(o)}_0 (x, \epsilon Y) + \frac {\partial \phi ^{(i)}}{\partial Y}. \quad \end{align}
We take this form so that we automatically match to the leading-order outer velocities,
$u^{(o)}_0$
and
$v^{(o)}_0$
, where any variations are introduced via
$\phi ^{(i)}$
. Substituting (5.11) into (5.3) gives
\begin{equation} \epsilon \frac {\partial ^2 \phi ^{(i)}}{\partial x^2} + \frac {1}{\epsilon } \frac {\partial ^2 \phi ^{(i)}}{\partial Y^2} = - 2 \kappa \sum _{i=1}^N p^{(i)}(x_i, 0) \Delta \left ( x-x_i, \, Y \right )\!, \end{equation}
since
\begin{equation} \frac {\partial {u}^{(o)}_0}{\partial x} (x, \epsilon Y) + \frac {1}{\epsilon } \frac {\partial {v}^{(o)}_0}{\partial Y} (x, \epsilon Y) = 0, \end{equation}
from (3.13) in the outer problem, and (5.5) becomes
\begin{equation} \frac {1}{\epsilon ^2} p^{(i)} + \frac {1}{2} \left [ \left ( u^{(o)}_0 (x, \epsilon Y) + \epsilon \frac {\partial \phi ^{(i)}}{\partial x} \right )^2 + \left ( v^{(o)}_0 (x, \epsilon Y) + \frac {\partial \phi ^{(i)}}{\partial Y} \right )^2 \right ] = \text{constant}. \end{equation}
The boundary conditions (5.6)–(5.8) become
\begin{align} \frac {\partial \phi ^{(i)}}{\partial x} = 0, \quad &\text{on} \quad {x} = 0, \\[-12pt]\nonumber \end{align}
\begin{align} \frac {\partial \phi ^{(i)}}{\partial Y} = 0, \quad &\text{on} \quad Y = 0, \\[-12pt]\nonumber \end{align}
\begin{align} p^{(i)} = {\mathcal{P}}_{\textit{out}}, \quad &\text{on} \quad x = 1. \end{align}
Before proceeding, we comment on the matching conditions (5.9) and (5.10) between the inner and outer solutions. Given our ansatz for
$u^{(i)}$
and
$v^{(i)}$
in (5.11), the leading-order matching is automatically satisfied. We also follow Van Dyke (Reference Van Dyke1975) to establish an additional matching condition between the inner and outer solutions, using
\begin{equation} { 1 \, \text{term outer} \, ( 2 \, \text{term inner} ) = 2 \, \text{term inner} \, (1 \, \text{term outer} ).} \end{equation}
Writing
$\phi ^{(i)} \sim \phi ^{(i)}_1 + \epsilon \phi ^{(i)}_2 + {\cdots}$
, this matching condition becomes
\begin{align} \lim _{Y \rightarrow + \infty } \frac {\partial \phi ^{(i)}_1}{\partial Y} (x,Y) = \frac {\partial \phi ^{(o)}_0}{\partial y} (x, 0) = -v^*(x). \end{align}
5.1. Method of multiple scales
The domain within the inner region is made up of a repeating structure of length
$\epsilon$
, surrounding each point sink. We describe each part of the repeating structure as a cell. We assume that the flow is slowly varying over each cell. We introduce a microscale variable,
$X$
, by letting
$x = \epsilon X$
, to capture both the behaviour over the domain length and over the length of each cell. We treat both the long scale,
$x$
, and short scale,
$X$
, as independent variables and assume that all dependent variables depend on
$x, X$
and
$Y$
in the inner region. Since we treat each of these variables as independent, spatial derivatives in
$x$
transform as
\begin{equation} \frac {\partial }{\partial x} \mapsto \frac {\partial }{\partial x} + \frac {1}{\epsilon } \frac {\partial }{\partial X}. \end{equation}
The sink in a particular cell is located at
$X=0$
. The addition of the microscale variable also gives us an additional degree of freedom. We remove this by imposing periodicity on the boundary walls such that
\begin{equation} \frac {\partial \phi ^{(i)}}{\partial X} \bigg \rvert _{X=-\frac {1}{2}} = \frac {\partial \phi ^{(i)}}{\partial X} \bigg \rvert _{X = \frac {1}{2}} = 0. \end{equation}
This retains the slow variation of
$u^{(i)}$
over each cell, where (5.11a
) becomes
\begin{align} u^{(i)} = u^{(o)}_0 (x, \epsilon Y) + \frac {\partial \phi ^{(i)}}{\partial X}+ \epsilon \frac {\partial \phi ^{(i)}}{\partial x}. \end{align}
We solve for
$\phi ^{(i)}$
in each periodic cell as a point-sink cell problem in a periodic semi-infinite half-strip, and match to the outer solution. The inner bulk equation (5.12) in the cell problem becomes
\begin{equation} \frac {\partial ^2 \phi ^{(i)}}{\partial X^2} + \frac {\partial ^2 \phi ^{(i)}}{\partial Y^2} + 2 \epsilon \frac {\partial ^2 \phi ^{(i)}}{\partial x \partial X} + \epsilon ^2 \frac {\partial ^2 \phi ^{(i)}}{\partial x^2} = - 2 \kappa p^{(i)}(x_i, 0, 0) \Delta (X, Y), \end{equation}
where we have used the fact that
$\epsilon \Delta (x-x_i,Y)=\varDelta (X,Y)$
.
5.2. Inner solution
We consider an asymptotic expansion of the inner pressure to be of the form
\begin{equation} p^{(i)} (x, X, Y) = p^{(i)}_0 (x, X, Y) + \epsilon p^{(i)}_1 (x, X, Y) + {\cdots} . \end{equation}
We first solve for the leading-order pressure in the inner region. The leading-order version of (5.14) implies that
$p^{(i)}_0$
is a constant and so, matching with the solution in the outer flow given by (4.6), we find that
\begin{equation} p^{(i)}_0 (x, X, Y) = \mathcal{P}_{\textit{out}}. \end{equation}
Thus, the leading-order pressure is constant everywhere in the flow. The problem for
$\phi ^{(i)}_1$
is given by
\begin{equation} \frac {\partial ^2 \phi ^{(i)}_1}{\partial X^2} + \frac {\partial ^2 \phi ^{(i)}_1}{\partial Y^2} = - 2 \kappa \mathcal{P}_{\textit{out}} \varDelta (X, Y), \end{equation}
with boundary and matching conditions
\begin{align} \frac {\partial \phi ^{(i)}_1}{\partial X} = 0, \quad &\text{on} \quad X = \pm \frac {1}{2}, \; Y \in (0,\infty ), \\[-12pt]\nonumber\end{align}
\begin{align} \frac {\partial \phi ^{(i)}_1}{\partial Y} = 0, \quad &\text{on} \quad Y=0, \; X \in \left [-\frac {1}{2}, \frac {1}{2} \right ]\!, \\[-12pt]\nonumber \end{align}
\begin{align} \frac {\partial \phi ^{(i)}_1}{\partial Y} \rightarrow -v^*(x), \quad &\text{as} \quad Y \rightarrow +\infty, \; X \in \left [-\frac {1}{2}, \frac {1}{2} \right ]\!. \end{align}
We use complex variable theory to find the solution for
$\phi ^{(i)}_1$
in the semi-infinite periodic strip with a point sink at the origin of strength
$2 \kappa \mathcal{P}_{\textit{out}}$
. Using the holomorphic conformal mapping
\begin{equation} \zeta = \xi + \mathrm{i} \eta = \sin {(\pi Z)}, \end{equation}
where
$Z = X + \mathrm{i} Y$
, the semi-infinite periodic strip domain is mapped to a semi-infinite half-plane,
$\eta = \textrm{Im} (\zeta ) \geqslant 0$
, as indicated in figure 5, where we denote
$\textrm{Re}$
and
$\textrm{Im}$
as the real and imaginary parts, respectively. Laplace’s equation still holds in the transformed domain, however, the periodic conditions become no-flux conditions on the half-plane boundary.
Conformal map of the semi-infinite half-strip inner region to the positive imaginary half-plane via the conformal map
$\zeta = \sin {(\pi Z)}$
.
The problem in the conformally mapped space is thus
\begin{equation} {\nabla} ^2_\zeta \phi ^{(i)}_1 = - 2 \kappa \mathcal{P}_{\textit{out}} \Delta {(\zeta )}, \end{equation}
with boundary and matching conditions
\begin{align} \frac {\partial \phi ^{(i)}_1}{\partial \eta } = 0, \quad &\text{on} \quad \textrm{Im} (\zeta ) = 0, \\[-12pt]\nonumber \end{align}
\begin{align} \frac {\partial \phi ^{(i)}_{1}}{\partial \eta } \rightarrow - v^*(x), \quad &\text{as} \quad \textrm{Im} (\zeta ) \rightarrow \infty . \end{align}
The solution to (5.31)–(5.33) is the well-known solution to Laplace’s equation with a point sink/source
\begin{equation} \phi _1^{(i)}(x,\zeta ) = C(x)- \frac {\kappa \mathcal{P}_{\textit{out}}}{\pi } \log (\zeta ), \end{equation}
where
$C(x)$
is a function of integration. Transforming back into
$(x,X,Y)$
coordinates, we have
\begin{equation} \phi _1^{(i)}(x,X,Y) = C(x)- \frac {\kappa \mathcal{P}_{\textit{out}}}{\pi } \textrm{Re} [ \log ( \sin { (\pi (X + \mathrm{i} Y) } ) ]. \end{equation}
The function of integration,
$C(x)$
, may be taken to be zero without loss of generality. Therefore, since
\begin{align} \frac {\partial \phi ^{(i)}_1}{\partial X} &= - \kappa \mathcal{P}_{\textit{out}} \textrm{Re} \left [ \cot { \left ( \pi \left ( X + \mathrm{i} Y \right ) \right )} \right ] \\[-12pt]\nonumber \end{align}
\begin{align} &= - \frac { \kappa \mathcal{P}_{\textit{out}} \sin {\left (\pi X \right )} \cos {\left ( \pi X \right )} }{ \cos ^2{\left (\pi X \right )} \sinh ^2{\left (\pi Y \right )} + \sin ^2{\left (\pi X \right )} \cosh ^2{\left (\pi Y \right )} }, \\[6pt] \nonumber\end{align}
and
\begin{align} \frac {\partial \phi ^{(i)}_1}{\partial Y} &= \kappa \mathcal{P}_{\textit{out}} \textrm{Im} \left [ \cot { \left ( \pi \left ( X + \mathrm{i} Y \right ) \right )} \right ] \\[-12pt]\nonumber \end{align}
\begin{align} &= - \frac { \kappa \mathcal{P}_{\textit{out}} \cosh {\left (\pi Y \right )} \sinh {\left ( \pi Y \right )} }{ \cos ^2{\left (\pi X \right )} \sinh ^2{\left (\pi Y \right )} + \sin ^2{\left (\pi X \right )} \cosh ^2{\left (\pi Y \right )} }, \end{align}
the leading-order
$x$
-component of the inner velocity is given by
\begin{align} u^{(i)}_0 (x, X, Y) = u^{(o)}_0 (x, 0) - \frac { \kappa \mathcal{P}_{\textit{out}} \sin {\left (\pi X \right )} \cos {\left ( \pi X \right )} }{ \cos ^2{\left (\pi X \right )} \sinh ^2{\left (\pi Y \right )} + \sin ^2{\left (\pi X \right )} \cosh ^2{\left (\pi Y \right )} }, \end{align}
and the leading-order
$y$
-component of the inner velocity is given by
\begin{align} v^{(i)}_0 (x, X, Y) = v^{(o)}_0 (x, 0) - \frac { \kappa \mathcal{P}_{\textit{out}} \cosh {\left (\pi Y \right )} \sinh {\left ( \pi Y \right )} }{ \cos ^2{\left (\pi X \right )} \sinh ^2{\left (\pi Y \right )} + \sin ^2{\left (\pi X \right )} \cosh ^2{\left (\pi Y \right )} }. \end{align}
Having found the solution in the inner region, we calculate that, as
$Y \rightarrow \infty$
,
\begin{align} \frac {\partial \phi ^{(i)}_1}{\partial Y} &\rightarrow - \kappa \mathcal{P}_{\textit{out}} \frac { e^{\pi Y} e^{\pi Y} }{ \cos ^2{\left (\pi X \right )} e^{2 \pi Y} + \sin ^2{\left (\pi X \right )} e^{2 \pi Y} } = - \kappa \mathcal{P}_{\textit{out}}. \end{align}
Therefore, via the matching condition (5.29), we find that
\begin{equation} v^* (x) = \kappa \mathcal{P}_{\textit{out}}, \end{equation}
and thus
\begin{equation} v^{(o)}_0 (x, 0) =-\kappa \mathcal{P}_{\textit{out}}. \end{equation}
Thus, the boundary layer smooths out the variation caused by the discrete sinks, and the outer flow simply experiences a spatially uniform flow of liquid out through the ‘effective’ bottom of the device.
6. Outer solution
Now that we have found
$v^*(x) = \kappa \mathcal{P}_{\textit{out}}$
, we return to the outer problem to find the leading-order outer velocities
$u^{(o)}_0$
and
$v^{(o)}_0$
. Solving (3.13) and (3.14), with boundary conditions (4.1), (4.2), (4.4) and (5.44), we find the simple solution
\begin{align} u^{(o)}_0 (x,y) = \frac {1 - \kappa \mathcal{P}_{\textit{out}} x}{\gamma }, \qquad v^{(o)}_0 (x,y) = -\kappa \mathcal{P}_{\textit{out}}\left (1 - \frac {y}{\gamma }\right )\!. \end{align}
We see that
$u_0^{(o)}$
is independent of
$y$
and so
$u^{(o)}_0 (x, \epsilon Y)=u^{(o)}_0 (x, 0)=u^{(o)}_0 (x, y)$
.
7. Composite solution and outflow flux
Using (5.40) and (6.1a
), we write the leading-order composite solution for
$u$
as
\begin{align} \hspace {-1mm} u_c (x, y) = \frac {1 - \kappa \mathcal{P}_{\textit{out}} x}{\gamma } - \dfrac { \kappa \mathcal{P}_{\textit{out}} \sin {\left [\dfrac {\pi }{\epsilon } \left ( x - \dfrac {\epsilon }{2} \right ) \right ]} \cos {\left [ \dfrac {\pi }{\epsilon } \left ( x - \dfrac {\epsilon }{2} \right ) \right ]} }{ \cos ^2{\left [\dfrac {\pi }{\epsilon } \left ( x - \dfrac {\epsilon }{2} \right ) \right ]} \sinh ^2{\left [\dfrac {\pi y}{\epsilon } \right ]} + \sin ^2{\left [\dfrac {\pi }{\epsilon } \left ( x - \dfrac {\epsilon }{2} \right ) \right ]} \cosh ^2{\left [\dfrac {\pi y}{\epsilon } \right ]} }, \end{align}
where we have taken care to correctly locate the sinks by setting
$X=(x-\epsilon /2)/\epsilon$
. Similarly, using (5.41) and (6.1b
), we find the leading-order composite solution for
$v$
is
\begin{align} \hspace {-1.5mm} v_c (x,y) = \kappa \mathcal{P}_{\textit{out}} \left ( \dfrac {y}{\gamma } -\dfrac { \cosh {\left [\dfrac {\pi y}{\epsilon } \right ]} \sinh {\left [ \dfrac {\pi y}{\epsilon } \right ]} }{ \cos ^2{\left [\dfrac {\pi }{\epsilon } \left ( x - \dfrac {\epsilon }{2} \right ) \right ]} \sinh ^2{\left [\dfrac {\pi y}{\epsilon } \right ]} + \sin ^2{\left [\dfrac {\pi }{\epsilon } \left ( x - \dfrac {\epsilon }{2} \right ) \right ]} \cosh ^2{\left [\dfrac {\pi y}{\epsilon } \right ]} } \right )\!. \end{align}
These explicit formulae for the velocities bypass the need to run computationally challenging numerical simulations.
We calculate the total flux,
$Q_T$
, through the lower boundary, using
\begin{equation} Q_T = - \int ^1_0 v^*(x) \, \mathrm{d}x = \frac {\delta ^3 {Re}}{12 \lambda } \mathcal{P}_{\textit{out}}, \end{equation}
while the flux through the main outlet is
$1 - Q_T$
. Hence, the dimensional flux,
$\hat {Q}_T$
, is
\begin{equation} \hat {Q}_T= \frac { h_2^3 N}{12 L_1 \mu } \hat {\mathcal{P}}_{\textit{out}}. \end{equation}
Thus we see that the flux out through the bottom boundary scales in the obvious way –
$N$
times the flux out of a single channel in the limit in which the pressure is constant.
8. Flow results
To validate our asymptotic prediction of the fluid flow and effective boundary condition, we compare with numerical solutions of the flow in the full branched domain.
8.1. Numerical results
We carry out numerical simulations in COMSOL in which we capture the nature of the high-Reynolds-number flow in the main channel and the Poiseuille flow in the branched channels. Since, in the regime of interest, the viscous boundary layers are much thinner than the
$\epsilon$
-layers and not a focus of this study, we impose free slip on the horizontal surfaces located at
$y=0$
and
$y=\gamma$
to ensure that the flow in the main channel emulates inviscid flow. We illustrate the solution structure and the relevant result comparisons in Appendix B. Thus, we solve (2.13)–(2.20), with (2.16) removed and (2.17) partially replaced by the free-slip condition on the main channel walls. We refine the mesh such that the total flux through the thin branched channels converges to a constant value for any smaller mesh size. We use a once refined free triangular, extremely fine mesh in the outer domain and two mapped distribution meshes along
$y=0$
and the walls of the branched channels. This achieves a convergent solution for the flow through the small branched channels, stable for any smaller mesh size.
8.2. Asymptotic and numerical comparison
We show the magnitude of the flow velocity resulting from our simulations and the streamlines of the flow, for an example set of parameters in which
$\mathcal{P}_{ {out}}=0.4$
,
$\textit{Re}=1000$
,
$\epsilon =0.04$
,
$\delta =0.1$
,
$\lambda =0.1$
and
$\gamma =0.5$
, in figure 6. We note that these values are on the edge of the asymptotic regime. However, this is the largest Reynolds number that we can achieve in COMSOL given the geometric complexities we are considering, but we also explore
$\epsilon = 0.1$
(which has fewer channels than in the Beko PLC device) which sits more squarely in this asymptotic regime. We see that the liquid flows across from left to right and down, as expected, with a dividing streamline separating the flow that exits through the branched channels from the flow that leaves through the main channel outlet. We see in figure 7(a) that the flow velocity is largest in the branched channels. We calculate the total flux through each of the branches,
$\sum _i Q_i^{\textit{branch}}$
to be
$0.315$
which corresponds to two times the height of the dividing streamline at
$x=0$
– this is because the inlet flux is
$q=1$
and the height of the main channel is
$\gamma = 0.5$
. For the same parameter values, the asymptotic prediction (7.3) gives
$Q = 1/3$
, which is an
${O}(\epsilon )$
difference.
Numerical solution for the magnitude of the flow velocity,
$|\boldsymbol{u}|$
, solved via the Navier–Stokes equations (2.13)–(2.14). We apply a slip condition on the main channel walls and no slip on the branched channel walls. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
. The black lines indicate streamlines and the red line indicates the dividing streamline.
Numerical solutions, zoomed in to individual branched channels, for (a) the magnitude of the flow velocity,
$|\boldsymbol{u}|$
, and (b) pressure,
$p$
, from figures 6 and 8, respectively, with the full colour range for
$p$
. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
. The black arrowed line indicates a particular streamline.
We show the corresponding pressure in figure 8. We see that the pressure is almost constant in the main channel, and then decreases to zero along the branched channels (see figure 7
b). Maybe surprisingly at a first glance, we see that the pressure increases from left to right. The simulations compare well with the leading-order asymptotic prediction,
$p=\mathcal{P}_{\textit{out}}$
, with deviations of
${O}(\epsilon )$
over the domain.
Numerical solution for the pressure,
$p$
, corresponding to figure 6, solved via the Navier–Stokes equations (2.13)–(2.14). Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
. The colour range for the pressure is restricted to
$[0.397,0.401]$
, and white otherwise, to highlight the asymptotic observation of constant pressure to leading order in
$\epsilon$
over the main channel. The branched channels are shown to illustrate their location.
We compare our asymptotic composite solution for
$v(x^*, y)$
, given by (7.2), with the numerical solution by plotting
$v$
versus
$y$
at the centre of a cell,
$x^* = x_i$
(figure 9
a) and at the right-hand edge of the cell,
$x_i + \epsilon /2$
, (figure 9
b), for
$\epsilon =0.04$
. We see that the asymptotic and numerical results agree well and that the agreement improves as
$\epsilon$
decreases whilst keeping
$\textit{Re}$
fixed, as expected (see figures 9 and 10). We clearly see the presence of the boundary layer at
$y=0$
, with the solution getting sharper as
$\epsilon$
decreases as seen in figure 11, highlighting that the leading-order outer flow is a good approximation for all
$y$
, as
$\epsilon \rightarrow 0$
.
Comparison of the leading-order asymptotic composite solution (7.2) in red dashed lines, with the numerical solution in solid black lines, for the
$y$
-component of velocity,
$v$
. The solutions are taken at (a)
$x = x_i$
, the centre of the branched channel/point sink and (b)
$x = x_i + \epsilon /2$
, the right-hand edge of the periodic unit cell. In this example, we consider the centremost cell, taking
$x_i = 0.5$
, and vary the cell width. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
Comparison of the leading-order asymptotic composite solution (7.2) in red dashed lines, with the numerical solution in solid black lines, for the
$y$
-component of velocity,
$v$
. The solutions are taken at (a)
$x = x_i$
, the centre of the branched channel/point sink and (b)
$x = x_i + \epsilon /2$
, the right-hand edge of the periodic unit cell. In this example, we consider the centremost cell, taking
$x_i = 0.55$
, and vary the cell width. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.1$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
Leading-order asymptotic composite solution (7.2) for
$v$
at
$x = x_i + \epsilon /2$
for
$\epsilon = 0.04, {} 0.1, \, 0.125, \, 1/6$
. The outer solution for the velocity (6.1) is indicated in purple. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
. For varying
$\epsilon$
, we remain in the correct limit,
$\epsilon \gg 1/\sqrt {\textit{Re}}$
.
We similarly compare the asymptotic solution for
$u(x^*, y)$
given by (7.1) with the numerical solution along the centre of a cell, at
$x^* = x_i$
(figure 12
a), and at the right-hand edge of the cell,
$x_i + \epsilon /2$
, (figure 12
b), for
$\epsilon =0.04$
. Once again, we see that the asymptotic and numerical results agree well, apart from in a region close to
$y=0$
. We notice a significant disagreement between the numerics and the asymptotics along the centre of a cell. In figure 12(a), we plot additional
$u(x^*, y)$
for
$x^*$
between the edge of the branched channel and the edge of the cell. For these additional values, the asymptotics agree reasonably well with the numerics, which suggests the disagreements are due to the finite size of the channel entrances. This error would reduce by decreasing the width of the channels in the numerical simulations. We note in passing, in figure 13(a), that
$u$
drops below zero around the right-hand corner of the branched channel entrance indicating a region of flow reversal.
Comparison of the leading-order asymptotic composite solution (7.1) in dashed lines, with the numerical solution in solid lines, for the
$x$
-component of velocity,
$u$
. Here, we only plot close to
$y=0$
as the solution is constant for larger
$y$
. The numerical solution in solid black lines and asymptotic solution in red dashed lines are taken at (a)
$x = x_i$
, the centre of the branched channel/point sink and (b)
$x = x_i + \epsilon /2$
, the right-hand edge of the periodic unit cell. In this example, we consider the centremost cell, taking
$x_i = 0.5$
. The additional solutions in (a) are taken at
$x=0.49$
(blue) and
$x=0.51$
(green). Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
Numerical solutions, zoomed in to individual branched channels, for (a) the
$x$
-component of the flow speed,
$u$
, and (b) the
$y$
-component of the flow speed,
$v$
. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
In figure 14(a), we compare the asymptotic solution for
$u(x,y^*)$
given by (7.1) with the numerical solution, for three values of
$y^*$
. We observe excellent agreement between the asymptotic and numerical solutions. We further see minimal
$y$
-dependence in both the numerical results and composite solution away from
$y=0$
. As we approach the wall, we see that, in the boundary layer, the asymptotic solution accurately captures the oscillations seen in the numerics. Our results indicate that the discrete behaviour of the branched channels is confined to the inner region and the outer flow only sees this as an effective boundary condition. We similarly compare the asymptotic and numerical solutions for
$v(x, y^*)$
given by (7.2) in figure 14, and again see excellent agreement for the three values of
$y^*$
.
Comparison of the leading-order asymptotic composite solution in red dashed lines, with the numerical solution in solid black lines, for (a) the
$x$
-component of velocity,
$u$
, given by (7.1) and (b) the
$y$
-component of velocity,
$v$
, given by (7.2). The solutions are taken at
$y = \epsilon /2$
,
$0.25$
and
$\gamma$
. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
In figure 15(a), we compare the flux through each branched channel found numerically with the leading-order flux found in (3.10), and observe an
${O}(\epsilon ^2)$
discrepancy. Integrating the effective boundary condition (7.3) over
$x$
, we find that the cumulative flux along the wall is given by
\begin{equation} Q^{\textit{cumulative}}(x) = \kappa \mathcal{P}_{\textit{out}} x. \end{equation}
In figure 15(b), we compare this with the cumulative flux through each branch found numerically, and again we see excellent agreement, albeit with an error that grows linearly since we are integrating values that have systematic and asymptotically small errors.
Comparison of the (a) flux through each branched channel, recorded at
$x = x_i$
, and (b) cumulative flux through each branched channel, recorded at
$x = \epsilon i$
, for
$i = 1, 2, {\cdots} , N$
. We plot the numerical solution of the flux (black dots), the asymptotic flux (red dashed lines) given by (3.10) and (8.1) and the flux (3.10) using the pressure found numerically at the top of each branched channel,
$p(x_i, 0)$
, (blue dots). Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
8.3. Solution effects due to angled branches
We now consider the effect of angling the branched channels. We denote the angle between the branched channel walls and the negative
$y$
-axis by
$\alpha$
, with
$\alpha \gt 0$
corresponding with channels that are angled in the overall flow direction. We consider the cases where either (a) the branched channel width stays constant as
$\alpha$
increases, or (b) the hole size stays constant as
$\alpha$
increases, as shown in figure 16.
Angled geometry with (a) constant branch width
$\delta \epsilon$
and hole width
$w$
, and (b) constant hole width
$\delta \epsilon$
. The angle,
$\alpha$
, is defined between the branched channels walls and the negative
$y$
-axis.
Before we consider each case, as in § 3.1, we draw attention to the behaviour around the entrance to the branched channels. We have assumed that the flow within the branched channels is unidirectional, however, at the entrance to the branched channel, there is a small region where the flow will not be fully developed. Once again, we neglect any complex behaviour in this region (see Appendix A). Under this assumption, and given that the resistance to the flow is along the channel, varying
$\alpha$
while keeping the branch width
$\delta \epsilon$
constant does not affect the boundary condition (4.4). Therefore the effective flux (7.3) stays the same. The only requirement is that the hole size,
$w$
, is sufficiently small when compared with the cell width, such that we may approximate the hole by a point sink. This means that we require
\begin{equation} w \ll \epsilon \quad \Rightarrow \quad \delta \ll \cos {\alpha }. \end{equation}
On the other hand, holding the hole size constant and varying
$\alpha$
changes the width of the branched channel to
$\delta \epsilon \cos {\alpha }$
, resulting in
\begin{equation} v^* = - \frac {\delta ^3 \textit{Re} \cos ^3{\alpha }}{12 \lambda } \mathcal{P}_{\textit{out}}, \end{equation}
and so the total flux through the branched channels is
\begin{equation} Q = \frac {\delta ^3 \textit{Re} \cos ^3{\alpha }}{12 \lambda } \mathcal{P}_{\textit{out}}, \end{equation}
i.e. the flux decreases as
$\alpha$
increases.
We compare these predictions with their numerical equivalents as
$\alpha$
is varied, in figure 17. We restrict our attention to angles
$\alpha \leqslant \pi /3$
, for which
$\delta = 0.1 \ll 0.5 \lt \cos {\alpha }$
. The asymptotic solutions match the numerical simulations well in both cases, with an
${O}(\epsilon )$
over-prediction, which decreases as
$\alpha$
approaches
$\pi /3$
. As we increase
$\alpha$
further, past
$\pi /3$
, we see an under-prediction in the constant width scenario, but continued agreement in the constant hole scenario. In both cases, we are ultimately restricted by
$\alpha \lt \pi /2$
. A similar result to figure 17(a) is found by Mao et al. (Reference Mao, Bischofberger and Hosoi2024), determining that the flux through the branched channel has minimal variation when changing the angle,
$\alpha$
, for a constant channel width.
Comparison of the total flux through the discrete branches,
$Q$
, found numerically (black solid) with the asymptotic prediction of the effective boundary condition (red dashed), for (a) constant branch channel width and (b) constant branch hole width. The asymptotic flux is given by (7.3) and (8.4), respectively. Numerical solutions are found similarly to figure 6, with the same parameter values, but with the addition of the respective angled branched channels. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
9. Particle model derivation
We now consider the behaviour of individual spherical particles moving through the device in the main channel,
$(\hat {x}, \hat {y}) \in [0,L] \times [0, h_1]$
. We suppose for simplicity that the only force acting on a particle is fluid drag, which has quadratic dependence on the velocity due to the high speed of the flow (Benchaita, Griffith & Rabinowicz Reference Benchaita, Griffith and Rabinowicz1983). We follow Vigolo et al. (Reference Vigolo, Griffiths, Radl and Stone2013) to establish a force balance for each particle, given by
${\hat {\boldsymbol{F}}} = m {\hat {\boldsymbol{a}}}$
, but significantly reduce the complexity by removing all forces other than drag. In doing this, we write
\begin{equation} \frac {4}{3} \rho _p \pi a^3 \frac {{d}^2 \hat {\boldsymbol{x}}_p}{\mathrm{d} \hat {t}^2} = \frac {C_D}{2} \rho _{\!f} \pi a^2 \, \bigg \rvert \, \hat {\boldsymbol{u}} - \frac {{\mathrm d}\hat {\boldsymbol{x}}_p}{{\mathrm d}\hat {t}} \bigg \rvert \left ( \hat {\boldsymbol{u}} - \frac {{\mathrm d}\hat {\boldsymbol{x}}_p}{{\mathrm d}\hat {t}} \right )\!, \end{equation}
where
$\hat {\boldsymbol{x}}_p = (\hat {x}_p (\hat {t}), \hat {y}_p (\hat {t}))$
is the particle position,
$a$
and
$\rho _p$
denote the particle radius and density, respectively,
$\hat {\boldsymbol{u}} = (\hat {u}, \hat {v})$
is the flow velocity and
$C_D$
is a drag coefficient (see, for example, Benchaita et al. Reference Benchaita, Griffith and Rabinowicz1983). We release a particle at the inlet,
$\hat {x}=0$
, with position
$\hat {y} = \hat {y}_0$
and assume that it has initial velocity given by the inlet flow velocity
\begin{align} \hat {\boldsymbol{x}}_p {(0)} &= (0, \hat {y}_0), \\[-12pt]\nonumber \end{align}
\begin{align} \frac {{\mathrm d}\hat {\boldsymbol{x}}_p}{{\mathrm d}\hat {t}} {(0)} &= \hat {\boldsymbol{u}} (0, \hat {y}_0). \end{align}
We assume particle collisions with the walls are perfectly elastic, and model the ricochet effect via a simple bounce condition, i.e. angle of incidence
$=$
angle of reflection, giving
\begin{equation} \left (\frac {{\mathrm d}\hat {x}_p}{{\mathrm d}\hat {t}}, \frac {{\mathrm d}\hat {y}_p}{{\mathrm d}\hat {t}} \right ) \bigg \rvert ^+ = \left (\frac {{\mathrm d}\hat {x}_p}{{\mathrm d}\hat {t}}, - \frac {{\mathrm d}\hat {y}_p}{{\mathrm d}\hat {t}} \right ) \bigg \rvert ^- \quad \text{on all walls,} \end{equation}
where
$\boldsymbol{\cdot }|^-$
means before the collision and
$\boldsymbol{\cdot }|^+$
means after the collision has occurred.
9.1. Non-dimensionalisation
We non-dimensionalise the model (9.1)–(9.4) using the scalings
\begin{equation} \hat {\boldsymbol{x}}_p = L \boldsymbol{x}_p, \qquad \hat {\boldsymbol{u}} = \frac {\mathcal{Q}}{L} \boldsymbol{u}, \qquad \hat {t} = \frac {L^2}{\mathcal{Q}} t, \end{equation}
to find that
\begin{equation} {\textit{St}} \frac {{d}^2 \boldsymbol{x}_p}{\mathrm{d} t^2} = \, \bigg \rvert \, \boldsymbol{u} - \frac {{\mathrm d}\boldsymbol{x}_p}{{\mathrm d}t} \bigg \rvert \left ( \boldsymbol{u} - \frac {{\mathrm d}\boldsymbol{x}_p}{{\mathrm d}t} \right )\!, \end{equation}
where
\begin{equation} \mbox{${St}$} = \frac {8}{3} \frac {a \rho _p}{C_D L \rho _{\!f}} \end{equation}
is the Stokes number. When
$\mbox{${St}$} = 0$
, particles are infinitely light tracer particles and follow the flow exactly. Otherwise,
$\mbox{${St}$} \gt 0$
represents a particle with some non-zero density.
The inlet/initial conditions (9.2) and (9.3) become
\begin{align} \boldsymbol{x}_p {(0)} &= (0, y_0), \\[-12pt]\nonumber \end{align}
\begin{align} \frac {{\mathrm d}\boldsymbol{x}_p}{{\mathrm d}t} {(0)} &= \boldsymbol{u}(0, y_0), \end{align}
where
$y_0 = \hat {y}_0/L$
, and the bounce condition (9.4) on the walls becomes
\begin{equation} \left (\frac {{\mathrm d}{x}_p}{{\mathrm d}{t}}, \frac {{\mathrm d}{y}_p}{{\mathrm d}{t}} \right ) \bigg \rvert ^+ = \left (\frac {{\mathrm d}{x}_p}{{\mathrm d}{t}}, - \frac {{\mathrm d}{y}_p}{{\mathrm d}{t}} \right ) \bigg \rvert ^- \!. \end{equation}
9.2. Suitable flow field
In the flow problem in §§ 2–7, we identified that there is an outer region away from the bottom wall, and an inner region of thickness
$\epsilon$
close to the bottom wall. We found that the flow profile is given by the composite solution
$\boldsymbol{u}_c = (u_c, v_c)$
shown in (7.1) and (7.2). The effect of the wall only becomes noticeable at a substantial distance into the boundary layer,
$y \approx \epsilon /2$
. A key question is whether the outer flow provides a sufficient enough description of the flow to accurately capture the particle trajectories, or whether the (computationally more expensive) composite solution is required. To answer this question, we solve (9.6)–(9.10) for several trajectories using (i) the outer flow, (ii) the composite flow, (iii) the full numerical flow and compare and contrast the results. We note that to ensure the flux out,
$Q$
, is the same in the asymptotic and numerical solutions, we must take
$\mathcal{P}_{\textit{out}} = 0.424$
in the numerical solution to compare with the asymptotic solution exactly, when
$\mathcal{P}_{\textit{out}} = 0.4$
, so that
$Q = 1/3$
in both. We explore two distinguished limits for the Stokes number,
$\mbox{${St}$} = {O}(1)$
and
$\mbox{${St}$}={O}(\epsilon )$
.
When simulating particle trajectories using the asymptotic results, we exploit the fact that the
$y$
-component of the velocity is an odd function in
$y$
, i.e.
$v_c (x, -y) = - v_c (x,y)$
. We therefore solve the governing equation (9.6) with inlet conditions (9.8) and (9.9) in the domain
$[0,1] \times [-\gamma, \gamma ]$
allowing the particle to pass through
$y=0$
unless the trajectory enters a branched channel. When a particle passes through the wall at
$y=0$
, the angle of incidence is equal to the angle of reflection. We then take the modulus of the
$y$
-component of the trajectory, which ensures the bounce condition is satisfied for
$y \geqslant 0$
. As noted, this would be unsuitable if a particle trajectory enters a branched channel – in our comparisons below, the specific values of
$y_0$
are chosen to ensure that each particle lands away from a branched channel on every bounce.
9.2.1.
$\boldsymbol{\mbox{${St}$} = {O}(1)}$
We see in figure 18, in the limit when
$\mbox{${St}$} = {O}(1)$
, that there is minimal difference in the particle trajectories when using either flow. In both cases of
$y_0$
, the particle rapidly escapes the boundary layer before exiting out of the end of the main channel. Therefore, in this limit the boundary layer contribution has minimal effect on any particle which starts outside of the boundary layer and so using the outer flow is suitable in this limit. We note that when a particle starts within the boundary layer, the bounce trajectory does not compare so well, and so the outer solution is not so suitable. We also see similar behaviour when we introduce particles into the full branched channel flow, indicated in black in figure 18 – they too rapidly escape the boundary layer and exit out of the end of the main channel, following a similar path to the particles in the composite solution.
Examples of possible bounces for
$\textit{St} = 1$
, with the flow field given by the outer solution
$\boldsymbol{u}^{(o)}_0 (x,y)$
(blue), the composite solution
$\boldsymbol{u}_c (x,y)$
(red) and the full numerical flow (black) from figure 6. We vary the initial position,
$y_0$
to show the minimal difference between trajectories when
$\mbox{${St}$} ={O}(1)$
, for particles initialise both inside and outside of the boundary layer indicated by black dashed lines. The location of the entrances to the branched channels are indicated by the semicircles on
$y=0$
. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
in the asymptotic solution and
$\mathcal{P}_{\textit{out}} = 0.424$
in the numerical simulation so that
$Q = 1/3$
in both. We also have
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
9.2.2.
$\boldsymbol{\mbox{${St}$} = {O}(\epsilon )}$
We see in figure 19 that the particle trajectories calculated using the outer flow, the composite flow and the numerical solution no longer follow the same path (although they do initially), and there are significantly more bounces than in the
$\mbox{${St}$}={O}(1)$
case. We see that, once a particle enters into the
$\epsilon$
-boundary layer, the Stokes number is small enough such that the particle motion is dominated by the flow and the particle remains within the boundary layer. Thus, in this limit, the composite solution should be used rather than the simpler outer flow. However, the agreement between these two flow fields and also the numerical solution will improve as
$\epsilon \rightarrow 0$
, as expected. Finally, in figure 19, we also observe that the particle trajectories in the composite flow seem disjoint or unresolved around the point sinks, however this behaviour is confirmed by similar behaviour in the full numerical solution of the flow and particles.
Examples of possible bounces for
$\mbox{${St}$} = 0.04$
, with the flow field given by the outer solution
$\boldsymbol{u}^{(o)}_0 (x,y)$
(blue), the composite solution
$\boldsymbol{u}_c (x,y)$
(red) and the full numerical flow (black) from figure 6. We vary the initial position,
$y_0$
to show the minimal difference between trajectories when
$\mbox{${St}$} = {O}(\epsilon )$
, for particles initialise both inside and outside of the boundary layer indicated by black dashed lines. The location of the entrances to the branched channels are indicated by the semicircles on
$y=0$
. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
in the asymptotic solution and
$\mathcal{P}_{\textit{out}} = 0.424$
in the numerical simulation so that
$Q = 1/3$
in both. We also have
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
In the results in § 10, due to the previous observations, we solve for the particle trajectories using the composite flow (7.1) and (7.2) only, to best mimic particles in the full branched channel flow for all values of
$\textit{St}$
.
10. Particle results
In branched channel filters, we are interested in diverting flow through the branched channels whilst retaining particles in the main channel flow. Since we have found an asymptotic prediction for the fluid flow through the branched channels, we now seek a similar result for particles. To quantify this, we define
\begin{align} \mathcal{K} = \frac {\text{Number of particles left through the branched channels}}{\text{Number of particles at the inlet}}. \end{align}
Ignoring particle–particle interactions, we release
$i=19\,999$
particles at slightly perturbed initial points
$(0, y_0^j)$
around
$i$
equispaced internal points (with step size
$1/20\,000$
), for
$j = 1, \ldots , i$
at
$t=0$
. The simulation terminates when all of the particles have left the system. We run 60 such simulations and average the value of
$\mathcal{K}$
over them (noting that averaging over completely random
$y_0$
positions achieves a similar result). We assume that the particles may bounce along the bottom wall,
$y=0$
, until they either leave through the end of the device or down a branched channel. Given that we have taken the limit
$\delta \rightarrow 0$
in the flow problem, we need to impose a ‘removal condition’ that takes account of the finite size of the branched channels (otherwise we will significantly underestimate the number of particles that exit through these channels). We assume that, if a particle collides with the wall within
$\delta \epsilon / 2$
either side of a point sink, we remove the particle from the simulation and assume that the particle will have gone down a branched channel (see figure 20 for an example particle removal via a branched channel). To calculate the proportion of particles,
$\mathcal{K}$
, we count the number of particles that go down a branched channel as described and divide by
$i$
.
Example trajectory in the composite solution for the flow (7.1) and (7.2), when hitting a point sink, or within a
$\delta \epsilon / 2$
radius on
$y=0$
. The point sinks are indicated by point markers on
$y=0$
. Here,
$St = 0.04$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\gamma = 0.5$
,
$Q = 1/3$
and
$y_0 = 0.07$
.
We start by calculating
$\mathcal{K}$
for two distinct limits,
$\mbox{${St}$} = 0$
and
$\mbox{${St}$} = \infty$
. These correspond, respectively, to infinitely light tracer particles and infinitely heavy particles that travel ballistically.
When
$\mbox{${St}$} = 0$
, particles follow the streamlines exactly, exiting the main channel either down one of the branched channels or out of the end of the main channel. No particle will deviate from the flow or bounce on the bottom wall. Hence, we have that
$\mathcal{K} = Q$
.
When
$\mbox{${St}$} = \infty$
, particles retain their initial velocity, and so follow the trajectory given by
\begin{equation} y\,= \,y_0 + \frac {v_c(0, y_0)}{u_c (0, y_0)} x. \end{equation}
Since the outer flow is a good approximation for the flow outside the
$\epsilon$
-boundary layer, the limiting trajectory for a particle, which we define as the dividing line between trajectories for particles that exit through the main channel without bouncing, and those that either exit through the branches or bounce then exit through the end, has the inlet position
$y_0 = \gamma Q / (1+Q)$
(see figure 21). Particles released above this point will never bounce for any
$\textit{St}$
, indicating that a large proportion of particles will simply flow out of the main channel outlet. One might naively think that, since the ratio of total hole size to wall length is
$\delta$
, a
$\delta$
proportion of these ballistic particles will enter the branched channels. This would be true if particles collide with the bottom wall at every
$x \in [0,1]$
, but it slightly overestimates
$\mathcal{K}$
at
$\mbox{${St}$} = \infty$
, since there are locations on the bottom wall close to the inlet which cannot be reached by any ballistic trajectory. We calculate from (7.1)–(7.2) that
\begin{equation} \frac {v_c(0, y_0)}{u_c (0, y_0)} = - \gamma{\kern-1pt}Q \left ( \tanh { \left ( \frac {\pi y_0}{\epsilon } \right )} - \frac {y_0}{\gamma } \right )\!, \end{equation}
Phase diagram for
$\mbox{${St}$} = \infty$
, with the limiting trajectory (red dashed) having inlet position
$\boldsymbol{x}_p = (0, y_0)$
and velocity
$\boldsymbol{u}(0, y_0)$
, for whether a particle will bounce.
which we use with (10.2) to find the intercept position on the bottom wall of a ballistic particle released from
$y_0$
at the inlet, given by
\begin{equation} x_{\textit{intercept}} = \frac {y_0}{\gamma{\kern-1pt}Q \left ( \tanh { \left ( \dfrac {\pi y_0}{\epsilon } \right )} - \dfrac {y_0}{\gamma } \right )}. \end{equation}
We calculate the range of
$y_0$
values which have an
$x$
-intercept value that corresponds with the location of the branched channel (see figure 22). We then divide this total amount by
$\gamma$
to normalise with the inlet height to give us the limiting value of
$\mathcal{K}$
when
$\mbox{${St}$} = \infty$
. For the parameter values chosen, this value is
$0.024$
(to
$3$
decimal places). We note that this value will change depending on the initial velocity of the particles, and is slightly smaller than
$\mathcal{K}=\delta Q / (1 + Q)$
– the value we would have found from using the outer flow. Thus, we see that for small
$\textit{St}$
, the filter allows the same fraction of particles out through the branches as fluid flow, while for large
$\textit{St}$
, the filter retains more particles in the main channel flow, as is preferable.
The
$x$
-intercept, calculated using (10.4), for a ballistic particle given an inlet position
$y_0$
, in the case where
$\mbox{${St}$}=\infty$
. We show the position of the branched channels by black horizontal lines, extended along to the
$y_0$
values for which the
$x$
-intercept falls at the centre of the branched channel. Note that
$y_0 \ne 0$
and so the first
$x$
-intercept is after the first two branched channels.
For a given
$0 \lt \mbox{${St}$} \lt \infty$
, we solve (9.6)–(9.10) numerically using Mathematica, and calculate
$\mathcal{K}$
from these simulations using (10.1). We show the dependence of
$\mathcal{K}$
on
$\textit{St}$
and the limiting behaviour when
$\mbox{${St}$} = 0$
and
$\mbox{${St}$} = \infty$
in figure 23. We see that the number of particles that travel through the branched channels decreases as
$\textit{St}$
increases, since more are deflected back into the main channel flow as we might expect. A similar conclusion is found by Divi et al. (Reference Divi, Strother and Paig-Tran2018), where particle retention in the main channel flow increases for increasing particle size and density. We further see the large
$\textit{St}$
approximation for
$\mathcal{K}$
readily approximates the actual
$\mathcal{K}$
for much smaller values of
$\textit{St}$
than we might have anticipated. The ideal regime is when the number of particles leaving through the branches is small compared with the flux – this allows this filter concept to remove some ‘clean’ water whilst retaining the majority of the particles flowing into the dead-end filter at the end of this simplified domain. However, the lowest fraction of particles that flow down the branched channels in this model is
$0.024$
. We note that the
$\mathcal{K}$
behaviour is non-monotonic, most visibly near
$\mbox{${St}$}=0.026$
. These are not numerical artefacts, but rather a dynamical result of the system. We comment further on this behaviour in Appendix C, where at
$\mbox{${St}$} = 0.26$
we see what we are calling a grazing point, or a type of separatrix, in the solution. But for the purpose of this study, we focus on the overall trend in
$\mathcal{K}$
.
The proportion of particles,
$\mathcal{K}$
, leaving through the main channel compared with the total at the inlet, calculated using (10.1), plotted as black points between (a)
$\mbox{${St}$} \in [0, 0.4]$
with a spacing of
$0.01$
and (b)
$\mbox{${St}$} \in [1,10]$
with a spacing of
$1$
. Here, we have taken
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\gamma = 0.5$
and
$Q = 1/3$
. The limiting values for
$\mathcal{K}$
are shown by red dashed lines. An explanation of the behaviour near
$\mbox{${St}$} = 0.26$
is given in Appendix C.
Finally, we define
\begin{align} \mathcal{R}=\frac {\text{Fraction of particles left through the main channel}}{\text{Fraction of particles left through the branched channels}} = \frac {1 - \mathcal{K}}{\mathcal{K}}. \end{align}
We calculate
\begin{align} \mbox{${St}$} = 0 &: \quad \mathcal{R} = 2, \\[-12pt]\nonumber \end{align}
\begin{align} \mbox{${St}$} = \infty &: \quad \mathcal{R} \approx 40.8. \end{align}
We show how
$\mathcal{R}$
varies with
$\textit{St}$
in figure 24. We see that
$\textit{St}$
varies from
$2$
to
$40.8$
(for
$Q=1/3$
), showing that it is possible to divert a reasonable fraction of the flow through the branched channels whilst retaining a very high proportion of the particles exiting through the end of the main channel, even when
$\textit{St}$
is relatively small. We see that it takes longer for
$\mathcal{R}$
to converge to the large
$\textit{St}$
asymptote than it does for
$\mathcal{K}$
; this is highlighted further in figure 25. Again, we highlight the unusual behaviour near
$\mbox{${St}$} = 0.26$
, and possibly other seemingly non-smooth points.
The proportion of particles,
$\mathcal{R}$
, leaving through the main channel rather than through the branched channels, calculated using (10.5), and plotted as black points between (a)
$\mbox{${St}$} \in [0, 0.4]$
with spacing
$0.01$
and (b)
$\mbox{${St}$} \in [1,10]$
with spacing
$1$
. We plot the corresponding values of
$\mathcal{R}$
for figure 23. Here, we have taken
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\gamma = 0.5$
and
$Q = 1/3$
. The limiting values for
$\mathcal{R}$
are shown by red dashed lines. An explanation of the behaviour near
$\mbox{${St}$} = 0.26$
is given in Appendix C.
The proportion of particles,
$\mathcal{R}$
, leaving through the main channel rather than through the branched channels, calculated using (10.5), and plotted as black points on a larger range of
$\textit{St}$
, on a
$\log$
scale, than figure 24. Here, we have taken
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\gamma = 0.5$
and
$Q = 1/3$
. The limiting value at
$\mbox{${St}$} = \infty$
for
$\mathcal{R}$
is shown by a red dashed line.
11. Conclusions and extensions
In this paper, we consider a high-Reynolds-number flow of a viscous liquid in a branched channel filter, comprising a series of T-junctions of width
$\epsilon$
. We prescribe a constant inlet flux, atmospheric pressure at the bottom of each branched channel and the pressure at the main channel outlet. In addition to the existence of a viscous boundary layer, of thickness
$1/\sqrt {\textit{Re}}$
, on the top and bottom walls of the main channel compartment, there is a layer of thickness
$\epsilon$
on the bottom wall where the flow adjusts to go down the branched channels.
To simplify the full problem, we first consider the flow in each separate branched channel, finding the flux,
$Q_i^{\textit{branch}}$
, as a function of the pressure at the top of the channel and the geometric and flow parameters. We take the limit as the channel width tends to zero and approximate the branch flow as point sinks with strength
$2 Q_i^{\textit{branch}}$
. In the main channel flow, we consider the regime where the branch separation is much larger than the viscous boundary layer, so that
$\epsilon \gg 1/\sqrt {\textit{Re}}$
, and for simplicity, neglect the viscous boundary layers completely. We then decompose the domain into two regions: one close to the wall where the discrete nature of the point sinks are apparent, and a region away from the wall, where the effects of individual sinks are smoothed out – the flow is inviscid in both regions. We find that the leading-order flow away from the wall has constant pressure.
We solve for the flow in an inner region close to
$y=0$
using multiscale asymptotics to capture the rapid oscillations caused by the point sinks. We utilise complex potential flow theory to find an explicit solution for the velocity. We match the solution to the outer flow, deriving an effective boundary condition that smooths out the discrete behaviour due to the branched channels. The dimensionless effective boundary condition is
$v = - \delta ^3 \textit{Re} \mathcal{P}_{\textit{out}}/12 \lambda$
on
$y=0$
, which is dependent on various design and operating parameters, but notably not on the T-junction width,
$\epsilon$
. Using this effective boundary condition, we solve for the leading-order outer flow velocity, and thus build a leading-order composite solution for the flow velocity that holds throughout the main compartment.
We validate our asymptotic solution by comparing with suitable numerical simulations for a large, finite
$\textit{Re}$
. We emulate inviscid main channel flow using a wall slip condition. We find good agreement between the asymptotic solution and numerical solution for an example set of parameters which improves as
$\epsilon$
decreases (i.e. as the number of channels increases). The
$\epsilon $
-boundary layer smooths out the discrete behaviour of the point sinks as required, and the effective boundary condition provides the main flow with the same necessary information as the full branched channel set-up. The explicit solution removes the need for computationally expensive simulations and thus has the potential to massively speed up the design process and optimisation of this device.
We explore how changing the branch angle affects the solution, for a constant branch width and constant hole size. We find that, for a constant branched channel width, which causes the branch hole size to alter with the angle,
$\alpha$
, the effective boundary condition remains constant for any angle provided the hole size,
$\delta \epsilon \sec (\alpha ) \ll \epsilon$
– a similar result is found numerically by Mao et al. (Reference Mao, Bischofberger and Hosoi2024). On the other hand, for a constant hole size, the width of the branch channel changes with angle,
$\alpha$
, and so the effective boundary condition depends on
$\alpha$
. When comparing our asymptotic prediction with the numerical solution, we once again see good agreement.
Having solved for the flow in the simplified main channel domain, we then solve for the trajectories of individual spherical particles placed in the flow, assuming a wall-bounce condition and the Stokes number,
$\textit{St}$
. We explore two distinguished limits between
$\textit{St}$
and
$\epsilon$
and conclude that simulating the particles in the leading-order outer flow does not mimic the effect of the branched channels well when
$\mbox{${St}$} = {O}(\epsilon )$
– we must therefore use the leading-order composite flow. To mimic particles flowing down the branched channels, a particle is removed from the simulation when it hits the bottom wall,
$y=0$
, at the location of a point sink or within a
$\delta \epsilon / 2$
distance away. We suppose that particles can bounce along the bottom wall an infinite number of times until they leave via the main outlet or via the bottom wall removal condition.
We show that for zero Stokes number, the fraction of particles which flow down the branched channels is the same as the flux through the wall,
$Q$
, since the particles follow the flow through the point sinks, and are therefore removed. For larger Stokes number, this limiting number falls to 0.024, slightly below the naive prediction of
$\delta Q/(1+Q)$
as the particles become ballistic, i.e.
$\mbox{${St}$} = \infty$
. We further numerically calculate the fraction of particles which we lose via the bottom wall for
$0 \lt \mbox{${St}$} \lt \infty$
. We find that given these particular rules for spherical particles in this branched channel filter concept, the fraction of particles that are removed through the branched channels is much less than the fluid flux through the wall via point sinks for
$\mbox{${St}$} \gt 0$
. Hence for particles to achieve the best trade-off between maximum fluid flux and minimal particles through the branched channels, spherical particles require a larger
$\textit{St}$
, which could correspond to flocculations of microfibres – a similar conclusion was found by Divi et al. (Reference Divi, Strother and Paig-Tran2018) for plankton particles.
There are various extensions to this work that should be considered. Firstly, it might be useful to study the thus-far neglected viscous boundary layers to see whether they have any influence on the outer flow, rather than neglecting them. This would introduce a non-zero vorticity near the branched channel entrances, which may change the dynamics of a particle’s ricochet trajectory through rotation. Secondly, one may consider the case where the viscous boundary layer is much larger than the
$\epsilon $
-boundary layer, and see how the multiple-scale process may be applied to find a suitable effective boundary condition in this regime. This would be another stepping stone towards our overarching aim of establishing effective boundary conditions for all regimes of laminar high-Reynolds-number flows, in addition to the low-Reynolds-number regime solutions given by Mao et al. (Reference Mao, Bischofberger and Hosoi2024) – which would negate the need for any numerics for the flow problem.
When deriving our effective boundary condition, we also assumed that we are given the outlet pressure,
$\hat {\mathcal{P}}_{\textit{out}}$
, as an operating parameter of the system (prescribed here as a constant, but in reality a function of time found by coupling the device to the downstream dead-end filter). Building and solving a model for the dead-end filter which incorporates clogging and caking of the filter surface, similarly to Köry et al. (Reference Köry, Krupp, Please and Griffiths2021), would provide further insight on this pressure in order to explore its influence on the flow in the ricochet device. Finally, further research might consider more realistic, rounded-corner, geometries for the entrances to the branched channels, as they are in a ricochet device.
There are also a range of extensions to our intentionally simple particle model: firstly, we have thus far only included a drag term in our force balance model and so we should investigate and incorporate other important forces causing this ricochet effect; we ought to identify the correct physics for the bouncing condition by introducing a coefficient of restitution. This might be by considering a smaller microscale structure and studying how the particles ricochet back into the flow before feeding this into our simplified domain; changing the inlet condition for the particles (considering the trajectory shape before entering this device); incorporating particle–particle interactions between these spherical particles; finally, and possibly most importantly, we have been modelling spherical particles in this paper. Although these spherical particles do mimic flocculations of microfibres, individual fibres are ultimately rod-like (and possibly flexible) and so the behaviours and forces would change. Li et al. (Reference Li, Bielinski, Lindner, du Roure and Delmotte2024) consider rod-like particles colliding and pole vaulting over various obstacles. The ricochet separation concept might have improved efficiency with this pole-vaulting effect.
The work done in this paper is a crucial step into understanding the operation of branched channel filters and their effective use in washing machines. Ultimately, it would be valuable to compare our results with those from experiments. Our work has the potential to play a valuable role in optimising Beko’s filtration devices and thus in preventing further microplastic pollution from entering our already heavily polluted oceans.
Acknowledgements
We would like to thank D. Booth, J. Chapman, M. Dalwadi and D. O’Kiely for helpful mathematical discussions, and G. Anderson, K. Pantelidis, S. Tarkuç and M. Toklu for insightful discussions on microfibre filtration at Beko PLC and Arçelik R&D.
Funding
T.F. is grateful to Beko PLC and EPSRC, grant reference number EP/W524311/1, for funding.
Declaration of interests
The authors report no conflict of interest.
Appendix A
As mentioned in §§ 3.1 and 8.3, there is a small region near the top of each branched channel where the flow develops into unidirectional, fully developed, flow. We show how the fully developed flow establishes in figure 26, where we present the velocity,
$v$
, down the channel with centre located at
$x=0.5$
, as a function of
$x$
, at various
$y$
positions. We see that the solution has become fully developed by
$y \approx -0.006$
, which is
${O}(\delta \epsilon )$
for the parameters chosen; the asymptotic prediction for the flow is also shown. We thus conclude that deviations from the asymptotic predictions are confined to an
${O}(\delta \epsilon )$
region near the entrance to the branched channel, and may thus be neglected.
Numerical solution (in solid lines) and asymptotic parabolic prediction (in red dashed lines) for the
$y$
-component of velocity,
$v$
, at the entrance to the centremost branched channel,
$x_i = 0.5$
, from figure 6. We plot various distances into this region to identify the transition length of the numerical solution into fully developed flow. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
Appendix B
In this appendix, we outline the solution procedure taken within this paper and highlight the comparisons we make between numerical simulations and asymptotic predictions, in the ideal and actual comparisons (see figure 27).
Flow solution structure. The full high-Reynolds-number structure is given by (a), applying several assumptions in (b) and (c), to find an effective boundary condition (d). We verify the asymptotics of (d) with the appropriate numerical model in (e), denoted by (
$\star$
). This relates to the original structure (a), via (
$\dagger$
), retaining the initial
$1 \ll \textit{Re} \lt \infty$
everywhere, whilst imposing an outer inviscid flow assumption via a wall slip condition.
Our model derivation begins with high-Reynolds-number flow in a branched domain (see figure 27
a). We explicitly solve for the flow in each branched channel and then simplify the structure by approximating each branched channel by a point sink (see figure 27
b) – the strength of each point sink depends on the finite
$\textit{Re}$
. We now suppose that the flow in the simplified main channel is inviscid everywhere, which allows us, via a multiple-scale analysis, to find an effective boundary condition on the bottom wall (see figure 27
c,d) and composite solutions for the velocities.
Since the viscous boundary layers are tiny in this limit, we emulate the flow in the main channel by imposing a slip condition on the walls as shown in figure 27(e), transforming the problem as indicated by (
$\dagger$
). The comparison between our asymptotic and numerical solution is indicated by (
$\star$
) in figure 27 enabling us to focus on the influence on the flow of the branched channels.
Appendix C
Here, we explore in more detail the trajectories that are used to make figures 23 and 24. We focus on particle-release positions within the
$\epsilon$
-boundary layer and we show in figure 28, as black points, the
$y_0$
locations in each simulation that result in particles going down a branched channel, for a given
$\textit{St}$
. We see that, as we increase
$\textit{St}$
, tongue-like regions of parameter space open up in which particles do not flow down branched channels. This is most striking for
$y_0 \lt 0.01$
. We associate this behaviour with particles hitting the bottom wall very close to the entrance to a branched channel, and we call this process grazing. The amount of grazing appears to increase as
$\textit{St}$
increases up to the final point near
$\mbox{${St}$} = 0.26$
.
The inlet positions
$y_0$
of particles that exit the device through a branched channel versus the Stokes number,
$\textit{St}$
. For each
$\textit{St}$
, we release
$i=19\,999$
particles at slightly perturbed initial points around
$i$
equispaced points, and run
$60$
separate simulations, plotting each point once. Here,
$\mathcal{P}_{\textit{out}} = 0.4$
,
$\textit{Re} = 1000$
,
$\epsilon = 0.04$
,
$\delta = 0.1$
,
$\lambda = 0.1$
and
$\gamma = 0.5$
.
Comparing figures 23 and 24 with figure 28, we notice that the bumps in
$\mathcal{K}$
align with the values of
$\textit{St}$
where there is a sharp switching or bifurcation-like behaviour in the
$y_0$
values that result in particles passing through the branched channels. We see that the large deviation in
$\mathcal{K}$
and
$\mathcal{R}$
at
$\mbox{${St}$} = 0.26$
is aligned with the steepest change in the
$y_0$
values in figure 28. We note that there are a finite number of rapidly changing
$y_0$
regions. In figure 29, we zoom in on the behaviour of
$\mathcal{K}$
and
$\mathcal{R}$
for
$\textit{St}$
near
$0.26$
. We see that, with this increased resolution, there are further non-monotonic regions, and that these also correspond to the regions of rapid
$y_0$
change in figure 28, with the slope here appearing to correlate with the magnitude of the non-monotonic deviations. We hypothesise that the phase-space behaviour discussed in this appendix will depend critically on the precise geometry at the branched channel entrance.
Zoomed in versions of figures 23 and 24 around
$\mbox{${St}$} \in [0.22,0.27]$
.
Open access
