Skip to content
Permalink
Browse files

rotorDiskSource: Enable opposite spin in rotorDisk, for fixed thrust …

…direction

Patch contributed by Robert Lee
Resolves patch request https://bugs.openfoam.org/view.php?id=3285
  • Loading branch information...
Henry Weller
Henry Weller committed Jun 4, 2019
1 parent 68e9c8e commit fcbebe9eb37462053879b1439cc1f3800a9a4b20
Showing with 28 additions and 35 deletions.
  1. +28 −35 src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
@@ -45,11 +45,11 @@ void Foam::fv::rotorDiskSource::calculate
const scalarField& V = mesh_.V();

// Logging info
scalar dragEff = 0.0;
scalar liftEff = 0.0;
scalar dragEff = 0;
scalar liftEff = 0;
scalar AOAmin = great;
scalar AOAmax = -great;
scalar powerEff = 0.0;
scalar powerEff = 0;

forAll(cells_, i)
{
@@ -61,42 +61,34 @@ void Foam::fv::rotorDiskSource::calculate

// Transform velocity into local cylindrical reference frame
vector Uc = cylindrical_->invTransform(U[celli], i);
// Uc.x(): radial direction.
// Uc.y(): drag direction.
// Uc.z(): lift / thrust direction.

// Transform velocity into local coning system
Uc = R_[i] & Uc;

// Set radial component of velocity to zero
Uc.x() = 0.0;
Uc.x() = 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;
scalar twist = 0;
scalar chord = 0;
label i1 = -1;
label i2 = -1;
scalar invDr = 0.0;
scalar invDr = 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;
}
const scalar alphaGeom = thetag[i] + twist;

// Effective angle of attack
scalar alphaEff = alphaGeom - atan2(-Uc.z(), Uc.y());
if (alphaEff > mathematical::pi)
{
alphaEff -= mathematical::twoPi;
}
if (alphaEff < -mathematical::pi)
{
alphaEff += mathematical::twoPi;
}
const int rotationSign = (omega_ > 0) ? 1 : -1;
const scalar alphaEff =
alphaGeom - atan2(-Uc.z(), rotationSign*Uc.y());

AOAmin = min(AOAmin, alphaEff);
AOAmax = max(AOAmax, alphaEff);
@@ -105,30 +97,32 @@ void Foam::fv::rotorDiskSource::calculate
const label profile1 = blade_.profileID()[i1];
const label profile2 = blade_.profileID()[i2];

scalar Cd1 = 0.0;
scalar Cl1 = 0.0;
scalar Cd1 = 0;
scalar Cl1 = 0;
profiles_[profile1].Cdl(alphaEff, Cd1, Cl1);

scalar Cd2 = 0.0;
scalar Cl2 = 0.0;
scalar Cd2 = 0;
scalar Cl2 = 0;
profiles_[profile2].Cdl(alphaEff, Cd2, Cl2);

scalar Cd = invDr*(Cd2 - Cd1) + Cd1;
scalar Cl = invDr*(Cl2 - Cl1) + Cl1;
const scalar Cd = invDr*(Cd2 - Cd1) + Cd1;
const scalar Cl = invDr*(Cl2 - Cl1) + Cl1;

// Apply tip effect for blade lift
scalar tipFactor = neg(radius/rMax_ - tipEffect_);
const scalar tipFactor = neg(radius/rMax_ - tipEffect_);

// Calculate forces perpendicular to blade
scalar pDyn = 0.5*rho[celli]*magSqr(Uc);
const scalar pDyn = 0.5*rho[celli]*magSqr(Uc);

const scalar f =
pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi;

scalar f = pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi;
vector localForce = vector(0.0, -f*Cd, tipFactor*f*Cl);
vector localForce = vector(0, rotationSign*-f*Cd, tipFactor*f*Cl);

// Accumulate forces
dragEff += rhoRef_*localForce.y();
liftEff += rhoRef_*localForce.z();
powerEff += rhoRef_ * localForce.y() * radius * omega_;
powerEff += rhoRef_*localForce.y()*radius*omega_;

// Transform force from local coning system into rotor cylindrical
localForce = invR_[i] & localForce;
@@ -193,8 +187,7 @@ void Foam::fv::rotorDiskSource::writeField

if (cells_.size() != values.size())
{
FatalErrorInFunction
<< abort(FatalError);
FatalErrorInFunction << abort(FatalError);
}

forAll(cells_, i)

0 comments on commit fcbebe9

Please sign in to comment.
You can’t perform that action at this time.