mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
695 lines
19 KiB
C
695 lines
19 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2011-2012 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 <http://www.gnu.org/licenses/>.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "rotorDiskSource.H"
|
|
#include "addToRunTimeSelectionTable.H"
|
|
#include "mathematicalConstants.H"
|
|
#include "trimModel.H"
|
|
#include "unitConversion.H"
|
|
#include "fvMatrices.H"
|
|
#include "syncTools.H"
|
|
|
|
using namespace Foam::constant;
|
|
|
|
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
|
|
|
|
namespace Foam
|
|
{
|
|
defineTypeNameAndDebug(rotorDiskSource, 0);
|
|
addToRunTimeSelectionTable(basicSource, rotorDiskSource, dictionary);
|
|
|
|
template<> const char* NamedEnum<rotorDiskSource::geometryModeType, 2>::
|
|
names[] =
|
|
{
|
|
"auto",
|
|
"specified"
|
|
};
|
|
|
|
const NamedEnum<rotorDiskSource::geometryModeType, 2>
|
|
rotorDiskSource::geometryModeTypeNames_;
|
|
|
|
template<> const char* NamedEnum<rotorDiskSource::inletFlowType, 3>::
|
|
names[] =
|
|
{
|
|
"fixed",
|
|
"surfaceNormal",
|
|
"local"
|
|
};
|
|
|
|
const NamedEnum<rotorDiskSource::inletFlowType, 3>
|
|
rotorDiskSource::inletFlowTypeNames_;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
|
|
|
|
void Foam::rotorDiskSource::checkData()
|
|
{
|
|
// set inflow type
|
|
switch (selectionMode())
|
|
{
|
|
case smCellSet:
|
|
case smCellZone:
|
|
case smAll:
|
|
{
|
|
// set the profile ID for each blade section
|
|
profiles_.connectBlades(blade_.profileName(), blade_.profileID());
|
|
switch (inletFlow_)
|
|
{
|
|
case ifFixed:
|
|
{
|
|
coeffs_.lookup("inletVelocity") >> inletVelocity_;
|
|
break;
|
|
}
|
|
case ifSurfaceNormal:
|
|
{
|
|
scalar UIn
|
|
(
|
|
readScalar(coeffs_.lookup("inletNormalVelocity"))
|
|
);
|
|
inletVelocity_ = -coordSys_.e3()*UIn;
|
|
break;
|
|
}
|
|
case ifLocal:
|
|
{
|
|
// do nothing
|
|
break;
|
|
}
|
|
default:
|
|
{
|
|
FatalErrorIn("void rotorDiskSource::checkData()")
|
|
<< "Unknown inlet velocity type" << abort(FatalError);
|
|
}
|
|
}
|
|
|
|
|
|
break;
|
|
}
|
|
default:
|
|
{
|
|
FatalErrorIn("void rotorDiskSource::checkData()")
|
|
<< "Source cannot be used with '"
|
|
<< selectionModeTypeNames_[selectionMode()]
|
|
<< "' mode. Please use one of: " << nl
|
|
<< selectionModeTypeNames_[smCellSet] << nl
|
|
<< selectionModeTypeNames_[smCellZone] << nl
|
|
<< selectionModeTypeNames_[smAll]
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
|
|
{
|
|
area_ = 0.0;
|
|
|
|
static const scalar tol = 0.8;
|
|
|
|
const label nInternalFaces = mesh_.nInternalFaces();
|
|
const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
|
|
const vectorField& Sf = mesh_.Sf();
|
|
const scalarField& magSf = mesh_.magSf();
|
|
|
|
vector n = vector::zero;
|
|
|
|
// calculate cell addressing for selected cells
|
|
labelList cellAddr(mesh_.nCells(), -1);
|
|
UIndirectList<label>(cellAddr, cells_) = identity(cells_.size());
|
|
labelList nbrFaceCellAddr(mesh_.nFaces() - nInternalFaces, -1);
|
|
forAll(pbm, patchI)
|
|
{
|
|
const polyPatch& pp = pbm[patchI];
|
|
|
|
if (pp.coupled())
|
|
{
|
|
forAll(pp, i)
|
|
{
|
|
label faceI = pp.start() + i;
|
|
label nbrFaceI = faceI - nInternalFaces;
|
|
label own = mesh_.faceOwner()[faceI];
|
|
nbrFaceCellAddr[nbrFaceI] = cellAddr[own];
|
|
}
|
|
}
|
|
}
|
|
|
|
// correct for parallel running
|
|
syncTools::swapBoundaryFaceList(mesh_, nbrFaceCellAddr);
|
|
|
|
// add internal field contributions
|
|
for (label faceI = 0; faceI < nInternalFaces; faceI++)
|
|
{
|
|
const label own = cellAddr[mesh_.faceOwner()[faceI]];
|
|
const label nbr = cellAddr[mesh_.faceNeighbour()[faceI]];
|
|
|
|
if ((own != -1) && (nbr == -1))
|
|
{
|
|
vector nf = Sf[faceI]/magSf[faceI];
|
|
|
|
if ((nf & axis) > tol)
|
|
{
|
|
area_[own] += magSf[faceI];
|
|
n += Sf[faceI];
|
|
}
|
|
}
|
|
else if ((own == -1) && (nbr != -1))
|
|
{
|
|
vector nf = Sf[faceI]/magSf[faceI];
|
|
|
|
if ((-nf & axis) > tol)
|
|
{
|
|
area_[nbr] += magSf[faceI];
|
|
n -= Sf[faceI];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// add boundary contributions
|
|
forAll(pbm, patchI)
|
|
{
|
|
const polyPatch& pp = pbm[patchI];
|
|
const vectorField& Sfp = mesh_.Sf().boundaryField()[patchI];
|
|
const scalarField& magSfp = mesh_.magSf().boundaryField()[patchI];
|
|
|
|
if (pp.coupled())
|
|
{
|
|
forAll(pp, j)
|
|
{
|
|
const label faceI = pp.start() + j;
|
|
const label own = cellAddr[mesh_.faceOwner()[faceI]];
|
|
const bool nbr = nbrFaceCellAddr[faceI - nInternalFaces];
|
|
const vector nf = Sfp[j]/magSfp[j];
|
|
|
|
if ((own != -1) && (nbr == -1) && ((nf & axis) > tol))
|
|
{
|
|
area_[own] += magSfp[j];
|
|
n += Sfp[j];
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
forAll(pp, j)
|
|
{
|
|
const label faceI = pp.start() + j;
|
|
const label own = cellAddr[mesh_.faceOwner()[faceI]];
|
|
const vector nf = Sfp[j]/magSfp[j];
|
|
|
|
if ((own != -1) && ((nf & axis) > tol))
|
|
{
|
|
area_[own] += magSfp[j];
|
|
n += Sfp[j];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if (correct)
|
|
{
|
|
reduce(n, sumOp<vector>());
|
|
axis = n/mag(n);
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::rotorDiskSource::createCoordinateSystem()
|
|
{
|
|
// construct the local rotor co-prdinate system
|
|
vector origin(vector::zero);
|
|
vector axis(vector::zero);
|
|
vector refDir(vector::zero);
|
|
|
|
geometryModeType gm =
|
|
geometryModeTypeNames_.read(coeffs_.lookup("geometryMode"));
|
|
|
|
switch (gm)
|
|
{
|
|
case gmAuto:
|
|
{
|
|
// determine rotation origin (cell volume weighted)
|
|
scalar sumV = 0.0;
|
|
const scalarField& V = mesh_.V();
|
|
const vectorField& C = mesh_.C();
|
|
forAll(cells_, i)
|
|
{
|
|
const label cellI = cells_[i];
|
|
sumV += V[cellI];
|
|
origin += V[cellI]*C[cellI];
|
|
}
|
|
reduce(origin, sumOp<vector>());
|
|
reduce(sumV, sumOp<scalar>());
|
|
origin /= sumV;
|
|
|
|
// determine first radial vector
|
|
vector dx1(vector::zero);
|
|
scalar magR = -GREAT;
|
|
forAll(cells_, i)
|
|
{
|
|
const label cellI = cells_[i];
|
|
vector test = C[cellI] - origin;
|
|
if (mag(test) > magR)
|
|
{
|
|
dx1 = test;
|
|
magR = mag(test);
|
|
}
|
|
}
|
|
reduce(dx1, maxMagSqrOp<vector>());
|
|
magR = mag(dx1);
|
|
|
|
// determine second radial vector and cross to determine axis
|
|
forAll(cells_, i)
|
|
{
|
|
const label cellI = cells_[i];
|
|
vector dx2 = C[cellI] - origin;
|
|
if (mag(dx2) > 0.5*magR)
|
|
{
|
|
axis = dx1 ^ dx2;
|
|
if (mag(axis) > SMALL)
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
reduce(axis, maxMagSqrOp<vector>());
|
|
axis /= mag(axis);
|
|
|
|
// correct the axis direction using a point above the rotor
|
|
{
|
|
vector pointAbove(coeffs_.lookup("pointAbove"));
|
|
vector dir = pointAbove - origin;
|
|
dir /= mag(dir);
|
|
if ((dir & axis) < 0)
|
|
{
|
|
axis *= -1.0;
|
|
}
|
|
}
|
|
|
|
coeffs_.lookup("refDirection") >> refDir;
|
|
|
|
// set the face areas and apply correction to calculated axis
|
|
// e.g. if cellZone is more than a single layer in thickness
|
|
setFaceArea(axis, true);
|
|
|
|
break;
|
|
}
|
|
case gmSpecified:
|
|
{
|
|
coeffs_.lookup("origin") >> origin;
|
|
coeffs_.lookup("axis") >> axis;
|
|
coeffs_.lookup("refDirection") >> refDir;
|
|
|
|
setFaceArea(axis, false);
|
|
|
|
break;
|
|
}
|
|
default:
|
|
{
|
|
FatalErrorIn("rotorDiskSource::createCoordinateSystem()")
|
|
<< "Unknown geometryMode " << geometryModeTypeNames_[gm]
|
|
<< ". Available geometry modes include "
|
|
<< geometryModeTypeNames_ << exit(FatalError);
|
|
}
|
|
}
|
|
|
|
coordSys_ = cylindricalCS("rotorCoordSys", origin, axis, refDir, false);
|
|
|
|
const scalar sumArea = gSum(area_);
|
|
const scalar diameter = Foam::sqrt(4.0*sumArea/mathematical::pi);
|
|
Info<< " Rotor gometry:" << nl
|
|
<< " - disk diameter = " << diameter << nl
|
|
<< " - disk area = " << sumArea << nl
|
|
<< " - origin = " << coordSys_.origin() << nl
|
|
<< " - r-axis = " << coordSys_.e1() << nl
|
|
<< " - psi-axis = " << coordSys_.e2() << nl
|
|
<< " - z-axis = " << coordSys_.e3() << endl;
|
|
}
|
|
|
|
|
|
void Foam::rotorDiskSource::constructGeometry()
|
|
{
|
|
const vectorField& C = mesh_.C();
|
|
|
|
forAll(cells_, i)
|
|
{
|
|
if (area_[i] > ROOTVSMALL)
|
|
{
|
|
const label cellI = cells_[i];
|
|
|
|
// position in (planar) rotor co-ordinate system
|
|
x_[i] = coordSys_.localPosition(C[cellI]);
|
|
|
|
// 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();
|
|
|
|
// blade flap angle [radians]
|
|
scalar beta =
|
|
flap_.beta0 - flap_.beta1c*cos(psi) - flap_.beta2s*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();
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
Foam::tmp<Foam::vectorField> Foam::rotorDiskSource::inflowVelocity
|
|
(
|
|
const volVectorField& U
|
|
) const
|
|
{
|
|
switch (inletFlow_)
|
|
{
|
|
case ifFixed:
|
|
case ifSurfaceNormal:
|
|
{
|
|
return tmp<vectorField>
|
|
(
|
|
new vectorField(mesh_.nCells(), inletVelocity_)
|
|
);
|
|
|
|
break;
|
|
}
|
|
case ifLocal:
|
|
{
|
|
return U.internalField();
|
|
|
|
break;
|
|
}
|
|
default:
|
|
{
|
|
FatalErrorIn
|
|
(
|
|
"Foam::tmp<Foam::vectorField> "
|
|
"Foam::rotorDiskSource::inflowVelocity"
|
|
"(const volVectorField&) const"
|
|
) << "Unknown inlet flow specification" << abort(FatalError);
|
|
}
|
|
}
|
|
|
|
return tmp<vectorField>(new vectorField(mesh_.nCells(), vector::zero));
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
|
|
|
|
Foam::rotorDiskSource::rotorDiskSource
|
|
(
|
|
const word& name,
|
|
const word& modelType,
|
|
const dictionary& dict,
|
|
const fvMesh& mesh
|
|
|
|
)
|
|
:
|
|
basicSource(name, modelType, dict, mesh),
|
|
rhoName_("none"),
|
|
rhoRef_(1.0),
|
|
omega_(0.0),
|
|
nBlades_(0),
|
|
inletFlow_(ifLocal),
|
|
inletVelocity_(vector::zero),
|
|
tipEffect_(1.0),
|
|
flap_(),
|
|
x_(cells_.size(), vector::zero),
|
|
R_(cells_.size(), I),
|
|
invR_(cells_.size(), I),
|
|
area_(cells_.size(), 0.0),
|
|
coordSys_(false),
|
|
rMax_(0.0),
|
|
trim_(trimModel::New(*this, coeffs_)),
|
|
blade_(coeffs_.subDict("blade")),
|
|
profiles_(coeffs_.subDict("profiles"))
|
|
{
|
|
read(dict);
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
|
|
|
Foam::rotorDiskSource::~rotorDiskSource()
|
|
{}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
void Foam::rotorDiskSource::calculate
|
|
(
|
|
const vectorField& U,
|
|
const scalarField& thetag,
|
|
vectorField& force,
|
|
const bool divideVolume,
|
|
const bool output
|
|
) const
|
|
{
|
|
const scalarField& V = mesh_.V();
|
|
const bool compressible = this->compressible();
|
|
tmp<volScalarField> trho(rho());
|
|
|
|
// logging info
|
|
scalar dragEff = 0.0;
|
|
scalar liftEff = 0.0;
|
|
scalar AOAmin = GREAT;
|
|
scalar AOAmax = -GREAT;
|
|
|
|
forAll(cells_, i)
|
|
{
|
|
if (area_[i] > ROOTVSMALL)
|
|
{
|
|
const label cellI = cells_[i];
|
|
|
|
const scalar radius = x_[i].x();
|
|
|
|
// velocity in local cylindrical reference frame
|
|
vector Uc = coordSys_.localVector(U[cellI]);
|
|
|
|
// transform from rotor cylindrical into local coning system
|
|
Uc = R_[i] & Uc;
|
|
|
|
// set radial component of velocity to zero
|
|
Uc.x() = 0.0;
|
|
|
|
// set blade normal component of velocity
|
|
Uc.y() = radius*omega_ - Uc.y();
|
|
|
|
// determine blade data for this radius
|
|
// i2 = index of upper radius bound data point in blade list
|
|
scalar twist = 0.0;
|
|
scalar chord = 0.0;
|
|
label i1 = -1;
|
|
label i2 = -1;
|
|
scalar invDr = 0.0;
|
|
blade_.interpolate(radius, twist, chord, i1, i2, invDr);
|
|
|
|
// flip geometric angle if blade is spinning in reverse (clockwise)
|
|
scalar alphaGeom = thetag[i] + twist;
|
|
if (omega_ < 0)
|
|
{
|
|
alphaGeom = mathematical::pi - alphaGeom;
|
|
}
|
|
|
|
// effective angle of attack
|
|
scalar alphaEff = alphaGeom - atan2(-Uc.z(), Uc.y());
|
|
|
|
AOAmin = min(AOAmin, alphaEff);
|
|
AOAmax = max(AOAmax, alphaEff);
|
|
|
|
// determine profile data for this radius and angle of attack
|
|
const label profile1 = blade_.profileID()[i1];
|
|
const label profile2 = blade_.profileID()[i2];
|
|
|
|
scalar Cd1 = 0.0;
|
|
scalar Cl1 = 0.0;
|
|
profiles_[profile1].Cdl(alphaEff, Cd1, Cl1);
|
|
|
|
scalar Cd2 = 0.0;
|
|
scalar Cl2 = 0.0;
|
|
profiles_[profile2].Cdl(alphaEff, Cd2, Cl2);
|
|
|
|
scalar Cd = invDr*(Cd2 - Cd1) + Cd1;
|
|
scalar Cl = invDr*(Cl2 - Cl1) + Cl1;
|
|
|
|
// apply tip effect for blade lift
|
|
scalar tipFactor = neg(radius/rMax_ - tipEffect_);
|
|
|
|
// calculate forces perpendicular to blade
|
|
scalar pDyn = 0.5*magSqr(Uc);
|
|
if (compressible)
|
|
{
|
|
pDyn *= trho()[cellI];
|
|
}
|
|
|
|
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;
|
|
|
|
// convert force to global cartesian co-ordinate system
|
|
force[cellI] = coordSys_.globalVector(localForce);
|
|
|
|
if (divideVolume)
|
|
{
|
|
force[cellI] /= V[cellI];
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
if (output)
|
|
{
|
|
reduce(AOAmin, minOp<scalar>());
|
|
reduce(AOAmax, maxOp<scalar>());
|
|
reduce(dragEff, sumOp<scalar>());
|
|
reduce(liftEff, sumOp<scalar>());
|
|
|
|
Info<< type() << " output:" << nl
|
|
<< " min/max(AOA) = " << radToDeg(AOAmin) << ", "
|
|
<< radToDeg(AOAmax) << nl
|
|
<< " Effective drag = " << dragEff << nl
|
|
<< " Effective lift = " << liftEff << endl;
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::rotorDiskSource::addSup(fvMatrix<vector>& eqn, const label fieldI)
|
|
{
|
|
dimensionSet dims = dimless;
|
|
if (eqn.dimensions() == dimForce)
|
|
{
|
|
coeffs_.lookup("rhoName") >> rhoName_;
|
|
dims.reset(dimForce/dimVolume);
|
|
}
|
|
else
|
|
{
|
|
coeffs_.lookup("rhoRef") >> rhoRef_;
|
|
dims.reset(dimForce/dimVolume/dimDensity);
|
|
}
|
|
|
|
volVectorField force
|
|
(
|
|
IOobject
|
|
(
|
|
"rotorForce",
|
|
mesh_.time().timeName(),
|
|
mesh_,
|
|
IOobject::NO_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
mesh_,
|
|
dimensionedVector("zero", dims, vector::zero)
|
|
);
|
|
|
|
const volVectorField& U = eqn.psi();
|
|
|
|
const vectorField Uin(inflowVelocity(U));
|
|
|
|
trim_->correct(Uin, force);
|
|
|
|
calculate(Uin, trim_->thetag(), force);
|
|
|
|
|
|
// add source to rhs of eqn
|
|
eqn -= force;
|
|
|
|
if (mesh_.time().outputTime())
|
|
{
|
|
force.write();
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::rotorDiskSource::writeData(Ostream& os) const
|
|
{
|
|
os << indent << name_ << endl;
|
|
dict_.write(os);
|
|
}
|
|
|
|
|
|
bool Foam::rotorDiskSource::read(const dictionary& dict)
|
|
{
|
|
if (basicSource::read(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;
|
|
|
|
coeffs_.lookup("nBlades") >> nBlades_;
|
|
|
|
inletFlow_ = inletFlowTypeNames_.read(coeffs_.lookup("inletFlowType"));
|
|
|
|
coeffs_.lookup("tipEffect") >> tipEffect_;
|
|
|
|
const dictionary& flapCoeffs(coeffs_.subDict("flapCoeffs"));
|
|
flapCoeffs.lookup("beta0") >> flap_.beta0;
|
|
flapCoeffs.lookup("beta1c") >> flap_.beta1c;
|
|
flapCoeffs.lookup("beta2s") >> flap_.beta2s;
|
|
flap_.beta0 = degToRad(flap_.beta0);
|
|
flap_.beta1c = degToRad(flap_.beta1c);
|
|
flap_.beta2s = degToRad(flap_.beta2s);
|
|
|
|
|
|
// create co-ordinate system
|
|
createCoordinateSystem();
|
|
|
|
// read co-odinate system dependent properties
|
|
checkData();
|
|
|
|
constructGeometry();
|
|
|
|
trim_->read(coeffs_);
|
|
|
|
if (debug)
|
|
{
|
|
writeField("thetag", trim_->thetag()(), true);
|
|
writeField("faceArea", area_, true);
|
|
}
|
|
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|