2.1. Statement of the problem

Consider the rarefied flow of argon gas through the straight and bent channels of the same height $h$ and axis length $L_{ch}=5h$ connecting the two reservoirs of equal size, see figure 2. The gas pressures in the upstream and downstream reservoirs are $p_{i}$ and $p_{e}$, respectively. The reservoir gas and wall temperatures are kept at $T_0=273\ \text {K}$. The resulting rarefied gas flow is characterised by the pressure ratio $\beta$ and the reference Knudsen number

(2.1a,b)\begin{equation} \textit{Kn}=\frac{\lambda}{L_{0}},\quad\lambda=\frac{\mu_{0}}{p_{0}}\sqrt{\frac{{\rm \pi} k_{B}T_{0}}{2m}}, \end{equation}

where $\mu _0$, $\lambda$, $m$ and $k_{B}$ are the gas viscosity at temperature $T_0$, the molecular mean free path, the molecular mass and the Boltzmann constant, respectively. The channel height and the inlet pressure are taken as the reference length $L_0=h$ and the reference pressure $p_0=p_i$, respectively. As a result, our reference Knudsen number is the inlet Knudsen number $\textit{Kn} = \textit{Kn}_i$. When the exit or mean Knudsen number ($\textit{Kn}_e$ or $\textit{Kn}_m$) is of interest, it can be calculated from the reference Knudsen number and pressure ratio $\beta$, i.e. $\textit{Kn}_e=\textit{Kn} / \beta$ or $\textit{Kn}_m = \textit{Kn}(1+\beta )/(2\beta )$, respectively. Alternatively, the gas rarefaction can be characterised by the rarefaction parameter $\delta$, which relates to $\textit{Kn}$ as $\delta ={\sqrt {{\rm \pi} }L_{0}}/{2\lambda }={\sqrt {{\rm \pi} }}/{2\textit{Kn}}$.

Figure 2.

Pressure-driven gas flow from the left-hand reservoir to the right-hand reservoir through (a) the straight channel and (b) the bent channel of the same height $h$ and axis length $L_{ch}=5h$. Note that only part of the reservoirs are illustrated here.

The effect of pressure ratio can be reflected by the Reynolds number $Re$ that is defined as

(2.2)\begin{equation} Re=\frac{\rho_0U_0L_0}{\mu_0}, \end{equation}

where $U_0$ is the characteristic flow speed, e.g. the average speed in the channel and $\rho _0=mp_0/(k_{B}T_0)$ is the reference density. The Mach number $Ma$ relates to the Knudsen number $\textit{Kn}$ and the Reynolds number $Re$ through the von Kármán relation

(2.3)\begin{equation} Ma=\sqrt{\frac{2}{{\rm \pi}\gamma}}\textit{Kn}Re, \end{equation}

where the specific heat ratio $\gamma$ is taken to be $5/3$ for monatomic gas.

2.2. The governing equation and boundary condition

The Boltzmann equation is the fundamental equation to describe the rarefied dynamics of monatomic gas. However, due to the complexity of its collision operator, the kinetic model (Shakhov Reference Shakhov1968) is often used to simulate rarefied gas flows. Without an external body force, the Shakhov model equation takes the following form:

(2.4)\begin{equation} \frac{\partial f}{\partial t} + \boldsymbol c \boldsymbol{\cdot} \boldsymbol{\nabla} f = \frac{ f^S-f}{\tau}. \end{equation}

Here $f = f(\boldsymbol x,\boldsymbol c, t)$ is the velocity distribution function of gas molecules with molecular velocity $\boldsymbol c = (c_x, c_y, c_z)$ at position $\boldsymbol x = (x, y, z)$ and time $t$, and the reference distribution function $f^S$ reads

(2.5)\begin{equation} f^{S}= \frac{n}{(2{\rm \pi} k_{B}T/m)^{3/2}}\exp\left( - \frac{m\xi^2}{2k_{B}T} \right) \left[1+\left(1-Pr\right)\frac{m\boldsymbol{\xi}\boldsymbol{\cdot}\boldsymbol q}{5n\left(k_{B}T\right)^{2}} \left(\frac{m\xi^{2}}{k_{B}T}-5\right)\right],\end{equation}

where $n$ and $T$ are the gas number density and temperature, respectively, $\boldsymbol \xi = \boldsymbol c - \boldsymbol u$ is the peculiar velocity, with $\boldsymbol u$ the macroscopic flow velocity and $\boldsymbol q=({m}/{2})\int \boldsymbol {\xi }\xi ^{2}f\,\text {d}\boldsymbol c$ is the heat flux. The Prandtl number $Pr$ is set to be $2/3$ for monatomic gas. Conservative flow variables are calculated as $\boldsymbol W \equiv (n, n\boldsymbol u, ne)^{\textrm {T}}=\int \boldsymbol ( 1, \boldsymbol c, c^2/2 )^{\textrm {T}} f \,\text {d}\boldsymbol c$ and the temperature is determined from the specific total energy $e = (u^2 + 3k_{B}T/m)/2$.

The relaxation time $\tau$ in the Shakhov model (2.4) is related to the dynamic viscosity $\mu$ and the local pressure $p$ as $\tau = \mu /p=\mu /(nk_{B}T)$. For gas molecules interacting through the inverse power law potential, the dynamic viscosity $\mu$ depends on the temperature $T$ as

(2.6)\begin{equation} \mu=\mu_{0}\left(\frac{T}{T_{0}}\right)^{\omega}, \end{equation}

where $\omega$ is the viscosity index. Argon gas with $m=6.63 \times 10^{-26}$ kg, $\mu _0=2.117 \times 10^{-5}$ Ns m$^{-2}$ and $\omega =0.81$ is considered in this paper.

To simulate the rarefied gas flow, the gas–surface interaction should be specified, which is modelled by the Maxwell diffuse–specular reflection, i.e.

(2.7)\begin{align} f\left(\boldsymbol{c} \mid \boldsymbol{c} \boldsymbol{\cdot} \boldsymbol{N} > 0 \right)&=\alpha_d n_s \left(\frac{m}{2{\rm \pi} k_{B}T_0}\right)^{3/2}\exp\left(-\frac{mc^2}{2k_{B}T_0}\right) \nonumber\\ &\quad + \left(1-\alpha_d\right)f\left[ \boldsymbol{c}-2\boldsymbol{N} (\boldsymbol{c} \boldsymbol{\cdot} \boldsymbol{N})\right],\quad (0\le\alpha_d\le 1), \end{align}

where $\boldsymbol {N},n_s,\alpha _d$ are the normal unit vector of the solid surface, the gas number density on the solid surface and the accommodation coefficient, respectively. In the diffuse–specular model, $\alpha _d$ portion of incident particles are reflected diffusively, whereas the remaining $(1-\alpha _d)$ portion of incident particles have a specular reflection. The diffuse boundary condition $\alpha _d=1$ at the solid surfaces is used in this study, except for the simulations in § 3.2 where diffuse–specular reflection is used. The gas number density on the solid surface is computed from the non-penetration condition, i.e. the number of gas molecules reflected from the wall is equal to those molecules approaching the same wall

(2.8)\begin{equation} n_s=-2\sqrt{\frac{\rm \pi}{2k_BT_0/m}}\int_{\boldsymbol{c} \boldsymbol{\cdot} \boldsymbol{N} < 0}\boldsymbol{c}\boldsymbol{\cdot} \boldsymbol{N} f \,{\mathrm{d}\boldsymbol{c}}. \end{equation}

At the free surfaces of reservoirs, molecules entering the computational domain follow the Maxwellian distribution with the local bulk velocity $\boldsymbol {u}$, pressure and temperature corresponding to each reservoir. The local velocity $\boldsymbol {u}$ at the free surfaces is extrapolated from that of the interior neighbour grid points.

2.3. The numerical methods

We employ the reference length $L_0$, temperature $T_0$, pressure $p_0$ and the most probable molecular speed $v_m=\sqrt {2k_{B}T_0/m}$ to normalise the following variables:

(2.9a–e)\begin{equation} \left(\tilde{X},\tilde{Y}\right)=\frac{\left(X,Y\right)}{L_0},\quad \left(\tilde{u},\tilde{v}\right)=\frac{\left(u,v\right)}{v_m},\quad \tilde{T}=\frac{T}{T_0},\quad \tilde{p}=\frac{p}{p_0},\quad \tilde{n}=\frac{n}{p_0/k_BT_0}, \end{equation}

where $u$ and $v$ are the velocity components of macroscopic velocity $\boldsymbol u$ along the coordinate axes $X$ and $Y$, respectively. Hereafter, the tildes on these dimensionless quantities are omitted for simplicity.

The DVM is one of the most commonly used deterministic approaches to solving the Boltzmann equation and its simplified models (Broadwell Reference Broadwell1964; Yang & Huang Reference Yang and Huang1995). It projects the continuous molecular velocity space $\boldsymbol c$ into a set of discrete velocities $\boldsymbol c^{(i)}$$(i=1,2,\ldots ,N_{c})$. As a result, the governing equation (2.4) is replaced by a system of $N_{c}$ independent equations. Here, we discretise this system in time by a time-implicit Godunov-type scheme (Yang & Huang Reference Yang and Huang1995; Titarev Reference Titarev2007)

(2.10)\begin{equation} \left.\begin{gathered} \left[\frac{1}{{\rm \Delta}{t^{(j)}}}+{\boldsymbol c^{(i)}}\boldsymbol{\cdot} \boldsymbol{\nabla}+\frac{1}{\tau^{(j)}}\right]{\rm \Delta} f^{(j)}=\text{RHS}^{(j)}, \\ \text{RHS}^{(j)}=\frac{1}{\tau^{(j)}}\left[f_{eq}^{(j)}-f^{(j)}\right]-{\boldsymbol c^{(i)}}\boldsymbol{\cdot} \boldsymbol{\nabla} f^{(j)}, \end{gathered}\right\} \end{equation}

where ${\rm \Delta} f^{(j)}=f^{(j+1)}-f^{(j)}$ needs to be determined at the time step ${\rm \Delta} {t^{(j)}}$ and $j$ is the time step index. Here, $\text {RHS}^{(j)}$ is the explicit part, and the spatial derivative is approximated by a third-order upwind scheme. This time-implicit scheme allows us to use a large time step to accelerate steady-state solution. Specifically, the classical Courant–Friedrichs–Lewy number of $10^6$ is chosen in this study. A few efficient implicit methods for solving the linearised Shakhov model have been developed by Titarev (Reference Titarev2012a,Reference Titarevb, Reference Titarev2013). The details of the DVM algorithm used in this work can be found in Ho et al. (Reference Ho, Li, Wu, Reese and Zhang2019). In the following simulations, the polar velocity grids of $N_c = N_{c_p} \times N_{\varphi }= 4\times 120$ and $4\times 40$ are chosen for $\textit{Kn}>0.1$ and $\textit{Kn} \le 0.1$, respectively. The number of discretised velocities in the angular coordinate $N_{\varphi }$ is uniformly spaced on $[0,2{\rm \pi} ]$, whereas the number of discretised velocities in the radial coordinate $N_{c_p}$ follows the half-range Gauss–Hermite abscissae. The molecular velocity is cut off at $2.3$ times the most probable molecular speed $v_m$ with $N_{c_p}=4$. In the cases of small pressure ratio, i.e. $\beta =0.75, 0.5$ and low Knudsen number $\textit{Kn} \le 0.03$, $N_{c_p}=8$ and the truncated molecular velocity of $3.7v_m$ are used due to relatively high Mach number. With these sets of the velocity grids, the numerical uncertainty is of the order of $0.1\,\%$ when $N_{c_p}$ is doubled (note: the truncated molecular velocity is also enlarged accordingly).

Different uniform spatial grids are used in DVM simulations based on the value of Knudsen number, i.e. $L10W4h160$, $L20W8h80$ and $L40W16h40$ are used for $\textit{Kn} \le 0.1$, $0.1 < \textit{Kn} \le 1$ and $1 < \textit{Kn}$, respectively. Here, for example, the spatial grid $L10W4h160$ denotes the reservoirs are of size of $L=10h$, $W=4h$ and the reference length (i.e. the channel height $h$) is resolved by $160$ uniform-cells, see figure 2.

Our simulations start from the global equilibrium state. The convergence criterion for the steady-state is checked every time step as follows:

(2.11)\begin{equation} E(t)=\frac{\sum \left|\boldsymbol{u}(t)-\boldsymbol{u}(t-{\rm \Delta} t)\right|}{\sum |\boldsymbol{u}(t)|}<10^{-6}. \end{equation}

The reduced mass flow rate $G$, normalised by $h^2(p_{i}-p_{e})/v_{m}L_{ch}$, reads

(2.12)\begin{equation} G=\frac{L_{ch}}{h}\frac{2p_0}{\left(p_{i}-p_{e}\right)} \int nu\,\text{d}Y=\frac{L_{ch}}{h}\frac{2p_0}{\left(p_{i}-p_{e}\right)}u_{avg}.\end{equation}

The average speed $u_{avg}=\int nu\,\text {d}Y$ is calculated at two vertical cross-sections inside the channel near the inlet and the exit to guarantee the mass conservation. The discrepancy of $u_{avg}$ near the inlet and the exit is less than $0.01\,\%$. To calculate the Reynolds number, the characteristic flow speed $U_0$ in (2.2) is determined by $U_0=u_{avg}v_m$.

In this study, DVM is used for solving the nonlinear Shakhov model, see (2.4) and (2.5) for all the values of pressure ratio $\beta$, except for $\beta =0.99999$. For $\beta =0.99999$, the DVM simulation converges very slowly and requires a very fine spatial grid at such low $\textit{Kn}$ (Valougeorgis & Naris Reference Valougeorgis and Naris2003; Wang et al. Reference Wang, Ho, Wu, Guo and Zhang2018). Therefore, the discontinuous Galerkin method (DGM) and general synthetic iterative scheme, which enables the fast convergence to the steady-state solution and retains asymptotic preserving natures of the Navier–Stokes equations, is employed to solve the linearised Shakhov equation for the case of pressure ratio $\beta =0.99999$ (Su et al. Reference Su, Wang, Zhang and Wu2019b, Reference Su, Zhu, Wang, Zhang and Wu2020b). The DVM and DGM solutions are compared with each other and with the other available data, see appendix A.

3.1. Influence of Reynolds number and Knudsen number

In our simulations, the Knudsen number $\textit{Kn}$ and the pressure ratio $\beta$ are chosen first, and the Reynolds number $Re$ is then obtained from the simulation results. The compressibility of the gas tends to suppress vortices generation which has been observed in supersonic flows, e.g. see Burggraf (Reference Burggraf1966). The compressibility effects may be insignificant as the Mach number $Ma$ calculated from the von Kármán relation (2.3) is up to $1.5\times 10^{-3}$, $1.5\times 10^{-2}$, $0.19$ and $0.25$ for all the considered cases of the pressure ratio $\beta =0.99999$, $0.99$, $0.75$ and $0.5$, respectively.

Figure 3 shows the streamlines and the pressure distribution inside the bent channel for $\textit{Kn}=0.01$, $0.02$ and $0.05$ with the pressure ratio $\beta =0.99$ and $0.5$. With the fixed pressure ratio $\beta$ and increasing $\textit{Kn}$, $Re$ becomes smaller so that the vortices at the bends shrink or even disappear. The strong influence of $Re$, i.e. the inertial effect from the mainstream, on the vortex formation is well known for continuum flows; however, with the presence of gas rarefaction effect in microflow, two unusual phenomena are observed in figure 3. First, the Reynolds number alone cannot characterise the flow separation in a microchannel. For example, from figures 3(c) and 3(f) we see that in the slip flow regime the vortices exist at $Re=0.68$ but not at a larger $Re$ of $4.72$. Second, as shown in figures 3(c) and 3(b), the vortices near the concave and convex corners occur at $Re$ as small as $0.68$ and $40.65$, respectively, while the corresponding critical Reynolds numbers (i.e. the smallest $Re$ where the vortex emerges) for air in the continuum flow regime, obtained from the Navier–Stokes equations with a no-slip boundary condition (Maharudrayya et al. Reference Maharudrayya, Jayanti and Deshpande2004), are $200$ and $100$; on the other hand, for water flow in a microchannel no vortex is observed experimentally when $Re<100$ (Xiong & Chung Reference Xiong and Chung2008). To analyse these contradictory observations, we need to look into the secondary flows at the bends.

Figure 3.

The pressure field and streamlines of rarefied argon gas flow inside the bent channel at different values of $\textit{Kn}$ and $Re$. Note: the rarefaction $\delta$ of $88.6, 44.3, 17.7$ corresponds to $\textit{Kn}=0.01, 0.02, 0.05$, respectively. (a) $\beta = 0.99, \textit{Kn}= 0.01, Re= 2.48$; (b) $\beta = 0.5, \textit{Kn}= 0.01, Re= 40.65$; (c) $\beta = 0.99, \textit{Kn}= 0.02, Re= 0.68$; (d) $\beta = 0.5, \textit{Kn}= 0.02, Re= 17.80$; (e) $\beta = 0.99, \textit{Kn}= 0.05, Re= 0.13$ and (f) $\beta = 0.50, \textit{Kn} = 0.05, Re= 4.72$.

To the best of authors’ knowledge, the concave vortex, which exists near the concave corners, is closely related to the Moffatt eddies originally found by the Stokes stream function for viscous fluid near deep sharp corners in Stokes flow (Moffatt Reference Moffatt1964). Contrary to the conventional vortex occurring at high $Re$, where the centrifugal force is balanced by the pressure force, the Moffatt viscous vortex is both driven and damped by the viscous force. By solving the stream function in the Stokes limit at the vicinity of a sharp corner, Moffatt found that a sequence of vortices generally appears if the angle of the sharp corner is smaller than the critical angle $146^\circ$. In the hydrodynamic regime, Moffatt eddies have been confirmed numerically by Biswas, Breuer & Durst (Reference Biswas, Breuer and Durst2004) in backward-facing step flows, where the primary vortex has been observed in a wide range of $Re$ from $10^{-4}$ to $10^2$. Interestingly, the vortex size is nearly constant when the Reynolds number increases up to unity, then becomes larger with a further increase of $Re$ (up to $100$). Also, the size ratio of the primary vortex to the secondary vortex obtained by their Navier–Stokes simulations for $Re=1$ is in good agreement with the Stokes approximation in the Moffatt analysis ($16.28$ for the right angle corner). This rapid reduction of the size of successive vortices leads to enormous difficulty in capturing a sequence of vortices experimentally. In gas microflow, the concave vortex was, surprisingly, observed in the numerical simulations with $Re$ of the order of unity (Agrawal et al. Reference Agrawal, Djenidi and Agrawal2009; White et al. Reference White, Borg, Scanlon and Reese2013). Without awareness of Moffatt's theory, the existence of a concave vortex at such low $Re$ is either attributed to gas rarefaction (Agrawal et al. Reference Agrawal, Djenidi and Agrawal2009) or not acknowledged (Wang & Li Reference Wang and Li2004; Sharipov & Graur Reference Sharipov and Graur2012; Liu et al. Reference Liu, Tang, Su, Wu and Zhang2018). We observe, for the first time, that the concave vortex in gas microflows does exist in the limit as the Reynolds number approaches zero, which agrees with Moffatt's theory. From our simulations, the smallest Reynolds number for a vortex to appear is found to be $0.32\times 10^{-3}$ for the case of pressure ratio $\beta =0.99999$ and $\textit{Kn}=0.03$, as shown in figure 4(a).

Figure 4.

Critical values of the Reynolds number for the emergence of (a) concave vortex (enlarged view of the pressure field and streamlines near the first concave corner of the upper wall) and (b) convex vortex. Note: the rarefaction $\delta$ of $88.6,29.5$ corresponds to $\textit{Kn}=0.01, 0.03$, respectively. (a) $Re= 0.32\times 10^{-3}, \textit{Kn}= 0.03, \beta = 0.99999$ and (b) $Re= 30.8, \textit{Kn}= 0.01, \beta = 0.75$.

The convex vortex, which appears in the downstream side of convex bent corners, is formed when the fluid streams cannot move closely along the sharp angle of the wall after the convex corner due to their relatively high momentum. The flow separation emerges from the convex corner point and extends beyond the vena contracta, i.e. the narrowest region of the stream after a sudden contraction of the duct (see the schematic diagram in figure 1). A similar effect is also observed in the slip flow through an orifice (Sharipov Reference Sharipov2004; Ho & Graur Reference Ho and Graur2014). This type of vortex was observed clearly from the experiments of water flow through a microchannel with bends at $Re \ge 300$ (Xiong & Chung Reference Xiong and Chung2008). To the best of the authors’ knowledge, this is the first study that detects this type of vortex for rarefied gas flow at a much smaller Reynolds number compared with continuum flow. From our simulations, the smallest Reynolds number for the emergence of a convex vortex is found to be $30.8$ in the case of pressure ratio $\beta =0.75$ and $\textit{Kn}=0.01$, in which a single convex vortex occurs near the second rectangular bend as shown in figure 4(b).

3.2. Influence of velocity slip

To understand the puzzling behaviour of concave and convex vortices for rarefied gas flows at low Reynolds number, we analyse the influence of velocity slip at the channel surface, which is one of the major rarefaction effects when $\textit{Kn}$ is small. Note that when the other parameters are fixed, the smaller the accommodation coefficient $\alpha _d$ is, the larger the slip velocity will be (Loyalka Reference Loyalka1968; Su et al. Reference Su, Wang, Liu and Wu2019a). Here, we change $\alpha _d$ and keep both $\textit{Kn}$ and $\beta$ fixed to allow significant variation of slip velocity $u_s$ at the wall while the mainstream velocity is only slightly affected. For example, in the following two cases the Reynolds numbers are changed by less than $5\,\%$ and $0.9\,\%$ when $\alpha _d$ is varied from $1$ to $0.6$.

The influence of slip velocity on the flow field near the concave corner is shown in figure 5, when $\textit{Kn}=0.01$ and $\beta =0.99$. It can be seen from figures 5(a) and 5(b) that, when $\alpha _d$ decreases (i.e. the slip velocity increases), the size of the vortex shrinks and its core centre moves towards the corner. Slip velocity profiles along the upper and lower walls, as shown in figures 5(c) and 5(d), demonstrate that no slip exists in the region near the concave corner (i.e. at $d/h\approx 2$ of the upper wall and $d/h\approx 3$ of the lower wall, where $d$ is the distance between the inlet of a microchannel and the considered point along the corresponding wall) and this non-slip region increases with $\alpha _d$ . This can be explained by the fact that, when $\alpha _d$ decreases, the induced larger slip velocity needs a longer distance to reduce to zero at the concave corner. As a consequence, the wall surface with no-slip flow becomes smaller and the Moffatt vortex is suppressed at a lower $\alpha _d$. A weak adverse pressure gradient is also found at the concave corner region, see figures 5(e) and 5(f).

Figure 5.

Vortex, slip velocity and adverse pressure gradient at the concave corner when $\textit{Kn}=0.01$, $\beta =0.99$. (a and b) Enlarged view of the pressure field and streamlines near the first concave corner (upper wall); (c and d) the slip velocity along the upper and lower walls; (e and f) the pressure along the upper and lower walls, and the channel axis. The distance $d$ is measured from the inlet of the microchannel along the upper wall, the lower wall or the channel axis, so that the concave corner is positioned at $d/h=2$ and $3$ for the upper and lower walls, respectively. The Reynolds numbers are $2.48$ and $2.60$ for $\alpha _d=1$ and $\alpha _d=0.6$, respectively. (a) $\alpha _d = 1$, the first concave corner; (b) $\alpha _d = 0.6$, the first concave corner; (c) $\alpha _d = 1$, slip velocity profile; (d) $\alpha _d = 0.6$, slip velocity profile; (e) $\alpha _d = 1$, pressure profile and (f) $\alpha _d = 0.6$, pressure profile.

The influence of slip velocity on the flow field at the convex corner is presented in figure 6, when $\textit{Kn}=0.01$ and $\beta =0.5$. In contrast to the concave vortex, the convex vortex is larger and its core centre moves away from the corner point when $\alpha _d$ decreases from 1 to 0.6, see figures 6(a) and 6(b). This can be explained by the fact that fluid streams, denoted as fluid streams ‘A’ in figure 1, departing from the upstream wall of the corner cannot follow a sharp rectangular turn and are separated from the downstream wall of the corner. These separations create a fluid region between the downstream wall of the corner and the separated streams ‘A’. As no fluid stream departing from the inlet can access this region, a recirculation occurs, resulting in the vortex in this region. Fluid streams ‘A’, through viscous force, drive the convex vortex. As the speed of fluid streams ‘A’ increases with reduction of $\alpha _d$ – see the slip velocity at the convex corner ($d/h=4$ for the upper wall and $d/h=1$ for the lower wall) shown in figures 6(c) and 6(d) – the convex vortex becomes larger. A larger vortex leads to a shortened vena contracta. When $\alpha _d$ decreases, the slip velocity of backward flow and separation length (from the corner point to the reattachment point) increases. The sudden jump and drop of slip velocity can also be seen at the convex corner, which is positively related to the abrupt drop of pressure profile along the wall before the convex corner shown in figures 6(e) and 6(f). An adverse pressure gradient along the wall after the convex corner is consistent with the pressure field near the convex corner as shown in figures 6(a) and 6(b). The adverse pressure gradient after the convex corner is significantly stronger than that at the concave corner ($d/h=2$ for the upper wall and $d/h=3$ for the lower wall), this can also be seen for the case shown in figure 5.

Figure 6.

Vortex, slip velocity and adverse pressure gradient at the convex corner when $\textit{Kn}=0.01$, $\beta =0.5$: (a and b) Enlarged view of the pressure field and streamlines near the second convex corner (upper wall); (c and d) the slip velocity along the upper and lower walls; (e and f) the pressure along the upper wall, the lower wall, and the channel axis. The convex corner is positioned at $d/h=4$ and $1$ for the upper and lower walls, respectively. The Reynolds numbers are $40.65$ and $41.01$ for $\alpha _d=1$ and $\alpha _d=0.6$, respectively. (a) $\alpha _d = 1$, the second convex corner; (b) $\alpha _d = 0.6$, the second convex corner; (c) $\alpha _d = 1$, slip velocity profile; (d) $\alpha _d = 0.6$, slip velocity profile; (e) $\alpha _d = 1$, pressure profile and (f) $\alpha _d = 0.6$, pressure profile.

3.3. Explanation of separation in rarefied gas flows

From the above numerical simulations, we know that the slip velocity suppresses flow separation at the concave corner and intensifies flow separation at the convex corner. This finding can explain the two unusual observations of rarefied flow separation presented in § 3.1. First, the momentum of bulk fluid is much larger in the case shown in figure 3(f) than that in figure 3(c), however, flow separation is not observed in figure 3(f). This is because the Moffatt vortex is suppressed by a significant increase of the slip velocity at the concave corner for $\textit{Kn}=0.05$, see figure 7. As the size of the Moffatt vortex is significantly reduced in rarefied flow, it is more difficult to capture in simulations and experiments. This explains why the Moffatt vortex was not detected in some numerical studies of rarefied gas flow in bent microchannels (Wang & Li Reference Wang and Li2004; Sharipov & Graur Reference Sharipov and Graur2012; Liu et al. Reference Liu, Tang, Su, Wu and Zhang2018), but was captured by others with much refined grids (Agrawal et al. Reference Agrawal, Djenidi and Agrawal2009; White et al. Reference White, Borg, Scanlon and Reese2013; Varade et al. Reference Varade, Agrawal, Prabhu and Pradeep2015). Second, the velocity slip in gas rarefied flows enhances the convex vortex leading to ‘early onset’ of flow separation in microsystems (with respect to the Reynolds number) compared with the flow in macrosystems, in which no-slip boundary condition is applied at the fluid–wall interface. The convex vortex in a rarefied gas can be found at Reynolds numbers as small as $30.8$, see figure 4(b).

Figure 7.

Slip velocity near the first concave corner (the upper wall, $d/h=2$) for the cases shown in figures 3(c) and 3(f). The distance $d$ is measured along the upper wall, starting from the inlet of the microchannel.

The attenuation of vortices near a bend with increasing $\textit{Kn}$ at fixed pressure ratio $\beta$, as shown in figure 3, can be interpreted by the reduction of $Re$, which is consistent with the findings in the continuum flow (Agrawal et al. Reference Agrawal, Djenidi and Agrawal2009; Varade et al. Reference Varade, Agrawal, Prabhu and Pradeep2015). This reason seems not sufficient in rarefied gas flow, especially for a concave vortex, which may occur even when $Re$ approaches zero, see figure 4(a). The underlying reason for the concave vortex and the additional reason for the convex vortex are the change of flow velocity due to the rarefaction effect characterised by the Knudsen number. From figure 8 we see that the slip velocity is also affected by the Knudsen number: at the concave corner, it increases with increasing $\textit{Kn}$, while at the convex corner it decreases with increasing $\textit{Kn}$. As a consequence, in agreement with the role of velocity slip on vortex formation as analysed in § 3.2, both types of vortices become smaller with increasing $\textit{Kn}$ at a fixed pressure ratio. Here, the concave and convex vortices start to disappear when the Knudsen number increases to $0.04$ and $0.01$, respectively.

Figure 8.

Influence of the Knudsen number on the slip velocity near (a) the first concave corner (the upper wall, $d/h=2$) for $\beta =0.99$, (b) the second convex corner (the upper wall, $d/h=4$) for $\beta =0.5$. The distance $d$ is measured along the upper wall, starting from the inlet of the microchannel.

The adverse pressure gradient along the channel walls, which can be conveniently detected by pressure sensors (Lee et al. Reference Lee, Wong and Zohar2001; Varade et al. Reference Varade, Agrawal, Prabhu and Pradeep2015), is usually used as an indicator for flow separation at bends. From our numerical data, adverse pressure gradients at both the concave and convex corners are reduced when the Knudsen number (Reynolds number) increases (decreases). The adverse pressure gradient along the concave wall exists in all the examined $Re$ and $\textit{Kn}$, while the adverse pressure gradient along the convex wall can only be found when $\textit{Kn}\le 0.1$, regardless of the pressure ratio $\beta$. However, the existence of adverse pressure gradient along the wall does not guarantee the flow separation as it may not lead to vortex generation. For example, although an adverse pressure gradient does occur at both the concave and convex corners for the case of pressure ratio $\beta =0.5$ and $\textit{Kn}=0.05$ as shown in figure 3(f), no flow separation is found for this case.