waves: Added waves library and setWaves utility

This addition allows for theoretical wave models to be utilised for
initialisation and as boundary conditions. Multiple models can be used
simultaneously, each with differing phases and orientations. If multiple
models are used the shapes and velocities are superimposed.

The wave models are specified in the velocity boundary condition. The
phase fraction boundary condition and the set utility both look up the
velocity condition in order to access the wave model. A velocity
boundary may be specified as follows:

    inlet
    {
        type            waveVelocity;
        origin          (0 0 0);
        direction       (1 0 0);
        speed           2;
        waves
        (
            Airy
            {
                length      300;
                amplitude   2.5;
                depth       150;
                phase       0;
                angle       0;
            }
        );
        scale           table ((1200 1) (1800 0));
        crossScale      constant 1;
    }

The alpha boundary only requires the type, unless the name of the
velocity field is non-standard, in which case a "U" entry will also be
needed. The setWaves utility does not require a dictionary file; non-
standard field names can be specified as command-line arguments.

Wave models currently available are Airy (1st order) and Stokes2 (second
order). If a depth is specified, and it is not too large, then shallow
terms will be included, otherwise the models assume that the liquid is
deep.

This work was supported by Jan Kaufmann and Jan Oberhagemann at DNV GL.

applications/utilities/preProcessing/setWaves/Make/files

@@ -0,0 +1,3 @@
+setWaves.C
+
+EXE = $(FOAM_APPBIN)/setWaves

applications/utilities/preProcessing/setWaves/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/waves/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lwaves

applications/utilities/preProcessing/setWaves/setWaves.C

@@ -0,0 +1,252 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    setWaves
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "levelSet.H"
+#include "timeSelector.H"
+#include "waveAlphaFvPatchScalarField.H"
+#include "waveVelocityFvPatchVectorField.H"
+#include "waveSuperposition.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void addWaves
+(
+    const waveSuperposition& waves,
+    const bool liquid,
+    volScalarField& alpha,
+    volVectorField& U
+)
+{
+    const scalar t = alpha.time().value();;
+    const fvMesh& mesh = alpha.mesh();
+    const pointMesh& pMesh = pointMesh::New(mesh);
+
+    // Height fields
+    const scalarField heightC(waves.height(t, mesh.cellCentres()));
+    const scalarField heightP(waves.height(t, mesh.points()));
+
+    // Velocity fields
+    const DimensionedField<vector, volMesh>
+        UGasC
+        (
+            IOobject("UGasC", mesh.time().timeName(), mesh),
+            mesh,
+            dimVelocity,
+            waves.UGas(t, mesh.cellCentres())
+        );
+    const DimensionedField<vector, pointMesh>
+        UGasP
+        (
+            IOobject("UGasP", mesh.time().timeName(), mesh),
+            pMesh,
+            dimVelocity,
+            waves.UGas(t, mesh.points())
+        );
+    const DimensionedField<vector, volMesh>
+        ULiquidC
+        (
+            IOobject("ULiquidC", mesh.time().timeName(), mesh),
+            mesh,
+            dimVelocity,
+            waves.ULiquid(t, mesh.cellCentres())
+        );
+    const DimensionedField<vector, pointMesh>
+        ULiquidP
+        (
+            IOobject("ULiquidP", mesh.time().timeName(), mesh),
+            pMesh,
+            dimVelocity,
+            waves.ULiquid(t, mesh.points())
+        );
+
+    // Convert from the level set to volume-averaged fields and sum up
+    alpha.ref() += levelSetFraction(mesh, heightC, heightP, !liquid);
+    U.ref() +=
+        levelSetAverage
+        (
+            mesh,
+            heightC,
+            heightP,
+            UGasC,
+            UGasP,
+            ULiquidC,
+            ULiquidP
+        );
+
+    // Now set the boundary fields
+    forAll(alpha.boundaryField(), patchi)
+    {
+        fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi];
+        fvPatchVectorField& Up = U.boundaryFieldRef()[patchi];
+
+        const fvPatch& patch = alphap.patch();
+
+        // Height fields
+        const scalarField heightF(waves.height(t, patch.Cf()));
+        const scalarField heightP(waves.height(t, patch.patch().localPoints()));
+
+        // Velocity fields
+        const vectorField UGasC(waves.UGas(t, mesh.cellCentres()));
+        const vectorField UGasP(waves.UGas(t, mesh.points()));
+        const vectorField ULiquidC(waves.ULiquid(t, mesh.cellCentres()));
+        const vectorField ULiquidP(waves.ULiquid(t, mesh.points()));
+
+        alphap == alphap + levelSetFraction(patch, heightF, heightP, !liquid);
+        Up == Up
+          + levelSetAverage
+            (
+                patch,
+                heightC,
+                heightP,
+                UGasC,
+                UGasP,
+                ULiquidC,
+                ULiquidP
+            );
+    }
+}
+
+
+int main(int argc, char *argv[])
+{
+    timeSelector::addOptions(false, false);
+
+    Foam::argList::addOption
+    (
+        "U",
+        "name",
+        "name of the velocity field, default is \"U\""
+    );
+
+    Foam::argList::addOption
+    (
+        "alpha",
+        "name",
+        "name of the volume fraction field, default is \"alpha\""
+    );
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    instantList timeDirs = timeSelector::selectIfPresent(runTime, args);
+
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+
+    forAll(timeDirs, timeI)
+    {
+        runTime.setTime(timeDirs[timeI], timeI);
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        mesh.readUpdate();
+
+        // Read the phase fraction and velocity fields
+        volScalarField alpha
+        (
+            IOobject
+            (
+                args.optionFound("alpha") ? args["alpha"] : "alpha",
+                runTime.timeName(),
+                mesh,
+                IOobject::MUST_READ
+            ),
+            mesh
+        );
+        volVectorField U
+        (
+            IOobject
+            (
+                args.optionFound("U") ? args["U"] : "U",
+                runTime.timeName(),
+                mesh,
+                IOobject::MUST_READ
+            ),
+            mesh
+        );
+
+        // Zero the fields
+        alpha = dimensionedScalar("0", alpha.dimensions(), 0);
+        U = dimensionedVector("0", U.dimensions(), vector::zero);
+        forAll(alpha.boundaryField(), patchi)
+        {
+            alpha.boundaryFieldRef()[patchi] == 0;
+            U.boundaryFieldRef()[patchi] == vector::zero;
+        }
+
+        // Loop the patches, looking for wave conditions
+        forAll(alpha.boundaryField(), patchi)
+        {
+            const fvPatchScalarField& alphap = alpha.boundaryField()[patchi];
+            const fvPatchVectorField& Up = U.boundaryField()[patchi];
+
+            const bool isWave = isA<waveAlphaFvPatchScalarField>(alphap);
+
+            if (isA<waveVelocityFvPatchVectorField>(Up) != isWave)
+            {
+                FatalErrorInFunction
+                    << "The alpha condition on patch " << Up.patch().name()
+                    << " is " << alphap.type() << " and the velocity condition"
+                    << " is " << Up.type() << ". Wave boundary conditions must"
+                    << " be set in pairs. If the alpha condition is "
+                    << waveAlphaFvPatchScalarField::typeName
+                    << " then the velocity condition must be "
+                    << waveVelocityFvPatchVectorField::typeName
+                    << " and vice-versa." << exit(FatalError);
+            }
+
+            if (isWave)
+            {
+                Info<< "Adding waves from patch " << Up.patch().name() << endl;
+                addWaves
+                (
+                    refCast<const waveVelocityFvPatchVectorField>(Up).waves(),
+                    refCast<const waveAlphaFvPatchScalarField>(alphap).liquid(),
+                    alpha,
+                    U
+                );
+            }
+        }
+
+        // Output
+        Info<< "Writing " << alpha.name() << nl;
+        alpha.write();
+        Info<< "Writing " << U.name() << nl << endl;
+        U.write();
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/Allwmake

@@ -78,6 +78,7 @@ functionObjects/Allwmake $targetType $*
 wmake $targetType sixDoFRigidBodyMotion
 wmake $targetType rigidBodyDynamics
 wmake $targetType rigidBodyMeshMotion
+wmake $targetType waves
 
 
 #------------------------------------------------------------------------------

src/finiteVolume/cfdTools/general/levelSet/levelSet.C

@@ -39,23 +39,29 @@ Foam::levelSetFraction
     const bool above
 )
 {
-    DimensionedField<scalar, volMesh> sum
+    tmp<DimensionedField<scalar, volMesh>> tResult
     (
-        IOobject
+        new DimensionedField<scalar, volMesh>
         (
-            "levelSetIntegral",
-            mesh.time().timeName(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar("0", dimVolume, 0)
+            IOobject
+            (
+                "levelSetFraction",
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar("0", dimless, 0)
+        )
     );
+    DimensionedField<scalar, volMesh>& result = tResult.ref();
 
-    forAll(sum, cI)
+    forAll(result, cI)
     {
         const List<tetIndices> cellTetIs =
             polyMeshTetDecomposition::cellTetIndices(mesh, cI);
 
+        scalar v = 0, r = 0;
+
         forAll(cellTetIs, cellTetI)
         {
             const tetIndices& tetIs = cellTetIs[cellTetI];
@@ -82,18 +88,22 @@ Foam::levelSetFraction
                     levelP[pIB]
                 };
 
+            v += cut::volumeOp()(tet);
+
             if (above)
             {
-                sum[cI] += tetCut(tet, level, cut::volumeOp(), cut::noOp());
+                r += tetCut(tet, level, cut::volumeOp(), cut::noOp());
             }
             else
             {
-                sum[cI] += tetCut(tet, level, cut::noOp(), cut::volumeOp());
+                r += tetCut(tet, level, cut::noOp(), cut::volumeOp());
             }
         }
+
+        result[cI] = r/v;
     }
 
-    return sum/mesh.V();
+    return tResult;
 }
 
 
@@ -105,12 +115,15 @@ Foam::tmp<Foam::scalarField> Foam::levelSetFraction
     const bool above
 )
 {
-    vectorField sum(patch.size(), vector::zero);
+    tmp<scalarField> tResult(new scalarField(patch.size(), 0));
+    scalarField& result = tResult.ref();
 
-    forAll(sum, fI)
+    forAll(result, fI)
     {
         const face& f = patch.patch().localFaces()[fI];
 
+        vector a = vector::zero, r = vector::zero;
+
         for(label eI = 0; eI < f.size(); ++ eI)
         {
             const edge e = f.faceEdge(eI);
@@ -130,18 +143,22 @@ Foam::tmp<Foam::scalarField> Foam::levelSetFraction
                     levelP[e[1]]
                 };
 
+            a += cut::areaOp()(tri);
+
             if (above)
             {
-                sum[fI] += triCut(tri, level, cut::areaOp(), cut::noOp());
+                r += triCut(tri, level, cut::areaOp(), cut::noOp());
             }
             else
             {
-                sum[fI] += triCut(tri, level, cut::noOp(), cut::areaOp());
+                r += triCut(tri, level, cut::noOp(), cut::areaOp());
             }
         }
+
+        result[fI] = a/magSqr(a) & r;
     }
 
-    return sum & patch.Sf()/sqr(patch.magSf());
+    return tResult;
 }
 
 // ************************************************************************* //