The Laminar Pipe Flow Example

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) 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
• gnuplot

Meshing

Geometric definition

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.
As the tutorial goes along with the steps of the original Cornell tutorial the necessary modifications of the files coming with the pitzDaily example will be explained. The example can be easily copied with the following terminal command:

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:

• we are not applying any scaling as we define the cordinates in meters
• vertex numbering starts with "0"
• in C++ (yes, this is the language OpenFOAM® was written in) "//" is the comment. Everything after these signs is not executed.
• we would like to have the top of the wedge cylindrical. Therefore the y coordinate of vertex #2-4 is not equal to 0.1. The coordinates belong to the coordinates on a circle with 2.5 deg from the z axis.

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:

user@machine:~$ octave

Calculation of the coordinates:

octave:1> cos(2.5*pi/180)*0.1
ans = 0.099905
octave:2> sin(2.5*pi/180)*0.1
ans = 0.0043619

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:

Please note:

• the bottom face is not a face but an edge, therefore the "corner" vertices are 0 and 1. See the numbering on the OpenFOAM® documentation page sketch here.
• we have 100 elements among the first defined direction (this is the axis direction defined by vortex #0 and #1) and 5 elements among the radial direction (direction of vertex #0 and #2; on the OpenFOAM® homepage wedge explanation sketch vertex #0 and #4 ).
• we do not apply a grading among any of the direction (simpleGrading (1 1 1)). Grading (continously changing element size - expansion) would be necessary towards large changes (with decreasing element sizes) in the calculated fields. This method is commonly used to decrease the discretisation errors during the calculation (e.g. truncation error). In the simpleGrading definition (1 1 1) stands for the ratio of the size of the last and first element in the block direction. Here the cells will have the same size.

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

Download the blockMeshDict file

Launch meshing

As the blockMeshDict is set the meshing can be launched with simply tiping:

user@machine:~$ blockMesh

in the console within the case directory. The ParaView postprocessor (with the OpenFOAM® reader) can be launched with:

user@machine:~$ paraFoam

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.
The field files contain (example: the U field file):

• the dimension of the actual field variable. As dimensions should be matching within the equations do not modify these values.
• the initial condition. The initial condition field can be distributed, not just constant.
• the boundary conditions. The boundary conditions above show the modified state of the U field file. The list of boundary regions and their names should match with the ones specified in the blockMeshDict file. Here their primitive or derived type should be listed with the actual value (List of primitive and derived type boundary conditions with keywords and values they need). According to the Cornell tutorial we specify a uniform inflow with 1 m/s and gradient 0 at the outlet for the velocity. The walls have a no-slip boundary condition with prescribed 0 m/s velocity.
• the pressure boundary condition is gradient 0 at the inlet and 0 at the outlet (the simpleFoam solver uses relative pressures, therefore the pressure level for the outlet and as initial condition can be 0).

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/m³): ν=μ/ρ. 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:

• linear interpolation for calculating the gradients, nu in the laplacian term > gradSchemes, laplacianSchemes
• upwind differenceing for the divergence term > divSchemes
• corrected scheme (explicit non-orthogonal correction) for the surface normal gradients > snGradSchemes, laplacianSchemes

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:

Important options:

• application simpleFoam > put the name of the solver which will used
• startFrom startTime > startTime has to be specified. There are other options (e.g. latestTime to restart an existing calculation from the last available iteration/time step).
• startTime 0 > we start the calculation new. Other (existing) steps could be used from here.
• stopAt endTime > endTime has to be specified. There are other options (e.g. writeNow will stop the solver at the next iteration and write result at that iteration even when originally there would be no result saving at that iteration).
• endTime 3000 > we will calculate 3000 iterations
• deltaT 1 > the same entry is used for steady and time dependent calculations. In transient simulations the timestep should be put here in [s], while for steady the value should be 1.
• writeControl timeStep > result will be written every time steps
• writeInterval 1000 > write result at every 1000 iterations

Calculation

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 &

Some explanation:

• the term "simpleFoam" itself will launch the simpleFoam solver. If the solver is launched simply like this than it will provide its output to the console window. The solver can be stopped with , but we will lose the output however it contains important information.
• the ">" sign in the Linux consoles is a redirection. In this case we redirect the output of the simpleFoam solver to a file called "log". If this file does not exist than it will be created, if it does exist than it will be overwritten (the console output can be appended to an existing file with ">>".).

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:

Monitor residuals

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.

Post Processing

Result files

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:

• rotate - left mouse button.
• zoom - middle mouse button
• pan -right mouse button

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:

Notes:

• output in xmgr format
• cellPoint interpolation to have cell center values only
• an axial and a radial profile is defined with their end points

Launch sampling:

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
user@machine:~$ octave
octave:1> a=load wallShearStress.xyzv;
octave:2> b=ones(rows(a),2);
octave:3> b(:,1)=a(:,1);
octave:4> b(:,2)=sqrt(a(:,4).^2+a(:,5).^2+a(:,6).^2);
octave:5> fid=fopen("wallShearStress.xv","w");
octave:6> for i=1:rows(b)
> fprintf(fid,"%g %g\n",b(i,1),b(i,2));
> end
octave:7> fclose(fid);
octave:8> quit

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".
Select the "Source" and "Destination" graphs the same and type the following into the "Evaluation" field:

y=y*2

This is what you should get:

The FLUENT result for comparison. Values at the developed flow region are similar.

Download the case directory

laminarPipe.tar.gz

Top  Home  Fő oldal  CFD en  CFD hu