waves: Moved mean velocity specification back into the wave models

With the inclusion of boundary layer modelling in the gas, the
separation of wave perturbation from and mean flow became less useful,
and potentially prevents further extension to support similar boundary
layer modelling in the liquid.

The mean velocity entry, UMean, is now needed in the
constant/waveProperties file rather than in the waveVelocity boundary
condition.

applications/utilities/preProcessing/setWaves/setWaves.C

@@ -91,12 +91,6 @@ int main(int argc, char *argv[])
     const word alphaName = setWavesDict.lookupOrDefault<word>("alpha", "alpha");
     const word UName = setWavesDict.lookupOrDefault<word>("U", "U");
     const bool liquid = setWavesDict.lookupOrDefault<bool>("liquid", true);
-    const dimensionedVector UMean
-    (
-        "UMean",
-        dimVelocity,
-        setWavesDict.lookup("UMean")
-    );
 
     // Get the wave models
     const waveSuperposition& waves = waveSuperposition::New(mesh);
@@ -172,33 +166,30 @@ int main(int argc, char *argv[])
             dimensionedVector("0", dimVelocity, vector::zero)
         );
 
-        // Offset
-        const vector offset = UMean.value()*t;
-
         // Cell centres and points
         const pointField& ccs = mesh.cellCentres();
         const pointField& pts = mesh.points();
 
         // Internal field
-        h.primitiveFieldRef() = waves.height(t, ccs + offset);
+        h.primitiveFieldRef() = waves.height(t, ccs);
-        hp.primitiveFieldRef() = waves.height(t, pts + offset);
+        hp.primitiveFieldRef() = waves.height(t, pts);
-        uGas.primitiveFieldRef() = waves.UGas(t, ccs + offset);
+        uGas.primitiveFieldRef() = waves.UGas(t, ccs);
-        uGasp.primitiveFieldRef() = waves.UGas(t, pts + offset);
+        uGasp.primitiveFieldRef() = waves.UGas(t, pts);
-        uLiq.primitiveFieldRef() = waves.ULiquid(t, ccs + offset);
+        uLiq.primitiveFieldRef() = waves.ULiquid(t, ccs);
-        uLiqp.primitiveFieldRef() = waves.ULiquid(t, pts + offset);
+        uLiqp.primitiveFieldRef() = waves.ULiquid(t, pts);
 
         // Boundary fields
         forAll(mesh.boundary(), patchj)
         {
             const pointField& fcs = mesh.boundary()[patchj].Cf();
-            h.boundaryFieldRef()[patchj] = waves.height(t, fcs + offset);
+            h.boundaryFieldRef()[patchj] = waves.height(t, fcs);
-            uGas.boundaryFieldRef()[patchj] = waves.UGas(t, fcs + offset);
+            uGas.boundaryFieldRef()[patchj] = waves.UGas(t, fcs);
-            uLiq.boundaryFieldRef()[patchj] = waves.ULiquid(t, fcs + offset);
+            uLiq.boundaryFieldRef()[patchj] = waves.ULiquid(t, fcs);
         }
 
         // Set the fields
         alpha == levelSetFraction(h, hp, !liquid);
-        U == UMean + levelSetAverage(h, hp, uGas, uGasp, uLiq, uLiqp);
+        U == levelSetAverage(h, hp, uGas, uGasp, uLiq, uLiqp);
 
         // Set the boundary fields
         forAll(mesh.boundary(), patchi)

etc/caseDicts/annotated/setWavesDict

@@ -24,8 +24,5 @@ U       U;
 // under the waves (true) or the gas over the waves (false)
 liquid  true;
 
-// The mean flow velocity over which to superimpose waves
-UMean   (1 0 0);
-
 
 // ************************************************************************* //

src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C

@@ -123,18 +123,13 @@ Foam::tmp<Foam::scalarField> Foam::waveAlphaFvPatchScalarField::alpha() const
 {
     const scalar t = db().time().timeOutputValue();
     const waveSuperposition& waves = waveSuperposition::New(db());
-    const waveVelocityFvPatchVectorField& Up =
-        refCast<const waveVelocityFvPatchVectorField>
-        (
-            patch().lookupPatchField<volVectorField, scalar>(UName_)
-        );
 
     return
         levelSetFraction
         (
             patch(),
-            waves.height(t, patch().Cf() + Up.offset()),
+            waves.height(t, patch().Cf()),
-            waves.height(t, patch().patch().localPoints() + Up.offset()),
+            waves.height(t, patch().patch().localPoints()),
             !liquid_
         );
 }
@@ -159,8 +154,8 @@ Foam::tmp<Foam::scalarField> Foam::waveAlphaFvPatchScalarField::alphan() const
         levelSetFraction
         (
             meshs,
-            waves.height(t, meshs.cellCentres() + Up.offset())(),
+            waves.height(t, meshs.cellCentres())(),
-            waves.height(t, meshs.points() + Up.offset())(),
+            waves.height(t, meshs.points())(),
             !liquid_
         )
     );

src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C

@@ -117,22 +117,17 @@ Foam::tmp<Foam::scalarField> Foam::wavePressureFvPatchScalarField::p() const
 {
     const scalar t = db().time().timeOutputValue();
     const waveSuperposition& waves = waveSuperposition::New(db());
-    const waveVelocityFvPatchVectorField& Up =
-        refCast<const waveVelocityFvPatchVectorField>
-        (
-            patch().lookupPatchField<volVectorField, scalar>(UName_)
-        );
 
     return
         levelSetAverage
         (
             patch(),
-            waves.height(t, patch().Cf() + Up.offset()),
+            waves.height(t, patch().Cf()),
-            waves.height(t, patch().patch().localPoints() + Up.offset()),
+            waves.height(t, patch().patch().localPoints()),
-            waves.pGas(t, patch().Cf() + Up.offset())(),
+            waves.pGas(t, patch().Cf())(),
-            waves.pGas(t, patch().patch().localPoints() + Up.offset())(),
+            waves.pGas(t, patch().patch().localPoints())(),
-            waves.pLiquid(t, patch().Cf() + Up.offset())(),
+            waves.pLiquid(t, patch().Cf())(),
-            waves.pLiquid(t, patch().patch().localPoints() + Up.offset())()
+            waves.pLiquid(t, patch().patch().localPoints())()
         );
 }
 
@@ -156,12 +151,12 @@ Foam::tmp<Foam::scalarField> Foam::wavePressureFvPatchScalarField::pn() const
         levelSetAverage
         (
             meshs,
-            waves.height(t, meshs.cellCentres() + Up.offset())(),
+            waves.height(t, meshs.cellCentres())(),
-            waves.height(t, meshs.points() + Up.offset())(),
+            waves.height(t, meshs.points())(),
-            waves.pGas(t, meshs.cellCentres() + Up.offset())(),
+            waves.pGas(t, meshs.cellCentres())(),
-            waves.pGas(t, meshs.points() + Up.offset())(),
+            waves.pGas(t, meshs.points())(),
-            waves.pLiquid(t, meshs.cellCentres() + Up.offset())(),
+            waves.pLiquid(t, meshs.cellCentres())(),
-            waves.pLiquid(t, meshs.points() + Up.offset())()
+            waves.pLiquid(t, meshs.points())()
         )
     );
 

src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C

@@ -42,7 +42,6 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
     phiName_("phi"),
     pName_("p"),
     inletOutlet_(true),
-    UMean_(nullptr),
     faceCellSubset_(nullptr),
     faceCellSubsetTimeIndex_(-1)
 {
@@ -63,7 +62,6 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
     pName_(dict.lookupOrDefault<word>("p", "p")),
     inletOutlet_(dict.lookupOrDefault<Switch>("inletOutlet", true)),
-    UMean_(Function1<vector>::New("UMean", dict)),
     faceCellSubset_(nullptr),
     faceCellSubsetTimeIndex_(-1)
 {
@@ -94,7 +92,6 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
     phiName_(ptf.phiName_),
     pName_(ptf.pName_),
     inletOutlet_(ptf.inletOutlet_),
-    UMean_(ptf.UMean_, false),
     faceCellSubset_(nullptr),
     faceCellSubsetTimeIndex_(-1)
 {}
@@ -109,7 +106,6 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
     phiName_(ptf.phiName_),
     pName_(ptf.pName_),
     inletOutlet_(ptf.inletOutlet_),
-    UMean_(ptf.UMean_, false),
     faceCellSubset_(nullptr),
     faceCellSubsetTimeIndex_(-1)
 {}
@@ -125,7 +121,6 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
     phiName_(ptf.phiName_),
     pName_(ptf.pName_),
     inletOutlet_(ptf.inletOutlet_),
-    UMean_(ptf.UMean_, false),
     faceCellSubset_(nullptr),
     faceCellSubsetTimeIndex_(-1)
 {}
@@ -165,16 +160,15 @@ Foam::tmp<Foam::vectorField> Foam::waveVelocityFvPatchVectorField::U() const
     const waveSuperposition& waves = waveSuperposition::New(db());
 
     return
-        UMean_->value(t)
+        levelSetAverage
-      + levelSetAverage
         (
             patch(),
-            waves.height(t, patch().Cf() + offset()),
+            waves.height(t, patch().Cf()),
-            waves.height(t, patch().patch().localPoints() + offset()),
+            waves.height(t, patch().patch().localPoints()),
-            waves.UGas(t, patch().Cf() + offset())(),
+            waves.UGas(t, patch().Cf())(),
-            waves.UGas(t, patch().patch().localPoints() + offset())(),
+            waves.UGas(t, patch().patch().localPoints())(),
-            waves.ULiquid(t, patch().Cf() + offset())(),
+            waves.ULiquid(t, patch().Cf())(),
-            waves.ULiquid(t, patch().patch().localPoints() + offset())()
+            waves.ULiquid(t, patch().patch().localPoints())()
         );
 }
 
@@ -190,16 +184,15 @@ Foam::tmp<Foam::vectorField> Foam::waveVelocityFvPatchVectorField::Un() const
 
     const vectorField Us
     (
-        UMean_->value(t)
+        levelSetAverage
-      + levelSetAverage
         (
             meshs,
-            waves.height(t, meshs.cellCentres() + offset())(),
+            waves.height(t, meshs.cellCentres())(),
-            waves.height(t, meshs.points() + offset())(),
+            waves.height(t, meshs.points())(),
-            waves.UGas(t, meshs.cellCentres() + offset())(),
+            waves.UGas(t, meshs.cellCentres())(),
-            waves.UGas(t, meshs.points() + offset())(),
+            waves.UGas(t, meshs.points())(),
-            waves.ULiquid(t, meshs.cellCentres() + offset())(),
+            waves.ULiquid(t, meshs.cellCentres())(),
-            waves.ULiquid(t, meshs.points() + offset())()
+            waves.ULiquid(t, meshs.points())()
         )
     );
 
@@ -297,7 +290,6 @@ void Foam::waveVelocityFvPatchVectorField::write
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "p", "p", pName_);
     writeEntryIfDifferent<Switch>(os, "inletOutlet", true, inletOutlet_);
-    UMean_->writeData(os);
 }
 
 

src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H

@@ -62,7 +62,6 @@ Usage
         phi         | Name of the flux field                 | no     | phi
         p           | Name of the pressure field             | no     | p
         inletOutlet | does the condition behave like inletOutlet | no | true
-        UMean       | velocity of the mean flow              | yes    |
         ramp        | ramping function for the mean flow speed | no   | None
     \endtable
 
@@ -74,7 +73,6 @@ Usage
         phi         phi;
         p           p;
         inletOutlet yes;
-        UMean       (2 0 0);
         ramp        constant 1;
     }
     \endverbatim
@@ -119,9 +117,6 @@ class waveVelocityFvPatchVectorField
         //- Act as an inlet/outlet patch?
         const Switch inletOutlet_;
 
-        //- Mean velocity
-        const autoPtr<Function1<vector>> UMean_;
-
         //- Mesh subset corresponding to the patch adjacent cells
         mutable autoPtr<fvMeshSubset> faceCellSubset_;
 
@@ -218,18 +213,6 @@ public:
             //- Return the current modelled velocity field in the neighbour cell
             tmp<vectorField> Un() const;
 
-            //- Return the mean velocity
-            inline vector UMean() const
-            {
-                return UMean_->value(db().time().timeOutputValue());
-            }
-
-            //- Return the distance offset for the wave models
-            inline vector offset() const
-            {
-                return - UMean_->integrate(0, db().time().timeOutputValue());
-            }
-
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 

src/waves/waveSuperpositions/waveAtmBoundaryLayerSuperposition/waveAtmBoundaryLayerSuperposition.C

@@ -77,19 +77,30 @@ Foam::tmp<Foam::vectorField> Foam::waveAtmBoundaryLayerSuperposition::UGas
             db().lookupObject<uniformDimensionedVectorField>("g").value()
         );
 
-    const scalar h0 = - gHat & origin();
+    const scalar h0 = - gHat & origin_;
 
-    atmBoundaryLayer atm
+    const vector UGasRefRel = UGasRef_ - UMean_->value(t);
-    (
-        normalised(UGasRef_),
-      - gHat,
-        mag(UGasRef_),
-        h0 + hRef_,
-        scalarField(p.size(), hWaveMax_ - hWaveMin_),
-        scalarField(p.size(), h0 + hWaveMin_)
-    );
 
-    return waveSuperposition::UGas(t, p) + atm.U(p);
+    const scalar magUGasRefRel = mag(UGasRefRel);
+
+    tmp<vectorField> tU = waveSuperposition::UGas(t, p);
+
+    if (magUGasRefRel > 0)
+    {
+        atmBoundaryLayer atm
+        (
+            UGasRefRel/magUGasRefRel,
+          - gHat,
+            magUGasRefRel,
+            h0 + hRef_,
+            scalarField(p.size(), hWaveMax_ - hWaveMin_),
+            scalarField(p.size(), h0 + hWaveMin_)
+        );
+
+        tU.ref() += atm.U(p);
+    }
+
+    return tU;
 }
 
 

src/waves/waveSuperpositions/waveAtmBoundaryLayerSuperposition/waveAtmBoundaryLayerSuperposition.H

@@ -26,24 +26,23 @@ Class
 
 Description
     An extension of waveSuperposition which adds an atmospheric boundary layer
-    model to the gas velocity. The user supplies a gas velocity relative to the
+    model to the gas velocity. The user supplies a gas velocity and a height
-    mean liquid velocity and a height above the wave coordinate system origin
+    above the wave coordinate system origin at which that velocity is reached.
-    at which that velocity is reached. Also needed are a maximum and minimum
+    Also needed are a maximum and minimum wave height which are used to set the
-    wave height which are used to set the surface roughness in the boundary
+    surface roughness in the boundary layer model. It is not trivial to
-    layer model. It is not trivial to determine these from an arbitrary
+    determine these from an arbitrary superposition of differently oriented
-    superposition of differently oriented wave models, so they are required as
+    wave models, so they are required as user inputs instead. For a pure
-    user inputs instead. For a pure sinusoidal wave, the maximum and minimum
+    sinusoidal wave, the maximum and minimum wave heights can be set to
-    wave heights can be set to positive and negative amplitude, respectively.
+    positive and negative amplitude, respectively.
 
 Usage
     \table
-        Property | Description                             | Req'd? | Default
+        Property | Description                              | Req'd? | Default
-        UGasRef  | The gas velocity relative to the \\
+        UGasRef  | The gas velocity at the reference height | yes    |
-                   liquid at the reference height          | yes    |
         hRef     | The reference height relative to the \\
-                   origin of the wave coordinate system    | yes    |
+                   origin of the wave coordinate system     | yes    |
-        hWaveMin | The minimum wave elevation              | yes    |
+        hWaveMin | The minimum wave elevation               | yes    |
-        hWaveMax | The maximum wave elevation              | yes    |
+        hWaveMax | The maximum wave elevation               | yes    |
     \endtable
 
     Example specification:

src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.C

@@ -48,6 +48,7 @@ namespace Foam
 
 void Foam::waveSuperposition::transformation
 (
+    const scalar t,
     const vectorField& p,
     tensor& axes,
     vectorField& xyz
@@ -64,7 +65,7 @@ void Foam::waveSuperposition::transformation
 
     axes = tensor(dSurfHat, - gHat ^ dSurfHat, - gHat);
 
-    xyz = axes & (p - origin_);
+    xyz = axes & (p - origin_ - UMean_->integrate(0, t));
 }
 
 
@@ -204,6 +205,7 @@ Foam::waveSuperposition::waveSuperposition(const objectRegistry& db)
     direction_(lookup("direction")),
     waveModels_(),
     waveAngles_(),
+    UMean_(Function1<vector>::New("UMean", *this)),
     scale_
     (
         found("scale")
@@ -252,7 +254,7 @@ Foam::tmp<Foam::scalarField> Foam::waveSuperposition::height
 {
     tensor axes;
     vectorField xyz(p.size());
-    transformation(p, axes, xyz);
+    transformation(t, p, axes, xyz);
 
     return
         xyz.component(2)
@@ -268,14 +270,14 @@ Foam::tmp<Foam::vectorField> Foam::waveSuperposition::ULiquid
 {
     tensor axes;
     vectorField xyz(p.size());
-    transformation(p, axes, xyz);
+    transformation(t, p, axes, xyz);
 
     if (heightAboveWave_)
     {
         xyz.replace(2, height(t, p));
     }
 
-    return velocity(t, xyz) & axes;
+    return UMean_->value(t) + (velocity(t, xyz) & axes);
 }
 
 
@@ -287,7 +289,7 @@ Foam::tmp<Foam::vectorField> Foam::waveSuperposition::UGas
 {
     tensor axes;
     vectorField xyz(p.size());
-    transformation(p, axes, xyz);
+    transformation(t, p, axes, xyz);
 
     axes = tensor(- axes.x(), - axes.y(), axes.z());
 
@@ -298,7 +300,7 @@ Foam::tmp<Foam::vectorField> Foam::waveSuperposition::UGas
 
     xyz.replace(2, - xyz.component(2));
 
-    return velocity(t, xyz) & axes;
+    return UMean_->value(t) + (velocity(t, xyz) & axes);
 }
 
 
@@ -310,7 +312,7 @@ Foam::tmp<Foam::scalarField> Foam::waveSuperposition::pLiquid
 {
     tensor axes;
     vectorField xyz(p.size());
-    transformation(p, axes, xyz);
+    transformation(t, p, axes, xyz);
 
     if (heightAboveWave_)
     {
@@ -329,7 +331,7 @@ Foam::tmp<Foam::scalarField> Foam::waveSuperposition::pGas
 {
     tensor axes;
     vectorField xyz(p.size());
-    transformation(p, axes, xyz);
+    transformation(t, p, axes, xyz);
 
     axes = tensor(- axes.x(), - axes.y(), axes.z());
 
@@ -358,6 +360,7 @@ void Foam::waveSuperposition::write(Ostream& os) const
             << nl << decrIndent << indent << token::END_BLOCK << nl;
     }
     os  << decrIndent << token::END_LIST << token::END_STATEMENT << nl;
+    UMean_->writeData(os);
     if (scale_.valid())
     {
         scale_->writeData(os);

src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.H

@@ -36,6 +36,7 @@ Usage
         origin     | origin of the wave coordinate system    | yes    |
         direction  | direction of the wave coordinate system | yes    |
         waves      | list of wave models to superimpose      | yes    |
+        UMean      | velocity of the mean flow               | yes    |
         scale      | scale factor in the direction           | no     | None
         crossScale | scale factor perpendicular to the direction | no | None
         heightAboveWave | use the height above the wave as the vertical \\
@@ -63,6 +64,7 @@ Usage
             angle       0;
         }
     );
+    UMean       (2 0 0);
     scale       table ((100 1) (200 0));
     crossScale  constant 1;
     heightAboveWave no;
@@ -92,7 +94,9 @@ class waveSuperposition
 :
     public IOdictionary
 {
-    // Private Data
+protected:
+
+    // Protected Data
 
         //- The origin of the wave coordinate system
         const vector origin_;
@@ -106,6 +110,9 @@ class waveSuperposition
         //- The angle relative to the direction at which the waves propagate
         scalarList waveAngles_;
 
+        //- Mean velocity
+        const autoPtr<Function1<vector>> UMean_;
+
         //- Scaling in the local x-direction
         const autoPtr<Function1<scalar>> scale_;
 
@@ -117,11 +124,12 @@ class waveSuperposition
         const Switch heightAboveWave_;
 
 
-    // Private Member Functions
+    // Protected Member Functions
 
         //- Get the transformation to actual coordinates
         void transformation
         (
+            const scalar t,
             const vectorField& p,
             tensor& axes,
             vectorField& xyz
@@ -193,20 +201,6 @@ public:
 
     // Member Functions
 
-        // Access
-
-            //- Return the origin of the wave coordinate system
-            const vector& origin() const
-            {
-                return origin_;
-            }
-
-            //- Return the direction of the wave coordinate system
-            const vector& direction() const
-            {
-                return direction_;
-            }
-
         //- Get the height above the waves at a given time and global positions
         virtual tmp<scalarField> height
         (

tutorials/multiphase/interFoam/RAS/DTCHullWave/0/U.orig

@@ -15,11 +15,11 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-UMean 1.668;
+#include        "$FOAM_CASE/constant/waveProperties";
 
 dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform (#neg $UMean 0 0);
+internalField   uniform ($UxMean 0 0);
 
 boundaryField
 {
@@ -29,14 +29,13 @@ boundaryField
     inlet
     {
         type            waveVelocity;
-        UMean           (#neg $UMean 0 0);
     }
 
     outlet
     {
         type            outletPhaseMeanVelocity;
         alpha           alpha.water;
-        UnMean          $UMean;
+        UnMean          #neg $UxMean;
         value           $internalField;
     }
 

tutorials/multiphase/interFoam/RAS/DTCHullWave/constant/waveProperties

@@ -30,6 +30,10 @@ waves
     }
 );
 
+UxMean          -1.668;
+
+UMean           ($UxMean 0 0);
+
 scale           table ((4 1) (12 0));
 
 

tutorials/multiphase/interFoam/RAS/DTCHullWave/system/setWavesDict

@@ -14,14 +14,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-U.orig
-{
-    #include "$FOAM_CASE/0/U.orig";
-}
-
 alpha   alpha.water;
 
-UMean   $U.orig.boundaryField.inlet.UMean;
-
 
 // ************************************************************************* //

tutorials/multiphase/interFoam/laminar/wave/0/U.orig

@@ -26,7 +26,6 @@ boundaryField
     left
     {
         type            waveVelocity;
-        UMean           (2 0 0);
     }
     right
     {

tutorials/multiphase/interFoam/laminar/wave/constant/waveProperties

@@ -30,6 +30,8 @@ waves
     }
 );
 
+UMean           (2 0 0);
+
 scale           table ((1200 1) (1800 0));
 
 crossScale      constant 1;

tutorials/multiphase/interFoam/laminar/wave/system/blockMeshDict

@@ -65,7 +65,7 @@ boundary
     }
     bottom
     {
-        type patch;
+        type wall;
         faces
         (
             (0 1 5 4)

tutorials/multiphase/interFoam/laminar/wave/system/setWavesDict

@@ -16,7 +16,5 @@ FoamFile
 
 alpha   alpha.water;
 
-UMean   (2 0 0);
-
 
 // ************************************************************************* //