• swak4Foam - Swiss Army Knife for OpenFOAM
  • Compilation and installation
  • Usage Example: Continuation of the spillway tutorial
    • 1: Inserting a surface elevation probe
    • 2: Creating a new and simpler blockMeshDict
    • 3: Modifying the boundary conditions
    • 4: Coupling the inlet boundary to the surface elevation probe
    • 5: Comparison between swak4Foam surface probe and Paraview integration
    • 6: Improving the velocity boundary condition
  • Download
  • Further information

swak4Foam - Swiss Army Knife for OpenFOAM

The third-party library swak4Foam is a useful tool witch enables you to do a lot of different tasks with OpenFOAM that otherwise require you to edit and recompile solvers, boundary conditions or other parts of the official OpenFOAM package.

Warning

As most of the rest of the internet, swak4Foam is provided as-is, and you should verify the functions and correctness of the results yourself! The tool is probably not designed with your usage in mind, and therefore, a lot of undiscovered bugs might be present. Please consult the swak4Foam wiki page for the latest announcements, updates and hints.

Compilation and installation

Before you try to compile and install swak4Foam, you should have a valid OpenFOAM installation with MPI. If you are using the HPC computers ve or vilje, this is already present with the module load openfoam and module load mpt commands. You should also have a decent compiler, on ve or vilje it is strongly recommended to use the compilers from Intel. These are loaded with module load intelcomp. On other systems use the recommended compiler for that system, witch often is GCC. This guide applies for OpenFOAM 2.x.x versions.

1. (if you are using a computer other than ve or vilje, for example your personal workstation, skip this step) Load the necessary modules:

   Terminal window module load openfoam module load mpt module load intelcomp

2. Go to the directory where the OpenFOAM build system expect to find third party and user added sources:

   Terminal window cd $WM_PROJECT_USER_DIR

3. Checkout the newest version of swak4Foam from SVN:

   Terminal window svn checkout svn://svn.code.sf.net/p/openfoam-extend/svn/trunk/Breeder_2.0/libraries/swak4Foam If OpenFOAM version newer than 2.0.1 use: svn checkout svn://svn.code.sf.net/p/openfoam-extend/svn/trunk/Breeder_2.0/libraries/swak4Foam/ swak4Foam_2.x

   This repository contains the newest stable version, and is not the development repository.
4. Enter the swak4Foam directory:

   Terminal window cd swak4Foam

5. Compile the sources:

   Terminal window nice ./Allwmake

   The nice command gives the compilation lower priority on the system (compared to for example ssh). This is good practice on shared machines as ve or vilje, since many people use the same login node at the same time. If you compile it on your personal computer using nice it is not important.

If no errors are present after the compilation, the installation is a success.

Usage Example: Continuation of the spillway tutorial

To demonstrate some of the features of swak4Foam we will extend and improve the spillway tutorial using swak4Foam.

The goal by using swak4Foam in this case is to:

1. Automatically measure the free surface elevation at runtime without the need to store and postprocess all the timesteps.
2. Being able of using a simpler blockMeshDict file consisting of only one block
3. Defining one inlet that is partially wet
4. Adjusting the inlet flow velocity and the volume fractionat the inlet to follow the evolution of the free surface, while managing a constant volume flow of water.

1: Inserting a surface elevation probe

First we need to insert a surface elevation probe to calculate the surface elevation while running OpenFOAM. This is done by loading some additional libraries from the swak4Foam package, and searching for

. To load the necessary libraries and define the probe functions, add the following lines to system/controlDict:

libs (
      "libOpenFOAM.so"  // OpenFOAM will (probably) crash if this library is not specified
      "libsimpleSwakFunctionObjects.so"
      "libswakFunctionObjects.so"
      "libgroovyBC.so"  // Not needed now, but we will need it later when altering the boundary conditions
     );
           
functions
{
    #include "sampledSets"
}

Then create a new file system/sampledSets containing:

fillHeight
{
    type swakExpression;
    valueType set;
    verbose true;
    setName surface;
                 
    set
    {
        type    uniform;
        axis    z;
        start   ( -10 0 0 );
        end     ( -10 0 8 );
        nPoints 200;
    }
                 
    expression "(alpha1 > 0.5) ? pos().z : 0";
    accumulations (max);
                 
    interpolate true;
    interpolationType cellPoint ;
}

This set is defined as a line from (-10, 0, 0) to (-10, 0, 8) with 200 points, and the values of the set are the z-coordinate of the point if

and 0 otherwise. Then it searches for the maximum value in this set, and records this as the free surface elevation. The result is appended to a file swakExpression_fillHeight/0/fillHeight at each timestep. The x ? y : z logical operator used above is a short and (very) useful form of if (x==true) then y else z and will be used frequently when dealing with swak4Foam.

Try to run a few seconds of simulation and verify that this works. You can plot the free surface elevation against time using gnuplot or similar.

2: Creating a new and simpler blockMeshDict

One of the goals of this work is to being able to use just one block in the blockMeshDict file, and then handle the "jump" in inlet conditions by groovyBC (part of the swak4Foam package). Simply define one single block of the same dimensions as before:

convertToMeters 1;
 
vertices
(
    (-20 -0.05 0 )
    ( 15 -0.05 0 )
    ( 15 -0.05 8 )
    (-20 -0.05 8 )
    (-20 -0.15 0 )
    ( 15 -0.15 0 )
    ( 15 -0.15 8 )
    (-20 -0.15 8 )
);
                 
blocks
(
    hex (0 1 2 3 4 5 6 7) (350 80 1) simpleGrading (1 1 1)
);
                 
edges
(
);
                 
boundary
(
    inlet
    {
        type patch;
        faces
        (
            (3 0 4 7)
        );
    }
                 
    outlet
    {
        type patch;
        faces
        (
            (2 1 5 6)
        );
    }
                 
    atmosphere
    {
        type patch;
        faces
        (
            (2 3 7 6)
        );
    }
                 
    bottomWall
    {
        type wall;
        faces
        (
            (0 1 5 4)
        );
    }
                 
    front
    {
        type empty;
        faces
        (
            (0 1 2 3)
        );
    }
                 
    back
    {
        type empty;
        faces
        (
            (5 4 7 6)
        );
    }
                 
);
                 
mergePatchPairs
(
);

3: Modifying the boundary conditions

First we need to replace the two previously defined boundaries inletAir and inletWater with a single new boundary inlet. For the turbulent properties k and omega you must make some assumptions like we did in the original spillway tutorial. The pressure conditions for _p_rgh_ will be the same.

The inlet conditions for U and alpha1 will be replaced by a new condition called groovyBC. Let us first make a simple boundary that will behave the same way that our previous inletAir and inletWater did:

inlet
{
    type                groovyBC;
    value               uniform (0.6 0 0); // Initial condition, will only be used until the first time step is complete
 
    variables           "Q=3.6;zSurf=6;";
 
    fractionExpression  "1"; // "1" here indicate that this is a Dirichlet boundary condition (0 is Neumann BC)
    valueExpression     "pos().z<=zSurf ? vector(Q/zSurf, 0, 0) : vector(0, 0, 0)";
}

inlet
{
    type                groovyBC;
    value               uniform 1; // Initail condition, will only be used until the first time step is complete
 
    variables           "zSurf=6;";
 
    fractionExpression  "1"; // "1" here indicate that this is a Dirichlet boundary condition (0 is Neumann BC)
    valueExpression     "pos().z<=zSurf ? 1 : 0";
}

Try to understand what the code above actually does before continuing. The key is the variable definitions (even tough the variables are constant in this case) and the valueExpression part, with the actual expression for the inlet value. Note the usage of the x ? y : z syntax as previously defined.

Try to run the siumulation and verify that you don't get any errors and that the boundary conditions behave as you expect. If you do get errors, please check that the syntax is right. GroovyBC might be a little picky when it comes to extra spaces etc. in the variable definitions and value expressions. You can use Paraview to measure the inlet volume flow and use that as a "validation" criteria.

4: Coupling the inlet boundary to the surface elevation probe

The last topic is to couple the boundary condition at the inlet to the free surface elevation probe that we defined previously. This is not difficult at all. Just change the variable definition of the free surface elevation from zSurf=6; to

zSurf{set'surface}=max(alpha1>0.5 ? pos().z : 0);

in the files with the boundary conditions for the velocity and volume fraction. Note that we are accessing the raw data from the line we defined in the system/sampledSets file, and not the surface elevation. That is recalculated with the max()-function the same way as previously explained.

Run the simulation and compare the result to the result you obtained previously!

The image above shows a comparison after 45 seconds of simulation time with the fixed inlet height on top and the dynamic boundary condition using swak and groovyBC on the bottom. You can clearly see that the streamlines in the upper image is somewhat distorted next to the inlet, and never really flatten out before they "bend over" the spillway. In the lower image the streamlines are almost perfectly horizontal when coming out from the inlet. The positive side of this is that it reduces the time it takes before the surface elevation converges, as the surface now behaves much more "smoothly" with less waves.

5: Comparison between swak4Foam surface probe and Paraview integration

Another version of finding the surface elevation was presented in the spillway tutorial when we integrated the volume fraction. The method presented here is a bit different, since we search for the point where

. But how different? When integrating in Paraview, the results can be stored in a text file, and using Gnuplot or Octave the two datasets can be plotted together:

As we can see, it seems like the "searching algorithm" presented here will be a lower threshold of the integration method. If we look at the data from the "searching algorithm" we can see that the difference between two "jumps" is almost exactly equal to 0.040 meters, witch happens to be the distance between the interpolation points along our sample line (8 meters divided by 200 pints = 0.04 meter/point)! Realizing this, we can increase the resolution on our sample line, rerun the simulation and get smoother data.

6: Improving the velocity boundary condition

If you study the plot above you see that there are clearly some oscillations or resonance in the free surface elevation. The cause of this is that the surface elevation sample set and the inlet does not match 1:1 with respect to cell numbers or their height. The logic implemented in the boundary condition is then not correct. To illustrate the error, imagine that the value of zSurf at a timestep is 6.14 meters, and that the cell height is 0.1 meter at the inlet. The face centers of the boundary will then be 0.05, 0.15, ..., 6.05, 6.15, 6.25, ... meters above the channel bottom. The result will be that the face with center at z=6.05 will become a water inlet, but the face with center at z=6.15 will not. The velocity at the inlet will however still be calculated based on the zSurf value of 6.14 meters: Ux = 3.6/6.14 = 0.586 m/s, given a target volume flow of 3.6 m^3/s. The real volume flow will however be the velocity calculated multiplied over the real inlet height: Qreal = 0.586 * 6.10 = 3.57 m^3/s. As the surface elevation changes, the volume flow changes due to discretization errors!

How to avoid it? The simple and easy solution is to make sure that the velocity calculated always will give the correct volume flow, even if it cannot be applied over the correct height. Adding another equation to the velocity inlet does the trick:

inlet
                {
                    type                groovyBC;
                    value               uniform (0.6 0 0); // Initial condition, will only be used until the first time step is complete
                 
                    variables           "Q=3.6;zSurf{set'surface}=max(alpha1>0.5 ? pos().z : 0);zInlet=max(pos().z<=zSurf ? pos().z : 0);";
                 
                    fractionExpression  "1"; // "1" here indicate that this is a Dirichlet boundary condition (0 is Neumann BC)
                    valueExpression     "pos().z<=zInlet ? vector(Q/zInlet, 0, 0) : vector(0, 0, 0)";
                }

What the code actually does now is that is finds the real inlet height zInlet first, and then calculates the inlet velocity based zInlet instead of the surface elevation zSurf. This gives a constant volume flow independent on the discretization of the inlet and sample line (but is it correct and exactly equal to 3.6 m^3/s? Think about that.).

The image above shows the volume flow through the inlet plotted against time. As we can see, using a high resolution line source together with the improved boundary condition logic makes the boundary a lot more stable. The final value of the flow into the domain in the last simulation is 3.5989 m^3/s, witch is just 0.03 percent off the target value of 3.6 m^3/s. When looking at these figures one must also take into account that this is a result of postprocessing, and the the writePrecision in this case is 6 digits (i.e. OpenFOAM only writes 6 significant digits to disk for prostprocessing). Increasing the write precision might have revealed that we are even closer to the target, but that is not investigated here.

As we can see from the plot above, with a surface elevation probe with 800 points and the improved boundary condition, the elevation signal is now also much more continuous. However, there is still some oscillations of high frequency on the surface. The reason for this is not properly investigated, however it might be waves due to the air-water interface:

Studying details in the flow velocity (pink line) around the surface shows us a significant velocity gradient between air and water. So wave-making due to the air-water interface cannot be ruled out. For all practical purposes however, this is not an issue and we will not deal with this any more. Perhaps you can experiment with some different boundary conditions for the non-water part of the inlet and the atmosphere to get rid of it?

Download

You can download all necessary case files in swak4Foam.tar.gz. Scripts to run and clean the case is also attached.

Further information

If you need more information about swak4Foam, please consult its wiki page or any of the wiki pages for groovyBC or funkysetFields. A good presentation is given at the 6th OpenFOAM workshop in 2011, where both slides and example case files can be downloaded.