diff --git a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C
index b568b6b89da4ca26ab55e845d4afa7fa3a0413e7..433407c8444bff7aca75ec00f1773016d3ac6bd6 100644
--- a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C
+++ b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C
@@ -509,132 +509,6 @@ Foam::fv::rotorDiskSource::~rotorDiskSource()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-template<class RhoFieldType>
-void Foam::fv::rotorDiskSource::calculate
-(
- const RhoFieldType& rho,
- const vectorField& U,
- const scalarField& thetag,
- vectorField& force,
- const bool divideVolume,
- const bool output
-) const
-{
- const scalarField& V = mesh_.V();
-
- // 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 = localAxesRotation_->transform(U[cellI], i);
-
- // 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());
- if (alphaEff > mathematical::pi)
- {
- alphaEff -= mathematical::twoPi;
- }
- if (alphaEff < -mathematical::pi)
- {
- alphaEff += mathematical::twoPi;
- }
-
- 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*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
- 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] = localAxesRotation_->invTransform(localForce, i);
-
- 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::fv::rotorDiskSource::addSup
(
fvMatrix<vector>& eqn,
diff --git a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
index 9be5e9e9a8985054abea2869e2f7f27a580ec721..722afcdb0ad9845dcc3cba4a3001a92a84fb9caa 100644
--- a/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
+++ b/src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -26,8 +26,136 @@ License
#include "rotorDiskSource.H"
#include "volFields.H"
+using namespace Foam::constant;
+
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
+template<class RhoFieldType>
+void Foam::fv::rotorDiskSource::calculate
+(
+ const RhoFieldType& rho,
+ const vectorField& U,
+ const scalarField& thetag,
+ vectorField& force,
+ const bool divideVolume,
+ const bool output
+) const
+{
+ const scalarField& V = mesh_.V();
+
+ // 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 = localAxesRotation_->transform(U[cellI], i);
+
+ // 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());
+ if (alphaEff > mathematical::pi)
+ {
+ alphaEff -= mathematical::twoPi;
+ }
+ if (alphaEff < -mathematical::pi)
+ {
+ alphaEff += mathematical::twoPi;
+ }
+
+ 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*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
+ 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] = localAxesRotation_->invTransform(localForce, i);
+
+ 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;
+ }
+}
+
+
template<class Type>
void Foam::fv::rotorDiskSource::writeField
(