Difference between revisions of "Contrib:BondMatt/LaminarPipeFlow"

From CAELinuxWiki
Jump to: navigation, search
 
(45 intermediate revisions by the same user not shown)
Line 1: Line 1:
 
='''Introduction'''=
 
='''Introduction'''=
The problem solved herein is that of a pipe with a uniform velocity field at the inlet and atmospheric pressure at the outlet. The pipe is sufficiently long that the viscous boundary layer is completely formed. A thin slice along the length of the pipe is studied to minimize the time to solve. The inlet velocity is very low (1 m/s) to ensure that flow is laminar. The medium is air and the diameter of the pipe is 0.2m. The Reynolds number, the ratio of inertial and viscous forces, is approximately 100. This is well below 2300, the Reynolds number magnitude commonly accepted as indicating the start of transition from laminar to turbulent flow. Side note: interesting tidbit: laminar flow can occur at extremely high Reynolds numbers (up to 100,000) if the pipe is very smooth and pipe vibrations are controlled. This problem is taken from a FLUENT tutorial from Cornell University, freely available on their website. Final desirables include plots of centerline velocity and wall skin-friction coefficient vs. axial position and the velocity profile at the outlet. Note that all screenshots are displayed at a reduced size and can be clicked on for larger images that may clarify small or cloudy text/fields.
+
The problem solved herein is that of a pipe with a uniform velocity field at the inlet and atmospheric pressure at the outlet. The pipe is sufficiently long that the viscous boundary layer is completely formed. A thin slice along the length of the pipe is studied to minimize the time to solve. The inlet velocity is very low (1 m/s) to ensure that flow is laminar. The medium is air and the diameter of the pipe is 0.2m. The Reynolds number, the ratio of inertial and viscous forces, is approximately 100. This is well below 2300, the Reynolds number magnitude commonly accepted as indicating the start of transition from laminar to turbulent flow. (Interesting tidbit: laminar flow can occur at extremely high Reynolds numbers (up to 100,000) if the pipe is very smooth and pipe vibrations are controlled.) This problem is taken from a FLUENT tutorial from Cornell University, freely available on their website. Final desirables include plots of centerline velocity and wall skin-friction coefficient vs. axial position and the velocity profile at the outlet. Note that all screenshots are displayed at a reduced size and can be clicked on for larger images that may clarify small or cloudy text/fields.
  
I encourage editing and revision of this article. I am not familiar with this software. This contribution is part of my attempts to familiarize myself with Salome and Code Saturne.
+
I encourage editing and revision of this article. I am not familiar with this software. This contribution is part of my attempts to familiarize myself with Salome, OpenFOAM, and Code Saturne.
  
 
='''Software Versions:'''=
 
='''Software Versions:'''=
Line 11: Line 11:
 
Code Saturne V1.4.a (GUI)
 
Code Saturne V1.4.a (GUI)
  
 +
OpenFOAM 2.0
  
 
='''Fluent Tutorial'''=
 
='''Fluent Tutorial'''=
Line 17: Line 18:
  
 
='''Geometry / Mesh - Salome'''=
 
='''Geometry / Mesh - Salome'''=
Creation of solid box: [[Media:lam-pipe-1a.ogg]]
 
  
The geometry is a box. The Cornell tutorial uses a 2-D box but Code Saturne can only work with 3-D geometry. The box represents a cross section of a pipe along its length cut in half to take advantage of the flow symmetry around the axis of the pipe. I created a box using the 'create a box' icon which brings up the same dialogue box if one goes to 'New Entity' and then 'Primitives' and then 'Box' in the task bar.
+
The geometry is a box. The Cornell tutorial uses a 2-D rectangle but Code Saturne can only work with 3-D geometry. The rectangle represents a cross section of a pipe along its length cut in half to take advantage of the flow symmetry around the axis of the pipe. I created a box using the 'create a box' icon which brings up the same dialogue box if one goes to 'New Entity' and then 'Primitives' and then 'Box' in the task bar.
  
  
Line 137: Line 137:
  
  
The next folder is 'Calculation control' and this is where I begin to become uncertain about what each value is. The first file is 'Steady management'. I leave the relaxation coefficient as is but change the 'Iterations number (restart included)' to 100. In the fluent tutorial the results converged to the desired residual after 46 iterations. With Fluent one specifies the maximum number of iterations without convergence to the desired residual. It looks like Saturne just runs a given number of iterations. It would be nice to see plots of residuals like you can get with Fluent. For a complex problem I can only guess that one would run the problem with an increasing number of iterations and do a quick analysis of the results each time to see if they are reasonable and improving. Perhaps there is a way of plotting residuals in the post-processor?
+
The next folder is 'Calculation control' and this is where I begin to become uncertain about what each value is. The first file is 'Steady management'. I leave the relaxation coefficient as is but change the 'Iterations number (restart included)' to 100. In the fluent tutorial the results converged to the desired residual after 46 iterations. With Fluent one specifies the maximum number of iterations without convergence to the desired residual. It looks like Saturne just runs a given number of iterations. It would be nice to see plots of residuals like you can get with Fluent. For a complex problem I can only guess that one would run the problem with an increasing number of iterations and do a quick analysis of the results each time to see if they are reasonable and improving. Perhaps there is a way of plotting residuals in the post-processor? The developers over at the EDF Code Saturne forums gave me some advice on the topic of residuals but I have not made any headway on that topic yet.
  
  
Line 149: Line 149:
  
  
The next folder is 'Numerical parameters' and the first file is 'Equation parameters'. Under the 'scheme' tab I changed the Scheme for each velocity component to 'Upwind'. This was done for several reasons. I am assuming that this refers to the differencing scheme used to evaluate the rate of change of each velocity component. In the Fluent tutorial a second-order upwind scheme is used. The Saturne scheme may not be second order but it is at least upwind. This is important because upwind schemes take into account the direction of the flow. The solution can become unstable if a scheme that does not consider direction of flow is used. Central differencing, in its simplest form, is more accurate than upwind differencing (in its simplest form). For the 'Equation parameters' file no changes were made under the 'solver' tab. In the 'global parameters' file the 'Pseudo-coupled velocity-pressure solver' box was selected. I did not want to check this box and have been attempting to change other settings (including the differencing scheme for the velocity components). This box was not checked in the flow over a step tutorial. Any thoughts?
+
The next folder is 'Numerical parameters' and the first file is 'Equation parameters'. Under the 'scheme' tab I changed the Scheme for each velocity component to 'Upwind'. This was done for several reasons. I am assuming that this refers to the differencing scheme used to evaluate the rate of change of each velocity component. In the Fluent tutorial a second-order upwind scheme is used. I think most journals require 2nd order differencing schemes. The Saturne scheme may not be second order but it is at least upwind. This is important because upwind schemes take into account the direction of the flow. The solution can become unstable if a scheme that does not consider direction of flow is used. Central differencing, in its simplest form, is more accurate than upwind differencing (in its simplest form). There is a whole lot of theory here that I do not claim to understand very well. From what I can recall off the top of my head, upwind differencing is best for certain discretizations. Therefore, without any other advice I would upwind difference as much as possible and use second order if accuracy was desired. For the 'Equation parameters' file no changes were made under the 'solver' tab. In the 'global parameters' file the 'Pseudo-coupled velocity-pressure solver' box was selected. I did not want to check this box and have been attempting to change other settings (including the differencing scheme for the velocity components). This box was not checked in the flow over a step tutorial.
  
  
Line 168: Line 168:
  
 
Save your case file. The last step (for the solver) is to click on the 'Code_Saturne batch running' icon. This simulation should run fairly quickly. I use the system monitor to keep an eye on progress but this one runs almost too fast.
 
Save your case file. The last step (for the solver) is to click on the 'Code_Saturne batch running' icon. This simulation should run fairly quickly. I use the system monitor to keep an eye on progress but this one runs almost too fast.
 
  
 
='''Post-Processing / Paraview'''=
 
='''Post-Processing / Paraview'''=
Line 189: Line 188:
  
  
The result is quite different from the Fluent tutorial results. Theoretically, the v component of the centerline velocity in the fully developed portion of the flow should be 2 m/s or twice that of the v component of the uniform inlet velocity field. Where did I go wrong? This is not a simulation of a cross section of a pipe. It is a simulation of a cross section of an infinitely wide rectangular conduit or flow between two infinitely wide 8m long plates with a distance of 0.2m between them. I was not clever enough to figure this out on my own, the developers of Code Saturne over at the EDF forum enlightened me. Kudos to them.
+
The result is quite different from the Fluent tutorial results. Theoretically, the v component of the centerline velocity in the fully developed portion of the flow should be 2 m/s or twice that of the v component of the uniform inlet velocity field. Where did I go wrong? This is not a simulation of a cross section of a pipe. It is a simulation of a cross section of an infinitely wide rectangular conduit or flow between two infinitely wide 8m long plates with a distance of 0.2m between them. I was not clever enough to figure this out on my own, the developers of Code Saturne over at the EDF forum enlightened me. Kudos to them. I believe OpenFOAM has the ability to take a 2D mesh and create a 3D flow domain with a variety of options. I have yet to run a simulation but it looks like there is an option to create a 5 degree segment of a cylinder.
  
  
Line 260: Line 259:
  
 
[[Image:centerline-velocity.jpg|800x500px]]
 
[[Image:centerline-velocity.jpg|800x500px]]
 +
 +
 +
='''Geometry/Mesh, boundary conditions, initial conditions - OpenFOAM'''=
 +
For some reason I thought OpenFOAM had 2D solvers akin to Fluent. I was not quite correct. It does have a solver for axisymmetric problems using a small wedge.
 +
 +
However, since I am just learning to use the software I start with a basic simulation: a rectangular prism.
 +
 +
I started with the cavity tutorial covered in the OpenFOAM manual. The most basic continuation that seemed reasonable was to make a rectangular prism with dimensions of 0.1 meters by 0.1 meters by 8 meters similar to that of the pipe in the previous Code Saturn models.
 +
 +
So the blockMeshDict looks like:
 +
 +
<code>
 +
vertices
 +
(
 +
    (0 0 0)
 +
    (0.1 0 0)
 +
    (0.1 0.1 0)
 +
    (0 0.1 0)
 +
    (0 0 8)
 +
    (0.1 0 8)
 +
    (0.1 0.1 8)
 +
    (0 0.1 8)
 +
);
 +
 +
blocks
 +
(
 +
    hex (0 1 2 3 4 5 6 7) (10 10 800) simpleGrading (1 1 1)
 +
);
 +
 +
edges
 +
(
 +
);
 +
 +
boundary
 +
(
 +
    fixedWalls
 +
    {
 +
        type wall;
 +
        faces
 +
        (
 +
            (3 7 6 2)
 +
        );
 +
    }
 +
    frontAndBack
 +
    {
 +
        type empty;
 +
        faces
 +
        (
 +
            (1 5 4 0)
 +
            (0 4 7 3)
 +
            (2 6 5 1)
 +
        );
 +
    }
 +
    inlet
 +
    {
 +
        type patch;
 +
        faces
 +
        (
 +
            (1 2 3 0)
 +
        );
 +
    }
 +
    outlet
 +
    {
 +
        type patch;
 +
        faces
 +
        (
 +
            (5 6 7 4)
 +
        );
 +
      }
 +
);
 +
);</code>
 +
 +
A logical next step would be initial conditions. In directory '0' the initial conditions for the pressure and velocity are set in two text files named 'p' and 'U'.
 +
 +
The text file 'p' might contain:
 +
 +
<code>
 +
boundaryField
 +
{
 +
    fixedWalls
 +
    {
 +
        type            zeroGradient;
 +
    }
 +
 +
    frontAndBack
 +
    {
 +
        type            empty;
 +
    }
 +
   
 +
    outlet
 +
    {
 +
        type            fixedValue;
 +
        value          uniform 0;
 +
    }
 +
 +
    inlet
 +
    {
 +
        type            zeroGradient;
 +
    }
 +
}
 +
</code>
 +
 +
Since there is no flow through or tangential to walls it makes sense that there should be no pressure gradient at any point on a wall.
 +
The 'frontAndBack' conditions are for planes of symmetry. There is a plane of symmetry condition in OpenFOAM but it is not used here. Consistent with other tutorials these surfaces were set to type 'empty' when creating the mesh. The outlet is a pressure outlet with a gauge pressure of 0 (atmospheric). I struggle to understand the zeroGradient condition at the inlet since in one direction there is clearly a pressure gradient driving the flow. In the other two principal directions the gradient certainly should be zero. Perhaps the velocity imposed in the text file named 'U' will overwrite the null pressure gradient?
 +
 +
And the main body of text file 'U':
 +
 +
<code>
 +
boundaryField
 +
{
 +
    inlet
 +
    {
 +
        type            fixedValue;
 +
        value          uniform (0 0 1);
 +
    }
 +
 +
    fixedWalls
 +
    {
 +
        type            fixedValue;
 +
        value          uniform (0 0 0);
 +
    }
 +
 +
    frontAndBack
 +
    {
 +
        type            empty;
 +
    }
 +
 +
    outlet
 +
    {
 +
        type            zeroGradient;
 +
    }
 +
}
 +
</code>
 +
 +
Nothing too shocking here.
 +
 +
='''Initial results - OpenFOAM'''=
 +
 +
Similar to the first simulations run with Code Saturn the axial velocity profile is shown below for a simulation with a termination time of 2 seconds and a time step of 0.005 seconds. The time step was calculated considering the Courant condition where the time for the fluid to flow completely through any direction in the smallest cell should be greater than the timestep. All cells are perfectly cubic with edge lengths of 0.1 meters divided by 10 or 0.01 meters. The velocity at the inlet is 1 m/s but once the velocity profile develops the maximum velocity will be approximately 1.5 m/s. The timestep could probably be set slightly larger at 0.006 or possibly even 0.0065 seconds but I rounded to the nearest 0.005 seconds.
 +
 +
[[File:OpenFOAM_results1.png]]
 +
 +
 +
snappyHexMesh was then used to obtain a mesh of a quarter cylinder. This required creating geometry, in ASCII STL file format, of the region of the blockmesh to be removed to obtain a quarter cylinder:
 +
 +
 +
[[File:Geometry_Boolean_quarter_cylinder.png]]
 +
 +
 +
Based upon previous results the domain was reduced to a pipe with a length of 4 meters. This is sufficient to full develop the velocity profile and reduces computational expense. snappyHexMesh is no speed demon when one does not know what they are doing.
 +
 +
 +
The snappyHexMesh dictionary was taken from some other tutorial. There are a few decent references out there including:
 +
http://www.hydroniumion.de/general/snappyhexmesh-tutorial/
 +
 +
 +
My initial snappyHexMesh dictionary and the STL file:
 +
[[File:snappyHexMeshDict.txt]]
 +
[[File:GeoSTL_half.stl.txt]]
 +
 +
 +
For both files the .txt extension must be removed.
 +
 +
 +
With the previous snappyHexMesh dictionary the resulting mesh looks like:
 +
 +
[[File:SnappeHexMeshNumeroUno.png]]
 +
 +
 +
The previous snappyHexMeshdict dictionary does not terminate normally with OpenFOAM v2.11. The following dictionary was revised to include missing keywords. The filename must be altered to remove the '2.txt'
 +
 +
 +
[[File:snappyHexMeshDict2.txt]]
 +
 +
 +
The following figure shows the mesh in relation to the STL geometry to which the mesh was snapped. To view this particular mesh the time was advanced in ParaView using the 'Next frame' button.
 +
 +
[[File:snappyHexMeshSnapped.png]]
 +
 +
 +
snappyHexMesh apparently creates a few different meshes (2 or 3?), sort of intermediate steps. It creates each of these meshes at different times. Upon running snappyHexMesh one should notice the creation of a few folders with names that will be associated with the time step specified in the control dictionary. To run a simulation I have been trying to change the start time to the time associated with the final mesh generated by snappyHexMesh. However, this leads to a few problems. The final mesh generated by snappyHexMesh has no boundary conditions associated with it. I would have hoped snappyHexMesh would automatically create new boundary conditions based upon the block mesh boundary conditions I specified. This is not a big problem, I just copy the files 'p' and 'U' from folder '0'. At this point I found that the wall boundary conditions in the block mesh (which I had named 'fixedWalls', consistent with the cavity tutorial) no longer have associated nodes or cell faces in the mesh. snappyHexMesh removed elements in this region and created a new group of cell faces and/or nodes named 'Geo_solid'.
 +
 +
My solution was to change the cell faces/nodes group 'fixedWalls' in my blockmesh to 'Geo_solid' and, correspondingly, change my boundary conditions to be consistent.
 +
 +
To run icoFoam I change the start time to the time at which snappyHexMesh generated the final mesh.
 +
 +
However, this does not work particularly well. When I run icoFoam the Courant number diverges very quickly:
 +
 +
 +
[[File:CourantNumberDivergence.png]]
 +
 +
 +
Running checkMesh there seems to be a problem. Note that checkMesh is indicating a problem for an intermediate mesh. It is not shown but there is also a problem with the final mesh.
 +
 +
 +
[[File:CheckMesh.png]]
 +
 +
 +
A quick search online seemed to indicate that this issue with checking the mesh is associated with the use of the 'empty' boundary condition type. A simple solution was to change the type from 'empty' to 'symmetryPlane' in the blockmesh dictionary and the initial/boundary conditions for pressure and velocity.
 +
 +
 +
='''Coarser mesh - OpenFOAM'''=
 +
 +
 +
The previous mesh has very small cells that may require a very small time step. Additionally, a logical next step was to fiddle with the snappyHexMesh dictionary and see what happens.
 +
 +
Decreasing the maxLocalCells and maxGlobalCells parameters each by an order of magnitude yielded the following mesh:
 +
 +
 +
[[File:PoorQualityCoarseMesh.png]]
 +
 +
 +
icoFoam will produce the following reasonable results with a time step of 0.0015 seconds:
 +
 +
 +
[[File:ResultsCoarseMesh.png]]
 +
 +
 +
I wanted to try to obtain a nicer looking mesh and found a very good resource:
 +
https://sites.google.com/site/snappywiki/snappyhexmesh/cylinder-case
 +
 +
 +
Changing level (1 2) to (0 0) under refinementSurfaces yielded the following mesh.
 +
This keyword, refinementSurfaces, controls cell splitting by proximity to a specified surface.
 +
The integers in brackets indicate minimum and maximum refinement (see section 5.4.3 of the OpenFOAM users manual).
 +
 +
 +
[[File:Coarse_hex_mesh.png]]
 +
 +
='''Skin friction coefficient - OpenFOAM'''=
 +
 +
 +
Another result in the Fluent tutorial, in addition to the axial velocity profile, is the skin friction coefficient:
 +
 +
[[File:skin_friction_formula.png]]
 +
 +
 +
OpenFOAM has a tool akin to foamCalc for calculating shear stress: wallShearStress
 +
In this case I used wallShearStress -latestTime to calculate shear stress at the termination time when a steady state condition had been reached.
 +
However, in order to even use this tool I had to modify the model. I believe there may be more than one way to approach this but this was my solution.
 +
 +
In the folder named 'constant' I opened the transportProperties dictionary and replaced what was there with:
 +
 +
transportModel  Newtonian;
 +
 +
rho            rho [ 1 -3 0 0 0 0 0 ] 1;
 +
 +
nu              nu [ 0 2 -1 0 0 0 0 ] 2e-03;
 +
 +
 +
These properties of the fluid were made consistent with the Fluent tutorial. Two other corrresponding additions to the model were the RASProperties and turbulenceProperties dictionaries also in the folder named constant. Both files may not be necessary. The RASProperties dictionary contains the following lines (see table 3.9 of the OpenFOAM users manual).
 +
 +
RASModel        laminar;
 +
 +
turbulence      off;
 +
 +
printCoeffs    off;
 +
 +
 +
The turbulenceProperties dictionary contains:
 +
 +
simulationType laminar;
 +
 +
After re-running the model with the new dictionaries 'wallShearStress -latestTime' can be run in the OpenFOAM terminal. Opening the results in Paraview (command: 'paraFoam'), jumping to the last time step, and plotting the wall shear stress over a line through (0.0707 0.0707 0.1) and (0.0707 0.0707 3.9) yields the following plot:
 +
 +
[[File:wallShearStressPlot.png]]
 +
 +
The results are a little different than those obtained with Fluent here: https://confluence.cornell.edu/download/attachments/118475269/SkinFrictionPlot_Full.png
 +
 +
However, the magnitude of the shear stress is of course half of the skin friction coefficient due to the one-half in the denominator and the value of unity for the density and velocity. If the points for the line plot are near the planes of symmetry the wall shear stress is under-predicted. The shear stress should just be the viscosity multiplied by the derivative of velocity with respect to distance from the wall assuming perfectly axial flow. Looking at velocity contour plots (magnitude and components) the only significant component is axial and the velocity is relatively constant at a specific axial position and radius (i.e. no radial flow). Perhaps the shear discrepancy with choice of line plot points is the result of how the numerical derivative is performed and the location? With no data outside of the mesh numerical derivatives at boundaries might use null values when unknown or switch to a different finite difference scheme.

Latest revision as of 01:23, 9 February 2013

Introduction

The problem solved herein is that of a pipe with a uniform velocity field at the inlet and atmospheric pressure at the outlet. The pipe is sufficiently long that the viscous boundary layer is completely formed. A thin slice along the length of the pipe is studied to minimize the time to solve. The inlet velocity is very low (1 m/s) to ensure that flow is laminar. The medium is air and the diameter of the pipe is 0.2m. The Reynolds number, the ratio of inertial and viscous forces, is approximately 100. This is well below 2300, the Reynolds number magnitude commonly accepted as indicating the start of transition from laminar to turbulent flow. (Interesting tidbit: laminar flow can occur at extremely high Reynolds numbers (up to 100,000) if the pipe is very smooth and pipe vibrations are controlled.) This problem is taken from a FLUENT tutorial from Cornell University, freely available on their website. Final desirables include plots of centerline velocity and wall skin-friction coefficient vs. axial position and the velocity profile at the outlet. Note that all screenshots are displayed at a reduced size and can be clicked on for larger images that may clarify small or cloudy text/fields.

I encourage editing and revision of this article. I am not familiar with this software. This contribution is part of my attempts to familiarize myself with Salome, OpenFOAM, and Code Saturne.

Software Versions:

Ubuntu 8.04 LTS

Salome V4.1.4

Code Saturne V1.4.a (GUI)

OpenFOAM 2.0

Fluent Tutorial

The Cornell University Fluent tutorial can be found here Cornell Tutorial


Geometry / Mesh - Salome

The geometry is a box. The Cornell tutorial uses a 2-D rectangle but Code Saturne can only work with 3-D geometry. The rectangle represents a cross section of a pipe along its length cut in half to take advantage of the flow symmetry around the axis of the pipe. I created a box using the 'create a box' icon which brings up the same dialogue box if one goes to 'New Entity' and then 'Primitives' and then 'Box' in the task bar.


Box.jpg


Dx is the width of the box in the x-direction, starting at the origin. My dimensions: Dx=0.001, Dy=8, and Dz=0.1 Units are meters. The next step is to create the groups of faces for boundary conditions. As in the Cornell tutorial, half the pipe is studied to take advantage of symmetry. There are three symmetry B.C.s, one inlet B.C., one outlet B.C., and a wall B.C.


Box2.jpg


Groups can be created by going to 'New Entity', 'Group', and then 'Create'. When creating a group the elements are selected by type (points, lines, planes, solids) and added with the 'add' button. Descriptive group names should be used to make applying B.C.s easier. Once all elements for a group are selected and the name is specified the 'Apply' button is clicked.


Groups.jpg


Groups of edges were also created by direction so that the same mesh that is used in the Cornell tutorials could be created. In other words lines running in the x-direction were grouped into a group named 'xlines' and so on. This allows one dimensional submeshes to be created with the number of intervals on the edges specified. These submeshes are 1-dimensional and use the 'wire discretization' algorithm with the 'Nb. Segments' (number of segments) hypothesis. The main mesh used the 'Hexahedron' algorithm (3D) and 'Quadrangle (Mapping)' algorithm (2D). For the main mesh no 1D algorithm is specified since all edges are discretized by the submeshes.


Mesh1-3D.jpg


Mesh1-2D.jpg


Mesh1-submesh1.jpg


Mesh1-submesh1-hypothesis.jpg


Mesh1-submesh2-hypothesis.jpg


Mesh1-submesh3-hypothesis.jpg


The geometry groups are used to create groups of nodes. One method is to go to 'Mesh' and then 'Create Groups from Geometry' in the taskbar. The mesh and the geometry group are selected and the 'Apply' button selected.


Mesh-groups.jpg


Solver - Code Saturne

The CFD study was started by launching the 'Code-Saturne-Wizard' and creating a case. This dialogue box can be used to remember/select the right type of mesh to export (MED). To export the mesh it is just a matter of going back to Salome and right clicking on the mesh in the tree and selecting 'Export to MED file'.


ExportMEDmesh.jpg


Launch-cfd-wizard.jpg


The GUI will launch and directories will be created. I specified that the directories should be created on the desktop but they can be anywhere. The next step is to create a new case. There is an icon provided in the navigation toolbar for this in the form of a piece of paper with a pencil. A new case could also be created by going to 'File' and then 'New file'. On the left hand side of the GUI window is a series of folders and files. Each file corresponds to an aspect of the solver including boundary conditions, output file type, fluid properties, and number of iterations, just to name a few. If you should decide to close the Saturne GUI it can be started by opening the directory that was created with the CAELinux wizard and opening the DATA folder. In this folder is a file entitled 'SatuneGUI'. If you double click it and select the 'run in terminal' option from the dialogue box that comes up, the Saturne GUI will appear. You will need to open the .xml file that was saved the last time you were working on that case. As with most computer based work saving often is always a good idea.


New-case.jpg


Restarting-GUI.jpg


I follow the folders from top to bottom going through each file. The first folder is 'Calculation environment' and the first file is 'Identity and paths'. I do not edit anything here. The next file is 'Solution domain'. There are several tabs for this file but the first and last are the only significant ones for this problem. In the 'meshes' tab you should confirm that your mesh is listed. Under the 'stand-alone running' tab the output file type is selected. I used the Ensight Gold format which is compatible with Paraview. I have been unable to get results with the MED format which is compatible with the post processor in Salome. The simulation appears to run just fine but there are no result files created. The next file is 'Definition of volume zones'. I did not edit anything in this file.


The next folder is 'Thermophysical models'. In the first file ('Calculation features') the flow type must be changed to 'steady flow'. This indicates that the flow does not change at each point with respect to time, ie. the velocity is constant at each point with respect to time. I did not edit the next file ('Mobile mesh'). The next file is 'Turbulence models'. 'No model (ie. laminar flow)' is selected for this problem.


Calc-features.jpg


Turb-model.jpg


The next (and last) file in the 'Thermophysical models' folder is 'Initialization'. Saturne uses the finite volume CFD solver technique which essentially just conserves mass, energy, and momentum over each element. Partial differential equations for each of these conservation laws (mass: continuity, momentum: navier stokes) are integrated over the faces of the element and the integrals are evaluated numerically. The result is a set of linear algebraic equations. These equations are solved using an indirect, iterative technique. In other words a solution is assumed and then all but one of the assumed or previously calculated values are used to re-evaluate each unknown. This is done until a convergent solution is reached. In this 'Initialization' file the initial guess for all velocity magnitudes is specified. I assume that the u and w components are 0 m/s (these are known) and the v component is 1 m/s. The actual v-component will vary from 0 - 2 m/s but 1 m/s is the average.


Initialization.jpg


The next folder is 'Physical properties'. The first file is 'Reference values' which is where atmospheric pressure is set. The next file is 'Fluid properties'. As per the Cornell tutorial (Cornell Fluent tutorial) the density is 1 kg per cubic meter and the dynamic viscosity is 2e-3 kg/m/s. We are not concerned with the temperature field so the specific heat magnitude is fine as is. I may have previously stated that the fluid is air but the viscosity is much too high (comparable to very cold water). The last file in this folder is 'Gravity, hydrostatic pressure'. I did not modify anything in this file.


Fluid-properties.jpg


The next folder is 'Additional scalars'. No files in this folder were modified.


The next folder is 'Boundary conditions'. The first file is 'Definition of boundary regions'. There are several steps that need to be taken with this first file. The first step is to click on the folder/magnifying glass icon annotated 'Import groups and references from Preprocessor listing'. A window will open and the file 'listenv.pre' should be selected. Five regions should then appear under 'Definition of boundary regions'. There are four columns: Label, Zone, Nature, and Localisation. The 'Localisation' column will contain the group names created in Salome and can be used to set the boundary condition types.


B-c-def1.jpg


B-c-def2.jpg


B-c-def3.jpg


B-c-def4.jpg


B-c-def5.jpg


B-c-def6.jpg


The next file in the 'Boundary conditions' folder is 'Dynamic variables bound. cond.' Two boundary conditions should appear: the inlet and the wall. The default settings for the wall are left as is. The velocity components for the inlet must be set: u = 0 m/s, v = 1 m/s, and w = 0 m/s. The next file in this folder has no settings that can modified for this problem.


Inlet-b-c.jpg


The next folder is 'Calculation control' and this is where I begin to become uncertain about what each value is. The first file is 'Steady management'. I leave the relaxation coefficient as is but change the 'Iterations number (restart included)' to 100. In the fluent tutorial the results converged to the desired residual after 46 iterations. With Fluent one specifies the maximum number of iterations without convergence to the desired residual. It looks like Saturne just runs a given number of iterations. It would be nice to see plots of residuals like you can get with Fluent. For a complex problem I can only guess that one would run the problem with an increasing number of iterations and do a quick analysis of the results each time to see if they are reasonable and improving. Perhaps there is a way of plotting residuals in the post-processor? The developers over at the EDF Code Saturne forums gave me some advice on the topic of residuals but I have not made any headway on that topic yet.


Steady-management.jpg


The second file in the 'Calculation control' folder is 'Output control'. I did not modify anything in this file. The last file in this folder is 'Solution control'. I unchecked the u and w velocity component boxes under 'Print in listing' and 'Post-processing'. I found that if either of the u-velocity component checkboxes were selected, the u velocity component would be the only one in the output. The u and w velocity components should be 0 m/s at all points anyways.


Solution-control.jpg


The next folder is 'Numerical parameters' and the first file is 'Equation parameters'. Under the 'scheme' tab I changed the Scheme for each velocity component to 'Upwind'. This was done for several reasons. I am assuming that this refers to the differencing scheme used to evaluate the rate of change of each velocity component. In the Fluent tutorial a second-order upwind scheme is used. I think most journals require 2nd order differencing schemes. The Saturne scheme may not be second order but it is at least upwind. This is important because upwind schemes take into account the direction of the flow. The solution can become unstable if a scheme that does not consider direction of flow is used. Central differencing, in its simplest form, is more accurate than upwind differencing (in its simplest form). There is a whole lot of theory here that I do not claim to understand very well. From what I can recall off the top of my head, upwind differencing is best for certain discretizations. Therefore, without any other advice I would upwind difference as much as possible and use second order if accuracy was desired. For the 'Equation parameters' file no changes were made under the 'solver' tab. In the 'global parameters' file the 'Pseudo-coupled velocity-pressure solver' box was selected. I did not want to check this box and have been attempting to change other settings (including the differencing scheme for the velocity components). This box was not checked in the flow over a step tutorial.


Equation-parameters.jpg


Global-parameters.jpg


The next folder is 'Calculation management'. The first file is 'Memory Management'. No changes were made to this file. The next file (which again was not changed) is 'Start/Restart'. The last file is 'Prepare batch calculation'. The first step with this file was to select the batch script file. This is done by clicking on the folder/magnifying glass icon annotated 'Select the batch script file'. The file is named 'runcase'. More icons will appear once this file is selected. The next step is to select the 'Advanced options' icon. This is not necessary but only suggested. The first box in the window that comes up allows the location of the temporary files directory to be selected. It can be useful to know where this is because there are a bunch of backed up files in there.


Select-batch-file.jpg


Temp-directory.jpg


Save your case file. The last step (for the solver) is to click on the 'Code_Saturne batch running' icon. This simulation should run fairly quickly. I use the system monitor to keep an eye on progress but this one runs almost too fast.

Post-Processing / Paraview

Paraview can be launched from the Applications/CAELinux menu. The Ensight Gold results file is loaded by going to File and then Open in the menu bar. The file is located in the directory created for the study in the RESU folder. For Ensight Gold result files there will be another directory along the lines of CHR.ENSIGHT.#.


OpenResults.jpg


I plotted the centerline velocity first. Since I was not interested in pressure I unchecked the 'Pressure' and 'total pressure' checkboxes under 'Properties - Cell/Point Array Status'. It is also necessary to change the 'Byte Order' from 'BigEndian' to 'LittleEndian'.


CHRcase.jpg


To plot the centerline velocity there is a Filter called 'Plot Over Line' in the 'Alphabetical' portion of the menu. To plot the centerline velocity the points for the line could be (0,0,0) and (0,8,0).


Centerline.jpg


The result is quite different from the Fluent tutorial results. Theoretically, the v component of the centerline velocity in the fully developed portion of the flow should be 2 m/s or twice that of the v component of the uniform inlet velocity field. Where did I go wrong? This is not a simulation of a cross section of a pipe. It is a simulation of a cross section of an infinitely wide rectangular conduit or flow between two infinitely wide 8m long plates with a distance of 0.2m between them. I was not clever enough to figure this out on my own, the developers of Code Saturne over at the EDF forum enlightened me. Kudos to them. I believe OpenFOAM has the ability to take a 2D mesh and create a 3D flow domain with a variety of options. I have yet to run a simulation but it looks like there is an option to create a 5 degree segment of a cylinder.


Geometry - Salome - Revised Geometry

When I re-did this simulation I went with a 45 degree segment of the pipe since this is very simple to construct in Salome. I approached creation of this geometry by creating a cylinder with a radius of 0.1m and a length of 8m and a box with dx=dy=0.1m and dz=8m. Then a boolean operation can be used to extract the intersecting volume. Unfortunately, this geometry is very difficult to mesh such that it represents the mesh in the Cornell Fluent tutorial. A coarse mesh that roughly corresponds can be created but the results are extremely erroneous (the boundaries become extremely distorted). A small slice might allow a more similar mesh and I may experiment with this at a later time. I have recently found a time saving alternative to creating the Boundary Condition groups in Salome by using the explode function to make 'sub-geometries' for each entity for all faces and edges. These 'sub-geometries' can then be used in creation of the mesh to control the local sizes. Once I have exploded the entity to get all of the edges and faces I systematically rename each one based upon the boundary conditions, geometry, or other relevance. In this case all faces were named after the boundary conditions that will be applied (symmetry, inlet, outlet, and wall). The edges were named after the two faces on either side.


Create-box.jpg


Create-cylinder.jpg


Boolean-common.jpg


Extract-faces.jpg


Extract-edges.jpg


Geometry-object-tree-before-renaming.jpg


Object-tree-after-renaming.jpg


Meshing - Salome - Revised Geometry

I struggled attempting to use hexahedron (box-like) elements which are more accurate and in the end resorted to using tetrahedral (pyramid shaped). I created the main mesh using boolean extracted 'Common' geometry, the 'Tetrahedron (Netgen)' 3D Algorithm, and the 'Triangle (Mefisto)' 2D Algorithm with the 'Length from edges' hypothesis. I may have approached the creation of submeshes somewhat differently for this revised mesh as I selected multiple edge groups if the desired number of intervals was the same. In retrospect I should have created groups of edges based upon the number of intervals but the difference in effort is pretty well negligible.


3DmeshAlgorithm.jpg


2DmeshAlgorithm.jpg


AxisSubmesh.jpg


Submesh2.jpg


I ended up using 300 intervals in the axial direction, 15 intervals for the straight edges in the x and y directions at the inlet and outlet, and 15 intervals along the curved edges at the inlet and outlet. I experimented with the mesh extensively to find the optimal interval numbers. The critical edges were those in the axial direction as a small number would result in very coarse meshes that did not approximate the curvature at the wall very well. It also often came down to a trade-off between oscillations in the fully developed centerline velocity and time to run simulation. Even using both cores of my CPU I found this simulation would take approximately 5 minutes to run. A smaller sector would reduce the number of elements for a given average element size which would reduce the time to run the simulation but I am unsure about the boundary conditions. I am not familiar with the periodicity boundary conditions and will need to spend some time familiarizing myself. One issue with this mesh is that it does not compare with the mesh used in the Fluent tutorials. The resulting mesh is shown below:


Mesh&tree.jpg


Solver - Code Saturne - Revised Geometry

The following settings were used: steady flow, no turbulence model, velocity initialization of (0,0,1) m/s, density of 1.0 kg/m^3, viscosity of 0.002 kg/m/s, inlet velocity of (0,0,1) m/s, 50 iterations, and velocity components u and v removed from results.


Calculation-features&turbModels.jpg


Initialization&fluidproperties.jpg


Boundary-conditions&inletBC.jpg


Steady-management&solution-control.jpg


Post Processing - Paraview - Revised Geometry

Using the same filter (plot over line) the centerline velocity was plotted. An interesting tidbit compliments of Code Saturne support at EDF: velocities come from cell centers so the velocity is not exactly from the centerline. This makes sense since the minimum velocities are greater than zero (should be zero at the wall). Generally the pressures and velocities are computed on separate grids in CFD simulations (incompressible flows?) since pressure drives the flow.


Centerline-velocity.jpg


Geometry/Mesh, boundary conditions, initial conditions - OpenFOAM

For some reason I thought OpenFOAM had 2D solvers akin to Fluent. I was not quite correct. It does have a solver for axisymmetric problems using a small wedge.

However, since I am just learning to use the software I start with a basic simulation: a rectangular prism.

I started with the cavity tutorial covered in the OpenFOAM manual. The most basic continuation that seemed reasonable was to make a rectangular prism with dimensions of 0.1 meters by 0.1 meters by 8 meters similar to that of the pipe in the previous Code Saturn models.

So the blockMeshDict looks like:

vertices (

   (0 0 0)
   (0.1 0 0)
   (0.1 0.1 0)
   (0 0.1 0)
   (0 0 8)
   (0.1 0 8)
   (0.1 0.1 8)
   (0 0.1 8)

);

blocks (

   hex (0 1 2 3 4 5 6 7) (10 10 800) simpleGrading (1 1 1)

);

edges ( );

boundary (

   fixedWalls
   {
       type wall;
       faces
       (
           (3 7 6 2)
       );
   }
   frontAndBack
   {
       type empty;
       faces
       (
           (1 5 4 0)
           (0 4 7 3)
           (2 6 5 1)
       );
   }
   inlet
   {
       type patch;
       faces
       (
           (1 2 3 0)
       );
    }
    outlet
    {
       type patch;
       faces
       (
           (5 6 7 4)
       );
     }

); );

A logical next step would be initial conditions. In directory '0' the initial conditions for the pressure and velocity are set in two text files named 'p' and 'U'.

The text file 'p' might contain:

boundaryField

{
   fixedWalls
   {
       type            zeroGradient;
   }
   frontAndBack
   {
       type            empty;
   }
   
   outlet
   {
       type            fixedValue;
       value           uniform 0;
   }
   inlet
   {
       type            zeroGradient;
   }

}

Since there is no flow through or tangential to walls it makes sense that there should be no pressure gradient at any point on a wall. The 'frontAndBack' conditions are for planes of symmetry. There is a plane of symmetry condition in OpenFOAM but it is not used here. Consistent with other tutorials these surfaces were set to type 'empty' when creating the mesh. The outlet is a pressure outlet with a gauge pressure of 0 (atmospheric). I struggle to understand the zeroGradient condition at the inlet since in one direction there is clearly a pressure gradient driving the flow. In the other two principal directions the gradient certainly should be zero. Perhaps the velocity imposed in the text file named 'U' will overwrite the null pressure gradient?

And the main body of text file 'U':

boundaryField {

   inlet
   {
       type            fixedValue;
       value           uniform (0 0 1);
   }
   fixedWalls
   {
       type            fixedValue;
       value           uniform (0 0 0);
   }
   frontAndBack
   {
       type            empty;
   }
   outlet
   {
       type            zeroGradient;
   }

}

Nothing too shocking here.

Initial results - OpenFOAM

Similar to the first simulations run with Code Saturn the axial velocity profile is shown below for a simulation with a termination time of 2 seconds and a time step of 0.005 seconds. The time step was calculated considering the Courant condition where the time for the fluid to flow completely through any direction in the smallest cell should be greater than the timestep. All cells are perfectly cubic with edge lengths of 0.1 meters divided by 10 or 0.01 meters. The velocity at the inlet is 1 m/s but once the velocity profile develops the maximum velocity will be approximately 1.5 m/s. The timestep could probably be set slightly larger at 0.006 or possibly even 0.0065 seconds but I rounded to the nearest 0.005 seconds.

OpenFOAM results1.png


snappyHexMesh was then used to obtain a mesh of a quarter cylinder. This required creating geometry, in ASCII STL file format, of the region of the blockmesh to be removed to obtain a quarter cylinder:


Geometry Boolean quarter cylinder.png


Based upon previous results the domain was reduced to a pipe with a length of 4 meters. This is sufficient to full develop the velocity profile and reduces computational expense. snappyHexMesh is no speed demon when one does not know what they are doing.


The snappyHexMesh dictionary was taken from some other tutorial. There are a few decent references out there including: http://www.hydroniumion.de/general/snappyhexmesh-tutorial/


My initial snappyHexMesh dictionary and the STL file: File:SnappyHexMeshDict.txt File:GeoSTL half.stl.txt


For both files the .txt extension must be removed.


With the previous snappyHexMesh dictionary the resulting mesh looks like:

SnappeHexMeshNumeroUno.png


The previous snappyHexMeshdict dictionary does not terminate normally with OpenFOAM v2.11. The following dictionary was revised to include missing keywords. The filename must be altered to remove the '2.txt'


File:SnappyHexMeshDict2.txt


The following figure shows the mesh in relation to the STL geometry to which the mesh was snapped. To view this particular mesh the time was advanced in ParaView using the 'Next frame' button.

SnappyHexMeshSnapped.png


snappyHexMesh apparently creates a few different meshes (2 or 3?), sort of intermediate steps. It creates each of these meshes at different times. Upon running snappyHexMesh one should notice the creation of a few folders with names that will be associated with the time step specified in the control dictionary. To run a simulation I have been trying to change the start time to the time associated with the final mesh generated by snappyHexMesh. However, this leads to a few problems. The final mesh generated by snappyHexMesh has no boundary conditions associated with it. I would have hoped snappyHexMesh would automatically create new boundary conditions based upon the block mesh boundary conditions I specified. This is not a big problem, I just copy the files 'p' and 'U' from folder '0'. At this point I found that the wall boundary conditions in the block mesh (which I had named 'fixedWalls', consistent with the cavity tutorial) no longer have associated nodes or cell faces in the mesh. snappyHexMesh removed elements in this region and created a new group of cell faces and/or nodes named 'Geo_solid'.

My solution was to change the cell faces/nodes group 'fixedWalls' in my blockmesh to 'Geo_solid' and, correspondingly, change my boundary conditions to be consistent.

To run icoFoam I change the start time to the time at which snappyHexMesh generated the final mesh.

However, this does not work particularly well. When I run icoFoam the Courant number diverges very quickly:


CourantNumberDivergence.png


Running checkMesh there seems to be a problem. Note that checkMesh is indicating a problem for an intermediate mesh. It is not shown but there is also a problem with the final mesh.


CheckMesh.png


A quick search online seemed to indicate that this issue with checking the mesh is associated with the use of the 'empty' boundary condition type. A simple solution was to change the type from 'empty' to 'symmetryPlane' in the blockmesh dictionary and the initial/boundary conditions for pressure and velocity.


Coarser mesh - OpenFOAM

The previous mesh has very small cells that may require a very small time step. Additionally, a logical next step was to fiddle with the snappyHexMesh dictionary and see what happens.

Decreasing the maxLocalCells and maxGlobalCells parameters each by an order of magnitude yielded the following mesh:


PoorQualityCoarseMesh.png


icoFoam will produce the following reasonable results with a time step of 0.0015 seconds:


ResultsCoarseMesh.png


I wanted to try to obtain a nicer looking mesh and found a very good resource: https://sites.google.com/site/snappywiki/snappyhexmesh/cylinder-case


Changing level (1 2) to (0 0) under refinementSurfaces yielded the following mesh. This keyword, refinementSurfaces, controls cell splitting by proximity to a specified surface. The integers in brackets indicate minimum and maximum refinement (see section 5.4.3 of the OpenFOAM users manual).


Coarse hex mesh.png

Skin friction coefficient - OpenFOAM

Another result in the Fluent tutorial, in addition to the axial velocity profile, is the skin friction coefficient:

Skin friction formula.png


OpenFOAM has a tool akin to foamCalc for calculating shear stress: wallShearStress In this case I used wallShearStress -latestTime to calculate shear stress at the termination time when a steady state condition had been reached. However, in order to even use this tool I had to modify the model. I believe there may be more than one way to approach this but this was my solution.

In the folder named 'constant' I opened the transportProperties dictionary and replaced what was there with:

transportModel Newtonian;

rho rho [ 1 -3 0 0 0 0 0 ] 1;

nu nu [ 0 2 -1 0 0 0 0 ] 2e-03;


These properties of the fluid were made consistent with the Fluent tutorial. Two other corrresponding additions to the model were the RASProperties and turbulenceProperties dictionaries also in the folder named constant. Both files may not be necessary. The RASProperties dictionary contains the following lines (see table 3.9 of the OpenFOAM users manual).

RASModel laminar;

turbulence off;

printCoeffs off;


The turbulenceProperties dictionary contains:

simulationType laminar;

After re-running the model with the new dictionaries 'wallShearStress -latestTime' can be run in the OpenFOAM terminal. Opening the results in Paraview (command: 'paraFoam'), jumping to the last time step, and plotting the wall shear stress over a line through (0.0707 0.0707 0.1) and (0.0707 0.0707 3.9) yields the following plot:

WallShearStressPlot.png

The results are a little different than those obtained with Fluent here: https://confluence.cornell.edu/download/attachments/118475269/SkinFrictionPlot_Full.png

However, the magnitude of the shear stress is of course half of the skin friction coefficient due to the one-half in the denominator and the value of unity for the density and velocity. If the points for the line plot are near the planes of symmetry the wall shear stress is under-predicted. The shear stress should just be the viscosity multiplied by the derivative of velocity with respect to distance from the wall assuming perfectly axial flow. Looking at velocity contour plots (magnitude and components) the only significant component is axial and the velocity is relatively constant at a specific axial position and radius (i.e. no radial flow). Perhaps the shear discrepancy with choice of line plot points is the result of how the numerical derivative is performed and the location? With no data outside of the mesh numerical derivatives at boundaries might use null values when unknown or switch to a different finite difference scheme.