2. A new perspective on the Hele-Shaw approximation

We consider the Navier–Stokes equations for an incompressible, homogeneous fluid confined within a uniform-thickness gap. The equations are non-dimensionalised using the gap height $h$ , so that the vertical coordinate $z$ lies in the interval $-1/2 \lt z \lt 1/2$ . In the Hele-Shaw limit, the in-plane velocity components have a parabolic profile across the gap. The flow can be approximated by $\boldsymbol{u}(x,y,z) \approx (3/2) ( 1-4z^2 )\, \boldsymbol{u}_0(x,y)$ , where $\boldsymbol{u}_{0}$ satisfies the no-slip boundary condition on $\boldsymbol{u}_i$ on any in-plane boundaries. Since $\int_{- {1}/{2}}^{{1}/{2}} (3/2) ( 1-4z^2 )\,{\textrm{d}}z=1$ , $u_0 (x,y)$ is the depth-averaged velocity. Since all velocity components must vanish at the boundaries $z = \pm 1/2$ , a natural choice for the higher-order vertical basis functions is the family of Gegenbauer (ultraspherical) polynomials with parameter $\alpha = -1/2$ ; see Belinsky (Reference Belinsky2000). The Gegenbauer polynomials $C_{n}^{(\alpha )}$ are orthogonal with respect to the weight function $(1-z^{2})^{\alpha - 1/2}$ . In particular, $C_{n}^{(1/2)}$ reduces to the Legendre polynomial of degree $n$ , and for $n \geqslant 2$ , one has $C_{n}^{(-{1}/{2})} (z)= -\int_{-1}^{z} C_{n-1}^{({1}/{2})} (z)\, {\textrm{d}} z$ . The first few polynomials are

(2.1) \begin{equation} C_{2}^{(-{1}/{2})} (z)=\frac {1}{2}\big( 1-z^{2}\big), \quad C_{3}^{(-{1}/{2})} (z)=\frac {1}{2}z\big( 1-z^{2}\big), \quad C_{4}^{(-{1}/{2})} (z)=\frac {1}{8}\big(5 z^2-1\big)\big(1-z^{2}\big). \end{equation}

Since $\lVert C_{n}^{(-{1}/{2})}(z)\rVert _{\infty } \sim n^{- {3}/{2}}$ , we choose the basis as $a_{n}\,C_{n+2}^{(- {1}/{2})} (2z)$ , where $a_{n}$ is a normalisation factor chosen so that the basis remains of order 1 in magnitude. In addition, by the identity $(1-z^{2})\, C_n^{ ({3}/{2})} (z)= (n+1) (n+2)\,C_{n+2}^{(-{1}/{2})} (z)$ , it is therefore equivalent to choose $a_{n}(1-4z^{2})\,C_{n}^{ ({3}/{2})} (2z)$ as the basis, since this differs from our earlier choice only by a multiplicative constant.

Having established an orthogonal family of basis functions that satisfy the no-slip condition at $z=\pm 1/2$ , we now represent the 3-D velocity field in this basis. Specifically, we expand

(2.2) \begin{equation} \boldsymbol{u} = \sum _{n=0}^\infty \boldsymbol{u}_n(x, y)\, a_n\, C_{n+2}^{(-{1}/{2})} (2z), \end{equation}

where $\boldsymbol{u}_n$ satisfies the no-slip boundary condition on any in-plane boundaries. The subsequent task in constructing the approximation is to specify the corresponding in-plane functions $\boldsymbol{u}_n(x,y)$ and derive the equations that they must satisfy. As a first step, we illustrate the procedure by applying the series expansion to a unidirectional shear flow (rectilinear flow).

2.1. Poiseuille flow

We first consider pressure-driven flow in a channel, aligned with the $x$ -axis, with uniform rectangular cross-section $(y,z) \in [-L/2,L/2]\times [-1/2,1/2]$ . Assuming an imposed pressure gradient $G$ , and scaling $p_{x}$ by $G$ , and velocities by $Gh^2/(12\mu )$ , the steady-state solution of the Navier–Stokes equations takes the form $\boldsymbol{u}=(u(y,z),0,0)$ , where $u$ satisfies a 2-D Poisson equation $(\partial _{y}^{2}+ \partial _{z}^{2}) u = - 12$ . In the rectangular channel, this equation admits an exact solution (Cornish Reference Cornish1928; Boussinesq Reference Boussinesq1868). Its series representation is

(2.3) \begin{equation} u_{\textit{exact}}=\frac {3}{2} \big ( 1-4z^{2} \big )-\frac {48}{\pi ^{3}}\sum _{m=0}^{\infty } \frac {(-1)^m}{(1+2m)^{3}}\cos \left ( (2m+1)\pi z \right )\frac {\cosh \left ( (2m+1)\pi y \right ) }{\cosh \left ( (2m+1)\pi {L}/{2} \right )}. \end{equation}

By substituting the leading-order forms of the components of $u\approx ({3}/{2}) (1-4z^{2})\, u_{0} (y)$ into the Poisson equation, we obtain the residual

(2.4) \begin{equation} r (y,z)=\frac {3}{2} \big ( 1-4z^{2} \big ) \partial _{y}^{2}u_{0} (y) -12 u_{0} (y) +12 \approx 0, \end{equation}

where the first two terms come from $\partial _{y}^{2}+ \partial _{z}^{2}$ , and the third term is the dimensionless pressure gradient. Including additional bases such as $C_{3}^{(- {1}/{2})} (z), C_{4}^{(- {1}/{2})} (z), \ldots$ in our expansion would allow us to eliminate the $z$ -dependent terms in (2.4), but at the expense of introducing a coupled system of equations for the unknown functions $u_i$ . Before considering the more complex coupled system, we first seek an optimal closure involving only $u_{0}$ , i.e. a single equation for $u_{0}$ with constant coefficients. Therefore, we assume that the optimal closure is of the form $a\,\partial _{y}^{2}u_{0} (y) -12 u_{0} (y) +12=0$ for some constant $a$ to be determined. (Alternatively, one may assume the more general form $a\, \partial _{y}^{2}u_{0} (y) -b\, u_{0} (y) +c=0$ , but, by applying the residual minimisation procedure introduced later, one would find that $b = 12$ and $c = 12$ , which reduces exactly to the same equation as above.) We can easily solve the equation and obtain an expression for $u_{0} (y)$ :

(2.5) \begin{equation} u_{0} (y)= 1-\text{sech}\left (\frac {\sqrt {3}\, L}{\sqrt {a}}\right ) \cosh \left (\frac {2 \sqrt {3}\, y}{\sqrt {a}}\right )\!. \end{equation}

To determine $a$ , we substitute it into (2.4), obtaining

(2.6) \begin{equation} r (y,z;a) =u^{\prime \prime}_{0}(y;a)\left ( \frac {3}{2} \big ( 1-4z^{2} \big )-a \right )\!. \end{equation}

We choose $a$ to minimise the second factor, integrating the squared residual against a weight function $w(z)$ :

(2.7) \begin{equation} E (a;w) \equiv \int _{-{1}/{2}}^{{1}/{2}} \left ( \frac {3}{2} \big ( 1-4z^{2} \big )-a \right )^{2} w(z)\,{\textrm{d}} z. \end{equation}

Choosing $w(z)=1$ yields $E(a;1)=a^2-2 a+ {6}/{5}$ , which is minimised by $a=1$ , thereby recovering the Darcy–Brinkman approximation, and its solution

(2.8) \begin{equation} u\approx \frac {3}{2} \big(1-4z^{2}\big)\, u_{0} (y)= \frac {3}{2} \big ( 1-4z^{2} \big ) \left ( 1- \cosh \big (2 \sqrt {3}\, y\big )\text{sech}\big (\sqrt {3}\, L\big ) \right )\!. \end{equation}

The key step to producing a better approximation is to observe that the space of functions spanned by Gegenbauer polynomials remains the natural space in which to optimise our approximation of (2.4). Since terms in (2.4) do not vanish at $z=\pm 1/2$ , we multiply the equation by a weight $w (z)= (1-4z^2)$ , ensuring that the product is contained within that space so that $r(z;a)\,w (z) = \sum _{n=0}^{\infty } \beta _n\, C_n^{(- {1}/{2})} (2z)$ , for some set of coefficients $\beta _m$ . Applying Parseval’s theorem to this series expansion gives

(2.9) \begin{equation} \begin{aligned} \kern-12pt E (a;w) &\equiv \!\int _{-1/2}^{1/2}r(z;a)^2\,\big(1-4z^2\big)\,{\textrm{d}}z= \!\int _{-1/2}^{1/2} \big(1-4z^2\big)^{-1} \left ( r(z;a)\,\big(1-4z^2\big) \right )^{2}{\textrm{d}}z\\ & = \sum _{n=0}^{\!\infty } \beta _n^2 \!\int _{-{1}/{2}}^{{1}/{2}}\big(1-4z^2\big)^{-1}\, C_n^{(-{1}/{2})} (2z)^{2}\,{\textrm{d}} z. \end{aligned} \end{equation}

Therefore, our optimal closure consists of minimising (2.7) with $w(z)=(1-4z^2)$ . Under this weight, we obtain $E(a;w)=2a^2/3-8 a/5+36/35$ , which is minimised when $a=6/5$ . This optimal closure leads to a different optimal profile:

(2.10) \begin{equation} u\approx \frac {3}{2} \big(1-4z^{2}\big)\, u_{0} (y) = \frac {3}{2}\big (1-4z^{2} \big ) \left (1-\cosh (\sqrt {10}\,y)\operatorname {sech}\left (L\sqrt {\frac {5}{2}}\right )\right )\! . \end{equation}

We provide two physical justifications for the parabolic weighting. First, the velocity field vanishes at the boundaries and attains its maximum at the centre, so it is natural to assign greater weight near the centre and less near the walls. Second, if the solution is independent of $y$ , then we have

(2.11) \begin{equation} \partial _{z}^{2}u = f, \quad u|_{z=\pm {1}/{2}}=0, \end{equation}

with a general forcing term $f (z)$ . The induced volume flux is then

(2.12) \begin{equation} F=\int _{-{1}/{2}}^{{1}/{2}} u(y)\,{\textrm{d}}z = \int _{-{1}/{2}}^{{1}/{2}} \frac {1}{8} \big ( 1-4z^{2} \big ) f (z)\,{\textrm{d}}z, \end{equation}

where the second equality follows from integration by parts. Thus when the error between the exact solution and its approximation is measured in terms of the resulting volume flux, the weight $1-4z^{2}$ naturally arises.

Both approximations derived above are compared against the exact solution. The choice $a=6/5$ consistently yields a more accurate approximation across the entire $y$ -domain (figure 2 a) and for all tested values of $L$ (figure 2 b). For example, when $L=1$ , the relative errors for $a=1$ and $a=6/5$ are 0.115 and 0.031, respectively. The accuracy of the $a=6/5$ approximation is somewhat surprising, given that $L=1$ does not satisfy the classical scale-separation assumption $L \gg 1$ . As $L$ increases beyond 1, the relative error decreases for both approximations, though the error measured in the $L^{\infty }$ norm does not vanish in the infinite aspect ratio limit. This is because there always exists a narrow region near the side walls $y = \pm L/2$ where the vertical velocity profile deviates from a parabola due to the no-slip boundary condition. However, as the aspect ratio grows, this region occupies a smaller fraction of the domain, outside which the pressure-driven flow in the channel approaches that in a two-parallel-plate geometry. When the error is measured in the $L^2$ norm instead, both approximations exhibit an error that decays as $L^{-1/2}$ in the infinite-aspect-ratio limit, as shown in figure 2(c).

Figure 2.

Validating a 2-D approximation for pressure-driven flow in a rectangular channel. (a) Velocity profile at $z=0$ for $L=1$ : red for (2.8), blue dashed for (2.10), orange dashed for higher-order approximation from (4.7), black for exact solution. The inset shows a zoomed view of the interval $z \in [0,0.1]$ . (b) Relative error in the $L^\infty$ norm $\lVert u - u_{\textit{app}} \rVert _{\infty } / \lVert u \rVert _{\infty }$ as a function of aspect ratio $L$ : red for (2.8), blue dashed for (2.10), orange dashed for (4.7). (c) Relative error in the $L^2$ norm $\lVert u - u_{\textit{app}} \rVert _{2} / \lVert u \rVert _{2}$ as a function of $L$ : red for (2.8), blue dashed for (2.10), orange dashed for (4.7), black dotted for $c L^{1/2}$ , illustrating asymptotic scaling. (d) Relative error in the volume flux as a function of $L$ : red for (2.8), blue dashed for (2.10), orange dashed for (4.7).

Another quantity of interest is the volume flux of the approximate flow solution

(2.13) \begin{equation} \begin{aligned} F &=\int_{-1}^{1}\int_{-{L}/{2}}^{{L}/{2}}u_{0} (y)\,Q_{0} (z)\,{\textrm{d}} y\, {\textrm{d}} z= L-\frac {\sqrt {a} \tanh({\sqrt {3}\, L}/{\sqrt {a}})}{\sqrt {3}}\\ &= L-\frac {\sqrt {a}}{\sqrt {3}}+ {O} ({\textrm e}^{-2\sqrt {3}\,L/\sqrt{a}}), \end{aligned} \end{equation}

which may be compared to the exact flux formula obtained from (2.3):

(2.14) \begin{equation} F=L-\frac {192}{\pi ^5}\sum \limits _{m=0}^{\infty } \frac {\tanh \left ( ({L\pi }/{2}) (2m+1) \right )}{(2m+1)^5}=L-\frac {186\, \zeta (5)}{\pi ^5}+{O} ({\textrm e}^{-L\pi } ), \end{equation}

where $\zeta (s)$ is the Riemann zeta function. Figure 2(d) shows the relative error in the volume flux computed using (2.8) and (2.10), demonstrating that the latter achieves consistently better accuracy. For reference, if side-wall effects are ignored, then the solution reduces to $u(y,z) = ({3}/{2})(1-4z^{2})$ , which yields flux $F=L$ . In terms of the volume flux, (2.13) shows that a channel of length $L$ with no-slip side walls is effectively equivalent to a channel of length $L - \sqrt {2/5}\approx L-0.63246$ with slip side walls. This agrees closely with the asymptotic exact flux, $L-( {186\, \zeta (5)}/{\pi ^5})\approx L-0.63025$ .

To probe the optimality of our closure relation more closely, at each aspect ratio, we numerically compute the value of $a$ that minimises

(2.15) \begin{equation} \frac {\lVert u_{\textit{exact}} (y,z) -\frac {3}{2}\big ( 1-4z^2\big ) u_{0}(y;a)\rVert _{\infty } }{\lVert u_{\textit{exact}} (y,z)\rVert _{\infty }}. \end{equation}

For aspect ratio 1, the value of $a$ is 1.196; for aspect ratio 3, it is 1.188; and for aspect ratio 5, it is 1.184. Hence the value $a=6/5=1.2$ that we obtained by minimising the residual closely approximates the optimal error-minimising values.

We remark that minimising (2.15) requires knowledge of the exact solution to $\Delta u = -12$ . For general or more complicated problems, this is impractical. If the exact solution were available, then constructing such an approximation would be of limited interest. By contrast, the $L^2$ minimisation in (2.7) does not require solving the governing equation in full, nor even computing $u_{0}(y;a)$ . Interestingly, although $u_{0}(y;a)$ depends on $a$ , the value of $a$ obtained by minimising $E(a;w)$ does not necessarily minimise the entire residual $r(y,z;a)$ in the $L^{\infty }$ norm or under other metrics. Nevertheless, it provides a good approximation to the solution. Because this approach yields a reasonably accurate result without requiring the solution of any differential equation, it is computationally attractive.

2.2. Method of weighted residuals

Our closure method follows the general framework of the method of weighted residuals (MWR) (Finlayson Reference Finlayson2013), which has been applied to derive reduced-order approximation for thin-gap fluid flows ranging from free-surface thin-film flows (Ruyer-Quil & Manneville Reference Ruyer-Quil and Manneville2000; Kalliadasis et al. Reference Kalliadasis, Ruyer-Quil, Scheid and Velarde2011) and Hele-Shaw flows (Ruyer-Quil Reference Ruyer-Quil2001) to the thermal convection in a vertical slot (Papanicolaou, Christov & Homsy Reference Papanicolaou, Christov and Homsy2009).

In the standard application of the MWR to obtain a reduced-order model for a differential equation $\mathcal{L}(u) = f$ , one first selects a set of basis functions $\left \{ \phi _n(z) \right \}_{n=1}^N$ and approximates the solution as $u(x,z) \approx \sum _{n=1}^N u_n(y)\, \phi _n(z)$ . Next, the residual is then defined as $r(y,z) = \mathcal{L}(u(y,z)) - f(y,z)$ . Finally, a set of weight functions $\left \{ w_n(z) \right \}_{n=1}^N$ is chosen, and the residual is projected onto these weights via $\int _\varOmega w_n(z)\, r(y,z)\,{\textrm{d}} z = 0$ , $ j = 1, \ldots , N$ , yielding a system of $N$ equations for the unknown coefficients $u_n(y)$ .

We can express both approximations derived in the previous section within the MWR framework. First, we take the weight functions to be $w_{n} = C_{n}^{(1/2)}(2z)$ . Owing to the orthogonality of the Gegenbauer polynomials, imposing $\int _\varOmega w_n(z)\, r(y,z)\,{\textrm{d}} z = 0$ is equivalent to expanding the residual in the $C_{n}^{(1/2)}$ basis and setting the corresponding coefficients to zero. Expanding the residual in (2.4) in terms of $C_{n}^{(1/2)}$ yields

(2.16) \begin{equation} \begin{aligned} r(y,z) & = \big ( \partial _{y}^{2}u_{0}(y) - 12 u_{0}(y) + 12 \big ) C_{0}^{(1/2)}(2z) - \partial _{y}^{2} u_{0}(y)\, C_{2}^{(1/2)}(2z),\\ C_{0}^{(1/2)}(2z) & =1, \quad C_{2}^{(1/2)}(2z)=\frac {1}{2}(12z^{2}-1). \end{aligned} \end{equation}

The coefficient of $C_{0}^{(1/2)}$ recovers the approximation (2.8) with $a = 1$ .

Next, we take $w_{n} =(1-4z^{2})\, C_{n}^{(3/2)}(2z)$ as the weight functions. Expanding the same residual in the $C_{n}^{(3/2)}$ basis yields

(2.17) \begin{equation} \begin{aligned} r(y,z) &= \left ( \frac {6}{5}\partial _{y}^{2} u_{0}(y) - 12 u_{0}(y) + 12 \right ) C_{0}^{(3/2)}(2z) - \frac {1}{5}\partial _{y}^{2} u_{0}(y)\, C_{2}^{(3/2)}(2z),\\ C_{0}^{(3/2)}(2z) & =1, \quad C_{2}^{(3/2)}(2z)=30z^{2}-\frac {3}{2}. \end{aligned} \end{equation}

The coefficient of $C_{0}^{(3/2)}$ corresponds exactly to the approximation (2.10) with $a = 6/5$ .

In this sense, the derivation in the previous subsection is fully consistent with the MWR formulation: the velocity field is expanded in $C_{n}^{(-1/2)}(2z)$ , while the residual is projected onto the space spanned by $C_{n}^{(3/2)}(2z)$ . The resulting procedure may therefore be interpreted as a Petrov–Galerkin method, since the trial and test spaces are constructed from distinct polynomial families.

Although the optimal closure equation can be obtained within the classical MWR framework by selecting appropriate basis functions and weights, our discussion in the previous subsection highlights two non-trivial points. First, we demonstrate an alternative to the third step in the MWR: instead of integrating the residual against multiple orthogonal weights, we propose a closure equation with unknown parameters, substitute it into the residual, and determine the parameters by minimising the weighted $L^{2}$ norm. This formulation provides greater flexibility in the closure form, and may lead to different relations in other contexts. Second, owing to the completeness of the basis, the systems obtained by projecting the residual onto Legendre polynomials $C_{n}^{(1/2)}$ , Gegenbauer polynomials $C_{n}^{(3/2)}$ , or even sine functions (as in Pepper et al. Reference Pepper, Roper, Ryu, Matsudaira and Stone2010), become asymptotically equivalent when many terms are retained in the velocity expansion. However, when only a few terms are used, the choice of weight strongly influences the approximation accuracy. In such cases, one must carefully select both the basis functions and the residual metric to capture the underlying physical behaviour of the system.

2.3. Stokes flow in coaxial flow devices

Under the assumption of Stokes flow, the nonlinear term $\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\!\boldsymbol{u}$ is again negligible. The optimal 2-D approximation replaces the operator $\varDelta$ by $a\,\varDelta _{H}-12$ , where $\varDelta _{H}\equiv \partial _{x}^{2}+ \partial _{y}^{2}$ , even when the flow is 3-D. The same choice of $a$ may be justified (residuals resembling (2.7) are obtained for each component of the fully 3-D Stokes flow). The residual analysis of the continuity equation yields $\partial _{x}u_{0}+\partial _{y}v_{0}=0$ .

We test the approximation on a thin-gap but non-rectilinear flow that is commonly realised in microfluidic devices. In coaxial flow, three fluid streams, flowing sufficiently slowly that ${\textit{Re}}=0$ , are combined and pass through a narrow aperture (figure 3(a); Anna, Bontoux & Stone Reference Anna, Bontoux and Stone2003), emerging as a inner fluid thread sandwiched between the outer fluids. In real coaxial flow devices, the surface tension of the inner fluid thread breaks it up into droplets. We consider a version of the coaxial flow experiment without viscosity contrasts or surface tension, focusing only on the ability of the optimal 2-D model to render the geometry of the interface between the inner and outer fluids. We used COMSOL Multiphysics 5.4 to solve both the 3-D Stokes equations and the 2-D approximate models, and to compare their streamlines with those from the optimal 2-D approximation. Simulations were performed using the geometry from Anna et al. (Reference Anna, Bontoux and Stone2003), while varying the mass flux ratio $Q$ between the inner and outer flows. At a representative value of $Q$ (figure 3 b), the optimal 2-D approximation with $a = 6/5$ closely approximates the interface in both the narrowing and diverging sections of the channel, whereas the classical Hele-Shaw approximation ( $a=0$ ) overpredicts the downstream width of the middle layer, and provides only qualitative agreement in shape.

Figure 3.

Coaxial flow. (a) Simulation geometry based on Anna et al. (Reference Anna, Bontoux and Stone2003). The outer inlets have mass flux $Q$ , the inner inlet has flux $1$ . Inlet and outlet lengths are extended by 20 to ensure fully developed flow. (b) Dividing streamlines for $Q=1.7$ : classical Hele-Shaw (red dashed), derived $a=6/5$ approximation (blue dotted), and 3-D simulation (black dashed).

The agreement between the 2-D approximation and the full 3-D computation in the far downstream region is expected for two reasons. First, the flow becomes rectilinear, resembling Poiseuille flow – a regime in which our optimal 2-D approximation is already known to be accurate (figure 2 b). Second, mass conservation between the inlets and the outlet constrains the position of the interfaces within the downstream Poiseuille flow. Indeed, it is unsurprising that the classical Hele-Shaw approximation will overpredict the width of the inner jet (figure 3 b), because it treats the downstream flow as uniform or plug-like in $y$ , omitting the slowdown near the side walls of the channel that prevents narrowing of the outer flows. However, the optimal 2-D approximation renders the interface accurately even through the narrow orifice separating the inlets and outlet. Although the orifice is relatively narrow (aspect ratio width to depth $\lt 0.4$ ), and conventional wisdom holds that Hele-Shaw-type models are accurate only when the channel is much wider (aspect ratio $\gtrsim 1$ ), the optimal 2-D approximation still performs well. Remarkably, the flow remains approximately 2-D even within the orifice: the ratio of out-of-plane to in-plane velocity components peaks at just 0.22. This modest degree of three-dimensionality, forced by the thin-gap geometry, suggests that the validity of our optimal 2-D model extends well beyond the classical limits typically associated with the Hele-Shaw approximation.

3. Thin-gap flows at moderate Reynolds numbers

After validating the accuracy of the approximation for Stokes flow (i.e. for thin-gap microfluidic devices), we can generalise it to the case of thin-gap flows at moderate Reynolds numbers. These conditions are often realised in inertial microfluidic devices that can be built using the same soft-lithography techniques, with similar dimensions to regular microfluidic devices, but are perfused with high-pressure pumps, to create approximately m s⁻¹ flows (see Di Carlo Reference Di Carlo2009).

For the finite Reynolds number flow, we incorporate an extra term into our closure equation:

(3.1) \begin{equation} \textit{Re}\, ( a\, \partial _{t} \boldsymbol{u}_{\!H} + b \boldsymbol{u}_{\!H}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\!H} \boldsymbol{u}_{\!H})= c\, \varDelta _{H} \boldsymbol{u}_{\!H}- \frac {d}{H^{2}} \boldsymbol{u}_{\!H} - \boldsymbol{\nabla} _{\!H} p, \!\!\quad\!\! \boldsymbol{\nabla} _{\!H} \boldsymbol{\cdot }\boldsymbol{u}_{\!H}=0, \!\!\quad\!\! \left . \boldsymbol{u}_{\!H}\right |_{\partial \varOmega }=0. \end{equation}

Here, $\textit{Re} = {LU}/{\nu }$ is the Reynolds number, with $L$ the characteristic length scale, $U$ the characteristic velocity, and $\nu$ the kinematic viscosity. The choice of $L$ influences the non-dimensionalisation. For example, taking $L={h}/{2}$ (half the channel depth) yields the vertical domain $z \in [-1,1]$ , while taking $L=h$ (the full depth, as in this work) rescales the domain to $z \in [-{1}/{2}, {1}/{2}]$ . To accommodate both conventions, we denote by $H$ the non-dimensionalised gap width.

Following the same weighted residual analysis, using the weight functions $w(z) = 1 -4 z^{2}$ and $w(z) = 1$ yields, respectively,

(3.2a) \begin{align} (a,b,c,d) & =\left (\frac 65, \frac {54}{35}, \frac 65, 12\right )\!, \end{align}

(3.2b) \begin{align} (a,b,c,d) & =\left (1, \frac 65, 1, 12\right )\!. \end{align}

Previous work has reported these and similar closure models, arrived at by different approximations and techniques. For example, Plouraboué & Hinch (Reference Plouraboué and Hinch2002) and Yuan & Azaiez (Reference Yuan and Azaiez2014) used $a = 1$ , $b = 1$ , $c = 1$ and $d = 12$ , while for the case of Stokes flow, Zeng, Yortsos & Salin (Reference Zeng, Yortsos and Salin2003) proposed $a = b = 0$ , $c = {12}/{\pi ^{2}} \approx 1.21$ and $d = 12$ . The closure equation with $w(z)=1$ coefficients in (3.2b ) was derived in Gondret & Rabaud (Reference Gondret and Rabaud1997) by assuming a parabolic flow profile, and averaging the Navier–Stokes equations across the gap.

Ruyer-Quil (Reference Ruyer-Quil2001) first demonstrated that inertial effects give rise to the coefficient ${54}/{35}$ in the depth-averaged equation obtained via the MWR. The present work differs in two important respects. First, the coefficient in Ruyer-Quil (Reference Ruyer-Quil2001) was derived by expanding the velocity field as a sequence of polynomials of the form $z^{j}(1-4z^{2})$ and matching coefficients of like powers of $z$ . This non-orthogonal polynomial basis can hinder the derivation of higher-order approximations, leading to more strongly coupled systems with reduced accuracy. Second, the approximation in Ruyer-Quil (Reference Ruyer-Quil2001) neglects in-plane viscous stresses. Despite this limitation, the model has been widely adopted in studies of fluid instabilities in Hele-Shaw geometries (Chevalier et al. Reference Chevalier, Amar, Bonn and Lindner2006). Martin, Rakotomalala & Salin (Reference Martin, Rakotomalala and Salin2002) add in-plane viscous stresses with an effective viscosity resulting in coefficients identical to those in (3.2a ). More recently, Li et al. (Reference Li, Song, Zhang, Yu and Xuan2021) performed an asymptotic analysis, expanding the pressure and velocity fields in powers of $h/L$ . By retaining terms up to $O(h^2/L^2)$ , discarding terms proportional to the square of the Reynolds number, and making several additional approximation assumptions, they obtained the same 2-D approximation as in (3.2a ). We do not, therefore, make any claim to be the first to arrive at this 2-D approximation, but note that optimal closure via the MWR provides a direct and physically transparent route to deriving the equation, and its higher-order approximations.

To compare weighted and unweighted residuals, we compare both 2-D approximations with 3-D simulations of flow in a centrifuge-on-a-chip (figure 1(c); Hur et al. Reference Hur, Mach and Di Carlo2011), which consists of a rectangular chamber connected to two narrow, elongated channels serving as the inlet and outlet. In our simulations, performed in COMSOL Multiphysics 6.3, the chamber occupies $(x,y,z) \in [-50,50] \times [-50,50] \times [-0.5,0.5]$ , while the inlet and outlet have unit square cross-sections. The inlet region of the flow resembles the channel expansion on which Li et al. (Reference Li, Song, Zhang, Yu and Xuan2021) has already validated the 2-D approximation (3.2a ), but the centrifuge-on-a-chip is a paradigmatic example of the design of inertial microfluidic devices to fractionate particles, such as cells in a blood sample, by size, so centrifuge-on-a-chip flows merit separate attention.

As the Reynolds number $\textit{Re} = Uh/\nu$ of the inlet flow is increased, a separation eddy forms at the centrifuge entrance, and this eddy plays a role in trapping large particles (Paiè et al. Reference Paiè, Che and Di Carlo2017). We measure the eddy size by interpolating the $v=0$ contour to the wall $x=-50$ . The 2-D approximation (3.2a ) (blue diamonds) models the growth of the 3-D eddy with $\textit{Re}$ (black circles) better than the unweighted residual model (red circles) (3.2b ). Though both approximations slightly overpredict the value of $\textit{Re}$ at which separation first occurs (17.5 compared to 16 for the fully 3-D simulation), the unweighted residual approximation underestimates the separation eddy size by a factor of 16–21 % over the entire $\textit{Re}$ range (figure 4), while the error from our 2-D approximation (3.2a ) never exceeds 12 %. Comparing flow fields (figure 4, right-hand panel), we see that the 2-D approximation from (3.2a ) also quantitatively matches the streamlines and velocities obtained in 3-D simulations.

Figure 4.

A 2-D approximation for finite $\textit{Re}$ flows in a centrifuge-on-a-chip geometry: comparison of the separation eddy size in 3-D simulation (black), 2-D approximation (blue, (3.2a )), second-order 2-D approximation (orange, (4.1)) and unweighted residual approximation (red, (3.2b )). The right-hand panel shows streamline patterns coloured by velocity magnitude (red for larger, blue for smaller), comparing the 2-D approximation (top) with the full 3-D simulation (bottom).

4. Higher-order corrections

Although the 2-D approximation (3.2a ) shows good accord with 3-D simulations, there is mild disagreement in flow features such as eddy size. The framework used to derive the approximation naturally lends itself to determining where and how the approximation breaks down. In forming our approximation, we have adopted the Hele-Shaw approximation of assuming a parabolic velocity profile across the gap, and no fluid velocity in the depth direction. Based on figure 4, we focus particularly on dissecting the features of moderate Reynolds number flows that cause these approximations to break down. We first studied whether the assumption of zero vertical velocity breaks down. Across the entire centrifuge-on-a-chip and for every value of $\textit{Re}$ that we simulated, the flow is nearly 2-D: the largest value of $|w|_{\infty }/|\boldsymbol{u}|_{\infty }$ that we observe is only 0.0555, obtained in the chamber outlet.

Next, we examine the validity of the parabolic-velocity assumption. Although a priori estimates for the breakdown of parabolicity can be obtained through linear perturbation analyses of the 3-D Navier–Stokes equations in simple geometries, such results remain limited in scope. For example, Plouraboué & Hinch (Reference Plouraboué and Hinch2002) analysed Kelvin–Helmholtz instability in a Hele-Shaw cell and showed that the velocity profile of an inhomogeneous fluid may lose its parabolic character at Reynolds numbers as low as $\textit{Re} \approx 10$ . Because the fluid in our study is homogeneous, the mechanism responsible for this breakdown does not directly apply, and the relevant Reynolds number threshold may be different.

Our interest is instead in an a posteriori assessment tailored to the Reynolds number and geometry of the microfuge. We quantify deviations from parabolicity by comparing the magnitudes of the first two terms in a Gegenbauer-polynomial expansion of the 3-D velocity field.

Specifically, we interpret the Hele-Shaw approximation as the leading term of an orthogonal polynomial expansion for the velocity. In this, we go beyond reproducing previously derived approximations (Martin et al. Reference Martin, Rakotomalala and Salin2002; Li et al. Reference Li, Song, Zhang, Yu and Xuan2021), generating higher-order corrections on the models for the first time. Including the next-order terms, the velocity field and pressure are approximated by

(4.1) \begin{equation} \begin{aligned} \boldsymbol{u}_{\!H} & =\frac {3}{2} \big (1-4z^2\big ) \left (\boldsymbol{u}_{0,H}+\left (30 z^2-\frac {3}{2}\right ) \boldsymbol{u}_{2,H}\right )\!,\\ w & =15 z \big (1-4z^2\big )^2w_{2}, \quad p=p_{0}+\left (30z^{2}-\frac {3}{2}\right ) p_2,\\ \boldsymbol{f}_{\!H} & =\boldsymbol{f}_{0,H}+\left (30z^{2}-\frac {3}{2}\right ) \boldsymbol{f}_{2,H}, \quad f_{z}=10zf_{2,z}, \end{aligned} \end{equation}

where $\boldsymbol{f}_H$ and $f_{z}$ are the body force terms in the horizontal direction and the vertical direction, included here for generality. In the expansion (4.1), the horizontal velocity component is taken to be even in $z$ , which is the most common situation in microfluidic applications. The expression for $\boldsymbol{u}_{\!H}$ follows from retaining both $C_{2}^{(-1/2)}(2z)$ and $C_{4}^{(-1/2)}(2z)$ , together with appropriate normalisation factors. Due to the no-slip boundary condition and the continuity equation, both $w$ and $\partial _{z}w$ vanish on the boundary, which determines the gapwise basis function that multiplies $w_{2}$ . Moreover, $n\geqslant 4$ , $C_n^{(- {3}/{2})}(2z)={-3}/{(2n-3)} ( C_n^{(- {1}/{2})}(2z)-C_{n-2}^{(- {1}/{2})}(2z) )$ and their derivatives vanish at $z=\pm ( {1}/{2})$ , making them a suitable choice of gapwise basis functions for expanding $w$ . Accordingly, the vertical momentum equation should be projected onto $C_n^{({5}/{2})}(2z)$ . The coefficient of $p_{2}$ is justified by examining the residual of the leading-order approximation (2.6). When $a = 6/5$ , the $z$ -dependent factor becomes

(4.2) \begin{equation} \frac {3}{2} \big (1-4 z^2\big )-a=-\frac {1}{5} \left (30 z^2-\frac {3}{2}\right )=-\frac {1}{5}C_{2}^{({3}/{2})} (2z). \end{equation}

This choice of basis for $p_{2}$ efficiently compensates for the dominant residual generated by the other fields. Furthermore, since $\partial _{z} ( C_n^{({3}/{2})}(z) )=3\,C_{n-1}^{ ({5}/{2})}(z)$ , the derivative $\partial _{z} p$ naturally produces the $C_{n}^{(5/2)}$ family required to balance the error in the vertical momentum equation.

Adopting the expansions in (4.1) and applying the same weighted residual analysis, the continuity equation yields

(4.3) \begin{equation} \begin{aligned} &\boldsymbol{\nabla} _{\!H} \boldsymbol{\cdot }\boldsymbol{u}_{0,H}=0, \quad \boldsymbol{\nabla} _{\!H} \boldsymbol{\cdot }\boldsymbol{u}_{2,H}=\frac {20}{3}w_{2}. \end{aligned} \end{equation}

The momentum equation in the horizontal direction gives

(4.4) \begin{equation} \begin{aligned} &\textit{Re}\! \left ( \frac {6}{5} \partial _{t} \left ( \boldsymbol{u}_{0,H} -\frac {3}{7} \boldsymbol{u}_{2,H} \right ) + \frac {24}{7} \left (\!-\boldsymbol{u}_{0,H}+ \frac {54}{11}\boldsymbol{u}_{2,H} \right ) w_{2} \right .\\ &\quad \left . {}+ \frac {54}{35}\left ( \boldsymbol{u}_{0,H}-\frac {2}{3}\boldsymbol{u}_{2,H} \right )\boldsymbol{\cdot }\boldsymbol{\nabla} _{\!H} \boldsymbol{u}_{0,H}+ \frac {36}{35} \left (\! -\boldsymbol{u}_{0,H}+\frac {24}{11} \boldsymbol{u}_{2,H} \right )\boldsymbol{\cdot }\boldsymbol{\nabla} _{\!H} \boldsymbol{u}_{2,H}\right )\\ &\quad{}= \frac {6}{5} \varDelta _{H} \left ( \boldsymbol{u}_{0,H}- \frac {3}{7} \boldsymbol{u}_{2,H} \right )- 12 \boldsymbol{u}_{0,H} - \boldsymbol{\nabla} _{\!H} p_{0}+ \boldsymbol{f}_{0,H},\\ &\textit{Re}\! \left (\frac {1}{5}\partial _{t} \left ( -\boldsymbol{u}_{0,H}+4 \boldsymbol{u}_{2,H} \right ) -\frac {8}{11}\left ( \boldsymbol{u}_{0,H}+ \frac {18}{13}\boldsymbol{u}_{2,H} \right ) w_{2} \right .\\ &\quad \left . {}+ \frac {2}{5}\left (\! -\boldsymbol{u}_{0,H}+\frac {24}{11}\boldsymbol{u}_{2,H} \right )\boldsymbol{\cdot }\boldsymbol{\nabla} _{\!H} \boldsymbol{u}_{0,H}+ \frac {48}{55} \left ( \boldsymbol{u}_{0,H}-\frac {9}{26} \boldsymbol{u}_{2,H} \right )\boldsymbol{\cdot }\boldsymbol{\nabla} _{\!H} \boldsymbol{u}_{2,H}\right )\\ &\quad{}= \frac {1}{5} \varDelta _{H}\! \left ( -\boldsymbol{u}_{0,H}+4\boldsymbol{u}_{2,H} \right )- 72\boldsymbol{u}_{2,H} - \boldsymbol{\nabla} _{\!H} p_{2}+\boldsymbol{f}_{2,H}, \\ \end{aligned} \end{equation}

and the momentum equation in the vertical direction gives

(4.5) \begin{equation} \begin{aligned} &\textit{Re}\partial _{t} w_{2} + \textit{Re}\,\frac {15}{13}\boldsymbol{u}_{0,H}\boldsymbol{\cdot }\boldsymbol{\nabla} _{\!H} w_{2} = \varDelta _{H} w_{2}- 44w_{2}- \frac {33}{4} p_{2}+\frac {11}{8}f_{2,z}. \\ \end{aligned} \end{equation}

We first apply the above equation to the unidirectional shear flow problem considered in § 2.1. We have $u=({3}/{2} )(1-4z^2 ) (u_{0}+ (30 z^2-({3}/{2}) ) u_{2} )$ , where $u_0$ and $u_{2}$ solve

(4.6) \begin{equation} \begin{aligned} 0 & = \frac {6}{5} \varDelta _{H} u_{0}- \frac {18}{35} \varDelta _{H} u_{2}- 12 u_{0} + 12,\\ 0 & = -\frac {18}{35} \varDelta _{H} u_{0}+ \frac {72}{35} \varDelta _{H} u_{2}- \frac {1296}{7} u_{2}. \\ \end{aligned} \end{equation}

If we set $u_{2}=0$ and keep only the first equation, then the system reduces to the closure equation discussed in § 2.1. By an appropriate change of variables, we can reduce (4.6) to two decoupled equations, leading to the solution

(4.7) \begin{equation} \begin{aligned} u_{i}(y) & = \sum _{j\in \{ 0,2 \}} a_{\textit{0j}}a_{\textit{ij}}\left ( 1-\text{sech}\left (\frac {L}{2 \sqrt {\lambda _j}}\right ) \cosh \left (\frac {y}{\sqrt {\lambda _j}}\right ) \right )\!,\quad i=0, 2,\\ \lambda _{0} & =\frac {1}{36} \left (2+\sqrt {\frac {19}{7}}\right )\approx 0.101,\quad \lambda _{2}=\frac {1}{36} \left (2-\sqrt {\frac {19}{7}}\right ) \approx 0.00979, \\ a_{00} & =- \sqrt {\frac {1}{2}+\frac {4 \sqrt {{7}/{19}}}{5}}\approx -0.992 ,\quad a_{02}= \sqrt {\frac {1}{2}-\frac {4 \sqrt {{7}/{19}}}{5}}\approx 0.120,\\ a_{20} & =\frac {\sqrt {7}}{6 \sqrt {3}} \sqrt {\frac {1}{2}-\frac {4 \sqrt {{7}/{19}}}{5}}\approx 0.0306, \quad a_{22}=\frac {\sqrt {7}}{6 \sqrt {3}} \sqrt {\frac {1}{2}+\frac {4 \sqrt {{7}/{19}}}{5}} \approx 0.253 . \end{aligned} \end{equation}

When $L = 1$ , the $L^{\infty }$ error of the second-order approximation is $0.0036$ , which is approximately one-tenth of the error associated with the leading-order approximation described in § 2.1. The second-order approximation, shown as the orange curve in figure 2(a), is visually indistinguishable from the exact solution. As illustrated in figures 2(b,c), the errors of the second-order approximation, measured in both the $L^{2}$ and $L^{\infty }$ norms, are substantially smaller than those of the leading-order approximation. Moreover, the $L^{\infty }$ norm of the correction term $( {3}/{2})(1 - 4 z^{2})(30 z^{2} - ({3}/{2}))\, u_{2}(y)$ is $0.0249$ , which is of the same order of magnitude as the error in the leading-order approximation. This indicates that the correction term provides a practical error estimate of the leading-order approximation.

The flux of this flow approximation (4.7) is

(4.8) \begin{equation} \begin{aligned} F & =\int_{-1}^{1}\int_{-{L}/{2}}^{{L}/{2}}\frac{3}{2} (1-4z^2) \left(u_{0}+\left(30 z^2-\frac{3}{2}\right) u_{2}\right)\mathrm{d} y\, \mathrm{d} z\\&=\sum_{j\in \{ 0,2 \}} a_{\textit{0j}}a_{\textit{0j}}\left ( L-2 \sqrt {\lambda _{j}}\tanh\left (\frac {L}{2 \sqrt {\lambda _j}}\right ) \right )\\ & = L -\frac {1}{5} \sqrt {\frac {1}{266} \big(9 \sqrt {7}+2618\big)}+ {O} \big({\textrm e}^{-{L}/{\sqrt {\lambda _{0}}}} \big). \end{aligned} \end{equation}

The constant term in the asymptotic expansion (4.8) obtained from the second-order approximation is $0.63029$ , matching the exact flux asymptotic expansion to four digits. In contrast, the leading-order approximation matches to only two digits. When $L = 1$ , the relative error of this flux approximation is $9.47 \times 10^{-5}$ . For comparison, the relative error of the leading-order flux formula (2.13) is $6.68 \times 10^{-3}$ , which is larger by a factor $70.6$ . The error of this second-order flux approximation over a range of aspect ratios $L$ is shown as the orange curve in figure 2(d), demonstrating uniformly higher accuracy across all tested domain sizes.

Next, we apply the second-order approximation (4.4) to the centrifuge-on-a-chip in figure 1(c). Equation (4.7) is used as the boundary condition at the inlet in the simulation of the 2-D approximate equation. To ensure that the flux is 1 at the inlet, we normalise (4.7) using (4.8).

Second-order corrections are especially valuable when the leading-order approximation suffers large errors, or equivalently, when the flow being approximated is highly non-parabolic, or has a cross-gap component. We analyse the centrifuge-on-a-chip to determine where these sources of error occur. In figure 5(a), we plot the magnitudes of $\boldsymbol{u}_{0}$ and $\boldsymbol{u}_{2}$ calculated by projecting 3-D simulation to the basis functions for $\boldsymbol{u}_{0}$ and $\boldsymbol{u}_{2}$ . The relative sizes of those components give a measure of how non-parabolic the flow is. In the top part of figure 5(b), we plot the 3-D simulated $w$ field. Taken together, the 3-D simulations reveal that the leading-order 2-D model becomes inaccurate near the junction between the camber and the outlet. In this region, the gapwise velocity profile departs significantly from the parabolic shape assumed by the leading-order theory, and the cross-gap velocity attains its largest magnitude over the entire domain. These effects indicate the breakdown of the leading-order approximation, and motivate the need for including second-order corrections.

Figure 5.

Simulation of the centrifuge-on-a-chip. The computational domain is $(x,y,z) \in [-5,5] \times [-5,5] \times [-0.5,0.5]$ . The simulation uses Reynolds number 100 and P2–P1 finite elements. (a) Heat maps of the $x$ -components of $\boldsymbol u_{0}$ (top) and $\boldsymbol u_{2}$ (bottom) obtained by projecting the 3-D field onto the corresponding basis functions. (b) Heat maps of $w_{2}$ from the 3-D simulation (top) and from the 2-D reduced model (bottom). (c) Transects of the $x$ -velocity near the outlet ( $x = 4.9$ , $y = 0.3$ ): blue dashed curve shows reconstruction using $\boldsymbol u_{0}$ only; orange dash-dotted curve shows reconstruction using $\boldsymbol u_{0}$ and $\boldsymbol u_{2}$ via (4.1); black solid curve shows the full 3-D result. (d) Transects of $w$ at $x = -1$ , $y = 0.25$ from the 3-D simulation (black) and from the 2-D reduced model (orange dashed). Vertical velocity units are scaled by $10^{-3}$ .

Including the next-order terms alleviates both deficiencies. The bottom part of figure 5(b), computed from the second-order 2-D equations, shows good overall agreement with the full 3-D simulation. Figures 5(c) and 5(d) present the in-plane and cross-gap velocity profiles as functions of $z$ at a representative point in the domain, demonstrating a close quantitative match between the second-order approximation and the 3-D results.

To quantitatively validate the approximation, we compare the produced separation eddy size in figure 4, finding that the extended 2-D approximation captures the eddy size, which is useful for predicting separation performance with almost imperceivable error.

What general benefit can we gain from higher-order approximations? First, in many microfluidic thin-gap flows, the vertical velocity remains small, making truly 3-D flows rare in thin-gap geometries. Conversely, truly 2-D flows occur only in special cases, such as Stokes flow in a channel with constant radius of curvature (Lauga, Stroock & Stone Reference Lauga, Stroock and Stone2004). Some applications, however, deliberately induce strong three-dimensionality, e.g. micromixing (Stroock et al. Reference Stroock, Dertinger, Ajdari, Mezic, Stone and Whitesides2002). The higher-order approximation (4.1) provides a non-trivial description of the vertical velocity component, offering valuable insight into these contexts. Second, higher-order terms provide a built-in mechanism for estimating and controlling approximation error. If the estimated error exceeds a chosen tolerance, then additional terms can be systematically added to improve accuracy.

While our analysis primarily focused on approximating fluid flows, so long as the domain has a thin dimension, the approach is flexible and broadly applicable. It can be readily extended to other physical quantities, such as solute concentration, electric fields and magnetic fields, arising in applications such as reactive solute transport (Salmon & Ajdari Reference Salmon and Ajdari2007), electrokinetic flows (Lin, Storey & Santiago Reference Lin, Storey and Santiago2008; Ding Reference Ding2023) and magnetohydrodynamic microfluidics (Qian & Bau Reference Qian and Bau2009). The method also allows for generalisation of boundary conditions, for example, incorporating no-stress boundaries to model bubble motion (Booth, Griffiths & Howell Reference Booth, Griffiths and Howell2025) or flows in microchannels with soft walls (Inamdar et al. Reference Inamdar, Wang and Christov2020).