From b87bc88b5efd510e74ded76bc3c44719e7f3f5e3 Mon Sep 17 00:00:00 2001 From: Vuko Vukcevic Date: Tue, 13 Jun 2023 14:24:31 +0200 Subject: [PATCH] Address MR comments --- .../fanMomentumSource/fanMomentumSource.C | 179 ++++++++++-------- .../fanMomentumSource/fanMomentumSource.H | 145 ++++++++------ 2 files changed, 183 insertions(+), 141 deletions(-) diff --git a/src/fvOptions/sources/derived/fanMomentumSource/fanMomentumSource.C b/src/fvOptions/sources/derived/fanMomentumSource/fanMomentumSource.C index 813fc7dc72..92bd6e558f 100644 --- a/src/fvOptions/sources/derived/fanMomentumSource/fanMomentumSource.C +++ b/src/fvOptions/sources/derived/fanMomentumSource/fanMomentumSource.C @@ -26,12 +26,8 @@ License \*---------------------------------------------------------------------------*/ #include "fanMomentumSource.H" -#include "fvMatrices.H" -#include "DimensionedField.H" #include "IFstream.H" #include "addToRunTimeSelectionTable.H" -#include "TableFile.H" -#include "turbulenceModel.H" #include "turbulentTransportModel.H" #include "turbulentFluidThermoModel.H" @@ -78,14 +74,14 @@ void Foam::fv::fanMomentumSource::writeProps } -void Foam::fv::fanMomentumSource::initializeUpstreamFaces() +void Foam::fv::fanMomentumSource::initUpstreamFaces() { // First calculate the centre of gravity of the cell zone const vectorField& cellCentres = mesh_.cellCentres(); const scalarField& cellVolumes = mesh_.cellVolumes(); - vector centreGravityCellZone = vector::zero; - scalar cellZoneVolume = 0.; + vector centreGravityCellZone(Zero); + scalar cellZoneVolume = 0; for (const label celli : cells_) { const scalar cellVolume = cellVolumes[celli]; @@ -96,7 +92,16 @@ void Foam::fv::fanMomentumSource::initializeUpstreamFaces() reduce(centreGravityCellZone, sumOp()); reduce(cellZoneVolume, sumOp()); - centreGravityCellZone /= max(cellZoneVolume, SMALL); + if (cellZoneVolume < SMALL) + { + FatalErrorInFunction + << "Cannot initialize upstream faces because the total cell zone" + << " volume is zero or negative: " << cellZoneVolume << "." << nl + << "Check whether there are cells in fan momentum source cell zone." + << abort(FatalError); + } + + centreGravityCellZone /= cellZoneVolume; // Collect faces upstream of the centre of gavity const faceZone& fZone = mesh_.faceZones()[surroundingFaceZoneID_]; @@ -109,7 +114,7 @@ void Foam::fv::fanMomentumSource::initializeUpstreamFaces() { if ( - (flowDir_ & (faceCentreGravity[facei] - centreGravityCellZone)) < 0. + (flowDir_ & (faceCentreGravity[facei] - centreGravityCellZone)) < 0 ) { labelPair patchFaceInfo(-1, -1); @@ -181,10 +186,26 @@ void Foam::fv::fanMomentumSource::initializeUpstreamFaces() } -template Foam::scalar -Foam::fv::fanMomentumSource::calculateFlowRate(FlowRateFunctor f) const +Foam::fv::fanMomentumSource::calcFlowRate +( + const surfaceScalarField& phi, + const surfaceScalarField& rhof +) const { + // Sanity check to make sure we're not left in an inconsistent state because + // here we're operating on primitive (non-dimensional) fields + if (isNull(rhof) != (phi.dimensions() == dimVolume/dimTime)) + { + FatalErrorInFunction + << "Incorrect usage of the function. " << nl + << "For incompressible flow, pass surfaceScalarField::null for rhof" + << " and volumetric flux for phi." << nl + << "For compressible flow, pass face-interpolated density as rhof" + << " and mass flux for phi." + << exit(FatalError); + } + // Calculate the flow rate through the upstream faces scalarList phif(upstreamPatchFaceInfo_.size()); @@ -194,13 +215,35 @@ Foam::fv::fanMomentumSource::calculateFlowRate(FlowRateFunctor f) const { const labelPair& patchFaceInfo = upstreamPatchFaceInfo_[i]; + if (isNull(rhof)) + { + // Incompressible variant (phi = volumetric flux) + phif[i] = + patchFaceInfo.first() < 0 + ? phi[patchFaceInfo.second()] + : phi.boundaryField()[patchFaceInfo]; + } + else + { + // Compressible variant (phi = mass flux) + phif[i] = + patchFaceInfo.first() < 0 + ? phi[patchFaceInfo.second()]/ + rhof.internalField()[patchFaceInfo.second()] + : phi.boundaryField()[patchFaceInfo]/ + rhof.boundaryField()[patchFaceInfo]; + } + // Sign of the flux needs to be flipped if this is an internal face // whose owner is found in the cell zone - phif[i] = + if + ( patchFaceInfo.first() < 0 - && cellsInZones_.found(owners[patchFaceInfo.second()]) - ? -f(patchFaceInfo) - : f(patchFaceInfo); + && cellsInZones_.found(owners[patchFaceInfo.second()]) + ) + { + phif[i] *= -1; + } } return gSum(phif); @@ -217,19 +260,23 @@ Foam::fv::fanMomentumSource::fanMomentumSource ) : fv::cellSetOption(sourceName, modelType, dict, mesh), - fanCurve_(Function1::New("fanCurve", dict)), + upstreamPatchFaceInfo_(), + cellsInZones_(), + fanCurvePtr_(Function1::New("fanCurve", coeffs_, &mesh)), + UName_(coeffs_.get("U")), flowDir_(coeffs_.get("flowDir")), thickness_(coeffs_.get("thickness")), gradPFan_(0.0), - upstreamPatchFaceInfo_(), - rho_(nullptr) + surroundingFaceZoneID_(-1), + rho_(coeffs_.getOrDefault("rho", SMALL)), + useRefRho_(coeffs_.found("rho")) { // Skip all the checks if the source term has been deactivated // because there are no selected cells. // Such a situation typically occurs for multiple fluid regions if (fv::option::isActive()) { - const word faceZoneName = coeffs_.get("faceZone"); + const word faceZoneName = coeffs_.getWord("faceZone"); surroundingFaceZoneID_ = mesh_.faceZones().findZoneID(faceZoneName); @@ -242,34 +289,29 @@ Foam::fv::fanMomentumSource::fanMomentumSource << exit(FatalError); } + flowDir_.normalise(); if (mag(flowDir_) < SMALL) { - FatalErrorInFunction - << "Detected zero-vector for flowDir. Check your settings." - << exit(FatalError); + FatalIOErrorInFunction(coeffs_) + << "'flowDir' cannot be a zero-valued vector." + << exit(FatalIOError); } - else + + if (thickness_ < VSMALL) { - flowDir_ /= mag(flowDir_); + FatalIOErrorInFunction(coeffs_) + << "'thickness' cannot be non-positive." + << exit(FatalIOError); } - if (thickness_ <= 0.) + if (coeffs_.found("fields")) { - FatalErrorInFunction - << "The thickness of the fan model region cannot be negative" - << " or null." - << exit(FatalError); - } - - - coeffs_.readEntry("fields", fieldNames_); - - if (fieldNames_.size() != 1) - { - FatalErrorInFunction - << "Source can only be applied to a single field. Current " - << "settings are:" << fieldNames_ << exit(FatalError); + IOWarningInFunction(coeffs_) + << "Found 'fields' entry, which will be ignored. This function" + << " object will operate only on field 'U' = " << UName_ + << endl; } + fieldNames_.resize(1, UName_); fv::option::resetApplied(); @@ -288,12 +330,14 @@ Foam::fv::fanMomentumSource::fanMomentumSource Info<< " Initial pressure gradient = " << gradPFan_ << nl << endl; - if (coeffs_.found("rho")) + if (rho_ < SMALL) { - rho_.reset(new scalar(coeffs_.get("rho"))); + FatalIOErrorInFunction(coeffs_) + << "'rho' cannot be zero or negative." + << exit(FatalIOError); } - initializeUpstreamFaces(); + initUpstreamFaces(); } } @@ -314,7 +358,8 @@ void Foam::fv::fanMomentumSource::addSup mesh_.time().timeName(), mesh_, IOobject::NO_READ, - IOobject::NO_WRITE + IOobject::NO_WRITE, + IOobject::NO_REGISTER ), mesh_, dimensionedVector(eqn.dimensions()/dimVolume, Zero) @@ -330,31 +375,20 @@ void Foam::fv::fanMomentumSource::addSup << abort(FatalError); } - if (!rho_.good() || rho_.ref() < VSMALL) + if (!useRefRho_) { FatalErrorInFunction - << "You called incompressible addSup without or with " - << "zero value reference density." + << "You called incompressible addSup without providing a " + << "reference density. Add 'rho' entry to the dictionary." << abort(FatalError); } - // Lambda function passed as argument to calculateFlowRate - auto getVolumetricFlowRate = - [&phi](const labelPair& patchFaceInfo) - { - return - ( - patchFaceInfo.first() < 0 - ? phi[patchFaceInfo.second()] - : phi.boundaryField()[patchFaceInfo] - ); - }; - - const scalar flowRate = calculateFlowRate(getVolumetricFlowRate); + const scalar flowRate = calcFlowRate(phi, surfaceScalarField::null()); // Pressure drop for this flow rate // if flow rate is negative, pressure is clipped at the static pressure - gradPFan_ = fanCurve_->value(max(flowRate, scalar(0)))/thickness_/rho_.ref(); + gradPFan_ = + fanCurvePtr_->value(max(flowRate, scalar(0)))/thickness_/rho_; // Create the source term UIndirectList(Su, cells_) = flowDir_*gradPFan_; @@ -380,7 +414,8 @@ void Foam::fv::fanMomentumSource::addSup mesh_.time().timeName(), mesh_, IOobject::NO_READ, - IOobject::NO_WRITE + IOobject::NO_WRITE, + IOobject::NO_REGISTER ), mesh_, dimensionedVector(eqn.dimensions()/dimVolume, Zero) @@ -396,27 +431,12 @@ void Foam::fv::fanMomentumSource::addSup << abort(FatalError); } - const surfaceScalarField rhof = fvc::interpolate(rho); - - // Lambda function passed as argument to calculateFlowRate - auto getVolumetricFlowRate = - [&phi, &rhof](const labelPair& patchFaceInfo) - { - return - ( - patchFaceInfo.first() < 0 - ? phi[patchFaceInfo.second()]/ - rhof.internalField()[patchFaceInfo.second()] - : phi.boundaryField()[patchFaceInfo]/ - rhof.boundaryField()[patchFaceInfo] - ); - }; - - const scalar flowRate = calculateFlowRate(getVolumetricFlowRate); + const surfaceScalarField rhof(fvc::interpolate(rho)); + const scalar flowRate = calcFlowRate(phi, rhof); // Pressure drop for this flow rate // if flow rate is negative, pressure is clipped at the static pressure - gradPFan_ = fanCurve_->value(max(flowRate, scalar(0)))/thickness_; + gradPFan_ = fanCurvePtr_->value(max(flowRate, scalar(0)))/thickness_; // Create the source term UIndirectList(Su, cells_) = flowDir_*gradPFan_; @@ -424,7 +444,6 @@ void Foam::fv::fanMomentumSource::addSup eqn += Su; writeProps(gradPFan_, flowRate); - } diff --git a/src/fvOptions/sources/derived/fanMomentumSource/fanMomentumSource.H b/src/fvOptions/sources/derived/fanMomentumSource/fanMomentumSource.H index 6e95b5f5c0..e1baa7572f 100644 --- a/src/fvOptions/sources/derived/fanMomentumSource/fanMomentumSource.H +++ b/src/fvOptions/sources/derived/fanMomentumSource/fanMomentumSource.H @@ -24,61 +24,73 @@ License along with OpenFOAM. If not, see . Class - Foam::fanMomentumSource + Foam::fv::fanMomentumSource + +Group + grpFvOptionsSources Description - This source term models the action of a fan on the flow. + This source term models the action of a fan on the flow. It calculates + flow rate through a set of upstream faces of encompassing the cell zone. The + set of upstream faces is automatically calculated based on the flow + direction and the surrounding face zone. Based on the flow rate, the + pressure gradient is calculated based on the fan pressure curve and the + thickness of the fan (measurement section). The pressure gradient is then + added to the momentum equation. Usage + Minimal example by using \c constant/fvOptions: + \verbatim + + { + // Mandatory entries + type fanMomentumSource; + fanCurve >; + flowDir ; + thickness ; + cellZone ; + faceZone ; + + // Optional entries + gradient ; + rho ; + + // Inherited entries + ... + } + \endverbatim + + where the entries mean: \table - Property | Description | Required | Default - fanCurve | Pressure drop vs flow-rate | yes | - fields | Name of operand field | yes | - flowDir | Direction of the flow through the fan| yes | - faceZone | Name of upstream faceZone | yes | - thickness | Thickness of the fan | yes | - \endtable + Property | Description | Type | Reqd | Deflt + type | Type name: fanMomentumSource | word | yes | - + fanCurve | Pressure drop vs flow-rate | Function1\| yes | - + flowDir | Direction of the flow through the fan | vector | yes | - + cellZone | Cell zone representing the fan | word | yes | - + faceZone | Face zone around the cell zone | word | yes | - + thickness | Thickness of the fan [m] | scalar | yes | - + gradient | Initial pressure gradient | scalar | no | - + rho | Reference density for incompressible flow | scalar | no | - + \endtable The inherited entries are elaborated in: - - \link cellSetOption.H \endlink - - Example of the source term specification: - \verbatim - fanModel - { - type fanMomentumSource; - active on; - selectionMode cellZone; - cellZone fanCellZone; - faceZone fanUpstreamFaceZone; - fanCurve - { - type table; - file "constant/fanCurve"; - } - fields (U); - flowDir (1.0 0.0 0.0); - thickness 0.2; - } + - \link fvOption.H \endlink + - \link cellSetOption.H \endlink + - \link Function1.H \endlink Note - The fan curve needs to provide for a pressure drop expressed in Pascal and - is specified as a function of the volumetric flow rate in (m^3/s) + - The fan curve needs to be provded for a pressure drop in [Pa] and + is specified as a function of the volumetric flow rate in [m^3/s]. SourceFiles fanMomentumSource.C - \*---------------------------------------------------------------------------*/ -#ifndef fanMomentumSource_H -#define fanMomentumSource_H +#ifndef Foam_fv_fanMomentumSource_H +#define Foam_fv_fanMomentumSource_H -#include "autoPtr.H" -#include "fvMesh.H" -#include "volFields.H" #include "cellSetOption.H" -#include "Function1.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -88,7 +100,7 @@ namespace fv { /*---------------------------------------------------------------------------*\ - Class fanMomentumSource Declaration + Class fanMomentumSource Declaration \*---------------------------------------------------------------------------*/ class fanMomentumSource @@ -97,30 +109,36 @@ class fanMomentumSource { // Private Data - //- Volumetric flow rate vs. pressure drop table - autoPtr> fanCurve_; - - //- Direction of flow through the fan - vector flowDir_; - - //- Thickness of the fan - const scalar thickness_; - - //- Pressure drop - scalar gradPFan_; - //- Local list of pairs of upstream patch IDs and upstream face IDs // For internal faces, the upstream patch ID is -1 List upstreamPatchFaceInfo_; - //- Id for the surrounding face zone - label surroundingFaceZoneID_; - //- Cells in zone labelHashSet cellsInZones_; + //- Volumetric flow rate [m3] vs. pressure drop [Pa] table + autoPtr> fanCurvePtr_; + + //- Name of the velocity field this function object operates on + word UName_; + + //- Direction of flow through the fan + vector flowDir_; + + //- Thickness of the fan, in [m] + scalar thickness_; + + //- Pressure drop + scalar gradPFan_; + + //- Id for the surrounding face zone + label surroundingFaceZoneID_; + //- Reference density for incompressible cases - autoPtr rho_; + scalar rho_; + + //- Whether to use reference density + bool useRefRho_; // Private Member Functions @@ -128,13 +146,18 @@ class fanMomentumSource //- Write the pressure gradient to file (for restarts etc) void writeProps(const scalar gradP, const scalar flowRate) const; - //- Collect all faces upstream of the centre of gravity of the cell - // zone - void initializeUpstreamFaces(); + //- Collect all faces upstream of + //- the centre of gravity of the cell zone + void initUpstreamFaces(); - //- Calculate the volumetric flow rate - template - scalar calculateFlowRate(FlowRateFunctor f) const; + //- Return the volumetric flow rate given flux and face-interpolated + // density. For incompressible cases, use surfaceScalarField::null() + // for rhof + scalar calcFlowRate + ( + const surfaceScalarField& phi, + const surfaceScalarField& rhof + ) const; public: