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

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

Steady flow, meaning $\frac{\partial}{\partial t}(\cdot) = 0$
Radial and swirl components of the fluid velocity are zero meaning $u_r = u_\phi = 0$
The flow is assumed to be axisymmetric, meaning $\frac{\partial}{\partial \phi}(\cdot) = 0$
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

On the boundary of the cylinder ( $r=R$ ) we want $u_z=0$
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