Commit 9c740ffc authored by mattijs's avatar mattijs
Browse files

Merge branch 'develop' of develop.openfoam.com:Development/OpenFOAM-plus into develop

parents 510b579d 6c3f04d0
......@@ -6,6 +6,8 @@ waveGenerationModels/base/regularWaveModel/regularWaveModel.C
waveGenerationModels/base/solitaryWaveModel/solitaryWaveModel.C
waveGenerationModels/derived/Boussinesq/BoussinesqWaveModel.C
waveGenerationModels/derived/cnoidal/cnoidalWaveModel.C
waveGenerationModels/derived/Grimshaw/GrimshawWaveModel.C
waveGenerationModels/derived/McCowan/McCowanWaveModel.C
waveGenerationModels/derived/StokesII/StokesIIWaveModel.C
waveGenerationModels/derived/StokesI/StokesIWaveModel.C
waveGenerationModels/derived/StokesV/StokesVWaveModel.C
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2017 IH-Cantabria
-------------------------------------------------------------------------------
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 "GrimshawWaveModel.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace waveModels
{
defineTypeNameAndDebug(Grimshaw, 0);
addToRunTimeSelectionTable
(
waveModel,
Grimshaw,
patch
);
}
}
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
Foam::scalar Foam::waveModels::Grimshaw::alfa
(
const scalar H,
const scalar h
) const
{
scalar eps = H/h;
return sqrt(0.75*eps)*(1.0 - 0.625*eps + (71.0/128.0)*eps*eps);
}
Foam::scalar Foam::waveModels::Grimshaw::eta
(
const scalar H,
const scalar h,
const scalar x,
const scalar y,
const scalar theta,
const scalar t,
const scalar X0
) const
{
const scalar eps = H/h;
const scalar eps2 = eps*eps;
const scalar eps3 = eps*eps2;
const scalar C = sqrt(mag(g_)*h)*sqrt(1.0 + eps - 0.05*eps2 - (3.0/70.0)*eps3);
const scalar ts = 3.5*h/sqrt(H/h);
const scalar xa = -C*t + ts - X0 + x*cos(theta) + y*sin(theta);
const scalar alfa = this->alfa(H, h);
const scalar s = (1.0)/(cosh(alfa*(xa/h)));
const scalar s2 = s*s;
const scalar q = tanh(alfa*(xa/h));
const scalar q2 = q*q;
return
h
*(
eps*s2
- 0.75*eps2*s2*q2
+ eps3*(0.625*s2*q2 - 1.2625*s2*s2*q2)
);
}
Foam::vector Foam::waveModels::Grimshaw::Uf
(
const scalar H,
const scalar h,
const scalar x,
const scalar y,
const scalar theta,
const scalar t,
const scalar X0,
const scalar z
) const
{
const scalar eps = H/h;
const scalar eps2 = eps*eps;
const scalar eps3 = eps*eps2;
const scalar C = sqrt(mag(g_)*h)*sqrt(1.0 + eps - 0.05*eps2 - (3.0/70.0)*eps3);
const scalar ts = 3.5*h/sqrt(eps);
const scalar xa = -C*t + ts - X0 + x*cos(theta) + y*sin(theta);
const scalar alfa = this->alfa(H, h);
const scalar s = (1.0)/(cosh(alfa*(xa/h)));
const scalar s2 = s*s;
const scalar s4 = s2*s2;
const scalar s6 = s2*s4;
const scalar zbyh = z/h;
const scalar zbyh2 = zbyh*zbyh;
const scalar zbyh4 = zbyh2*zbyh2;
scalar outa = eps*s2 - eps2*(-0.25*s2 + s4 + zbyh2*(1.5*s2 - 2.25*s4));
scalar outb = 0.475*s2 + 0.2*s4 - 1.2*s6;
scalar outc = zbyh2*(-1.5*s2 - 3.75*s4 + 7.5*s6);
scalar outd = zbyh4*(-0.375*s2 + (45.0/16.0)*s4 - (45.0/16.0)*s6);
scalar u = sqrt(mag(g_)*h)*(outa - eps3*(outb + outc + outd));
outa = eps*s2 - eps2*(0.375*s2 + 2*s4 + zbyh2*(0.5*s2 - 1.5*s4));
outb = (49.0/640.0)*s2 - 0.85*s4 - 3.6*s6;
outc = zbyh2*((-13.0/16.0)*s2 -(25.0/16.0)*s4 + 7.5*s6);
outd = zbyh4*(-0.075*s2 -1.125*s4 - (27.0/16.0)*s6);
const scalar w = sqrt(mag(g_)*h)*(outa - eps3*(outb + outc + outd));
const scalar v = u*sin(waveAngle_);
u *= cos(waveAngle_);
return vector(u, v, w);
}
void Foam::waveModels::Grimshaw::setLevel
(
const scalar t,
const scalar tCoeff,
scalarField& level
) const
{
forAll(level, paddlei)
{
const scalar eta =
this->eta
(
waveHeight_,
waterDepthRef_,
xPaddle_[paddlei],
yPaddle_[paddlei],
waveAngle_,
t,
x0_
);
level[paddlei] = waterDepthRef_ + tCoeff*eta;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::waveModels::Grimshaw::Grimshaw
(
const dictionary& dict,
const fvMesh& mesh,
const polyPatch& patch,
const bool readFields
)
:
solitaryWaveModel(dict, mesh, patch, false)
{
if (readFields)
{
readDict(dict);
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::waveModels::Grimshaw::~Grimshaw()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::waveModels::Grimshaw::readDict(const dictionary& overrideDict)
{
if (solitaryWaveModel::readDict(overrideDict))
{
return true;
}
return false;
}
void Foam::waveModels::Grimshaw::setVelocity
(
const scalar t,
const scalar tCoeff,
const scalarField& level
)
{
forAll(U_, facei)
{
// Fraction of geometry represented by paddle - to be set
scalar fraction = 1;
// Height - to be set
scalar z = 0;
setPaddlePropeties(level, facei, fraction, z);
if (fraction > 0)
{
const label paddlei = faceToPaddle_[facei];
const vector Uf = this->Uf
(
waveHeight_,
waterDepthRef_,
xPaddle_[paddlei],
yPaddle_[paddlei],
waveAngle_,
t,
x0_,
z
);
U_[facei] = fraction*Uf*tCoeff;
}
}
}
void Foam::waveModels::Grimshaw::info(Ostream& os) const
{
solitaryWaveModel::info(os);
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2017 IH-Cantabria
-------------------------------------------------------------------------------
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::Grimshaw
Description
Grimshaw wave model
\*---------------------------------------------------------------------------*/
#ifndef waveModels_Grimshaw_H
#define waveModels_Grimshaw_H
#include "solitaryWaveModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace waveModels
{
/*---------------------------------------------------------------------------*\
Class Grimshaw Declaration
\*---------------------------------------------------------------------------*/
class Grimshaw
:
public solitaryWaveModel
{
protected:
// Protected Member Functions
//- Wave height
virtual scalar eta
(
const scalar H,
const scalar h,
const scalar x,
const scalar y,
const scalar theta,
const scalar t,
const scalar X0
) const;
//-
virtual scalar alfa
(
const scalar H,
const scalar h
) const;
//- Wave velocity
virtual vector Uf
(
const scalar H,
const scalar h,
const scalar x,
const scalar y,
const scalar theta,
const scalar t,
const scalar X0,
const scalar z
) const;
//- Set the water level
virtual void setLevel
(
const scalar t,
const scalar tCoeff,
scalarField& level
) const;
//- Calculate the wave model velocity
virtual void setVelocity
(
const scalar t,
const scalar tCoeff,
const scalarField& level
);
public:
//- Runtime type information
TypeName("Grimshaw");
//- Constructor
Grimshaw
(
const dictionary& dict,
const fvMesh& mesh,
const polyPatch& patch,
const bool readFields = true
);
//- Destructor
virtual ~Grimshaw();
// Public Member Functions
//- Read from dictionary
virtual bool readDict(const dictionary& overrideDict);
//- Info
virtual void info(Ostream& os) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace waveModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\/ M anipulation | Copyright (C) 2017 IH-Cantabria
-------------------------------------------------------------------------------
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 "McCowanWaveModel.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace waveModels
{
defineTypeNameAndDebug(McCowan, 0);
addToRunTimeSelectionTable
(
waveModel,
McCowan,
patch
);
}
}
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
Foam::scalar Foam::waveModels::McCowan::eta
(
const scalar H,
const scalar h,
const scalar x,
const scalar y,
const scalar theta,
const scalar t,
const scalar X0
) const
{
vector vec = this->mn(H, h);
scalar mm = vec[0];
scalar nn = vec[1];
scalar C = sqrt(((mag(g_)*h)/mm)*tan(mm));
scalar ts = 3.5*h/sqrt(H/h);
scalar Xa = -C*t + ts - X0 + x*cos(theta) + y*sin(theta);
scalar xin = 0.5*H;
scalar etas = newtonRapsonF2(xin, H, h, Xa, mm, nn);
return etas;
}
Foam::vector Foam::waveModels::McCowan::mn
(
const scalar H,
const scalar h
) const
{
// m
scalar xin = 1;
scalar m = newtonRapsonF1(xin, H, h);
// n
scalar c1 = sin(m + (1.0 + (2.0*H/(3.0*h))));
scalar n = (2.0/3.0)*sqr(c1);
return vector(m, n, n);
}
Foam::scalar Foam::waveModels::McCowan::newtonRapsonF1
(
const scalar x0,
const scalar H,
const scalar h
) const
{
label N = 10000;
scalar eps = 1.e-5;
scalar maxval = 10000.0;
label iter = 1;
scalar x = x0;
scalar residual = 0;
while (iter <= N)
{
// f
scalar a = x + 1.0 + 2.0*H/(3.0*h);
scalar b = 0.5*x*(1.0 + H/h);
scalar c = 0.5*x*(1.0 + h/H);
scalar c1 = sin(a);
scalar fx = (2.0/3.0)*sqr(c1) - x*H/(h*tan(b));
residual = mag(fx);
if (residual < eps)
{
return x;
}
else if ((iter > 1) && (residual > maxval))
{
FatalErrorInFunction
<< "Newton-Raphson iterations diverging: "
<< "iterations = " << iter
<< ", residual = " << residual
<< exit(FatalError);
}
// f-prime
scalar c2 = 1.0/tan(c);
scalar c3 = 1.0/sin(b);
scalar fprime = (4.0/3.0)*c1*cos(a) - c2*h/H - b*sqr(c3);
x -= fx/fprime;
iter++;
}
WarningInFunction
<< "Failed to converge in " << iter << " iterations. Residual = "
<< residual << nl << endl;
return x;
}
Foam::scalar Foam::waveModels::McCowan::newtonRapsonF2
(
const scalar x0,
const scalar H,
const scalar h,
const scalar xa,
const scalar m,
const scalar n
) const
{
label N = 10000;
scalar eps = 1.e-5;
scalar maxval = 10000;
label iter = 1;
scalar x = x0;
scalar residual = 0;
while (iter <= N)
{
// f
scalar a = m*(1.0 + x/h);
scalar c1 = cos(a);
scalar c2 = sin(a);
scalar fx = x - (h*n/m*(c2/(c1 + cosh(m*xa/h))));
residual = mag(fx);
if (residual < eps)
{
return x;
}
else if ((iter > 1) && (residual > maxval))