18.5.2 Navier-Stokes Equation for the Entrance Flow

18.5.2.1 Derivation

General Version of the Navier-Stokes Equation. The first thing we need is the modified Navier-Stokes equation. We neglect changes with respect to time, as the entrance effects are not time-dependent, but only dependent on z, which is why we can set $\frac{\partial \overrightarrow{v}}{\partial t}=0$. The flow is again driven by a pressure along the z-axis and there are no volume forces involved. Taken together, these assumptions will simplify Eq. 11.40 to

$\rho \left(\overrightarrow{v}·\overrightarrow{\nabla }\right)\overrightarrow{v}=-\frac{dp}{dz}+\eta \mathrm{\Delta }\overrightarrow{v}$

Eq. 18.65 is the Navier-Stokes equation that we need to solve. It is in general form, i.e., it can also be applied for systems which are not radially symmetric.

Neglecting the Convection Term.Eq. 18.65 is still a rather complex equation and we would have a hard time finding a general solution to it. We will therefore make another assumption that is in line with the assumption made for Poiseuille flow. We will consider only flows at a very low Reynolds number Re, i.e., we neglect momentum transport by convection. In most microfluidic systems, this is a valid assumption and it will make our problem significantly easier to solve. For low Reynolds numbers we can drop the entire left side of Eq. 18.65 and are left with

$\begin{array}{lll}0& =& -\frac{dp}{dz}+\eta \mathrm{\Delta }\overrightarrow{v}\\ \eta \mathrm{\Delta }\overrightarrow{v}& =& \frac{dp}{dz}\end{array}$

This equation should look familiar - it is the fundamental equation for pressure-driven Poiseuille flow (see Eq. 15.14) when neglecting effects of the volume forces.

Cylindrical Coordinates. We now convert Eq. 18.66 into cylindrical coordinates. For this we need the Laplace (see Eq. 7.99) operator in cylindrical coordinates that once applied to Eq. 18.66, results in the following equation

$\frac{\eta }{r}\frac{\partial }{\partial r}\left(r\frac{\partial v}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^{2}v}{\partial {\phi }^2}+\frac{{\partial }^{2}v}{\partial z^2}=\frac{dp}{dz}$

Please note that we have used the scalar Laplacian as we have converted the velocity vector $\overrightarrow{v}$ to the scalar v_z, because there is no flow in radial and azimuthal direction, thus v_φ = v_r = 0. This reduces the vector notation to a single equation of v_z. There are no changes of v_z in azimuthal direction, but now we need to consider changes along the z-axis. This is in contrast to the fully developed Hagen-Poiseuille flow which we derived in section 16.3 where we assumed a fully developed flow profile and thus neglected changes along the z-axis, i.e., $\frac{\partial v_x}{\partial z}=0$. See section 16.3.2 for details on this assumption. Taking together we simplify Eq. 18.67 to

$\begin{array}{lll}\frac{\eta }{r}\frac{\partial }{\partial r}\left(r\frac{\partial v_z}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2{v}_z}{\partial {\phi }^2}+\frac{{\partial }^2{v}_z}{\partial x^2}& =& \frac{dp}{dz}\\ \frac{\eta }{r}\frac{\partial }{\partial r}\left(r\frac{\partial v_z}{\partial r}\right)+\frac{{\partial }^2{v}_z}{\partial z^2}& =& \frac{dp}{dz}\end{array}$

Eq. 18.68 is as far as we can go with the simplifications. This is the equation that we will need to solve.

18.5.2.2 Transient and Steady-State Solution

We now apply our solution consisting of a steady state and transient solution as given by Eq. 18.64 to Eq. 18.68. In doing so we find

$\begin{array}{lll}\frac{\eta }{r}\frac{\partial }{\partial r}\left(r\frac{\partial v_{z,\mathrm{steady}‐\mathrm{state}+}{v}_{z,\mathrm{transient}}}{\partial r}\right)+\frac{{\partial }^2{v}_{z,\mathrm{steady}‐\mathrm{state}+}{v}_{z,\mathrm{transient}}}{\partial z^2}& =& \frac{dp}{dz}\\ \frac{\eta }{r}\frac{\partial }{\partial r}\left(r\frac{\partial v_{z,\mathrm{steady}‐\mathrm{state}}}{\partial r}+r\frac{\partial v_{z,\mathrm{transient}}}{\partial r}\right)+\frac{{\partial }^2{v}_{z,\mathrm{steady}‐\mathrm{state}}}{\partial z^2}+\frac{{\partial }^2{v}_{z,\mathrm{transient}}}{\partial z^2}& =& \frac{dp}{dz}\end{array}$

At this point we note that the steady-state solution is not a function of z therefore $\frac{{\partial }^2{v}_z,\mathrm{steady}‐\mathrm{state}}{\partial z^2}=0$. Furthermore we already know the steady-state solution to the Poiseuille flow that we derived in section 16.3 and that is given by Eq. 18.63. For Eq. 18.69 we need the derivative with respect to the radius that is given by

$\frac{\partial v_{z,\mathrm{steady}‐\mathrm{state}}}{\partial r}=-4v_{\mathrm{a}\mathrm{v}}\frac{r}{R^2}$

We therefore rewrite Eq. 18.69 to

$\begin{array}{lll}\frac{\eta }{r}\left(\frac{\partial \left(r{v}_{z,\mathrm{steady}‐\mathrm{state}}\right)}{\partial r}+\frac{\partial \left(r{v}_{z,\mathrm{transient}}\right)}{\partial r}\right)+\frac{{\partial }^2{v}_{z,\mathrm{transient}}}{\partial z^2}& =& \frac{dp}{dz}\\ \frac{\eta }{r}\left(\frac{\partial -4v_{\mathrm{a}\mathrm{v}}\frac{r^2}{R^2}}{\partial r}+\frac{\partial \left(r{v}_{z,\mathrm{transient}}\right)}{\partial r}\right)+\frac{{\partial }^2{v}_{z,\mathrm{transient}}}{\partial z^2}& =& \frac{dp}{dz}\\ \frac{\eta }{r}\left(-8v_{\mathrm{a}\mathrm{v}}\frac{r}{R^2}+\frac{\partial \left(r{v}_{z,\mathrm{transient}}\right)}{\partial r}\right)+\frac{{\partial }^2{v}_{z,\mathrm{transient}}}{\partial z^2}& =& \frac{dp}{dz}\\ -\frac{8\eta v_{\mathrm{a}\mathrm{v}}}{R^2}+\frac{\eta }{r}\left(\frac{\partial \left(r{v}_{z,\mathrm{transient}}\right)}{\partial r}\right)+\frac{{\partial }^2{v}_{z,\mathrm{transient}}}{\partial z^2}& =& \frac{dp}{dz}\end{array}$

From Eq. 16.34 we know that

$v_{\mathrm{a}\mathrm{v}}=-\frac{R^2}{8\eta }\frac{dp}{dz}$

which simplifies Eq. 18.70

$\begin{array}{lll}\frac{8\eta }{R^2}\frac{R^2}{8\eta }\frac{dp}{dz}+\frac{\eta }{r}\left(\frac{\partial \left(r{v}_{z,\mathrm{transient}}\right)}{\partial r}\right)+\frac{{\partial }^2{v}_{z,\mathrm{transient}}}{\partial z^2}& =& \frac{dp}{dz}\\ \frac{1}{r}\left(\frac{\partial \left(r{v}_{z,\mathrm{transient}}\right)}{\partial r}\right)+\frac{{\partial }^2{v}_{z,\mathrm{transient}}}{\partial z^2}& =& 0\\ \frac{1}{r}\left(\frac{\partial v_{z,\mathrm{transient}}}{\partial r}+r\frac{{\partial }^2{v}_{z,\mathrm{transient}}}{\partial r^2}\right)+\frac{{\partial }^2{v}_{z,\mathrm{transient}}}{\partial z^2}& =& 0\end{array}$

where we now, once again, have a homogeneous PDE. Please note that this is expected, as the steady-state solution must satisfy the Navier-Stokes equation (see Eq. 18.65) and this solution is effectively derived from this inhomogeneous PDE. This procedure is identical to the procedure we have taken when deriving the solution to the time-dependent transient problem in section 18.3. We now proceed to solving Eq. 18.71 again with a substitution and separation of variables approach.

18.5.3 Substitution and Separation of Variables

18.5.3.1 Approach

For Eq. 18.71 we assume a solution of the form

$v_{z,\mathrm{transient}}\left(r,z\right)=A\left(r\right)Z\left(z\right)$

which after inserting into Eq. 18.71 yields the following equation

$\begin{array}{lll}\frac{1}{r}\left(Z\frac{\partial A}{\partial r}+Xr\frac{{\partial }^{2}A}{\partial r^2}\right)+A\frac{{\partial }^{2}Z}{\partial z^2}& =& 0\\ \eta \frac{1}{A}\left(\frac{1}{r}\frac{\partial A}{\partial r}+\frac{{\partial }^{2}A}{\partial r^2}\right)& =& -\frac{1}{Z}\frac{{\partial }^{2}Z}{\partial z^2}=-\lambda \end{array}$

which yields the two ODEs

$\frac{1}{A}\left(\frac{1}{r}\frac{dA}{dr}+\frac{d^{2}A}{d{r}^2}\right)=-\lambda$

$\frac{1}{Z}\frac{d^{2}Z}{d{z}^2}=\lambda$

18.5.3.2 Solution for A (r)

The first ODE will look familiar by now if rewritten to

$\begin{array}{lll}\frac{dA}{dr}r+\frac{d^{2}A}{d{r}^2}{r}^2& =& -\lambda A{r}^2\\ r^2\frac{d^{2}A}{d{r}^2}+r\frac{dA}{dr}+\left(\lambda r^2-0\right)A& =& 0\end{array}$

which is a Bessel differential equation according to Eq. 3.44 with $\alpha =\sqrt{\lambda }$ and ν = 0 with the solution according to Eq. 3.46

$A\left(r\right)=c_0{J}_0\left(\sqrt{\lambda }r\right)$

18.5.3.3 Boundary Conditions

Obviously, we also need to apply the boundary conditions which are given by

$A\left(R\right)=c_0{J}_0\left(\sqrt{\lambda }R\right)\stackrel{!}{=}0\to \lambda ={\left(\frac{{\lambda }_i}{R}\right)}^2$

where λ_i is the i^th root of the Bessel function J_ν (see section 3.2.4.3). We therefore find the general solution for A (r) by

$A\left(r\right)={\displaystyle \sum _i^{\infty }{c}_i}{J}_0\left({\lambda }_i\frac{r}{R}\right)$

18.5.3.4 Solution for Z (z)

The general solution for Z (z) from Eq. 18.73 can be simply taken from Tab. 8.2 as

$Z\left(z\right)=c_1{\mathrm{e}}^{\frac{{\lambda }_{\mathrm{i}}}{\mathit{R}}z}+c_2{\mathrm{e}}^{-\frac{{\lambda }_{\mathrm{i}}}{\mathit{R}}z}$

We have one initial condition that will solve for one of the unknowns of Eq. 18.75. The second we can derive by stating that we expect the effects of z on the flow profile to diminish as z increases. This is due to the fact that the flow profile will eventually become fully-developed and thus independent of z. As for z → ∞ the first term of Eq. 18.75 will become large and we know that $c_1\stackrel{!}{=}0$ in order for the solution to be physically meaningful.

18.5.3.5 General Solution

We now have our general solution given by

$v_{z,\mathrm{transient}}\left(r,z\right)=A\left(r\right)Z\left(z\right)={\displaystyle \sum _i^{\infty }{c}_i}{J}_0\left({\lambda }_i\frac{r}{R}\right){\mathrm{e}}^{-\frac{{\lambda }_{\mathrm{i}}}{R}z}$

for which we obviously still have to determine the missing constants c_i. These we can derive by the initial conditions.

18.5.3.6 Applying Initial Conditions

We know that for z = 0 the overall velocity profile must be equal to the constant v_av. Therefore

$\begin{array}{lll}{v}_z\left(r,0\right)\hfill & =\hfill & v_{z,\mathrm{steady}‐\mathrm{state}}\left(r\right)+v_{z,\mathrm{transient}}\left(r,0\right)\hfill \\ \hfill & =\hfill & 2v_{\mathrm{a}\mathrm{v}}\left(1-{\left(\frac{r}{R}\right)}^2\right)+{\displaystyle \sum _i^{\infty }{c}_i}{J}_0\left({\lambda }_i\frac{r}{R}\right){\mathrm{e}}^{-\frac{{\lambda }_i}{\mathrm{R}}0}\hfill \\ \hfill & =\hfill & 2v_{\mathrm{a}\mathrm{v}}\left(1-{\left(\frac{r}{R}\right)}^2\right)+{\displaystyle \sum _i^{\infty }{c}_i}{J}_0\left({\lambda }_i\frac{r}{R}\right)\hfill \\ \hfill & \stackrel{!}{=}\hfill & v_{\mathrm{a}\mathrm{v}}\hfill \end{array}$

from which we derive

${\displaystyle \sum _i^{\infty }{c}_i}{J}_0\left({\lambda }_i\frac{r}{R}\right)=v_{\mathrm{a}\mathrm{v}}\left(2{\left(\frac{r}{R}\right)}^2-1\right)$

We now have to compare coefficients in order to derive c_i.

18.5.3.7 Fourier-Bessel Series Expansion

In order to find c_i we have to perform a Fourier-Bessel series expansions (see section 4.4). We know that the coefficients are determined by Eq. 4.51 as

$c_i=v_{\mathrm{a}\mathrm{v}}\frac{2}{R^2{J}_1^2\left({\lambda }_i\right)}{\displaystyle {\int }_0^{R}r}{J}_0\left({\lambda }_i\frac{r}{R}\right)\left(2{\left(\frac{r}{R}\right)}^2-1\right)dr$

Fortunately we have already derived the two terms of the integral during an introductory example in section 4.4.2.1. Using Eq. 4.58 we can immediately write down the first integral as

$2{\displaystyle {\int }_0^{R}r}{J}_0\left(\frac{{\lambda }_i}{R}r\right){\left(\frac{r}{R}\right)}^{2}dr=2\frac{R}{\frac{{\lambda }_{\mathit{i}}}{R}}{J}_1\left(\frac{{\lambda }_i}{R}R\right)-\frac{4}{{\left(\frac{{\lambda }_i}{R}\right)}^2}{J}_2\left(\frac{{\lambda }_i}{R}R\right)$

Using Eq. 4.57 we can write down the integral for the second term as

$-{\displaystyle {\int }_0^{R}r}{J}_0\left(\frac{{\lambda }_i}{R}r\right)dr=-\frac{R}{\frac{{\lambda }_{\mathit{i}}}{R}}{J}_1\left(\frac{{\lambda }_i}{R}R\right)$

Summing up Eq. 18.80 and Eq. 18.81 we find

$\begin{array}{lll}{\displaystyle {\int }_0^{R}r}{J}_0\left({\lambda }_i\frac{r}{R}\right)\left(2{\left(\frac{r}{R}\right)}^2-1\right)dr& =& 2\frac{R}{\frac{{\lambda }_i}{R}}{J}_1\left(\frac{{\lambda }_i}{R}R\right)-\frac{4}{{\left(\frac{{\lambda }_{\mathit{i}}}{R}\right)}^2}{J}_2\left(\frac{{\lambda }_i}{R}R\right)-\frac{R}{\frac{{\lambda }_i}{R}}{J}_1\left(\frac{{\lambda }_i}{R}R\right)\\ & =& \frac{R}{\frac{{\lambda }_i}{R}}{J}_1\left({\lambda }_i\right)-\frac{4}{{\left(\frac{{\lambda }_{\mathit{i}}}{R}\right)}^2}{J}_2\left({\lambda }_i\right)\end{array}$

Here we can apply the recurrence relation for expression the Bessel function of order 2, as Bessel functions of lower order according to Eq. 3.60 from which we find for ν = 1

$J_2\left({\lambda }_i\right)=\frac{2}{{\lambda }_i}\left(J_1\left({\lambda }_i\right)-J_0\left({\lambda }_i\right)\right)$

in which case we can rewrite Eq. 18.82 to

$\begin{array}{lll}{\displaystyle {\int }_0^{R}r}{J}_0\left({\lambda }_i\frac{r}{R}\right)\left(2{\left(\frac{r}{R}\right)}^2-1\right)dr& =& \frac{R}{\frac{{\lambda }_i}{R}}{J}_1\left({\lambda }_i\right)-\frac{8}{\frac{{\lambda }_i^3}{R^2}}\left(J_1\left({\lambda }_i^2\right)-J_0\left({\lambda }_i\right)\right)\\ & =& \frac{R^2}{{\lambda }_i^3}\left({\lambda }_i{J}_1\left({\lambda }_i^2\right)-8J_1\left({\lambda }_i\right)\right)\\ & =& \frac{R^2{J}_1\left({\lambda }_i\right)}{{\lambda }_i^3}\left({\lambda }_i^2-8\right)\end{array}$

where we have used the fact that λ_i are the roots of J₀ and thus J₀ (λ_i) = 0. We can now write down the missing coefficients c_i according to Eq. 18.79 as

$\begin{array}{lll}{c}_i& =& v_{\mathrm{a}\mathrm{v}}\frac{2}{R^2{J}_1^2\left({\lambda }_i\right)}\frac{R^2}{{\lambda }_i^3}\left({\lambda }_i^2{J}_1\left({\lambda }_i\right)-8J_1\left({\lambda }_i\right)\right)\\ & =& v_{\mathrm{a}\mathrm{v}}\frac{2}{{\lambda }_i^3{J}_1\left({\lambda }_i\right)}\left({\lambda }_i^2-8\right)\end{array}$

Eq. 18.84 gives us the coefficients we are seeking.