The Laminar Pipe Flow Example
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) 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 original tutorial can be found here.
The laminar pipe flow example introduces a 2D axi-symmetric mesh which is really 2D in FLUENT. In OpenFOAM® the mesh for the same purpose is a one layer thick "wedge" mesh (3D).
I decided to use the simpleFoam solver therefore I took the pitzDaily incompressible example from the tutorials. The simpleFoam solver is an incompressible RANS turbulent code using the SIMPLE algorythm.
cp -rv $FOAM_TUTORIALS/incompressible/simpleFoam/pitzDaily <location_where_you_want_to_place_the_example>
We will need one wedge block with 6 vertices. The axis will be perpendicular with the x axis. To generate the mesh We will use the blockMesh meshing utility of OpenFOAM®. It requires the blockMeshDict file where the vertices, blocks and boundary interfaces are defined. The blockMeshDict is placed here within the case directory:
The geometry is a pipe with D=0.2m diameter and L=8m length, therefore the vertice definition is the following:
Note the followings:
I usually do the necessary calculations with Octave which is a Matlab clone. One can calculate the coordinates needed in the console after launching Octave:
Calculation of the coordinates:
According to the Cornell tutorial the mesh should have 5 elements in radial direction and 100 elements in axial direction. Definition of the block therefore looks like this:
In The edges section curved edges (spline, arc) can be defined.
First the edge type should be defined followed by the starting and end point numbers (these points should are defined at the "Vertices" section). These are followed by the coordinates of a third point (in case of an arc or spline). One of the available options is the polySpline where several points can be used which will be crossed. (See all options here.)
Finally the boundary interfaces should be defined with their corner points. Use a continous list (e.g. do not put points of diagonal corners after each other). In the patches section therefore the base patch type, the patch name and the block (defined by 4 vertices for each) should be listed:
Note that in case of the wedge boundary condition the two side should be in separate patches! (On the contrary both parts of the "2D boundary condition" - empty - can be in one patch.)
As the blockMeshDict is set the meshing can be launched with simply tiping:
in the console within the case directory. The ParaView postprocessor (with the OpenFOAM® reader) can be launched with:
This will load the case. First click on "Apply", than set "Surface With Edges" in the drop down menu.
The boundary patches can be viewed by activating them on the left side in the mesh parts window. Changes should be Applied again.
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):
We are calculating laminar flow therefore the turbulence related fields (epsilon, k, nut, nuTilda and R field files) can be deleted.
Summary of all the boundary conditions:
Setting viscosity (Transport properties - transportProperties)
Viscosity (with many others) can be set in the transportProperties file. This file is located in the constant directory:
Set the followings in the file:
simpleFoam uses the kinematic viscosity which can be obtained from the given dinamic viscosity (μ=2e-3 kg/m/s) and density (ρ=1 kg/m3): ν=μ/ρ. The fluid type stays Newtonian.
Disable turbulence (RASProperties)
We are about to calculate laminar flow, so the turbulence settings should be modified. The related file is the RASProperties:
In this file one of the available turbulence models should be specified. See the list of turbulence models here under the "RAS turbulence models for incompressible fluids" keyword. One other option is laminar, what we need here:
Review the solver settings (fvSolution)
The fvSolution solver settings file is located in the system directory:
The file contains the settings of linear solvers, number of non-orthogonal corrector loops (in our case the mesh is rather nice hence no corrector loops are required), and under-relaxation factors.
Not expecting problems during the calculation the under-relaxation factors can be also left untuched.
Review the settings of numerical schemes (fvSchemes)
The fvSchemes settings file is located in the system directory:
Descritpion of the available options can be found here. The basic set up of the pitzDaily example uses:
Available options can be anytime checked by entering something which is not availble for sure (e.g. put div(phi,U) Gauss a; instead of div(phi,U) Gauss upwind. The solver will list the valid options.)!
Setting solver running options (controlDict)
In OpenFOAM® time and data control is done via the controlDict file. It is placed in the system directory:
Launch the solver
As we have finished with all the pre-processing work we can start the solver now. Type in the following command to the console (in the case directory):
user@machine:~$ simpleFoam > log &
The continously written log file can be monitored with:
user@machine:~$ tail -200f log
This comand will display the last 200 lines of the file called "log". The "f" option forces the tail program to re-display the last lines whenever it is changeing. Using redirection and the tail Linux program we can still monitor the log file while it is beeing written.
An example piece of the solver output:
The residuals are written into the log file. These values can be monitored during the run with gnuplot using a nice script I have found on the cfd-online OpenFOAM forum here. The script extracts the numbers so it is readable for gnuplot. The script I present here is not working with all the solvers: if the name of a variables changes, or we add a corrector loop (and therefore a variable is calculated more than once per iteration) the script will fail. But of course it can be adapted. Moreover the sript will look for the file called "log".
Launch gnuplot with the script like this:
user@machine:~$ gnuplot Residuals-lam
Download the Residuals-lam gnuplot script.
Here is the output:
We can see, that the residuals are not decreasing anymore after iteration 700. So 1000 iteration would have been enough for this calculation.
Once the solver writes an output it creates a new directory for it into the case directory. In our case (due to the writeInterval setting of 1000) we have 3 new directories:
They contain the new field files. The internalField section (we used this the internalField for initialisation in the 0 field) contains now the newly calculated values. Result field U:
Check the result fields in ParaView
Results can be loaded to the previously used ParaView. The program would start using the command "paraview", but OpenFOAM® comes with a script (paraFoam) which loads the actual case where it is executed (as we see before when checking the mesh). To see all available options of paraFoam type:
user@machine:~$ paraFoam -help
To load the last results set (iteration 3000) type:
user@machine:~$ paraFoam -latestTime
Once the ParaView window appears click on "Apply" (as shown before). One may find a good description about the usage of ParaView at the OpenFOAM® documentation page. I will give here only a limited set of information, for details pleas consult the previous link.
Navigation is the viewer window:
The toolbars are shown here. One may change to the last iteration using the arrovs in the "VCR Controls" toolbar. A field (U or p) can be loaded in the drop-down menu of the "Active Variable Controls". The velocity distribution viewed from the inlet side of the pipe:
Velocity profiles can be visualized with the "Plot Over Line" filter from the "Filters" menu, but the profiles can be provided from OpenFOAM® directly.
Exporting xy data from OpenFOAM
We can export graphs - axial and radial velocity profiles - from OpenFOAM® using the sample utility (The full list of utilities). The set up is done in the sampleDict file which should be placed in the system directory. An example can be found in the distribution under the applications directory, which can be copied by:
user@machine:~$ cp $FOAM_APP/utilities/postProcessing/sampling/sample/sampleDict ./system/
The sampleDict modified for our needs:
user@machine:~$ sample -latestTime
The output will be placed in the newly created sets directory:
The new files can be opened with xmgrace (To change the appearance of a graph, just double click on it and set what you wish.):
user@machine:~$ xmgrace sets/3000/axial_profile_U.agr
user@machine:~$ xmgrace sets/3000/radial_profile_U.agr
One can see that the results are quite similar to the ones obtained from FLUENT according to the original Cornell tutorial (axial profile; radial profile > be careful, the axes are exchanged in the radial profile).
Skin Friction Coefficient
According to the definition of the skin friction coefficient we the only unknown for the evaluaton is the wall shear stress. This can be calculated from the saved fields with the wallShearStress utility.
user@machine:~$ wallShearStress -latestTime
will create a new field in the directory of the last iteration:
wallShearStress is a surface field (it can be visualized in ParaView only when the wall patch is enabled). We need its values in a graph. This can get the walues of this field again with the sample utility. The modified sampleDict:
I usually create links for the sampleDict file:
user@machine:~$ ln -s system/sampleDict2 system/sampleDict
Naming the two sampleDict files "sampleDict1" and "sampleDict2" sampling can be done one after the other by changing the link (once it points to "sampleDict1" than to "sampleDict2").
The new output is a surface output it will be placed to the newly creates surfaces directory. The file is written in raw mode, so it contains coordinate components and wallShearStress components. Removing the header lines (let's call the new file "wallShearStress.xyzv") the file can be loaded in to Octave and the required distribution can be calculated from by calculating the magnitude of the wall shear stress (for the distance We will keep the x coordinate):
user@machine:~$ cd surfaces/3000
With these commands we have saved the magnitude of wall shear stress to the surface/3000/wallShearStress.xv file. It can now be loaded to xmgrace:
user@machine:~$ xmgrace surface/3000/wallShearStress.xv
This curve should be multipied by 2 (reference velocity and reference density is equal to 1). Within xmgrace go to "Data" > "Transformation" > "Evaluate expression".
This is what you should get:
The FLUENT result for comparison. Values at the developed flow region are similar.
Download the case directory