waves: Improved handling of reversed flow and additional wave models

The wave library from dev has been copied over into 5.x. It is still
backwards compatible. This change has fixed a failure in the DTCHullWave
tutorial.

src/waves/Make/files

@@ -2,6 +2,8 @@ waveModels/waveModel/waveModel.C
 waveModels/waveModel/newWaveModel.C
 waveModels/Airy/Airy.C
 waveModels/Stokes2/Stokes2.C
+waveModels/Stokes5/Stokes5.C
+waveModels/solitary/solitary.C
 
 waveSuperposition/waveSuperposition.C
 

src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C

@@ -38,6 +38,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
 )
 :
     directionMixedFvPatchVectorField(p, iF),
+    phiName_("phi"),
     waves_(db())
 {
     refValue() = Zero;
@@ -54,6 +55,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
 )
 :
     directionMixedFvPatchVectorField(p, iF),
+    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
     waves_(db(), dict)
 {
     if (dict.found("value"))
@@ -80,6 +82,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
 )
 :
     directionMixedFvPatchVectorField(ptf, p, iF, mapper),
+    phiName_(ptf.phiName_),
     waves_(ptf.waves_)
 {}
 
@@ -90,6 +93,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
 )
 :
     directionMixedFvPatchVectorField(ptf),
+    phiName_(ptf.phiName_),
     waves_(ptf.waves_)
 {}
 
@@ -101,6 +105,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
 )
 :
     directionMixedFvPatchVectorField(ptf, iF),
+    phiName_(ptf.phiName_),
     waves_(ptf.waves_)
 {}
 
@@ -132,15 +137,33 @@ void Foam::waveVelocityFvPatchVectorField::updateCoeffs()
         return;
     }
 
-    refValue() = U();
+    const vectorField UWave(U());
 
     const scalarField& phip =
-        patch().lookupPatchField<surfaceScalarField, scalar>("phi");
+        patch().lookupPatchField<surfaceScalarField, scalar>(phiName_);
 
-    const symmTensorField nnp(sqr(patch().nf()));
+    const scalarField out(pos0(phip));
 
-    // Fix all components if inflow, just normal if outflow
+    // Where inflow, fix all velocity components to values specified by the
-    valueFraction() = (1 - pos0(phip))*I + pos0(phip)*nnp;
+    // wave model.
+    refValue() = (1 - out)*UWave;
+    valueFraction() = (1 - out)*symmTensor::I;
+
+    // Where outflow, set the normal component of the velocity to a value
+    // consistent with phi, but scale it to get the volumentic flow rate
+    // specified by the wave model. Tangential components are extrapolated.
+    const scalar QPhip = gSum(out*phip);
+    const scalar QWave = gSum(out*(UWave & patch().Sf()));
+    const vectorField nBySf(patch().Sf()/sqr(patch().magSf()));
+    if (QPhip > VSMALL)
+    {
+        refValue() += out*(QWave/QPhip)*phip*nBySf;
+    }
+    else
+    {
+        refValue() += out*QWave*nBySf;
+    }
+    valueFraction() += out*sqr(patch().nf());
 
     directionMixedFvPatchVectorField::updateCoeffs();
     directionMixedFvPatchVectorField::evaluate();
@@ -153,6 +176,7 @@ void Foam::waveVelocityFvPatchVectorField::write
 ) const
 {
     directionMixedFvPatchVectorField::write(os);
+    writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     waves_.write(os);
 }
 

src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H

@@ -42,6 +42,7 @@ Usage
         waves       | list of wave models to superimpose     | yes    |
         scale       | scale factor along the mean flow direction | no | None
         crossScale  | scale factor across the mean flow direction | no | None
+        phi         | Name of the flux field                 | no     | phi
     \endtable
 
     Example of the boundary condition specification:
@@ -100,6 +101,9 @@ class waveVelocityFvPatchVectorField
 {
     // Private data
 
+        //- Name of the flux field
+        const word phiName_;
+
         //- Wave superposition
         const waveSuperposition waves_;
 

src/waves/waveModels/Airy/Airy.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Airy.H"
+#include "mathematicalConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -38,15 +39,85 @@ namespace waveModels
 }
 
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::waveModels::Airy::k() const
+{
+    return 2*Foam::constant::mathematical::pi/length_;
+}
+
+
+Foam::scalar Foam::waveModels::Airy::celerity() const
+{
+    return sqrt(g()/k()*tanh(k()*depth()));
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::waveModels::Airy::angle
+(
+    const scalar t,
+    const scalar u,
+    const scalarField& x
+) const
+{
+    return phase_ + k()*(x - (u + celerity())*t);
+}
+
+
+bool Foam::waveModels::Airy::deep() const
+{
+    return k()*depth() > log(GREAT);
+}
+
+
+Foam::tmp<Foam::vector2DField> Foam::waveModels::Airy::vi
+(
+    const label i,
+    const scalar t,
+    const scalar u,
+    const vector2DField& xz
+) const
+{
+    const scalarField x(xz.component(0));
+    const scalarField z(xz.component(1));
+
+    const scalarField phi(angle(t, u, x));
+    const scalarField kz(- k()*z);
+
+    if (deep())
+    {
+        return i*exp(- mag(kz))*zip(cos(i*phi), sin(i*phi));
+    }
+    else
+    {
+        const scalar kd = k()*depth();
+        const scalarField kdz(max(scalar(0), kd - mag(kz)));
+        return i*zip(cosh(i*kdz)*cos(i*phi), sinh(i*kdz)*sin(i*phi))/sinh(kd);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+Foam::waveModels::Airy::Airy(const Airy& wave)
+:
+    waveModel(wave),
+    length_(wave.length_),
+    phase_(wave.phase_),
+    depth_(wave.depth_)
+{}
+
+
 Foam::waveModels::Airy::Airy
 (
     const objectRegistry& db,
     const dictionary& dict
 )
 :
-    waveModel(db, dict)
+    waveModel(db, dict),
+    length_(readScalar(dict.lookup("length"))),
+    phase_(readScalar(dict.lookup("phase"))),
+    depth_(dict.lookupOrDefault<scalar>("depth", log(2*GREAT)/k()))
 {}
 
 
@@ -76,22 +147,21 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Airy::velocity
     const vector2DField& xz
 ) const
 {
-    const scalarField x(xz.component(0));
+    const scalar ka = k()*amplitude(t);
-    const scalarField z(xz.component(1));
 
-    const scalarField phi(angle(t, u, x));
+    return celerity()*ka*vi(1, t, u, xz);
-    const scalarField kz(- k()*mag(z - elevation(t, u, x)));
+}
-    const scalar wa = omega(u)*amplitude(t);
 
-    if (shallow())
+
+void Foam::waveModels::Airy::write(Ostream& os) const
+{
+    waveModel::write(os);
+
+    os.writeKeyword("length") << length_ << token::END_STATEMENT << nl;
+    os.writeKeyword("phase") << phase_ << token::END_STATEMENT << nl;
+    if (!deep())
     {
-        const scalar kh = k()*depth();
+        os.writeKeyword("depth") << depth_ << token::END_STATEMENT << nl;
-        const scalarField khz((kh + kz)*pos0(kh + kz));
-        return wa*zip(cosh(khz)*cos(phi), sinh(khz)*sin(phi))/sinh(kh);
-    }
-    else
-    {
-        return wa*exp(kz)*zip(cos(phi), sin(phi));
     }
 }
 

src/waves/waveModels/Airy/Airy.H

@@ -60,6 +60,49 @@ class Airy
 :
     public waveModel
 {
+    //- Private data
+
+        //- Peak-to-peak length [m]
+        const scalar length_;
+
+        //- Phase offset [rad]
+        const scalar phase_;
+
+        //- Depth [m]
+        const scalar depth_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- The angular wavenumber [rad/m]
+        scalar k() const;
+
+        //- The wave celerity [m/s]
+        scalar celerity() const;
+
+        //- Angle of the oscillation [rad]
+        tmp<scalarField> angle
+        (
+            const scalar t,
+            const scalar u,
+            const scalarField& x
+        ) const;
+
+        //- Return whether shallow and intermdiate effects are to be omitted
+        bool deep() const;
+
+        //- Return the non-dimensionalised i-th harmonic of the velocity
+        tmp<vector2DField> vi
+        (
+            const label i,
+            const scalar t,
+            const scalar u,
+            const vector2DField& xz
+        ) const;
+
+
 public:
 
     //- Runtime type information
@@ -68,6 +111,9 @@ public:
 
     // Constructors
 
+        //- Construct a copy
+        Airy(const Airy& wave);
+
         //- Construct from a database and a dictionary
         Airy(const objectRegistry& db, const dictionary& dict);
 
@@ -84,6 +130,26 @@ public:
 
     // Member Functions
 
+        // Access
+
+            //- Get the length
+            scalar length() const
+            {
+                return length_;
+            }
+
+            //- Get the phase
+            scalar phase() const
+            {
+                return phase_;
+            }
+
+            //- Get the depth
+            scalar depth() const
+            {
+                return depth_;
+            }
+
         //- Get the wave elevation at a given time, mean velocity and local
         //  coordinates. Local x is aligned with the mean velocity.
         virtual tmp<scalarField> elevation
@@ -102,6 +168,9 @@ public:
             const scalar u,
             const vector2DField& xz
         ) const;
+
+        //- Write
+        virtual void write(Ostream& os) const;
 };
 
 

src/waves/waveModels/Stokes2/Stokes2.C

@@ -65,18 +65,20 @@ Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes2::elevation
     const scalarField& x
 ) const
 {
-    const scalar kh = k()*depth();
+    const scalar kd = k()*depth(), ka = k()*amplitude(t);
-    const scalarField correction(k()*sqr(amplitude(t))*cos(2*angle(t, u, x))/2);
 
-    if (shallow())
+    const scalar T = deep() ? 1 : tanh(kd);
+
+    const scalar B22 = (3/sqr(T) - 1)/T/4;
+
+    if (debug)
     {
-        const scalar factor((3/sqr(tanh(kh)) - 1)/tanh(kh)/2);
+        Info<< "B22 = " << B22 << endl;
-        return Airy::elevation(t, u, x) + factor*correction;
-    }
-    else
-    {
-        return Airy::elevation(t, u, x) + correction;
     }
+
+    return
+        Airy::elevation(t, u, x)
+      + (1/k())*sqr(ka)*B22*cos(2*angle(t, u, x));
 }
 
 
@@ -87,29 +89,19 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Stokes2::velocity
     const vector2DField& xz
 ) const
 {
-    if (shallow())
+    const scalar kd = k()*depth(), ka = k()*amplitude(t);
-    {
-        const scalarField x(xz.component(0));
-        const scalarField z(xz.component(1));
 
-        const scalarField phi(angle(t, u, x));
+    const scalar A22ByA11 = deep() ? 0 : 0.375/pow3(sinh(kd));
-        const scalarField kz(- k()*mag(z - elevation(t, u, x)));
+
-        const scalar kh = k()*depth();
+    if (debug)
-        const scalarField khz((kh + kz)*pos0(kh + kz));
+    {
-        const scalar kwaa = k()*omega(u)*sqr(amplitude(t));
+        const scalar A11 = 1/sinh(kd);
+        Info<< "A22 = " << A22ByA11*A11 << endl;
+    }
 
     return
         Airy::velocity(t, u, xz)
-          + 6*kwaa*zip
+      + celerity()*sqr(ka)*A22ByA11*vi(2, t, u, xz);
-            (
-                cosh(2*khz)*cos(2*phi),
-                sinh(2*khz)*sin(2*phi)
-            )/pow4(sinh(kh));
-    }
-    else
-    {
-        return Airy::velocity(t, u, xz);
-    }
 }
 
 

src/waves/waveModels/Stokes5/Stokes5.C

@@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "Stokes5.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveModels
+{
+    defineTypeNameAndDebug(Stokes5, 0);
+    addToRunTimeSelectionTable(waveModel, Stokes5, objectRegistry);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::waveModels::Stokes5::Stokes5
+(
+    const objectRegistry& db,
+    const dictionary& dict
+)
+:
+    Stokes2(db, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::waveModels::Stokes5::~Stokes5()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes5::elevation
+(
+    const scalar t,
+    const scalar u,
+    const scalarField& x
+) const
+{
+    const scalar kd = k()*depth(), ka = k()*amplitude(t);
+
+    const scalar S = deep() ? 0 : 1/cosh(2*kd), T = deep() ? 1 : tanh(kd);
+
+    const scalar B31 =
+      - 3.0/8/pow3(1 - S)
+       *(
+            1 + 3*S + 3*sqr(S) + 2*pow3(S)
+        );
+
+    const scalar B42 =
+        1.0/6/T/(3 + 2*S)/pow4(1 - S)
+       *(
+            6 - 26*S - 182*sqr(S) - 204*pow3(S) - 25*pow4(S) + 26*pow5(S)
+        );
+
+    const scalar B44 =
+        1.0/24/T/(3 + 2*S)/pow4(1 - S)
+       *(
+            24 + 92*S + 122*sqr(S) + 66*pow3(S) + 67*pow4(S) + 34*pow5(S)
+        );
+
+    const scalar B53 =
+        9.0/128/(3 + 2*S)/(4 + S)/pow6(1 - S)
+       *(
+            132 + 17*S - 2216*sqr(S) - 5897*pow3(S) - 6292*pow4(S)
+          - 2687*pow5(S) + 194*pow6(S) + 467*S*pow6(S) + 82*sqr(pow4(S))
+        );
+
+    const scalar B55 =
+        5.0/384/(3 + 2*S)/(4 + S)/pow6(1 - S)
+       *(
+            300 + 1579*S + 3176*sqr(S) + 2949*pow3(S) + 1188*pow4(S)
+          + 675*pow5(S) + 1326*pow6(S) + 827*S*pow6(S) + 130*sqr(pow4(S))
+        );
+
+    if (debug)
+    {
+        Info<< "B31 = " << B31 << endl
+            << "B42 = " << B42 << endl
+            << "B44 = " << B44 << endl
+            << "B53 = " << B53 << endl
+            << "B55 = " << B55 << endl;
+    }
+
+    const scalarField phi(angle(t, u, x));
+
+    return
+        Stokes2::elevation(t, u, x)
+      + (1/k())
+       *(
+           pow3(ka)*B31*(cos(phi) - cos(3*phi))
+         + pow4(ka)*(B42*cos(2*phi) + B44*cos(4*phi))
+         + pow5(ka)*(- (B53 + B55)*cos(phi) + B53*cos(3*phi) + B55*cos(5*phi))
+        );
+}
+
+
+Foam::tmp<Foam::vector2DField> Foam::waveModels::Stokes5::velocity
+(
+    const scalar t,
+    const scalar u,
+    const vector2DField& xz
+) const
+{
+    const scalar kd = k()*depth(), ka = k()*amplitude(t);
+
+    const scalar S = deep() ? 0 : 1/cosh(2*kd);
+    const scalar SByA11 = deep() ? 0 : S*sinh(kd);
+
+    const scalar A31ByA11 =
+        1.0/8/pow3(1 - S)
+       *(
+          - 4 - 20*S + 10*sqr(S) - 13*pow3(S)
+        );
+
+    const scalar A33ByA11 =
+        1.0/8/pow3(1 - S)
+       *(
+          - 2*sqr(S) + 11*pow3(S)
+        );
+
+    const scalar A42ByA11 =
+        SByA11/24/pow5(1 - S)
+       *(
+            12 - 14*S - 264*sqr(S) - 45*pow3(S) - 13*pow4(S)
+        );
+
+    const scalar A44ByA11 =
+        SByA11/48/(3 + 2*S)/pow5(1 - S)
+       *(
+            10*sqr(S) - 174*pow3(S) + 291*pow4(S) + 278*pow5(S)
+        );
+
+    const scalar A51ByA11 =
+        1.0/64/(3 + 2*S)/(4 + S)/pow6(1 - S)
+       *(
+         - 1184 + 32*S + 13232*sqr(S) + 21712*pow3(S) + 20940*pow4(S)
+         + 12554*pow5(S) - 500*pow6(S) - 3341*S*pow6(S) - 670*sqr(pow4(S))
+        );
+
+    const scalar A53ByA11 =
+        1.0/32/(3 + 2*S)/pow6(1 - S)
+       *(
+            4*S + 105*sqr(S) + 198*pow3(S) - 1376*pow4(S) - 1302*pow5(S)
+          - 117*pow6(S) + 58*S*pow6(S)
+        );
+
+    const scalar A55ByA11 =
+        1.0/64/(3 + 2*S)/(4 + S)/pow6(1 - S)
+       *(
+          - 6*pow3(S) + 272*pow4(S) - 1552*pow5(S) + 852*pow6(S)
+          + 2029*S*pow6(S) + 430*sqr(pow4(S))
+        );
+
+    if (debug)
+    {
+        const scalar A11 = 1/sinh(kd);
+        Info<< "A31 = " << A31ByA11*A11 << endl
+            << "A33 = " << A33ByA11*A11 << endl
+            << "A42 = " << A42ByA11*A11 << endl
+            << "A44 = " << A44ByA11*A11 << endl
+            << "A51 = " << A51ByA11*A11 << endl
+            << "A53 = " << A53ByA11*A11 << endl
+            << "A55 = " << A55ByA11*A11 << endl;
+    }
+
+    const vector2DField v1(vi(1, t, u, xz)), v3(vi(3, t, u, xz));
+
+    return
+        Stokes2::velocity(t, u, xz)
+      + celerity()
+       *(
+            pow3(ka)*(A31ByA11*v1 + A33ByA11*v3)
+          + pow4(ka)*(A42ByA11*vi(2, t, u, xz) + A44ByA11*vi(4, t, u, xz))
+          + pow5(ka)*(A51ByA11*v1 + A53ByA11*v3 + A55ByA11*vi(5, t, u, xz))
+        );
+}
+
+
+// ************************************************************************* //

src/waves/waveModels/Stokes5/Stokes5.H

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::waveModels::Stokes5
+
+Description
+    Fifth-order wave model.
+
+    Reference:
+    \verbatim
+        "A Fifth Order Stokes Theory for Steady Waves"
+        J D Fenton
+        Journal of Waterway, Port, Coastal, and Ocean Engineering (1985),
+        Volume 111, Issue 2, Pages 216-234
+    \endverbatim
+
+SourceFiles
+    Stokes5.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Stokes5_H
+#define Stokes5_H
+
+#include "Stokes2.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Stokes5 Declaration
+\*---------------------------------------------------------------------------*/
+
+class Stokes5
+:
+    public Stokes2
+{
+public:
+
+   //- Runtime type information
+    TypeName("Stokes5");
+
+
+    // Constructors
+
+        //- Construct from a database and a dictionary
+        Stokes5(const objectRegistry& db, const dictionary& dict);
+
+        //- Construct a clone
+        virtual autoPtr<waveModel> clone() const
+        {
+            return autoPtr<waveModel>(new Stokes5(*this));
+        }
+
+
+    //- Destructor
+    virtual ~Stokes5();
+
+
+    // Member Functions
+
+        //- Get the wave elevation at a given time, mean velocity and local
+        //  coordinates. Local x is aligned with the mean velocity.
+        virtual tmp<scalarField> elevation
+        (
+            const scalar t,
+            const scalar u,
+            const scalarField& x
+        ) const;
+
+        //- Get the wave velocity at a given time, mean velocity and local
+        //  coordinates. Local x is aligned with the mean velocity, and z with
+        //  negative gravity.
+        virtual tmp<vector2DField> velocity
+        (
+            const scalar t,
+            const scalar u,
+            const vector2DField& xz
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace waveModels
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/waves/waveModels/solitary/solitary.C

@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "solitary.H"
+#include "mathematicalConstants.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveModels
+{
+    defineTypeNameAndDebug(solitary, 0);
+    addToRunTimeSelectionTable(waveModel, solitary, objectRegistry);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::waveModels::solitary::k(const scalar t) const
+{
+    return sqrt(0.75*amplitude(t)/pow3(depth()));
+}
+
+
+Foam::scalar Foam::waveModels::solitary::alpha(const scalar t) const
+{
+    return amplitude(t)/depth();
+}
+
+
+Foam::scalar Foam::waveModels::solitary::celerity(const scalar t) const
+{
+    return sqrt(g()*depth()/(1 - alpha(t)));
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::waveModels::solitary::parameter
+(
+    const scalar t,
+    const scalar u,
+    const scalarField& x
+) const
+{
+    return k(t)*(x - offset_ - (u + celerity(t))*t);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::waveModels::solitary::Pi
+(
+    const scalar t,
+    const scalar u,
+    const scalarField& x
+) const
+{
+    const scalar clip = log(GREAT);
+
+    return 1/sqr(cosh(max(- clip, min(clip, parameter(t, u, x)))));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::waveModels::solitary::solitary(const solitary& wave)
+:
+    waveModel(wave),
+    offset_(wave.offset_),
+    depth_(wave.depth_)
+{}
+
+
+Foam::waveModels::solitary::solitary
+(
+    const objectRegistry& db,
+    const dictionary& dict
+)
+:
+    waveModel(db, dict),
+    offset_(readScalar(dict.lookup("offset"))),
+    depth_(readScalar(dict.lookup("depth")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::waveModels::solitary::~solitary()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField> Foam::waveModels::solitary::elevation
+(
+    const scalar t,
+    const scalar u,
+    const scalarField& x
+) const
+{
+    return amplitude(t)*Pi(t, u, x);
+}
+
+
+Foam::tmp<Foam::vector2DField> Foam::waveModels::solitary::velocity
+(
+    const scalar t,
+    const scalar u,
+    const vector2DField& xz
+) const
+{
+    const scalar A = alpha(t);
+    const scalarField Z(max(scalar(0), 1 - mag(xz.component(1)/depth())));
+    const scalarField P(Pi(t, u, xz.component(0)));
+
+    return
+        celerity(t)
+       *zip
+        (
+            A/4
+           *(
+                (4 + 2*A - 6*A*sqr(Z))*P
+              + (- 7*A + 9*A*sqr(Z))*sqr(P)
+            ),
+          - A*Z*depth()*k(t)*tanh(parameter(t, u, xz.component(0)))
+           *(
+                (2 + A - A*sqr(Z))*P
+              + (- 7*A + 3*A*sqr(Z))*sqr(P)
+            )
+        );
+}
+
+
+void Foam::waveModels::solitary::write(Ostream& os) const
+{
+    waveModel::write(os);
+
+    os.writeKeyword("offset") << offset_ << token::END_STATEMENT << nl;
+    os.writeKeyword("depth") << depth_ << token::END_STATEMENT << nl;
+}
+
+
+// ************************************************************************* //

src/waves/waveModels/solitary/solitary.H

@@ -0,0 +1,174 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::waveModels::solitary
+
+Description
+    Solitary wave model.
+
+    Reference:
+    \verbatim
+        "Water Wave Mechanics For Engineers And Scientists"
+        R G Dean and R A Dalrymple
+        Pages 314-317
+    \endverbatim
+
+SourceFiles
+    solitary.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solitary_H
+#define solitary_H
+
+#include "waveModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class solitary Declaration
+\*---------------------------------------------------------------------------*/
+
+class solitary
+:
+    public waveModel
+{
+    //- Private data
+
+        //- Offset [m]
+        const scalar offset_;
+
+        //- Depth [m]
+        const scalar depth_;
+
+
+    // Private Member Functions
+
+        //- The wavenumber [1/m]
+        scalar k(const scalar t) const;
+
+        //- The dimensionless amplitude [1]
+        scalar alpha(const scalar t) const;
+
+        //- The wave celerity [m/s]
+        scalar celerity(const scalar t) const;
+
+        //- The evolution parameter [1]
+        //  This is analogous to the oscillation angle of a periodic wave
+        tmp<scalarField> parameter
+        (
+            const scalar t,
+            const scalar u,
+            const scalarField& x
+        ) const;
+
+        //- The dimensionless elevation [1]
+        tmp<scalarField> Pi
+        (
+            const scalar t,
+            const scalar u,
+            const scalarField& x
+        ) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("solitary");
+
+
+    // Constructors
+
+        //- Construct a copy
+        solitary(const solitary& wave);
+
+        //- Construct from a database and a dictionary
+        solitary(const objectRegistry& db, const dictionary& dict);
+
+        //- Construct a clone
+        virtual autoPtr<waveModel> clone() const
+        {
+            return autoPtr<waveModel>(new solitary(*this));
+        }
+
+
+    //- Destructor
+    virtual ~solitary();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Get the offset
+            scalar offset() const
+            {
+                return offset_;
+            }
+
+            //- Get the depth
+            scalar depth() const
+            {
+                return depth_;
+            }
+
+        //- Get the wave elevation at a given time, mean velocity and local
+        //  coordinates. Local x is aligned with the mean velocity.
+        virtual tmp<scalarField> elevation
+        (
+            const scalar t,
+            const scalar u,
+            const scalarField& x
+        ) const;
+
+        //- Get the wave velocity at a given time, mean velocity and local
+        //  coordinates. Local x is aligned with the mean velocity, and z with
+        //  negative gravity.
+        virtual tmp<vector2DField> velocity
+        (
+            const scalar t,
+            const scalar u,
+            const vector2DField& xz
+        ) const;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace waveModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/waves/waveModels/waveModel/waveModel.C

@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "waveModel.H"
-#include "mathematicalConstants.H"
 #include "Time.H"
 #include "uniformDimensionedFields.H"
 
@@ -37,65 +36,23 @@ namespace Foam
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-Foam::scalar Foam::waveModel::k() const
-{
-    return 2*Foam::constant::mathematical::pi/length_;
-}
-
-
-Foam::scalar Foam::waveModel::sigma() const
-{
-    const uniformDimensionedVectorField& g =
-        db_.lookupObject<uniformDimensionedVectorField>("g");
-
-    return sqrt(mag(g.value())*k()*tanh(k()*depth()));
-}
-
-
-Foam::scalar Foam::waveModel::omega(const scalar u) const
-{
-    return sigma() + k()*u;
-}
-
-
-Foam::tmp<Foam::scalarField> Foam::waveModel::angle
-(
-    const scalar t,
-    const scalar u,
-    const scalarField& x
-) const
-{
-    return k()*x - omega(u)*t;
-}
-
-
-bool Foam::waveModel::shallow() const
-{
-    return k()*depth() < log(GREAT);
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::waveModel::waveModel(const waveModel& wave)
 :
     db_(wave.db_),
-    length_(wave.length_),
+    amplitude_(wave.amplitude_, false)
-    amplitude_(wave.amplitude_, false),
-    phase_(wave.phase_),
-    depth_(wave.depth_)
 {}
 
 
-Foam::waveModel::waveModel(const objectRegistry& db, const dictionary& dict)
+Foam::waveModel::waveModel
+(
+    const objectRegistry& db,
+    const dictionary& dict
+)
 :
     db_(db),
-    length_(readScalar(dict.lookup("length"))),
+    amplitude_(Function1<scalar>::New("amplitude", dict))
-    amplitude_(Function1<scalar>::New("amplitude", dict)),
-    phase_(readScalar(dict.lookup("phase"))),
-    depth_(dict.lookupOrDefault<scalar>("depth", GREAT))
 {}
 
 
@@ -107,12 +64,15 @@ Foam::waveModel::~waveModel()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+Foam::scalar Foam::waveModel::g() const
+{
+    return mag(db_.lookupObject<uniformDimensionedVectorField>("g").value());
+}
+
+
 void Foam::waveModel::write(Ostream& os) const
 {
-    os.writeKeyword("length") << length_ << token::END_STATEMENT << nl;
     amplitude_->writeData(os);
-    os.writeKeyword("phase") << phase_ << token::END_STATEMENT << nl;
-    os.writeKeyword("depth") << depth_ << token::END_STATEMENT << nl;
 }
 
 

src/waves/waveModels/waveModel/waveModel.H

@@ -60,43 +60,9 @@ class waveModel
         //- Reference to the database
         const objectRegistry& db_;
 
-        //- Peak-to-peak length [m]
-        const scalar length_;
-
         //- Peak-to-mean amplitude [m]
         autoPtr<Function1<scalar>> amplitude_;
 
-        //- Phase offset [rad]
-        const scalar phase_;
-
-        //- Depth [m]
-        const scalar depth_;
-
-
-protected:
-
-    // Protected Member Functions
-
-        //- The angular wavenumber [rad/m]
-        scalar k() const;
-
-        //- The intrinsic angular frequency [rad/s]
-        scalar sigma() const;
-
-        //- The observed angular frequency [rad/s]
-        scalar omega(const scalar u) const;
-
-        //- Angle of the oscillation [rad]
-        tmp<scalarField> angle
-        (
-            const scalar t,
-            const scalar u,
-            const scalarField& x
-        ) const;
-
-        //- Return whether shallow effects are to be included
-        bool shallow() const;
-
 
 public:
 
@@ -153,29 +119,14 @@ public:
 
         // Access
 
-            //- Get the length
-            scalar length() const
-            {
-                return length_;
-            }
-
             //- Get the amplitude
             scalar amplitude(const scalar t) const
             {
                 return amplitude_->value(t);
             }
 
-            //- Get the phase
+        //- Get the (scalar) value of gravity.
-            scalar phase() const
+        scalar g() const;
-            {
-                return phase_;
-            }
-
-            //- Get the depth
-            scalar depth() const
-            {
-                return depth_;
-            }
 
         //- Get the wave elevation at a given time, mean velocity and local
         //  coordinates. Local x is aligned with the mean velocity.