Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
Henry
2012-04-16 11:47:34 +01:00
48 changed files with 1583 additions and 319 deletions

View File

@ -152,6 +152,16 @@ dimensioned<Type>::dimensioned
}
template <class Type>
dimensioned<Type>::dimensioned
()
:
name_("undefined"),
dimensions_(dimless),
value_(pTraits<Type>::zero)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template <class Type>

View File

@ -107,6 +107,9 @@ public:
//- Construct from an Istream with a given name and dimensions
dimensioned(const word&, const dimensionSet&, Istream&);
//- Null constructor
dimensioned();
//- Construct from dictionary, with default value.
static dimensioned<Type> lookupOrDefault
(

View File

@ -180,7 +180,7 @@ bool Foam::mergePoints
const bool verbose,
labelList& pointMap,
List<Type>& newPoints,
const Type& origin = Type::zero
const Type& origin
)
{
label nUnique = mergePoints

View File

@ -163,6 +163,22 @@ public:
return TableBase<Type>::integrate(x1, x2);
}
//- Return dimensioned constant value
virtual dimensioned<Type> dimValue(const scalar x) const
{
return TableBase<Type>::dimValue(x);
}
//- Integrate between two values and return dimensioned type
virtual dimensioned<Type> dimIntegrate
(
const scalar x1,
const scalar x2
) const
{
return TableBase<Type>::dimIntegrate(x1, x2);
}
// I/O

View File

@ -34,9 +34,27 @@ Foam::CompatibilityConstant<Type>::CompatibilityConstant
)
:
DataEntry<Type>(entryName),
value_(pTraits<Type>::zero)
value_(pTraits<Type>::zero),
dimensions_(dimless)
{
dict.lookup(entryName) >> value_;
Istream& is(dict.lookup(entryName));
token firstToken(is);
if (firstToken.isWord())
{
token nextToken(is);
if (nextToken == token::BEGIN_SQR)
{
is.putBack(nextToken);
is >> dimensions_;
is >> value_;
}
}
else
{
is.putBack(firstToken);
is >> value_;
}
}
@ -47,10 +65,10 @@ Foam::CompatibilityConstant<Type>::CompatibilityConstant
)
:
DataEntry<Type>(cnst),
value_(cnst.value_)
value_(cnst.value_),
dimensions_(cnst.dimensions_)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template<class Type>
@ -78,6 +96,23 @@ Type Foam::CompatibilityConstant<Type>::integrate
}
template<class Type>
Foam::dimensioned<Type> Foam::CompatibilityConstant<Type>::
dimValue(const scalar x) const
{
return dimensioned<Type>("dimensionedValue", dimensions_, value_);
}
template<class Type>
Foam::dimensioned<Type> Foam::CompatibilityConstant<Type>::dimIntegrate
(
const scalar x1, const scalar x2
) const
{
return dimensioned<Type>("dimensionedValue", dimensions_, (x2-x1)*value_);
}
// * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
#include "CompatibilityConstantIO.C"

View File

@ -42,6 +42,7 @@ SourceFiles
#define CompatibilityConstant_H
#include "DataEntry.H"
#include "dimensionSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -68,6 +69,9 @@ class CompatibilityConstant
//- Constant value
Type value_;
//- The dimension set
dimensionSet dimensions_;
// Private Member Functions
@ -111,6 +115,16 @@ public:
//- Integrate between two values
Type integrate(const scalar x1, const scalar x2) const;
//- Return dimensioned constant value
dimensioned<Type> dimValue(const scalar) const;
//- Integrate between two values and return dimensioned type
dimensioned<Type> dimIntegrate
(
const scalar x1,
const scalar x2
) const;
// I/O

View File

@ -31,12 +31,27 @@ template<class Type>
Foam::Constant<Type>::Constant(const word& entryName, const dictionary& dict)
:
DataEntry<Type>(entryName),
value_(pTraits<Type>::zero)
value_(pTraits<Type>::zero),
dimensions_(dimless)
{
Istream& is(dict.lookup(entryName));
word entryType(is);
is >> value_;
token firstToken(is);
if (firstToken.isWord())
{
token nextToken(is);
if (nextToken == token::BEGIN_SQR)
{
is.putBack(nextToken);
is >> dimensions_;
is >> value_;
}
}
else
{
is.putBack(firstToken);
is >> value_;
}
}
@ -44,7 +59,8 @@ template<class Type>
Foam::Constant<Type>::Constant(const Constant<Type>& cnst)
:
DataEntry<Type>(cnst),
value_(cnst.value_)
value_(cnst.value_),
dimensions_(cnst.dimensions_)
{}
@ -71,6 +87,22 @@ Type Foam::Constant<Type>::integrate(const scalar x1, const scalar x2) const
}
template<class Type>
Foam::dimensioned<Type> Foam::Constant<Type>::dimValue(const scalar x) const
{
return dimensioned<Type>("dimensionedValue", dimensions_, value_);
}
template<class Type>
Foam::dimensioned<Type> Foam::Constant<Type>::dimIntegrate
(
const scalar x1, const scalar x2
) const
{
return dimensioned<Type>("dimensionedValue", dimensions_, (x2-x1)*value_);
}
// * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
#include "ConstantIO.C"

View File

@ -41,6 +41,7 @@ SourceFiles
#define Constant_H
#include "DataEntry.H"
#include "dimensionSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -67,6 +68,9 @@ class Constant
//- Constant value
Type value_;
//- The dimension set
dimensionSet dimensions_;
// Private Member Functions
@ -107,6 +111,16 @@ public:
//- Integrate between two values
Type integrate(const scalar x1, const scalar x2) const;
//- Return dimensioned constant value
dimensioned<Type> dimValue(const scalar) const;
//- Integrate between two values and return dimensioned type
dimensioned<Type> dimIntegrate
(
const scalar x1,
const scalar x2
) const;
// I/O

View File

@ -76,6 +76,22 @@ Type Foam::DataEntry<Type>::value(const scalar x) const
}
template<class Type>
Type Foam::DataEntry<Type>::integrate(const scalar x1, const scalar x2) const
{
notImplemented
(
"Type Foam::DataEntry<Type>::integrate"
"("
"const scalar, "
"const scalar"
") const"
);
return pTraits<Type>::zero;
}
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::DataEntry<Type>::value
(
@ -93,22 +109,6 @@ Foam::tmp<Foam::Field<Type> > Foam::DataEntry<Type>::value
}
template<class Type>
Type Foam::DataEntry<Type>::integrate(const scalar x1, const scalar x2) const
{
notImplemented
(
"Type Foam::DataEntry<Type>::integrate"
"("
"const scalar, "
"const scalar"
") const"
);
return pTraits<Type>::zero;
}
template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::DataEntry<Type>::integrate
(
@ -127,6 +127,90 @@ Foam::tmp<Foam::Field<Type> > Foam::DataEntry<Type>::integrate
}
template<class Type>
Foam::dimensioned<Type> Foam::DataEntry<Type>::dimValue(const scalar x) const
{
notImplemented
(
"dimensioned<Type> Foam::DataEntry<dimensioned<Type> >::dimValue"
"(const scalar) const"
);
return dimensioned<Type>("zero", dimless, pTraits<Type>::zero);
}
template<class Type>
Foam::dimensioned<Type> Foam::DataEntry<Type>::dimIntegrate
(
const scalar x1,
const scalar x2
) const
{
notImplemented
(
"dimensioned<Type> Foam::DataEntry<Type>::dimIntegrate"
"("
"const scalar, "
"const scalar"
") const"
);
return dimensioned<Type>("zero", dimless, pTraits<Type>::zero);
}
template<class Type>
Foam::tmp<Foam::Field<Foam::dimensioned<Type> > >
Foam::DataEntry<Type>::dimValue
(
const scalarField& x
) const
{
tmp<Field<dimensioned<Type> > > tfld
(
new Field<dimensioned<Type> >
(
x.size(),
dimensioned<Type>("zero", dimless, pTraits<Type>::zero)
)
);
Field<dimensioned<Type> >& fld = tfld();
forAll(x, i)
{
fld[i] = this->dimValue(x[i]);
}
return tfld;
}
template<class Type>
Foam::tmp<Foam::Field<Foam::dimensioned<Type> > >
Foam::DataEntry<Type>::dimIntegrate
(
const scalarField& x1,
const scalarField& x2
) const
{
tmp<Field<dimensioned<Type> > > tfld
(
new Field<dimensioned<Type> >(x1.size())
);
Field<dimensioned<Type> >& fld = tfld();
forAll(x1, i)
{
fld[i] = this->dimIntegrate(x1[i], x2[i]);
}
return tfld;
}
// * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
#include "DataEntryIO.C"

View File

@ -41,6 +41,7 @@ SourceFiles
#include "dictionary.H"
#include "Field.H"
#include "dimensionedType.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -141,23 +142,49 @@ public:
virtual void convertTimeBase(const Time& t);
// Evaluation
public:
//- Return value as a function of (scalar) independent variable
virtual Type value(const scalar x) const;
//- Return value as a function of (scalar) independent variable
virtual tmp<Field<Type> > value(const scalarField& x) const;
// Evaluation
//- Integrate between two (scalar) values
virtual Type integrate(const scalar x1, const scalar x2) const;
//- Return value as a function of (scalar) independent variable
virtual Type value(const scalar x) const;
//- Integrate between two (scalar) values
virtual tmp<Field<Type> > integrate
(
const scalarField& x1,
const scalarField& x2
) const;
//- Return value as a function of (scalar) independent variable
virtual tmp<Field<Type> > value(const scalarField& x) const;
//- Integrate between two (scalar) values
virtual Type integrate(const scalar x1, const scalar x2) const;
//- Integrate between two (scalar) values
virtual tmp<Field<Type> > integrate
(
const scalarField& x1,
const scalarField& x2
) const;
//- Return dimensioned type
virtual dimensioned<Type> dimValue(const scalar x) const;
//- Return dimensioned type as a function of (scalar)
virtual tmp<Field<dimensioned<Type> > > dimValue
(
const scalarField& x
) const;
//- Integrate between two scalars and returns a dimensioned type
virtual dimensioned<Type> dimIntegrate
(
const scalar x1,
const scalar x2
) const;
//- Integrate between two scalars and returns list of dimensioned type
virtual tmp<Field<dimensioned<Type> > > dimIntegrate
(
const scalarField& x1,
const scalarField& x2
) const;
// I/O

View File

@ -41,7 +41,15 @@ Foam::autoPtr<Foam::DataEntry<Type> > Foam::DataEntry<Type>::New
word DataEntryType;
if (firstToken.isWord())
{
DataEntryType = firstToken.wordToken();
// Dimensioned type default compatibility
if (firstToken.wordToken() == entryName)
{
DataEntryType = "CompatibilityConstant";
}
else
{
DataEntryType = firstToken.wordToken();
}
}
else
{

View File

@ -36,6 +36,12 @@ Foam::Table<Type>::Table(const word& entryName, const dictionary& dict)
Istream& is(dict.lookup(entryName));
word entryType(is);
token firstToken(is);
is.putBack(firstToken);
if (firstToken == token::BEGIN_SQR)
{
is >> this->dimensions_;
}
is >> this->table_;
TableBase<Type>::check();

View File

@ -30,7 +30,7 @@ Description
in the form, e.g. for an entry \<entryName\> that is (scalar, vector):
\verbatim
<entryName> table
<entryName> table [0 1 0 0 0] //dimension set optional
(
0.0 (1 2 3)
1.0 (4 5 6)
@ -129,6 +129,22 @@ public:
return TableBase<Type>::integrate(x1, x2);
}
//- Return dimensioned constant value
virtual dimensioned<Type> dimValue(const scalar x) const
{
return TableBase<Type>::dimValue(x);
}
//- Integrate between two values and return dimensioned type
virtual dimensioned<Type> dimIntegrate
(
const scalar x1,
const scalar x2
)
{
return TableBase<Type>::dimIntegrate(x1, x2);
}
// I/O

View File

@ -39,7 +39,8 @@ Foam::TableBase<Type>::TableBase(const word& name, const dictionary& dict)
dict.lookupOrDefault<word>("outOfBounds", "clamp")
)
),
table_()
table_(),
dimensions_(dimless)
{}
@ -48,7 +49,8 @@ Foam::TableBase<Type>::TableBase(const TableBase<Type>& tbl)
:
name_(tbl.name_),
boundsHandling_(tbl.boundsHandling_),
table_(tbl.table_)
table_(tbl.table_),
dimensions_(tbl.dimensions_)
{}
@ -330,6 +332,11 @@ Type Foam::TableBase<Type>::value(const scalar x) const
i++;
}
Info <<
(xDash - table_[i].first())/(table_[i+1].first() - table_[i].first())
* (table_[i+1].second() - table_[i].second())
+ table_[i].second() << endl;
// Linear interpolation to find value
return Type
(
@ -403,6 +410,28 @@ Type Foam::TableBase<Type>::integrate(const scalar x1, const scalar x2) const
}
template<class Type>
Foam::dimensioned<Type> Foam::TableBase<Type>::
dimValue(const scalar x) const
{
return dimensioned<Type>("dimensionedValue", dimensions_, this->value(x));
}
template<class Type>
Foam::dimensioned<Type> Foam::TableBase<Type>::dimIntegrate
(
const scalar x1, const scalar x2
) const
{
return dimensioned<Type>
(
"dimensionedValue",
dimensions_,
this->integrate(x2, x1)
);
}
// * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * * * //
#include "TableBaseIO.C"

View File

@ -37,6 +37,7 @@ SourceFiles
#include "DataEntry.H"
#include "Tuple2.H"
#include "dimensionSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -87,6 +88,9 @@ protected:
//- Table data
List<Tuple2<scalar, Type> > table_;
//- The dimension set
dimensionSet dimensions_;
// Protected Member Functions
@ -138,6 +142,16 @@ public:
//- Integrate between two (scalar) values
virtual Type integrate(const scalar x1, const scalar x2) const;
//- Return dimensioned constant value
virtual dimensioned<Type> dimValue(const scalar x) const;
//- Integrate between two values and return dimensioned type
virtual dimensioned<Type> dimIntegrate
(
const scalar x1,
const scalar x2
) const;
// I/O

View File

@ -37,6 +37,11 @@ Foam::TableFile<Type>::TableFile(const word& entryName, const dictionary& dict)
const dictionary coeffs(dict.subDict(type() + "Coeffs"));
coeffs.lookup("fileName") >> fName_;
if (coeffs.found("dimensions"))
{
coeffs.lookup("dimensions") >> this->dimensions_;
}
fileName expandedFile(fName_);
IFstream is(expandedFile.expand());

View File

@ -31,8 +31,9 @@ Description
<entryName> tableFile;
tableFileCoeffs
{
fileName dataFile; // name of data file
outOfBounds clamp; // optional out-of-bounds handling
dimensions [0 0 1 0 0]; // optional dimensions
fileName dataFile; // name of data file
outOfBounds clamp; // optional out-of-bounds handling
}
\endverbatim
@ -145,6 +146,22 @@ public:
return TableBase<Type>::integrate(x1, x2);
}
//- Return dimensioned constant value
virtual dimensioned<Type> dimValue(const scalar x) const
{
return TableBase<Type>::dimValue(x);
}
//- Integrate between two values and return dimensioned type
virtual dimensioned<Type> dimIntegrate
(
const scalar x1,
const scalar x2
)
{
return TableBase<Type>::dimIntegrate(x1, x2);
}
// I/O

View File

@ -42,11 +42,19 @@ Foam::polynomial::polynomial(const word& entryName, const dictionary& dict)
:
DataEntry<scalar>(entryName),
coeffs_(),
canIntegrate_(true)
canIntegrate_(true),
dimensions_(dimless)
{
Istream& is(dict.lookup(entryName));
word entryType(is);
token firstToken(is);
is.putBack(firstToken);
if (firstToken == token::BEGIN_SQR)
{
is >> this->dimensions_;
}
is >> coeffs_;
if (!coeffs_.size())
@ -85,7 +93,8 @@ Foam::polynomial::polynomial
:
DataEntry<scalar>(entryName),
coeffs_(coeffs),
canIntegrate_(true)
canIntegrate_(true),
dimensions_(dimless)
{
if (!coeffs_.size())
{
@ -125,7 +134,8 @@ Foam::polynomial::polynomial(const polynomial& poly)
:
DataEntry<scalar>(poly),
coeffs_(poly.coeffs_),
canIntegrate_(poly.canIntegrate_)
canIntegrate_(poly.canIntegrate_),
dimensions_(poly.dimensions_)
{}
@ -180,4 +190,26 @@ Foam::scalar Foam::polynomial::integrate(const scalar x1, const scalar x2) const
}
Foam::dimensioned<Foam::scalar> Foam::polynomial::dimValue
(
const scalar x
) const
{
return dimensioned<scalar>("dimensionedValue", dimensions_, value(x));
}
Foam::dimensioned<Foam::scalar> Foam::polynomial::dimIntegrate
(
const scalar x1, const scalar x2
) const
{
return dimensioned<scalar>
(
"dimensionedValue",
dimensions_,
integrate(x1, x2)
);
}
// ************************************************************************* //

View File

@ -30,7 +30,7 @@ Description
describes y = x^2 + 2x^3
\verbatim
<entryName> polynomial
<entryName> polynomial [0 0 1 0 0] // optional dimensions
(
(1 2)
(2 3)
@ -47,6 +47,7 @@ SourceFiles
#include "DataEntry.H"
#include "Tuple2.H"
#include "dimensionSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -79,6 +80,9 @@ class polynomial
//- Flag to indicate whether poly can be integrated
bool canIntegrate_;
//- The dimension set
dimensionSet dimensions_;
// Private Member Functions
@ -129,6 +133,16 @@ public:
//- Integrate between two (scalar) values
scalar integrate(const scalar x1, const scalar x2) const;
//- Return dimensioned constant value
dimensioned<scalar> dimValue(const scalar) const;
//- Integrate between two values and return dimensioned type
dimensioned<scalar> dimIntegrate
(
const scalar x1,
const scalar x2
) const;
// I/O

View File

@ -240,10 +240,15 @@ Foam::dynamicRefineFvMesh::refine
}
}
// // Remove the stored tet base points
// tetBasePtIsPtr_.clear();
// // Remove the cell tree
// cellTreePtr_.clear();
// Update fields
updateMesh(map);
// Move mesh
/*
pointField newPoints;

View File

@ -15,6 +15,10 @@ basicSource/rotorDiskSource/profileModel/profileModel.C
basicSource/rotorDiskSource/profileModel/profileModelList.C
basicSource/rotorDiskSource/profileModel/lookup/lookupProfile.C
basicSource/rotorDiskSource/profileModel/series/seriesProfile.C
basicSource/rotorDiskSource/trimModel/trimModel/trimModel.C
basicSource/rotorDiskSource/trimModel/trimModel/trimModelNew.C
basicSource/rotorDiskSource/trimModel/fixed/fixedTrim.C
basicSource/rotorDiskSource/trimModel/targetForce/targetForceTrim.C
basicSource/actuationDiskSource/actuationDiskSource.C
basicSource/radialActuationDiskSource/radialActuationDiskSource.C
@ -25,4 +29,4 @@ $(interRegion)/constantHeatTransfer/constantHeatTransfer.C
$(interRegion)/tabulatedHeatTransfer/tabulatedHeatTransfer.C
$(interRegion)/variableHeatTransfer/variableHeatTransfer.C
LIB = $(FOAM_LIBBIN)/libfieldSources
LIB = $(FOAM_LIBBIN)/libfieldSources

View File

@ -26,8 +26,8 @@ License
#include "rotorDiskSource.H"
#include "addToRunTimeSelectionTable.H"
#include "mathematicalConstants.H"
#include "trimModel.H"
#include "unitConversion.H"
#include "geometricOneField.H"
#include "fvMatrices.H"
#include "syncTools.H"
@ -123,6 +123,8 @@ void Foam::rotorDiskSource::checkData()
void Foam::rotorDiskSource::setFaceArea(vector& axis, const bool correct)
{
area_ = 0.0;
static const scalar tol = 0.8;
const label nInternalFaces = mesh_.nInternalFaces();
@ -299,12 +301,15 @@ void Foam::rotorDiskSource::createCoordinateSystem()
reduce(axis, maxMagSqrOp<vector>());
axis /= mag(axis);
// axis direction is somewhat arbitrary - check if user needs
// needs to reverse
bool reverse(readBool(coeffs_.lookup("reverseAxis")));
if (reverse)
// correct the axis direction using a point above the rotor
{
axis *= -1.0;
vector pointAbove(coeffs_.lookup("pointAbove"));
vector dir = pointAbove - origin;
dir /= mag(dir);
if ((dir & axis) < 0)
{
axis *= -1.0;
}
}
coeffs_.lookup("refDirection") >> refDir;
@ -380,9 +385,6 @@ void Foam::rotorDiskSource::constructGeometry()
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 [radians]
alphag_[i] = trim_.alphaC + trim_.A*cos(psi) + trim_.B*sin(psi);
}
}
@ -444,13 +446,12 @@ Foam::rotorDiskSource::rotorDiskSource
inletVelocity_(vector::zero),
tipEffect_(1.0),
flap_(),
trim_(),
trim_(trimModel::New(*this, coeffs_)),
blade_(coeffs_.subDict("blade")),
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),
rMax_(0.0)
@ -467,34 +468,179 @@ Foam::rotorDiskSource::~rotorDiskSource()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::rotorDiskSource::calculate
(
const vectorField& U,
const scalarField& alphag,
vectorField& force,
const bool divideVolume,
const bool output
) const
{
tmp<volScalarField> trho;
if (rhoName_ != "none")
{
trho = mesh_.lookupObject<volScalarField>(rhoName_);
}
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 = coordSys_.localVector(U[cellI]);
// apply correction in local system due to coning
Uc = R_[i] & Uc;
// set radial component of velocity to zero
Uc.x() = 0.0;
// remove blade linear velocity from blade normal component
Uc.y() -= radius*omega_;
// determine blade data for this radius
// i1 = index of upper 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);
// effective angle of attack
scalar alphaEff =
mathematical::pi + atan2(Uc.z(), Uc.y()) - (alphag[i] + twist);
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*magSqr(Uc);
if (trho.valid())
{
pDyn *= trho()[cellI];
}
scalar f = pDyn*chord*nBlades_*area_[i]/mathematical::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();
liftEff += localForce.z();
// convert force to global cartesian co-ordinate system
force[cellI] = coordSys_.globalVector(localForce);
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::rotorDiskSource::addSup(fvMatrix<vector>& eqn, const label fieldI)
{
// add source to rhs of eqn
const volVectorField& U = eqn.psi();
dimensionSet dims = dimless;
if (eqn.dimensions() == dimForce)
{
coeffs_.lookup("rhoName") >> rhoName_;
const volScalarField& rho =
mesh_.lookupObject<volScalarField>(rhoName_);
eqn -= calculateForces
(
rho.internalField(),
inflowVelocity(U),
dimForce/dimVolume
);
dims.reset(dimForce/dimVolume);
}
else
{
eqn -= calculateForces
(
oneField(),
inflowVelocity(U),
dimForce/dimVolume/dimDensity
);
dims.reset(dimForce/dimVolume/dimDensity);
}
volVectorField force
(
IOobject
(
"rotorForce",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedVector("zero", dims, vector::zero)
);
const volVectorField& U = eqn.psi();
const vectorField Uin(inflowVelocity(U));
trim_->correct(Uin, force);
calculate(Uin, trim_->thetag(), force);
// add source to rhs of eqn
eqn -= force;
if (mesh_.time().outputTime())
{
force.write();
}
}
@ -527,12 +673,10 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
flapCoeffs.lookup("beta1") >> flap_.beta1;
flapCoeffs.lookup("beta2") >> flap_.beta2;
flap_.beta0 = degToRad(flap_.beta0);
flap_.beta1 = degToRad(flap_.beta1);
flap_.beta2 = degToRad(flap_.beta2);
const dictionary& trimCoeffs(coeffs_.subDict("trimCoeffs"));
trimCoeffs.lookup("alphaC") >> trim_.alphaC;
trimCoeffs.lookup("A") >> trim_.A;
trimCoeffs.lookup("B") >> trim_.B;
trim_.alphaC = degToRad(trim_.alphaC);
trim_->read(coeffs_);
checkData();
@ -542,7 +686,7 @@ bool Foam::rotorDiskSource::read(const dictionary& dict)
if (debug)
{
writeField("alphag", alphag_, true);
writeField("alphag", trim_->thetag()(), true);
writeField("faceArea", area_, true);
}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,13 +37,16 @@ Description
fieldNames (U); // names of fields on which to apply source
rhoName rho; // density field if compressible case
nBlades 3; // number of blades
tip effect 0.96; // normalised radius above which lift = 0
tipEffect 0.96; // normalised radius above which lift = 0
inletFlowType local; // inlet flow type specification
geometryMode auto; // geometry specification
refDirection (-1 0 0); // reference direction
// - used as reference for psi angle
trimModel fixed; // fixed || targetForce
flapCoeffs
{
@ -51,16 +54,12 @@ Description
beta1 0; // lateral flapping coeff
beta2 0; // longitudinal flapping coeff
}
trimCoeffs
{
alphac 15; // collective pitch angle [deg]
A 0; // lateral cyclic coeff
B 0; // longitudinal cyclic coeff
}
blade
{
...
}
profiles
{
...
@ -102,6 +101,9 @@ SourceFiles
namespace Foam
{
// Forward declaration of classes
class trimModel;
/*---------------------------------------------------------------------------*\
Class rotorDiskSource Declaration
\*---------------------------------------------------------------------------*/
@ -140,13 +142,6 @@ protected:
scalar beta2; // longitudinal flapping coeff
};
struct trimData
{
scalar alphaC; // collective pitch angle
scalar A; // lateral cyclic coeff
scalar B; // longitudinal cyclic coeff
};
// Protected data
@ -172,8 +167,8 @@ protected:
//- Blade flap coefficients [rad/s]
flapData flap_;
//- Blad trim coefficients
trimData trim_;
//- Trim model
autoPtr<trimModel> trim_;
//- Blade data
bladeModel blade_;
@ -181,7 +176,8 @@ protected:
//- Profile data
profileModelList profiles_;
//- Cell centre positions in local rotor frame (Cartesian x, y, z)
//- Cell centre positions in local rotor frame
// (Cylindrical r, theta, z)
List<point> x_;
//- Rotation tensor for flap angle
@ -190,9 +186,6 @@ protected:
//- Inverse rotation tensor for flap angle
List<tensor> invR_;
//- Geometric angle of attack [deg]
List<scalar> alphag_;
//- Area [m2]
List<scalar> area_;
@ -220,15 +213,6 @@ protected:
//- Return the inlet flow field
tmp<vectorField> inflowVelocity(const volVectorField& U) const;
//- Calculate forces
template<class RhoType>
tmp<volVectorField> calculateForces
(
const RhoType& rho,
const vectorField& U,
const dimensionSet& dims
);
//- Helper function to write rotor values
template<class Type>
void writeField
@ -264,6 +248,29 @@ public:
// Member Functions
// Access
//- Return the cell centre positions in local rotor frame
// (Cylindrical r, theta, z)
inline const List<point>& x() const;
//- Return the rotor co-ordinate system (r, theta, z)
inline const cylindricalCS& coordSys() const;
// Evaluation
//- Calculate forces
void calculate
(
const vectorField& U,
const scalarField& alphag,
vectorField& force,
const bool divideVolume = true,
const bool output = true
) const;
// Source term addition
//- Source term to fvMatrix<vector>
@ -286,6 +293,10 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "rotorDiskSourceI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "rotorDiskSourceTemplates.C"
#endif

View File

@ -0,0 +1,41 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "rotorDiskSource.H"
const Foam::List<Foam::point>& Foam::rotorDiskSource::x() const
{
return x_;
}
const Foam::cylindricalCS& Foam::rotorDiskSource::coordSys() const
{
return coordSys_;
}
// ************************************************************************* //

View File

@ -24,13 +24,8 @@ License
\*---------------------------------------------------------------------------*/
#include "rotorDiskSource.H"
#include "addToRunTimeSelectionTable.H"
#include "mathematicalConstants.H"
#include "unitConversion.H"
#include "volFields.H"
using namespace Foam::constant::mathematical;
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
template<class Type>
@ -81,149 +76,4 @@ void Foam::rotorDiskSource::writeField
}
template<class RhoType>
Foam::tmp<Foam::volVectorField> Foam::rotorDiskSource::calculateForces
(
const RhoType& rho,
const vectorField& U,
const dimensionSet& dims
)
{
tmp<volVectorField> tForce
(
new volVectorField
(
IOobject
(
"rotorForce",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedVector("zero", dims, vector::zero)
)
);
vectorField& force = tForce().internalField();
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 = coordSys_.localVector(U[cellI]);
// apply correction in local system due to coning
Uc = R_[i] & Uc;
// set radial component of velocity to zero
Uc.x() = 0.0;
// remove blade linear velocity from blade normal component
Uc.y() -= radius*omega_;
// velocity magnitude
scalar magUc = mag(Uc);
// determine blade data for this radius
// i1 = index of upper 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);
// effective angle of attack
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);
// 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 = 1.0;
if (radius/rMax_ > tipEffect_)
{
tipFactor = 0.0;
}
// 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();
liftEff += localForce.z();
// convert force to global cartesian co-ordinate system
force[cellI] = coordSys_.globalVector(localForce);
force[cellI] /= V[cellI];
}
}
if (mesh_.time().outputTime())
{
tForce().write();
}
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;
return tForce;
}
// ************************************************************************* //

View File

@ -0,0 +1,96 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "fixedTrim.H"
#include "addToRunTimeSelectionTable.H"
#include "unitConversion.H"
#include "mathematicalConstants.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(fixedTrim, 0);
addToRunTimeSelectionTable(trimModel, fixedTrim, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::fixedTrim::fixedTrim(const rotorDiskSource& rotor, const dictionary& dict)
:
trimModel(rotor, dict, typeName),
thetag_(rotor.cells().size(), 0.0)
{
read(dict);
}
// * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * //
Foam::fixedTrim::~fixedTrim()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::fixedTrim::read(const dictionary& dict)
{
trimModel::read(dict);
scalar theta0 = degToRad(readScalar(coeffs_.lookup("theta0")));
scalar theta1c = degToRad(readScalar(coeffs_.lookup("theta1c")));
scalar theta1s = degToRad(readScalar(coeffs_.lookup("theta1s")));
const List<vector>& x = rotor_.x();
forAll(thetag_, i)
{
scalar psi = x[i].y();
if (psi < 0)
{
psi += mathematical::twoPi;
}
thetag_[i] = theta0 + theta1c*cos(psi) + theta1s*sin(psi);
}
}
Foam::tmp<Foam::scalarField> Foam::fixedTrim::thetag() const
{
return tmp<scalarField>(thetag_);
}
void Foam::fixedTrim::correct(const vectorField& U, vectorField& force)
{
// do nothing
}
// ************************************************************************* //

View File

@ -0,0 +1,95 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::fixedTrim
Description
Fixed trim coefficients
SourceFiles
fixedTrim.C
\*---------------------------------------------------------------------------*/
#ifndef fixedTrim_H
#define fixedTrim_H
#include "trimModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class fixedTrim Declaration
\*---------------------------------------------------------------------------*/
class fixedTrim
:
public trimModel
{
protected:
// Protected data
//- Geometric angle of attack [rad]
scalarField thetag_;
public:
//- Run-time type information
TypeName("fixedTrim");
//- Constructor
fixedTrim(const rotorDiskSource& rotor, const dictionary& dict);
//- Destructor
virtual ~fixedTrim();
// Member functions
//- Read
void read(const dictionary& dict);
//- Return the geometric angle of attack [rad]
virtual tmp<scalarField> thetag() const;
//- Correct the model
virtual void correct(const vectorField& U, vectorField& force);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,234 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "targetForceTrim.H"
#include "addToRunTimeSelectionTable.H"
#include "unitConversion.H"
#include "mathematicalConstants.H"
using namespace Foam::constant;
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(targetForceTrim, 0);
addToRunTimeSelectionTable(trimModel, targetForceTrim, dictionary);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
Foam::vector Foam::targetForceTrim::calcForce
(
const vectorField& U,
const scalarField& thetag,
vectorField& force
) const
{
rotor_.calculate(U, thetag, force, false, false);
const labelList& cells = rotor_.cells();
const vectorField& C = rotor_.mesh().C();
const vector& origin = rotor_.coordSys().origin();
const vector& rollAxis = rotor_.coordSys().e1();
const vector& pitchAxis = rotor_.coordSys().e2();
const vector& yawAxis = rotor_.coordSys().e3();
vector f(vector::zero);
forAll(cells, i)
{
label cellI = cells[i];
vector moment = force[cellI]^(C[cellI] - origin);
f[0] += force[cellI] & yawAxis;
f[1] += moment & pitchAxis;
f[2] += moment & rollAxis;
}
reduce(f, sumOp<vector>());
return f;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::targetForceTrim::targetForceTrim
(
const rotorDiskSource& rotor,
const dictionary& dict
)
:
trimModel(rotor, dict, typeName),
calcFrequency_(-1),
target_(vector::zero),
theta_(vector::zero),
nIter_(50),
tol_(1e-8),
relax_(1.0),
dTheta_(degToRad(0.1))
{
read(dict);
}
// * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * //
Foam::targetForceTrim::~targetForceTrim()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::targetForceTrim::read(const dictionary& dict)
{
trimModel::read(dict);
const dictionary& targetDict(coeffs_.subDict("target"));
target_[0] = readScalar(targetDict.lookup("fThrust"));
target_[1] = readScalar(targetDict.lookup("mPitch"));
target_[2] = readScalar(targetDict.lookup("mRoll"));
const dictionary& pitchAngleDict(coeffs_.subDict("pitchAngles"));
theta_[0] = degToRad(readScalar(pitchAngleDict.lookup("theta0Ini")));
theta_[1] = degToRad(readScalar(pitchAngleDict.lookup("theta1cIni")));
theta_[2] = degToRad(readScalar(pitchAngleDict.lookup("theta1sIni")));
coeffs_.lookup("calcFrequency") >> calcFrequency_;
coeffs_.readIfPresent("nIter", nIter_);
coeffs_.readIfPresent("tol", tol_);
coeffs_.readIfPresent("relax", relax_);
if (coeffs_.readIfPresent("dTheta", dTheta_))
{
dTheta_ = degToRad(dTheta_);
}
}
Foam::tmp<Foam::scalarField> Foam::targetForceTrim::thetag() const
{
const List<vector>& x = rotor_.x();
tmp<scalarField> ttheta(new scalarField(x.size()));
scalarField& t = ttheta();
forAll(t, i)
{
scalar psi = x[i].y();
if (psi < 0)
{
psi += mathematical::twoPi;
}
t[i] = theta_[0] + theta_[1]*cos(psi) + theta_[2]*sin(psi);
}
return ttheta;
}
void Foam::targetForceTrim::correct(const vectorField& U, vectorField& force)
{
if (rotor_.mesh().time().timeIndex() % calcFrequency_ == 0)
{
// iterate to find new pitch angles to achieve target force
scalar err = GREAT;
label iter = 0;
tensor J(tensor::zero);
while ((err > tol_) && (iter < nIter_))
{
// cache initial theta vector
vector theta0(theta_);
// set initial values
vector old = calcForce(U, thetag(), force);
// construct Jacobian by perturbing the pitch angles
// by +/-(dTheta_/2)
for (label pitchI = 0; pitchI < 3; pitchI++)
{
theta_[pitchI] -= dTheta_/2.0;
vector f0 = calcForce(U, thetag(), force);
theta_[pitchI] += dTheta_;
vector f1 = calcForce(U, thetag(), force);
vector ddTheta = (f1 - f0)/dTheta_;
J[pitchI + 0] = ddTheta[0];
J[pitchI + 3] = ddTheta[1];
J[pitchI + 6] = ddTheta[2];
theta_ = theta0;
}
// calculate the change in pitch angle vector
vector dt = inv(J) & (target_ - old);
// update pitch angles
vector thetaNew = theta_ + relax_*dt;
// update error
err = mag(thetaNew - theta_);
// update for next iteration
theta_ = thetaNew;
iter++;
}
if (iter == nIter_)
{
WarningIn
(
"void Foam::targetForceTrim::correct"
"("
"const vectorField&, "
"vectorField&"
")"
) << "Trim routine not converged in " << iter
<< " iterations, max residual = " << err << endl;
}
else
{
Info<< type() << ": converged in " << iter
<< " iterations" << endl;
}
Info<< " new pitch angles:" << nl
<< " theta0 = " << radToDeg(theta_[0]) << nl
<< " theta1c = " << radToDeg(theta_[1]) << nl
<< " theta1s = " << radToDeg(theta_[2]) << nl
<< endl;
}
}
// ************************************************************************* //

View File

@ -0,0 +1,126 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::targetForceTrim
Description
Target force trim coefficients
SourceFiles
targetForceTrim.C
\*---------------------------------------------------------------------------*/
#ifndef targetForceTrim_H
#define targetForceTrim_H
#include "trimModel.H"
#include "tensor.H"
#include "vector.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class targetForceTrim Declaration
\*---------------------------------------------------------------------------*/
class targetForceTrim
:
public trimModel
{
protected:
// Protected data
//- Number of iterations between calls to 'correct'
label calcFrequency_;
//- Target force [N]
vector target_;
//- Pitch angles (collective, roll, pitch) [rad]
vector theta_;
//- Maximum number of iterations in trim routine
label nIter_;
//- Convergence tolerance
scalar tol_;
//- Under-relaxation coefficient
scalar relax_;
//- Perturbation angle used to determine jacobian
scalar dTheta_;
// Protected member functions
//- Calculate the rotor forces
vector calcForce
(
const vectorField& U,
const scalarField& alphag,
vectorField& force
) const;
public:
//- Run-time type information
TypeName("targetForceTrim");
//- Constructor
targetForceTrim(const rotorDiskSource& rotor, const dictionary& dict);
//- Destructor
virtual ~targetForceTrim();
// Member functions
//- Read
void read(const dictionary& dict);
//- Return the geometric angle of attack [rad]
virtual tmp<scalarField> thetag() const;
//- Correct the model
virtual void correct(const vectorField& U, vectorField& force);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,67 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "trimModel.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(trimModel, 0);
defineRunTimeSelectionTable(trimModel, dictionary);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::trimModel::trimModel
(
const rotorDiskSource& rotor,
const dictionary& dict,
const word& name
)
:
rotor_(rotor),
name_(name),
coeffs_(dictionary::null)
{
read(dict);
}
// * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * //
Foam::trimModel::~trimModel()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::trimModel::read(const dictionary& dict)
{
coeffs_ = dict.subDict(name_ + "Coeffs");
}
// ************************************************************************* //

View File

@ -0,0 +1,135 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::trimModel
Description
Trim model base class
SourceFiles
trimModel.C
\*---------------------------------------------------------------------------*/
#ifndef trimModel_H
#define trimModel_H
#include "rotorDiskSource.H"
#include "dictionary.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class trimModel Declaration
\*---------------------------------------------------------------------------*/
class trimModel
{
protected:
// Protected data
//- Reference to the rotor source model
const rotorDiskSource& rotor_;
//- Name of model
const word name_;
//- Coefficients dictionary
dictionary coeffs_;
public:
//- Run-time type information
TypeName("trimModel");
// Declare runtime constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
trimModel,
dictionary,
(
const rotorDiskSource& rotor,
const dictionary& dict
),
(rotor, dict)
);
// Constructors
//- Construct from components
trimModel
(
const rotorDiskSource& rotor,
const dictionary& dict,
const word& name
);
// Selectors
//- Return a reference to the selected trim model
static autoPtr<trimModel> New
(
const rotorDiskSource& rotor,
const dictionary& dict
);
//- Destructor
virtual ~trimModel();
// Member functions
//- Read
virtual void read(const dictionary& dict);
//- Return the geometric angle of attack [rad]
virtual tmp<scalarField> thetag() const = 0;
//- Correct the model
virtual void correct(const vectorField& U, vectorField& force) = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,59 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "trimModel.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::autoPtr<Foam::trimModel> Foam::trimModel::New
(
const rotorDiskSource& rotor,
const dictionary& dict
)
{
const word modelType(dict.lookup(typeName));
Info<< " Selecting " << typeName << " " << modelType << endl;
dictionaryConstructorTable::iterator cstrIter =
dictionaryConstructorTablePtr_->find(modelType);
if (cstrIter == dictionaryConstructorTablePtr_->end())
{
FatalErrorIn
(
"trimModel::New(const rotorDiskSource&, const dictionary&)"
) << "Unknown " << typeName << " type "
<< modelType << nl << nl
<< "Valid " << typeName << " types are:" << nl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<trimModel>(cstrIter()(rotor, dict));
}
// ************************************************************************* //

View File

@ -240,8 +240,7 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
),
origin_(dict_.lookup("origin")),
axis_(dict_.lookup("axis")),
omega_(dict_.lookup("omega")),
Omega_("Omega", omega_*axis_)
omega_(DataEntry<scalar>::New("omega", dict_))
{
if (dict_.found("patches"))
{
@ -256,7 +255,6 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
axis_ = axis_/mag(axis_);
Omega_ = omega_*axis_;
excludedPatchLabels_.setSize(excludedPatchNames_.size());
@ -309,7 +307,9 @@ void Foam::MRFZone::addCoriolis
vectorField& ddtUc = ddtU.internalField();
const vectorField& Uc = U.internalField();
const vector& Omega = Omega_.value();
const scalar t = mesh_.time().timeOutputValue();
const vector& Omega = omega_->value(t)*axis_;
forAll(cells, i)
{
@ -331,7 +331,9 @@ void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const
vectorField& Usource = UEqn.source();
const vectorField& U = UEqn.psi();
const vector& Omega = Omega_.value();
const scalar t = mesh_.time().timeOutputValue();
const vector& Omega = omega_->value(t)*axis_;
forAll(cells, i)
{
@ -357,7 +359,9 @@ void Foam::MRFZone::addCoriolis
vectorField& Usource = UEqn.source();
const vectorField& U = UEqn.psi();
const vector& Omega = Omega_.value();
const scalar t = mesh_.time().timeOutputValue();
const vector& Omega = omega_->value(t)*axis_;
forAll(cells, i)
{
@ -371,9 +375,11 @@ void Foam::MRFZone::relativeVelocity(volVectorField& U) const
{
const volVectorField& C = mesh_.C();
const vector& origin = origin_.value();
const vector& origin = origin_;
const vector& Omega = Omega_.value();
const scalar t = mesh_.time().timeOutputValue();
const vector& Omega = omega_->value(t)*axis_;
const labelList& cells = mesh_.cellZones()[cellZoneID_];
@ -411,9 +417,11 @@ void Foam::MRFZone::absoluteVelocity(volVectorField& U) const
{
const volVectorField& C = mesh_.C();
const vector& origin = origin_.value();
const vector& origin = origin_;
const vector& Omega = Omega_.value();
const scalar t = mesh_.time().timeOutputValue();
const vector& Omega = omega_->value(t)*axis_;
const labelList& cells = mesh_.cellZones()[cellZoneID_];
@ -481,9 +489,12 @@ void Foam::MRFZone::absoluteFlux
void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
{
const vector& origin = origin_.value();
const vector& origin = origin_;
const scalar t = mesh_.time().timeOutputValue();
const vector& Omega = omega_->value(t)*axis_;
const vector& Omega = Omega_.value();
// Included patches
forAll(includedFaces_, patchi)
@ -511,7 +522,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const MRFZone& MRF)
os << token::BEGIN_BLOCK << incrIndent << nl;
os.writeKeyword("origin") << MRF.origin_ << token::END_STATEMENT << nl;
os.writeKeyword("axis") << MRF.axis_ << token::END_STATEMENT << nl;
os.writeKeyword("omega") << MRF.omega_ << token::END_STATEMENT << nl;
MRF.omega_->writeData(os);
if (MRF.excludedPatchNames_.size())
{

View File

@ -49,6 +49,7 @@ SourceFiles
#include "fvMatricesFwd.H"
#include "fvMatrices.H"
#include "DataEntry.H"
#include "autoPtr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -86,10 +87,14 @@ class MRFZone
//- Excluded faces (per patch) that do not move with the MRF
labelListList excludedFaces_;
const dimensionedVector origin_;
dimensionedVector axis_;
dimensionedScalar omega_;
dimensionedVector Omega_;
//- Origin of the axis
const vector origin_;
//- Axis vector
vector axis_;
//- Angular velocty (rad/sec)
autoPtr<DataEntry<scalar> > omega_;
// Private Member Functions

View File

@ -41,8 +41,10 @@ void Foam::MRFZone::relativeRhoFlux
const surfaceVectorField& Cf = mesh_.Cf();
const surfaceVectorField& Sf = mesh_.Sf();
const vector& origin = origin_.value();
const vector& Omega = Omega_.value();
const vector& origin = origin_;
const scalar t = mesh_.time().timeOutputValue();
const vector& Omega = omega_->value(t)*axis_;
const vectorField& Cfi = Cf.internalField();
const vectorField& Sfi = Sf.internalField();
@ -92,8 +94,10 @@ void Foam::MRFZone::absoluteRhoFlux
const surfaceVectorField& Cf = mesh_.Cf();
const surfaceVectorField& Sf = mesh_.Sf();
const vector& origin = origin_.value();
const vector& Omega = Omega_.value();
const vector& origin = origin_;
const scalar t = mesh_.time().timeOutputValue();
const vector& Omega = omega_->value(t)*axis_;
const vectorField& Cfi = Cf.internalField();
const vectorField& Sfi = Sf.internalField();

View File

@ -431,7 +431,7 @@ void Foam::mappedPatchBase::calcMapping() const
tmp<pointField> patchPoints(facePoints(patch_));
// Get offsetted points
const pointField offsettedPoints = samplePoints(patchPoints());
const pointField offsettedPoints(samplePoints(patchPoints()));
// Do a sanity check

View File

@ -153,14 +153,14 @@ Type Foam::fieldValues::faceSource::processSameTypeValues
}
case opAreaAverage:
{
const scalarField magSf = mag(Sf);
const scalarField magSf(mag(Sf));
result = sum(values*magSf)/sum(magSf);
break;
}
case opAreaIntegrate:
{
const scalarField magSf = mag(Sf);
const scalarField magSf(mag(Sf));
result = sum(values*magSf);
break;
@ -177,7 +177,7 @@ Type Foam::fieldValues::faceSource::processSameTypeValues
}
case opCoV:
{
const scalarField magSf = mag(Sf);
const scalarField magSf(mag(Sf));
Type meanValue = sum(values*magSf)/sum(magSf);

View File

@ -138,7 +138,7 @@ Foam::labelList Foam::springRenumber::renumber
<< " average force:" << average(mag(sumForce)) << endl;
// Determine displacement.
scalarField displacement = deltaT*sumForce;
scalarField displacement(deltaT*sumForce);
//Pout<< "Displacement :" << nl
// << " min : " << min(displacement) << nl

View File

@ -157,7 +157,7 @@ void Foam::ensightSetWriter<Type>::write
<< "coordinates" << nl;
for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
{
const scalarField fld = valueSets[setI]->component(cmpt);
const scalarField fld(valueSets[setI]->component(cmpt));
forAll(fld, i)
{
if (mag(fld[i]) >= scalar(floatScalarVSMALL))
@ -293,7 +293,7 @@ void Foam::ensightSetWriter<Type>::write
cmpt++
)
{
const scalarField fld = fieldVals[trackI].component(cmpt);
const scalarField fld(fieldVals[trackI].component(cmpt));
forAll(fld, i)
{
if (mag(fld[i]) >= scalar(floatScalarVSMALL))

View File

@ -317,7 +317,8 @@ void LaunderSharmaKE::correct()
- fvm::Sp(fvc::ddt(rho_) + fvc::div(phi_), epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
==
C1_*G*epsilon_/k_ + fvm::SuSp((C3_ - 2.0/3.0*C1_)*rho_*divU, epsilon_)
C1_*G*epsilon_/k_
- fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)
- fvm::Sp(C2_*f2()*rho_*epsilon_/k_, epsilon_)
//+ 0.75*1.5*flameKproduction*epsilon_/k_
+ E

View File

@ -22,9 +22,9 @@ FoamFile
// Fixed patches (by default they 'move' with the MRF zone)
nonRotatingPatches ();
origin origin [0 1 0 0 0 0 0] (0 0 0);
axis axis [0 0 0 0 0 0 0] (0 0 1);
omega omega [0 0 -1 0 0 0 0] 1047.2;
origin (0 0 0);
axis (0 0 1);
omega constant 1047.2;
}
)

View File

@ -22,9 +22,9 @@ FoamFile
// Fixed patches (by default they 'move' with the MRF zone)
nonRotatingPatches ();
origin origin [0 1 0 0 0 0 0] (0 0 0);
axis axis [0 0 0 0 0 0 0] (0 0 1);
omega omega [0 0 -1 0 0 0 0] 104.72;
origin (0 0 0);
axis (0 0 1);
omega 104.72;
}
)

View File

@ -22,9 +22,9 @@ FoamFile
// Fixed patches (by default they 'move' with the MRF zone)
nonRotatingPatches ();
origin origin [0 1 0 0 0 0 0] (0 0 0);
axis axis [0 0 0 0 0 0 0] (0 0 1);
omega omega [0 0 -1 0 0 0 0] 104.72;
origin (0 0 0);
axis (0 0 1);
omega constant 104.72;
}
)

View File

@ -22,9 +22,9 @@ FoamFile
// Fixed patches (by default they 'move' with the MRF zone)
nonRotatingPatches ();
origin origin [0 1 0 0 0 0 0] (0 0 0);
axis axis [0 0 0 0 0 0 0] (0 0 1);
omega omega [0 0 -1 0 0 0 0] 6.2831853;
origin (0 0 0);
axis (0 0 1);
omega constant 6.2831853;
}
)

View File

@ -22,9 +22,9 @@ FoamFile
// Fixed patches (by default they 'move' with the MRF zone)
nonRotatingPatches ();
origin origin [0 1 0 0 0 0 0] (0 0 0);
axis axis [0 0 0 0 0 0 0] (0 0 1);
omega omega [0 0 -1 0 0 0 0] 6.2831853;
origin (0 0 0);
axis (0 0 1);
omega constant 6.2831853;
}
)

View File

@ -22,9 +22,9 @@ FoamFile
// Fixed patches (by default they 'move' with the MRF zone)
nonRotatingPatches ();
origin origin [0 1 0 0 0 0 0] (0 0 0);
axis axis [0 0 0 0 0 0 0] (0 0 1);
omega omega [0 0 -1 0 0 0 0] 10.472;
origin (0 0 0);
axis (0 0 1);
omega constant 10.472;
}
)

View File

@ -22,9 +22,9 @@ FoamFile
// Fixed patches (by default they 'move' with the MRF zone)
nonRotatingPatches ();
origin origin [0 1 0 0 0 0 0] (0 0 0);
axis axis [0 0 0 0 0 0 0] (0 0 1);
omega omega [0 0 -1 0 0 0 0] 10.472;
origin (0 0 0);
axis (0 0 1);
omega constant 10.472;
}
)

View File

@ -22,9 +22,9 @@ FoamFile
// Fixed patches (by default they 'move' with the MRF zone)
nonRotatingPatches ();
origin origin [0 1 0 0 0 0 0] (0 0 0);
axis axis [0 0 0 0 0 0 0] (0 0 1);
omega omega [0 0 -1 0 0 0 0] 10.472;
origin (0 0 0);
axis (0 0 1);
omega constant 10.472;
}
)