Commit 7c35ea52 authored by Henry's avatar Henry
Browse files

rotorDiskSource: Corrected coordinate transformations

Now the rotor axis may be specified in other than the z-direction
Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1290
parent 2366712a
......@@ -69,14 +69,14 @@ namespace Foam
void Foam::fv::rotorDiskSource::checkData()
{
// set inflow type
// Set inflow type
switch (selectionMode())
{
case smCellSet:
case smCellZone:
case smAll:
{
// set the profile ID for each blade section
// Set the profile ID for each blade section
profiles_.connectBlades(blade_.profileName(), blade_.profileID());
switch (inletFlow_)
{
......@@ -96,7 +96,7 @@ void Foam::fv::rotorDiskSource::checkData()
}
case ifLocal:
{
// do nothing
// Do nothing
break;
}
default:
......@@ -137,7 +137,7 @@ void Foam::fv::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
vector n = vector::zero;
// calculate cell addressing for selected cells
// Calculate cell addressing for selected cells
labelList cellAddr(mesh_.nCells(), -1);
UIndirectList<label>(cellAddr, cells_) = identity(cells_.size());
labelList nbrFaceCellAddr(mesh_.nFaces() - nInternalFaces, -1);
......@@ -157,10 +157,10 @@ void Foam::fv::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
}
}
// correct for parallel running
// Correct for parallel running
syncTools::swapBoundaryFaceList(mesh_, nbrFaceCellAddr);
// add internal field contributions
// Add internal field contributions
for (label faceI = 0; faceI < nInternalFaces; faceI++)
{
const label own = cellAddr[mesh_.faceOwner()[faceI]];
......@@ -189,7 +189,7 @@ void Foam::fv::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
}
// add boundary contributions
// Add boundary contributions
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
......@@ -262,7 +262,7 @@ void Foam::fv::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
void Foam::fv::rotorDiskSource::createCoordinateSystem()
{
// construct the local rotor co-prdinate system
// Construct the local rotor co-prdinate system
vector origin(vector::zero);
vector axis(vector::zero);
vector refDir(vector::zero);
......@@ -274,7 +274,7 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
{
case gmAuto:
{
// determine rotation origin (cell volume weighted)
// Determine rotation origin (cell volume weighted)
scalar sumV = 0.0;
const scalarField& V = mesh_.V();
const vectorField& C = mesh_.C();
......@@ -288,7 +288,7 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
reduce(sumV, sumOp<scalar>());
origin /= sumV;
// determine first radial vector
// Determine first radial vector
vector dx1(vector::zero);
scalar magR = -GREAT;
forAll(cells_, i)
......@@ -304,7 +304,7 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
reduce(dx1, maxMagSqrOp<vector>());
magR = mag(dx1);
// determine second radial vector and cross to determine axis
// Determine second radial vector and cross to determine axis
forAll(cells_, i)
{
const label cellI = cells_[i];
......@@ -321,7 +321,7 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
reduce(axis, maxMagSqrOp<vector>());
axis /= mag(axis);
// correct the axis direction using a point above the rotor
// Correct the axis direction using a point above the rotor
{
vector pointAbove(coeffs_.lookup("pointAbove"));
vector dir = pointAbove - origin;
......@@ -345,7 +345,7 @@ void Foam::fv::rotorDiskSource::createCoordinateSystem()
)
);
// set the face areas and apply correction to calculated axis
// 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);
......@@ -405,20 +405,20 @@ void Foam::fv::rotorDiskSource::constructGeometry()
{
const label cellI = cells_[i];
// position in (planar) rotor co-ordinate system
// Position in (planar) rotor co-ordinate system
x_[i] = coordSys_.localPosition(C[cellI]);
// cache max radius
// Cache max radius
rMax_ = max(rMax_, x_[i].x());
// swept angle relative to rDir axis [radians] in range 0 -> 2*pi
// Swept angle relative to rDir axis [radians] in range 0 -> 2*pi
scalar psi = x_[i].y();
// blade flap angle [radians]
// 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
// Determine rotation tensor to convert from planar system into the
// rotor cone system
scalar c = cos(beta);
scalar s = sin(beta);
......@@ -601,7 +601,7 @@ bool Foam::fv::rotorDiskSource::read(const dictionary& dict)
coeffs_.lookup("fieldNames") >> fieldNames_;
applied_.setSize(fieldNames_.size(), false);
// read co-ordinate system/geometry invariant properties
// Read co-ordinate system/geometry invariant properties
scalar rpm(readScalar(coeffs_.lookup("rpm")));
omega_ = rpm/60.0*mathematical::twoPi;
......@@ -620,10 +620,10 @@ bool Foam::fv::rotorDiskSource::read(const dictionary& dict)
flap_.beta2s = degToRad(flap_.beta2s);
// create co-ordinate system
// Create co-ordinate system
createCoordinateSystem();
// read co-odinate system dependent properties
// Read co-odinate system dependent properties
checkData();
constructGeometry();
......
......@@ -43,7 +43,7 @@ void Foam::fv::rotorDiskSource::calculate
{
const scalarField& V = mesh_.V();
// logging info
// Logging info
scalar dragEff = 0.0;
scalar liftEff = 0.0;
scalar AOAmin = GREAT;
......@@ -57,19 +57,19 @@ void Foam::fv::rotorDiskSource::calculate
const scalar radius = x_[i].x();
// velocity in local cylindrical reference frame
vector Uc = localAxesRotation_->transform(U[cellI], i);
// Transform velocity into local cylindrical reference frame
vector Uc = localAxesRotation_->invTransform(U[cellI], i);
// transform from rotor cylindrical into local coning system
// Transform velocity into local coning system
Uc = R_[i] & Uc;
// set radial component of velocity to zero
// Set radial component of velocity to zero
Uc.x() = 0.0;
// set blade normal component of velocity
// Set blade normal component of velocity
Uc.y() = radius*omega_ - Uc.y();
// determine blade data for this radius
// 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;
......@@ -78,14 +78,14 @@ void Foam::fv::rotorDiskSource::calculate
scalar invDr = 0.0;
blade_.interpolate(radius, twist, chord, i1, i2, invDr);
// flip geometric angle if blade is spinning in reverse (clockwise)
// 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
// Effective angle of attack
scalar alphaEff = alphaGeom - atan2(-Uc.z(), Uc.y());
if (alphaEff > mathematical::pi)
{
......@@ -99,7 +99,7 @@ void Foam::fv::rotorDiskSource::calculate
AOAmin = min(AOAmin, alphaEff);
AOAmax = max(AOAmax, alphaEff);
// determine profile data for this radius and angle of attack
// Determine profile data for this radius and angle of attack
const label profile1 = blade_.profileID()[i1];
const label profile2 = blade_.profileID()[i2];
......@@ -114,24 +114,24 @@ void Foam::fv::rotorDiskSource::calculate
scalar Cd = invDr*(Cd2 - Cd1) + Cd1;
scalar Cl = invDr*(Cl2 - Cl1) + Cl1;
// apply tip effect for blade lift
// Apply tip effect for blade lift
scalar tipFactor = neg(radius/rMax_ - tipEffect_);
// calculate forces perpendicular to blade
// Calculate forces perpendicular to blade
scalar pDyn = 0.5*rho[cellI]*magSqr(Uc);
scalar f = pDyn*chord*nBlades_*area_[i]/radius/mathematical::twoPi;
vector localForce = vector(0.0, -f*Cd, tipFactor*f*Cl);
// accumulate forces
// Accumulate forces
dragEff += rhoRef_*localForce.y();
liftEff += rhoRef_*localForce.z();
// convert force from local coning system into rotor cylindrical
// Transform force from local coning system into rotor cylindrical
localForce = invR_[i] & localForce;
// convert force to global cartesian co-ordinate system
force[cellI] = localAxesRotation_->invTransform(localForce, i);
// Transform force into global Cartesian co-ordinate system
force[cellI] = localAxesRotation_->transform(localForce, i);
if (divideVolume)
{
......
......@@ -103,9 +103,10 @@ void Foam::axesRotation::calcTransform
}
}
// the global->local transformation
// Global->local transformation
Rtr_ = Rtr;
// the local->global transformation
// Local->global transformation
R_ = Rtr.T();
}
......@@ -298,7 +299,8 @@ void Foam::axesRotation::operator=(const dictionary& dict)
}
else if (dict.found("axis") || dict.found("direction"))
{
// let it bomb if only one of axis/direction is defined
// Both "axis" and "direction" are required
// If one is missing the appropriate error message will be generated
order = e3e1;
dict.lookup("axis") >> axis1;
dict.lookup("direction") >> axis2;
......
......@@ -88,6 +88,7 @@ class axesRotation
const axisOrder& order = e3e1
);
public:
//- Runtime type information
......@@ -134,7 +135,7 @@ public:
//- Update the rotation for a list of cells
virtual void updateCells(const polyMesh&, const labelList&)
{
// do nothing
// Do nothing
}
//- Return local-to-global transformation tensor
......@@ -196,14 +197,14 @@ public:
) const;
//- Transform vectorField using transformation tensorField and return
// symmetrical tensorField
// symmetric tensorField
virtual tmp<symmTensorField> transformVector
(
const vectorField& st
) const;
//- Transform vector using transformation tensor and return
// symmetrical tensor
// symmetric tensor
virtual symmTensor transformVector(const vector& st) const;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment