Poiseuille Flow: Analytical Solution Derivation and Finite Difference Validation

Two plates, gap $2h$. Pressure drops only in $x$. Fully developed means $v=0$ and $u=u(y)$. A round pipe carries the same name, Poiseuille flow, but the formulas below are the plane channel. If a code misses this parabola, the viscous term or the no-slip wall is wrong.

Move the −dp/dx or μ slider and only the printed $u_{\max}$ changes. The drawn curve does not.

The curve on the canvas is $u/u_{\max} = 1 - (y/h)^2$. The Re slider only shifts the corner label. With $h=1$, the peak readout is the −dp/dx value divided by $2\mu$.

Problem Setup#

Geometry and Boundary Conditions#

Consider a fully developed laminar flow between two infinite parallel plates.

• Channel half-width: $h$ (plates are located at $y = \pm h$)
• Flow direction: $x$ (assumed infinitely long)
• Boundary conditions: no-slip at both walls $\Rightarrow$ $u(\pm h) = 0$
• Driving force: constant pressure gradient $\dfrac{dp}{dx} = -G < 0$ ($G > 0$ is a constant)

Governing Equation#

Since the flow is fully developed, $\partial/\partial x = 0$ and $v = 0$. The $x$-direction Navier–Stokes equation is:

$\rho \left( u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} \right) = -\frac{dp}{dx} + \mu \frac{\partial^2 u}{\partial y^2}$

Applying the fully developed conditions ($\partial u/\partial x = 0$, $v = 0$):

$\mu \frac{d^2 u}{dy^2} = \frac{dp}{dx} = -G$

Analytical Derivation#

Step 1: Simplify to an ODE#

$\frac{d^2 u}{dy^2} = -\frac{G}{\mu} \equiv -C \quad (C > 0)$

Step 2: First Integration#

$\frac{du}{dy} = -C y + A_1$

Step 3: Second Integration#

$u(y) = -\frac{C}{2} y^2 + A_1 y + A_2$

Step 4: Apply Boundary Conditions#

$u(+h) = 0$:

$u(-h) = 0$:

(1) $-$ (2): $0 = 2 A_1 h \Rightarrow A_1 = 0$ (Symmetry)

(1) $+$ (2): $0 = -C h^2 + 2A_2 \Rightarrow A_2 = \dfrac{C h^2}{2}$

Final Analytical Solution#

$\boxed{u(y) = \frac{G}{2\mu}\bigl(h^2 - y^2\bigr)}$

Maximum velocity at the channel center ($y=0$):

$u_{\max} = \frac{G h^2}{2\mu}$

Physical Interpretation#

Velocity Profile#

The analytical solution is a parabolic distribution. Wall friction is transmitted uniformly throughout the flow cross-section, resulting in a maximum at the center and zero at the walls.

The dimensionless velocity $u/u_{\max} = 1 - (y/h)^2$ is a universal profile, independent of the channel half-width $h$ or the viscosity $\mu$.

Wall Shear Stress#

$\tau_w = \mu \left.\frac{du}{dy}\right|_{y=h} = \mu \cdot \frac{G}{\mu}(-h) = -G h$

(Magnitude: $\tau_w = G h$) The larger the pressure gradient and the wider the channel, the greater the wall friction.

Volumetric Flow Rate#

$Q = \int_{-h}^{h} u \, dy = \frac{G}{2\mu} \int_{-h}^{h} (h^2 - y^2) \, dy = \frac{G}{2\mu} \cdot \frac{4h^3}{3} = \frac{2 G h^3}{3\mu}$

This is known as the 2D version of the Hagen-Poiseuille law. Since the flow rate is proportional to $h^3$, doubling the channel width increases the flow rate eightfold.

Code defaults $h=G=\mu=1$ give $u_{\max}=1/2$ and $Q=2/3$. Divide by the height $2h$ and the mean is $1/3$, two thirds of the peak. Integrate $1-(y/h)^2$ from $-h$ to $h$ and divide by $2h\,u_{\max}$: the same $2/3$. $h$ and $\mu$ cancel.

Double only $h$: $u_{\max}$ grows by 4 ($h^2$), $Q$ by 8 ($h^3$), wall stress $Gh$ by 2. Treat the stress like the flow-rate ratio and the wall value is four times too high. The laminar pipe that locks the kinetic-energy factor at 2.0000 is the coefficient in front of velocity head.

Visualizing the Analytical Solution with Python#

import numpy as np
import matplotlib.pyplot as plt
 
# Parameters
h   = 1.0    # Channel half-width [m]
G   = 1.0    # Pressure gradient [Pa/m]
mu  = 1.0    # Dynamic viscosity [Pa·s]
 
# Grid
y = np.linspace(-h, h, 400)
u_exact = G / (2 * mu) * (h**2 - y**2)
u_max   = G * h**2 / (2 * mu)
 
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
 
# Velocity Profile (Horizontal)
axes[0].plot(u_exact / u_max, y / h, "k-", lw=2)
axes[0].fill_betweenx(y / h, u_exact / u_max, alpha=0.15)
axes[0].axvline(0, color="gray", lw=0.8, ls="--")
axes[0].set_xlabel(r"$u / u_{\max}$")
axes[0].set_ylabel(r"$y / h$")
axes[0].set_title("Poiseuille Velocity Profile")
axes[0].set_xlim(-0.05, 1.1)
 
# Shear Stress Distribution
tau = mu * np.gradient(u_exact, y)
axes[1].plot(tau / G / h, y / h, "r-", lw=2, label=r"$\tau / \tau_w$")
axes[1].axvline(0, color="gray", lw=0.8, ls="--")
axes[1].set_xlabel(r"$\tau / \tau_w$")
axes[1].set_ylabel(r"$y / h$")
axes[1].set_title("Shear Stress Distribution (Linear)")
axes[1].legend()
 
plt.tight_layout()
plt.savefig("poiseuille_profile.png", dpi=150)
plt.show()

Numerical Validation with Finite Difference Method#

We discretize $d^2u/dy^2 = -G/\mu$ using a second-order central difference scheme. On a uniform grid $\Delta y = 2h/(N-1)$:

$\frac{u_{i+1} - 2u_i + u_{i-1}}{(\Delta y)^2} = -\frac{G}{\mu}$

Boundary conditions: $u_0 = u_{N-1} = 0$.

Writing this as a tridiagonal system of linear equations:

$\begin{pmatrix} -2 & 1 \\ 1 & -2 & 1 \\ & \ddots & \ddots & \ddots \\ & & 1 & -2 \end{pmatrix} \mathbf{u} = -\frac{G (\Delta y)^2}{\mu} \mathbf{1}$

Python Implementation#

import numpy as np
 
def solve_poiseuille_fd(N, h=1.0, G=1.0, mu=1.0):
    """
    N : Number of internal grid points (excluding boundaries)
    returns (y, u_num, u_exact, L2_error)
    """
    dy = 2 * h / (N + 1)
    y_inner = np.linspace(-h + dy, h - dy, N)
 
    # Tridiagonal coefficient matrix
    diag  = -2 * np.ones(N)
    upper =  np.ones(N - 1)
    lower =  np.ones(N - 1)
    A = (np.diag(diag) + np.diag(upper, k=1) + np.diag(lower, k=-1))
 
    rhs = -G * dy**2 / mu * np.ones(N)
 
    u_num = np.linalg.solve(A, rhs)
 
    # Full grid including boundaries
    y_full  = np.concatenate([[-h], y_inner, [h]])
    u_full  = np.concatenate([[0.0], u_num, [0.0]])
    u_exact = G / (2 * mu) * (h**2 - y_full**2)
 
    L2 = np.sqrt(np.mean((u_full - u_exact)**2))
    return y_full, u_full, u_exact, L2
 
# Grid convergence test
Ns = [4, 8, 16, 32, 64, 128, 256]
errors = []
for N in Ns:
    _, _, _, err = solve_poiseuille_fd(N)
    errors.append(err)
 
# Convergence order
for i in range(1, len(Ns)):
    ratio = errors[i-1] / errors[i]
    order = np.log2(ratio)
    print(f"N={Ns[i]:4d}  L2={errors[i]:.3e}  Convergence Order={order:.2f}")

Convergence Results#

“

Note: Since Poiseuille flow is a quadratic polynomial, the second-order central difference scheme yields results at the level of machine precision. The actual error does not scale with $\Delta y^2$ but saturates at the floating-point rounding error level.

At low $N$ the parabola already sits on the formula. That is enough to see whether the viscous term is alive. The order column cannot make that call.

The discretization line writes $\Delta y = 2h/(N-1)$. The table does not. At $N=4,8,16,32,64$ the spacings are 0.4000, 0.2222, 0.1176, 0.0606, 0.0308. With $h=1$ those equal $2/(N+1)$: $2/5=0.4$, $2/9 \approx 0.2222$, $2/17 \approx 0.1176$, $2/33 \approx 0.0606$, $2/65 \approx 0.0308$. Both walls are dropped from the unknowns, so the spacing is $2h/(N+1)$.

The five $L_2$ entries are $2.78\times 10^{-16}$, $1.67\times 10^{-16}$, $2.22\times 10^{-16}$, $2.78\times 10^{-16}$, $2.22\times 10^{-16}$. Successive base-2 logs are about 0.7, then negative, negative, then about 0.3. Not ~2.0. A quadratic in $y$ has zero third derivative, so the second-order truncation error is zero. This column does not measure order.

Comparison Plot of Results#

import matplotlib.pyplot as plt
 
fig, ax = plt.subplots(figsize=(6, 5))
 
colors = plt.cm.viridis(np.linspace(0.2, 0.9, 4))
Ns_plot = [4, 8, 16, 64]
 
for N, c in zip(Ns_plot, colors):
    y, u_num, _, _ = solve_poiseuille_fd(N)[:3]
    _, _, u_exact, _ = solve_poiseuille_fd(N)
    ax.plot(u_num, y, "o--", color=c, ms=4, label=f"FD N={N}")
 
# Analytical solution
y_fine = np.linspace(-1, 1, 400)
ax.plot(0.5 * (1 - y_fine**2), y_fine, "k-", lw=2, label="Analytical Solution")
 
ax.set_xlabel(r"$u$ (Non-dimensional)")
ax.set_ylabel(r"$y/h$")
ax.set_title("Comparison: FD vs Analytical Solution")
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig("poiseuille_fd_vs_exact.png", dpi=150)
plt.show()

OpenFOAM Setup Hints#

When setting up Poiseuille flow with OpenFOAM's icoFoam, the key parameters are:

0/U — Velocity Boundary Conditions#

boundaryField
{
    walls
    {
        type    noSlip;         // No-slip at wall
    }
    inlet
    {
        type    fixedValue;
        value   uniform (1 0 0); // Or parabolicVelocity
    }
    outlet
    {
        type    zeroGradient;
    }
}

0/p — Pressure Boundary Conditions#

boundaryField
{
    inlet   { type fixedValue; value uniform 1; }
    outlet  { type fixedValue; value uniform 0; }
    walls   { type zeroGradient; }
}

constant/transportProperties#

nu              nu [ 0 2 -1 0 0 0 0 ] 0.01;  // Adjust kinematic viscosity

Entrance length and centerline speed#

1. Fully Developed Condition: The channel length must be sufficiently longer than the entrance development length ($L_{dev} \approx 0.05 \, Re \cdot h$).
2. Grid Independence: A grid count of 16 or more in the $y$-direction ensures second-order accuracy.
3. Validation Metrics: Compare the centerline velocity $u_c = G h^2 / (2\mu)$ and wall shear stress $\tau_w = G h$.

Stop reading the order column once L2 is 1e-16?#

At $N=4$, $L_2$ is already $2.78\times 10^{-16}$. Pushing the grid to 64 only jiggles that column inside the $10^{-16}$ band. The printed ~2.0 did not come from those ratios. On a quadratic the central scheme has zero truncation error, so signing off the order from this table is the wrong call. Whether $\mu \nabla^2 \mathbf{u}$ is alive shows up in the centerline speed and the wall stress.

Check $Gh^2/(2\mu)$ and $Gh$. Order that actually splits lives on a solution with a discontinuity. The same code reporting 0.79 and 2.00 at once is where that column stops meaning the scheme order. Stuck on roundoff, a finer grid here still does not produce an order.

“

References

• Batchelor, G. K. An Introduction to Fluid Dynamics. Cambridge, 1967.
• Ferziger, J. H., Perić, M. Computational Methods for Fluid Dynamics. Springer, 2002.
• OpenFOAM Documentation — icoFoam tutorial: cavity