Documentation

In this article a brief introduction is given on how to implement a Numerical Wave Tank (NWT) in OpenFOAM to simulate regular water waves. NWTs are used to simulate regular and irregular waves to study water free surface dynamics, sediment transport, ﬂuid-structure interaction due to wave loads and more.

In this case, a regular linear airy wave is simulated in the created NWT. The NWT’s dimension are (L = 10, W = 1.2 and H=0.6) m. The simulated wave height is 2 cm and the wave length is 3 m.

For ease of use, this numerical simulation is implemented via the blueCFD software (a Windows version of OpenFOAM). This software can be downloaded and installed as an executable ﬁle from http://bluecfd.github.io/Core/Downloads/. After installation is completed, the user can use the blueCFD core terminal to execute OpenFOAM commands.

OpenFOAM is a Finite Volume based open source code that includes several solvers dedicated to solve various ﬂuid ﬂow problems. The choice between these solvers depends on the user application of interest.

For this NWT application we will use the interFOAM solver (Damian (2012)) that is able to solve for multiphase ﬂows using the Volume Of Fluid technique for ﬂuid-ﬂuid interface capturing. The solver is incompressible and isothermal. interFoam solves for the coupled PDE’s, continuity equation 1, Navier-Stoke’s equations 2 and the phase fraction equation 4). Where $u$ is the ﬂow velocity. $p$ and $ρ$ are the thermodynamic pressure and mixture density (5). $T$ is the viscous stress tensor that is a function of the eﬀective dynamic viscosity and the general rate of strain tensor $D$ . $fb$ represents body forces, including the gravitational force and the surface tension. The solver uses the PIMPLE algorithm for velocity-pressure coupling (Holzmann (2016)).

∇ ⋅u = 0

(1)

∂(ρu) +∇ ⋅(ρuu ) = − ∇p + ∇ ⋅T + ρfb ∂t

(2)

T = 2μD

(3)

∂α-+ ∇ ⋅αu = 0 ∂t

(4)

ρ = αρw + (1− α )ρa

(5)

OpenFOAM simulation cases should have a certain folder/ﬁle structure.

Figure 1 shows the folder/ﬁle structure for the NWT case. folder (0) contains the ﬁles that represents the boundary and initial conditions for the velocity ﬁeld, pressure ﬁeld and the phase fraction ﬁeld. Folder (system) contains the ﬁles that control the numerical schemes such as fvSchemes and fvSolution. It also includes the controlDict ﬁle where the user can indicate the solver of interest such as (interFoam) and the simulation time in addition to other control parameters. The system folder includes the blockMeshDict ﬁle that is used to generate the numerical grid for the simulation case.

The constant folder includes additional information for the solver, not related to boundary condition or numerical schemes. The ﬁles in the constant folder can include a g ﬁle that determines the direction of the gravitational force, (transportProperties) includes the material properties such as viscosity and density of the ﬂuid and ﬁle waveProperties contains the simulated wave properties such as wave amplitude and wave length. Figure 2 shows the the blueCFD terminal on the root folder location and the sub folders/ﬁles. (The blueCFD terminal can be opened on the root folder location by right-clicking on the root folder and choose ”Open with blueCFD core terminal”).Note: The folder structure for the example presented in this article is attached as NWT.zip.

Figure 1: Folder structure for an interFoam solver case

Figure 2: blueCFD terminal showing the main folders and ﬁles in the root folder

To generate a simple computational domain for the NWT with cuboid shape, the blockMesh utility is used. With blockMesh, the user can assign Cartesian coordinates for all 8 vertices in the cuboid shape as well as deﬁning certain boundary names and types for faces constructed from these vertices. The numbering system for the blockMesh utility is shown in ﬁgure 3.

Figure 3: blockMesh cuboid vertex numbering, This image is adopted form (Greenshields (2018))

Code 1 shows the blockMesh utility dictionary ﬁle structure. First, the min/max dimension variables for the NWT are listed, then the vertices Cartesian coordinates $(x1,x2,x3)$ are deﬁned based on the numbering system in ﬁgure 3. Then the cuboid shape is constructed by the hex object with the numbered vertices followed by the number of elements in $(x1,x2,x3)$ directions, here we used 45,15 and 45 elements in $x1$ , $x2$ and $x3$ directions respectively.

The boundary object deﬁnes the name, type and location of each boundary on the computational domain outer surface. In this case we deﬁne ﬁve boundaries for the NWT (Inlet, Outlet, sides, bottom and top). Each boundary is constructed by four vertices as shown in in ﬁgure 3.

Typing blockMesh in the blueCFD terminal makes OpenFOAM generate the computational domain (the higher number of elements the longer time it takes to construct the grid).

Now, the user is ready to type paraFoam in the terminal to show the grid in the ParaView post processing software as shown in ﬁgure 4

Code 1: blockMeshDict to generate a simple grid for a cuboid domain

  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*--------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      blockMeshDict;
}

// * * * * * * * * * * * * * * * * * *
convertToMeters 1;
//  Tank dimensions
X_MIN   0  ;
X_MAX   10  ;
Y_MIN   -.6  ;
Y_MAX  .6  ;
Z_MIN   -0.3 ;
Z_MAX   .3  ;
// Vertices locations for the Tank
vertices
(
    ($X_MIN $Y_MIN $Z_MIN)
    ($X_MAX $Y_MIN $Z_MIN)
    ($X_MAX $Y_MAX $Z_MIN)
    ($X_MIN $Y_MAX $Z_MIN)
    ($X_MIN $Y_MIN $Z_MAX)
    ($X_MAX $Y_MIN $Z_MAX)
    ($X_MAX $Y_MAX $Z_MAX)
    ($X_MIN $Y_MAX $Z_MAX)
);
// Build the block Mesh based on the indices numbering
blocks
(hex (0 1 2 3 4 5 6 7) (45 15 45)  simpleGrading (1 1 1));
edges
(
);
defaultPatch
{
    name frontAndBack;
    type empty;}
// Give a name for each side
boundary
(
    inlet
    { type patch;
        faces
        ((0 4 7 3) );}
    outlet
    {
        type wall;
        faces
        ((2 6 5 1));}
    bottom
    {
        type wall;
        faces
        ((1 2 3 0)); }
    top
    {
        type patch;
        faces
        ( (7 4 5 6));}
      sides
    {
        type wall;
        faces
        ((4 5 1 0)
            (7 3 2 6));}    );
mergePatchPairs();
// ***********************************

Figure 4: Grid generated by the blockMesh utility

In this section we deﬁne the initial and boundary conditions IC/BC for the interFoam solver when solving the set of equations (1, 2 and 4). The IC/BC are deﬁned for the velocity ﬁeld u in the U ﬁle, pressure ﬁeld p in the p_rgh ﬁle and the phase fraction ﬁeld $α$ in the alpha.water ﬁle.

In addition, we need to deﬁne two more ﬁles , waveProperties where the incident wave parameters are deﬁned and setFields where the IC for the water/air phase fraction $α$ are deﬁned. This will essentially ﬁll the tank with water and deﬁne the mean water level.

The velocity BCs are shown in code 2. For the Inlet boundary, the type is wave velocity which will be deﬁned in code 5. noSlip BC is applied to the tank side, bottom and outlet walls.

A pressureInletOutletVelocity BC is applied to the top wall, this BC is usually used for boundaries that are open to the atmosphere and where inﬂow or outﬂow may take place. When inﬂow occurs the velocity is calculated by interpolating the neighbouring cells, for outﬂow, the pressure gradient is set to zero at that boundary.

The attached test case folder contains the BC called 0.org, it is customary to keep the original folder untouched and to make a copy of it calling it 0 where the phase fraction BC is set.

More information about the diﬀerent OpenFOAM boundary conditions can be found in (Greenshields (2018))

Code 2: U. Velocity boundary condition

  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*-----------------------------------------
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "0";
    object      U;
}
// * * * * * * * * * * * * *  * * * * * * * * * * *
dimensions      [0 1 -1 0 0 0 0];
internalField   uniform (0 0 0);
boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"
    inlet
    {type            waveVelocity;}
    outlet
    {type            noSlip;}
    top
    {type            pressureInletOutletVelocity;
        value           uniform (0 0 0); }
    bottom
    {type            noSlip;}
     sides
    {type            noSlip; }}
// **************************************************** //

Code 3 shows the pressure ﬁeld BC, for all boundaries zero gradient is enforced except for the top boundary where a total pressure BC is applied that represents atmospheric datum pressure.

Code 3: p_rgh. Pressure boundary condition

  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*------------------------------------------
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      p_rgh;}
// * * * * * * * * * * * * * * * * * *
dimensions      [1 -1 -2 0 0 0 0];
internalField   uniform 0;
boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"
    "(inlet|outlet|bottom|sides)"
    {
        type            fixedFluxPressure;
        value           uniform 0;}
    top
    {
        type            totalPressure;
        p0              uniform 0;}}
// ***************//

Code 4 shows the alpha.water BC ﬁle for the phase fraction $α$ . For the outlet, sides and bottom boundaries a zeroGradient BC is used. For the inlet boundary, waveAlpha BC is used that is linked to the waveProperties ﬁle discussed later in code 5. inletOutlet BC is applied for the top boundary, this is a mixed boundary condition which applies zeroGradient for outﬂow condition and user speciﬁed value for inﬂow condition.

Code 4: alpha.water. Phase fraction boundary condition

  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*--------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      alpha.water;
}
// * * * * *  * * * * * * * * * * * * * * * * //
dimensions      [0 0 0 0 0 0 0];
internalField   uniform 0;
boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"
    "(outlet|bottom|sides)"
    {
        type            zeroGradient;}
    inlet
    {
        type            waveAlpha;
        U               U;
        inletOutlet     true;}
    top
    {
        type            inletOutlet;
        inletValue      uniform 0;
        value           uniform 0; }}
// ***************************************** //

In code 5 the wave properties are deﬁned by the wave amplitude, length and the water current speed UMean. In addition, the user can deﬁne the wave direction (angle). The wave origin and direction are speciﬁed with respect to the computational domain coordinate system generated by the blockMesh utility. Here we are simulating an (Airy) wave that obeys linear wave theory. Other higher order nonlinear wave models are also included in the OpenFOAM wave library such as Stoke’s waves.

At the end of the waveProperties ﬁle, scale object is deﬁned to linearly damp the generated waves in the range 7-9 m of the tank length in the wave propagation direction.

Code 5: waveProperties. Wave properties deﬁnition

  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*-----------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      waveProperties;
}
// * * * * * *  * * * * * * * * * * * * * * * //
origin          (0 0 0);
direction       (1 0 0);
waves
(
    Airy
    //Stokes2
    {
        length      3;
        amplitude   0.01;
        phase       0;
        angle       0; });
UMean           (0 0 0);   // current speed
scale           table ((7 1) (9 0));
crossScale      constant 1;
// ************************************* //

Code 6 shows the IC for the phase fraction α in the computational domain. First, the whole domain α has value of zero that represents air. Then using the boxToCell object, all cells in the in the box dimension of interest can be assigned new value of α = 1 that represent water. The z_max value in the box object is set to zero which puts the undisturbed mean water free surface level at zero on the z axis.

The user can run the setFields command in the blueCFD terminal by typing setFields. To view the α ﬁeld run paraFoam in blueCFD terminal to open ParaView as shown in ﬁgure 5.

Code 6: setFileds. Deﬁning the initial mean water level in the NWT

  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*-----------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      setFieldsDict;
}
// * * * * * *  * * * * * * * * * * * * //
defaultFieldValues
( volScalarFieldValue alpha.water 0
);
regions
(
    // Set cell values
    // (does zerogradient on boundaries)
    boxToCell
    {  box (-999 -999 -999) (999 999 0);
        fieldValues
        (
            volScalarFieldValue alpha.water 1
        );});
// ******************************* //

Figure 5: Undisturbed phase fraction α. α=1 is water, α = 0 is air

The user can control the major parameters of the case using the controlDict. These parameters include solver type (interFoam), simulation endtime, data saving and writing intervals and the maximum allowable Courant number for simulation stability control. Also in controlDict, the user can call the libraries used in the simulation such as the wave generation library libwaves.so. In addition, other post post processing and, output controls and data logging at certain locations or surfaces can be achieved. In this example we used the interfaceHeight function to monitor the water free surface elevation at certain location in the NWT.

Other simulation control ﬁles are placed in the system folder such as fvSchemes and fvSoluiton, these two ﬁles can control the numerical schemes such as the spatio-temporal numerical diﬀerentiation in the interFoam solver or the parameters of the velocity-pressure coupling algorithm.

To start the interFoam solver, run interFoam command in the blueCFD terminal. The user will notice new folders are being created in the root folder while the simulation is running. These have names corresponding to the time stamps being written. A postProcessing folder is also generated which includes the interfaceHeight output containing the wave elevation at the wave probe point.

Code 7: ControlDict.

  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*----------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * ** * * * * * * * * * * * * * * //
application     interFoam;
startFrom       latestTime;
startTime       0;
stopAt          endTime;
endTime         80;
deltaT          0.01;
writeControl    adjustableRunTime;
writeInterval   0.1;
purgeWrite      0;
writeFormat     ascii;
writePrecision  6;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable yes;
adjustTimeStep  yes;
maxCo           1.5;
maxAlphaCo      1.5;
maxDeltaT       0.1;
libs
( "libwaves.so");

functions
{
    interfaceHeight1
    { type            interfaceHeight;
        libs            ("libfieldFunctionObjects.so");
        locations       (( 5 0 -1) );
        alpha           alpha.water;
    }}
// ********************************************** //

After the solver ﬁnishes, the user can view the results in ParView by typing paraFoam in the blueCFD terminal. Figure 6 shows the pressure contour for the water free surface at isosurface α = 0.5.

Figure 7 shows the free surface elevation at the virtual wave probe point (The output of the interfaceHeight function deﬁned in the controlDict).

Figure 6: Pressure contour at the water free surface

Figure 7: Free surface water elevation at the virtual wave probe point

The user can automate the case running process by running the Allrun script in the blueCFD terminal as shown in code 8. The Allrun script ﬁrst makes a copied directory for the initial condition 0 folder, then runs (blockMesh, setFields,interFoam) in order. To run the Allrun script, simply type ./Allrun in the blueCFD terminal. To delete all results and data logging of the simulation in an eﬃcient and easy way, the user can run the Allclean script shown in code 9.

Code 8: Allrun script

. $WM_PROJECT_DIR/bin/tools/RunFunctions
cp -r 0.org 0 # To copy the 0.org folder to 0 folder
runApplication blockMesh
runApplication setFields
runApplication interfoam

Code 9: Allclean script

# Source tutorial clean functions
. $WM_PROJECT_DIR/bin/tools/CleanFunctions
# Remove surface
rm -f constant/triSurface/DTC-scaled.stl.gz > /dev/null 2>&1
rm -rf constant/extendedFeatureEdgeMesh > /dev/null 2>&1
rm -f constant/triSurface/DTC-scaled.eMesh > /dev/null 2>&1
cleanCase
rm -r 0

Damian, S. M. (2012). Description and utilization of interfoam multiphase solver. International Center for Computational Methods in Engineering, pages 1–64.

Greenshields, C. (2018). Openfoam–the openfoam foundation–user guide.

Holzmann, T. (2016). Mathematics, numerics, derivations and openfoam®. Loeben, Germany: Holzmann CFD.

Creating a Numerical wave tank on OpenFOAM BlueCFD

1 Introduction and problem statement

2 Software installation

3 interfoam solver

4 OpenFOAM case folder structure

5 Mesh creation

6 Boundary and initial conditions

7 OpenFOAM simulation control

8 Results and post processing

9 Allrun/Allclean scripts

References