Planar Poiseuille Flow Tutorial

This tutorial documents the analytical solution implemented for incompressible planar Poiseuille flow with either constant-viscosity or power-law rheology. The solution gives the streamwise velocity and temperature profiles across a planar channel.

The implementation assumes the exact pressure is zero because periodic boundary conditions are used in the streamwise direction.


1. Problem setup

Consider fully developed planar Poiseuille flow between two parallel plates. The cross-channel coordinate is \(y\), and the channel walls are located at

\[ y = -H, \qquad y = H, \]

where

\[ H = \frac{h_{\max} - h_{\min}}{2}. \]

The flow is driven by a uniform streamwise momentum source \(S_u\). The velocity field is

\[\begin{split} \mathbf{u}(y) = \begin{bmatrix} u(y) \\ 0 \end{bmatrix} \end{split}\]

in two dimensions, or

\[\begin{split} \mathbf{u}(y) = \begin{bmatrix} u(y) \\ 0 \\ 0 \end{bmatrix} \end{split}\]

in three dimensions.

The wall temperatures are prescribed as

\[ T(-H) = T_l, \qquad T(H) = T_u. \]

The exact Lagrange pressure is

\[ p(y) = 0. \]

2. Parameters

Symbol

Code parameter

Description

\(h_{\min}\)

H min

Lower geometric bound used to define the channel height

\(h_{\max}\)

H max

Upper geometric bound used to define the channel height

\(H\)

computed

Channel half-height, \(H=(h_{\max}-h_{\min})/2\)

\(T_l\)

Lower wall temperature

Lower-wall temperature

\(T_u\)

Upper wall temperature

Upper-wall temperature

\(S_u\)

Momentum source

Streamwise momentum source

\(\nu\)

Kinematic viscosity or Flow consistency index

Constant kinematic viscosity or power-law consistency index

\(n\)

Flow behavior index

Power-law index; \(n=1\) for a Newtonian fluid

\(k\)

Thermal conductivity

Thermal conductivity

For a constant-property Newtonian fluid, the implementation sets

\[ n = 1. \]

For a power-law fluid, \(n\) is read from the input as the flow behavior index.


3. Analytical velocity solution

The streamwise velocity profile is

\[ u(y) = \frac{n}{n+1} \left(\frac{S_u}{\nu}\right)^{1/n} H^{(n+1)/n} \left[ 1 - \left(\frac{|y|}{H}\right)^{(n+1)/n} \right]. \]

The transverse velocity components are zero:

\[ v(y) = 0, \]

and, in three dimensions,

\[ w(y) = 0. \]

Newtonian limit

For \(n=1\), the velocity reduces to the standard parabolic Poiseuille profile:

\[ u(y) = \frac{1}{2} \frac{S_u}{\nu} \left( H^2 - y^2 \right). \]

4. Analytical temperature solution

The analytical temperature profile is

\[\begin{split} T(y) = \frac{n^2}{(3n+1)(2n+1)} \frac{\nu}{k} \left(\frac{S_u}{\nu}\right)^{(n+1)/n} \left[ H^{(3n+1)/n} - |y|^{(3n+1)/n} \right] \\ + \frac{1}{2} \frac{y}{H} \left(T_u - T_l\right) + \frac{1}{2} \left(T_u + T_l\right). \end{split}\]

This expression contains two contributions:

  1. A viscous-heating contribution caused by the imposed momentum source.

  2. A linear conductive contribution satisfying the imposed wall temperatures.

The boundary values are recovered directly:

\[ T(-H) = T_l, \qquad T(H) = T_u. \]

Newtonian limit

For \(n=1\),

\[ T(y) = \frac{1}{20} \frac{\nu}{k} \left(\frac{S_u}{\nu}\right)^2 \left( H^4 - y^4 \right) + \frac{1}{2} \frac{y}{H} \left(T_u - T_l\right) + \frac{1}{2} \left(T_u + T_l\right). \]

5. Python script for plotting the analytical solution

Save the following as plot_planar_poiseuille_analytical_solution.py and run

python plot_planar_poiseuille_analytical_solution.py

The script generates velocity and temperature plots and saves them as PNG files.

import numpy as np
import matplotlib.pyplot as plt


def planar_poiseuille_solution(
    y,
    h_min=-1.0,
    h_max=1.0,
    T_l=300.0,
    T_u=310.0,
    S_u=1.0,
    nu=1.0,
    k=1.0,
    n=1.0,
):
    """
    Analytical solution for planar Poiseuille flow with viscous heating.

    Parameters
    ----------
    y : array_like
        Cross-channel coordinate.
    h_min, h_max : float
        Geometric bounds used to define the channel half-height,
        H = (h_max - h_min) / 2.
    T_l, T_u : float
        Lower- and upper-wall temperatures.
    S_u : float
        Streamwise momentum source.
    nu : float
        Kinematic viscosity for n = 1, or power-law consistency index for n != 1.
    k : float
        Thermal conductivity.
    n : float
        Flow behavior index. Use n = 1 for a Newtonian fluid.

    Returns
    -------
    u : ndarray
        Streamwise velocity.
    T : ndarray
        Temperature.
    p : ndarray
        Exact Lagrange pressure, equal to zero.
    """
    y = np.asarray(y, dtype=float)
    H = 0.5 * (h_max - h_min)

    if H <= 0.0:
        raise ValueError("The channel half-height H must be positive.")
    if nu <= 0.0:
        raise ValueError("nu must be positive.")
    if k <= 0.0:
        raise ValueError("k must be positive.")
    if n <= 0.0:
        raise ValueError("n must be positive.")
    if np.any(np.abs(y) > H):
        raise ValueError("All y values must satisfy |y| <= H.")

    abs_y = np.abs(y)

    u = (
        n / (n + 1.0)
        * (S_u / nu) ** (1.0 / n)
        * H ** ((n + 1.0) / n)
        * (1.0 - (abs_y / H) ** ((n + 1.0) / n))
    )

    T = (
        n**2
        / ((3.0 * n + 1.0) * (2.0 * n + 1.0))
        * nu
        / k
        * (S_u / nu) ** ((n + 1.0) / n)
        * (H ** ((3.0 * n + 1.0) / n) - abs_y ** ((3.0 * n + 1.0) / n))
    )
    T += 0.5 * y * (T_u - T_l) / H + 0.5 * (T_u + T_l)

    p = np.zeros_like(y)

    return u, T, p


def main():
    # Example parameters. Modify these to match your input deck.
    h_min = -1.0
    h_max = 1.0
    T_l = 300.0
    T_u = 310.0
    S_u = 1.0
    nu = 1.0
    k = 1.0
    n = 1.0

    H = 0.5 * (h_max - h_min)
    y = np.linspace(-H, H, 500)

    u, T, p = planar_poiseuille_solution(
        y,
        h_min=h_min,
        h_max=h_max,
        T_l=T_l,
        T_u=T_u,
        S_u=S_u,
        nu=nu,
        k=k,
        n=n,
    )

    plt.figure()
    plt.plot(u, y)
    plt.xlabel("Streamwise velocity, u(y)")
    plt.ylabel("Cross-channel coordinate, y")
    plt.title("Analytical velocity profile")
    plt.grid(True)
    plt.tight_layout()
    plt.savefig("analytical_velocity_profile.png", dpi=300)

    plt.figure()
    plt.plot(T, y)
    plt.xlabel("Temperature, T(y)")
    plt.ylabel("Cross-channel coordinate, y")
    plt.title("Analytical temperature profile")
    plt.grid(True)
    plt.tight_layout()
    plt.savefig("analytical_temperature_profile.png", dpi=300)

    plt.show()


if __name__ == "__main__":
    main()


6. Notes for verification

Useful checks are:

\[ u(-H) = u(H) = 0, \]
\[ T(-H) = T_l, \qquad T(H) = T_u, \]
\[ p(y) = 0. \]

For \(n=1\), the velocity profile should be parabolic and the viscous-heating contribution to temperature should be quartic in \(y\).


7. Running this case with VERTEX-CFD on the Helios CPU platform

The VERTEX-CFD input file for this analytical solution is located at

vertex-cfd/examples/inputs/incompressible/incompressible_2d_planar_poiseuille.xml

These instructions assume that VERTEX-CFD has already been built on the Helios CPU platform and that the Helios setup environment is available at

/software/build/vertex-cfd/setup-env.sh

Runs may be performed from your home directory.

7.1. Stage the input file

Create a run directory and copy the planar Poiseuille input file into it:

mkdir -p ~/planar_poiseuille
cd ~/planar_poiseuille

cp /path/to/vertex-cfd/examples/inputs/incompressible/incompressible_2d_planar_poiseuille.xml .

If you are working from a clone of the VERTEX-CFD repository, replace /path/to/vertex-cfd with the path to your source tree.

The copied input file should be named

incompressible_2d_planar_poiseuille.xml

The submission script below sets

BASEFILE=incompressible_2d_planar_poiseuille

so the input file is resolved as

INPUTFILE=${BASEFILE}.xml

7.2. Create a SLURM submission script

Create a file named

touch run-planar-poiseuille-on-helios-cpu.sh

and copy the following contents into it:

#!/bin/bash
#SBATCH -N 1
#SBATCH --ntasks-per-node 2
#SBATCH --cpus-per-task 1
#SBATCH --output output%j.txt
#SBATCH --error error%j.txt
#SBATCH -J planar-poiseuille
#SBATCH -t 1-00:00
#SBATCH --mail-type=ALL
#SBATCH --mail-user=YOUR-EMAIL@ornl.gov

# Load the VERTEX-CFD Helios CPU environment.
source /software/build/vertex-cfd/setup-env.sh
module list

export OMP_PROC_BIND=true
export OMP_PLACES=threads

# Input file for this verification case.
# This assumes the job is submitted from the directory containing the XML file.
BASEFILE=incompressible_2d_planar_poiseuille
INPUTFILE=${BASEFILE}.xml

# Update this path to point to your VERTEX-CFD executable.
# For example, this may be in your build or install directory.
VERTEXCFD=/PATH/TO/BUILD/DIRECTORY/src/vertexcfd

mpirun $VERTEXCFD --i=$INPUTFILE

# Recombine decomposed Exodus output files.
# Each MPI rank writes one file during the simulation.
EPU=/software/build/vertex-cfd/trilinos/16.0.0/bin/epu
$EPU -auto ${BASEFILE}_solution.exo.${SLURM_NTASKS}.000

# Remove decomposed files after successful recombination.
rm ${BASEFILE}_solution.exo.${SLURM_NTASKS}.*

Before submitting, edit the following entries as needed:

  • #SBATCH -N 2: number of nodes.

  • #SBATCH --ntasks-per-node 185: MPI ranks per node.

  • #SBATCH --mail-user=YOUR-EMAIL@ornl.gov: your email address.

  • VERTEXCFD=/PATH/TO/BUILD/DIRECTORY/src/vertexcfd: path to your VERTEX-CFD executable.

For a small verification run, you may want to reduce the node count and MPI rank count to match the size of the mesh and the resources available on Helios.

7.3. Submit the job

Submit the job using

sbatch run-planar-poiseuille-on-helios-cpu.sh

During the run, VERTEX-CFD writes decomposed Exodus solution files, one per MPI rank. After the solver finishes, the script recombines them using Trilinos epu:

/software/build/vertex-cfd/trilinos/16.0.0/bin/epu -auto incompressible_2d_planar_poiseuille_solution.exo.${SLURM_NTASKS}.000

The script then removes the decomposed per-rank files:

rm incompressible_2d_planar_poiseuille_solution.exo.${SLURM_NTASKS}.*

If the job exits before recombination, you can manually rerun the epu command after checking how many decomposed files were written.

7.4. Expected output

For this case, the main recombined solution file should have a name similar to

incompressible_2d_planar_poiseuille_solution.exo

depending on the output settings in the XML input file.

You can inspect the solution in ParaView, for example on a Helios visualization node, and compare the numerical fields against the analytical profiles documented above:

\[ u(y) = \frac{n}{n+1} \left(\frac{S_u}{\nu}\right)^{1/n} H^{(n+1)/n} \left[ 1 - \left(\frac{|y|}{H}\right)^{(n+1)/n} \right], \]
\[\begin{split} T(y) = \frac{n^2}{(3n+1)(2n+1)} \frac{\nu}{k} \left(\frac{S_u}{\nu}\right)^{(n+1)/n} \left[ H^{(3n+1)/n} - |y|^{(3n+1)/n} \right] \\ + \frac{1}{2} \frac{y}{H} \left(T_u - T_l\right) + \frac{1}{2} \left(T_u + T_l\right). \end{split}\]

The exact Lagrange pressure for this periodic setup is

\[ p(y)=0. \]