2.1. Physical domain

The geometry of the problem is represented in figure 1. The analysis herein develops a solution for the two scenarios represented in figure 1, wherein the channel sidewalls are straight and oriented at a constant angle of convergence or divergence, $2\varTheta$ . Previous studies on converging planar flows (Duryodhan et al. Reference Duryodhan, Singh and Agrawal2013; Singhal & Ansari Reference Singhal and Ansari2016; Goli et al. Reference Goli, Saha and Agrawal2019) have remarked on the difficulty of selection of characteristic length scales, since the width of the channel is not constant for $\varTheta \neq 0$ . To address this difficulty, the approach here develops a solution along an arc with radius of curvature $R$ , where $2R\varTheta$ is the length of the arc bounded by the channel sidewalls, as shown in figure 1(c). Such an arc can be constructed at any arbitrary location within the domain of interest, as long as its radius coincides with lines tangent to the sidewalls at the points of intersection. The mathematical treatment here adopts a cylindrical coordinate system centred at the point of intersection for these two tangent lines, with planar channel boundaries located at $\pm H$ along the $z$ -axis. With these geometric definitions established, $R\varTheta$ and $H$ emerge as the two natural length scales for mathematical analysis. In the discussion that follows, it will be shown that the velocity and pressure fields lend themselves naturally to an asymptotic solution along this arc, using this cylindrical coordinate system.

The three-dimensional geometry of the present problem is reminiscent of Jeffery–Hamel flow, which describes two-dimensional flow emanating from a point source, or converging towards a point sink, at the intersection between two angled walls (Panton Reference Panton2013, Chapter 14.7). Both problems comprise diverging or converging flow in the vicinity of straight sidewalls oriented along a constant angle. However, the wall boundaries at $z=\pm H$ in figure 1, which are not present in the case of Jeffery–Hamel flow, are highly influential for the flow regime and the form of the solution. Beyond the obvious difference in dimensionality, the characteristic length scale for the flow, and therefore, the balance between inertial and viscous forces, is determined significantly by the distance between these wall boundaries. Furthermore, whereas Jeffery–Hamel flow can be solved exactly for arbitrary $\textit{Re}$ , such solutions are not readily extensible to the case shown in figure 1. The present work aims to apply several simplifying assumptions, detailed in the section that follows, that permit an analytical solution for flow regimes that occur regularly in the context of practical microfluidic applications.

2.2. Flow regime and assumptions

The present analysis assumes incompressible, steady flow, and the viscosity of the fluid is assumed to be Newtonian within the range of shear rate achieved within the channel. The velocity field is assumed to be independent of the pressure profiles far upstream and downstream from the arc located at $r=R$ , $-\varTheta \leq \theta \leq \varTheta$ , as illustrated in figure 1(c). The ratio between the channel depth, $H$ , and the sidewall arc length illustrated in figure 1, $R\varTheta$ , is assumed to be small. This ratio is the local inverse aspect ratio of the channel cross-sectional arc, denoted as the small geometric parameter $\alpha \ll 1$ , and is defined as follows:

(2.1) \begin{equation}\alpha =\frac{H}{R\varTheta }.\end{equation}

Additionally, inertial forces are assumed to be negligible compared with viscous forces. Validating this assumption requires a mathematical definition of $\textit{Re}$ . The following is adopted for the present analysis:

(2.2) \begin{equation}\textit{Re}=\frac{2\rho Q/\left(H+R\varTheta \right)}{{\rm \mu} },\end{equation}

where $Q$ is the volumetric flow rate through the channel, $\rho$ is the fluid density and ${\rm \mu}$ is the dynamic viscosity. This definition of $\textit{Re}$ implies the average radial velocity component as the characteristic velocity, and as the characteristic length, the hydraulic diameter of the cross-sectional arc illustrated in figure 1, $D_{h}$ , which is defined as follows:

(2.3) \begin{equation}D_{h}=\frac{2HR\varTheta }{H+R\varTheta }.\end{equation}

For the special case of $\varTheta \rightarrow 0$ and $R\varTheta \rightarrow W$ , which equates to laminar steady flow through a straight rectangular channel of constant width $W$ and height $H$ , the characteristic velocity reverts to the average unidirectional flow velocity. The constant coefficient included in the definition of $\textit{Re}$ ensures that, for this same special case, $D_{h}$ reverts to its conventional definition of hydraulic diameter for a constant rectangular cross-section

(2.4) \begin{equation}\left.D_{h}\right| _{R\varTheta \rightarrow W,\varTheta \rightarrow 0}=\frac{2HW}{H+W}=\frac{4A_{h}}{P_{h}},\end{equation}

where $A_{h}$ is the area, and $P_{h}$ is the wetted perimeter, of the rectangular channel cross-section. Finally, it is worth remarking that, in contrast to many of the computational studies listed in table 1, the definition of ${\textit{Re}}$ given by (2.2) considers only the local cross-section arc, rather than a larger section of the channel along a range of $r$ . The justification for this distinction is the assumption of fully developed flow, which precludes effects sufficiently far upstream and downstream from influencing flow at the cross-section arc.

It is important to appreciate that $\alpha$ scales proportionally to $r^{-1}$ along the length of the channel according to the definition in (2.1). Furthermore, as $\alpha \rightarrow 0$ , $H+R\varTheta \rightarrow R\varTheta$ and thus, according to (2.2), ${\textit{Re}}$ also scales approximately with $r^{-1}$ provided that $\alpha \ll 1$ . In other words, if $\alpha \ll 1$ and ${\textit{Re}}\ll 1$ at the narrowest channel width under consideration, then the assumption of small $\alpha$ and small ${\textit{Re}}$ also holds for wider sections of the channel. It is, of course, possible to conceive various scenarios that fail to satisfy either one or both of these assumptions. However, from a practical standpoint – that is, in a functional microfluidic device – it is not unusual for both of these conditions to be satisfied.

Finally, a point of caution is also warranted regarding the assumption of fully developed flow in relation to the sidewall angle, $\varTheta$ . Practical implementations of the geometry shown in figure 1 presumably require an inlet and outlet somewhere along the channel, each displaced sufficiently from the cross-sectional arc of interest that developing flow effects can be neglected. As $\varTheta \rightarrow \pi /2$ , the inlet or outlet associated with the widening end of the channel requires an increasingly large distance from the cross-section arc in order for the assumption of fully developed flow to apply, unless the pressure field solution of the following analysis is purposefully contrived upstream and downstream of the region of interest by some alternative means.

2.3. Governing equations and boundary conditions

Given the stated assumptions, governing equations for continuity and momentum are, respectively,

(2.5a) \begin{align}&\qquad\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \end{align}

(2.5b) \begin{align}& \textbf{0}=-\boldsymbol{\nabla }p+{\rm \mu} {\nabla} ^{2}\boldsymbol{u}, \end{align}

where $\boldsymbol{u}$ denotes the vector representation of the velocity field, and p is the scalar pressure field.

In the cylindrical coordinate system, the governing equations are expressed as follows:

(2.6a) \begin{align}&\qquad\qquad\quad\frac{1}{r}\frac{\partial }{\partial r}\left(ru_{r}\right)+\frac{1}{r}\frac{\partial }{\partial \theta }\left(u_{\theta }\right)+\frac{\partial }{\partial z}\left(u_{z}\right)=0,\end{align}

(2.6b) \begin{align}&\,\,\frac{\partial p}{\partial r}={\rm \mu} \left[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_{r}}{\partial r}\right)+\frac{1}{r^{2}}\frac{\partial ^{2}u_{r}}{\partial \theta ^{2}}+\frac{\partial ^{2}u_{r}}{\partial z^{2}}-\frac{u_{r}}{r^{2}}-\frac{2}{r^{2}}\frac{\partial u_{\theta }}{\partial \theta }\right]\!,\end{align}

(2.6c) \begin{align}&\frac{1}{r}\frac{\partial p}{\partial \theta }={\rm \mu} \left[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_{\theta }}{\partial r}\right)+\frac{1}{r^{2}}\frac{\partial ^{2}u_{\theta }}{\partial \theta ^{2}}+\frac{\partial ^{2}u_{\theta }}{\partial z^{2}}-\frac{u_{\theta }}{r^{2}}+\frac{2}{r^{2}}\frac{\partial u_{r}}{\partial \theta }\right]\!,\end{align}

(2.6d) \begin{align}&\qquad\qquad\frac{\partial p}{\partial z}={\rm \mu} \left[\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_{z}}{\partial r}\right)+\frac{1}{r^{2}}\frac{\partial ^{2}u_{z}}{\partial \theta ^{2}}+\frac{\partial ^{2}u_{z}}{\partial z^{2}}\right]\!, \end{align}

where u _r , u _θ and u _z are the components of $\boldsymbol{u}$ along unit vectors e _r , e _θ and e _z , respectively.

Asymptotic approximations of the velocity and pressure fields are developed around the small parameter $\alpha$ from (2.1). To represent (2.6a )–(2.6d ) in dimensionless form, the following scaling is used:

(2.7a) \begin{align}&\qquad\qquad\qquad\qquad r=\frac{H}{\alpha }\hat{r}, \theta =\hat{\theta }, z=H\hat{z}, \end{align}

(2.7b) \begin{align}&\left(u_{r},u_{\theta },u_{z}\right)=\left(\frac{Q}{R{\unicode[Arial]{x0398}} H}\hat{u}_{r},\frac{Q}{R{\unicode[Arial]{x0398}} H}\hat{u}_{\theta },\frac{Q}{\left(R{\unicode[Arial]{x0398}} \right)^{2}}\hat{u}_{z}\right)=\left(\frac{\alpha Q}{H^{2}}\right)\left(\hat{u}_{r},\hat{u}_{\theta },\alpha \hat{u}_{z}\right)\!, \end{align}

(2.7c) \begin{align}&\qquad\qquad\qquad\qquad\qquad p=\left(\frac{{\rm \mu} Q}{H^{3}}\right)\hat{p}. \end{align}

Based on (2.7a )–(2.7c ), $H$ , $\alpha Q/H^{2}$ and ${\rm \mu} Q/H^{3}$ are, respectively, the length scale, velocity scale and pressure scale adopted by the present analysis. The dimensionless governing equations in the cylindrical coordinate system are given by (2.8a )–(2.8d ), with hats on all dimensionless variables dropped hereafter for brevity:

(2.8a) \begin{align}&\qquad\qquad\quad\,\,\frac{1}{r}\frac{\partial }{\partial r}\left(ru_{r}\right)+\frac{1}{r}\frac{\partial }{\partial \theta }\left(u_{\theta }\right)+\frac{\partial }{\partial z}\left(u_{z}\right)=0, \end{align}

(2.8b) \begin{align}&\frac{\partial p}{\partial r}=\alpha ^{2}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_{r}}{\partial r}\right)+\alpha ^{2}\frac{1}{r^{2}}\frac{\partial ^{2}u_{r}}{\partial \theta ^{2}}+\frac{\partial ^{2}u_{r}}{\partial z^{2}}-\alpha ^{2}\frac{u_{r}}{r^{2}}-\alpha ^{2}\frac{2}{r^{2}}\frac{\partial u_{\theta }}{\partial \theta }, \end{align}

(2.8c) \begin{align}&\frac{1}{r}\frac{\partial p}{\partial \theta }=\alpha ^{2}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_{\theta }}{\partial r}\right)+\alpha ^{2}\frac{1}{r^{2}}\frac{\partial ^{2}u_{\theta }}{\partial \theta ^{2}}+\frac{\partial ^{2}u_{\theta }}{\partial z^{2}}-\alpha ^{2}\frac{u_{\theta }}{r^{2}}+\alpha ^{2}\frac{2}{r^{2}}\frac{\partial u_{r}}{\partial \theta }, \end{align}

(2.8d) \begin{align}&\qquad\qquad \frac{\partial p}{\partial z}=\alpha ^{4}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_{z}}{\partial r}\right)+\alpha ^{4}\frac{1}{r^{2}}\frac{\partial ^{2}u_{z}}{\partial \theta ^{2}}+\alpha ^{2}\frac{\partial ^{2}u_{z}}{\partial z^{2}}. \end{align}

For the channel shown in figure 1(c), the cross-section arc is located at dimensionless coordinate $r=1/\varTheta$ . The limit in which $\varTheta \rightarrow 0$ represents the special case of flow through a straight rectangular duct, whose velocity field is known (Panton Reference Panton2013, Chapter 11.2). In the dimensionless polar cylindrical coordinate system, the remaining channel walls are located at $z=\pm 1$ and $\theta = {\pm}\varTheta$ . Symmetry of the solution about $z=0$ and $\theta =0$ can be inferred based on the geometry of the channel. The velocity field governed by the equations above is therefore subject to the following boundary conditions:

(2.9a) \begin{align}&\qquad\qquad\qquad \left.u_{r}\right| _{\theta =\varTheta }=\left.u_{z}\right| _{\theta =\varTheta }=\left.u_{r}\right| _{z=1}=\left.u_{\theta }\right| _{z=1}=0, \end{align}

(2.9b) \begin{align}&\qquad\qquad\qquad\qquad\qquad \left.u_{\theta }\right| _{\theta =\varTheta }= \left.u_{z}\right| _{z=1}=0, \end{align}

(2.9c) \begin{align}&\left.\frac{\partial u_{r}}{\partial \theta }\right| _{\theta =0}=\left.u_{\theta }\right| _{\theta =0}=\left.\frac{\partial u_{z}}{\partial \theta }\right| _{\theta =0}=\left.\frac{\partial u_{r}}{\partial z}\right| _{z=0}=\left.\frac{\partial u_{\theta }}{\partial z}\right| _{z=0}=\left.u_{z}\right| _{z=0}=0, \end{align}

where (2.9a ) represents no-slip boundary conditions, (2.9b ) represents no-flux boundary conditions and (2.9c ) represents symmetry about $\theta =0$ and $z=0$ ,

The velocity field is also subject to the following flow rate condition:

(2.9d)

where $-1$ on the right-hand side of (2.9d ) corresponds to flow through a planar channel with converging sidewalls, and $1$ corresponds to diverging sidewalls.

Finally, the pressure field must be evaluated relative to its value at a reference point, which may be located outside the cross-section arc. This is formalised as the following:

(2.9e) \begin{equation}\left.p\right| _{\left(r,\theta , z\right)=\left(r_{p},0,0\right)}=P,\end{equation}

where $P$ denotes the pressure at point $(r,\theta ,z)=(r_{P},0,0)$ .

2.4. Asymptotic expansion representation

Setting $\alpha =0$ reduces the order of (2.8a )–(2.8d ) such that the $\theta$ -axis boundary conditions cannot all be satisfied simultaneously. Therefore, this analysis seeks a composite expansion incorporating an inner velocity field $\boldsymbol{V}\equiv (V_{r}, V_{\theta } ,V_{z})$ , which is valid near $\theta =\varTheta$ , and an outer velocity field $\boldsymbol{U}\equiv (U_{r}, U_{\theta } ,U_{z})$ , which is valid far from the $\theta$ -axis boundaries. In physical terms, $\boldsymbol{V}$ represents the boundary layers of the flow at the sidewalls of the channel. $\boldsymbol{U}$ , which represents outer flow beyond the influence of these boundary layers, is nevertheless influenced by the varying channel width by virtue of continuity.

For the outer velocity solution, an asymptotic expansion of the following form is assumed:

(2.10) \begin{equation} (U_{r}, U_{\theta },U_{z},p)= (U_{r,0}, U_{\theta ,0} ,U_{z,0},p_{0})+\alpha (U_{r,1}, U_{\theta ,1} ,U_{z,1},p_{1})+O(\alpha ^{2}).\end{equation}

The notation $O(\alpha ^{n})$ denotes terms that are of order $\alpha ^{n}$ . Substituting $(u_{r}, u_{\theta },u_{z},p)=(U_{r}, U_{\theta },U_{z},p)$ in (2.8a )–(2.8d ) and separating terms of the same order into separate equations yields the following asymptotic expansion formulation of the governing equations for the outer velocity:

(2.11a) \begin{align}&\frac{1}{r}\frac{\partial }{\partial r}(rU_{r,n})+\frac{1}{r}\frac{\partial }{\partial \theta }(U_{\theta ,n})+\frac{\partial }{\partial z}(U_{z,n})=0, \end{align}

(2.11b) \begin{align}&\qquad\qquad \qquad\,\, \frac{\partial p_{n}}{\partial r}=\frac{\partial ^{2}U_{r,n}}{\partial z^{2}}, \end{align}

(2.11c) \begin{align}&\,\,\,\,\quad\qquad\qquad \frac{1}{r}\frac{\partial p_{n}}{\partial \theta }=\frac{\partial ^{2}U_{\theta ,n}}{\partial z^{2}}, \end{align}

(2.11d) \begin{align}&\quad\qquad\qquad\quad\,\,\, \frac{\partial p_{n}}{\partial z}=0, \end{align}

subject to the following boundary conditions:

(2.12a) \begin{align}\left.U_{r,n}\right| _{z=1}=\left.U_{\theta ,n}\right| _{z=1}&=0, \end{align}

(2.12b) \begin{align} \left.U_{z,n}\right| _{z=1} & =0, \end{align}

(2.12c) \begin{align}\left.\frac{\partial U_{\theta ,n}}{\partial z}\right| _{z=0}=\left.\frac{\partial U_{r,n}}{\partial z}\right| _{z=0}=\left.U_{z,n}\right| _{z=0} & =0, \end{align}

(2.12d) \begin{align}\left.p_{0}\right| _{\left(r,\theta , z\right)=\left(r_{p},0,0\right)}& =P, \end{align}

(2.12e) \begin{align} \left.p_{1}\right| _{\left(r,\theta , z\right)=\left(r_{p},0,0\right)} & =0, \end{align}

where (2.12a ) represents no slip, (2.12b ) represents no flux, (2.12c ) represents symmetry about $z=0$ and (2.12d ) and (2.12e ) express pressure at a known reference point. These equations are valid for $n=0$ , which produces the $O(1)$ governing equations, and $n=1$ , which produces the $O(\alpha )$ governing equations, for the outer region.

To study the velocity field near $\theta =\varTheta$ , a scaled inner variable $\xi =(\theta -\varTheta )/\alpha$ is introduced, and an asymptotic expansion of the velocity field of the following form is assumed:

(2.13) \begin{equation}(V_{r}, V_{\theta } ,V_{z},p)=(V_{r,0}, V_{\theta ,0} ,V_{z,0},p_{0})+\alpha (V_{r,1}, V_{\theta ,1} ,V_{z,1},p_{1})+O(\alpha ^{2}). \end{equation}

Here, it is worth noting that the governing equations, expressed generally as (2.8a )–(2.8d ), and for the outer regions as (2.11a )–(2.11d ), have so far contained only powers of $\alpha ^{2}$ among their constant coefficients. This might suggest $\alpha ^{2}$ as the appropriate choice of small parameter for the asymptotic expansion, were it not for the coefficient powers of  $\alpha$ that appear when the governing equations are expressed in terms of the scaled inner variable. The inner variable scaling elevates the order of terms containing derivatives of  $\theta$ in (2.8a )–(2.8d ), and even produces two terms of $O(\alpha ^{-1})$ , which warrants the following order $O(\alpha ^{-1})$ inner equations for continuity and $\theta$ -momentum, respectively:

(2.14a) \begin{align}\frac{1}{r}\frac{\partial }{\partial \xi }(V_{\theta ,0})& =0, \end{align}

(2.14b) \begin{align}\frac{1}{r}\frac{\partial p_{0}}{\partial \xi }& =0. \end{align}

The $O(\alpha ^{-1})$ expressions for conservation of $r$ -momentum and $z$ -momentum are both trivial.

Inner expansion terms of $O(1)$ and $O(\alpha )$ can be organised under the following expression of (2.8a )–(2.8d ):

(2.15a) \begin{align}&\frac{1}{r}\frac{\partial }{\partial r}(rV_{r,n})+\frac{1}{r}\frac{\partial }{\partial \xi }(V_{\theta ,n+1})+\frac{\partial }{\partial z}(V_{z,n})=0, \end{align}

(2.15b) \begin{align}&\qquad\frac{\partial p_{n}}{\partial r}=\frac{1}{r^{2}}\frac{\partial ^{2}V_{r,n}}{\partial \xi ^{2}}+\frac{\partial ^{2}V_{r,n}}{\partial z^{2}}-\frac{2}{r^{2}}\frac{\partial V_{\theta ,n-1}}{\partial \xi }, \end{align}

(2.15c) \begin{align}&\,\,\frac{1}{r}\frac{\partial p_{n+1}}{\partial \xi }=\frac{1}{r^{2}}\frac{\partial ^{2}V_{\theta ,n}}{\partial \xi ^{2}}+\frac{\partial ^{2}V_{\theta ,n}}{\partial z^{2}}+\frac{2}{r^{2}}\frac{\partial V_{r,n-1}}{\partial \xi }, \end{align}

(2.15d) \begin{align}&\qquad\frac{\partial p_{n}}{\partial z}=0, \end{align}

where $n=0$ produces the $O(1)$ governing equations, and $n=1$ produces the $O(\alpha )$ governing equations, for the inner region.

Here, a minor issue arises. The $z$ -velocity components appear in exactly one governing equation, (2.15a ), which expresses continuity. However, the $\theta$ -velocity component which appears in the same equation, also fails to appear in the momentum equations of the same order. As a result, the $O(\alpha ^{n})$ governing equations, as formulated above, afford fewer equations than needed to resolve the unknown velocity components. This issue can be remedied by resorting to the $O(\alpha ^{n+1})$ conservation of $\theta$ -momentum and $(\alpha ^{n+2})$ conservation of $z$ -momentum equations. This pair of higher-order equations can be combined into a new equation which contains $V_{\theta ,n+1}$ and $V_{z,n}$ , and crucially, none of the other velocity components. This additional set of equations is expressed as follows:

(2.16) \begin{equation}\frac{1}{r^{2}}\frac{\partial ^{3}V_{\theta ,n+1}}{\partial \xi ^{2}\partial z}+\frac{\partial ^{3}V_{\theta ,n+1}}{\partial z^{3}}+\frac{2}{r^{2}}\frac{\partial ^{2}V_{r,n}}{\partial \xi \partial z}=\frac{1}{r^{3}}\frac{\partial ^{3}V_{z,n}}{\partial \xi ^{3}}+\frac{1}{r}\frac{\partial ^{3}V_{z,n}}{\partial \xi \partial z^{2}}.\end{equation}

In a physical sense, (2.16) reflects the requirement that pressure field terms and corresponding velocity field terms must be self-consistent. Although the derivation of (2.16) requires resorting to higher-order governing equations, this approach does not require explicit definitions for the higher-order terms that appear in those equations, and is therefore is permissible according to standard perturbation methods. Detailed intermediate steps in the derivation of (2.16) are provided in Appendix A.

The velocity and pressure in the inner region are subject to the following boundary conditions:

(2.17a) \begin{align}&\left.V_{r,n}\right| _{\xi =0}=\left.V_{z,n}\right| _{\xi =0}=\left.V_{r,n}\right| _{z=1}= \left.V_{\theta ,n+1}\right| _{z=1}=0, \end{align}

(2.17b) \begin{align}&\qquad\qquad\quad\left.V_{\theta ,n+1}\right| _{\xi =0}= \left.V_{z,n}\right| _{z=1}=0, \end{align}

(2.17c) \begin{align}&\qquad\qquad\left.\frac{\partial V_{r,n}}{\partial z}\right| _{z=0}=\left.\frac{\partial V_{\theta ,n+1}}{\partial z}\right| _{z=0}=0, \end{align}

where (2.17a ) expresses no slip, (2.17b ) expresses no flux and (2.17c ) expresses symmetry about $z=0$ .

2.5. Computational simulation

Validation of the analytical solution of the converging slot domain is performed using an implementation of the semi-implicit method pressure linked equation (SIMPLE) solver. Originally developed by Patankar and Spalding (Patankar & Spalding Reference Patankar and Spalding1972; Patankar Reference Patankar1980), SIMPLE is a finite volume formulation of the problem that utilises a staggered grid, with pressure computed at cell centres, and velocity computed at boundary facets. In this method, a pressure guess is used to generate a corresponding velocity field; the pressure is corrected to account for the local residuals of continuity; and the corrected pressure and a corresponding velocity correction are used to generate the next iteration in each successive improvement to solution accuracy. The influence of the pressure field on the velocity is ‘semi-implicit’ in the sense that, at each mesh element, contributions from neighbouring elements are accounted for in the calculation of velocity from pressure, and ignored in the calculation of a pressure correction from velocity. The termination condition for SIMPLE is based on the error calculation for continuity, in which these neighbour contributions are not ignored, thus ensuring accuracy of the converged solution.

The majority of assumptions outlined in § 2.2 are also used for the computational model, summarised as follows:

1. (i) Incompressible, steady flow.

2. (ii) Newtonian viscosity.

3. (iii) Negligible inertial terms (i.e. ${\textit{Re}}\ll 1$ ).

4. (iv) $\alpha \ll 1$ .

In place of the final assumption used in the analytical formulation – that of fully developed flow – the pressure states at the upstream and downstream domain boundaries of the computational model are specified directly. To facilitate comparison between the analytical and computational results, the dimensionless representation of continuity from (2.8a ), and the dimensionless representation of conservation of momentum from (2.8b )–(2.8d ), are used as the governing equations for the model.

The validation study of this work uses a structured cylindrical mesh, spanning from $r=1/2\varTheta$ to $r=3/2\varTheta$ , from $\theta =0$ to $\theta =\varTheta$ and from $z=0$ to $z=1$ for each case, with a scaling factor of $\alpha ^{-1}$ along the $\theta$ -dimension applied to the $\theta =\varTheta$ boundary cells. For the transition between the scaled boundary cells and interior cells, a growth ratio of 1.3 for cell width along $\theta$ is used. The purpose of this scaling factor is to resolve the velocity field in the near vicinity of the sidewall, in sufficient detail to enable comparison with the analytical solution. In total, six distinct meshes, summarised table 2, are considered. An example of the mesh domain is plotted in figure 2, for the first case tabulated in table 2.

Table 2. Summary of mesh parameters for each validation case reported in this study.

Figure 2. The mesh used for case 1 of table 2, with $\alpha =0.2$ and $\varTheta =\pi /4$ , visualised in (a) isometric view and (b) from the top down.

The velocity boundary conditions used for the computational study are the same as those for the analytical formulation, as expressed in (2.9a )–(2.9c ). However, rather than assuming fully developed flow within the computational domain, the following pressure conditions are used for the inlet and outlet boundaries:

(2.18a) \begin{align}&\left.p\right| _{r=1/2\varTheta }=\left.\left(p_{0}+\alpha p_{1}\right)\right| _{r=1/2\varTheta }\!, \end{align}

(2.18b) \begin{align}&\left.p\right| _{r=3/2\varTheta }=\left.\left(p_{0}+\alpha p_{1}\right)\right| _{r=3/2\varTheta }\!, \end{align}

where $p_{0}$ and $p_{1}$ represent, respectively, the $O(1)$ and order $O(\alpha )$ terms in the asymptotic expansion defined in (2.10), whose solution form is derived later in § 3.3 of this work. Equation (2.18a ) expresses the pressure condition at the channel outlet (for converging flow) and (2.18b ) expresses the pressure condition at the channel inlet (again, for converging flow). These pressure conditions are adopted to facilitate direct comparison between the analytical solution and the computational validation results.

3.1. Outer solution

The following form of the $O(1)$ and $O(\alpha )$ pressure field expansions can be found to satisfy the governing equations expressed as (2.11a ) and the boundary conditions given by (2.12a ):

(3.19) \begin{align}p_{n}-P_{n}&=C_{n,0}\ln \left(\frac{r}{r_{P}}\right)+\sum _{m=1}^{\infty }\left(-C_{n,m,1}r^{-m}+C_{n,m,2}r^{m}\right)\cos (m\theta )\nonumber\\&\quad -\!\sum _{m=1}^{\infty }\!\left(-C_{n,m,1}{r_{P}}^{-m}+C_{n,m,2}{r_{P}}^{m}\right)+\!\sum _{m=1}^{\infty }\!\left(-C_{n,m,3}r^{-m}\!+C_{n,m,4}r^{m}\right)\sin (m\theta ),\end{align}

where $P_{0}=P$ and $P_{1}=0$ express the reference pressure condition, and $C_{n,0}$ , $ C_{n,m,1}$ , $C_{n,m,2}$ , $ C_{n,m,3}$ and $ C_{n,m,4}$ are constants to be determined from the flow rate condition and by matching between inner and outer solutions.

Taking partial derivatives of (3.19), substituting into (2.11a ) and twice integrating with respect to $z$ gives the following form of the outer velocity field:

(3.20a) \begin{align} U_{r,n}& = -\frac{1}{r}\left\{C_{n,0}+\sum _{m=1}^{\infty }m\left(C_{n,m,1}r^{-m}+C_{n,m,2}r^{m}\right)\cos (m\theta )\right .\nonumber\\&\quad \left .+\sum _{m=1}^{\infty }m\left(C_{n,m,3}r^{-m}+C_{n,m,4}r^{m}\right)\sin (m\theta )\right\}\left(\frac{1-z^{2}}{2}\right)\!, \end{align}

(3.20b) \begin{align} U_{\theta ,n} &= \frac{1}{r}\left\{\sum _{m=1}^{\infty }m\left(-C_{n,m,1}r^{-m}+C_{n,m,2}r^{m}\right)\sin (m\theta )\right .\nonumber\\& \quad \left .+\sum _{m=1}^{\infty }n\left(-C_{n,m,3}r^{-m}+C_{n,m,4}r^{m}\right)\cos (m\theta )\right\}\left(\frac{1-z^{2}}{2}\right)\!, \end{align}

(3.20c) \begin{align} U_{z,n} &= 0. \end{align}

3.2. Inner solution

Inner velocity terms must be determined starting with the $O(\alpha ^{-1})$ contributions that appear in (2.14). In conjunction with the no-flux boundary condition at $\xi =0$ , (2.14a ) requires that $V_{\theta ,0}=0$ . Equation (2.14b ) implies, helpfully, that the $O(1)$ pressure terms are inherited from the outer flow. Furthermore, a re-examination of (2.15c ) with these results reveals that the $O(\alpha )$ pressure terms are inherited from the outer region as well.

Taking a separation of variables approach, the following solution form can be found to satisfy (2.15b ) and the boundary conditions in (2.17) for $V_{r,0}$ and $V_{r,1}$ :

(3.21a) \begin{equation}V_{r,n}=\left(\left.\frac{\partial p_{n}}{\partial r}\right| _{\theta =\varTheta }\right)\left\{\sum _{m=1}^{\infty }E_{m}\exp \left\{\delta _{m}r\xi \right\}\cos \left(\delta _{m}z\right)-\left(\frac{1-z^{2}}{2}\right)\right\}\!,\end{equation}

where $n=0$ corresponds to the $O(1)$ solution, $n=1$ corresponds to the $O(\alpha )$ solution, $\delta _{m}=({1}/{2})(2m-1)\pi$ and $E_{m}=(-1)^{m-1}2/\delta _{m}^{3}$ .

The remaining non-zero velocity components are $V_{\theta ,1}$ , $V_{z,0}$ and $V_{z,1}$ , which appear in (2.15a ) and (2.16). These equations are satisfied by the following solution forms of the velocity $\theta$ - and $z$ -components:

(3.21b) \begin{align}V_{\theta ,1}&=-\frac{1}{r}\left(\left.\frac{\partial p_{0}}{\partial r}\right| _{\theta =\varTheta }\right)\sum _{m=1}^{\infty }\left\{E_{m}r\xi \exp \left\{\delta _{m}r\xi \right\}\cos \left(\delta _{m}z\right)\vphantom{\frac{1}{\delta _{m}^{4}}}\right.\nonumber\\&\quad\left . +\left(1-\exp \left\{\delta _{m}r\xi \right\}\right)\left[FE_{m}\cos \left(\delta _{m}z\right)-\left(F-\frac{1}{\delta _{m}}\right)\frac{1}{\delta _{m}^{4}}g'(z)\right]\right\}\!, \end{align}

(3.21c) \begin{align}V_{z,n}&=-\frac{1}{r}\left(\left.\frac{\partial p_{n}}{\partial r}\right| _{\theta =\varTheta }\right)\sum _{m=1}^{\infty }\left(F-\frac{1}{\delta _{m}}\right)\exp \left\{\delta _{m}r\xi \right\}\left[E_{n}\sin \left(\delta _{m}z\right)-\frac{1}{\delta _{m}^{3}}g(z)\right]\!. \end{align}

The coefficients $\delta _{m}$ , $E_{m}$ , $F$ and $g(z)$ that appear in (3.21a )–(3.21c ) are defined as follows:

(3.22a) \begin{align}\delta _{m} & = \frac{1}{2}\left(2m-1\right)\pi, \end{align}

(3.22b) \begin{align} E_{m}& = \frac{\left(-1\right)^{m-1}2}{\delta _{m}^{3}}, \end{align}

(3.22c) \begin{align} F & = \sum _{m=1}^{\infty }\left(\frac{6}{\delta _{m}^{5}}\right), \end{align}

(3.22d) \begin{align}g(z) & = \frac{\sum _{m=1}^{\infty }\left[\left(F\delta _{m}-1\right)E_{m}\sin \left(\delta _{m}z\right)/\delta _{m}\right]}{\sum _{m=1}^{\infty }\left[\left(F\delta _{m}-1\right)/\delta _{m}^{4}\right]}. \end{align}

Intermediate steps used to arrive at these solution forms are provided in Appendix A.

3.3. Composite solution

To construct a composite solution, matching between the outer solutions and each of the inner solutions is required. Matching between $\boldsymbol{U}$ and $\boldsymbol{V}$ can be accomplished using van Dyke’s rule, which asserts

(3.23) \begin{equation}\lim _{\alpha \rightarrow 0} \left\{\boldsymbol{U}\left(\theta =\alpha \xi +\varTheta \right)\right\}=\lim _{\alpha \rightarrow 0} \left\{\boldsymbol{V}\left(\xi =\left(\theta -\varTheta \right)/\alpha \right)\right\}\!.\end{equation}

Matching the $O(1)$ outer and inner solutions requires that constant coefficients $C_{0,m,1}$ , $C_{0,m,2}$ , $C_{0,m,3}$ and $C_{0,m,4}$ that appear in (3.19) and (3.20a )–(3.20c ) must be zero for all $m$ . Similarly, matching the $O(\alpha )$ solutions requires that $C_{1,m,1}=0$ for all $m\neq 1$ , $C_{1,1,1}=[C_{0,0}/\sin (\varTheta )]\sum _{n=1}^{\infty }(6/\delta _{n}^{5})$ and $C_{1,m,2}=C_{1,m,3}=C_{1,m,4}=0$ for all $m$ . Alternatively, these same results can be produced by matching by means of an intermediate variable. The constant coefficients $C_{0,0}$ and $C_{0,1}$ , which remain unknown following the matching procedure, must be determined based on the flow rate condition.

After matching, a composite velocity can be constructed from the superposition of inner and outer solutions, subtracting the redundant contributions from the inner and outer velocity components, called the overlap. Here, it is necessary to include inner velocity contributions at both sidewall boundaries: the boundary at $\theta =\varTheta$ , which corresponds to $\boldsymbol{V}$ , and the boundary at $\theta =-\varTheta$ , which corresponds to $\boldsymbol{V}$ reflected about the plane of symmetry at $\theta =0$ . Similarly, it is necessary to subtract the overlap corresponding to each boundary. The overlap corresponding to the intermediate region near $\theta =\varTheta$ is defined as $\lim _{\alpha \rightarrow 0} \{\boldsymbol{U}(\theta =\alpha \xi +\varTheta )\}$ , or equivalently, $\lim _{\alpha \rightarrow 0} \{\boldsymbol{V}(\xi =(\theta -\varTheta )/\alpha )\}$ , while the overlap corresponding to the opposite sidewall boundary can be defined as the same overlap reflected about $\theta =0$ . After subtracting the overlap and collecting terms, the composite solutions can be expressed as follows:

(3.24a) \begin{align}& p=P+ (C_{0,0}+\alpha C_{1,0})\ln \left(\frac{r}{r_{P}}\right)-\frac{\alpha C_{0,0}F}{\sin (\varTheta )}\left(\frac{\cos \left(\theta \right)}{r}-\frac{1}{r_{p}}\right)+O(\alpha ^{2}), \end{align}

(3.24b) \begin{align}& u_{r}=-\frac{1}{r}\left[C_{0,0}+\alpha C_{1,0}+\alpha \frac{C_{0,0}F}{r}\frac{\cos \left(\theta \right)}{\sin (\varTheta )}\right]\left(\frac{1-z^{2}}{2}\right)+\frac{1}{r}\left[C_{0,0}+\alpha C_{1,0}\vphantom{\frac{1}{r}}\right .\nonumber\\&\qquad\left .+\alpha \frac{C_{0,0}F}{r}\cot (\varTheta )\right]\left[2\sum _{n=1}^{\infty }E_{n}\exp \left\{-\delta _{n}\frac{r\varTheta }{\alpha }\right\}\cosh \left\{\delta _{n}\frac{r\theta }{\alpha }\right\}\cos \left(\delta _{n}z\right)\right]+O(\alpha ^{2}), \end{align}

(3.24c) \begin{align}& u_{\theta }=-\frac{\alpha }{r^{2}}C_{0,0}F\frac{\sin (\theta)}{\sin (\varTheta )}\left(\frac{1-z^{2}}{2}\right)+\frac{2}{r^{2}}C_{0,0}\sum _{n=1}^{\infty }\left\{\exp \left\{-\delta _{n}\frac{r\varTheta }{\alpha }\right\}\left(r\varTheta \sinh \left\{\delta _{n}\frac{r\theta }{\alpha }\right\}\right . \right .\nonumber\\&\qquad\left .\left . -\,r\theta \cosh \left\{\delta _{n}\frac{r\theta }{\alpha }\right\}\right)\left[E_{n}\cos \left(\delta _{n}z\right)\right]\right\}+\frac{2\alpha }{r^{2}}C_{0,0}\sum _{n=1}^{\infty }\left\{\exp \left\{-\delta _{n}\frac{r\varTheta }{\alpha }\right\}\sinh \left\{\delta _{n}\frac{r\theta }{\alpha }\right\}\right .\nonumber\\&\qquad\times \left .\left[FE_{n}\cos \left(\delta _{n}z\right)-\left(F-\frac{1}{\delta _{n}}\right)\frac{1}{\delta _{n}^{4}}g'(z)\right]\right\}+O(\alpha ^{2}), \end{align}

(3.24d) \begin{align}& u_{z}=-\frac{2}{r^{2}}\left(C_{0,0}+\alpha C_{1,0}+\alpha \frac{C_{0,0}F}{r}\cot (\varTheta )\right)\boldsymbol{\cdot }\sum _{n=1}^{\infty }\exp \left\{-\delta _{n}\frac{r\varTheta }{\alpha }\right\}\cosh \left\{\delta _{n}\frac{r\theta }{\alpha }\right\}\nonumber\\&\qquad\times\left(F-\frac{1}{\delta _{n}}\right)\left[E_{n}\sin \left(\delta _{n}z\right)-\frac{1}{\delta _{n}^{3}}g(z)\right]+O(\alpha ^{2}). \end{align}

Finally, the flow rate condition from (2.9d ) must be applied to determine the constant coefficients $C_{0,0}$ and $C_{1,0}$ . This flow rate condition can be expressed in terms of the asymptotic expansion as follows:

(3.25)

where $-1+O(\alpha ^{2})$ on the right-hand side of (3.25) denotes flow through a planar channel with converging sidewalls, and $1+O(\alpha ^{2})$ denotes flow through a planar channel with diverging sidewalls. Evaluating this integral using the expressions in (3.24a )–(3.24d ) provides the following definitions for $C_{0,0}$ and $C_{1,0}$ :

(3.26a) \begin{align}&\pm C_{0,0}=\frac{3}{4\varTheta }\left[1-\alpha \varTheta \cot (\varTheta )F\right]^{-1}, \end{align}

(3.26b) \begin{align}&\qquad \pm (C_{0,0}+\alpha C_{1,0})=\frac{3}{4\varTheta }, \end{align}

where a positive sign on the left-hand sides of (3.26a )–(3.26b ) denotes converging sidewalls, and a negative sign denotes diverging sidewalls. Intermediate steps in the application of the flow rate condition are provided in Appendix B.

With these definitions, the pressure field and pressure gradients along the arc located at $r=1/\varTheta$ can be expressed as follows:

(3.27a) \begin{align} \pm p & =P-\frac{3}{4\varTheta }\left[\ln \left(r_{P}\varTheta \right)+\beta \left(\cos \left(\theta \right)-\frac{1}{r_{p}\varTheta }\right)\right]\!, \end{align}

(3.27b) \begin{align} \pm \frac{\partial p}{\partial r}& =\frac{3}{4}\left[1+\beta \cos \left(\theta \right)\right] \end{align}

(3.27c) \begin{align} \pm \frac{1}{r}\frac{\partial p}{\partial \theta }& =\frac{3}{4}\beta \sin \left(\theta \right)\!, \end{align}

(3.27d) \begin{align}\pm \frac{\partial p}{\partial z}&=0, \end{align}

where $P$ is a known reference pressure at $r=r_{P}$ , and the signs of the left-hand sides of (3.27a )–(3.27c ) are positive for channels with converging sidewalls and negative for channels with diverging sidewalls. Additionally, the parameter $\beta$ expresses the magnitude of variation of the pressure gradient along $\theta$ , in terms of $\varTheta$ and $\alpha$ , and is defined as follows:

(3.28) \begin{equation}\beta =\frac{\alpha \varTheta F\csc (\varTheta )}{1-\alpha \varTheta F\cot (\varTheta )}.\end{equation}

After substituting these definitions into (3.24a )–(3.24c ), the $O(\alpha )$ expansion for the velocity field along the arc located at $r=1/\varTheta$ can be expressed as follows:

(3.29a) \begin{align}\pm u_{r}&=-\frac{3}{4}\left[1+\beta \cos \left(\theta \right)\right]\left(\frac{1-z^{2}}{2}\right)+\frac{3}{4}\left[1+\beta \cos (\varTheta )\right]S_{r1}+O(\alpha ^{2}), \end{align}

(3.29b) \begin{align}\pm u_{\theta }&=-\frac{3}{4}\beta \sin \left(\theta \right)\left(\frac{1-z^{2}}{2}\right)+\frac{3}{4}\beta \left(\frac{\sin (\varTheta )}{F}\right)\left(S_{\theta 1}+S_{\theta 2}\right)+O(\alpha ^{2}), \end{align}

(3.29c) \begin{align} \pm u_{z}&=-\frac{3}{4}\varTheta \left[1+\beta \cos (\varTheta )\right]S_{z1}+O(\alpha ^{2}), \end{align}

where the signs of the left-hand sides of (3.29a )–(3.29c ) are positive for channels with converging sidewalls and negative for channels with diverging sidewalls, and $S_{r1}$ , $S_{r2}$ , $S_{\theta 1}$ , $S_{z1}$ and $S_{z2}$ are functions involving infinite summations, defined as follows:

(3.30a) \begin{align} S_{r1}\left(\theta ,z\right)& =2\sum _{m=1}^{\infty }E_{m}\exp \{-\delta _{m}\alpha ^{-1}\}\cosh \left\{\frac{\delta _{m}}{\alpha }\frac{\theta }{\varTheta }\right\}\cos \left(\delta _{m}z\right), \end{align}

(3.30b) \begin{align}S_{\theta 1}\left(\theta ,z\right)& =2\sum _{m=1}^{\infty }\!\left\{\frac{1}{\alpha }\exp \! \left\{-\frac{\delta _{m}}{\alpha }\right\}\!\left(\!\sinh \!\left\{\frac{\delta _{m}}{\alpha }\frac{\theta }{\varTheta }\right\}-\frac{\theta }{\varTheta }\cosh \! \left\{\frac{\delta _{m}}{\alpha }\frac{\theta }{\varTheta }\right\} \right)\left[\!\frac{E_{m}}{F}\cos \!\left(\delta _{m}z\right)\right]\right\}\!, \end{align}

(3.30c) \begin{align}S_{\theta 2}\left(\theta ,z\right)& =2\sum _{m=1}^{\infty }\left\{\exp \left\{-\frac{\delta _{m}}{\alpha }\right\}\sinh \left\{\frac{\delta _{m}}{\alpha }\frac{\theta }{\varTheta }\right\}\left[E_{m}\cos \left(\delta _{m}z\right)-\left(\frac{1}{\delta _{m}^{4}}-\frac{1}{F\delta _{m}^{5}}\right)g'(z)\right]\right\}\!, \end{align}

(3.30d) \begin{align}S_{z1}\left(\theta ,z\right)& =2\sum _{m=1}^{\infty }\exp \left\{-\frac{\delta _{m}}{\alpha }\right\}\cosh \left\{\frac{\delta _{m}}{\alpha }\frac{\theta }{\varTheta }\right\}\left(F-\frac{1}{\delta _{m}}\right)\left[E_{m}\sin \left(\delta _{m}z\right)-\frac{1}{\delta _{m}^{3}}g(z)\right]\!, \end{align}

where the coefficients $\delta _{m}$ , $E_{m}$ , $F$ and $g(z)$ that appear in (3.30a )–(3.30d ) follow the same definitions as previously given in (3.22a )–(3.22d ).

It is worth remarking that the infinite summation terms in (3.30a )–(3.30d ) can be computed numerically and truncated to reasonable accuracy after only a few terms in each series.

3.4. Composite solution plots

The behaviour of the $O(\alpha )$ composite expansion in the intermediate region, where the solution from (3.29a ) and (3.30) transitions between inner and outer approximations, is illustrated in figure 3. For three illustrative values of $\alpha$ , figure 3 overlays the inner, outer and composite expansions within the domain bounded by $-3\leq \xi \leq 0$ . This range for the scaled coordinate $\xi$ is sufficiently large to encompass significant contributions from both the inner and outer velocity expansions. Respectively, figures 3(a)–3(c) each comprise separate subplots for each velocity component, alongside a graphical illustration of the plot range along $\xi$ , the scaled inner $\theta$ coordinate.

Figure 3. Matching between inner and outer velocity field solutions plotted along the intermediate region for (a) $\alpha =0.2$ , (b) $\alpha =0.1$ and (c) $\alpha =0.05$ . The channel cross-section geometry is plotted on the leftmost sub-figure in each row of sub-figures, with a dashed arc representing the portion of the channel along which the velocity components are plotted. The $r$ -velocity components are plotted in sub-figure (i) of each row, $\theta$ -velocity components are plotted in sub-figure (ii) of each row, and non-zero $z$ -velocity components are plotted in sub-figure (iii) of each row.

A few notable features of the perturbation solution can be observed in these plots. First, the size of the boundary layer at the sidewall scales with $\alpha$ , which offers further confirmation for the choice of scaled inner variable. This is also reflected in the mathematical forms of the inner velocity components, whose largest terms are of order $O(\alpha )$ , and which diminish to zero as $\xi \rightarrow 0$ . Second, $u_{r}$ is the dominant velocity component within the channel, particularly for small $\alpha$ , which underscores the convenience of the cylindrical coordinate system for analysing the flow problem. This aspect of the problem formulation neatly delineates between a primary flow component and the other, secondary, flow components. Third, among the three velocity components, the $\theta$ -velocity exhibits the most diverse and interesting behaviour in the transition between the inner and outer regions. In this sense, the $\theta$ -velocity in the intermediate region can be thought of as reconciling between an outer radial flow, and a more complex three-dimensional flow near the sidewall boundaries, where the radial component exemplifies the outer flow regime, and the $z$ -component exemplifies the inner flow regime.

The form of the $O(\alpha )$ perturbation solution for pressure is plotted in figure 4. The radial and circumferential gradients are shown in figures 4(a) and 4(b) respectively, while the gradient along $z$ is not plotted, since it is zero to mathematical order $O(\alpha )$ . These gradients are cosine and sine curves centred at $\theta =0$ , which follows directly from their mathematical expressions in (3.27b ) and (3.27c ). The highest possible peak amplitude across the channel arc, for either pressure gradient component, is $\beta$ , and is achieved as $\varTheta \rightarrow \pi /2$ .

Figure 4. Pressure gradients in a converging channel. Radial and circumferential components are plotted against $\theta$ in (a) and (b), respectively for several values of $\beta$ . The value of $\beta$ , the magnitude of gradient variation along $\theta$ , is plotted versus $\varTheta$ for various $\alpha$ in (c), and versus $\alpha$ for various $\varTheta$ in (d).

The value of $\beta$ , as defined in (3.28), is determined entirely by the channel geometry, and it is therefore worth remarking on the relationship between $\alpha$ , $\varTheta$ and the variation in pressure across the channel. Figure 4(c) plots the relationship between $\beta$ and $\alpha$ , which exhibits a power-law scaling for fixed $\varTheta$ , while figure 4(d) presents a weak, nonlinear relationship between $\beta$ and $\varTheta$ for various values of $\alpha$ . However, this does not necessarily equate to a weak influence of $\varTheta$ on the variation of pressure across the channel; here, it is important to recall that $\varTheta$ defines the valid bounds of $\theta$ in figures 4(a) and 4(b), and the pressure variation across the channel is therefore influenced by $\varTheta$ quite strongly. The value of $\alpha$ is also influential on shape of the pressure field, in the sense that pressure variation along the cross-section arc diminishes with $\alpha$ according to a power-law relationship.

3.5. Streamlines

Figure 5 plots the streamlines for the $O(\alpha )$ perturbation solution for converging flow, averaged through the channel depth ( $z)$ . Three illustrative cases are represented here: $\varTheta =\pi /12$ in figure 5(a) corresponds to a channel encroaching upon the special case of a straight rectangular cross-section, $\varTheta =\pi /4$ in figure 5(b) represents a case of moderately large $\varTheta$ and $\varTheta =7\pi /16$ in figure 5(c) approaches the special case of planar flow at an abrupt contraction. For all three cases, $\alpha =0.2$ .

Figure 5. Streamlines of the $O(\alpha )$ perturbation solution for velocity, averaged across $z$ , in the converging planar channel from $r=1/(2\varTheta )$ to $r=3/(2\varTheta )$ , with $\alpha =0.2$ . Streamlines spanning the entire channel width ( $-\varTheta \lt \theta \lt \varTheta$ ) are plotted in (a)–(c). Streamlines near the sidewall at $\theta =\varTheta$ are plotted in (d)–(f), with a log scale used for the scaled inner coordinate $\xi$ . The dotted curved lines in (a)–(c) denote the vertical axis range for each corresponding plot among (d)–(f). Three separate cases of sidewall angle are represented: (a), (d) $\varTheta =\pi /12$ , (b), (e) $\varTheta =\pi /4$ and (c), (f) $\varTheta =7\pi /16$ .

At first glance, the shapes of these streamlines are straightforward throughout the majority of the physical domain, exhibiting some measure of complexity only in very near proximity to the channel sidewalls, and closely conforming to the overall channel convergence angle everywhere else. Similar plots generated with smaller $\alpha$ exhibit streamline shapes that conform even more strongly to the sidewall geometry. This conformal behaviour is well known among Hele-Shaw flows generally, and therefore expected as $\alpha \rightarrow 0$ . Near the sidewalls, the streamlines curve gradually inward along the channel circumference, along $\theta$ , as flow converges. Only in the nearest vicinity of the sidewall do the streamline shapes appear to reverse this trend; however, this range of proximity to the wall is well below $O(\alpha ^{2})$ , which, arguably, puts it beyond the stated resolution of the perturbation solution as originally formulated.

3.6. Pressure and velocity contours

The order $O(\alpha )$ pressure field solution is plotted in figure 6 for various permutations of channel geometry. Here, the influence of $\alpha$ can be seen by comparing across the rows of panels, and the influence of $\varTheta$ , can be seen by comparing across columns. The dotted line in each contour plot denotes the cross-section arc at $r\varTheta =1$ , and serves as a useful point of reference to judge visually the distortion of the pressure field. This distortion is, perhaps not surprisingly, most extreme for large channel sidewall angles. Comparing the pressure plots across each row, it can be seen that the effect of decreasing $\alpha$ is to drive the pressure field towards conformance with the overall channel geometry, such that the pressure contours become orthogonal to the streamline shapes plotted in figure 5. In figure 6(i), $\alpha =0.05$ is sufficiently small to ensure this conformance even for a planar channel geometry approaching that of abrupt contraction.

Figure 6. Contour plots of the $O(\alpha )$ perturbation solution for pressure plotted across the physical domain from $r=1/(2\varTheta )$ to $r=3/(2\varTheta )$ for nine different cases: (a)–(c) $\alpha =0.2$ , (d)–(e) $\alpha =0.1$ and (g)–(i) $\alpha =0.05$ ; (a), (d), (g) $\varTheta =\pi /6$ , (b), (e), (h) $\varTheta =\pi /3$ and (c), (f), (i) $\varTheta =7\pi /16$ .

Figure 7. Contour plots of the $O(\alpha )$ perturbation solution for velocity. The value of $u_{r}$ is plotted for $\alpha =0.2$ , $\alpha =0.1$ and $\alpha =0.05$ in (a), (b) and (c), respectively. The value of $u_{\theta }$ is plotted for $\alpha =0.2$ , $\alpha =0.1$ and $\alpha =0.05$ in (d), (e) and (f), respectively. Lastly, $u_{z}$ is plotted for $\alpha =0.20$ , $\alpha =0.1$ and $\alpha =0.05$ in (g), (h) and (i), respectively.

Similar visualisations of the $O(\alpha )$ solution can also be constructed for the velocity components. These are shown in figure 7 for three representative values of $\alpha$ , and in figure 8(a–c) for three representative values of $\varTheta$ . All cases represented in these plots share several important features. First, as observed in the matching visualisations from figure 3, the dominant velocity component is oriented along $r$ , with secondary flow along $\theta$ in both the inner and outer regions. Second, the domain of influence for $u_{z}$ is highly localised to the regions near the channel sidewalls, whereas $u_{\theta }$ spans both inner and outer regions. And third, the primary and secondary flow components differ by roughly one to two orders of magnitude, which ensures that, for the outer regions of the physical domain, most flow behaviour is captured by the primary flow component, $u_{r}$ , alone.

Figure 8. The $r$ -velocity contours at $\alpha =0.1$ compared for (a) $\varTheta =\pi /6$ , (b) $\varTheta =\pi /3$ and (c) $\varTheta =7\pi /16$ . (d) The value of $\beta$ , the magnitude of variation of the pressure gradient along $\theta$ , normalised against its minimum possible value for each $\alpha$ and plotted against $\varTheta$ , with normalised $\beta$ values corresponding to panels (a)–(c) indicated. (e) Value of 1+ $\beta$ plotted without normalisation against the full range of possible $\varTheta$ .

The influence of channel inverse aspect ratio, $\alpha$ , on the velocity field can be seen by comparing the three cases represented in figure 7. Decreasing $\alpha$ across figure 7(a–c) – that is, approaching pure Hele-Shaw flow in the converging channel – drives $u_{r}$ toward a two-dimensional parabolic profile along $z$ . For $u_{\theta }$ , decreasing $\alpha$ across figure 7(d–f) has the effect of shifting the outer edges of velocity contours toward the channel sidewalls, with the interior contour edges gradually tapering towards zero at the centre of the channel. Finally, comparing figure 7(g–i) shows that, as $\alpha \rightarrow 0$ , the velocity contribution of $u_{z}$ is increasingly relegated to the near vicinity of the sidewalls, while the magnitude of $u_{\theta }$ diminishes across the entire width of the channel.

In contrast to $\alpha$ , varying $\varTheta$ influences the shapes of the velocity contour plots with relative subtlety. Figure 8(a–c) plots contours of $u_{r}$ across the channel cross-section arc for three representative cases of $\varTheta$ , with $\alpha =0.1$ for all three cases. The lack of variation among the contour shapes in figure 8(a–c) reflects the form of the mathematical expression from (3.29a )–(3.29c ), wherein $u_{r}$ , $u_{\theta }$ and $u_{z}$ all scale in some manner proportionally with $\beta$ , which encapsulates the geometry of the channel. Therefore, here it is helpful to scrutinise the relationship between $\varTheta$ and $\beta$ , in order to understand the influence of $\varTheta$ on the $O(\alpha )$ velocity field solution. This relationship is plotted in figure 8(d–e), with $\beta$ normalised to its minimum value as $\varTheta \rightarrow 0$ in figure 8(d), and the 1+ $\beta$ scaling value plotted without normalisation in figure 8(e). The general point conveyed by this visualisation is that the direct influence of the value of $\varTheta$ on $\beta$ is weak. Even for large $\alpha$ , the range of variation for $\beta$ across $\varTheta$ falls within one half of the minimum value of $\beta$ for the same value of $\alpha$ . In turn, the narrow range of possible $\beta$ for a given value of $\alpha$ is reflected as minimal variation among the velocity contour plots of figure 8(a–c).

3.7. Physical interpretation

Here, some additional remarks are warranted regarding the physical meaning of the results shown in figures 4–8. First, much of the behaviour evident in visualisations of the pressure and velocity solutions can be discerned directly from the mathematical forms of the $O(\alpha )$ solutions themselves, (3.29a )–(3.29c ). These expressions convey that the pressure field is two-dimensional to mathematical order $\alpha$ ; that the velocity field has one primary flow component oriented along $r$ , which is constant along the cross-section arc to leading order; that the velocity field also comprises two secondary flow components, one of which, $u_{z}$ , is of $O(\alpha )$ and confined to the regions in near proximity of the channel sidewalls; and that continuity balance within the inner region necessitates an additional $O(\alpha )\ \theta$ -velocity component, which matches to an $O(\alpha )$ contribution in the outer region of the channel.

The dimensionless parameters describing the channel geometry are $\alpha$ and $\varTheta$ , and the visualisations of §§ 3.5 through 3.6 elucidate several aspects of their influence on both the pressure and velocity fields. Decreasing $\alpha$ drives the solution towards the case of pure Hele-Shaw flow, where flow in the channel is parabolic along $z$ and oriented along the primary flow direction, $r$ . In the context of the planar converging or diverging channel, this effectively orients the pressure field orthogonally to the chosen cylindrical coordinate system, and forces the pressure field towards a two-dimensional configuration, which the perturbation model predicts accurately to order $\alpha$ . Additionally, the relative importance of the two secondary flow components diminishes as $\alpha \rightarrow 0$ , in manners distinct to each component: $u_{\theta }$ diminishes in magnitude, while $u_{z}$ diminishes in terms of the size of its domain of influence relative to the full span of the channel.

With respect to increasing $\varTheta$ , the perturbation analysis predicts a warping of the pressure field that becomes increasingly significant as $\alpha \rightarrow 1$ . The velocity field solution, however, appears largely insensitive to $\varTheta$ when plotted across the channel cross-section arc in terms of dimensionless cylindrical coordinates. In physical terms, this insensitivity reflects the fact that the majority of variation in velocity behaviour effected by $\varTheta$ is captured by the problem formulation itself. The constituent parts of the problem formulation – in particular, the orientation of the coordinate system, and the dimensional scaling of the domain – describe the influence of $\varTheta$ accurately, at least to mathematical order $O(\alpha )$ .

4.1. Case results

Here, numerical results derived using the SIMPLE methodology described in § 2.5 are presented as a basis for validation of the perturbation analysis results from § 3. The preceding discussion has structured many observations in terms of the relationships between pressure and velocity, and the geometric parameters $\alpha$ and $\varTheta$ . Therefore, the computational results are organised into one set of cases spanning a range of $\alpha$ with constant $\varTheta$ , shown in figure 9, and a second set of cases spanning a range of $\varTheta$ with constant $\alpha$ , shown in figure 10. The selection of these cases is tailored towards illustrating where the shortcomings of the perturbation analysis become apparent, while still highlighting its accuracy to the expected mathematical order.

Figure 9. Pressure and velocity results of the numerical simulation for three illustrative cases of converging planar channel flow. (a)–(c) Overlays of the streamlines from $r=1/(2\varTheta )$ to $r=3/(2\varTheta )$ overlaid on contour plots of the pressure field. Contours of velocity along $r=1/\varTheta$ are shown for (d)–(f) $u_{r}$ , (g)–(i) $u_{\theta }$ and (j)–(l) $u_{z}$ . Each case is organised as a column in the grid of panels: $\alpha =0.2$ is shown in (a), (d), (g), (j); $\alpha =0.1$ is shown in (b), (e), (h), (k); and $\alpha =0.05$ is shown in (c), (f), (i), (l). Here, $\varTheta =\pi /4$ in all three cases.

Figure 10. Pressure and velocity results of the numerical simulation for three illustrative cases of converging planar channel flow. (a)–(c) Overlays of the streamlines from $r=1/(2\varTheta )$ to $r=3/(2\varTheta )$ overlaid on contour plots of the pressure field. Contours of velocity along $r=1/\varTheta$ are shown for (d)–(f) $u_{r}$ , (g)–(i) $u_{\theta }$ and (j)–(l) $u_{z}$ . Each case is organised as a column in the grid of panels: $\varTheta =\pi /6$ is shown in (a), (d), (g), (j); $\varTheta =\pi /3$ is shown in (b), (e), (h), (k); and $\varTheta =7\pi /16$ is shown in (c), (f), (i), (l). Here, $\alpha =0.1$ in all three cases.

Comparing the streamline and pressure plots for the three cases in figure 9(a–c), the computational model confirms both the expected streamline shapes and pressure contour results from the discussions around figures 5 and 6, respectively. The primary velocity component, $u_{r}$ , in figure 9(d–f), and as well as the secondary velocity $u_{\theta }$ in figure 9(g–i) are also in close agreement with the analytical results shown in figure 7(a–f). The results for $u_{z}$ in figure 9(j–l), however, illustrate well how the expected $O(\alpha )$ accuracy of the perturbation analysis can become consequential. All three results, when compared with the identical cases computed using the perturbation solution in figure 7(g–i), suggest that the analytical solution slightly under-predicts the size of the domain of influence for $u_{z}$ , although the general shape of this velocity component near the sidewalls is the same for both computational and analytical results. Furthermore, for the case of $\alpha =0.2$ , the numerical simulation predicts significant $u_{z}$ contributions in the outer regions of the channel domain. These differences do not appear to be an artefact of the mesh selection, and repeating the simulations with as little as half the mesh grid count used in this study produces numerical results that are visually identical to those in figure 9.

The validation results in figure 10(a–c), similarly, agree with the pressure contour plots in figure 6(d–f), and depict streamline shapes conforming closely to the converging geometry of the channel, also in agreement with the perturbation analysis results. The velocity plots, however, do exhibit significant differences relative to the perturbation solution results in figures 7 and 8. The contour shapes for $u_{r}$ are influenced noticeably by the choice of $\varTheta$ , becoming more focused inward towards the centre of the channel with increasing $\varTheta$ . Here, $u_{\theta }$ does not differ remarkably in shape between the three cases, but does differ in magnitude, commensurate with the corresponding differences in $u_{r}$ . The differences among $u_{z}$ are perhaps the most illustrative, suggesting an under-estimation on the part of the perturbation solution in the span of $u_{z}$ near the sidewalls, and for increasingly large $\varTheta$ , the emergence of $u_{z}$ contributions in the outer regions that do not appear in the $O(\alpha )$ perturbation solution. Interestingly, the computational results do confirm that, as predicted by the perturbation solution, $u_{z}$ contributions confined to the regions near the sidewall span the same fraction of the channel angle regardless of the choice of $\varTheta$ .

4.2. Validity of analytical solution

The numerical results reported in figures 9 and 10 corroborate the predictions of the perturbation solution of (3.29a ) to a large degree. The evident differences between the numerical and analytical results are both reasonably within the advertised $O(\alpha )$ accuracy, and helpful for illustrating one additional point of physical significance in the model formulation. Specifically, the non-zero $u_{z}$ contributions in the outer regions evoke terms of $O(\alpha ^{2})$ and smaller that appear in (2.8b )–(2.8c ), which were discarded in formulating the $O(1)$ and $O(\alpha )$ governing equations for the outer regions. These terms are sufficiently large to influence the pressure field to $O(\alpha ^{2})$ . Were the perturbation analysis to be extended to encompass further corrections, these corrections would account for said terms as well.

The validity of the perturbation analysis, like any analysis, extends only as far as the simplifying assumptions hold true. Considering that one of the key assumptions for this model is the restriction of $\alpha \ll 1$ , it is worth emphasising that the cases of figures 9 and 10 comprise relatively large $\alpha$ . Here, the numerical results reported in figures 9 and 10 are helpful in illustrating specific examples of what ‘relatively large’ may mean. In the context of the perturbation analysis developed in this work, even $\alpha \sim 10^{-1}$ may be large enough to produce obvious inaccuracies, especially for abrupt sidewall angles for the channel. But this does not preclude the perturbation model from practical application among its intended range of problems: low- ${\textit{Re}}$ flow in planar microchannels, with channel sidewalls converging at an arbitrary finite angle, and with small channel depth relative to the width between sidewalls.