1. Introduction
The flow of a gas in a microchannel usually falls into the slip or the transition regime. The gaseous flow in a microchannel is dependent on non-dimensional numbers: the Knudsen (
$Kn$), Mach (
$Ma$) and Reynolds (
$Re$) numbers, out of which only two are independent while the third number can be calculated from 
$Kn = \sqrt {({{\rm \pi} \gamma }/{2})}({Ma}/{Re})$, where 
$\gamma$ is described as the ratio of specific heat of the gas. The Knudsen number becomes relatively large at the microscale because of the small characteristic length scale. The Navier–Stokes equations are no longer valid for flow in microchannels (Singh, Dongari & Agrawal Reference Singh, Dongari and Agrawal2014), micro-devices (Gad-el Hak Reference Gad-el Hak1999) or in hypersonic flows of space vehicles in low Earth orbits (Ivanov & Gimelshein Reference Ivanov and Gimelshein1998). The Knudsen number for such flows lies in the continuum-transition regime (Agarwal, Yun & Balakrishnan Reference Agarwal, Yun and Balakrishnan2001; Agrawal, Kushwaha & Jadhav Reference Agrawal, Kushwaha and Jadhav2020).
To model such flows, extended hydrodynamics equations have been proposed over the last fifty years. One of these models is the Burnett equation derived from the Chapman–Enskog expansion of the Boltzmann equation. Burnett (Reference Burnett1935) derived the original Burnett equation. However, there are a few limitations (Garcia-Colin, Velasco & Uribe Reference Garcia-Colin, Velasco and Uribe2008) of the Burnett equations, and Bobylev (Reference Bobylev1982) reported one of them – known as Bobylev's instability. Due to a very cumbersome form of the original Burnett equations and the limitations with their stability, several variations to the original form of the Burnett equations have since been proposed in the literature. For example, Chapman & Cowling (Reference Chapman and Cowling1970) presented the conventional Burnett equations; Zhong, MacCormack & Chapman (Reference Zhong, MacCormack and Chapman1993) presented the augmented Burnett equations for application to hypersonic flows, while Balakrishnan and Agarwal came up with the Bhatnagar–Gross–Krook (BGK)–Burnett equations (Balakrishnan Reference Balakrishnan1999). Recently, Singh, Jadhav & Agrawal (Reference Singh, Jadhav and Agrawal2017) derived the OBurnett equations, a stable set of Burnett equations for rarefied flows using the Onsager-reciprocity principle. Despite the complex nature of the Burnett equations, attempts have been made to solve the various forms of these equations.
The improvements to the original Burnett equations over the past few decades have made possible a comparison of the available forms of the Burnett equations. Yun, Agarwal & Balakrishnan (Reference Yun, Agarwal and Balakrishnan1998) compared the BGK–Burnett equations and the augmented Burnett equations with the Navier–Stokes solutions. They reported that, for very high Knudsen number (
$Kn=0.1$), the difference between the augmented Burnett equations is significant, whereas, the solutions of the BGK–Burnett equations are close to the Navier–Stokes equation. Fiscko & Chapman (Reference Fiscko and Chapman1989) compared the direct simulation Monte Carlo (DSMC) results to the solution of the Navier–Stokes, Burnett and super-Burnett equations for a one-dimensional shock structure. The solutions were presented for a hard-sphere gas, argon and a Maxwellian gas from 
$Ma = 1.3$ to 50. They found that solution of the Burnett equations is in much closer agreement to both the Monte Carlo and the experimental results as compared to that from the Navier–Stokes equations. However, they found that the conventional Burnett equations suffer from instability for a very fine grid. Zhong et al. (Reference Zhong, MacCormack and Chapman1993) corrected this instability by adding higher-order terms from the super-Burnett equations to the Burnett equations. Recently, Jadhav & Agrawal (Reference Jadhav and Agrawal2020) proposed Grad's second problem, a new problem to test various forms of the Burnett equations and other higher-order equations. The solution of the conventional Burnett equations was found to compare well with the available DSMC data. The analytical solutions of the Burnett equations in cylindrical coordinates were given by Singh & Agrawal (Reference Singh and Agrawal2014). Agrawal et al. (Reference Agrawal, Kushwaha and Jadhav2020) discussed various forms of the Burnett equations, the challenges involved in solving the cumbersome Burnett equations and the available solutions obtained using analytical and numerical approaches in the past decades to the various forms of the Burnett equations.
The aim of this paper is to derive an analytical solution of the conventional Burnett equations for the plane Poiseuille flow problem. The study of gaseous flow in microchannels has relevance to the design of microdevices, which find numerous applications. Therefore, an ample amount of work has been done in the past to describe the flow characteristics in a microchannel. Aoki, Takata & Nakanishi (Reference Aoki, Takata and Nakanishi2002) solved the BGK–Burnett equations for flow between parallel plates driven by an external body force. They presented an asymptotic analysis as well as a numerical study. Jadhav, Singh & Agrawal (Reference Jadhav, Singh and Agrawal2017) solved the same problem using the OBurnett equations; a good agreement of these new equations against benchmark data was noted. Cercignani & Daneri (Reference Cercignani and Daneri1963) solved the Poiseuille flow of a rarefied gas numerically and showed that the volume flow rate versus pressure obtained using the BGK model fits well with the experimental results. Xu & Li (Reference Xu and Li2004) and Bao & Lin (Reference Bao and Lin2008) investigated the transition flows in microchannels by solving the Burnett equations numerically. The DSMC technique was employed to investigate the flows in microelectromechanical devices to obtain insights which could facilitate the design and optimization of the mass flow rate (Piekos & Breuer Reference Piekos and Breuer1996; Cheng & Liao Reference Cheng and Liao2000). Singh et al. (Reference Singh, Dongari and Agrawal2014) provided an analytical solution using the augmented Burnett equations to investigate low-speed isothermal microscale gas flows. They also presented the results for normalized mass flow rate, friction factor and axial velocity profile. A perturbation based solution for gaseous flow in a long microchannel using a two-dimensional conventional Burnett equation was presented by Rath, Singh & Agrawal (Reference Rath, Singh and Agrawal2018). The study was an analytical investigation, which provided the solutions for the pressure distribution and velocity distribution inside a microchannel. These solutions were compared to the DSMC solutions obtained by Xu & Li (Reference Xu and Li2004). There have been experimental investigations carried out to find the pressure distribution in gaseous slip flow in a long microchannel (Pong et al. Reference Pong, Ho, Liu and Tai1994), to find the value of the tangential momentum accommodation coefficient (Ewart et al. Reference Ewart, Perrier, Graur and Méolans2007; Hemadri, Bhandarkar & Agrawal Reference Hemadri, Bhandarkar and Agrawal2018) and to show the existence of choking in a microchannel (Garg & Agrawal Reference Garg and Agrawal2019), among other studies.
In the present work, we solve a two-dimensional conventional Burnett equation which involves solving the partial differential equation analytically. The primary outcome is the derivation of a closed-form analytical expression for the pressure distribution, obtained by solving the cross-stream momentum equation inside a long microchannel in terms of the characteristics of the gas and the microchannel through which the flow is happening. The streamwise momentum equation has also been solved analytically for the velocity distribution in a microchannel. The solution was compared with numerical data and perturbation based analysis results obtained by Rath et al. (Reference Rath, Singh and Agrawal2018) and found to be in close agreement. The novelty in the work is in obtaining a closed-form analytical expression for pressure that can be used in future analysis.
2. Governing equations
The governing equations for flow in a microchannel are given by the standard mass and momentum equations
\begin{align} &\begin{bmatrix} 1 & 0 & 0 \\ 0 & \rho u & 0 \\ 0 & 0 & \rho u\end{bmatrix} \begin{bmatrix} \dfrac{\partial(\rho u)}{\partial x} \\ \dfrac{\partial u}{\partial x} \\ \dfrac{\partial v}{\partial x} \end{bmatrix} + \begin{bmatrix} 1 & 0 & 0 \\ 0 & \rho v & 0 \\ 0 & 0 & \rho v \end{bmatrix} \begin{bmatrix} \dfrac{\partial(\rho v)}{\partial y} \\ \dfrac{\partial u}{\partial y} \\ \dfrac{\partial v}{\partial y} \end{bmatrix} + \begin{bmatrix} 0 \\ \dfrac{\partial(p+\sigma_{xx})}{\partial x} \\ \dfrac{\partial(\sigma_{xy})}{\partial x} \end{bmatrix}\nonumber\\ &\quad + \begin{bmatrix} 0 \\ \dfrac{\partial(\sigma_{xy})}{\partial y}\\ \dfrac{\partial(p+\sigma_{yy})}{\partial y}\end{bmatrix}=0 \end{align}
where 
$u$ and 
$v$ represent the streamwise and cross-streamwise velocities, respectively, 
$x$ and 
$y$ are the directions defined in figure 1, 
$\sigma _{ij}$ is the viscous stress tensor, 
$p$ is the pressure and 
$T$ is the absolute temperature. The pressure and temperature are related by the ideal gas law, i.e. 
$p=\rho R T$, where 
$R$ is the ideal gas constant, and 
$\rho$ is the density of ideal gas.
Microchannel schematic showing the 
$x$–
$y$ axis and streamwise velocity for steady two-dimensional and isothermal flow.
2.1. Stress and heat flux terms
The stress and heat flux terms are divided into Navier–Stokes- and Burnett-order terms. The stress and heat flux terms are therefore given as
\begin{gather} \sigma_{ij}=\sigma_{ij}^{N-S}+\sigma_{ij}^{B-E}, \end{gather}
\begin{gather}q_{i}=q_{i}^{N-S}+q_{i}^{B-E}. \end{gather}The combined Burnett and Navier–Stokes stress and heat flux terms are substituted in (2.1).
The Navier–Stokes terms on the right-hand side of (2.2) and (2.3) are given as
\begin{gather} \sigma_{ij}^{N-S}=2\mu\left[\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{1}{3}\delta_{ij}\frac{\partial u_k}{\partial x_k}\right], \end{gather}
\begin{gather}q_i^{N-S}={-}\kappa\frac{\partial T}{\partial x_i}. \end{gather}The Burnett terms in the right-hand side of (2.2) and (2.3) in the conventional Burnett equations for hard-sphere gas molecules are given below. However, these expressions are general for all interatomic potentials and the effect of the type of molecule is expressed solely through the Burnett coefficients
\begin{align} \sigma_{ij}^{B-E}&= k_1\dfrac{\mu^2}{p} \dfrac{\partial u_k}{\partial x_k} \frac{\partial u_{\langle i}}{\partial x_{j\rangle }} + k_2 \dfrac{\mu^2}{p}\left[ - \left(\dfrac{\partial}{\partial x_{\langle i}} \left[\dfrac{1}{\rho} \dfrac{\partial p}{\partial x_{j\rangle }}\right] + \dfrac{\partial u_k}{\partial x_{\langle i}} \dfrac{\partial u_{j\rangle }}{\partial x_k} + 2 \dfrac{\partial u_k}{\partial x_{\langle i}} S_{j\rangle k} \right) \right] \nonumber\\ &\quad + k_3 \dfrac{\mu^2}{\rho T} \dfrac{\partial^2 T}{\partial x_{\langle i} \partial x_{j\rangle }} + k_4 \dfrac{\mu^2}{\rho^2 R T^2} \dfrac{\partial T}{\partial x_{\langle i}} \dfrac{\partial p}{\partial x_{j\rangle }} + k_5 \dfrac{\mu^2}{\rho T^2} \dfrac{\partial T}{\partial x_{\langle i}} \dfrac{\partial T }{\partial x_{j\rangle }} + k_6 \frac{\mu^2}{p}\frac{\partial u_{k}}{\partial x_{\langle i }}\frac{\partial u_{ j\rangle }}{\partial x_{k }}, \end{align}
\begin{align} q_{i}^{B-E}&=\theta_1\frac{\mu ^2}{\rho T}\frac{\partial u_k}{\partial x_k}\frac{\partial T}{\partial x_i} + \theta_2\frac{\mu ^2}{\rho T}\left[-\frac{2}{3}\frac{\partial}{\partial x_i}\left(T\frac{\partial u_k}{\partial x_k}\right)\right] +\theta_2\frac{\mu ^2}{\rho T}\left({-}2\frac{\partial u_k}{\partial x_i}\frac{\partial T}{\partial x_k}\right)\nonumber\\ &\quad +\theta_3\frac{\mu ^2}{\rho p}\left(\frac{\partial u_{\langle i}}{\partial x_{k\rangle }}\frac{\partial p}{\partial x_k}\right)+\theta_4\frac{\mu ^2}{\rho}\left(\frac{\partial}{\partial x_k}\frac{\partial u_{\langle i}}{\partial x_{k\rangle }}\right) + 3\theta_5\frac{\mu ^2}{\rho T}\left(\frac{\partial u_{\langle i}}{\partial x_{k\rangle }}\frac{\partial T}{\partial x_k}\right), \end{align}
where 
$S_{ij}$ is the symmetric and trace free part of the velocity gradient tensor and is defined as
\begin{equation} S_{ij}=\frac{\partial u\langle i}{\partial x j\rangle}=\frac{1}{2}\frac{\partial u_i}{\partial x_j}+\frac{1}{2}\frac{\partial u_j}{\partial x_i}-\frac{1}{3}\frac{\partial u_k}{\partial x_k}\delta_{ij}, \end{equation}
where 
$\delta _{ij}$ is the Kronecker delta. Although the original Burnett equations were given by Burnett (Reference Burnett1935), these were modified to derive the conventional Burnett equations. Specifically, Chapman & Cowling (Reference Chapman and Cowling1970) replaced the material derivative terms using the Euler equations and also obtained the above mentioned Burnett stress and heat flux terms.
3. Solution
The solution is obtained under the following assumptions: the flow is assumed to be isothermal, two-dimensional, with no external body force applied to it. The gaseous flow is considered to have a low Reynolds number and a low Mach number, while the Knudsen number is of order unity. Arkilic, Schmidt & Breuer (Reference Arkilic, Schmidt and Breuer1997) also employed a similar set of conditions to perform a two-dimensional analysis of gaseous flow in a long microchannel using the Navier–Stokes equation. Similarly, Rath et al. (Reference Rath, Singh and Agrawal2018) solved the gaseous flow in a long microchannel using two-dimensional Burnett equations under the same set of assumptions.
The governing equations are the equations of mass and momentum containing the terms of the Navier–Stokes and the Burnett order, and thus (2.4) and (2.6) have to be substituted in (2.1). The full form of the governing equations obtained after substituting the stress and heat flux terms are rather cumbersome; these equations are given in Rath et al. (Reference Rath, Singh and Agrawal2018) but not reproduced here for brevity. The further steps involved in the derivation are as follows: first, the variables involved in the governing equation are non-dimensionalized as
\begin{gather} x^{*}=\frac{x}{L},\quad y^{*}=\frac{y}{H},\quad u^{*}=\frac{u}{U_{out}}, \quad p^{*}=\frac{p}{p_{out}},\quad {\rho}^{*}=\frac{\rho}{\rho_{out}},\quad v^*=\frac{v}{U_{out}}, \end{gather}
where 
$L$ and 
$H$ are respectively the length and height of the microchannel, and ‘out’ denotes the reference location (taken at the exit of the microchannel). Therefore, 
$p_{out}$ is the reference pressure used to normalize the pressure; similarly, 
$u_{out}$ is used to normalize the velocity. Next, since this problem concerns gaseous flow in a long microchannel, a smallness parameter 
$\epsilon ={H}/{L}$ has been employed. This helps us to identify terms with a small magnitude; such terms can then be discarded.
Hence, after implementing the aforementioned assumptions, and upon neglecting the higher-order terms of the small parameter ‘
$\epsilon$’, the following equations are obtained:
\begin{gather} \epsilon\frac{\partial (\rho u)}{\partial x} + \frac{\partial (\rho v)}{\partial y} = 0, \end{gather}
\begin{gather} -\frac{2k_2}{3 p^2}\frac{\epsilon}{Re}\frac{\partial^3 p}{\partial x \partial^2 y} - \frac{2k_2}{p^4}\frac{\epsilon}{Re}\frac{\partial p}{\partial x}\left(\frac{\partial p}{\partial y} \right)^2 + \frac{k_2}{3 p^3}\frac{\epsilon}{Re}\frac{\partial^2 p}{\partial x\partial y} \frac{\partial ^2p}{\partial y^2}\nonumber\\ \hspace{-4.7pc}+\frac{7 k_2}{3p^3}\frac{\epsilon}{Re}\frac{\partial p}{\partial y}\frac{\partial ^2p}{\partial x \partial y} + \frac{Re}{\gamma}\frac{\epsilon}{Ma^2}\frac{\partial p}{\partial x} = \frac{\partial^2 u}{\partial y^2}, \end{gather}
\begin{gather} \frac{8k_2}{3Re}p\frac{\partial p}{\partial y}\frac{\partial ^2p}{\partial y^2} - \frac{2k_2}{3Re}p^2\frac{\partial^3 p}{\partial y^3}-\frac{2k_2}{Re}\left(\frac{\partial p}{\partial y}\right)^3 +\frac{Re}{Ma^2} \frac{p^4}{\gamma}\frac{\partial p}{\partial y} = 0. \end{gather}
Note that all the terms in the equation have already been non-dimensionalized; however, for simplicity, the variables involved are written after dropping ‘
$*$’. The term 
$\gamma$ in (3.3) and (3.4) represents the ratio of specific heat capacities of the gas, as already mentioned above. The momentum equations involve two variables of interest, i.e. velocity ‘
$u$’ and pressure ‘
$p$’. We employ an analytical approach to find these variables.
3.1. Solution for pressure field
For obtaining the final solution of the pressure field from (3.4), multiple arrangements have been performed during the entire solution process, as explained below. In the first step, we re-arrange (3.4) by dividing the equation by 
$p^4$ to obtain
\begin{equation} \frac{8\hat{k}_2}{3p^3}\frac{\partial p}{\partial y}\frac{\partial ^2p}{\partial y^2} - \frac{2\hat{k}_2}{3p^2}\frac{\partial^3 p}{\partial y^3}-\frac{2\hat{k}_2}{p^4}\left(\frac{\partial p}{\partial y}\right)^3 + \frac{1}{\gamma}\frac{\partial p}{\partial y} = 0, \end{equation}where,
\begin{equation} \hat{k}_2=\frac{k_2 Ma^2}{Re^2}. \end{equation}Note that (3.5) is a non-standard nonlinear differential equation. Such equations need special mathematical techniques, involving combining of some terms in the equation to make it simpler to deal with it.
Hence, in order to solve (3.5) we split the left-hand side of the equation into two parts: one part consists of all the terms independent of 
${1}/{\gamma }$, while the second part contains terms with 
${1}/{\gamma }$. Accordingly,
\begin{equation} \left.\begin{gathered} (a)=\frac{8\hat{k}_2}{3p^3}\frac{\partial p}{\partial y}\frac{\partial ^2p}{\partial y^2} - \frac{2\hat{k}_2}{3p^2}\frac{\partial^3 p}{\partial y^3}-\frac{2\hat{k}_2}{p^4}\left(\frac{\partial p}{\partial y}\right)^3\\ (b)=\frac{\partial p}{\partial y}. \end{gathered}\right\} \end{equation}
It is evident that part (a) of (3.5) has a term consisting of a third-order derivative, a second-order derivative term of 
$p$ and the first-order term of third degree, respectively. The presence of a third-order partial differential equation with a constant coefficient 
$\hat {k}_2$ multiplied by all three terms of part (a) gives us a hint and motivation to exploit the equation to reduce the terms in (3.5). The expression is independent of 
$x$, and the arrangement of the derivatives of 
$p$ shows symmetry and encourages us to use the differential formula with two terms to retain the third term. Finally, by using a trial and error method, we came up with an expression denoted by 
$F$ and given as
\begin{equation} F = \frac{2\hat{k}_2}{3}\left[\frac{1}{p^3}\left(\frac{\partial p}{\partial y}\right)^2-\frac{1}{p^2}\frac{\partial^2 p}{\partial y^2}\right]. \end{equation}
Taking the 
$y$-derivative of (3.8) gives part (a) of (3.5) as shown below,
\begin{equation} \frac{\partial}{\partial y}\left(\frac{2\hat{k}_2}{3}\left(\frac{1}{p^3}\left(\frac{\partial p}{\partial y}\right)^2-\frac{1}{p^2}\frac{\partial^2 p}{\partial y^2}\right)\right) =\frac{8\hat{k}_2}{3p^3}\frac{\partial p}{\partial y}\frac{\partial ^2p}{\partial y^2} - \frac{2\hat{k}_2}{3p^2}\frac{\partial^3 p}{\partial y^3}-\frac{2\hat{k}_2}{p^4}\left(\frac{\partial p}{\partial y}\right)^3. \end{equation}
Adding the second part, i.e. 
$({{1}/{\gamma })({\partial p}/{\partial y}})$ to both sides of (3.9) gives,
\begin{align} &\frac{\partial}{\partial y}\left[\frac{2\hat{k}_2}{3}\left(\frac{1}{p^3}\left(\frac{\partial p}{\partial y}\right)^2-\frac{1}{p^2}\frac{\partial^2 p}{\partial y^2}\right)\right] + \frac{1}{\gamma}\frac{\partial p}{\partial y} = \frac{8\hat{k}_2}{3p^3}\frac{\partial p}{\partial y}\frac{\partial ^2p}{\partial y^2} - \frac{2\hat{k}_2}{3p^2}\frac{\partial^3 p}{\partial y^3}\nonumber\\ &\quad -\frac{2\hat{k}_2}{p^4}\left(\frac{\partial p}{\partial y}\right)^3+ \frac{1}{\gamma}\frac{\partial p}{\partial y}. \end{align}Hence, (3.5) can also be written as
\begin{equation} \frac{\partial F}{\partial y} + \frac{1}{\gamma}\frac{\partial p}{\partial y} = 0, \end{equation}
which is in a readily integrable form. Integrating (3.11) with respect to 
$y$ gives,
\begin{equation} \frac{2\hat{k}_2}{3}\left[\frac{1}{p^3}\left(\frac{\partial p}{\partial y}\right)^2-\frac{1}{p^2}\frac{\partial^2 p}{\partial y^2}\right] =c(x)- \frac{p}{\gamma}, \end{equation}
where 
$c$ is the integration constant and is a function of 
$x$. The resulting equation (3.12) is a second-order differential equation. Solving (3.12) for 
$p(x,y)$, generates two more integration constants which are a function of 
$x$. The analytical expression for 
$p(x,y)$ is given by
\begin{gather} p(x,y) = \frac{8c_1(x)^2c_2(x)\gamma^2\hat{k}_2^2e^{c_1(x)y}}{36c(x)^2\gamma^4\hat{k}_2^2-24c_1(x)^2\gamma^3\hat{k}_2^3 + 12c(x)c_2(x)\gamma^2\hat{k}_2e^{c_1(x)y} + c_2(x)^2e^{2c_1(x)y}}. \end{gather}Writing (3.13) in a slightly simpler way gives the following closed-form analytical expression:
\begin{equation} p(x,y)=\frac{8c_1(x)^2c_2(x)\gamma^2\hat{k}_2^2e^{c_1(x)y}}{\left(6c(x)\gamma^2\hat{k}_2+c_2(x)e^{c_1(x)y}\right)^2-24c_1(x)^2\gamma^3\hat{k}_2^3}. \end{equation}
The above expression is an exact analytical solution for the pressure field for fluid flow in a long microchannel, to which different boundary conditions can be applied for the purpose of evaluating the constants of integration involved in the above expression. The novelty in this solution is coming up with the form of 
$F$, which allowed us to obtain the above analytical solution. In a subsequent subsection, we will show the applicability of the solution obtained above.
3.2. Solution for velocity field
The velocity equation has been obtained from the non-dimensional 
$x$-momentum equation (3.3), which is derived from (2.1) after the non-dimensionalization of the field variables. Once again, upon simplifying and substituting (3.3) gives
\begin{gather} {-}\frac{2\hat{k}_2}{3 p^2}\frac{\partial^3 p}{\partial x \partial y^2} - \frac{2\hat{k}_2}{p^4}\frac{\partial p}{\partial x}\left(\frac{\partial p}{\partial y}\right)^2 + \frac{\hat{k}_2}{3 p^3}\frac{\partial p}{\partial x}\frac{\partial ^2p}{\partial y^2} + \frac{7 \hat{k}_2}{3p^3}\frac{\partial p}{\partial y}\frac{\partial ^2p}{\partial x \partial y} + \frac{1}{\gamma}\frac{\partial p}{\partial x} = \frac{\partial^2 u}{\partial y^2}\frac{Ma^2}{\epsilon Re}. \end{gather}The procedure to solve (3.15) takes into account the form of (3.14). Therefore, (3.14) is written in an alternate form
\begin{equation} p(x,y)= \frac{ f e^{c_1 y}}{-F_3 + (F_4 + c_2e^{c_1 y})^2}, \end{equation}where,
\begin{equation} \left.\begin{gathered} f= 8c_1(x)^2\gamma^2 \hat{k}_2^2c_2(x),\\ F_3 = 24c_1(x)^2\gamma^3\hat{k}_2^3,\\ F_4 = 6 c(x)\gamma^2\hat{k}_2. \end{gathered}\right\} \end{equation}Substituting (3.16) in (3.15) and then integrating the equation yields
\begin{equation} \frac{\textrm{d}u}{\textrm{d} y}= \frac{1}{(3f(X)^2)\gamma Re Ma^2}\int{I_{Term}}\,{\textrm{d} y} +3 F_5(x) f(x)^2 Re \gamma Ma^2, \end{equation}
where 
$F_5$ is the constant of integration and 
$I_{Term}$ is given by (3.19)
\begin{align} I_{Term}&= (([[-(24\gamma)\hat{k}_2 Ma^2c_1(x)^2({-}F_4(x)^{2/3} + F_3(x))c_2(x)^2 + (6Re^2)f(x)^2]F_4(x)'F_4(x)\nonumber\\ &\quad -(4c_1(x)^2)c_2(x)Ma^2\hat{k}_2\gamma(2F_4(x)^2 + F_3(x))({-}F_4(x)^2 + F_3(x))c_2(x)'\nonumber\\ &\quad + ((6\gamma)\hat{k}_2Ma^2c_1(x)^2(F_4(x)^2/3 + F_3(x))c_2(x)^2 - (3Re^2)f(x)^2)F_3(x)']f(x)\nonumber\\ &\quad - (8f(x))(\hat{k}_2c_1(x)\gamma Ma^2(c_1(x)y F_3(x) - (29F_4(x)^2)/8 + (13F_3(x))/8)c_2(x)^2 \nonumber\\ &\quad -(3y)Re^2f(x)^2/8)({-}F_4(x)^2 + F_3(x))c_1(x)' - 4(\gamma\hat{k}_2Ma^2c_1(x)^2({-}2F_4(x)^2 \nonumber\\ &\quad + F_3(x))c_2(x)^2 - (3Re^2)f(x)^2/4)({-}F_4(x)^2 + F_3(x))f(x)')\exp(c_1(x)y) \nonumber\\ &\quad + ((((-(20\gamma)Ma^2(-(3F_4(x)^2)/5 + F_3(x))\hat{k}_2c_1(x)^2c_2(x)^2 + (6Re^2)f(x)^2)F_4(x)' \nonumber\\ &\quad+ (4c_1(x)^2)c_2(x)^2F_3(x)'Ma^2\hat{k}_2\gamma F_4(x))c_2(x)+ (-(12\gamma)\hat{k}_2Ma^2c_1(x)^2({-}F_4(x)^2 \nonumber\\ &\quad+ F_3(x))c_2(x)^2 + (6Re^2)f(x)^2)c_2(x)'F_4(x))f(x) + ((12\gamma)\hat{k}_2Ma^2c_1(x)^2({-}F_4(x)^2 \nonumber\\ &\quad + F_3(x))c_2(x)^2 - (6Re^2)f(x)^2)c_2(x)f(x)'F_4(x) + (36c_1(x))c_2(x)^3Ma^2\hat{k}_2\gamma f(x)\nonumber\\ &\quad \times F_4(x)({-}F_4(x)^2 + F_3(x))c_1(x)')\exp((2c_1(x))y) + (((-(12\gamma)Ma^2(-(2F_4(x)^2)/3 \nonumber\\ &\quad +F_3(x))\hat{k}_2c_1(x)^2c_2(x)^2 +(6Re^2)f(x)^2)c_2(x)'+(2c_1(x)^2)c_2(x)^3Ma^2\hat{k}_2\nonumber\\ &\quad \times \gamma((4F_4(x)')F_4(x)+ F_3(x)'))f(x) + (-(8\hat{k}_2)c_1(x)\gamma Ma^2(c_1(x)y F_3(x) \nonumber\\ &\quad +(29F_4(x)^2)/8 - (13F_3(x))/8)c_2(x)^2 + (3y)Re^2f(x)^2)f(x)c_2(x)c_1(x)' \nonumber\\ &\quad +((4\gamma)\hat{k}_2Ma^2c_1(x)^2({-}2F_4(x)^2 + F_3(x))c_2(x)^2 -(3Re^2)f(x)^2)c_2(x)f(x)')\nonumber\\ &\quad \times c_2(x)\exp((3c_1(x))y) - (2\gamma)Ma^2c_2(x)^4((({-}c_2(x)F_4(x)' - c_2(x)'F_4(x))f(x) \nonumber\\ &\quad +c_2(x)f(x)'F_4(x))c_1(x)+(7c_2(x))c_1(x)'f(x)F_4(x))\hat{k}_2c_1(x)\exp((4c_1(x))y) \nonumber\\ &\quad -(2\gamma)Ma^2({-}F_4(x)^2 + F_3(x))^2((({-}c_2(x)F_4(x)' - c_2(x)'F_4(x))f(x) \nonumber\\ &\quad +c_2(x)f(x)'F_4(x))c_1(x) + (7c_2(x))c_1(x)'f(x)F_4(x))\hat{k}_2c_1(x) \nonumber\\ &\quad +(3c_1(x))c_1(x)'Ma^2\hat{k}_2\gamma f(x)({-}F_4(x)^2 + F_3(x))^3\exp({-}c_1(x)y) \nonumber\\ &\quad -(3\hat{k}_2)\exp((5c_1(x))y)c_1(x)'c_2(x)^6c_1(x)f(x)\gamma Ma^2)\epsilon1/(-(2c_2(x)^2)({-}3F_4(x)^2\nonumber\\ &\quad+ F_3(x))\exp((2c_1(x))y) + (4\exp((3c_1(x))y))c_2(x)^3F_4(x)+ \exp((4c_1(x))y)c_2(x)^4 \nonumber\\ &\quad- 4(\exp(c_1(x)y)F_4(x)c_2(x) + F_4(x)^2/4 - F_3(x)/4)({-}F_4(x)^2 + F_3(x))). \end{align}Equation (3.18) can be integrated again to obtain the analytical expression for velocity.
4. Results
In this section, we present our results and benchmark them with the available solutions in the literature.
4.1. Solution for pressure
In order to find the constants present in (3.14), the required boundary conditions are obtained as follows:
(i) Due to the symmetry of the flow, the normal pressure gradient at the centreline is zero,
(4.1)
\begin{equation} {\left.\left(\frac{\partial p}{\partial y}\right)\right|_{y=0}=0}. \end{equation}(ii) Pressure at the centreline
${(p|_{y=0}=p_c)}$ is obtained from DSMC data.(iii) The second-order curvature of pressure is also obtained from DSMC data, and given as
$a$,
(4.2)\begin{equation} {\left.\frac{\partial^2 p}{\partial y^2}\right|_{y=0}=a}. \end{equation}
Using the above boundary conditions the constants of integration are deduced as follows:
\begin{equation} \left.\begin{gathered} c(0) ={-}\frac{2a\gamma \hat{k}_2 - 3p_c^3}{3\gamma p_c^2},\\ c_1(0)=\sqrt{\frac{3p_c^3-4a\gamma \hat{k}_2}{2\gamma \hat{k}_2 p_c}},\\ c_2(0)={-}\frac{4\gamma^2 \hat{k}_2^2a}{p_c^2}. \end{gathered}\right\} \end{equation}Substituting the above constants for the expression of pressure, we get the closed-form solution for pressure in its simplest form.
Figure 2 shows the pressure distribution with the wall-normal position inside a long microchannel. The Knudsen number is taken as 0.1 and 
$p_{out}=6.05 \times 10^{-4}$ Pa, as used by Zheng, Garcia & Alder (Reference Zheng, Garcia and Alder2002). The two boundary conditions of second-order curvature (
${{\partial ^2 p}/{\partial y^2}}$) and centreline pressure (
$p_c$), are obtained from Xu (Reference Xu2003) for calculating the variation of pressure across the microchannel. The analytical pressure obtained from (3.14) matches with the DSMC solutions obtained by Zheng et al. (Reference Zheng, Garcia and Alder2002) and Xu (Reference Xu2003) for monoatomic gases having 
$\gamma =5/3$. Zheng et al. (Reference Zheng, Garcia and Alder2002) and Xu (Reference Xu2003) used the super-Burnett and BGK–Burnett equation solutions and compared against the numerical solution from DSMC. The result obtained from the present analytical approach has also been compared against the perturbation solution of Rath et al. (Reference Rath, Singh and Agrawal2018) in the same figure, and the solution overlaps with that of Rath et al. (Reference Rath, Singh and Agrawal2018).
Pressure distribution plotted against the wall-normal position in a long microchannel. The Knudsen number is 0.1.
The green line in the figure corresponds to the Navier–Stokes solution; as we know, the cross-stream Navier–Stokes equation simplifies to 
$\textrm {d}p/{\textrm {d} y}=0$. Therefore, according to the Navier–Stokes equation, the pressure should not vary across the cross-section. This limitation of the Navier–Stokes equation is also pointed out in the literature (Mansour, Baras & Garcia Reference Mansour, Baras and Garcia1997; Uribe & Garcia Reference Uribe and Garcia1999; Aoki et al. Reference Aoki, Takata and Nakanishi2002), that is, the Navier–Stokes equations cannot produce a non-uniform variation of pressure across the microchannel even after modifying the equation in different ways. However, pressure does vary with in the cross-stream direction at high Knudsen numbers as evident from the DSMC data in figure 2.
The DSMC results obtained from Zheng et al. (Reference Zheng, Garcia and Alder2002), Xu & Li (Reference Xu and Li2004) use a hard-sphere fluid with particle mass 
$m=1$ and diameter 
$d=1$, the thermal wall length is 1860 and the reference fluid speed 
$u_{out}=1$, speed of sound = 0.91 and 
$\gamma =5/3$ for a monoatomic gas. The reference pressure is 
$6.05 \times 10^{-4}$, pressure inlet and pressure outlet are respectively 
$(\frac {1}{2})p_{ref}$ and 
$(\frac {3}{2})p_{ref}$. The Knudsen number is 0.1 at the 
$x=0$.
Further comparison has also been carried out using the DSMC results of Mansour et al. (Reference Mansour, Baras and Garcia1997). Table 1 presents the non-dimentionalized parameters used for plotting the pressure and velocity profiles. The simulated fluid used by Mansour et al. (Reference Mansour, Baras and Garcia1997) is dilute gas with particle mass 
$m=1$ and diameter 
$d=1$. They have an acceleration-driven flow; therefore, the speed of sound is 0.9129, which gives 
$Ma=u/c = 1.548$, and the Knudsen number is 0.1. The thermal wall length is 1860. The reference fluid speed is 1. Figure 3 presents the comparison of the analytical pressure with the DSMC results of Mansour et al. (Reference Mansour, Baras and Garcia1997). The deviation between the analytical result and the DSMC result is plotted in the error percentage graph shown in figure 4. The error percentage is calculated as 
$(p_{ana}-p_{DSMC})/p_{DSMC}\times 100$. The maximum error recorded is less than 
$0.7\,\%$. The discrepancy at the wall is observed because we have a pressure-driven flow, whereas, Mansour et al. (Reference Mansour, Baras and Garcia1997) have an acceleration-driven flow.
Pressure distribution plotted against the wall-normal position in a long microchannel. The Knudsen number is 0.1 and 
$Ma=1.5480$.
Error in the percentage distribution between the analytical pressure and DSMC results obtained by Mansour et al. (Reference Mansour, Baras and Garcia1997) with respect to the normal-wall position. The Knudsen number is 0.1 and 
$Ma=1.5480$.
Table presenting the values of the parameters.
The value of 
$p_c$ is the non-dimensionalized value of pressure at the centreline; 
$a$ is the curvature obtained by curve fitting with the obtained DSMC data points; 
$p_{ref}$ is the reference pressure.
4.2. Solution for velocity
The expression for velocity is difficult to obtain because the 
$x$-momentum equation is lengthy and cumbersome compared to the 
$y$-momentum equation. However, we solved (3.15) by substituting the pressure expression obtained in (3.14) and integrating it once. The boundary conditions required to find the constants of integration are as follows:
(i) The symmetrical condition at the centreline
${({\partial u}/{\partial y}|_{y=0}=0)}$.(ii) Maxwell slip boundary condition at the wall
(4.4)\begin{equation} { u|_{wall}=\left(\frac{2-\sigma}{\sigma}\right)Kn \left.\frac{\partial u}{\partial y}\right|_{wall}}. \end{equation}
The velocity has been benchmarked using the above boundary conditions with the DSMC solution obtained previously and the perturbation solution proposed by Rath et al. (Reference Rath, Singh and Agrawal2018).
Figure 5 shows the velocity distribution along the wall-normal direction, which has been obtained after solving (3.18) numerically. The value of the tangential momentum coefficient, 
$\sigma$ employed in the calculations is unity. The analytical solution presented here corresponds to the flow with Mach number, 
$Ma=1.0989$, and Reynolds number, 
$Re=17.7$. The configurations of the microchannel used have been taken from Xu & Li (Reference Xu and Li2004). A comparison of results is made against the DSMC solution obtained by Xu (Reference Xu2003) and Zheng et al. (Reference Zheng, Garcia and Alder2002) and the perturbation solution presented by Rath et al. (Reference Rath, Singh and Agrawal2018). Figure 5 reveals that the result shown in the current paper is consistent with the previously reported results across the entire cross-section of the microchannel.
Velocity distribution plotted against the wall-normal position in a long microchannel, all non-dimensionalized. The Knudsen number is 0.1.
5. Discussion
Fluid flow in a long channel has always been an area of interest in both the academic and industrial fields. With the advancement in technology, manufacturing of long microchannels using micro-machining has become a reality, which led to extensive research on understanding the physics of gas flow in a microchannel in recent years. This extensive research has provided experimental results (Pong et al. Reference Pong, Ho, Liu and Tai1994; Arkilic et al. Reference Arkilic, Schmidt and Breuer1997), analytical solutions of the Navier–Stokes equations (Arkilic et al. Reference Arkilic, Schmidt and Breuer1997; Dongari, Agrawal & Agrawal Reference Dongari, Agrawal and Agrawal2007) and numerical solutions of the Burnett equation (Zheng et al. Reference Zheng, Garcia and Alder2002; Xu & Li Reference Xu and Li2004) for gas flow in a microchannel. The applicability of the analytical solutions of the Navier–Stokes equations, even after using different slip models, is limited to relatively small Knudsen number (Agrawal et al. Reference Agrawal, Kushwaha and Jadhav2020). This limitation of the Navier–Stokes equations led to the derivation of the Burnett equation, with a potentially larger envelop of Knudsen number (Agarwal et al. Reference Agarwal, Yun and Balakrishnan2001). However, the complexity of the Burnett equations and the highly nonlinear characteristics acted as an obstacle to progress in obtaining the analytical solution of these equations in the past. So although limited numerical solutions for flow in microchannels are available, analytical solutions are missing from the literature due to the cumbersome nature of the Burnett equations. Singh et al. (Reference Singh, Dongari and Agrawal2014) obtained a first analytical solution to the conventional Burnett equations for any configuration. They proposed and utilized a novel iterative approach, where physical insights about the flow were employed, along with the mathematical procedure to obtain the solution. Their novel approach was later employed to obtain solution for gas flow in a microtube by Singh & Agrawal (Reference Singh and Agrawal2014). Although Rath et al. (Reference Rath, Singh and Agrawal2018) derived an analytical expression for fluid flow in a long microchannel, their solution was based upon the perturbation method and is therefore still an approximation to the governing equations. In the current study, we have derived an exact analytical expression for both the pressure and velocity fields of a gas flowing in a long microchannel by using the conventional Burnett equations and approximations from Arkilic et al. (Reference Arkilic, Schmidt and Breuer1997). Hence, the result presented is a significant breakthrough in the area of analytical solution of the Burnett equations.
The solution presented here opens the possibility for further interesting investigations as we now have a closed-form analytical expression for the pressure distribution inside a long microchannel. For example, the variation of flow-induced stresses and entropy generation in microchannels of different sizes can be estimated using the expression provided here. A study which arrives at the optimum size of the microchannel by minimizing entropy production is available in Pallavi & Mahulikar (Reference Pallavi and Mahulikar2018). The other advantage is that it is significantly more convenient to benchmark data against the analytical solution presented here as compared to DSMC. Further, these results might change the perception about the Burnett equations in the community (Garcia-Colin et al. Reference Garcia-Colin, Velasco and Uribe2008), and are expected to inspire additional studies to find the analytical solution for other variants of the Burnett equations using similar or alternative approaches.
6. Conclusions
In this work, the conventional Burnett equations have been solved analytically to obtain the pressure and velocity fields for gas flow in a microchannel. The analytical expression presented in the current manuscript is valid for low Reynolds and Mach numbers. The obtained expressions for pressure and velocity have been validated against DSMC data and theoretical results available in the literature.
Our analysis has brought out the complete solution of the Burnett equations within an analytical framework, thereby putting the solution on a par with that obtained from the Navier–Stokes equations. The derived solution can be advantageously utilized, for example, in deriving solutions for heat transfer, entropy generation and undertaking other analyses in the future.
Funding
This work is funded by the Department of Atomic Energy, under its prestigious ‘DAE-SRC Outstanding Award’ scheme (No. 21/01/2015-BRNS/35038).
Declaration of interests
The authors report no conflict of interest.
Open access
