Skip to content
Snippets Groups Projects
Commit 1e0358c4 authored by Henry's avatar Henry Committed by Andrew Heather
Browse files

rotorDiskSource: move templated functions into the template .C file

parent 0412aaf4
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSource.C
+ 0
− 126
View file @ 1e0358c4
...@@ -509,132 +509,6 @@ Foam::fv::rotorDiskSource::~rotorDiskSource() ...@@ -509,132 +509,6 @@ Foam::fv::rotorDiskSource::~rotorDiskSource()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * 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 void Foam::fv::rotorDiskSource::addSup
( (
fvMatrix<vector>& eqn, fvMatrix<vector>& eqn,
......
src/fvOptions/sources/derived/rotorDiskSource/rotorDiskSourceTemplates.C
+ 129
− 1
View file @ 1e0358c4
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
...@@ -26,8 +26,136 @@ License ...@@ -26,8 +26,136 @@ License
#include "rotorDiskSource.H" #include "rotorDiskSource.H"
#include "volFields.H" #include "volFields.H"
using namespace Foam::constant;
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * 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> template<class Type>
void Foam::fv::rotorDiskSource::writeField void Foam::fv::rotorDiskSource::writeField
( (
......
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