Hagen-Poiseuille flow

Here, I am going to go over the solution to fully developed laminar pipe flow. This is a canonical problem and provides an exact solution to the Navier-Stokes equations. It is often referred to as the Hagen-Poiseuille flow problem.

The problem is depicted in the figure below:

Fig: Hagen-Poiseuille Flow Problem Definition

In this problem, we are examining laminar flow through a pipe. The problem states that the flow in the pipe is being driven by a constant pressure gradient in the axial direction (dP/dz=constant). We assume that the flow is purely axial (v_r = v_\theta =0), steady state (\partial / \partial t  =0), incompressible (ρ=constant), axisymmetric (\partial / \partial \theta  =0). We also neglect gravity.

First we start with the axial Navier-Stokes momentum equation in cylindrical coordinates:

\rho v_z \frac{\partial v_z}{\partial z} = -\frac{dp}{dz}+\frac{\mu}{r} \frac{d}{dr}\left(r \frac{dv_z}{dr}\right)

By using our assumptions we can reduce this equation to give us:

\frac{dp}{dz}=\frac{\mu}{r} \frac{d}{dr}\left(r \frac{dv_z}{dr}\right)

Because we know that dP/dz is a constant, this function is easily integrated twice. The first integration:


And the second integration:

v_z=\frac{r^2}{4\mu}\frac{dP}{dz}+A \ln(r)+B

Now we apply our boundary conditions!

No slip condition at r=R : v_z=0=\frac{dP}{dz}\frac{R^2}{4\mu}+A\ln(r)+B

Finite axial velocity at r=0: v_z=0+A \ln(0)+B

Since ln(0) is a discontinuity, we know that in order for these equations to be satisfied, A must be equal to zero (ie A=0). Then combining these two equations, we get that:


Thus our final solution is:


And there we have our answer!

Couette Flow (Moving Wall)

In this post I am going to go through the solution to the moving wall Couette flow problem.

An illustration of the problem is given below:

Fig 1: Illustration of Couette Flow

In this problem, the fluid between two parallel plates is being driven by the motion of the top plate. Here, we assume that the flow is axial (v=w=0), incompressible (ρ=constant), fully developed (u only a function of y). We also assume that gravity can be neglected and that no pressure gradient is present (dp/dx=0).

First we start with the Navier-Stokes momentum equation in the x-direction:

\rho\left(u\frac{\partial u}{\partial x}+ v \frac{\partial u}{\partial y}\right) = \frac{\partial p}{\partial x} + \mu\left( \frac{\partial^2 u}{\partial y^2}+\frac{\partial^2 u}{\partial y^2}\right)

By applying our assumptions listed above, you should be able to see that the equation simply becomes:

\frac{d^2 u}{dy^2}=0

To solve the problem we must integrate this equation and solve using the boundary conditions defined by the problem. The integration results in:


Our boundary conditions come from the problem and our super smart knowledge of the no-slip condition ;). Ie.

@ y=h, u=V

@ y=0, u=0

After subbing in we get two equations, with two unknowns (the integration constants)



Therefore we can see that B=0 and A=V/h. This leads to the final solution of this plane Couette flow problem:


Now we have shown that the velocity profile in this case is the linear profile above. We can also calculate the shear stress:

\tau=\mu \frac{du}{dy} = \mu \frac{V}{h}

Any thoughts or questions please respond in the comments!



Definition of the Speed of Sound

Deriving the speed of sound is a main component of any undergraduate textbook in gas dynamics and is very important to understanding more complicated concepts! Here, I will give an honest attempt to explain and provide definitions of what it is and why it is important .I will go through the derivation of the speed of sound mathematically.

First, I’ll attempt to explain what the speed of sound is. Most simply put, it is the speed at which sound propagates in a medium. But in fact, it is not just sound. All disturbances and motions are communicated throughout a medium at the speed of sound. For example, ever notice how streamlines around a vehicle begin to bend far upstream from the vehicle itself? This occurs because the presence of the vehicle has been communicated upstream (at the speed of sound) and the gas is forced to move out of the way. The speed of sound is dependent on the properties of the fluid being described, and the fluid itself. For instance the speed of sound in air at standard conditions is ~343 m/s whereas the speed of sound in water is ~1482 m/s. Additionally, and very importantly, the speed of sound represents the coupling of the pressure and density fields in a gas. In an incompressible fluid, the density is constant and therefore not a function of pressure. So we can call these fields uncoupled. But in real flows this is not the case. In fact, it is this couple of the pressure and density fields that make important and interesting changes to the conservation of momentum equation! (More on that in another post)

So why is the speed of sound important? Well let’s think about the vehicle example again. The streamlines upstream are able to bend in response to the oncoming vehicle. This is because the vehicle is travelling less than the speed of sound. But what if it were travelling faster than the speed of sound? Now the streamlines do not see what is coming and are suddenly struck by the car. This results in a shock-wave, and a very different aerodynamic scenario. In my related post regarding (subsonic vs. supersonic flow) I will re-enforce this discussion mathematically.

Now let’s make use of an example to define the speed of sound mathematically. The most straight forward approach to defining the speed of sound is with a piston example:

Fig: Schematic of the piston example

In the scenario depicted above, a piston is suddenly accelerated to a velocity of dV. This then causes a pressure increase of dp, and sends a compression wave (at the speed of sound) down the channel. The gas in the channel that has not yet been reached by the wave sits at the initial conditions of zero velocity, and initial pressure and density p_o and \rho_o respectively. As the wave moves down stream it causes the velocity of the gas to increase to dV and the pressure and density to increase by dp and dρ, respectively.

Remember: In an incompressible fluid, a piston motion would instantaneously cause all of the fluid in the tube to change velocity to dV. But recall that in a compressible fluid the motion of the piston does NOT do this. It sends a message downstream at the speed of sound telling the gas that the piston has moved. The wave then causes the fluid downstream to respond so as to ensure that the governing laws of fluid mechanics are held.

We are going to use this scenario and solve for the speed of sound, a. In order to analyze this problem we are going to perform a control volume analysis of the compression wave itself. To accomplish this, we are using a control volume that is in the frame of reference of the wave (ie moving to the right at velocity a). This is shown in the following figure:

Fig: Control Volume Schematic

Notice that when we put the control volume in the reference frame of the wave how the problem is simplified? We now have a simple control volume with an inflow and outflow and we need only apply conservation of mass and momentum to solve the problem.

First we start with conservation of mass (Note: A-> cross sectional area of channel):


\left(\rho_o+d\rho\right)\left(a-dV\right)A-\rho_o a A = 0

If we expand, and simplify by removing the high order terms we get:

a d\rho-\rho_o dV=0

Now we can add the momentum equation:

\Sigma F = \dot{m}V_{out}-\dot{m}V_{in}

p_oA-\left(p_o+dp\right)A = \rho_o a A \left[\left(a-dV\right)-a\right]

We can simplify this to:

dp=\rho_oa dV

This equation, and the simplified conservation of mass are then easily combined to achieve:


This is the definition of the speed of sound! However, we must recognize that a sound wave is an isentropic process. And this is not always the case. So the correct definition of the speed of sound uses partial derivatives and a subscript s to denote that it is for an isentropic wave.

a^2=\left(\frac{\partial p}{\partial \rho}\right)_s

So then how do we calculate the speed of sound for a medium! Well, this comes from simply recalling that in an isentropic process the relation between pressure and density is:

\frac{p}{\rho^\gamma} = \textnormal{constant}

From this we can show that:

\frac{\partial p}{\partial \rho} = \frac{\gamma p}{\rho}

And then finally by using the ideal gas law we end up at

a= \sqrt{\gamma R T}

The above relation is the most commonly used formula for the speed of sound.


Thanks for reading! Hopefully this provided some clarification and assistance to those who needed it! As always, please comment if you notice any mistakes, or even little semantic errors. It is important to me to be completely correct!

Some Useful References

… Any Undergraduate Gas Dynamics textbook but my preferred ones are:

[1] Anderson, J. D. (1990). Modern compressible flow: with historical perspective (Vol. 12). McGraw-Hill.

[2] John, J. E. A., & Keith, T, G. (2006) Gas Dynamics, 3rd Edition, Pearson Prentice Hall

[4] Zucrow, M. J., & Hoffman, J. D. (1976). Gas dynamics. New York: Wiley, 1976, 1.

Turbulent Zero Pressure Gradient Flat Plate – simpleFOAM, komegaSST

An excellent test case, and case to familiarize yourself with some of the turbulence models available in OpenFOAM is a 2D flat plate with zero pressure gradient. I will solve this problem using the solver simpleFoam and the komegaSST model.

Here are the sections of this post:

  1. Quick Overview: kω-SST (komegaSST) Boundary Conditions
  2. Case set-up and mesh
  3. Results
    • u+ vs. y+
    • Coefficient of friction
  4. Conclusions and useful resources

Download the case file here: kOmega Flat Plate Tutorial File

Quick Overview: kω-SST (komegaSST) Boundary Conditions

In this section I will describe the boundary set-up for komegaSST  where no wall functions are implemented. This requires that the y+ along the wall is less than or equal to one. For the komegaSST turbulence model the boundary conditions are as follows:

At the wall:

  • ω (omega) – specific dissipation rate
    • BC type: fixedValue
    • BC value: \omega_{wall}=10\frac{6\nu_\infty}{\beta_1\left(\Delta y_{wall}\right)^2}
  • k – turbulent kinetic energy
    • BC type: fixedValue
    • BC value: 0
  • nut – turbulent viscosity
    • BC type: fixedValue
    • BC value: 0

In the free-stream:

  • ω (omega) – specific dissipation rate
    • BC type: fixedValue
    • BC value: \frac{U_\infty}{L} < \omega_{\infty} < \frac{10 U_\infty}{L}
  • k – turbulent kinetic energy
    • BC type: fixedValue
    • BC value:\frac{10^{-5}U_\infty^2}{Re_L} < k_{\infty} < \frac{0.1U_\infty^2}{Re_L}
  • nut – turbulent viscosity
    • BC type: calculated
    • BC value: 0 (this is just an initial value)

where \beta_1=0.075, and \Delta y_{wall} is the normal distance from the wall to the first cell center.

Case set-up and mesh

Free-Stream Properties

For this case I have followed a similar set up to the 2D flat plate case used on the NASA turbulence modelling resource website. By doing this it gives me something to compare to! However, I am not going to use the exact set-up on that website. Since we are using simpleFoam, I am going to set up this case the way I prefer which is to use a velocity (U) of 1 m/s and scale all other properties accordingly. The simulation properties that I used are :

  • U∞=1 m/s
  • ν=4e-7 m2/s
  • L=2 m

These correspond to a Reynolds number at L=2m of 5 million .

Grid Generation

The grid used was generated in blockMesh. High inflation was used in the boundary layer region in order to achieve the desired y+ value of less than 1. For more details on grid generation using blockMesh see the OpenFOAM manual!


Boundary Conditions

For the incompressible solver simpleFoam, the minimum boundary conditions required for a simulation are p and U. However, if the simulation is a RANS simulation additional boundary conditions are required. For the kω-SST model we need to have a boundary condition on k and ω as well. The boundary conditions I defined in the zero (0) folder can be found in the attached tutorial file.

The only boundary condition that really needs any comment is omega. We calculate omega using \omega_{wall}=10\frac{6\nu_\infty}{\beta_1\left(\Delta y_{wall}\right)^2}. In our case the wall distance to the first cell center is 5.48316E-06. Using our free stream viscosity of 4E-7 this gives a value of omega at the wall of 10643595.39.

Tip for fvSolution

If you find that the results you are getting are wrong, it could be that the residuals for the different properties are too high! Certain properties converge before others and therefore you need to ensure that they all converge to a sufficiently low value!


First we compare the coefficient of friction to the .dat file available from the .DAT available from the NASA Turbulence modelling resource. NOTE: We simulated at Reynolds number of 5 million whereas the NASA setup is at 10 million. So the x coordinate in the following is plots is rescaled and in fact we are only using half of the data from the .DAT file!

Coefficient of Friction

We can see from the figure that the coefficient of friction from our simulation matches the expected data closely! Hurray!

Next let’s compare the u+ vs. y+ profile to the  universal profile for turbulent boundary layers. Recall that u+ is the velocity normalized by the friction velocity (u_*), and y+ is given by the following equation:

y^+ = \frac{u_* y}{\nu}

The u+ is vs y+ is plotted here:

u+ vs y+

We can see from the figure that our solution is pretty good! The y+ value of the first node is located around a y+ of approximately 0.5, the viscous sublayer matches very closely and the log law layer is not significantly off!


In this post we simulated a zero pressure gradient flat plate at a Reynolds number of 5 million. We compared the results for shear stress to the NASA turbulence modelling resource expected results and showed close agreement. Then the u+ vs y+ profile was compared to the universal law of the wall and again the results were okay!

Some useful references

The NASA Turbulence Modelling Resource, http://turbmodels.larc.nasa.gov/



Alternate Set-up: Oblique Shock in rhoCentralFOAM

I have already done a post on supersonic flow over 23 degree wedge. In that post, we were intentionally creating a situation where the shock wave would be detached from the obstacle. The set-up in that case is certainly appropriate when the shock is attached as well. However it is possible to use an even simpler set up for an attached oblique shock.

Continue reading Alternate Set-up: Oblique Shock in rhoCentralFOAM

Numerical Solution of the Compressible Laminar Boundary Layer Equations

In this post I go over the numerical solution to the compressible boundary layer equations. This is very useful when a quick estimate of shear stress, wall heat flux, or boundary layer height if necessary.

Continue reading Numerical Solution of the Compressible Laminar Boundary Layer Equations