Calculation Tutorials 
 Turbulent Pipe Flow 

The Turbulent Pipe Flow Example - Low-Re

Foreword

The purpose of the Cornell tutorials with OpenFOAM® is to give a few easy examples which are well known for the CFD users. The original tutorials were written by Rajesh Bhaskaran. They give the user a practical guide on software usage (FLUENT / ANSYS) and the bacground is also explained well. Therefore it is recommended to read the particular Cornell tutorial first. In my tutorials I show the same steps using OpenFOAM® and other open software so that the reader can produce comparable results. I made these calculations on a Linux operating system (Ubuntu 11.04) using OpenFOAM-2.0.0, Octave, gnuplot and xmgrace.
The tutorials are built for reading one after the other, therefore once something is explained it is expected as known afterwards.

The original tutorial can be found here.

Software Prequisites:

  • OpenFoam® 2.0.0
  • Octave
  • Grace or xmgrace (I am using Grace-5.1.22)
  • gnuplot

Turbulent pipe cases

The case will be calculated with a Low-Re turbulence model.
The dimensions of the domain does not change compared to the laminar pipe flow example. The case will be turbulent because the viscosity is decreased (2e-5 instead of 2e-3) resulting in a higher Re number:
Re = v*D/ν = 1*0.2/2e-5 = 10000
I will modify the high-Re case so one can start with copiing it under a new case directory.

Meshing

A similar mesh can be used as the one in the turbulent High-Re case finally. The mesh will be 8 m long again as in the laminar case. According to the Cornell tutorial we need high resolution in the radial direction. The Cornell tutorial use 30 cells in radial direction and a bias of 10 is applied. This means that ratio of cell length of the first and last cell in the radial direction is 10 (towards the side of the pipe cells are getting smaller).
Let's copy the full case as we used in the High-Re case:

cp -rv <turbulent_High-Re_case_directory> <turbulent_Low-Re_case_directory>

Change the blockMeshDict file:
hex (0 1 1 0 2 3 4 5) (100 1 30) simpleGrading (1 1 0.1))

Run the meshing again:

blockMesh

As in the High-Re tutorial do not forget to remove the default patch.

Set up the case

Boundary & Initial Conditions

The boundary and initial conditions will be the same for the flow field variables (velocity and pressure), but the turbulence boundary conditions need some changes.

A usual set up would use a small fixed value for k at the walls (small BUT NOT 0, e.g. 1e-10) and zeroGradient for epsilon . For nut we will use the nutUSpaldingWallFunction wall function. According to the description in its header file: "Wall function boundary condition for walls, based on velocity, using Spaldings law to give a continuous nut profile to the wall (y+ = 0)". The path of the header file is: /opt/openfoam200/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions/nutUSpaldingWallFunction/nutUSpaldingWallFunctionFvPatchScalarField.H. This wall function therefore can be used at low y+ values. The usage of nutLowReWallFunction should give similar results as this is the standard low-Re wall function.

See the changes in the k,epsilon and nut files:



Initialize the velocity field to the inlet value (according to the shown turbulence fields).

Turbulence boundary conditions at the inlet can be specified using the definitions in the Fluent Help:

  • Inlet velocity: V=1 m/s
  • Reynolds number: Re=V*D/ν=10000
  • Turbulence intensity: I=0.16*Re^(-1/8)=0.050596 - I use the assumption instead of the given value of 0.01.
  • Turbulent length scale: is ~7% of the diameter, lsc=0.07*D=0.014 m
  • Turbulent Kinetic Energy: k=3/2*(V*I)2=3.84e-3 m2/s2
  • Turbulent Dissipation Rate: ε=Cμ3/4*k3/2/lsc=2.79e-3 m2/s3, where Cμ=0.09.

These values should be put to the inlet value and to the value at the wall functions also.

Summary of all the boundary conditions:

Boundary OF type Velocity - derived type / value Pressure - derived type / value
inlet patch fixedValue; / value uniform (1 0 0); zeroGradient; / -
outlet patch zeroGradient; / - fixedValue; / value uniform 0;
wall wall fixedValue; / value uniform (0 0 0); zeroGradient; / -
axi_symm-f wedge wedge; / - wedge; / -
axi_symm-r wedge wedge; / - wedge; / -

Boundary OF type k - derived type / value epsilon - derived type / value
inlet patch fixedValue; / value uniform 3.84e-3; fixedValue; / value uniform 2.79e-3;
outlet patch zeroGradient; / - zeroGradient; / -
wall wall fixedValue; / value uniform 1e-10; zeroGradient; / -
axi_symm-f wedge wedge; / - wedge; / -
axi_symm-r wedge wedge; / - wedge; / -

Boundary OF type nut - derived type / value
inlet patch calculated; / value uniform 0;
outlet patch calculated; / value uniform 0;
wall wall nutUSpaldingWallFunction; / value uniform 0;
axi_symm-f wedge wedge; / -
axi_symm-r wedge wedge; / -

Setting viscosity (Transport properties - transportProperties)

We shall change now the value of viscosity in the transportProperties file (located in the constant directory):

Highslide JS

Set the followings in the file (keep the fluid type Newtonian):

Set turbulence model (RASProperties)

We will use a low-Re k-epsilon turbulence model: LaunderSharmaKE (see available turbulence models here > RAS turbulence models for incompressible fluids). Modify the turbulence model in the RASProperties file:

Highslide JS

Launder Sharma description in its header file: "Launder and Sharma low-Reynolds k-epsilon turbulence model for incompressible flows.". Header file path: /opt/openfoam200/src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.H.

Solver settings (fvSolution) and numerical schemes (fvSchemes)

fvSolution: we keep the file as it was used in the high-Re case. The changes are:

  • Solver for p equation is AMG.
  • Decreased solver tolerances:

fvSchemes: according to a change in the solver code in OpenFoam® 2.0.0 one term should be changed in the fvSchemes file:

Setting solver running options (controlDict)

Change only the endTime and writeInterval in the controlDict file (both to 1500):

Highslide JS

Calculation

Launch the solver

Launch the solver the same way as in the laminar flow example:

user@machine:~$ simpleFoam > log &

Monitor the log file:

user@machine:~$ tail -200f log

Monitor residuals

Use the same gnuplot script as in the High-Re example.

Residual history for this calculation:

Highslide JS

Download the Residuals-turb gnuplot script.

Post Processing

Axial velocity profile

Start with checking the axial velocity profile. Exporting the profile using the sample utility (use the sampleDict1 file from the laminar flow example files) one will get the following (as datasets were written to the format of Grace, it is enough to double-click on the axial_profile_U.agr in the file manager e.g. in nautilus when using Ubuntu):

Highslide JS

We can see that the axial velocity finally converged to the value of 1.25 m/s. The comparable Fluent result can be checked here. The converged axial velocity obtained with ANSYS 12 is 1.19 m/s (according to the chart shown).

Check the Skin friction coefficient

We should check the height of the near wall cells now. The non-dimensional measure for this is the y+ (definition @ cfd-online). Using a High-Re turbulence model y+ should be ~=30. The yPlus field is not written by default, we should generate it using the yPlusRAS utility:

user@machine:~$ yPlusRAS -latestTime

We will need the wallShearStress for checking the skin friction coefficient (definition @ cfd-online), therefore create it now:

user@machine:~$ wallShearStress -latestTime

One may find the definition of skin friction coefficient

The following new files are created:

Highslide JS

Load the yPlus plot the same way as in the High-Re tutorial:

Highslide JS

We can see that y+ is in the appropriate range.

Create the skin friction coefficient now from the wallShearStress plot the same way as in the High-Re tutorial. This is what you should get:

Highslide JS

The value at the developed flow range (end of pipe) is ~0.00867. The corresponding Fluent result is close to 0.01 (according to the graph shown).

Radial velocity profile

The radial profile is already created together with the axial profile when using the sampleDict1 file with the sample utility. Create the radial velocity profile the same way as in the High-Re tutorial:

Highslide JS

Comparing it with the Fluent radial profile we can see that the max velocity is predicted to be higher with OpenFOAM®. In addition the profile is not as flat in the center as in the Fluent result.

Non-dimensional velocity profile

Let's check the velocity profile in a non dimensional form which gives a good posibility of validation against measurements. For reference please read through this wikipedia article.
We shall calculate the non dimensional velocity and distance from the wall with the following formulae:





Start to create the new plot from the existing radial velocity plot (shall be placed here: /sets/1500/radial_profile_U.agr if created with sampleDict1). It is the plot already posted before. Apply the following modifications to the plot:

  1. move the data with -0.1 among the x axis
    Data > Transformations > Geometric transforms > Translate X: -0.1
  2. mirror the data to the y axis
    Geometric transforms > Scale X: -1.0, in the meanwile change "Translate X" back to 0!
  3. apply transformation on the y values (note: y values are the velocity values - u+!)
    Data > Transformations > Evaluate expression
    select "Set" and "Source" to the available data set. Type to the formula window: "y=y/0.065849" > Apply
  4. apply transformation on the x values (note: x values are the distances!)
    Type to the formula window: "x=x*3292.4" > Apply
  5. set x axis as logarithmic
    Plot > Axis properties > Edit: X axis, Scale: Logarithmic > Apply
  6. create reference data set for the viscous sublayer
    Edit > Data sets > Edit > Create new > By formula > Start at: 0.1, Stop at: 1000, Length: 4000, X=$t, Y=$t > Accept
  7. create reference data set for the log-law layer
    Edit > Data sets > Edit > Create new > By formula > Start at: 0.1, Stop at: 1000, Length: 4000, X=$t, Y=1/0.41*ln($t)+5 (note: constants κ=0.41 and C+=5.0 are applied to this equation!)
  8. do some formatting on the data sets
    double-click on one of the data graps - Set Appearance pops up > apply your settings

You should get something similar:

Highslide JS

Download the case directory

turbulentPipeLowRe.tar.gz


Version: 1.1, 2012-01-28

  1. TOP
  2. Foreword
  3. Software Prequisites
  4. Turbulent pipe cases
  5. Meshing
  6. Set up the case
    1. Boundary & Initial Conditions
    2. Setting viscosity
    3. Set turbulence model
    4. Solver settings, numerical schemes
    5. Solver running options
  7. Calculation
    1. Launch the solver
    2. Monitor residuals
  8. Post Processing
    1. Axial velocity profile
    2. Skin friction coefficient
    3. Radial velocity profile
    4. Non-dimensional velocity profile
  9. Downloads
  10. Navigation to other pages
navigation