From 1078234f1871453cc36de841bea81899fb01ba5a Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Fri, 13 Jun 2025 14:33:48 +0200 Subject: [PATCH] ENH: Added new bladeForces function object Calculates thrust, drag, torque and lift/drag/pressure coefficients of single or multiple blades (eg, propeller, turbine blades) This function object differs from the propellerInfo and forces function objects in that all forces are calculated within the local cylindrical coordinate system, which yields thrust and drag values within the expected reference frame. The output comprises: - coefficients per radial bin - area-weighted total coefficients - integrated forces and torque For convenient post-processing, the results are also written in a VTK (.vtp) output format. All surface results are registered internally, which makes them available for other function objects. --- src/functionObjects/forces/Make/files | 1 + .../forces/bladeForces/bladeForces.H | 406 +++++ .../forces/bladeForces/bladeForces.cxx | 1464 +++++++++++++++++ .../forces/bladeForces/bladeForcesDict | 69 + .../pimpleFoam/RAS/propeller1/0.orig/U | 58 + .../pimpleFoam/RAS/propeller1/0.orig/epsilon | 46 + .../pimpleFoam/RAS/propeller1/0.orig/k | 46 + .../pimpleFoam/RAS/propeller1/0.orig/nut | 45 + .../pimpleFoam/RAS/propeller1/0.orig/p | 43 + .../pimpleFoam/RAS/propeller1/Allclean | 12 + .../pimpleFoam/RAS/propeller1/Allrun | 22 + .../pimpleFoam/RAS/propeller1/Allrun.post | 13 + .../pimpleFoam/RAS/propeller1/Allrun.pre | 28 + .../RAS/propeller1/constant/dynamicMeshDict | 34 + .../propeller1/constant/transportProperties | 22 + .../propeller1/constant/turbulenceProperties | 29 + .../RAS/propeller1/system/bladeForces | 88 + .../RAS/propeller1/system/blockMeshDict | 87 + .../RAS/propeller1/system/controlDict | 61 + .../RAS/propeller1/system/decomposeParDict | 27 + .../RAS/propeller1/system/fvSchemes | 58 + .../RAS/propeller1/system/fvSolution | 79 + .../RAS/propeller1/system/snappyHexMeshDict | 466 ++++++ .../system/surfaceFeatureExtractDict | 53 + .../system/surfaceFeatureExtractDict.defaults | 15 + .../resources/geometry/propeller1/hub.obj.gz | Bin 0 -> 12965 bytes .../geometry/propeller1/propeller.obj.gz | Bin 0 -> 1530851 bytes .../geometry/propeller1/transmission.obj.gz | Bin 0 -> 4779 bytes 28 files changed, 3272 insertions(+) create mode 100644 src/functionObjects/forces/bladeForces/bladeForces.H create mode 100644 src/functionObjects/forces/bladeForces/bladeForces.cxx create mode 100644 src/functionObjects/forces/bladeForces/bladeForcesDict create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/0.orig/U create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/0.orig/epsilon create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/0.orig/k create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/0.orig/nut create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/0.orig/p create mode 100755 tutorials/incompressible/pimpleFoam/RAS/propeller1/Allclean create mode 100755 tutorials/incompressible/pimpleFoam/RAS/propeller1/Allrun create mode 100755 tutorials/incompressible/pimpleFoam/RAS/propeller1/Allrun.post create mode 100755 tutorials/incompressible/pimpleFoam/RAS/propeller1/Allrun.pre create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/constant/dynamicMeshDict create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/constant/transportProperties create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/constant/turbulenceProperties create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/system/bladeForces create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/system/blockMeshDict create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/system/controlDict create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/system/decomposeParDict create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/system/fvSchemes create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/system/fvSolution create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/system/snappyHexMeshDict create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/system/surfaceFeatureExtractDict create mode 100644 tutorials/incompressible/pimpleFoam/RAS/propeller1/system/surfaceFeatureExtractDict.defaults create mode 100644 tutorials/resources/geometry/propeller1/hub.obj.gz create mode 100644 tutorials/resources/geometry/propeller1/propeller.obj.gz create mode 100644 tutorials/resources/geometry/propeller1/transmission.obj.gz diff --git a/src/functionObjects/forces/Make/files b/src/functionObjects/forces/Make/files index 4d8d3dcd63..8a9bc848a8 100644 --- a/src/functionObjects/forces/Make/files +++ b/src/functionObjects/forces/Make/files @@ -1,5 +1,6 @@ forces/forces.C forceCoeffs/forceCoeffs.C +bladeForces/bladeForces.cxx propellerInfo/propellerInfo.C LIB = $(FOAM_LIBBIN)/libforces diff --git a/src/functionObjects/forces/bladeForces/bladeForces.H b/src/functionObjects/forces/bladeForces/bladeForces.H new file mode 100644 index 0000000000..c0d3819d2e --- /dev/null +++ b/src/functionObjects/forces/bladeForces/bladeForces.H @@ -0,0 +1,406 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2024-2025 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +Class + Foam::functionObjects::bladeForces + +Group + grpForcesFunctionObjects + +Description + Computes forces and coefficients over a given list of patches by integrating + pressure and viscous forces and moments. + + Forces and moments are output in their constituent components + within the user-defined cylindrical coordinate system: + - thrust + - drag + - torque + + Operands: + \table + Operand | Type | Location + input | - | - + output file | dat | postProcessing/\/\/\s + \endtable + + where \c \s: + \verbatim + force.dat | Forces (thrust drag torque) + C_d.dat | Cd values (per bin) + C_l.dat | Cl values (per bin) + C_p.dat | Cp values (per bin) + \endverbatim + +Usage + Minimal example by using \c system/controlDict.functions: + \verbatim + + { + // Mandatory entries + type bladeForces; + libs (forces); + patches (); + + // Optional entries + writeFields ; + fieldsInterval ; + useNamePrefix ; + + // Mandatory entries + origin (0 0 0); + axis (1 0 0); + n 25; + + // Conditional optional entries + p ; + U ; + rho ; + rhoInf ; // enabled if rho=rhoInf + pRef ; + Uref ; + + // Inherited entries + ... + } + \endverbatim + + where the entries mean: + \table + Property | Description | Type | Reqd | Deflt + type | Type name: bladeForces | word | yes | - + libs | Library name: forces | word | yes | - + patches | Names of operand patches | wordRes | yes | - + writeFields | Flag to write surface with fields | bool | no | false + fieldsInterval | Frequency of writeInterval for writeFields | int | no | 0 + useNamePrefix | Flag to include prefix for field names | bool | no | false + outputName | Name for registered surface and VTP output | word | no | + origin | Origin of cylindrical coordinate system | vector | yes | - + axis | Axis of cylindrical coordinate system | vector | yes | - + n | Rotation speed [rev/sec] | scalar | yes | - + rpm | Rotation speed [rev/min] | scalar | no | - + p | Name of pressure field | word | no | p + U | Name of velocity field | word | no | U + rho | Name of density field | word | no | rho + rhoInf | Value of reference density | scalar | cndtnl | 1 + pRef | Value of reference pressure | scalar | cndtnl | 0 + Uref | Magnitude of inlet reference axial velocity | scalar | yes | - + + radius | The blade outer radius | scalar | no | 1 + nRadial | Divisions in radial direction | label | no | 10 + lefthand | Using a left-hand blade | bool | false | false + nearCellValue | Patch velocity extrapolated from fluid | bool | no | false + \endtable + + Experimental entries (may be removed in the future): + \table + Property | Description | Type | Reqd | Deflt + geometricVelocity | Patch velocity based on position | bool | no | false + mag.thrust | Ignore sign for thrust values | bool | false | false + mag.drag | Ignore sign for drag values | bool | false | false + \endtable + + The inherited entries are elaborated in: + - \link functionObject.H \endlink + - \link writeFile.H \endlink + - \link coordinateSystem.H \endlink + +Note + - For incompressible cases, set \c rho to \c rhoInf. + You will then be required to provide a \c rhoInf + value corresponding to the constant freestream density. + - \c writeControl and \c writeInterval entries of function + object do control when to output force/coefficients files and fields. + +SourceFiles + bladeForces.cxx + +\*---------------------------------------------------------------------------*/ + +#ifndef Foam_functionObjects_bladeForces_H +#define Foam_functionObjects_bladeForces_H + +#include "fvMeshFunctionObject.H" +#include "writeFile.H" +#include "cylindricalCS.H" +#include "volFieldsFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace functionObjects +{ + +/*---------------------------------------------------------------------------*\ + Class bladeForces Declaration +\*---------------------------------------------------------------------------*/ + +class bladeForces +: + public fvMeshFunctionObject, + public writeFile +{ + // Private Data + + // Results + + //- Sum of thrust forces (axial) + scalar sumThrust_; + + //- Sum of drag forces (tangential) + scalar sumDrag_; + + //- Sum of torque moments (tangential) + scalar sumTorque_; + + //- Overall blade area (single side) + scalar totalArea_; + + //- Overall drag coefficient (area-averaged of bin coefficients) + scalar totalCd_; + + //- Overall lift coefficient (area-averaged of bin coefficients) + scalar totalCl_; + + //- Overall pressure coefficient (area-averaged of bin coefficients) + scalar totalCp_; + + //- Blade area (single side) per radial band + scalarList bandArea_; + + //- Drag coefficient per radial band + scalarList bandCd_; + + //- Lift coefficient per radial band + scalarList bandCl_; + + //- Pressure coefficient per radial band + scalarList bandCp_; + + + // Geometric Information + + //- Selected operand patches + labelList patchIDs_; + + //- Max (or reference) blade radius + scalar refRadius_; + + //- Number of divisions in radial direction + label nRadialDiv_; + + //- Rotational speed (revolutions per second) + scalar revPerSec_; + + //- Cylindrical coordinate system used for force and torque + coordSystem::cylindrical cylCoord_; + + + // Read from dictionary + + //- Reference axial velocity + scalar Uref_; + + //- Reference pressure + scalar pRef_; + + //- Reference density (for incompressible) + scalar rhoRef_; + + //- Name of pressure field (default: "p") + word pName_; + + //- Name of velocity field (default: "U") + word UName_; + + //- Name of density field (default: "rho") + word rhoName_; + + //- Name for registred surface and VTP output. + // Default is the functionObject::name() + word outputName_; + + //- Internal counter, incremented when write() is called. + label writeCounter_; + + //- Write surface/fields interval. + // A value <= 1 means write surface/fields at each write(). + label fieldsInterval_; + + //- Flag of initialisation (internal) + bool initialised_; + + //- Flag to write force and moment fields + bool writeFields_; + + //- Using a left-hand blade (flips orientation of drag value) + bool lefthand_; + + //- Extrapolate velocity to patch + bool nearCellValue_; + + //- Blade speed based on position [experimental] + bool useGeometricVelocity_; + + //- Ignore sign for thrust values [experimental] + bool useMagThrust_; + + //- Ignore sign for drag values [experimental] + bool useMagDrag_; + + + // File streams + + //- File stream for forces, torque and summary + autoPtr forceFilePtr_; + + //- File stream for areas (in radial bins) + autoPtr areaFilePtr_; + + //- File stream for coefficient of drag (in radial bins) + autoPtr CdFilePtr_; + + //- File stream for coefficient of lift (in radial bins) + autoPtr ClFilePtr_; + + //- File stream for coefficient of pressure (in radial bins) + autoPtr CpFilePtr_; + + + // Private Member Functions + + //- Initialise containers and fields + void initialise(); + + //- Write surface (VTK format) and fields + void writeSurface() const; + + + // Evaluation + + //- Return the effective stress (viscous + turbulent) for patch + tmp devRhoReff + ( + const tensorField& gradUp, + const label patchi + ) const; + + //- Return dynamic viscosity field + tmp mu() const; + + //- Return rho if specified otherwise rhoRef + tmp rho() const; + + //- Return rho if specified otherwise rhoRef for patch + tmp rho(const label patchi) const; + + //- Return rhoRef if the pressure field is dynamic (i.e. p/rho), + //- otherwise return 1 + scalar rho(const volScalarField& p) const; + + + // I-O + + //- Create the integrated-data files + void createIntegratedDataFiles(); + + //- Write integrated data to files + void writeIntegratedDataFiles(); + + +public: + + //- Runtime type information + TypeName("bladeForces"); + + + // Constructors + + //- Construct from Time and dictionary + bladeForces + ( + const word& name, + const Time& runTime, + const dictionary& dict, + const bool readFields = true + ); + + //- Construct from objectRegistry and dictionary + bladeForces + ( + const word& name, + const objectRegistry& obr, + const dictionary& dict, + const bool readFields = true + ); + + //- No copy construct + bladeForces(const bladeForces&) = delete; + + //- No copy assignment + void operator=(const bladeForces&) = delete; + + + //- Destructor + virtual ~bladeForces() = default; + + + // Member Functions + + //- The integrated thrust force (axial) + scalar thrust() const noexcept { return sumThrust_; } + + //- The integrated drag force (tangential) + scalar drag() const noexcept { return sumDrag_; } + + //- The integrated torque + scalar torque() const noexcept { return sumTorque_; } + + //- Calculate forces, torque, coefficients + void calculate(); + + + //- Read the dictionary + virtual bool read(const dictionary& dict); + + //- Execute the function object + virtual bool execute(); + + //- Write to data files/fields and to streams + virtual bool write(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace functionObjects +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/functionObjects/forces/bladeForces/bladeForces.cxx b/src/functionObjects/forces/bladeForces/bladeForces.cxx new file mode 100644 index 0000000000..5914463b5c --- /dev/null +++ b/src/functionObjects/forces/bladeForces/bladeForces.cxx @@ -0,0 +1,1464 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2024-2025 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 . + +\*---------------------------------------------------------------------------*/ + +#include "bladeForces.H" +#include "fvcGrad.H" +#include "turbulentTransportModel.H" +#include "turbulentFluidThermoModel.H" +#include "cartesianCS.H" +#include "unitConversion.H" +#include "addToRunTimeSelectionTable.H" +#include "polySurfaceFields.H" +// #include "foamVtkSurfaceWriter.H" +// #include "surfaceWriter.H" +#include "vtkSurfaceWriter.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace functionObjects +{ + defineTypeNameAndDebug(bladeForces, 0); + addToRunTimeSelectionTable(functionObject, bladeForces, dictionary); +} +} + + +// Somewhat experimental... +#undef Foam_bladeForces_CalculateChord +#undef Foam_bladeForces_DebugAngles + +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Reduce with sum operation on list - MUST have identical length on all ranks! +// Uses en-mass reduction +template +static void reduceList_sum +( + UList& list, + const int tag = UPstream::msgType(), + const label comm = UPstream::worldComm +) +{ + Foam::reduce(list.data(), list.size(), sumOp(), tag, comm); +} + + +#ifdef Foam_bladeForces_CalculateChord + +// Reduce min/max values - MUST have identical length on all ranks! +// Uses en-mass reduction +static void reduceMinMax +( + UList& limits, + List& work, // work array + const int tag = UPstream::msgType(), + const label comm = UPstream::worldComm +) +{ + const label len = limits.size(); + + work.resize_nocopy(len); + + for (label i = 0; i < len; ++i) + { + work[i] = limits[i].min(); + } + + Foam::reduce(work.data(), work.size(), minOp(), tag, comm); + + for (label i = 0; i < len; ++i) + { + limits[i].min() = work[i]; + work[i] = limits[i].max(); + } + + Foam::reduce(work.data(), work.size(), maxOp(), tag, comm); + + for (label i = 0; i < len; ++i) + { + limits[i].max() = work[i]; + } +} + +#endif /* Foam_bladeForces_CalculateChord */ + +} // End namespace Foam + +// Cannot debug blade rotation angles without chord calculations +#ifndef Foam_bladeForces_CalculateChord +#undef Foam_bladeForces_DebugAngles +#endif + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +void Foam::functionObjects::bladeForces::initialise() +{ + if (initialised_) + { + return; + } + + const bool isCompressible = (rhoName_ != "rhoInf"); + + if + ( + !foundObject(UName_) + || !foundObject(pName_) + || (isCompressible && !foundObject(rhoName_)) + ) + { + FatalErrorInFunction + << "Could not find velocity '" << UName_ + << "' or pressure '" << pName_ << '\''; + + if (isCompressible) + { + FatalError + << " or density '" << rhoName_ << '\''; + } + + FatalError + << " in database" << endl + << exit(FatalError); + } + + initialised_ = true; +} + + +Foam::tmp +Foam::functionObjects::bladeForces::devRhoReff +( + const tensorField& gradUp, + const label patchi +) const +{ + typedef incompressible::turbulenceModel icoModel; + typedef compressible::turbulenceModel cmpModel; + + if (const auto* turbp = cfindObject(icoModel::propertiesName)) + { + const auto& turb = *turbp; + + return -rho(patchi)*turb.nuEff(patchi)*devTwoSymm(gradUp); + } + + if (const auto* turbp = cfindObject(cmpModel::propertiesName)) + { + const auto& turb = *turbp; + + return -turb.muEff(patchi)*devTwoSymm(gradUp); + } + + if (const auto* thermop = cfindObject(fluidThermo::dictName)) + { + const auto& thermo = *thermop; + + return -thermo.mu(patchi)*devTwoSymm(gradUp); + } + + if (const auto* props = cfindObject("transportProperties")) + { + const auto& laminarT = *props; + + return -rho(patchi)*laminarT.nu(patchi)*devTwoSymm(gradUp); + } + + if (const auto* props = cfindObject("transportProperties")) + { + const dimensionedScalar nu("nu", dimViscosity, *props); + + return -rho(patchi)*nu.value()*devTwoSymm(gradUp); + } + + // Failed to resolve any model + FatalErrorInFunction + << "No valid model for viscous stress calculation" + << exit(FatalError); + + return nullptr; +} + + +Foam::tmp Foam::functionObjects::bladeForces::mu() const +{ + if (const auto* thermop = cfindObject(basicThermo::dictName)) + { + const auto& thermo = *thermop; + + return thermo.mu(); + } + + if (const auto* props = cfindObject("transportProperties")) + { + const auto& laminarT = *props; + + return rho()*laminarT.nu(); + } + + if (const auto* props = cfindObject("transportProperties")) + { + const dimensionedScalar nu("nu", dimViscosity, *props); + + return rho()*nu; + } + + // Failed to resolve any model + FatalErrorInFunction + << "No valid model for dynamic viscosity calculation" + << exit(FatalError); + + return nullptr; +} + + +Foam::tmp Foam::functionObjects::bladeForces::rho() const +{ + if (rhoName_ == "rhoInf") + { + return volScalarField::New + ( + "rho", + mesh_, + dimensionedScalar(dimDensity, rhoRef_) + ); + } + + return lookupObject(rhoName_); +} + + +Foam::tmp +Foam::functionObjects::bladeForces::rho(const label patchi) const +{ + if (rhoName_ == "rhoInf") + { + return tmp::New + ( + mesh_.boundary()[patchi].size(), + rhoRef_ + ); + } + + const auto& rho = lookupObject(rhoName_); + return rho.boundaryField()[patchi]; +} + + +Foam::scalar Foam::functionObjects::bladeForces::rho +( + const volScalarField& p +) const +{ + if (p.dimensions() == dimPressure) + { + return 1; + } + + if (rhoName_ != "rhoInf") + { + FatalErrorInFunction + << "Dynamic pressure is expected but kinematic is provided." + << exit(FatalError); + } + + return rhoRef_; +} + + +void Foam::functionObjects::bladeForces::createIntegratedDataFiles() +{ + if (!forceFilePtr_) + { + forceFilePtr_ = newFileAtStartTime("force"); + auto& os = forceFilePtr_(); + + writeHeader(os, "Force"); + writeHeader(os, ""); + writeCommented(os, "Time"); + writeTabbed(os, "thrust"); + writeTabbed(os, "drag"); + writeTabbed(os, "torque"); + writeTabbed(os, "Cd"); + writeTabbed(os, "Cl"); + writeTabbed(os, "Cp"); + writeTabbed(os, "area"); + os << endl; + } + + const auto writeCoeffsHeader = + [&](auto& os, const auto& longName, const auto& shortName) + { + writeHeader(os, longName); + writeHeader(os, ""); + writeCommented(os, "number of radial bins:"); + os << ' ' << nRadialDiv_ << nl; + + // General formatting (not fixed/scientific) for the bin limits: + const auto oldFlags = os.flags(); + os.unsetf(std::ios_base::floatfield); + + writeCommented(os, "bin lower limits:"); + for (label i = 0; i < nRadialDiv_; ++i) + { + os << ' ' << float(i)/nRadialDiv_; + } + os << nl; + writeCommented(os, "bin upper limits:"); + for (label i = 1; i <= nRadialDiv_; ++i) + { + os << ' ' << float(i)/nRadialDiv_; + } + os << nl; + writeHeader(os, ""); + writeCommented(os, "Time"); + writeTabbed(os, shortName); + os << endl; + + os.setf(oldFlags); + }; + + + if (!areaFilePtr_) + { + areaFilePtr_ = newFileAtStartTime("area"); + auto& os = areaFilePtr_(); + + writeCoeffsHeader(os, "Area", "area(per bin)"); + } + + if (!CdFilePtr_) + { + CdFilePtr_ = newFileAtStartTime("C_d"); + auto& os = CdFilePtr_(); + + writeCoeffsHeader(os, "Drag coefficient", "Cd(per bin)"); + } + + if (!ClFilePtr_) + { + ClFilePtr_ = newFileAtStartTime("C_l"); + auto& os = ClFilePtr_(); + + writeCoeffsHeader(os, "Lift coefficient", "Cl(per bin)"); + } + + if (!CpFilePtr_) + { + CpFilePtr_ = newFileAtStartTime("C_p"); + auto& os = CpFilePtr_(); + + writeCoeffsHeader(os, "Pressure coefficient", "Cp(per bin)"); + } +} + + +void Foam::functionObjects::bladeForces::writeIntegratedDataFiles() +{ + { + auto& os = forceFilePtr_(); + + writeCurrentTime(os); + + writeValue(os, (useMagThrust_ ? Foam::mag(sumThrust_) : sumThrust_)); + writeValue(os, (useMagDrag_ ? Foam::mag(sumDrag_) : sumDrag_)); + writeValue(os, sumTorque_); + writeValue(os, totalCd_); + writeValue(os, totalCl_); + writeValue(os, totalCp_); + writeValue(os, totalArea_); + os << endl; + } + + { + auto& os = areaFilePtr_(); + + writeCurrentTime(os); + + for (const auto& val : bandArea_) + { + writeValue(os, val); + } + os << endl; + } + + { + auto& os = CdFilePtr_(); + + writeCurrentTime(os); + + for (const auto& val : bandCd_) + { + writeValue(os, val); + } + os << endl; + } + + { + auto& os = ClFilePtr_(); + + writeCurrentTime(os); + + for (const auto& val : bandCl_) + { + writeValue(os, val); + } + os << endl; + } + + { + auto& os = CpFilePtr_(); + + writeCurrentTime(os); + + for (const auto& val : bandCp_) + { + writeValue(os, val); + } + os << endl; + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::functionObjects::bladeForces::bladeForces +( + const word& name, + const Time& runTime, + const dictionary& dict, + bool readFields +) +: + fvMeshFunctionObject(name, runTime, dict), + writeFile(mesh_, name), + + sumThrust_(0), + sumDrag_(0), + sumTorque_(0), + + totalArea_(0), + totalCd_(0), + totalCl_(0), + totalCp_(0), + + refRadius_(1), + nRadialDiv_(10), + revPerSec_(1), + + Uref_(1), + pRef_(0), + rhoRef_(1), + pName_("p"), + UName_("U"), + rhoName_("rho"), + outputName_("blade"), // overwritten with name() or user-specified + writeCounter_(0), + fieldsInterval_(0), + initialised_(false), + writeFields_(false), + lefthand_(false), + nearCellValue_(false), + useGeometricVelocity_(false), + useMagThrust_(false), + useMagDrag_(false) +{ + if (readFields) + { + read(dict); + Log << endl; + } +} + + +Foam::functionObjects::bladeForces::bladeForces +( + const word& name, + const objectRegistry& obr, + const dictionary& dict, + bool readFields +) +: + fvMeshFunctionObject(name, obr, dict), + writeFile(mesh_, name), + + sumThrust_(0), + sumDrag_(0), + sumTorque_(0), + + totalArea_(0), + totalCd_(0), + totalCl_(0), + totalCp_(0), + + refRadius_(1), + nRadialDiv_(10), + revPerSec_(1), + + Uref_(1), + pRef_(0), + rhoRef_(1), + pName_("p"), + UName_("U"), + rhoName_("rho"), + outputName_("blade"), // overwritten with name() or user-specified + writeCounter_(0), + fieldsInterval_(0), + initialised_(false), + writeFields_(false), + lefthand_(false), + nearCellValue_(false), + useGeometricVelocity_(false), + useMagThrust_(false), + useMagDrag_(false) +{ + if (readFields) + { + read(dict); + Log << endl; + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::functionObjects::bladeForces::read(const dictionary& dict) +{ + const polyBoundaryMesh& pbm = mesh_.boundaryMesh(); + + if (!fvMeshFunctionObject::read(dict) || !writeFile::read(dict)) + { + return false; + } + + initialised_ = false; + + // Reset to defaults + Uref_ = 1; + pRef_ = 0; + rhoRef_ = 1; + outputName_ = this->name(); + writeCounter_ = 0; + fieldsInterval_ = 0; + writeFields_ = false; + lefthand_ = false; + nearCellValue_ = false; + useGeometricVelocity_ = false; + useMagThrust_ = false; + useMagDrag_ = false; + + + Info<< type() << ' ' << name() << ':' << nl; + + // Can also use pbm.indices(wordRes), but no warnings... + patchIDs_ = pbm.patchSet(dict.get("patches")).sortedToc(); + + // With "excludePatches": + // { + // wordRes selectPatches, blockPatches; + // dict.readEntry("patches", selectPatches); + // dict.readIfPresent("excludePatches", blockPatches) + // + // if (!selectPatches.empty() && !excludePatches.empty()) + // { + // // Avoid default "select all" if both are empty + // patchIDs_ = pbm.indices(selectPatches, excludePatches); + // } + // else + // { + // // Can also use pbm.indices(selectPatches), but no warnings... + // patchIDs_ = pbm.patchSet(selectPatches).sortedToc(); + // } + // } + + refRadius_ = dict.getOrDefault("radius", 1); + nRadialDiv_ = dict.getOrDefault