Poiseuille Flow: Analytical Solution Derivation and Finite Difference Validation

Poiseuille flow is a fully developed laminar flow driven by a constant pressure gradient between two parallel plates or inside a circular pipe. Since an exact analytical solution exists, it appears in textbooks worldwide as the very first test case for CFD code validation.

Today, we will derive the analytical solution for channel flow (between 2D plates) step-by-step and verify the convergence order by writing a simple finite difference (FD) code.

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$:

0 = -\frac{C}{2} h^2 + A_1 h + A_2 \tag{1}

$u(-h) = 0$:

0 = -\frac{C}{2} h^2 - A_1 h + A_2 \tag{2}

(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.

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.

This is precisely why Poiseuille flow is an ideal CFD validation case — it provides near-perfect answers even at low grid resolutions, allowing for a quick check of the code's fundamental accuracy.

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

Validation Checklist#

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$.

Conclusion#

Poiseuille flow may seem simple, but this benchmark allows you to verify quite a few things:

• Whether the second-order central difference scheme is correctly implemented.
• Whether wall no-slip and pressure boundary conditions are handled properly.
• Whether the viscous stress tensor ($\mu \nabla^2 \mathbf{u}$) is accurately calculated in the code.

The principle of validating starting from problems with analytical solutions (V&V) is a shared culture across all CFD fields, including aerospace, nuclear, and meteorology. In the next post, let's challenge ourselves with unsteady analytical solution validation using Stokes' First Problem (impulsively started plate).

“

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

Adjust pressure gradient and viscosity to watch the parabola steepen.