From d0c566217220f7d902cf466fb9bac0722d316fbe Mon Sep 17 00:00:00 2001
From: Gabriel Barajas <>
Date: Sat, 7 Dec 2019 09:19:26 +0000
Subject: [PATCH] INT: Initial IHCantabria update for waveMaker BC
Adds support for paddles to generate 3-D waves
---
.../waveMakerPointPatchVectorField.C | 269 ++++++++++++++----
.../waveMakerPointPatchVectorField.H | 55 +++-
2 files changed, 262 insertions(+), 62 deletions(-)
diff --git a/src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.C b/src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.C
index 301c1e2b21f..609ba79d533 100644
--- a/src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.C
+++ b/src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.C
@@ -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);
}
diff --git a/src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.H b/src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.H
index 2b0030c92ed..42108523eae 100644
--- a/src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.H
+++ b/src/waveModels/derivedPointPatchFields/waveMaker/waveMakerPointPatchVectorField.H
@@ -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:
--
GitLab