diff --git a/src/waves/Make/files b/src/waves/Make/files index 985d87616..4d4f0361a 100644 --- a/src/waves/Make/files +++ b/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 diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C index 48e39a555..f9f0d7f30 100644 --- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C +++ b/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("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("phi"); + patch().lookupPatchField(phiName_); - const symmTensorField nnp(sqr(patch().nf())); + const scalarField out(pos0(phip)); - // Fix all components if inflow, just normal if outflow - valueFraction() = (1 - pos0(phip))*I + pos0(phip)*nnp; + // Where inflow, fix all velocity components to values specified by the + // 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(os, "phi", "phi", phiName_); waves_.write(os); } diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H index d0ce0e4c5..1e3bec546 100644 --- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H +++ b/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_; diff --git a/src/waves/waveModels/Airy/Airy.C b/src/waves/waveModels/Airy/Airy.C index 7a639951e..bea8cb985 100644 --- a/src/waves/waveModels/Airy/Airy.C +++ b/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::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::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("depth", log(2*GREAT)/k())) {} @@ -76,22 +147,21 @@ Foam::tmp Foam::waveModels::Airy::velocity const vector2DField& xz ) const { - const scalarField x(xz.component(0)); - const scalarField z(xz.component(1)); + const scalar ka = k()*amplitude(t); - const scalarField phi(angle(t, u, x)); - const scalarField kz(- k()*mag(z - elevation(t, u, x))); - const scalar wa = omega(u)*amplitude(t); + return celerity()*ka*vi(1, t, u, xz); +} - 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(); - 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)); + os.writeKeyword("depth") << depth_ << token::END_STATEMENT << nl; } } diff --git a/src/waves/waveModels/Airy/Airy.H b/src/waves/waveModels/Airy/Airy.H index 761614546..1ee70783d 100644 --- a/src/waves/waveModels/Airy/Airy.H +++ b/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 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 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 elevation @@ -102,6 +168,9 @@ public: const scalar u, const vector2DField& xz ) const; + + //- Write + virtual void write(Ostream& os) const; }; diff --git a/src/waves/waveModels/Stokes2/Stokes2.C b/src/waves/waveModels/Stokes2/Stokes2.C index 958734e1c..3145f66b9 100644 --- a/src/waves/waveModels/Stokes2/Stokes2.C +++ b/src/waves/waveModels/Stokes2/Stokes2.C @@ -65,18 +65,20 @@ Foam::tmp Foam::waveModels::Stokes2::elevation const scalarField& x ) const { - const scalar kh = k()*depth(); - const scalarField correction(k()*sqr(amplitude(t))*cos(2*angle(t, u, x))/2); + const scalar kd = k()*depth(), ka = k()*amplitude(t); - 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); - return Airy::elevation(t, u, x) + factor*correction; - } - else - { - return Airy::elevation(t, u, x) + correction; + Info<< "B22 = " << B22 << endl; } + + return + Airy::elevation(t, u, x) + + (1/k())*sqr(ka)*B22*cos(2*angle(t, u, x)); } @@ -87,29 +89,19 @@ Foam::tmp Foam::waveModels::Stokes2::velocity const vector2DField& xz ) const { - if (shallow()) - { - const scalarField x(xz.component(0)); - const scalarField z(xz.component(1)); + const scalar kd = k()*depth(), ka = k()*amplitude(t); - const scalarField phi(angle(t, u, x)); - const scalarField kz(- k()*mag(z - elevation(t, u, x))); - const scalar kh = k()*depth(); - const scalarField khz((kh + kz)*pos0(kh + kz)); - const scalar kwaa = k()*omega(u)*sqr(amplitude(t)); + const scalar A22ByA11 = deep() ? 0 : 0.375/pow3(sinh(kd)); - return - Airy::velocity(t, u, xz) - + 6*kwaa*zip - ( - cosh(2*khz)*cos(2*phi), - sinh(2*khz)*sin(2*phi) - )/pow4(sinh(kh)); - } - else + if (debug) { - return Airy::velocity(t, u, xz); + const scalar A11 = 1/sinh(kd); + Info<< "A22 = " << A22ByA11*A11 << endl; } + + return + Airy::velocity(t, u, xz) + + celerity()*sqr(ka)*A22ByA11*vi(2, t, u, xz); } diff --git a/src/waves/waveModels/Stokes5/Stokes5.C b/src/waves/waveModels/Stokes5/Stokes5.C new file mode 100644 index 000000000..237b69822 --- /dev/null +++ b/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 . + +\*---------------------------------------------------------------------------*/ + +#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::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::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)) + ); +} + + +// ************************************************************************* // diff --git a/src/waves/waveModels/Stokes5/Stokes5.H b/src/waves/waveModels/Stokes5/Stokes5.H new file mode 100644 index 000000000..19302c37f --- /dev/null +++ b/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 . + +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 clone() const + { + return autoPtr(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 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 velocity + ( + const scalar t, + const scalar u, + const vector2DField& xz + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam +} // End namespace waveModels + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/waves/waveModels/solitary/solitary.C b/src/waves/waveModels/solitary/solitary.C new file mode 100644 index 000000000..5143c4034 --- /dev/null +++ b/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 . + +\*---------------------------------------------------------------------------*/ + +#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::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::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::waveModels::solitary::elevation +( + const scalar t, + const scalar u, + const scalarField& x +) const +{ + return amplitude(t)*Pi(t, u, x); +} + + +Foam::tmp 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; +} + + +// ************************************************************************* // diff --git a/src/waves/waveModels/solitary/solitary.H b/src/waves/waveModels/solitary/solitary.H new file mode 100644 index 000000000..5b7880f0b --- /dev/null +++ b/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 . + +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 parameter + ( + const scalar t, + const scalar u, + const scalarField& x + ) const; + + //- The dimensionless elevation [1] + tmp 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 clone() const + { + return autoPtr(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 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 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 + +// ************************************************************************* // diff --git a/src/waves/waveModels/waveModel/waveModel.C b/src/waves/waveModels/waveModel/waveModel.C index a461d13e7..572a32102 100644 --- a/src/waves/waveModels/waveModel/waveModel.C +++ b/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("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::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), - phase_(wave.phase_), - depth_(wave.depth_) + amplitude_(wave.amplitude_, false) {} -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::New("amplitude", dict)), - phase_(readScalar(dict.lookup("phase"))), - depth_(dict.lookupOrDefault("depth", GREAT)) + amplitude_(Function1::New("amplitude", dict)) {} @@ -107,12 +64,15 @@ Foam::waveModel::~waveModel() // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // +Foam::scalar Foam::waveModel::g() const +{ + return mag(db_.lookupObject("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; } diff --git a/src/waves/waveModels/waveModel/waveModel.H b/src/waves/waveModels/waveModel/waveModel.H index e972464c0..cccafe0d6 100644 --- a/src/waves/waveModels/waveModel/waveModel.H +++ b/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> 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 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 - scalar phase() const - { - return phase_; - } - - //- Get the depth - scalar depth() const - { - return depth_; - } + //- Get the (scalar) value of gravity. + scalar g() const; //- Get the wave elevation at a given time, mean velocity and local // coordinates. Local x is aligned with the mean velocity.