diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C index 3b300e08ae..80e12aae66 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.C @@ -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), diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H index 7565fa17ce..ae2d8a2a47 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSource.H @@ -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 R_; + //- Inverse rotation tensor for flap angle + List invR_; + //- Geometric angle of attack [deg] List alphag_; diff --git a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C index acadb0ec9d..f4d02efd54 100644 --- a/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C +++ b/src/finiteVolume/cfdTools/general/fieldSources/basicSource/rotorDiskSource/rotorDiskSourceTemplates.C @@ -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::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::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::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();