Law of Hagen-Poiseuille

Module: tutorials.06_fluid.01_HagenPouseille_Stationary.run

Section author: Elias Karabelas <elias.karabelas@medunigraz.at>

This example demonstrates a simple application for fluid dynamics in a straight cylindrical pipe.

Law of Hagen-Poiseuille

The law of Hagen-Poiseuille is a physical law that gives the pressure drop \Delta p for an incompressible Newtonian fluid in the laminar regime flowing through a long cylindrical pipe of constant cross section. This geometric setup is depicted in figure Fig. 166.

../../_images/HagenPouseilleSetup.png

Fig. 166 Geometry for the Hagen-Poiseuille law.

The law states that

\Delta p = \frac{8\mu LQ}{\pi R^4}

where

  • \Delta p is the pressure drop (in Pascal)
  • L is the length of the pipe (in meter)
  • \mu is the dynamic viscosity of the fluid (in Pascal seconds)
  • Q is the volumetric flow rate (in cubic meter per seconds)
  • R is the pipe radius (in meter)

Derivation of Hagen-Poiseuille’s Law

For the derivation we will use the Navier-Stokes equations

\rho \partial_t \vec u + \rho \nabla \vec u \vec u - \mu \Delta \vec u + \nabla p &= \vec f,\\
\nabla \cdot \vec u &= 0.

First we assume that \vec f = \vec 0 and \rho constant. To derive Hagen-Poiseuille’s law we will rewrite the Navier-Stokes equations in cylindrical coordinates

(65)\rho \left(\frac{\partial u_r}{\partial t} + u_r \frac{\partial u_r}{\partial r} + \frac{u_{\phi}}{r}
 \frac{\partial u_r}{\partial \phi} + u_z \frac{\partial u_r}{\partial z} - \frac{u_{\phi}^2}{r}\right) &=
 -\frac{\partial p}{\partial r} + \mu \left(\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial
 u_r}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u_r}{\partial \phi^2} + \frac{\partial^2 u_r}{\partial z^2} - \frac{u_r}{r^2} -
     \frac{2}{r^2}\frac{\partial u_\phi}{\partial \phi} \right),\\
 \rho \left(\frac{\partial u_{\phi}}{\partial t} + u_r \frac{\partial u_{\phi}}{\partial r} +
    \frac{u_{\phi}}{r} \frac{\partial u_{\phi}}{\partial \phi} + u_z \frac{\partial u_{\phi}}{\partial z} + \frac{u_r u_{\phi}}{r}\right) &= -\frac{1}{r}\frac{\partial p}{\partial \phi} + \mu \left(\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_{\phi}}{\partial r}\right) +
    \frac{1}{r^2}\frac{\partial^2 u_{\phi}}{\partial \phi^2} + \frac{\partial^2 u_{\phi}}{\partial z^2} +
    \frac{2}{r^2}\frac{\partial u_r}{\partial \phi}-\frac{u_{\phi}}{r^2}\right),\\
 \rho \left(\frac{\partial u_z}{\partial t} + u_r \frac{\partial u_z}{\partial r} + \frac{u_{\phi}}{r} \frac{\partial u_z}{\partial \phi} +
    u_z \frac{\partial u_z}{\partial z}\right) &= -\frac{\partial p}{\partial z} + \mu \left(\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial u_z}{\partial r}\right) +
    \frac{1}{r^2}\frac{\partial^2 u_z}{\partial \phi^2} + \frac{\partial^2 u_z}{\partial z^2}\right),\\
    \frac{1}{r}\frac{\partial}{\partial r}\left(r u_r\right) +
    \frac{1}{r}\frac{\partial \left(u_\phi\right)}{\partial \phi} + \frac{\partial \left(
    u_z\right)}{\partial z} &= 0.

Additionally the following is assumed

  1. Steady flow, meaning \frac{\partial}{\partial t}(\cdot) = 0
  2. Radial and swirl components of the fluid velocity are zero meaning u_r = u_\phi = 0
  3. The flow is assumed to be axisymmetric, meaning \frac{\partial}{\partial \phi}(\cdot) = 0
  4. The flow is fully developed, meaning \frac{\partial}{\partial z} u_z = 0

First with this assumptions the continuity equation (fourth line of (65)) is trivially fulfilled. Further it follows that

(66)\frac{\partial}{\partial r} p &= 0,\\
\frac{1}{r}\frac{\partial}{\partial \phi} p &= 0,\\
\mu\left(\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial}{\partial r} u_z\right)\right) -
\frac{\partial}{\partial z}p &= 0.

From this we can deduce that p = p(z). Plugging this into the last line of (66) we can solve for u_z by twice integrating and get

u_z(r,z) = \frac{r^2}{4 \mu} \frac{dp}{dz} + c_1 \ln r + c_2

with some arbitrary integration constants c_1,c_2. To get a closed representation we will use the following observations

  1. On the boundary of the cylinder (r=R) we want u_z=0
  2. For r=0 we want u_z to be finite

The first obervation yields c_2 = -\frac{R^2}{4\mu} and the second yields that c_1=0. Putting all together we get a closed expression for u_z

u_z(r,z) = -\frac{1}{4\mu}\left(R^2 - r^2\right) \frac{dp}{dz}.

As a last assumption, let $p$ be decreasing linearly from z=0 to z=L yielding \frac{dp}{dz} =
\frac{\Delta p}{L}. Now we can calculate the volumetric flux Q through the cross section of the cylinder as

Q &= -2\pi \int_{0}^{R} u_z(r) r dr\\
  &= \frac{\Delta p \pi R^4}{8\mu L}

which gives the Law of Hagen-Poiseuille.

Numeric Verification

For the experiments we choose a cylindrical pipe with R=0.05m,~L=0.5m. We will vary the pressure drop \Delta p as well as the viscosity \mu. For simplicity we assume that \rho=1.0. We have three different discretizations (coarse, medium, fine). The coarse discretization has a maximal edge length of h_\text{coarse}=0.2R, the medium has h_\text{medium}=0.5h_\text{coarse}, and the fine has h_\text{fine}=0.25h_\text{coarse}.

  \Delta p [Pa] \mu [Pa s] \mathrm{Re}
Experiment 1 1.0 0.01 0.625
Experiment 2 5.0 0.01 3.125
Experiment 3 10.0 0.01 6.25
Experiment 4 15.0 0.01 9.38
Experiment 5 1.0 0.001 62.5
Experiment 6 5.0 0.001 312.5
Experiment 7 10.0 0.001 625.0
Experiment 8 15.0 0.001 937.5

The results for the coarse discretization are depicted in Table Tab. 22, the results for the medium discretization are depicted in Table Tab. 23, and the results for the fine discretization are depcited in Table Tab. 24. One can clearly see the influence of the Reynolds number on the accuracy of the numeric solution. For \mathrm{Re}\approx 1000 in a cylindrical pipe the assumptions for the law of Hagen-Poiseuille are not longer valid.

Usage

To run your own validation experiment just type in the command

./run.py --discretization TYPE --mu VALM --rho VALR --pressureDrop VALP --np NP

Here TYPE can be either coarse, medium or fine. This relates to the meshes used. The value VALM is the value for the dynamic viscosity \mu of the fluid in Pascal seconds. Value VALR is the value for the fluid density \rho in kilogram per cubic meters. VALP denotes the value for the pressure drop \Delta p given in Pascal. At the inflow a value of \frac{3 \Delta p}{2} will be prescribed and at the outflow a value of \frac{\Delta p}{2}. Last, NP stands for the number of processors.

The above code uses the standard meshes provided with this run script. This cylindrical mesh has the default length L=0.5[m], and default radius R=5[cm]. You can also generate your own mesh. This is achieved via the flags:

./run.py --discretization TYPE --mu VALM --rho VALR --pressureDrop VALP --np NP --generate 1 --radius R --length L
--radiusfactor HR

Here, R denotes the radius of your cylinder in centimeter, L denotes the length of your cylinder in meter and HR gives the factor for calculating the average edgelength in the mesh. It should be between zero and one. The edgelength is calculated as h_e=HR*0.01*R. To use this functionalty you need to have meshtool in your search path as well as gmsh. The software package gmsh can be downloaded from here. Depending on the chosen parameters the generation of the meshes can take some time.

Coarse

  Q_\text{numeric} Q_\text{HP} \lvert\frac{Q_\text{numeric}-Q_\text{HP}}{Q_\text{HP}}\rvert
coarse (1) 4.8157e-04 4.9087e-04 0.0190
coarse (2) 2.4077e-03 2.4544e-03 0.0190
coarse (3) 4.8157e-03 4.9087e-03 0.0190
coarse (4) 7.2243e-03 7.3631e-03 0.0188
coarse (5) 4.5198e-03 4.9087e-03 0.0792
coarse (6) 2.0044e-02 2.4544e-02 0.1833
coarse (7) 3.6429e-02 4.9087e-02 0.2579
coarse (8) 5.0822e-02 7.3631e-02 0.3098

Medium

  Q_\text{numeric} Q_\text{HP} \lvert\frac{Q_\text{numeric}-Q_\text{HP}}{Q_\text{HP}}\rvert
medium (1) 4.8869e-04 4.9087e-04 0.0044
medium (2) 2.4435e-03 2.4544e-03 0.0044
medium (3) 4.8875e-03 4.9087e-03 0.0043
medium (4) 7.3321e-03 7.3631e-03 0.0043
medium (5) 4.8394e-03 4.9087e-03 0.0141
medium (6) 2.3317e-02 2.4544e-02 0.05
medium (7) 4.5841e-02 4.9087e-02 0.0661
medium (8) 6.7741e-02 7.3631e-02 0.08

Fine

  Q_\text{numeric} Q_\text{HP} \lvert\frac{Q_\text{numeric}-Q_\text{HP}}{Q_\text{HP}}\rvert
fine (1) 4.9044e-04 4.9087e-04 0.0009
fine (2) 2.4523e-03 2.4544e-03 0.0009
fine (3) 4.9046e-03 4.9087e-03 0.0008
fine (4) 7.3573e-03 7.3631e-03 0.0008
fine (5) 4.9081e-03 4.9087e-03 0.0001
fine (6) 2.3535e-02 2.4544e-02 0.0411
fine (7) 4.6091e-02 4.9087e-02 0.0610
fine (8) 6.9728e-02 7.3631e-02 0.0530