 Calculation Tutorials   Turbulent Pipe Flow

# The Turbulent Pipe Flow Example - High-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-1.7 and other open software so that the reader can produce comparable results. I made these calculations on a Linux operating system (Ubuntu 10.04) using OpenFOAM-1.7, 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 1.7.1
• Octave
• Grace or xmgrace (I am using Grace-5.1.22)
• gnuplot

## Turbulent pipe cases

The case will be calculated with a High-Re turbulence model. An other turbulent pipe flow tutorial calculates the same example with a Low-Re model. One may find information on the two calculation types here and here.
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

## Meshing

A similar mesh can be used as in the laminar case. The only difference is that we will use a higher resolution in radial direction (10 cells). First copy the full laminar case:

cp -rv <laminar_case_directory> <turbulent_case_directory>

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

Run the meshing again:

blockMesh

The mesher will create a default patch with no elements in it. Delete this patch and change the number of patches to 5 in the constant/polyMesh/boundary file: ## Set up the case

### Boundary & Initial Conditions

The boundary and initial conditions can be set within the field files in the directory of the first iteration/timestep (counting time/iteration starts with 0). This time we will need the turbulent field variables as well, so they should be copied new from the tutorials if they have been deleted in the laminar flow tutorial. Copy the k field first:

cp -rv \$FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily/0/k <turbulent_case_directory>/0/

Repeat this for the epsilon and nut fields also.
The boundary fields should be edited to match with the patch names defined in the boundary file in <turbulent_case_directory>/constant/polyMesh.

When using a High-Re model a wall function should be used. The main idea is that the boundary layer is not resolved by the mesh (we do not have to use a really fine mesh near the wall), but the boundary layer is calculated by a semi-empirical formula in one cell: the boundary layer is in the first cell. The type of wall function to be used can be defined in the nut field file. First We will use the so called standard wall function which is called nutWallFunction in OpenFOAM®.

Additionally we have to tell the solver that we use wall functions for k and epsilon.
See the k,epsilon and nut files:   Initialize the velocity field to the inlet value (according to the shown turbulence fields).

This is similar case to the simpleFoam pitzDaily tutorials where wall functions are used also.
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: is given I=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=1.5e-4 m2/s2
• Turbulent Dissipation Rate: ε=Cμ3/4*k3/2/lsc=2.1562e-5 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 1.5e-4; fixedValue; / value uniform 2.1562e-5;
wall wall kqRWallFunction; / value uniform 1.5e-4; epsilonWallFunction; / value uniform 2.1562e-5;
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 nutWallFunction; / 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): Set the followings in the file (keep the fluid type Newtonian): ### Set turbulence model (RASProperties)

We will use the k-epsilon turbulence model (called kEpsilon in OpenFOAM®, see available turbulence models here) therefore the RASProperties file should be modified:  ### Solver settings (fvSolution) and numerical schemes (fvSchemes)

fvSolution:

• Lets change the linear solver type for the p equation to AMG. This will speed up the calculation.
• Decrease the tolerances to obtain better convergence for all the used solvers: fvSchemes: keep it unchanged.

### Setting solver running options (controlDict)

Change only the endTime and writeInterval in the controlDict file (both to 250):  ## 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

A similar gnuplot script can be used as in the laminar flow example, but it should be changed to monitor the k and epsilon residuals as well: Residual history for this calculation: ## 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): We can not be sure if that axial velocity is converged until the end of the pipe. Therefore launch the case on a longer mesh. Lets change the mesh length to 16m. Follow these steps:

• change all occurences of the 8m coordinate to 16m in the constant/polyMesh/blockMeshDict file
• re-generate the mesh
user@machine:~\$ blockMesh
• edit the constant/polyMesh/boundary file and delete the default patch. Do not forget to change the number of patches to 6. This will result in a coarser mesh.
• re-run the case:
user@machine:~\$ simpleFoam > log &
• change the sampling coordinates in system/sampleDict1 and sampleDict2. Run the sample utility again using both sampleDict files:
user@machine:~\$ rm -fv system/sampleDict
user@machine:~\$ ln -s system/sampleDict1 system/sampleDict
user@machine:~\$ sample -latestTime
user@machine:~\$ rm -fv system/sampleDict
user@machine:~\$ ln -s system/sampleDict1 system/sampleDict
user@machine:~\$ sample -latestTime

These sampling steps are included in a shell script. Se it in the case package at the downloads section.

This is the new axial velocity plot: We can see that the axial velocity finally converged to the value of 1.2 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: Modify the sampleDict2 file from the laminar example files to have also the yPlus sampled:

: This is finally the yPlus plot: We can see that y+ is too low. Therefore we should recalculate the case with a lower radial mesh resolution. Follow these steps:

• copy the case to have the curent results available for comparison
• edit the blockMeshDict file and change the z resolution to 6 (hex (0 1 1 0 2 3 4 5) (200 1 6) simpleGrading (1 1 1))
• (move to the newly created directory)
• create the mesh new:
user@machine:~\$ blockMesh
• start the solver (previously copied results will be overwritten!):
user@machine:~\$ simpleFoam > log &
• create the yPlus and wallShearStress fields new:
user@machine:~\$ yPlusRAS -latestTime
user@machine:~\$ wallShearStress -latestTime
• run the sample utility again using both sampleDict files y+ is closer to 30 with the rough mesh!

Create the skin friction coefficient now from the wallShearStress plot:

• load the wallShearStress_wall.raw file as block data. Load column 1 as x and column 4 as y
• multiply the curve with -2 to have absolute values and to calculate the skin friction coefficient from wallShearStress (ρ=1 kg/m3, we have to multiply by 2 to have the division by 1/2): Data > Transformations > Geometric transformations > Select the curve on the Apply to set panel > set Scale Y as -2 > Accept

This is what you should get: The value at the developed flow range (end of pipe) is ~0.00734 for the both meshes. The corresponding Fluent result is close to 0.01 (according to the graph shown).

The radial profile is already created together with the axial profile when using the sampleDict1 file with the sample utility. As sample created the output the default formating is: velocity on the y axis and distance on the x axis. We can change this to have an easier comparison with the Fluent results:

• Data > Transformations > Geometric transformations > Select the curve on the Apply to set panel > set Rotation (degrees) as 90 and Scale X as -1 > Accept Comparing it with the Fluent radial profile we can see that the max velocity is better captured on the rough mesh (200x6). Please note that the Fluent calcultion was done with a Low-Re model, therefore the distribution is finer due to the higher number of points displayed. 