54 { rotationMode::MRF,
"MRF" },
77 if (
dict.readIfPresent(
"rpm",
n_))
80 dict.readIfPresent(
"n",
n_);
84 n_ =
dict.get<scalar>(
"n");
91 const auto* MRFZones =
92 mesh_.cfindObject<IOMRFZoneList>(
"MRFProperties");
97 <<
"Unable to find MRFProperties in the database. "
98 <<
"Is this an MRF case?"
101 const auto& mrf = MRFZones->MRFZoneList::getFromName(
MRFName_);
103 dict.readIfPresent(
"originOffset", origin);
104 origin += mrf.origin();
122 if (
dict.readIfPresent(
"alphaAxis", alphaAxis))
124 alphaAxis.normalise();
127 new coordSystem::cylindrical(origin, axis, alphaAxis)
140 switch (rotationMode_)
142 case rotationMode::SPECIFIED:
147 case rotationMode::MRF:
149 const auto* MRFZones =
150 mesh_.cfindObject<IOMRFZoneList>(
"MRFProperties");
155 <<
"Unable to find MRFProperties in the database. "
156 <<
"Is this an MRF case?"
160 const auto& mrf = MRFZones->MRFZoneList::getFromName(MRFName_);
183 if (writePropellerPerformance_ && !propellerPerformanceFilePtr_)
185 propellerPerformanceFilePtr_ =
186 newFileAtStartTime(
"propellerPerformance");
187 auto&
os = propellerPerformanceFilePtr_();
190 writeHeaderValue(
os,
"CofR", coordSysPtr_->origin());
191 writeHeaderValue(
os,
"Radius", radius_);
192 writeHeaderValue(
os,
"Axis", coordSysPtr_->e3());
196 writeCommented(
os,
"Time");
197 writeTabbed(
os,
"n");
198 writeTabbed(
os,
"URef");
199 writeTabbed(
os,
"J");
200 writeTabbed(
os,
"KT");
201 writeTabbed(
os,
"10*KQ");
202 writeTabbed(
os,
"eta0");
206 if (writeWakeFields_)
222 <<
"Unable to find velocity field " << UName_
223 <<
" . Available vector fields are: "
245 label nPoint = nRadius*nTheta;
254 const label nFace = nRadius*nTheta;
256 points.resize_nocopy(nPoint);
260 const scalar zCoord = 0;
262 for (
int radiusi = 0; radiusi <= nRadius; ++radiusi)
264 if (r1 < SMALL && radiusi == 0)
266 points[pointi++] = origin;
270 const scalar r = r1 + radiusi*(r2 - r1)/nRadius;
272 for (label i = 0; i < nTheta; ++i)
287 const List<label> ptIDs(
identity(nTheta));
291 label pointOffset0 = 0;
292 label radiusOffset = 0;
293 DynamicList<label> facePts(4);
294 for (
int radiusi = 0; radiusi < nRadius; ++radiusi)
296 if (r1 < SMALL && radiusi == 0)
302 for (label thetai = 1; thetai <= nTheta; ++thetai)
308 facePts.append(thetai);
309 facePts.append(1 + ptIDs.fcIndex(thetai - 1));
311 faces[facei++] = face(facePts);
316 for (label thetai = 0; thetai < nTheta; ++thetai)
320 label offset = pointOffset0 + (radiusi-radiusOffset)*nTheta;
323 facePts.append(offset + ptIDs.fcIndex(thetai - 1));
324 facePts.append(offset + ptIDs.fcIndex(thetai));
327 facePts.append(offset + nTheta + ptIDs.fcIndex(thetai));
328 facePts.append(offset + nTheta + ptIDs.fcIndex(thetai - 1));
330 faces[facei++] =
face(facePts);
344 const scalar r1 = sampleDiskDict.getScalar(
"r1");
345 const scalar r2 = sampleDiskDict.getScalar(
"r2");
347 nTheta_ = sampleDiskDict.getLabel(
"nTheta");
348 nRadial_ = sampleDiskDict.getLabel(
"nRadial");
350 setSampleDiskGeometry
364 if (sampleDiskDict.readIfPresent(
"surfaceWriter", writerType))
378 surfaceWriterPtr_->useTimeDir(
true);
382 sampleDiskDict.getOrDefault(
"errorOnPointNotFound",
false);
390 if (!writeWakeFields_)
399 const auto& meshCells = mesh_.cells();
400 const auto& meshFaces = mesh_.faces();
401 const auto& meshPoints = mesh_.points();
407 for (
const label facei : meshCells[celli])
409 for (
const label fpi : meshFaces[facei])
411 if (bb.contains(meshPoints[fpi]))
420 treeCellIDs.append(celli);
426 indexedOctree<treeDataCell>
tree
435 cellIds_.setSize(points_.size(), -1);
436 pointMask_.setSize(points_.size(),
false);
439 (void)mesh_.tetBasePtIs();
441 const auto& treeData =
tree.shapes();
448 label treeCelli =
tree.findInside(
pos);
452 reduce(proci, maxOp<label>());
454 pointMask_[pointi] = treeCelli != -1;
461 ? treeData.objectIndex(treeCelli)
467 if (errorOnPointNotFound_)
470 <<
"Position " <<
pos <<
" not found in mesh"
476 <<
"Position " <<
pos <<
" not found in mesh"
491 if (
field.size() != points_.size())
494 <<
"Inconsistent field sizes: input:" <<
field.size()
495 <<
" points:" << points_.size()
500 scalar sumFieldArea = 0;
502 for (
const face&
f : faces_)
505 scalar faceValue = 0;
506 for (
const label pti :
f)
509 if (!pointMask_[pti])
514 faceValue +=
field[pti];
519 scalar
area =
f.mag(points_);
525 return sumFieldArea/(sumArea + ROOTVSMALL);
537 if (!surfaceWriterPtr_)
544 surfaceWriterPtr_->isPointData(
true);
545 surfaceWriterPtr_->beginTime(time_);
546 surfaceWriterPtr_->open
550 (baseFileDir() /
name() /
"surfaces" /
"disk"),
553 surfaceWriterPtr_->nFields(4);
554 surfaceWriterPtr_->write(
"U",
U);
555 surfaceWriterPtr_->write(
"Ur", Ur);
565 if (!writePropellerPerformance_)
571 setRotationalSpeed();
573 const vector sumForce = forceEff();
574 const vector sumMoment = momentEff();
576 const scalar diameter = 2*radius_;
577 const scalar URef = URefPtr_->value(time_.timeOutputValue());
578 const scalar j = -URef/
mag(n_ + ROOTVSMALL)/diameter;
579 const scalar denom = rhoRef_*
sqr(n_)*
pow4(diameter);
580 const scalar kt = (sumForce & coordSysPtr_->e3())/denom;
582 -
sign(n_)*(sumMoment & coordSysPtr_->e3())/(denom*diameter);
587 auto&
os = propellerPerformanceFilePtr_();
589 writeCurrentTime(
os);
602 <<
" Revolutions per second, n : " << n_ <<
nl
603 <<
" Reference velocity, URef : " << URef <<
nl
604 <<
" Advance coefficient, J : " << j <<
nl
605 <<
" Thrust coefficient, Kt : " << kt <<
nl
606 <<
" Torque coefficient, 10*Kq : " << 10*kq <<
nl
607 <<
" Efficiency, etaO : " << etaO <<
nl
613 setResult(
"URef", URef);
630 auto&
os = wakeFilePtr_();
632 const pointField propPoints(coordSysPtr_->localPosition(points_));
634 mag(propPoints[1][0] - propPoints[0][0]) < SMALL ? 0 : 1;
635 const scalar rMax = propPoints.last()[0];
637 const scalar UzMean = meanSampleDiskField(
U.component(2));
639 writeHeaderValue(
os,
"Time", time_.timeOutputValue());
640 writeHeaderValue(
os,
"Reference velocity", URef);
641 writeHeaderValue(
os,
"Direction", coordSysPtr_->e3());
642 writeHeaderValue(
os,
"Wake", 1 - UzMean/URef);
644 writeCommented(
os,
"r/R");
645 writeTabbed(
os,
"alpha");
646 writeTabbed(
os,
"(x y z)");
647 writeTabbed(
os,
"(Ur Utheta Uz)");
650 for (label thetai = 0; thetai < nTheta_; ++thetai)
652 const scalar deg = 360*thetai/scalar(nTheta_);
654 for (label radiusi = 0; radiusi <= nRadial_; ++radiusi)
656 label pointi = radiusi*nTheta_ + thetai + offset;
658 if (radiusi == 0 && offset == 1)
664 if (pointMask_[pointi])
666 const scalar rR = propPoints[radiusi*nTheta_][0]/rMax;
669 << points_[pointi] <<
tab <<
U[pointi] <<
nl;
689 auto&
os = axialWakeFilePtr_();
691 const pointField propPoints(coordSysPtr_->localPosition(points_));
693 mag(propPoints[1][0] - propPoints[0][0]) < SMALL ? 0 : 1;
694 const scalar rMax = propPoints.last()[0];
696 writeHeaderValue(
os,
"Time", time_.timeOutputValue());
699 for (label radiusi = 0; radiusi <= nRadial_; ++radiusi)
701 label pointi = radiusi*nTheta_;
702 scalar r = propPoints[pointi][0];
703 os <<
tab <<
"r/R=" << r/rMax;
707 for (label thetai = 0; thetai < nTheta_; ++thetai)
709 os << 360*thetai/scalar(nTheta_);
711 for (label radiusi = 0; radiusi <= nRadial_; ++radiusi)
713 label pointi = radiusi*nTheta_ + thetai + offset;
715 if (radiusi == 0 && offset == 1)
721 if (pointMask_[pointi])
723 os <<
tab << 1 -
U[pointi][2]/URef;
727 os <<
tab <<
"undefined";
742 if (!writeWakeFields_)
748 if (
mag(URef) < ROOTSMALL)
751 <<
"Magnitude of reference velocity should be greater than zero"
759 const vectorField UrDisk(coordSysPtr_->localVector(UDisk));
762 writeSampleDiskSurface(UDisk, UrDisk, URef);
774 const Type& defaultValue
778 auto&
field = tfield.ref();
787 const label celli = cellIds_[pointi];
789 if (cellIds_[pointi] != -1)
791 const point& position = points_[pointi];
792 field[pointi] = interpolator().interpolate(position, celli);
816 rotationMode_(rotationMode::SPECIFIED),
819 writePropellerPerformance_(true),
820 propellerPerformanceFilePtr_(nullptr),
821 writeWakeFields_(true),
822 surfaceWriterPtr_(nullptr),
826 errorOnPointNotFound_(false),
830 interpolationScheme_(
"cell"),
831 wakeFilePtr_(nullptr),
832 axialWakeFilePtr_(nullptr),
864 radius_ =
dict.getScalar(
"radius");
865 URefPtr_.reset(Function1<scalar>::New(
"URef",
dict, &mesh_));
866 rotationMode_ = rotationModeNames_.get(
"rotationMode",
dict);
868 writePropellerPerformance_ =
869 dict.get<
bool>(
"writePropellerPerformance");
871 writeWakeFields_ =
dict.get<
bool>(
"writeWakeFields");
872 if (writeWakeFields_)
874 dict.readIfPresent(
"interpolationScheme", interpolationScheme_);
876 dict.readIfPresent(
"nanValue", nanValue_);
891 setCoordinateSystem(dict_);
893 if (writeWakeFields_)
895 setSampleDiskSurface(dict_);
905 if (writeWakeFields_)
913 coordSysPtr_->localVector
922 const scalar UzMean = meanSampleDiskField(UDisk.component(2));
924 setResult(
"UzMean", UzMean);
950 updateSampleDiskCells();
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
void clear() noexcept
Clear the addressed list, i.e. set the size to zero.
void append(const T &val)
Copy append an element to the end of this list.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
tmp< Field< cmptType > > component(const direction) const
Return a component field of the field.
Generic GeometricField class.
static const this_type & null() noexcept
List of MRF zones with IO functionality. MRF zones are specified by a list of dictionary entries,...
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
void resize_nocopy(const label len)
Adjust allocated size of list without necessarily.
static void listCombineReduce(UList< T > &values, CombineOp cop, const int tag=UPstream::msgType(), const int communicator=UPstream::worldComm)
Forwards to Pstream::listReduce with an in-place cop.
static void listReduce(UList< T > &values, BinaryOp bop, const int tag=UPstream::msgType(), const int communicator=UPstream::worldComm)
Reduce list elements (list must be equal size on all ranks), applying bop to each list element.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
void size(const label n)
Older name for setAddressableSize.
T & last()
Access last element of the list, position [size()-1].
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
static int myProcNo(const label communicator=worldComm)
Rank of this process in the communicator (starting from masterNo()). Negative if the process is not a...
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
static Form uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
void inflate(const scalar factor)
Expand box by factor*mag(span) in all dimensions.
A cylindrical coordinate system (r-theta-z). The coordinate system angle theta is always in radians.
Base class for coordinate system specification, the default coordinate system type is cartesian .
point globalPosition(const point &local) const
From local coordinate position to global (cartesian) position.
virtual const point & origin() const
Return origin.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T, or return the given default value. FatalIOError if it is found and the number of...
label getLabel(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Same as get<label>(const word&, keyType::option).
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Find an entry if present, and assign to T val. FatalIOError if it is found and the number of tokens i...
scalar getScalar(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Same as get<scalar>(const word&, keyType::option).
A face is a list of labels corresponding to mesh vertices.
Abstract base-class for Time/database function objects.
Computes forces and moments over a given list of patches by integrating pressure and viscous forces a...
virtual void calcForcesMoments()
Calculate forces and moments.
scalar rhoRef_
Reference density needed for incompressible calculations.
word UName_
Name of velocity field.
forces(const word &name, const Time &runTime, const dictionary &dict, const bool readFields=true)
Construct from name, Time and dictionary.
virtual vector forceEff() const
Return the total force.
autoPtr< coordinateSystem > coordSysPtr_
Coordinate system used when evaluating forces and moments.
virtual vector momentEff() const
Return the total moment.
virtual bool read(const dictionary &)
Read the function-object dictionary.
const fvMesh & mesh_
Reference to the fvMesh.
Calculates propeller performance and wake field properties.
word MRFName_
Name of MRF zone (if applicable).
label nTheta_
Number of surface divisions in theta direction.
boolList pointMask_
List of participating points (parallel reduced).
autoPtr< OFstream > propellerPerformanceFilePtr_
Propeller performance file.
void UpdateMesh(const mapPolyMesh &mpm)
labelList cellIds_
Surface point sample cell IDs.
void writeWake(const vectorField &U, const scalar URef)
Write the wake text file.
void writeWakeFields(const scalar URef)
Write the wake fields.
void createFiles()
Create output files.
rotationMode rotationMode_
Rotation mode.
void setSampleDiskSurface(const dictionary &dict)
Set the sample surface based on dictionary settings.
void writeAxialWake(const vectorField &U, const scalar URef)
Write the axial wake text file.
autoPtr< OFstream > wakeFilePtr_
Wake field file.
void writeSampleDiskSurface(const vectorField &U, const vectorField &Ur, const scalar URef)
Write the sample surface.
faceList faces_
Surface faces.
const volVectorField & U() const
Return the velocity field.
bool initialised_
Initialised flag.
dictionary dict_
Copy of dictionary used during construction.
propellerInfo(const propellerInfo &)=delete
No copy construct.
autoPtr< OFstream > axialWakeFilePtr_
Axial wake field file.
void writePropellerPerformance()
Write the wake fields.
void movePoints(const polyMesh &mesh)
Update for changes of mesh.
tmp< Field< Type > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &psi, const Type &defaultValue) const
Interpolate from the mesh onto the sample surface.
autoPtr< Function1< scalar > > URefPtr_
Reference velocity.
bool writePropellerPerformance_
Flag to write performance data.
bool errorOnPointNotFound_
Flag to raise an error if the sample point is not found in the mesh. Default = false to enable....
void setCoordinateSystem(const dictionary &dict)
Set the coordinate system.
scalar meanSampleDiskField(const scalarField &field) const
Return the area average of a field.
void setRotationalSpeed()
Set the rotational speed.
pointField points_
Surface points.
scalar radius_
Propeller radius.
void setSampleDiskGeometry(const coordinateSystem &coordSys, const scalar r1, const scalar r2, const scalar nTheta, const label nRadius, faceList &faces, pointField &points) const
Set the faces and points for the sample surface.
scalar nanValue_
Default value when a sample point is not found; default = scalar::min.
word interpolationScheme_
Interpolation scheme.
label nRadial_
Number of surface divisions in radial direction.
autoPtr< surfaceWriter > surfaceWriterPtr_
Surface writer.
scalar n_
Propeller speed (revolutions per second).
virtual bool execute()
Execute the function-object operations.
virtual bool write()
Write the function-object results.
void updateSampleDiskCells()
Set the sample cells corresponding to the sample points.
bool writeWakeFields_
Flag to write wake fields.
static const Enum< rotationMode > rotationModeNames_
virtual bool read(const dictionary &)
Read the function-object dictionary.
Reads fields from the time directories and adds them to the mesh database for further post-processing...
virtual const objectRegistry & obr() const
The region or sub-region registry being used.
void setResult(const word &entryName, const Type &value)
Add result.
const Time & time_
Reference to the time database.
const Time & time() const
Return time database.
fileName baseFileDir() const
Return the base directory for output.
virtual autoPtr< OFstream > newFileAtStartTime(const word &name) const
Return autoPtr to a new file using the simulation start time.
virtual void writeTabbed(Ostream &os, const string &str) const
Write a tabbed string to stream.
void writeHeaderValue(Ostream &os, const string &property, const Type &value) const
Write a (commented) header property and value pair.
virtual void writeBreak(Ostream &os) const
Write a break marker to the stream.
virtual void writeCommented(Ostream &os, const string &str) const
Write a commented string to stream.
virtual void writeCurrentTime(Ostream &os) const
Write the current time to stream.
virtual bool writeToFile() const
Flag to allow writing to file.
Non-pointer based hierarchical recursive searching.
static autoPtr< interpolation< Type > > New(const word &interpolationType, const GeometricField< Type, fvPatchField, volMesh > &psi)
Return a reference to the specified interpolation scheme.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Registry of regIOobjects.
A traits class, which is primarily used for primitives and vector-space.
Mesh consisting of general polyhedral cells.
static autoPtr< surfaceWriter > New(const word &writeType)
Select construct a surfaceWriter.
static dictionary formatOptions(const dictionary &dict, const word &formatName, const word &entryName="formatOptions")
Same as fileFormats::getFormatOptions.
A class for managing temporary objects.
Standard boundBox with extra functionality for use in octree.
bool contains(const vector &dir, const point &) const
Contains point (inside or on edge) and moving in direction.
Encapsulation of data needed to search in/for cells. Used to find the cell containing a point (e....
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const volScalarField & psi
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
#define DebugInfo
Report an information message using Foam::Info.
#define WarningInFunction
Report a warning using Foam::Warning.
constexpr scalar twoPi(2 *M_PI)
const wordList area
Standard area field types (scalar, vector, tensor, etc).
Function objects are OpenFOAM utilities to ease workflow configurations and enhance workflows.
dimensionedScalar pos(const dimensionedScalar &ds)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
GeometricField< vector, fvPatchField, volMesh > volVectorField
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const NameMatchPredicate &selectedFields, DynamicList< regIOobject * > &storedObjects)
Read the selected GeometricFields of the templated type and store on the objectRegistry.
dimensionedScalar sign(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static void writeHeader(Ostream &os, const word &fieldName)
List< face > faceList
List of faces.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void reduce(T &value, BinaryOp bop, const int tag=UPstream::msgType(), const int communicator=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce).
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
FlatOutput::OutputAdaptor< Container, Delimiters > flatOutput(const Container &obj, Delimiters delim)
Global flatOutput() function with specified output delimiters.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
dimensionedScalar pow4(const dimensionedScalar &ds)
errorManip< error > abort(error &err)
Field< vector > vectorField
Specialisation of Field<T> for vector.
IOerror FatalIOError
Error stream (stdout output on all processes), with additional 'FOAM FATAL IO ERROR' header text and ...
vector point
Point is a vector.
static constexpr const zero Zero
Global zero (0).
bool interpolate(const vector &p1, const vector &p2, const vector &o, vector &n, scalar l)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
vectorField pointField
pointField is a vectorField.
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a).
constexpr char tab
The tab '\t' character(0x09).
static bool initialised_(false)
Tree tree(triangles.begin(), triangles.end())
#define forAll(list, i)
Loop across all elements in list.