Plane Poiseuille Flow

January 2025
Fluid-Dynamics Physics

Plane Poiseuille flow  external link is the name given to fluid flowing in a 2 dimensional pipe. It is a useful system for testing CFD codes because it can be solved exactly in both steady-state and for the ‘ramp-up’ case, allowing us to test a static and transient (laminar) flow. Further, it allows testing of boundaries because the system contains non-slip boundaries at the walls, and (if driving with a pressure gradient) pressure boundary conditions at the in- and outlet.

Contents

Problem

Flow Geometry

Plane Poiseuille flow is flow in a 2 dimensional pipe. It can also be viewed as a slice through a 3 dimensional flow between parallel infinite planes. For our analysis let’s say the pipe has a length \(L,\) width \(W\) and is oriented along the \(x\)-axis, so that the plates are at \(y=0\) and \(y=W.\)

Flow Geometry
Flow Geometry

Pressure Drop vs Body Force

The flow can either be driven by a pressure drop \(\Delta P\) over the pipe length \(L\), or by an external body force such as gravity. The two views are equivalent and can be interconverted simply. To derive the conversion let’s imagine a slice of fluid, perpendicular to the pipe’s length, with width \(\delta x\). The mass in such a slice is \(m = \rho W \delta x\). Given the total pressure drop over the whole length \(L\) is \(\Delta P\) this tell us the pressure drop over the slice is \(\delta P = \delta x \Delta P / L\) thus the net force is

$$F = A\,\delta P = W \delta x \frac{\Delta P}{L}$$

and the net acceleration on the fluid slice given by

$$g = \frac{F}{m} = \frac{1}{\rho}\frac{\Delta P}{L}$$

As expected this ultimately does not depend on the width of the pipe, since the pressure is assumed to be applied evenly across the width. Nor does it depend on the slice width, since there was nothing special about the slice; the width is arbirary. Further, we didn’t even say where along the pipe’s length the slice was located. By symmetry it is clear that cannot make a difference.

Example

For water undergoing gravitational acceleration in a pipe we have $\rho = 1000 \,\text{kg}\,\text{m}^{-3}$ and $g = 9.81 \,\text{ms}^{-2}$, which corresponds to a pressure drop of $9.81 \,\text{kPa/m}$. The same calculation for air, however, gives at pressure drop of only $0.98 \,\text{kPa/m}$.

Steady-state Solution

The steady-state solution is a solution to the time-independent Navier-Stokes equation for a Newtonian fluid

$$0 = -\nabla p + \mu\nabla^2\mathbf u + \rho\mathbf g$$

Symmetry means that the solution cannot change along the pipe. This means that any gradients wrt $x$ must be zero. An immediate corollary of this is that we must also have $u_y = 0$ everywhere.

Proof

The flow is incompressible and we know gradients wrt $x$ are zero, so

$$0 = \nabla\cdot\mathbf u = \partial_x u_x + \partial_y u_y = \partial_y u_y$$

Additionally, the no-slip boundary conditions assert that $u_y(0) = u_y(W) = 0$. Combining these two facts shows that $u_y = 0$ everywhere.

Because only \(u_x\) is non-zero we can ignore \(u_y\) and henceforth we will just write \(u_x = u\) for brevity.

From above we know that we can write our driving force as either a pressure gradient or a constant body force density; wlog. let’s write it using the body force.

So we have

$$0 = \mu \partial_{yy} u + \rho g \tag{1}\label{eq:tins}$$

This implies the solution is quadratic in $y$. The Ansatz $u(y) = Ay(W - y)$, for some constant $A$, is quadratic and consistent with the no-slip boundary conditions. This gives

$$\partial_{yy}u = -2A$$

Thus

$$A = \frac{g \rho}{2\mu} = \frac{g}{2\nu}$$

giving

$$\boxed{ u(y) = \frac{g}{2\nu}y(W-y) }$$

Properties

The maximum velocity occurs along the midline $y = W/2$

$$u^\mathrm{max} = \frac{1}{8}\frac{g}{\nu}W^2$$

Sanity Check

Let’s check the units make sense. We should get $[u] = LT^{-1}$:

$$[u] = \Big[\frac{g\rho}{\mu} W^2\Big] = \frac{LT^{-2} \,ML^{-3}}{ML^{-1}T^{-1}}L^2 = LT^{-1}$$

It checks out.

We can use this to calculate the Reynold’s number  external link for the system

$$\mathrm{Re} = \frac{v^\mathrm{max} W}{\nu} = \frac{1}{8}g \nu^{-2}W^3$$

The volumetric flow rate is

$$Q = \int_0^W \!\!\! u(y) \,dy = \frac{2}{3}W u^\mathrm{max}$$

Time Dependent Solution

Next we want to know how the fluid behaves if we initially start it from rest and let it accelerate. We expect at long times the solution will approach the time-independent solution above, but what does it do before then?

Because we are not at steady state we are now interested in the usual time-dependent Navier-Stokes equation

$$ \rho \frac{\mathrm D u}{\mathrm D t} = \mu \frac{\partial^2 u}{\partial y^2} + \rho g $$

As above, we know from symmetry that the solution must have \(u_y(y, t) = 0\) and cannot depend on \(x\). From the \(x\)-translational symmetry we see that the material time derivative ends up being equal to the partial time derivative:

$$ \frac{\mathrm D u}{\mathrm D t} = \frac{\partial u}{\partial t} + \cancelto{0}{\frac{\partial u}{\partial x}}\frac{\partial x}{\partial t} = \frac{\partial u}{\partial t} $$

Using this, and dividing by $\rho$ gives

$$ \partial_t u = \nu \partial_{yy} u + g \tag{2}\label{eq:tdns1} $$

The solution method is relatively simple, but there is one slight trick. It turns out to be necessary to consider not the velocity \(u(y, t)\) itself, but the difference between it and the steady-state solution. Denote the steady-state solution as \(\bar u(y)\) and the difference as

$$v(y, t) = u(y, t) - \bar u(y)$$

If we sub \(u = \bar u + v\) into the Navier-Stokes equation \eqref{eq:tdns1} we find

$$\partial_t v = \nu \partial_{yy} v \tag{3}\label{eq:tdns2}$$
Proof
$$\begin{align} \partial_t u &= \nu \partial_{yy} u + g \\ \partial_t (v + \bar u) &= \nu \partial_{yy} (v + \bar u) + g & \qquad\text{Substitution} \\ \partial_t v + \cancel{\partial_t \bar u} &= \nu (\partial_{yy} v + \partial_{yy} \bar u) + g & \qquad\text{Linearity} \\ \partial_t v &= \nu \partial_{yy} v + \underbrace{\nu\partial_{yy} \bar u + g}_{=0} & \qquad\text{Soln. to \eqref{eq:tins}} \\ \end{align}$$

Given the solution is defined on a finite interval one may naturally think of performing a Fourier series  external link expansion. Specifically, our solution is defined on the interval \([0, W]\) which suggests we can use the half-range Fourier expansion  external link . Moreover, the no-slip boundary conditions tell us \(v(0, t) = v(W, t) = 0\) therefore we choose the odd continuation so that we get a decomposition in terms of sines, which we know will be zero at \(0\) and \(W\).

The fundamental wavelength of our expansion is \(2W\) so define the wavenumber of each harmonic as \(k_n = \pi n / W\) for brevity. The expansion then reads

$$v(y, t) = \frac{a_0(t)}{2} + \sum_{n=1}^\infty a_n(t)\sin(k_n y)$$

"(Odd) Half-range Fourier Series"
(Odd) Half-range Fourier Series

The Fourier expansion lets us decouple the spatial and temporal dependences, which plays very nicely with the derivatives appearing in \eqref{eq:tdns1} because of their linearity. Substituting in the Fourier expansion into \eqref{eq:tdns1} gives

$$\begin{align} \frac{1}{2}\partial_t a_0 +\vphantom{c} \sum_{n=1}^\infty \big(\partial_t a_n\big)\,\sin(k_n y) &= \nu \sum_{n=1}^\infty a_n \partial_{yy}\sin(k_n y) \\ &= \nu \sum_{n=1}^\infty a_n (-k_n^2) \sin(k_n y) \end{align}$$

where we have directly calculated the derivative wrt \(y\) on the right-hand side. Conveniently the second derivatives of sine is itself (with a minus sign), which means we can equate terms in the sum on both sides. This gives

$$\begin{align} \partial_t a_0 &= 0 \\ \partial_t a_n &= - \nu k_n^2 a_n \qquad n \geq 1 \end{align}$$

Nice! We’ve converted our PDE into an infinite series of ODEs. Thankfully these are nice and simple, and we get

$$\begin{align} a_0 &= A_0 \\ a_n &= A_n e^{- \nu k_n^2 t} \qquad n \geq 1 \end{align}$$

where \(A_0\) and \(A_n\) are integration constants. To pin down the values of the constants we can look at the initial and “final” values.

In the limit \(t\rightarrow\infty\) the coefficients \(a_n\) decay to zero exponentially leaving only \(A_0\). We know we must recover \(\bar u(y)\) so the difference must be zero and thus \(A_0 = 0\). The initial value is more interesting. Our initial conditions are \(u(y, 0) = 0\) so we have

$$v(y, 0) = u(y, 0) - \bar u(y) = -\bar u(y)$$

Also, setting \(t=0\) gives us \(a_n(0) = A_n\), thus we get

$$\sum_{n=1}^\infty A_n \sin(k_n y) = -\frac{g}{2\nu} y (W - y)$$

Notice this is just the Fourier expansion of (minus) the steady-state solution; no time dependence in sight. We do not need to worry about the constant here because we are using the odd expansion and so it is zero. The (non-constant) coefficients of the half-range Fourier expansion are given by

$$ A_n = \frac{2}{W} \int_0^W -\frac{g}{2\nu} y (W - y) \sin(k_n y)\,dy $$

You could solve this integral by hand, or you could ask your favourite symbolic algebra system. For example, I used the SymPy  external link Python package:

import sympy as sp

y, W, g, nu = sp.symbols("y W g nu", positive=True)
n           = sp.symbols("n",        positive=True, integer=True)

u = (g / (2 * nu)) * y * (W - y)
k = sp.pi * n / W
a_n = (2 / W) * sp.integrate(-u * sp.sin(k * y), (y, 0, W))

print(sp.latex(sp.simplify(a_n)))

which tells me

$$A_n = \frac{g}{\nu} \frac{2 W^2 ((-1)^n - 1)}{\pi^3 n^3} = \frac{g}{\nu} \frac{2 ((-1)^n - 1)}{k_n^3 W}$$

Notice that if \(n\) is even then \(A_n = 0\) so we can rewrite the sum to just exclude these terms. Let \(n=2m+1\), with \(m\) running from zero, to give

$$A_m = -\frac{g}{\nu}\frac{4}{k_{2m+1}^3 W}$$

And there we have it, the final answer is thus

$$\boxed{ u(y, t) = \frac{g}{2\nu} y (W - y) - \sum_{m=0}^\infty \frac{g}{\nu} \frac{4}{k_{2m+1}^3 W} \exp(-\nu k_{2m+1}^2 t) \sin(k_{2m+1} y) }$$

It looks somewhat daunting but it’s not so bad. It’s easy enough to calculate numerically because the terms in the expansion decay in magnitude as \(O(m^{-3})\), so we only need a few terms to get a good approximation.

Sanity Checks

We can verify easily enough that this is indeed a solution to \eqref{eq:tdns1} simply by plugging it in, and clearly it meets the boundary conditions and the “final condition”. The initial condition is less obvious but we can see numerically it does satisfy these.

One thing we may expect is that exactly at time zero all parts of the bulk of the fluid should start accelerating with \(g\). This is because the viscous effects have had no time to propagate (diffuse) from the boundaries. For general \(y\) it seems tricky to show that the acceleration at time zero evalutates to \(g\), but for the special value of \(y=W/2\) we can show this is the case.

Proof

First notice

$$\sin\Big(k_{2m+1} \frac W 2\Big) = \sin\Big(\frac{(2m+1)\pi}{W} \frac W 2 \Big) = \sin\Big(m + \frac{1}{2}\Big) \equiv (-1)^m$$

Then differentiating \(u\) and setting \(t=0\) gives

$$\partial_t u(W/2, 0) = g \sum_{m=0}^\infty \frac{4}{k_{2m+1} W} (-1)^m = g \frac{4}{\pi} \sum_{m=0}^\infty \frac{(-1)^m}{(2m+1)}$$

As we hope, the sum equals \(\pi/4\) (again, I asked SymPy) and thus, as expected,

$$\partial_t u(W/2, 0) = g$$

Properties

Characteristic Time

So what does the solution look like? We can see that it tends to the steady-state solution at long times, but more specifically we approach this because the Fourier components decay exponentially. The half-life of each component is

$$\tau_m = \frac{\ln(2)}{\nu k_{2m+1}^2}= \tau_0 (2m+1)^{-2} \qquad\text{with}\qquad\tau = \frac{\ln(2) W^2}{\nu\pi^2}$$

As it is the weights decrease quickly for the higher frequency components, additionally they decay much faster, so at even moderate times the solution is dominated by the low frequency components.

The form of \(\tau_m\) shows us that changing \(W\) and \(\nu\) change the rate at which we approach the steady-state solution, but does not otherwise change the shape of the solution. This is because of the way \(\tau\) factors out for all \(\tau_m\), so changing either \(W\) or \(\nu\) will change all half-lives by the same factor. For this reason is is natural to plot the time evolution in terms of \(\tau_0\) which is a characteristic time for our solution. The following plot shows the time evolution over a ’long’ time as \(u\) starts to approach \(\bar u\):

“Time-dependent Solution”
Analytical Solution (Long times)

Plug Flow

At short times the bulk of the fluid has not yet felt the effects of the no-slip boundaries. As such, away from the boundaries, the fluid moves as a plug  external link . Over time the viscous effects spread inwards from the boundaries causing deviation from this and letting the fluid reach its terminal profile. The following plot shows this flow feature at short times:

“Time-dependent Solution”
Analytical Solution (Short times)

Because the flow is incompressible, at any non-zero time, no matter how small, there is technically some viscous effect felt everywhere in the domain. Intuitively you can think about this being because viscosity is diffusive, and the analytical solution for the diffusion equation  external link gives a non-zero value everywhere instantly. We can see that our solution also suggests this because the Fourier components will have always decayed some non-zero amount for any \(t > 0\).

Flow Rate

The time dependence of the flow rate is easy enough to calculate by integrating our solution wrt \(y\). If we calculate the flow rate relative to the steady-state flow rate we find

$$\frac{Q(t)}{\bar Q} = 1 - \frac{96}{\pi^4} \sum_{m=0}^\infty (2 m + 1)^{-4} \exp(-\nu k_{2 m + 1}^2 t)$$

where \(\bar Q\) is the steady-state flow-rate from above.

Given the speed with which the terms drop off, this is reasonably well approximated (to within 2%) by just the first term

$$\frac{Q(t)}{\bar Q} \approx 1 - \frac{96}{\pi^4} \exp\Big(-\frac{\log(2)}{\tau} t\Big)$$

This will overestimate the flow rate at early times, for example it gives a small positive value for \(t=0\) when we know it should be zero. This tells us that the flow rate expnentially approaches its steady-state value as time goes on.

Fluid Motion from the Boltzmann Equation

Latest Posts