Skip to content

GitLab

Commit d0c56621 authored by Gabriel Barajas's avatar Gabriel Barajas Committed by Andrew Heather
Browse files

INT: Initial IHCantabria update for waveMaker BC 

Adds support for paddles to generate 3-D waves
parent 4818af3b
Hide whitespace changes
Inline Side-by-side
src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.C
View file @ d0c56621
......@@ -33,6 +33,10 @@ License
#include "Time.H"
#include "gravityMeshObject.H"
#include "polyMesh.H"
#include "surfaceFields.H"
#include "volFields.H"
// * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * * //
const Foam::Enum<Foam::waveMakerPointPatchVectorField::motionTypes>
......@@ -89,6 +93,48 @@ Foam::scalar Foam::waveMakerPointPatchVectorField::timeCoeff
}
void Foam::waveMakerPointPatchVectorField::initialiseGeometry()
{
// Global patch extents
const vectorField& Cp = this->patch().localPoints();
const vectorField CpLocal(Cp);
boundBox bb(CpLocal, true);
const scalar xMin = bb.min().x();
const scalar xMax = bb.max().x();
const scalar yMin = bb.min().y();
const scalar yMax = bb.max().y();
zSpan_ = bb.max().z() - bb.min().z();
zMinGb_ = bb.min().z();
reduce(zMinGb_, minOp<scalar>());
// Global x, y positions of the paddle centres
xPaddle_.setSize(nPaddle_, 0);
yPaddle_.setSize(nPaddle_, 0);
const scalar xMid = xMin + 0.5*(xMax - xMin);
const scalar paddleDy = (yMax - yMin)/scalar(nPaddle_);
for (label paddlei = 0; paddlei < nPaddle_; ++paddlei)
{
xPaddle_[paddlei] = xMid;
yPaddle_[paddlei] = paddlei*paddleDy + yMin + 0.5*paddleDy;
}
// Local face centres
x_ = this->patch().localPoints().component(0);
y_ = this->patch().localPoints().component(1);
z_ = this->patch().localPoints().component(2);
// Local point-to-paddle addressing
pointToPaddle_.setSize(this->patch().size(), -1);
forAll(pointToPaddle_, ppi)
{
pointToPaddle_[ppi] = floor((y_[ppi] - yMin)/(paddleDy+0.01*paddleDy));
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::waveMakerPointPatchVectorField::waveMakerPointPatchVectorField
......@@ -105,10 +151,11 @@ Foam::waveMakerPointPatchVectorField::waveMakerPointPatchVectorField
wavePeriod_(0),
waveHeight_(0),
wavePhase_(0),
waveLength_(0),
waveAngle_(0),
startTime_(0),
rampTime_(1),
secondOrder_(false)
secondOrder_(false),
nPaddle_(0)
{}
......@@ -127,7 +174,7 @@ Foam::waveMakerPointPatchVectorField::waveMakerPointPatchVectorField
wavePeriod_(dict.get<scalar>("wavePeriod")),
waveHeight_(dict.get<scalar>("waveHeight")),
wavePhase_(dict.get<scalar>("wavePhase")),
waveLength_(this->waveLength(initialDepth_, wavePeriod_)),
waveAngle_(dict.get<scalar>("waveAngle")),
startTime_
(
dict.lookupOrDefault<scalar>
......@@ -137,7 +184,8 @@ Foam::waveMakerPointPatchVectorField::waveMakerPointPatchVectorField
)
),
rampTime_(dict.get<scalar>("rampTime")),
secondOrder_(dict.lookupOrDefault<bool>("secondOrder", false))
secondOrder_(dict.lookupOrDefault<bool>("secondOrder", false)),
nPaddle_(readScalar(dict.lookup("nPaddle")))
{
// Create the co-ordinate system
if (mag(n_) < SMALL)
......@@ -151,6 +199,13 @@ Foam::waveMakerPointPatchVectorField::waveMakerPointPatchVectorField
gHat_ = (g() - n_*(n_&g()));
gHat_ /= mag(gHat_);
waveAngle_ *= constant::mathematical::pi/180;
Info << " Initialise geometry... " << endl;
initialiseGeometry();
waterDepthRef_.setSize(nPaddle_, -1);
if (!dict.found("value"))
{
updateCoeffs();
......@@ -174,10 +229,11 @@ Foam::waveMakerPointPatchVectorField::waveMakerPointPatchVectorField
wavePeriod_(ptf.wavePeriod_),
waveHeight_(ptf.waveHeight_),
wavePhase_(ptf.wavePhase_),
waveLength_(ptf.waveLength_),
waveAngle_(ptf.waveAngle_),
startTime_(ptf.startTime_),
rampTime_(ptf.rampTime_),
secondOrder_(ptf.secondOrder_)
secondOrder_(ptf.secondOrder_),
nPaddle_(ptf.nPaddle_)
{}
......@@ -195,10 +251,11 @@ Foam::waveMakerPointPatchVectorField::waveMakerPointPatchVectorField
wavePeriod_(ptf.wavePeriod_),
waveHeight_(ptf.waveHeight_),
wavePhase_(ptf.wavePhase_),
waveLength_(ptf.waveLength_),
waveAngle_(ptf.waveAngle_),
startTime_(ptf.startTime_),
rampTime_(ptf.rampTime_),
secondOrder_(ptf.secondOrder_)
secondOrder_(ptf.secondOrder_),
nPaddle_(ptf.nPaddle_)
{}
......@@ -211,86 +268,174 @@ void Foam::waveMakerPointPatchVectorField::updateCoeffs()
return;
}
if (firstTime==0)
{
// Set the reference water depth
if (initialDepth_ != 0 )
{
forAll(waterDepthRef_,paddlei)
{
waterDepthRef_[paddlei] = initialDepth_;
}
}
else
{
FatalErrorInFunction
<< "initialDepth is not set. Please update "
<< abort(FatalError);
}
Info << " WaterDepth at the wavepaddles = " << waterDepthRef_ << endl;
firstTime = 1;
}
const scalar t = db().time().value() - startTime_;
const scalar waveK = constant::mathematical::twoPi/waveLength_;
const scalar sigma = constant::mathematical::twoPi/wavePeriod_;
scalarField waveLength_(nPaddle_,-1);
const scalar kh = waveK*initialDepth_;
scalarField waveK(nPaddle_,-1);
scalarField waveKx(nPaddle_,-1);
scalarField waveKy(nPaddle_,-1);
forAll (waveK, labelP)
{
waveLength_[labelP] = waveLength (waterDepthRef_[labelP], wavePeriod_);
waveK[labelP] = constant::mathematical::twoPi/waveLength_[labelP];
waveKx[labelP] = waveK[labelP]*cos(waveAngle_);
waveKy[labelP] = waveK[labelP]*sin(waveAngle_);
}
const scalar sigma_ = (2.0*constant::mathematical::pi) / wavePeriod_;
switch (motionType_)
{
case motionTypes::flap:
{
const scalar m1 =
4*sinh(kh)/(sinh(2*kh) + 2*kh)*(sinh(kh) + (1 - cosh(kh))/kh);
scalar motionX = 0.5*waveHeight_/m1*sin(sigma*t);
if (secondOrder_)
{
motionX +=
sqr(waveHeight_)/(16*initialDepth_)
*(3*cosh(kh)/pow3(sinh(kh)) - 2/m1)
*sin(2*sigma*t);
}
const pointField& points = patch().localPoints();
const scalarField dz(-(points & gHat_) - initialDepth_);
scalarField motionX (patch().localPoints().size(),-1);
forAll (points, labelP)
{
const label paddlei = pointToPaddle_[labelP];
scalar phaseTot = waveKx[paddlei]*xPaddle_[paddlei] + waveKy[paddlei]*yPaddle_[paddlei];
if (secondOrder_)
{
scalar m1_ = ((4.0*sinh(waveK[paddlei]*waterDepthRef_[paddlei]))
/ (sinh(2.0*waveK[paddlei]*waterDepthRef_[paddlei])+2.0*waveK[paddlei]*waterDepthRef_[paddlei]))
* ( sinh(waveK[paddlei]*waterDepthRef_[paddlei]) + 1.0/(waveK[paddlei]*waterDepthRef_[paddlei])
* (1.0-cosh(waveK[paddlei]*waterDepthRef_[paddlei])));
motionX[labelP] = waveHeight_/(2.0*m1_)*sin(phaseTot - sigma_*t)
+ pow(waveHeight_,2)/(16.0*waterDepthRef_[paddlei])*(3.0*cosh(waveK[paddlei]*waterDepthRef_[paddlei])
/ pow(sinh(waveK[paddlei]*waterDepthRef_[paddlei]),3) -2.0/m1_)*sin(phaseTot - 2.0*sigma_*t);
motionX[labelP] = timeCoeff(t)
* (1.0+((points[labelP].component(2)-zMinGb_)-waterDepthRef_[paddlei])
/ (waterDepthRef_[paddlei])) * motionX[labelP];
}
else
{
scalar waveBoardStroke_ = (sinh(2.0*waveK[paddlei]*waterDepthRef_[paddlei])
+ 2.0*waveK[paddlei]*waterDepthRef_[paddlei])
/ (4.0*sinh(waveK[paddlei]*waterDepthRef_[paddlei]))
* (1.0 / ( sinh(waveK[paddlei]*waterDepthRef_[paddlei])
+ (1.0-cosh(waveK[paddlei]*waterDepthRef_[paddlei]))
/ (waveK[paddlei]*waterDepthRef_[paddlei]) ) ) * waveHeight_;
motionX[labelP] = timeCoeff(t)
* (1.0+((points[labelP].component(2)-zMinGb_)-waterDepthRef_[paddlei])
/ (waterDepthRef_[paddlei]))*(waveBoardStroke_/2.0)*sin(phaseTot-sigma_*t);
}
}
Field<vector>::operator=
(
n_*timeCoeff(t)*motionX*(1 + dz/initialDepth_)
n_*motionX
);
break;
}
case motionTypes::piston:
{
const scalar m1 = 2*(cosh(2*kh) - 1)/(sinh(2*kh) + 2*kh);
scalar motionX = 0.5*waveHeight_/m1*sin(sigma*t);
if (secondOrder_)
{
motionX +=
sqr(waveHeight_)
/(32*initialDepth_)*(3*cosh(kh)
/pow3(sinh(kh)) - 2/m1);
}
Field<vector>::operator=(n_*timeCoeff(t)*motionX);
const pointField& points = patch().localPoints();
scalarField motionX (patch().localPoints().size(),-1);
forAll (points, labelP)
{
const label paddlei = pointToPaddle_[labelP];
scalar phaseTot = waveKx[paddlei]*xPaddle_[paddlei] + waveKy[paddlei]*yPaddle_[paddlei];
if (secondOrder_)
{
const scalar m1_ = 2.0*(cosh(2.0*waveK[paddlei]*waterDepthRef_[paddlei]) - 1.0)
/ (sinh(2.0*waveK[paddlei]*waterDepthRef_[paddlei])
+ 2.0*waveK[paddlei]*waterDepthRef_[paddlei]);
motionX[labelP] = timeCoeff(t) * (waveHeight_/(2.0*m1_)
* sin(phaseTot - sigma_*t) + pow(waveHeight_,2)
/ (32.0*waterDepthRef_[paddlei])*(3.0*cosh(waveK[paddlei]*waterDepthRef_[paddlei])
/ pow(sinh(waveK[paddlei]*waterDepthRef_[paddlei]),3)-2.0/m1_)
* sin(phaseTot - 2.0*sigma_*t));
}
else
{
scalar waveBoardStroke_ = (sinh(2.0*waveK[paddlei]*waterDepthRef_[paddlei])
+ 2.0*waveK[paddlei]*waterDepthRef_[paddlei])
/ (2.0*(cosh(2.0*waveK[paddlei]*waterDepthRef_[paddlei])
- 1.0)) * waveHeight_;
motionX[labelP] = timeCoeff(t)*(waveBoardStroke_/2.0)
* sin(phaseTot-sigma_*t);
}
}
Field<vector>::operator=
(
n_*motionX
);
break;
}
case motionTypes::solitary:
{
const scalar kappa = sqrt(0.75*waveHeight_/(pow3(initialDepth_)));
const scalar waveCelerity =
sqrt(mag(g())*(initialDepth_ + waveHeight_));
const scalar stroke = sqrt(16.0*waveHeight_*initialDepth_/3.0);
const scalar hr = waveHeight_/initialDepth_;
wavePeriod_ = (2.0/(kappa*waveCelerity))*(3.8 + hr);
const scalar tSolitary = -0.5*wavePeriod_ + t;
// Newton-Rapshon
scalar theta1 = 0;
scalar theta2 = 0;
scalar er = 10000;
const scalar error = 0.001;
while (er > error)
{
theta2 =
theta1
- (theta1 - kappa*waveCelerity*tSolitary + hr*tanh(theta1))
/(1.0 + hr*(1.0/cosh(theta1))*(1.0/cosh(theta1)));
const pointField& points = patch().localPoints();
scalarField motionX (patch().localPoints().size(),-1);
forAll (points, labelP)
{
const label paddlei = pointToPaddle_[labelP];
scalar kappa = sqrt(0.75*waveHeight_/(pow3(waterDepthRef_[paddlei])));
scalar waveCelerity = sqrt(mag(g())*(waterDepthRef_[paddlei] + waveHeight_));
scalar stroke = sqrt(16.0*waveHeight_*waterDepthRef_[paddlei]/3.0);
scalar hr = waveHeight_/waterDepthRef_[paddlei];
wavePeriod_ = (2.0/(kappa*waveCelerity))*(3.8 + hr);
scalar tSolitary = -0.5*wavePeriod_ + t;
// Newton-Rapshon
scalar theta1 = 0;
scalar theta2 = 0;
scalar er = 10000;
const scalar error = 0.001;
while (er > error)
{
theta2 =
theta1
- (theta1 - kappa*waveCelerity*tSolitary + hr*tanh(theta1))
/(1.0 + hr*(1.0/cosh(theta1))*(1.0/cosh(theta1)));
er = mag(theta1 - theta2);
theta1 = theta2;
}
}
scalar motionX =
waveHeight_/(kappa*initialDepth_)*tanh(theta1) + 0.5*stroke;
motionX[labelP] =
waveHeight_/(kappa*waterDepthRef_[paddlei])*tanh(theta1) + 0.5*stroke;
}
Field<vector>::operator=(n_*motionX);
......@@ -318,9 +463,11 @@ void Foam::waveMakerPointPatchVectorField::write(Ostream& os) const
os.writeEntry("wavePeriod", wavePeriod_);
os.writeEntry("waveHeight", waveHeight_);
os.writeEntry("wavePhase", wavePhase_);
os.writeEntry("waveAngle", waveAngle_);
os.writeEntry("startTime", startTime_);
os.writeEntry("rampTime", rampTime_);
os.writeEntry("secondOrder", secondOrder_);
os.writeEntry("nPaddle", nPaddle_);
writeEntry("value", os);
}
......
src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.H
View file @ d0c56621
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2018-2019 IH-Cantabria
Copyright (C) 2018 OpenCFD Ltd.
Copyright (C) 2018
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -136,6 +136,9 @@ class waveMakerPointPatchVectorField
//- Wave phase
scalar wavePhase_;
//wave angle
scalar waveAngle_;
//- Wave length
scalar waveLength_;
......@@ -148,6 +151,54 @@ class waveMakerPointPatchVectorField
//- On/off second order calculation switch
scalar secondOrder_;
// number wavepaddles
scalar nPaddle_;
//- Rotation tensor from global to local system
tensor Rgl_;
//- Rotation tensor from local to global system
tensor Rlg_;
//- Paddle x co-ordinates / [m]
scalarField xPaddle_;
//- Paddle y co-ordinates / [m]
scalarField yPaddle_;
//- Addressing from point patch index to paddle index
labelList pointToPaddle_;
//- Addressing from patch face index to paddle index
labelList faceToPaddle_;
//- Patch face centre x co-ordinates / [m]
scalarField x_;
//- Patch face centre y co-ordinates / [m]
scalarField y_;
//- Patch face centre z co-ordinates / [m]
scalarField z_;
//- Overall (point) span in z-direction / [m]
scalar zSpan_;
//- Minimum z (point) height per patch face / [m]
scalarField zMin_;
//- Global Minimum z (point) / [m]
scalar zMinGb_;
//- Maximum z (point) height per patch face / [m]
scalarField zMax_;
// calculated water depth at the patch
scalarField waterDepthRef_;
//
scalar firstTime = 0;
// Protected Member Functions
......@@ -160,6 +211,8 @@ class waveMakerPointPatchVectorField
//- Return the time scaling coefficient
virtual scalar timeCoeff(const scalar t) const;
//- Initialise
virtual void initialiseGeometry();
public:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment