ENH: Updates to rotorDiskSource momentum source

This commit is contained in:
andy 2011-12-07 16:28:10 +00:00
parent 492cf1e760
commit 3c0a78fb1b
3 changed files with 47 additions and 22 deletions

View File

@ -354,25 +354,33 @@ void Foam::rotorDiskSource::constructGeometry()
{
const label cellI = cells_[i];
// position in rotor co-ordinate system
// position in (planar) rotor co-ordinate system
x_[i] = coordSys_.localPosition(C[cellI]);
// cache max radius
rMax_ = max(rMax_, x_[i].x());
// determine swept angle relative to rDir axis
scalar psi = x_[i].y() - rDir.y();
// swept angle relative to rDir axis [radians] in range 0 -> 2*pi
scalar psi = x_[i].y();
if (psi < 0)
{
psi += mathematical::twoPi;
}
// blade flap angle
// blade flap angle [radians]
scalar beta = flap_.beta0 - flap_.beta1*cos(psi) - flap_.beta2*sin(psi);
// determine rotation tensor to convert into the rotor cone plane
scalar c = cos(-beta);
scalar s = sin(-beta);
R_[i] = tensor(1, 0, 0, 0, c, s, 0, -s, c);
// 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);
// geometric angle of attack - not including twist
alphag_[i] = trim_.alphaC - trim_.A*cos(psi) - trim_.B*sin(psi);
// geometric angle of attack - not including twist [radians]
alphag_[i] = trim_.alphaC + trim_.A*cos(psi) + trim_.B*sin(psi);
}
}
@ -439,6 +447,7 @@ Foam::rotorDiskSource::rotorDiskSource
profiles_(coeffs_.subDict("profiles")),
x_(cells_.size(), vector::zero),
R_(cells_.size(), I),
invR_(cells_.size(), I),
alphag_(cells_.size(), 0.0),
area_(cells_.size(), 0.0),
coordSys_(false),

View File

@ -131,17 +131,18 @@ public:
protected:
// Helper structures to encapsulate flap and trim data
// Note: all input in degrees (converted to radians internally)
struct flapData
{
scalar beta0; // coning angle [deg]
scalar beta0; // coning angle
scalar beta1; // lateral flapping coeff
scalar beta2; // longitudinal flapping coeff
};
struct trimData
{
scalar alphaC; // collective pitch angle [deg]
scalar alphaC; // collective pitch angle
scalar A; // lateral cyclic coeff
scalar B; // longitudinal cyclic coeff
};
@ -186,6 +187,9 @@ protected:
//- Rotation tensor for flap angle
List<tensor> R_;
//- Inverse rotation tensor for flap angle
List<tensor> invR_;
//- Geometric angle of attack [deg]
List<scalar> alphag_;

View File

@ -29,7 +29,7 @@ License
#include "unitConversion.H"
#include "volFields.H"
using namespace Foam::constant;
using namespace Foam::constant::mathematical;
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
@ -124,11 +124,11 @@ Foam::tmp<Foam::volVectorField> Foam::rotorDiskSource::calculateForces
const scalar radius = x_[i].x();
// apply correction due to flap in cartesian frame
vector Uc = R_[i] & U[cellI];
// velocity in local cylindrical reference frame
vector Uc = coordSys_.localVector(U[cellI]);
// velocity in local reference frame
Uc = coordSys_.localVector(Uc);
// apply correction in local system due to coning
Uc = R_[i] & Uc;
// set radial component of velocity to zero
Uc.x() = 0.0;
@ -149,7 +149,16 @@ Foam::tmp<Foam::volVectorField> Foam::rotorDiskSource::calculateForces
blade_.interpolate(radius, twist, chord, i1, i2, invDr);
// effective angle of attack
scalar alphaEff = alphag_[i] + twist - atan(Uc.z()/Uc.y());
scalar alphaEff = pi + atan2(Uc.z(), Uc.y()) - (alphag_[i] + twist);
if (alphaEff > pi)
{
alphaEff -= twoPi;
}
if (alphaEff < -pi)
{
alphaEff += twoPi;
}
AOAmin = min(AOAmin, alphaEff);
AOAmax = max(AOAmax, alphaEff);
@ -175,10 +184,13 @@ Foam::tmp<Foam::volVectorField> Foam::rotorDiskSource::calculateForces
tipFactor = 0.0;
}
// calculate forces
const scalar pDyn = 0.5*rho[cellI]*sqr(magUc);
const scalar f = pDyn*chord*nBlades_*area_[i]/(mathematical::twoPi);
const vector localForce(0.0, f*Cd, tipFactor*f*Cl);
// calculate forces perpendicular to blade
scalar pDyn = 0.5*rho[cellI]*sqr(magUc);
scalar f = pDyn*chord*nBlades_*area_[i]/twoPi;
vector localForce = vector(0.0, f*Cd, tipFactor*f*Cl);
// convert force from local coning system into rotor cylindrical
localForce = invR_[i] & localForce;
// accumulate forces
dragEff += localForce.y();