From 5a80198bef39013b51b600c4a0c836d04230992f Mon Sep 17 00:00:00 2001
From: andy <andy>
Date: Wed, 23 May 2012 15:11:55 +0100
Subject: [PATCH] ENH: Updates to rotor momentum source
---
.../rotorDiskSource/bladeModel/bladeModel.C | 49 ++++----
.../profileModel/lookup/lookupProfile.C | 37 +++---
.../rotorDiskSource/rotorDiskSource.C | 116 ++++++++----------
.../rotorDiskSource/rotorDiskSource.H | 3 +
.../trimModel/fixed/fixedTrim.C | 5 -
.../trimModel/targetForce/targetForceTrim.C | 5 -
6 files changed, 98 insertions(+), 117 deletions(-)
diff --git a/src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C b/src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C
index a1f76d3e993..af12f6b0ce2 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C
+++ b/src/fieldSources/basicSource/rotorDiskSource/bladeModel/bladeModel.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -47,33 +47,34 @@ void Foam::bladeModel::interpolateWeights
scalar& ddx
) const
{
- scalar x = -GREAT;
- label nElem = values.size();
-
i2 = 0;
- while ((x < xIn) && (i2 < nElem))
- {
- x = values[i2];
- i2++;
- }
+ label nElem = values.size();
- if (i2 == 0)
+ if (nElem == 1)
{
i1 = i2;
ddx = 0.0;
return;
}
- else if (i2 == values.size())
- {
- i2 = values.size() - 1;
- i1 = i2;
- ddx = 0.0;
- return;
- }
else
{
- i1 = i2 - 1;
- ddx = (xIn - values[i1])/(values[i2] - values[i1]);
+ while ((values[i2] < xIn) && (i2 < nElem))
+ {
+ i2++;
+ }
+
+ if (i2 == nElem)
+ {
+ i2 = nElem - 1;
+ i1 = i2;
+ ddx = 0.0;
+ return;
+ }
+ else
+ {
+ i1 = i2 - 1;
+ ddx = (xIn - values[i1])/(values[i2] - values[i1]);
+ }
}
}
@@ -120,14 +121,8 @@ Foam::bladeModel::bladeModel(const dictionary& dict)
}
else
{
- FatalErrorIn
- (
- "Foam::bladeModel::bladeModel"
- "("
- "const dictionary&, "
- "const word&"
- ")"
- ) << "No blade data specified" << exit(FatalError);
+ FatalErrorIn("Foam::bladeModel::bladeModel(const dictionary&)")
+ << "No blade data specified" << exit(FatalError);
}
}
diff --git a/src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C b/src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C
index 1dda6b403b3..93e54bc8efe 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C
+++ b/src/fieldSources/basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C
@@ -49,33 +49,34 @@ void Foam::lookupProfile::interpolateWeights
scalar& ddx
) const
{
- scalar x = -GREAT;
- label nElem = values.size();
-
i2 = 0;
- while ((x < xIn) && (i2 < nElem))
- {
- x = values[i2];
- i2++;
- }
+ label nElem = values.size();
- if (i2 == 0)
+ if (nElem == 1)
{
i1 = i2;
ddx = 0.0;
return;
}
- else if (i2 == nElem)
- {
- i2 = nElem - 1;
- i1 = i2;
- ddx = 0.0;
- return;
- }
else
{
- i1 = i2 - 1;
- ddx = (xIn - values[i1])/(values[i2] - values[i1]);
+ while ((values[i2] < xIn) && (i2 < nElem))
+ {
+ i2++;
+ }
+
+ if (i2 == nElem)
+ {
+ i2 = nElem - 1;
+ i1 = i2;
+ ddx = 0.0;
+ return;
+ }
+ else
+ {
+ i1 = i2 - 1;
+ ddx = (xIn - values[i1])/(values[i2] - values[i1]);
+ }
}
}
diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C
index 20161147df0..6014b623fb9 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C
+++ b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C
@@ -67,6 +67,7 @@ namespace Foam
void Foam::rotorDiskSource::checkData()
{
+ // set inflow type
switch (selectionMode())
{
case smCellSet:
@@ -129,11 +130,10 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
const label nInternalFaces = mesh_.nInternalFaces();
const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
- const vectorField& Sf = mesh_.Sf().internalField();
- const scalarField& magSf = mesh_.magSf().internalField();
+ const vectorField& Sf = mesh_.Sf();
+ const scalarField& magSf = mesh_.magSf();
vector n = vector::zero;
- label nFace = 0;
// calculate cell addressing for selected cells
labelList cellAddr(mesh_.nCells(), -1);
@@ -158,7 +158,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
// correct for parallel running
syncTools::swapBoundaryFaceList(mesh_, nbrFaceCellAddr);
-
// add internal field contributions
for (label faceI = 0; faceI < nInternalFaces; faceI++)
{
@@ -173,7 +172,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
{
area_[own] += magSf[faceI];
n += Sf[faceI];
- nFace++;
}
}
else if ((own == -1) && (nbr != -1))
@@ -184,7 +182,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
{
area_[nbr] += magSf[faceI];
n -= Sf[faceI];
- nFace++;
}
}
}
@@ -210,7 +207,6 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
{
area_[own] += magSfp[j];
n += Sfp[j];
- nFace++;
}
}
}
@@ -222,11 +218,10 @@ void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
const label own = cellAddr[mesh_.faceOwner()[faceI]];
const vector nf = Sfp[j]/magSfp[j];
- if ((own != -1) && ((nf & axis) > 0.8))
+ if ((own != -1) && ((nf & axis) > tol))
{
area_[own] += magSfp[j];
n += Sfp[j];
- nFace++;
}
}
}
@@ -359,32 +354,30 @@ void Foam::rotorDiskSource::constructGeometry()
forAll(cells_, i)
{
- const label cellI = cells_[i];
+ if (area_[i] > ROOTVSMALL)
+ {
+ const label cellI = cells_[i];
- // position in (planar) rotor co-ordinate system
- x_[i] = coordSys_.localPosition(C[cellI]);
+ // position in (planar) rotor co-ordinate system
+ x_[i] = coordSys_.localPosition(C[cellI]);
- // cache max radius
- rMax_ = max(rMax_, x_[i].x());
+ // cache max radius
+ rMax_ = max(rMax_, x_[i].x());
- // swept angle relative to rDir axis [radians] in range 0 -> 2*pi
- scalar psi = x_[i].y();
- if (psi < 0)
- {
- psi += mathematical::twoPi;
- }
+ // swept angle relative to rDir axis [radians] in range 0 -> 2*pi
+ scalar psi = x_[i].y();
- // blade flap angle [radians]
- scalar beta = flap_.beta0 - flap_.beta1*cos(psi) - flap_.beta2*sin(psi);
-
- // determine rotation tensor to convert from planar system into the
- // rotor cone system
- scalar cNeg = cos(-beta);
- scalar sNeg = sin(-beta);
- R_[i] = tensor(cNeg, 0.0, -sNeg, 0.0, 1.0, 0.0, sNeg, 0.0, cNeg);
- scalar cPos = cos(beta);
- scalar sPos = sin(beta);
- invR_[i] = tensor(cPos, 0.0, -sPos, 0.0, 1.0, 0.0, sPos, 0.0, cPos);
+ // blade flap angle [radians]
+ scalar beta =
+ flap_.beta0 - flap_.beta1*cos(psi) - flap_.beta2*sin(psi);
+
+ // determine rotation tensor to convert from planar system into the
+ // rotor cone system
+ scalar c = cos(beta);
+ scalar s = sin(beta);
+ R_[i] = tensor(c, 0, -s, 0, 1, 0, s, 0, c);
+ invR_[i] = R_[i].T();
+ }
}
}
@@ -440,6 +433,7 @@ Foam::rotorDiskSource::rotorDiskSource
:
basicSource(name, modelType, dict, mesh),
rhoName_("none"),
+ rhoRef_(1.0),
omega_(0.0),
nBlades_(0),
inletFlow_(ifLocal),
@@ -477,13 +471,15 @@ void Foam::rotorDiskSource::calculate
const bool output
) const
{
- tmp<volScalarField> trho;
- if (rhoName_ != "none")
- {
- trho = mesh_.lookupObject<volScalarField>(rhoName_);
- }
-
const scalarField& V = mesh_.V();
+ const bool compressible = rhoName_ != "none";
+
+ tmp<volScalarField> trho
+ (
+ compressible
+ ? mesh_.lookupObject<volScalarField>(rhoName_)
+ : volScalarField::null()
+ );
// logging info
@@ -503,17 +499,17 @@ void Foam::rotorDiskSource::calculate
// velocity in local cylindrical reference frame
vector Uc = coordSys_.localVector(U[cellI]);
- // apply correction in local system due to coning
+ // transform from rotor cylindrical into local coning system
Uc = R_[i] & Uc;
// set radial component of velocity to zero
Uc.x() = 0.0;
- // remove blade linear velocity from blade normal component
- Uc.y() -= radius*omega_;
+ // set blade normal component of velocity
+ Uc.y() = radius*omega_ - Uc.y();
// determine blade data for this radius
- // i1 = index of upper bound data point in blade list
+ // i2 = index of upper radius bound data point in blade list
scalar twist = 0.0;
scalar chord = 0.0;
label i1 = -1;
@@ -529,17 +525,7 @@ void Foam::rotorDiskSource::calculate
}
// effective angle of attack
- scalar alphaEff =
- mathematical::pi + atan2(Uc.z(), Uc.y()) - alphaGeom;
-
- if (alphaEff > mathematical::pi)
- {
- alphaEff -= mathematical::twoPi;
- }
- if (alphaEff < -mathematical::pi)
- {
- alphaEff += mathematical::twoPi;
- }
+ scalar alphaEff = alphaGeom - atan2(-Uc.z(), Uc.y());
AOAmin = min(AOAmin, alphaEff);
AOAmax = max(AOAmax, alphaEff);
@@ -564,21 +550,21 @@ void Foam::rotorDiskSource::calculate
// calculate forces perpendicular to blade
scalar pDyn = 0.5*magSqr(Uc);
- if (trho.valid())
+ if (compressible)
{
pDyn *= trho()[cellI];
}
- scalar f = pDyn*chord*nBlades_*area_[i]/mathematical::twoPi;
- vector localForce = vector(0.0, f*Cd, tipFactor*f*Cl);
+ scalar f = pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi;
+ vector localForce = vector(0.0, -f*Cd, tipFactor*f*Cl);
+
+ // accumulate forces
+ dragEff += rhoRef_*localForce.y();
+ liftEff += rhoRef_*localForce.z();
// convert force from local coning system into rotor cylindrical
localForce = invR_[i] & localForce;
- // accumulate forces
- dragEff += localForce.y();
- liftEff += localForce.z();
-
// convert force to global cartesian co-ordinate system
force[cellI] = coordSys_.globalVector(localForce);
@@ -616,6 +602,7 @@ void Foam::rotorDiskSource::addSup(fvMatrix<vector>& eqn, const label fieldI)
}
else
{
+ coeffs_.lookup("rhoRef") >> rhoRef_;
dims.reset(dimForce/dimVolume/dimDensity);
}
@@ -666,6 +653,8 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
coeffs_.lookup("fieldNames") >> fieldNames_;
applied_.setSize(fieldNames_.size(), false);
+
+ // read co-ordinate system/geometry invariant properties
scalar rpm(readScalar(coeffs_.lookup("rpm")));
omega_ = rpm/60.0*mathematical::twoPi;
@@ -683,14 +672,17 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
flap_.beta1 = degToRad(flap_.beta1);
flap_.beta2 = degToRad(flap_.beta2);
- trim_->read(coeffs_);
-
- checkData();
+ // create co-ordinate system
createCoordinateSystem();
+ // read co-odinate system dependent properties
+ checkData();
+
constructGeometry();
+ trim_->read(coeffs_);
+
if (debug)
{
writeField("alphag", trim_->thetag()(), true);
diff --git a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H
index 8e26425e107..a2f856988da 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H
+++ b/src/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H
@@ -148,6 +148,9 @@ protected:
//- Name of density field
word rhoName_;
+ //- Reference density if rhoName = 'none'
+ scalar rhoRef_;
+
//- Rotational speed [rad/s]
// Positive anti-clockwise when looking along -ve lift direction
scalar omega_;
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C b/src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C
index d5f4843d587..28b208c7c94 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C
+++ b/src/fieldSources/basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C
@@ -71,11 +71,6 @@ void Foam::fixedTrim::read(const dictionary& dict)
forAll(thetag_, i)
{
scalar psi = x[i].y();
- if (psi < 0)
- {
- psi += mathematical::twoPi;
- }
-
thetag_[i] = theta0 + theta1c*cos(psi) + theta1s*sin(psi);
}
}
diff --git a/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C b/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C
index 409015d885f..f22226f84c1 100644
--- a/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C
+++ b/src/fieldSources/basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C
@@ -142,11 +142,6 @@ Foam::tmp<Foam::scalarField> Foam::targetForceTrim::thetag() const
forAll(t, i)
{
scalar psi = x[i].y();
- if (psi < 0)
- {
- psi += mathematical::twoPi;
- }
-
t[i] = theta_[0] + theta_[1]*cos(psi) + theta_[2]*sin(psi);
}
--
GitLab