momentum.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2018-2021 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "momentum.H"
29#include "fvMesh.H"
30#include "volFields.H"
31#include "cellSet.H"
32#include "cylindricalRotation.H"
33#include "emptyPolyPatch.H"
35
36// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37
38namespace Foam
39{
40namespace functionObjects
41{
44}
45}
46
47
48// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
49
50void Foam::functionObjects::momentum::purgeFields()
51{
52 obr_.checkOut(scopedName("momentum"));
53 obr_.checkOut(scopedName("angularMomentum"));
54 obr_.checkOut(scopedName("angularVelocity"));
55}
56
57
58template<class GeoField>
61(
62 const word& baseName,
63 const dimensionSet& dims,
64 bool registerObject
65) const
66{
67 return
69 (
71 (
72 scopedName(baseName),
73 time_.timeName(),
74 mesh_,
77 registerObject
78 ),
79 mesh_,
81 );
82}
83
84
85void Foam::functionObjects::momentum::calc()
86{
87 initialise();
88
89 // Ensure volRegion is properly up-to-date.
90 // Purge old fields if anything has changed (eg, mesh size etc)
92 {
93 purgeFields();
94 }
95
96 // When field writing is not enabled we need our local storage
97 // for the momentum and angular velocity fields
98 autoPtr<volVectorField> tmomentum, tAngularMom, tAngularVel;
99
100
101 // The base fields required
102 const auto& U = lookupObject<volVectorField>(UName_);
103 const auto* rhoPtr = findObject<volScalarField>(rhoName_);
104
105 const dimensionedScalar rhoRef("rho", dimDensity, rhoRef_);
106
107 // For quantities such as the mass-averaged angular velocity,
108 // we would calculate the mass per-cell here.
109
110 // tmp<volScalarField::Internal> tmass =
111 // (
112 // rhoPtr
113 // ? (mesh_.V() * (*rhoPtr))
114 // : (mesh_.V() * rhoRef)
115 // );
116
117
118 // Linear momentum
119 // ~~~~~~~~~~~~~~~
120
121 auto* pmomentum = getObjectPtr<volVectorField>(scopedName("momentum"));
122
123 if (!pmomentum)
124 {
125 tmomentum = newField<volVectorField>("momentum", dimVelocity*dimMass);
126 pmomentum = tmomentum.get(); // get(), not release()
127 }
128 auto& momentum = *pmomentum;
129
130 if (rhoPtr)
131 {
132 momentum.ref() = (U * mesh_.V() * (*rhoPtr));
133 }
134 else
135 {
136 momentum.ref() = (U * mesh_.V() * rhoRef);
137 }
138 momentum.correctBoundaryConditions();
139
140
141 // Angular momentum
142 // ~~~~~~~~~~~~~~~~
143
144 auto* pAngularMom =
145 getObjectPtr<volVectorField>(scopedName("angularMomentum"));
146
147 if (hasCsys_ && !pAngularMom)
148 {
149 tAngularMom =
150 newField<volVectorField>("angularMomentum", dimVelocity*dimMass);
151 pAngularMom = tAngularMom.get(); // get(), not release()
152 }
153 else if (!pAngularMom)
154 {
155 // Angular momentum not requested, but alias to normal momentum
156 // to simplify logic when calculating the summations
157 pAngularMom = pmomentum;
158 }
159 auto& angularMom = *pAngularMom;
160
161
162 // Angular velocity
163 // ~~~~~~~~~~~~~~~~
164
165 auto* pAngularVel =
166 getObjectPtr<volVectorField>(scopedName("angularVelocity"));
167
168 if (hasCsys_)
169 {
170 if (!pAngularVel)
171 {
172 tAngularVel =
173 newField<volVectorField>("angularVelocity", dimVelocity);
174 pAngularVel = tAngularVel.get(); // get(), not release()
175 }
176 auto& angularVel = *pAngularVel;
177
178
179 // Global to local
180
181 angularVel.primitiveFieldRef() =
182 csys_.invTransform(mesh_.cellCentres(), U.internalField());
183
184 angularVel.correctBoundaryConditions();
185
186 if (rhoPtr)
187 {
188 angularMom.ref() = (angularVel * mesh_.V() * (*rhoPtr));
189 }
190 else
191 {
192 angularMom.ref() = (angularVel * mesh_.V() * rhoRef);
193 }
194
195 angularMom.correctBoundaryConditions();
196 }
197
198
199 // Integrate the selection
200
201 sumMomentum_ = Zero;
202 sumAngularMom_ = Zero;
203
205 {
206 for (label celli=0; celli < mesh_.nCells(); ++celli)
207 {
208 sumMomentum_ += momentum[celli];
209 sumAngularMom_ += angularMom[celli];
210 }
211 }
212 else
213 {
214 for (const label celli : cellIDs())
215 {
216 sumMomentum_ += momentum[celli];
217 sumAngularMom_ += angularMom[celli];
218 }
219 }
220
221 reduce(sumMomentum_, sumOp<vector>());
222 reduce(sumAngularMom_, sumOp<vector>());
223}
224
225
226// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
227
229{
230 if (!writeToFile() || writtenHeader_)
231 {
232 return;
233 }
234
235 if (hasCsys_)
236 {
237 writeHeader(os, "Momentum, Angular Momentum");
238 writeHeaderValue(os, "origin", csys_.origin());
239 writeHeaderValue(os, "axis", csys_.e3());
240 }
241 else
242 {
243 writeHeader(os, "Momentum");
244 }
245
247 {
249 (
250 os,
251 "Selection " + regionTypeNames_[regionType_]
252 + " = " + regionName_
253 );
254 }
255
256 writeHeader(os, "");
257 writeCommented(os, "Time");
258 writeTabbed(os, "(momentum_x momentum_y momentum_z)");
259
260 if (hasCsys_)
261 {
262 writeTabbed(os, "(momentum_r momentum_rtheta momentum_axis)");
263 }
264
265 writeTabbed(os, "volume");
266 os << endl;
267
268 writtenHeader_ = true;
269}
270
271
273{
274 if (initialised_)
275 {
276 return;
277 }
278
279 if (!foundObject<volVectorField>(UName_))
280 {
282 << "Could not find U: " << UName_ << " in database"
283 << exit(FatalError);
284 }
285
286
287 const auto* pPtr = cfindObject<volScalarField>(pName_);
288
289 if (pPtr && pPtr->dimensions() == dimPressure)
290 {
291 // Compressible - rho is mandatory
292
293 if (!foundObject<volScalarField>(rhoName_))
294 {
296 << "Could not find rho:" << rhoName_
297 << exit(FatalError);
298 }
299 }
300
301 initialised_ = true;
302}
303
304
306{
307 if (log)
308 {
309 Info<< type() << " " << name() << " write:" << nl;
310
311 Info<< " Sum of Momentum";
312
314 {
315 Info<< ' ' << regionTypeNames_[regionType_]
316 << ' ' << regionName_;
317 }
318
319 Info<< " (volume " << volRegion::V() << ')' << nl
320 << " linear : " << sumMomentum_ << nl;
321
322 if (hasCsys_)
323 {
324 Info<< " angular : " << sumAngularMom_ << nl;
325 }
326
327 Info<< endl;
328 }
329
330 if (writeToFile())
331 {
332 writeCurrentTime(os);
333
334 os << tab << sumMomentum_;
335
336 if (hasCsys_)
337 {
338 os << tab << sumAngularMom_;
339 }
340
341 os << tab << volRegion::V() << endl;
342 }
343}
344
345
346// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
347
349(
350 const word& name,
351 const Time& runTime,
352 const dictionary& dict,
353 bool readFields
354)
355:
358 writeFile(mesh_, name, typeName, dict),
359 sumMomentum_(Zero),
360 sumAngularMom_(Zero),
361 UName_(),
362 pName_(),
363 rhoName_(),
364 rhoRef_(1.0),
365 csys_(),
366 hasCsys_(false),
367 writeMomentum_(false),
368 writeVelocity_(false),
369 writePosition_(false),
370 initialised_(false)
371{
372 if (readFields)
373 {
374 read(dict);
375 Log << endl;
376 }
377}
378
379
381(
382 const word& name,
383 const objectRegistry& obr,
384 const dictionary& dict,
385 bool readFields
386)
387:
390 writeFile(mesh_, name, typeName, dict),
391 sumMomentum_(Zero),
392 sumAngularMom_(Zero),
393 UName_(),
394 pName_(),
395 rhoName_(),
396 rhoRef_(1.0),
397 csys_(),
398 hasCsys_(false),
399 writeMomentum_(false),
400 writeVelocity_(false),
401 writePosition_(false),
402 initialised_(false)
403{
404 if (readFields)
405 {
406 read(dict);
407 Log << endl;
408 }
409}
410
411
412// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
413
415{
419
420 initialised_ = false;
421
422 Info<< type() << " " << name() << ":" << nl;
423
424 // Optional entries U and p
425 UName_ = dict.getOrDefault<word>("U", "U");
426 pName_ = dict.getOrDefault<word>("p", "p");
427 rhoName_ = dict.getOrDefault<word>("rho", "rho");
428 rhoRef_ = dict.getOrDefault<scalar>("rhoRef", 1.0);
429 hasCsys_ = dict.getOrDefault("cylindrical", false);
430
431 if (hasCsys_)
432 {
434 }
435
436 writeMomentum_ = dict.getOrDefault("writeMomentum", false);
437 writeVelocity_ = dict.getOrDefault("writeVelocity", false);
438 writePosition_ = dict.getOrDefault("writePosition", false);
439
440 Info<<"Integrating for selection: "
441 << regionTypeNames_[regionType_]
442 << " (" << regionName_ << ")" << nl;
443
444 if (writeMomentum_)
445 {
446 Info<< " Momentum fields will be written" << endl;
447
448 mesh_.objectRegistry::store
449 (
450 newField<volVectorField>("momentum", dimVelocity*dimMass)
451 );
452
453 if (hasCsys_)
454 {
455 mesh_.objectRegistry::store
456 (
457 newField<volVectorField>("angularMomentum", dimVelocity*dimMass)
458 );
459 }
460 }
461
462 if (hasCsys_)
463 {
464 if (writeVelocity_)
465 {
466 Info<< " Angular velocity will be written" << endl;
467
468 mesh_.objectRegistry::store
469 (
470 newField<volVectorField>("angularVelocity", dimVelocity)
471 );
472 }
473
474 if (writePosition_)
475 {
476 Info<< " Angular position will be written" << endl;
477 }
478 }
479
480 return true;
481}
482
483
485{
486 calc();
487
488 if (Pstream::master())
489 {
490 writeFileHeader(file());
491
492 writeValues(file());
493
494 Log << endl;
495 }
496
497 // Write state/results information
498 setResult("momentum_x", sumMomentum_[0]);
499 setResult("momentum_y", sumMomentum_[1]);
500 setResult("momentum_z", sumMomentum_[2]);
501
502 setResult("momentum_r", sumAngularMom_[0]);
503 setResult("momentum_rtheta", sumAngularMom_[1]);
504 setResult("momentum_axis", sumAngularMom_[2]);
505
506 return true;
507}
508
509
511{
512 if (writeMomentum_ || (hasCsys_ && (writeVelocity_ || writePosition_)))
513 {
514 Log << "Writing fields" << nl;
515
516 const volVectorField* fieldPtr;
517
518 fieldPtr = findObject<volVectorField>(scopedName("momentum"));
519 if (fieldPtr) fieldPtr->write();
520
521 fieldPtr = findObject<volVectorField>(scopedName("angularMomentum"));
522 if (fieldPtr) fieldPtr->write();
523
524 fieldPtr = findObject<volVectorField>(scopedName("angularVelocity"));
525 if (fieldPtr) fieldPtr->write();
526
527 if (hasCsys_ && writePosition_)
528 {
529 // Clunky, but currently no simple means of handling
530 // component-wise conversion and output
531
532 auto cyl_r = newField<volScalarField>("cyl_r", dimLength);
533 auto cyl_t = newField<volScalarField>("cyl_theta", dimless);
534 auto cyl_z = newField<volScalarField>("cyl_z", dimLength);
535
536 // Internal
537 {
538 const auto& pts = mesh_.cellCentres();
539 const label len = pts.size();
540
541 UList<scalar>& r = cyl_r->primitiveFieldRef(false);
542 UList<scalar>& t = cyl_t->primitiveFieldRef(false);
543 UList<scalar>& z = cyl_z->primitiveFieldRef(false);
544
545 for (label i=0; i < len; ++i)
546 {
547 point p(csys_.localPosition(pts[i]));
548
549 r[i] = p.x();
550 t[i] = p.y();
551 z[i] = p.z();
552 }
553 }
554
555 // Boundary
556 const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
557
558 forAll(pbm, patchi)
559 {
560 if (isA<emptyPolyPatch>(pbm[patchi]))
561 {
562 continue;
563 }
564
565 const auto& pts = pbm[patchi].faceCentres();
566 const label len = pts.size();
567
568 UList<scalar>& r = cyl_r->boundaryFieldRef(false)[patchi];
569 UList<scalar>& t = cyl_t->boundaryFieldRef(false)[patchi];
570 UList<scalar>& z = cyl_z->boundaryFieldRef(false)[patchi];
571
572 for (label i=0; i < len; ++i)
573 {
574 point p(csys_.localPosition(pts[i]));
575
576 r[i] = p.x();
577 t[i] = p.y();
578 z[i] = p.z();
579 }
580 }
581
582 cyl_r->write();
583 cyl_t->write();
584 cyl_z->write();
585 }
586 }
587
588 return true;
589}
590
591
593{
595 purgeFields(); // Mesh motion makes calculated fields dubious
596}
597
598
600{
602 purgeFields(); // Mesh motion makes calculated fields dubious
603}
604
605
606// ************************************************************************* //
#define Log
Definition: PDRblock.C:35
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
labelList cellIDs
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
virtual bool read()
Re-read model coefficients if they have changed.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:94
const complexVectorField & newField()
Definition: UOprocess.C:95
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
T * get() noexcept
Return pointer to managed object without nullptr checking.
Definition: autoPtr.H:152
A cylindrical coordinate system (r-theta-z). The coordinate system angle theta is always in radians.
Definition: cylindricalCS.H:74
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Dimension set for the base types, which can be used to implement rigorous dimension checking for alge...
Definition: dimensionSet.H:109
Generic dimensioned Type class.
Abstract base-class for Time/database function objects.
word scopedName(const word &name) const
Return a scoped (prefixed) name.
Specialization of Foam::functionObject for an Foam::fvMesh, providing a reference to the Foam::fvMesh...
Computes the natural logarithm of an input volScalarField.
Definition: log.H:230
Computes linear/angular momentum, reporting integral values and optionally writing the fields.
Definition: momentum.H:252
void initialise()
Initialise the fields.
Definition: momentum.C:272
void writeValues(Ostream &os)
Write momentum data.
Definition: momentum.C:305
virtual void writeFileHeader(Ostream &os)
Output file header information.
Definition: momentum.C:228
virtual bool execute()
Calculate and report the integral momentum.
Definition: momentum.C:484
virtual bool write()
Write the momentum, possibly angular momentum and velocity.
Definition: momentum.C:510
virtual bool read(const dictionary &)
Read the momentum data.
Definition: momentum.C:414
Reads fields from the time directories and adds them to the mesh database for further post-processing...
Definition: readFields.H:158
const objectRegistry & obr_
Reference to the region objectRegistry.
Volume (cell) region selection class.
Definition: volRegion.H:116
bool useAllCells() const noexcept
Use all cells, not the cellIDs.
Definition: volRegionI.H:31
scalar V() const
Return total volume of the selected region.
Definition: volRegionI.H:59
bool update()
Update the cached values as required.
Definition: volRegion.C:250
Base class for writing single files from the function objects.
Definition: writeFile.H:120
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:162
void movePoints()
Update for new mesh geometry.
void updateMesh()
Update for new mesh topology.
Registry of regIOobjects.
bool checkOut(regIOobject *io) const
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
virtual bool write(const bool valid=true) const
Write using setting from DB.
splitCell * master() const
Definition: splitCell.H:113
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
U
Definition: pEqn.H:72
volScalarField & p
engineTime & runTime
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
Namespace for OpenFOAM.
const dimensionSet dimPressure
const dimensionSet dimless
Dimensionless.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:52
static void writeHeader(Ostream &os, const word &fieldName)
const dimensionSet dimVelocity
messageStream Info
Information stream (stdout output on master, null elsewhere)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: MSwindows.C:598
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
const dimensionSet dimDensity
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:51
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
constexpr char tab
The tab '\t' character(0x09)
Definition: Ostream.H:52
Info<< "Reading mechanical properties\n"<< endl;IOdictionary mechanicalProperties(IOobject("mechanicalProperties", runTime.constant(), mesh, IOobject::MUST_READ_IF_MODIFIED, IOobject::NO_WRITE));const dictionary &rhoDict(mechanicalProperties.subDict("rho"));word rhoType(rhoDict.get< word >("type"));autoPtr< volScalarField > rhoPtr
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333