Loading...
Searching...
No Matches
ITHACAforces.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ██╗████████╗██╗ ██╗ █████╗ ██████╗ █████╗ ███████╗██╗ ██╗
3 ██║╚══██╔══╝██║ ██║██╔══██╗██╔════╝██╔══██╗ ██╔════╝██║ ██║
4 ██║ ██║ ███████║███████║██║ ███████║█████╗█████╗ ██║ ██║
5 ██║ ██║ ██╔══██║██╔══██║██║ ██╔══██║╚════╝██╔══╝ ╚██╗ ██╔╝
6 ██║ ██║ ██║ ██║██║ ██║╚██████╗██║ ██║ ██║ ╚████╔╝
7 ╚═╝ ╚═╝ ╚═╝ ╚═╝╚═╝ ╚═╝ ╚═════╝╚═╝ ╚═╝ ╚═╝ ╚═══╝
8
9 * In real Time Highly Advanced Computational Applications for Finite Volumes
10 * Copyright (C) 2017 by the ITHACA-FV authors
11-------------------------------------------------------------------------------
12
13License
14 This file is part of ITHACA-FV
15
16 ITHACA-FV is free software: you can redistribute it and/or modify
17 it under the terms of the GNU Lesser General Public License as published by
18 the Free Software Foundation, either version 3 of the License, or
19 (at your option) any later version.
20
21 ITHACA-FV is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU Lesser General Public License for more details.
25
26 You should have received a copy of the GNU Lesser General Public License
27 along with ITHACA-FV. If not, see <http://www.gnu.org/licenses/>.
28
29\*---------------------------------------------------------------------------*/
30
31#if OPENFOAM >= 1812
32
33#include "ITHACAforces18.H"
34#include "fvcGrad.H"
35#include "porosityModel.H"
36#include "turbulentTransportModel.H"
37#include "turbulentFluidThermoModel.H"
38#include "addToRunTimeSelectionTable.H"
39#include "cartesianCS.H"
40
41// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42
43namespace Foam
44{
45namespace functionObjects
46{
47defineTypeNameAndDebug(ITHACAforces, 0);
48addToRunTimeSelectionTable(functionObject, ITHACAforces, dictionary);
49}
50}
51
52
53// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
54
55Foam::word Foam::functionObjects::ITHACAforces::fieldName(
56 const word& name) const
57{
58 return this->name() + ":" + name;
59}
60
61
62void Foam::functionObjects::ITHACAforces::createFiles()
63{
64 // Note: Only possible to create bin files after bins have been initialised
65 if (writeToFile() && !forceFilePtr_.valid())
66 {
67 forceFilePtr_ = createFile("force");
68 writeIntegratedHeader("Force", forceFilePtr_());
69 momentFilePtr_ = createFile("moment");
70 writeIntegratedHeader("Moment", momentFilePtr_());
71
72 if (nBin_ > 1)
73 {
74 forceBinFilePtr_ = createFile("forceBin");
75 writeBinHeader("Force", forceBinFilePtr_());
76 momentBinFilePtr_ = createFile("momentBin");
77 writeBinHeader("Moment", momentBinFilePtr_());
78 }
79
80 if (localSystem_)
81 {
82 localForceFilePtr_ = createFile("localForce");
83 writeIntegratedHeader("Force", localForceFilePtr_());
84 localMomentFilePtr_ = createFile("localMoment");
85 writeIntegratedHeader("Moment", localMomentFilePtr_());
86
87 if (nBin_ > 1)
88 {
89 localForceBinFilePtr_ = createFile("localForceBin");
90 writeBinHeader("Force", localForceBinFilePtr_());
91 localMomentBinFilePtr_ = createFile("localMomentBin");
92 writeBinHeader("Moment", localMomentBinFilePtr_());
93 }
94 }
95 }
96}
97
98
99void Foam::functionObjects::ITHACAforces::writeIntegratedHeader
100(
101 const word& header,
102 Ostream& os
103) const
104{
105 writeHeader(os, header);
106 writeHeaderValue(os, "CofR", coordSys_.origin());
107 writeHeader(os, "");
108 writeCommented(os, "Time");
109 writeTabbed(os, "(total_x total_y total_z)");
110 writeTabbed(os, "(pressure_x pressure_y pressure_z)");
111 writeTabbed(os, "(viscous_x viscous_y viscous_z)");
112
113 if (porosity_)
114 {
115 writeTabbed(os, "(porous_x porous_y porous_z)");
116 }
117
118 os << endl;
119}
120
121
122void Foam::functionObjects::ITHACAforces::writeBinHeader
123(
124 const word& header,
125 Ostream& os
126) const
127{
128 writeHeader(os, header + " bins");
129 writeHeaderValue(os, "bins", nBin_);
130 writeHeaderValue(os, "start", binMin_);
131 writeHeaderValue(os, "delta", binDx_);
132 writeHeaderValue(os, "direction", binDir_);
133 vectorField binPoints(nBin_);
134 writeCommented(os, "x co-ords :");
135 forAll(binPoints, pointi)
136 {
137 binPoints[pointi] = (binMin_ + (pointi + 1) * binDx_) * binDir_;
138 os << tab << binPoints[pointi].x();
139 }
140
141 os << nl;
142 writeCommented(os, "y co-ords :");
143 forAll(binPoints, pointi)
144 {
145 os << tab << binPoints[pointi].y();
146 }
147
148 os << nl;
149 writeCommented(os, "z co-ords :");
150 forAll(binPoints, pointi)
151 {
152 os << tab << binPoints[pointi].z();
153 }
154
155 os << nl;
156 writeHeader(os, "");
157 writeCommented(os, "Time");
158
159 for (label j = 0; j < nBin_; j++)
160 {
161 const word jn(Foam::name(j) + ':');
162 os << tab << jn << "(total_x total_y total_z)"
163 << tab << jn << "(pressure_x pressure_y pressure_z)"
164 << tab << jn << "(viscous_x viscous_y viscous_z)";
165
166 if (porosity_)
167 {
168 os << tab << jn << "(porous_x porous_y porous_z)";
169 }
170 }
171
172 os << endl;
173}
174
175
176
177void Foam::functionObjects::ITHACAforces::initialise()
178{
179 if (initialised_)
180 {
181 return;
182 }
183
184 if (directForceDensity_)
185 {
186 if (!foundObject<volVectorField>(fDName_))
187 {
188 FatalErrorInFunction
189 << "Could not find " << fDName_ << " in database"
190 << exit(FatalError);
191 }
192 }
193 else
194 {
195 if
196 (
197 !foundObject<volVectorField>(UName_)
198 || !foundObject<volScalarField>(pName_)
199 )
200 {
201 FatalErrorInFunction
202 << "Could not find U: " << UName_ << " or p:" << pName_
203 << " in database"
204 << exit(FatalError);
205 }
206
207 if (rhoName_ != "rhoInf" && !foundObject<volScalarField>(rhoName_))
208 {
209 FatalErrorInFunction
210 << "Could not find rho:" << rhoName_
211 << exit(FatalError);
212 }
213 }
214
215 initialiseBins();
216 initialised_ = true;
217}
218
219
220void Foam::functionObjects::ITHACAforces::initialiseBins()
221{
222 if (nBin_ > 1)
223 {
224 const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
225 // Determine extents of patches
226 binMin_ = GREAT;
227 scalar binMax = -GREAT;
228
229 for (const label patchi : patchSet_)
230 {
231 const polyPatch& pp = pbm[patchi];
232 scalarField d(pp.faceCentres() & binDir_);
233 binMin_ = min(min(d), binMin_);
234 binMax = max(max(d), binMax);
235 }
236
237 // Include porosity
238 if (porosity_)
239 {
240 const HashTable<const porosityModel*> models =
241 obr_.lookupClass<porosityModel>();
242 const scalarField dd(mesh_.C() & binDir_);
243 forAllConstIter(HashTable<const porosityModel*>, models, iter)
244 {
245 const porosityModel& pm = * iter();
246 const labelList& cellZoneIDs = pm.cellZoneIDs();
247
248 for (const label zonei : cellZoneIDs)
249 {
250 const cellZone& cZone = mesh_.cellZones()[zonei];
251 const scalarField d(dd, cZone);
252 binMin_ = min(min(d), binMin_);
253 binMax = max(max(d), binMax);
254 }
255 }
256 }
257
258 reduce(binMin_, minOp<scalar>());
259 reduce(binMax, maxOp<scalar>());
260 // Slightly boost binMax so that region of interest is fully
261 // within bounds
262 binMax = 1.0001 * (binMax - binMin_) + binMin_;
263 binDx_ = (binMax - binMin_) / scalar(nBin_);
264 // Create the bin points used for writing
265 binPoints_.setSize(nBin_);
266 forAll(binPoints_, i)
267 {
268 binPoints_[i] = (i + 0.5) * binDir_ * binDx_;
269 }
270 }
271
272 // Allocate storage for forces and moments
273 forAll(force_, i)
274 {
275 force_[i].setSize(nBin_, vector::zero);
276 moment_[i].setSize(nBin_, vector::zero);
277 }
278}
279
280
281void Foam::functionObjects::ITHACAforces::resetFields()
282{
283 force_[0] = Zero;
284 force_[1] = Zero;
285 force_[2] = Zero;
286 moment_[0] = Zero;
287 moment_[1] = Zero;
288 moment_[2] = Zero;
289
290 if (writeFields_)
291 {
292 volVectorField& force =
293 lookupObjectRef<volVectorField>(fieldName("force"));
294 force == dimensionedVector(force.dimensions(), Zero);
295 volVectorField& moment =
296 lookupObjectRef<volVectorField>(fieldName("moment"));
297 moment == dimensionedVector(moment.dimensions(), Zero);
298 }
299}
300
301
302Foam::tmp<Foam::volSymmTensorField>
303Foam::functionObjects::ITHACAforces::devRhoReff() const
304{
305 typedef compressible::turbulenceModel cmpTurbModel;
306 typedef incompressible::turbulenceModel icoTurbModel;
307
308 if (foundObject<cmpTurbModel>(cmpTurbModel::propertiesName))
309 {
310 const cmpTurbModel& turb =
311 lookupObject<cmpTurbModel>(cmpTurbModel::propertiesName);
312 return turb.devRhoReff();
313 }
314 else if (foundObject<icoTurbModel>(icoTurbModel::propertiesName))
315 {
316 const incompressible::turbulenceModel& turb =
317 lookupObject<icoTurbModel>(icoTurbModel::propertiesName);
318 return rho() * turb.devReff();
319 }
320 else if (foundObject<fluidThermo>(fluidThermo::dictName))
321 {
322 const fluidThermo& thermo =
323 lookupObject<fluidThermo>(fluidThermo::dictName);
324 const volVectorField& U = lookupObject<volVectorField>(UName_);
325 return -thermo.mu() * dev(twoSymm(fvc::grad(U)));
326 }
327 else if (foundObject<transportModel>("transportProperties"))
328 {
329 const transportModel& laminarT =
330 lookupObject<transportModel>("transportProperties");
331 const volVectorField& U = lookupObject<volVectorField>(UName_);
332 return -rho() * laminarT.nu() * dev(twoSymm(fvc::grad(U)));
333 }
334 else if (foundObject<dictionary>("transportProperties"))
335 {
336 const dictionary& transportProperties =
337 lookupObject<dictionary>("transportProperties");
338 dimensionedScalar nu("nu", dimViscosity, transportProperties);
339 const volVectorField& U = lookupObject<volVectorField>(UName_);
340 return -rho() * nu * dev(twoSymm(fvc::grad(U)));
341 }
342 else
343 {
344 FatalErrorInFunction
345 << "No valid model for viscous stress calculation"
346 << exit(FatalError);
347 return volSymmTensorField::null();
348 }
349}
350
351
352Foam::tmp<Foam::volScalarField> Foam::functionObjects::ITHACAforces::mu() const
353{
354 if (foundObject<fluidThermo>(basicThermo::dictName))
355 {
356 const fluidThermo& thermo =
357 lookupObject<fluidThermo>(basicThermo::dictName);
358 return thermo.mu();
359 }
360 else if (foundObject<transportModel>("transportProperties"))
361 {
362 const transportModel& laminarT =
363 lookupObject<transportModel>("transportProperties");
364 return rho() * laminarT.nu();
365 }
366 else if (foundObject<dictionary>("transportProperties"))
367 {
368 const dictionary& transportProperties =
369 lookupObject<dictionary>("transportProperties");
370 dimensionedScalar nu
371 (
372 "nu",
373 dimViscosity,
374 transportProperties.lookup("nu")
375 );
376 return rho() * nu;
377 }
378 else
379 {
380 FatalErrorInFunction
381 << "No valid model for dynamic viscosity calculation"
382 << exit(FatalError);
383 return volScalarField::null();
384 }
385}
386
387
388Foam::tmp<Foam::volScalarField> Foam::functionObjects::ITHACAforces::rho() const
389{
390 if (rhoName_ == "rhoInf")
391 {
392 return tmp<volScalarField>::New
393 (
394 IOobject
395 (
396 "rho",
397 mesh_.time().timeName(),
398 mesh_
399 ),
400 mesh_,
401 dimensionedScalar("rho", dimDensity, rhoRef_)
402 );
403 }
404
405 return (lookupObject<volScalarField>(rhoName_));
406}
407
408
409Foam::scalar Foam::functionObjects::ITHACAforces::rho(const volScalarField& p)
410const
411{
412 if (p.dimensions() == dimPressure)
413 {
414 return 1.0;
415 }
416
417 if (rhoName_ != "rhoInf")
418 {
419 FatalErrorInFunction
420 << "Dynamic pressure is expected but kinematic is provided."
421 << exit(FatalError);
422 }
423
424 return rhoRef_;
425}
426
427
428void Foam::functionObjects::ITHACAforces::applyBins
429(
430 const vectorField& Md,
431 const vectorField& fN,
432 const vectorField& fT,
433 const vectorField& fP,
434 const vectorField& d
435)
436{
437 if (nBin_ == 1)
438 {
439 force_[0][0] += sum(fN);
440 force_[1][0] += sum(fT);
441 force_[2][0] += sum(fP);
442 moment_[0][0] += sum(Md ^ fN);
443 moment_[1][0] += sum(Md ^ fT);
444 moment_[2][0] += sum(Md ^ fP);
445 }
446 else
447 {
448 scalarField dd((d & binDir_) - binMin_);
449 forAll(dd, i)
450 {
451 label bini = min(max(floor(dd[i] / binDx_), 0), force_[0].size() - 1);
452 force_[0][bini] += fN[i];
453 force_[1][bini] += fT[i];
454 force_[2][bini] += fP[i];
455 moment_[0][bini] += Md[i] ^ fN[i];
456 moment_[1][bini] += Md[i] ^ fT[i];
457 moment_[2][bini] += Md[i] ^ fP[i];
458 }
459 }
460}
461
462
463void Foam::functionObjects::ITHACAforces::addToFields
464(
465 const label patchi,
466 const vectorField& Md,
467 const vectorField& fN,
468 const vectorField& fT,
469 const vectorField& fP
470)
471{
472 if (!writeFields_)
473 {
474 return;
475 }
476
477 volVectorField& force =
478 lookupObjectRef<volVectorField>(fieldName("force"));
479 vectorField& pf = force.boundaryFieldRef()[patchi];
480 pf += fN + fT + fP;
481 volVectorField& moment =
482 lookupObjectRef<volVectorField>(fieldName("moment"));
483 vectorField& pm = moment.boundaryFieldRef()[patchi];
484 pm += Md;
485}
486
487
488void Foam::functionObjects::ITHACAforces::addToFields
489(
490 const labelList& cellIDs,
491 const vectorField& Md,
492 const vectorField& fN,
493 const vectorField& fT,
494 const vectorField& fP
495)
496{
497 if (!writeFields_)
498 {
499 return;
500 }
501
502 volVectorField& force =
503 lookupObjectRef<volVectorField>(fieldName("force"));
504 volVectorField& moment =
505 lookupObjectRef<volVectorField>(fieldName("moment"));
506 forAll(cellIDs, i)
507 {
508 label celli = cellIDs[i];
509 force[celli] += fN[i] + fT[i] + fP[i];
510 moment[celli] += Md[i];
511 }
512}
513
514
515void Foam::functionObjects::ITHACAforces::writeIntegratedForceMoment
516(
517 const string& descriptor,
518 const vectorField& fm0,
519 const vectorField& fm1,
520 const vectorField& fm2,
521 autoPtr<OFstream>& osPtr
522) const
523{
524 vector pressure = sum(fm0);
525 vector viscous = sum(fm1);
526 vector porous = sum(fm2);
527 vector total = pressure + viscous + porous;
528 Log << " Sum of " << descriptor.c_str() << nl
529 << " Total : " << total << nl
530 << " Pressure : " << pressure << nl
531 << " Viscous : " << viscous << nl;
532
533 if (porosity_)
534 {
535 Log << " Porous : " << porous << nl;
536 }
537
538 if (writeToFile())
539 {
540 Ostream& os = osPtr();
541#if OPENFOAM > 1906
542 writeCurrentTime(os);
543#else
544 writeTime(os);
545#endif
546 os << tab << total
547 << tab << pressure
548 << tab << viscous;
549
550 if (porosity_)
551 {
552 os << tab << porous;
553 }
554
555 os << endl;
556 }
557}
558
559
560void Foam::functionObjects::ITHACAforces::writeForces()
561{
562 Log << type() << " " << name() << " write:" << nl;
563 writeIntegratedForceMoment
564 (
565 "forces",
566 force_[0],
567 force_[1],
568 force_[2],
569 forceFilePtr_
570 );
571 writeIntegratedForceMoment
572 (
573 "moments",
574 moment_[0],
575 moment_[1],
576 moment_[2],
577 momentFilePtr_
578 );
579
580 if (localSystem_)
581 {
582 writeIntegratedForceMoment
583 (
584 "local forces",
585 coordSys_.localVector(force_[0]),
586 coordSys_.localVector(force_[1]),
587 coordSys_.localVector(force_[2]),
588 localForceFilePtr_
589 );
590 writeIntegratedForceMoment
591 (
592 "local moments",
593 coordSys_.localVector(moment_[0]),
594 coordSys_.localVector(moment_[1]),
595 coordSys_.localVector(moment_[2]),
596 localMomentFilePtr_
597 );
598 }
599
600 Log << endl;
601}
602
603
604void Foam::functionObjects::ITHACAforces::writeBinnedForceMoment
605(
606 const List<Field<vector >> & fm,
607 autoPtr<OFstream>& osPtr
608) const
609{
610 if ((nBin_ == 1) || !writeToFile())
611 {
612 return;
613 }
614
615 List<Field<vector >> f(fm);
616 if (binCumulative_)
617 {
618 for (label i = 1; i < f[0].size(); i++)
619 {
620 f[0][i] += f[0][i - 1];
621 f[1][i] += f[1][i - 1];
622 f[2][i] += f[2][i - 1];
623 }
624 }
625
626 Ostream& os = osPtr();
627#if OPENFOAM > 1906
628 writeCurrentTime(os);
629#else
630 writeTime(os);
631#endif
632 forAll(f[0], i)
633 {
634 vector total = f[0][i] + f[1][i] + f[2][i];
635 os << tab << total
636 << tab << f[0][i]
637 << tab << f[1][i];
638 if (porosity_)
639 {
640 os << tab << f[2][i];
641 }
642 }
643
644 os << nl;
645}
646
647
648void Foam::functionObjects::ITHACAforces::writeBins()
649{
650 writeBinnedForceMoment(force_, forceBinFilePtr_);
651 writeBinnedForceMoment(moment_, momentBinFilePtr_);
652 if (localSystem_)
653 {
654 List<Field<vector >> lf(3);
655 List<Field<vector >> lm(3);
656 lf[0] = coordSys_.localVector(force_[0]);
657 lf[1] = coordSys_.localVector(force_[1]);
658 lf[2] = coordSys_.localVector(force_[2]);
659 lm[0] = coordSys_.localVector(moment_[0]);
660 lm[1] = coordSys_.localVector(moment_[1]);
661 lm[2] = coordSys_.localVector(moment_[2]);
662 writeBinnedForceMoment(lf, localForceBinFilePtr_);
663 writeBinnedForceMoment(lm, localMomentBinFilePtr_);
664 }
665}
666
667
668// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
669
670Foam::functionObjects::ITHACAforces::ITHACAforces
671(
672 const word& name,
673 const Time& runTime,
674 const dictionary& dict,
675 bool readFields
676)
677 :
678 fvMeshFunctionObject(name, runTime, dict),
679 writeFile(mesh_, name),
680 force_(3),
681 moment_(3),
682 forceFilePtr_(),
683 momentFilePtr_(),
684 forceBinFilePtr_(),
685 momentBinFilePtr_(),
686 localForceFilePtr_(),
687 localMomentFilePtr_(),
688 localForceBinFilePtr_(),
689 localMomentBinFilePtr_(),
690 patchSet_(),
691 pName_(word::null),
692 UName_(word::null),
693 rhoName_(word::null),
694 directForceDensity_(false),
695 fDName_(""),
696 rhoRef_(VGREAT),
697 pRef_(0),
698 coordSys_(),
699 localSystem_(false),
700 porosity_(false),
701 nBin_(1),
702 binDir_(Zero),
703 binDx_(0.0),
704 binMin_(GREAT),
705 binPoints_(),
706 binCumulative_(true),
707 writeFields_(false),
708 initialised_(false)
709{
710 if (readFields)
711 {
712 read(dict);
713 Log << endl;
714 }
715}
716
717
718Foam::functionObjects::ITHACAforces::ITHACAforces
719(
720 const word& name,
721 const objectRegistry& obr,
722 const dictionary& dict,
723 bool readFields
724)
725 :
726 fvMeshFunctionObject(name, obr, dict),
727 writeFile(mesh_, name),
728 force_(3),
729 moment_(3),
730 forceFilePtr_(),
731 momentFilePtr_(),
732 forceBinFilePtr_(),
733 momentBinFilePtr_(),
734 localForceFilePtr_(),
735 localMomentFilePtr_(),
736 localForceBinFilePtr_(),
737 localMomentBinFilePtr_(),
738 patchSet_(),
739 pName_(word::null),
740 UName_(word::null),
741 rhoName_(word::null),
742 directForceDensity_(false),
743 fDName_(""),
744 rhoRef_(VGREAT),
745 pRef_(0),
746 coordSys_(),
747 localSystem_(false),
748 porosity_(false),
749 nBin_(1),
750 binDir_(Zero),
751 binDx_(0.0),
752 binMin_(GREAT),
753 binPoints_(),
754 binCumulative_(true),
755 writeFields_(false),
756 initialised_(false)
757{
758 if (readFields)
759 {
760 read(dict);
761 Log << endl;
762 }
763
764 /*
765 // Turn off writing to file
766 writeToFile_ = false;
767
768 forAll(force_, i)
769 {
770 force_[i].setSize(nBin_, vector::zero);
771 moment_[i].setSize(nBin_, vector::zero);
772 }
773 */
774}
775
776
777// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
778
779Foam::functionObjects::ITHACAforces::~ITHACAforces()
780{}
781
782
783// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
784
785bool Foam::functionObjects::ITHACAforces::read(const dictionary& dict)
786{
787 fvMeshFunctionObject::read(dict);
788 writeFile::read(dict);
789 initialised_ = false;
790 Info << type() << " " << name() << ":" << nl;
791 directForceDensity_ = dict.lookupOrDefault("directForceDensity", false);
792 const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
793 patchSet_ = pbm.patchSet(dict.get<wordRes>("patches"));
794
795 if (directForceDensity_)
796 {
797 // Optional entry for fDName
798 fDName_ = dict.lookupOrDefault<word>("fD", "fD");
799 }
800 else
801 {
802 // Optional entries U and p
803 pName_ = dict.lookupOrDefault<word>("p", "p");
804 UName_ = dict.lookupOrDefault<word>("U", "U");
805 rhoName_ = dict.lookupOrDefault<word>("rho", "rho");
806
807 // Reference density needed for incompressible calculations
808 if (rhoName_ == "rhoInf")
809 {
810 dict.readEntry("rhoInf", rhoRef_);
811 }
812
813 // Reference pressure, 0 by default
814 pRef_ = dict.lookupOrDefault<scalar>("pRef", 0.0);
815 }
816 coordSys_.clear();
817 localSystem_ = false;
818
819 // Centre of rotation for moment calculations
820 // specified directly, from coordinate system, or implicitly (0 0 0)
821 if (!dict.readIfPresent<point>("CofR", coordSys_.origin()))
822 {
823 // The 'coordinateSystem' sub-dictionary is optional,
824 // but enforce use of a cartesian system.
825 if (dict.found(coordinateSystem::typeName_()))
826 {
827 // New() for access to indirect (global) coordinate system
828 coordSys_ =
829 coordinateSystem::New
830 (
831 obr_, dict, coordinateSystem::typeName_()
832 );
833 }
834 else
835 {
836 coordSys_ = coordSystem::cartesian(dict);
837 }
838
839 localSystem_ = true;
840 }
841 dict.readIfPresent("porosity", porosity_);
842
843 if (porosity_)
844 {
845 Info << " Including porosity effects" << endl;
846 }
847 else
848 {
849 Info << " Not including porosity effects" << endl;
850 }
851 if (dict.found("binData"))
852 {
853 const dictionary& binDict(dict.subDict("binData"));
854 binDict.readEntry("nBin", nBin_);
855
856 if (nBin_ < 0)
857 {
858 FatalIOErrorInFunction(dict)
859 << "Number of bins (nBin) must be zero or greater"
860 << exit(FatalIOError);
861 }
862 else if (nBin_ == 0)
863 {
864 // Case of no bins equates to a single bin to collect all data
865 nBin_ = 1;
866 }
867 else
868 {
869 binDict.readEntry("cumulative", binCumulative_);
870 binDict.readEntry("direction", binDir_);
871 binDir_.normalise();
872 }
873 }
874 writeFields_ = dict.lookupOrDefault("writeFields", false);
875
876 if (writeFields_)
877 {
878 Info << " Fields will be written" << endl;
879 volVectorField* forcePtr
880 (
881 new volVectorField
882 (
883 IOobject
884 (
885 fieldName("force"),
886 time_.timeName(),
887 mesh_,
888 IOobject::NO_READ,
889 IOobject::NO_WRITE
890 ),
891 mesh_,
892 dimensionedVector(dimForce, Zero)
893 )
894 );
895 mesh_.objectRegistry::store(forcePtr);
896 volVectorField* momentPtr
897 (
898 new volVectorField
899 (
900 IOobject
901 (
902 fieldName("moment"),
903 time_.timeName(),
904 mesh_,
905 IOobject::NO_READ,
906 IOobject::NO_WRITE
907 ),
908 mesh_,
909 dimensionedVector(dimForce * dimLength, Zero)
910 )
911 );
912 mesh_.objectRegistry::store(momentPtr);
913 }
914 return true;
915}
916
917
918void Foam::functionObjects::ITHACAforces::calcForcesMoment()
919{
920 initialise();
921 resetFields();
922
923 if (directForceDensity_)
924 {
925 const volVectorField& fD = lookupObject<volVectorField>(fDName_);
926 const surfaceVectorField::Boundary& Sfb = mesh_.Sf().boundaryField();
927
928 for (const label patchi : patchSet_)
929 {
930 vectorField Md
931 (
932 mesh_.C().boundaryField()[patchi] - coordSys_.origin()
933 );
934 scalarField sA(mag(Sfb[patchi]));
935 // Normal force = surfaceUnitNormal*(surfaceNormal & forceDensity)
936 vectorField fN
937 (
938 Sfb[patchi] / sA
939 * (
940 Sfb[patchi] & fD.boundaryField()[patchi]
941 )
942 );
943 // Tangential force (total force minus normal fN)
944 vectorField fT(sA * fD.boundaryField()[patchi] - fN);
945 // Porous force
946 vectorField fP(Md.size(), Zero);
947 addToFields(patchi, Md, fN, fT, fP);
948 applyBins(Md, fN, fT, fP, mesh_.C().boundaryField()[patchi]);
949 }
950 }
951 else
952 {
953 const volScalarField& p = lookupObject<volScalarField>(pName_);
954 const surfaceVectorField::Boundary& Sfb = mesh_.Sf().boundaryField();
955 tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
956 const volSymmTensorField::Boundary& devRhoReffb
957 = tdevRhoReff().boundaryField();
958 // Scale pRef by density for incompressible simulations
959 scalar pRef = pRef_ / rho(p);
960
961 for (const label patchi : patchSet_)
962 {
963 vectorField Md
964 (
965 mesh_.C().boundaryField()[patchi] - coordSys_.origin()
966 );
967 vectorField fN
968 (
969 rho(p) * Sfb[patchi] * (p.boundaryField()[patchi] - pRef)
970 );
971 vectorField fT(Sfb[patchi] & devRhoReffb[patchi]);
972 vectorField fP(Md.size(), Zero);
973 addToFields(patchi, Md, fN, fT, fP);
974 applyBins(Md, fN, fT, fP, mesh_.C().boundaryField()[patchi]);
975 }
976 }
977
978 if (porosity_)
979 {
980 const volVectorField& U = lookupObject<volVectorField>(UName_);
981 const volScalarField rho(this->rho());
982 const volScalarField mu(this->mu());
983 const HashTable<const porosityModel*> models =
984 obr_.lookupClass<porosityModel>();
985
986 if (models.empty())
987 {
988 WarningInFunction
989 << "Porosity effects requested, but no porosity models found "
990 << "in the database"
991 << endl;
992 }
993
994 forAllConstIters(models, iter)
995 {
996 // Non-const access required if mesh is changing
997 porosityModel& pm = const_cast<porosityModel&>(* iter());
998 vectorField fPTot(pm.force(U, rho, mu));
999 const labelList& cellZoneIDs = pm.cellZoneIDs();
1000
1001 for (const label zonei : cellZoneIDs)
1002 {
1003 const cellZone& cZone = mesh_.cellZones()[zonei];
1004 const vectorField d(mesh_.C(), cZone);
1005 const vectorField fP(fPTot, cZone);
1006 const vectorField Md(d - coordSys_.origin());
1007 const vectorField fDummy(Md.size(), Zero);
1008 addToFields(cZone, Md, fDummy, fDummy, fP);
1009 applyBins(Md, fDummy, fDummy, fP, d);
1010 }
1011 }
1012 }
1013
1014 Pstream::listCombineGather(force_, plusEqOp<vectorField>());
1015 Pstream::listCombineGather(moment_, plusEqOp<vectorField>());
1016 Pstream::listCombineScatter(force_);
1017 Pstream::listCombineScatter(moment_);
1018}
1019
1020
1021Foam::vector Foam::functionObjects::ITHACAforces::forceEff() const
1022{
1023 return sum(force_[0]) + sum(force_[1]) + sum(force_[2]);
1024}
1025
1026
1027Foam::vector Foam::functionObjects::ITHACAforces::momentEff() const
1028{
1029 return sum(moment_[0]) + sum(moment_[1]) + sum(moment_[2]);
1030}
1031
1032
1033bool Foam::functionObjects::ITHACAforces::execute()
1034{
1035 calcForcesMoment();
1036
1037 if (Pstream::master())
1038 {
1039 createFiles();
1040 writeForces();
1041 writeBins();
1042 Log << endl;
1043 }
1044
1045 // Write state/results information
1046 setResult("normalForce", sum(force_[0]));
1047 setResult("tangentialForce", sum(force_[1]));
1048 setResult("porousForce", sum(force_[2]));
1049 setResult("normalMoment", sum(moment_[0]));
1050 setResult("tangentialMoment", sum(moment_[1]));
1051 setResult("porousMoment", sum(moment_[2]));
1052 return true;
1053}
1054
1055
1056bool Foam::functionObjects::ITHACAforces::write()
1057{
1058 if (writeFields_)
1059 {
1060 lookupObject<volVectorField>(fieldName("force")).write();
1061 lookupObject<volVectorField>(fieldName("moment")).write();
1062 }
1063
1064 return true;
1065}
1066
1067
1068Foam::vector Foam::functionObjects::ITHACAforces::forcePressure() const
1069{
1070 return sum(force_[0]);
1071}
1072
1073Foam::vector Foam::functionObjects::ITHACAforces::forceTau() const
1074{
1075 return sum(force_[1]);
1076}
1077
1078
1079
1080// ************************************************************************* //
1081
1082
1083#else
1084#include "ITHACAforces.H"
1085#include "fvcGrad.H"
1086#include "porosityModel.H"
1087#include "turbulentTransportModel.H"
1088#include "turbulentFluidThermoModel.H"
1089#include "addToRunTimeSelectionTable.H"
1090
1091// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
1092
1093namespace Foam
1094{
1101}
1102
1103
1104// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
1105
1106Foam::wordList Foam::functionObjects::ITHACAforces::createFileNames
1107(
1108 const dictionary& dict
1109) const
1110{
1111 DynamicList<word> names(1);
1112 const word forceType(dict.lookup("type"));
1113 // Name for file(MAIN_FILE=0)
1114 names.append(forceType);
1115
1116 if (dict.found("binData"))
1117 {
1118 const dictionary& binDict(dict.subDict("binData"));
1119 label nb = readLabel(binDict.lookup("nBin"));
1120
1121 if (nb > 0)
1122 {
1123 // Name for file(BINS_FILE=1)
1124 names.append(forceType + "_bins");
1125 }
1126 }
1127
1128 return names;
1129}
1130
1131
1132void Foam::functionObjects::ITHACAforces::writeFileHeader(const label i)
1133{
1134 switch (fileID(i))
1135 {
1136 case MAIN_FILE:
1137 {
1138 // force data
1139 writeHeader(file(i), "Forces");
1140 writeHeaderValue(file(i), "CofR", coordSys_.origin());
1141 writeCommented(file(i), "Time");
1142 const word forceTypes("(pressure viscous porous)");
1143 file(i)
1144 << "forces" << forceTypes << tab
1145 << "moments" << forceTypes;
1146
1147 if (localSystem_)
1148 {
1149 file(i)
1150 << tab
1151 << "localForces" << forceTypes << tab
1152 << "localMoments" << forceTypes;
1153 }
1154
1155 break;
1156 }
1157
1158 case BINS_FILE:
1159 {
1160 // bin data
1161 writeHeader(file(i), "Force bins");
1162 writeHeaderValue(file(i), "bins", nBin_);
1163 writeHeaderValue(file(i), "start", binMin_);
1164 writeHeaderValue(file(i), "delta", binDx_);
1165 writeHeaderValue(file(i), "direction", binDir_);
1166 vectorField binPoints(nBin_);
1167 writeCommented(file(i), "x co-ords :");
1168 forAll(binPoints, pointi)
1169 {
1170 binPoints[pointi] = (binMin_ + (pointi + 1) * binDx_) * binDir_;
1171 file(i) << tab << binPoints[pointi].x();
1172 }
1173
1174 file(i) << nl;
1175 writeCommented(file(i), "y co-ords :");
1176 forAll(binPoints, pointi)
1177 {
1178 file(i) << tab << binPoints[pointi].y();
1179 }
1180
1181 file(i) << nl;
1182 writeCommented(file(i), "z co-ords :");
1183 forAll(binPoints, pointi)
1184 {
1185 file(i) << tab << binPoints[pointi].z();
1186 }
1187
1188 file(i) << nl;
1189 writeCommented(file(i), "Time");
1190 const word binForceTypes("[pressure,viscous,porous]");
1191
1192 for (label j = 0; j < nBin_; j++)
1193 {
1194 const word jn('(' + Foam::name(j) + ')');
1195 const word f("forces" + jn + binForceTypes);
1196 const word m("moments" + jn + binForceTypes);
1197 file(i) << tab << f << tab << m;
1198 }
1199
1200 if (localSystem_)
1201 {
1202 for (label j = 0; j < nBin_; j++)
1203 {
1204 const word jn('(' + Foam::name(j) + ')');
1205 const word f("localForces" + jn + binForceTypes);
1206 const word m("localMoments" + jn + binForceTypes);
1207 file(i) << tab << f << tab << m;
1208 }
1209 }
1210
1211 break;
1212 }
1213
1214 default:
1215 {
1216 FatalErrorInFunction
1217 << "Unhandled file index: " << i
1218 << abort(FatalError);
1219 }
1220 }
1221
1222 file(i) << endl;
1223}
1224
1225
1226void Foam::functionObjects::ITHACAforces::initialise()
1227{
1228 if (initialised_)
1229 {
1230 return;
1231 }
1232
1233 if (directForceDensity_)
1234 {
1235 if (!obr_.foundObject<volVectorField>(fDName_))
1236 {
1237 FatalErrorInFunction
1238 << "Could not find " << fDName_ << " in database."
1239 << exit(FatalError);
1240 }
1241 }
1242 else
1243 {
1244 if
1245 (
1246 !obr_.foundObject<volVectorField>(UName_)
1247 || !obr_.foundObject<volScalarField>(pName_)
1248 )
1249 {
1250 FatalErrorInFunction
1251 << "Could not find " << UName_ << ", " << pName_
1252 << exit(FatalError);
1253 }
1254
1255 if
1256 (
1257 rhoName_ != "rhoInf"
1258 && !obr_.foundObject<volScalarField>(rhoName_)
1259 )
1260 {
1261 FatalErrorInFunction
1262 << "Could not find " << rhoName_
1263 << exit(FatalError);
1264 }
1265 }
1266
1267 initialised_ = true;
1268}
1269
1270
1271Foam::tmp<Foam::volSymmTensorField>
1272Foam::functionObjects::ITHACAforces::devRhoReff() const
1273{
1274 typedef compressible::turbulenceModel cmpTurbModel;
1275 typedef incompressible::turbulenceModel icoTurbModel;
1276
1277 if (obr_.foundObject<cmpTurbModel>(cmpTurbModel::propertiesName))
1278 {
1279 const cmpTurbModel& turb =
1280 obr_.lookupObject<cmpTurbModel>(cmpTurbModel::propertiesName);
1281 return turb.devRhoReff();
1282 }
1283 else if (obr_.foundObject<icoTurbModel>(icoTurbModel::propertiesName))
1284 {
1285 const incompressible::turbulenceModel& turb =
1286 obr_.lookupObject<icoTurbModel>(icoTurbModel::propertiesName);
1287 return rho() * turb.devReff();
1288 }
1289 else if (obr_.foundObject<fluidThermo>(fluidThermo::dictName))
1290 {
1291 const fluidThermo& thermo =
1292 obr_.lookupObject<fluidThermo>(fluidThermo::dictName);
1293 const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
1294 return -thermo.mu() * dev(twoSymm(fvc::grad(U)));
1295 }
1296 else if
1297 (
1298 obr_.foundObject<transportModel>("transportProperties")
1299 )
1300 {
1301 const transportModel& laminarT =
1302 obr_.lookupObject<transportModel>("transportProperties");
1303 const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
1304 return -rho() * laminarT.nu() * dev(twoSymm(fvc::grad(U)));
1305 }
1306 else if (obr_.foundObject<dictionary>("transportProperties"))
1307 {
1308 const dictionary& transportProperties =
1309 obr_.lookupObject<dictionary>("transportProperties");
1310 dimensionedScalar nu(transportProperties.lookup("nu"));
1311 const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
1312 return -rho() * nu * dev(twoSymm(fvc::grad(U)));
1313 }
1314 else
1315 {
1316 FatalErrorInFunction
1317 << "No valid model for viscous stress calculation"
1318 << exit(FatalError);
1319 return volSymmTensorField::null();
1320 }
1321}
1322
1323
1324Foam::tmp<Foam::volScalarField> Foam::functionObjects::ITHACAforces::mu() const
1325{
1326 if (obr_.foundObject<fluidThermo>(basicThermo::dictName))
1327 {
1328 const fluidThermo& thermo =
1329 obr_.lookupObject<fluidThermo>(basicThermo::dictName);
1330 return thermo.mu();
1331 }
1332 else if
1333 (
1334 obr_.foundObject<transportModel>("transportProperties")
1335 )
1336 {
1337 const transportModel& laminarT =
1338 obr_.lookupObject<transportModel>("transportProperties");
1339 return rho() * laminarT.nu();
1340 }
1341 else if (obr_.foundObject<dictionary>("transportProperties"))
1342 {
1343 const dictionary& transportProperties =
1344 obr_.lookupObject<dictionary>("transportProperties");
1345 dimensionedScalar nu
1346 (
1347 "nu",
1348 dimViscosity,
1349 transportProperties.lookup("nu")
1350 );
1351 return rho() * nu;
1352 }
1353 else
1354 {
1355 FatalErrorInFunction
1356 << "No valid model for dynamic viscosity calculation"
1357 << exit(FatalError);
1358 return volScalarField::null();
1359 }
1360}
1361
1362
1363Foam::tmp<Foam::volScalarField> Foam::functionObjects::ITHACAforces::rho() const
1364{
1365 if (rhoName_ == "rhoInf")
1366 {
1367 return tmp<volScalarField>
1368 (
1369 new volScalarField
1370 (
1371 IOobject
1372 (
1373 "rho",
1374 mesh_.time().timeName(),
1375 mesh_
1376 ),
1377 mesh_,
1378 dimensionedScalar("rho", dimDensity, rhoRef_)
1379 )
1380 );
1381 }
1382 else
1383 {
1384 return (obr_.lookupObject<volScalarField>(rhoName_));
1385 }
1386}
1387
1388
1389Foam::scalar Foam::functionObjects::ITHACAforces::rho(const volScalarField& p)
1390const
1391{
1392 if (p.dimensions() == dimPressure)
1393 {
1394 return 1.0;
1395 }
1396 else
1397 {
1398 if (rhoName_ != "rhoInf")
1399 {
1400 FatalErrorInFunction
1401 << "Dynamic pressure is expected but kinematic is provided."
1402 << exit(FatalError);
1403 }
1404
1405 return rhoRef_;
1406 }
1407}
1408
1409
1410void Foam::functionObjects::ITHACAforces::applyBins
1411(
1412 const vectorField& Md,
1413 const vectorField& fN,
1414 const vectorField& fT,
1415 const vectorField& fP,
1416 const vectorField& d
1417)
1418{
1419 if (nBin_ == 1)
1420 {
1421 force_[0][0] += sum(fN);
1422 force_[1][0] += sum(fT);
1423 force_[2][0] += sum(fP);
1424 moment_[0][0] += sum(Md ^ fN);
1425 moment_[1][0] += sum(Md ^ fT);
1426 moment_[2][0] += sum(Md ^ fP);
1427 }
1428 else
1429 {
1430 scalarField dd((d & binDir_) - binMin_);
1431 forAll(dd, i)
1432 {
1433 label bini = min(max(floor(dd[i] / binDx_), 0), force_[0].size() - 1);
1434 force_[0][bini] += fN[i];
1435 force_[1][bini] += fT[i];
1436 force_[2][bini] += fP[i];
1437 moment_[0][bini] += Md[i] ^ fN[i];
1438 moment_[1][bini] += Md[i] ^ fT[i];
1439 moment_[2][bini] += Md[i] ^ fP[i];
1440 }
1441 }
1442}
1443
1444
1445void Foam::functionObjects::ITHACAforces::writeForces()
1446{
1447 Log << type() << " " << name() << " write:" << nl
1448 << " sum of forces:" << nl
1449 << " pressure : " << sum(force_[0]) << nl
1450 << " viscous : " << sum(force_[1]) << nl
1451 << " porous : " << sum(force_[2]) << nl
1452 << " sum of moments:" << nl
1453 << " pressure : " << sum(moment_[0]) << nl
1454 << " viscous : " << sum(moment_[1]) << nl
1455 << " porous : " << sum(moment_[2])
1456 << endl;
1457 writeTime(file(MAIN_FILE));
1458 file(MAIN_FILE) << tab << setw(1) << '('
1459 << sum(force_[0]) << setw(1) << ' '
1460 << sum(force_[1]) << setw(1) << ' '
1461 << sum(force_[2]) << setw(3) << ") ("
1462 << sum(moment_[0]) << setw(1) << ' '
1463 << sum(moment_[1]) << setw(1) << ' '
1464 << sum(moment_[2]) << setw(1) << ')'
1465 << endl;
1466
1467 if (localSystem_)
1468 {
1469 vectorField localForceN(coordSys_.localVector(force_[0]));
1470 vectorField localForceT(coordSys_.localVector(force_[1]));
1471 vectorField localForceP(coordSys_.localVector(force_[2]));
1472 vectorField localMomentN(coordSys_.localVector(moment_[0]));
1473 vectorField localMomentT(coordSys_.localVector(moment_[1]));
1474 vectorField localMomentP(coordSys_.localVector(moment_[2]));
1475 writeTime(file(MAIN_FILE));
1476 file(MAIN_FILE) << tab << setw(1) << '('
1477 << sum(localForceN) << setw(1) << ' '
1478 << sum(localForceT) << setw(1) << ' '
1479 << sum(localForceP) << setw(3) << ") ("
1480 << sum(localMomentN) << setw(1) << ' '
1481 << sum(localMomentT) << setw(1) << ' '
1482 << sum(localMomentP) << setw(1) << ')'
1483 << endl;
1484 }
1485}
1486
1487
1488void Foam::functionObjects::ITHACAforces::writeBins()
1489{
1490 if (nBin_ == 1)
1491 {
1492 return;
1493 }
1494
1495 List<Field<vector >> f(force_);
1496 List<Field<vector >> m(moment_);
1497
1498 if (binCumulative_)
1499 {
1500 for (label i = 1; i < f[0].size(); i++)
1501 {
1502 f[0][i] += f[0][i - 1];
1503 f[1][i] += f[1][i - 1];
1504 f[2][i] += f[2][i - 1];
1505 m[0][i] += m[0][i - 1];
1506 m[1][i] += m[1][i - 1];
1507 m[2][i] += m[2][i - 1];
1508 }
1509 }
1510 writeTime(file(BINS_FILE));
1511 forAll(f[0], i)
1512 {
1513 file(BINS_FILE)
1514 << tab << setw(1) << '('
1515 << f[0][i] << setw(1) << ' '
1516 << f[1][i] << setw(1) << ' '
1517 << f[2][i] << setw(3) << ") ("
1518 << m[0][i] << setw(1) << ' '
1519 << m[1][i] << setw(1) << ' '
1520 << m[2][i] << setw(1) << ')';
1521 }
1522
1523 if (localSystem_)
1524 {
1525 List<Field<vector >> lf(3);
1526 List<Field<vector >> lm(3);
1527 lf[0] = coordSys_.localVector(force_[0]);
1528 lf[1] = coordSys_.localVector(force_[1]);
1529 lf[2] = coordSys_.localVector(force_[2]);
1530 lm[0] = coordSys_.localVector(moment_[0]);
1531 lm[1] = coordSys_.localVector(moment_[1]);
1532 lm[2] = coordSys_.localVector(moment_[2]);
1533
1534 if (binCumulative_)
1535 {
1536 for (label i = 1; i < lf[0].size(); i++)
1537 {
1538 lf[0][i] += lf[0][i - 1];
1539 lf[1][i] += lf[1][i - 1];
1540 lf[2][i] += lf[2][i - 1];
1541 lm[0][i] += lm[0][i - 1];
1542 lm[1][i] += lm[1][i - 1];
1543 lm[2][i] += lm[2][i - 1];
1544 }
1545 }
1546 forAll(lf[0], i)
1547 {
1548 file(BINS_FILE)
1549 << tab << setw(1) << '('
1550 << lf[0][i] << setw(1) << ' '
1551 << lf[1][i] << setw(1) << ' '
1552 << lf[2][i] << setw(3) << ") ("
1553 << lm[0][i] << setw(1) << ' '
1554 << lm[1][i] << setw(1) << ' '
1555 << lm[2][i] << setw(1) << ')';
1556 }
1557 }
1558 file(BINS_FILE) << endl;
1559}
1560
1561
1562// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1563
1564Foam::functionObjects::ITHACAforces::ITHACAforces
1565(
1566 const word& name,
1567 const Time& runTime,
1568 const dictionary& dict
1569)
1570 :
1571 fvMeshFunctionObject(name, runTime, dict),
1572 logFiles(obr_, name),
1573 force_(3),
1574 moment_(3),
1575 patchSet_(),
1576 pName_(word::null),
1577 UName_(word::null),
1578 rhoName_(word::null),
1579 directForceDensity_(false),
1580 fDName_(""),
1581 rhoRef_(VGREAT),
1582 pRef_(0),
1583 coordSys_(),
1584 localSystem_(false),
1585 porosity_(false),
1586 nBin_(1),
1587 binDir_(Zero),
1588 binDx_(0.0),
1589 binMin_(GREAT),
1590 binPoints_(),
1591 binCumulative_(true),
1592 initialised_(false)
1593{
1594 read(dict);
1595 resetNames(createFileNames(dict));
1596}
1597
1598
1599Foam::functionObjects::ITHACAforces::ITHACAforces
1600(
1601 const word& name,
1602 const objectRegistry& obr,
1603 const dictionary& dict
1604)
1605 :
1606 fvMeshFunctionObject(name, obr, dict),
1607 logFiles(obr_, name),
1608 force_(3),
1609 moment_(3),
1610 patchSet_(),
1611 pName_(word::null),
1612 UName_(word::null),
1613 rhoName_(word::null),
1614 directForceDensity_(false),
1615 fDName_(""),
1616 rhoRef_(VGREAT),
1617 pRef_(0),
1618 coordSys_(),
1619 localSystem_(false),
1620 porosity_(false),
1621 nBin_(1),
1622 binDir_(Zero),
1623 binDx_(0.0),
1624 binMin_(GREAT),
1625 binPoints_(),
1626 binCumulative_(true),
1627 initialised_(false)
1628{
1629 read(dict);
1630 resetNames(createFileNames(dict));
1631}
1632
1633
1634// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1635
1638
1639
1640// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1641
1642bool Foam::functionObjects::ITHACAforces::read(const dictionary& dict)
1643{
1644 fvMeshFunctionObject::read(dict);
1645 initialised_ = false;
1646 Log << type() << " " << name() << ":" << nl;
1647 directForceDensity_ = dict.lookupOrDefault("directForceDensity", false);
1648 const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
1649 patchSet_ = pbm.patchSet(wordReList(dict.lookup("patches")));
1650
1651 if (directForceDensity_)
1652 {
1653 // Optional entry for fDName
1654 fDName_ = dict.lookupOrDefault<word>("fD", "fD");
1655 }
1656 else
1657 {
1658 // Optional entries U and p
1659 pName_ = dict.lookupOrDefault<word>("p", "p");
1660 UName_ = dict.lookupOrDefault<word>("U", "U");
1661 rhoName_ = dict.lookupOrDefault<word>("rho", "rho");
1662
1663 // Reference density needed for incompressible calculations
1664 if (rhoName_ == "rhoInf")
1665 {
1666 dict.lookup("rhoInf") >> rhoRef_;
1667 }
1668
1669 // Reference pressure, 0 by default
1670 pRef_ = dict.lookupOrDefault<scalar>("pRef", 0.0);
1671 }
1672
1673 coordSys_.clear();
1674
1675 // Centre of rotation for moment calculations
1676 // specified directly, from coordinate system, or implicitly (0 0 0)
1677 if (!dict.readIfPresent<point>("CofR", coordSys_.origin()))
1678 {
1679 coordSys_ = coordinateSystem(obr_, dict);
1680 localSystem_ = true;
1681 }
1682
1683 dict.readIfPresent("porosity", porosity_);
1684
1685 if (porosity_)
1686 {
1687 Log << " Including porosity effects" << endl;
1688 }
1689 else
1690 {
1691 Log << " Not including porosity effects" << endl;
1692 }
1693
1694 if (dict.found("binData"))
1695 {
1696 const dictionary& binDict(dict.subDict("binData"));
1697 binDict.lookup("nBin") >> nBin_;
1698
1699 if (nBin_ < 0)
1700 {
1701 FatalIOErrorInFunction(dict)
1702 << "Number of bins (nBin) must be zero or greater"
1703 << exit(FatalIOError);
1704 }
1705 else if ((nBin_ == 0) || (nBin_ == 1))
1706 {
1707 nBin_ = 1;
1708 forAll(force_, i)
1709 {
1710 force_[i].setSize(1);
1711 moment_[i].setSize(1);
1712 }
1713 }
1714
1715 if (nBin_ > 1)
1716 {
1717 binDict.lookup("direction") >> binDir_;
1718 binDir_ /= mag(binDir_);
1719 binMin_ = GREAT;
1720 scalar binMax = -GREAT;
1721 forAllConstIter(labelHashSet, patchSet_, iter)
1722 {
1723 label patchi = iter.key();
1724 const polyPatch& pp = pbm[patchi];
1725 scalarField d(pp.faceCentres() & binDir_);
1726 binMin_ = min(min(d), binMin_);
1727 binMax = max(max(d), binMax);
1728 }
1729
1730 reduce(binMin_, minOp<scalar>());
1731 reduce(binMax, maxOp<scalar>());
1732 // slightly boost binMax so that region of interest is fully
1733 // within bounds
1734 binMax = 1.0001 * (binMax - binMin_) + binMin_;
1735 binDx_ = (binMax - binMin_) / scalar(nBin_);
1736 // create the bin points used for writing
1737 binPoints_.setSize(nBin_);
1738 forAll(binPoints_, i)
1739 {
1740 binPoints_[i] = (i + 0.5) * binDir_ * binDx_;
1741 }
1742
1743 binDict.lookup("cumulative") >> binCumulative_;
1744 // allocate storage for forces and moments
1745 forAll(force_, i)
1746 {
1747 force_[i].setSize(nBin_);
1748 moment_[i].setSize(nBin_);
1749 }
1750 }
1751 }
1752
1753 if (nBin_ == 1)
1754 {
1755 // allocate storage for forces and moments
1756 force_[0].setSize(1);
1757 force_[1].setSize(1);
1758 force_[2].setSize(1);
1759 moment_[0].setSize(1);
1760 moment_[1].setSize(1);
1761 moment_[2].setSize(1);
1762 }
1763
1764 return true;
1765}
1766
1767
1768void Foam::functionObjects::ITHACAforces::calcForcesMoment()
1769{
1770 initialise();
1771 force_[0] = Zero;
1772 force_[1] = Zero;
1773 force_[2] = Zero;
1774 moment_[0] = Zero;
1775 moment_[1] = Zero;
1776 moment_[2] = Zero;
1777
1778 if (directForceDensity_)
1779 {
1780 const volVectorField& fD = obr_.lookupObject<volVectorField>(fDName_);
1781 const surfaceVectorField::Boundary& Sfb =
1782 mesh_.Sf().boundaryField();
1783 forAllConstIter(labelHashSet, patchSet_, iter)
1784 {
1785 label patchi = iter.key();
1786 vectorField Md
1787 (
1788 mesh_.C().boundaryField()[patchi] - coordSys_.origin()
1789 );
1790 scalarField sA(mag(Sfb[patchi]));
1791 // Normal force = surfaceUnitNormal*(surfaceNormal & forceDensity)
1792 vectorField fN
1793 (
1794 Sfb[patchi] / sA
1795 * (
1796 Sfb[patchi] & fD.boundaryField()[patchi]
1797 )
1798 );
1799 // Tangential force (total force minus normal fN)
1800 vectorField fT(sA * fD.boundaryField()[patchi] - fN);
1801 //- Porous force
1802 vectorField fP(Md.size(), Zero);
1803 applyBins(Md, fN, fT, fP, mesh_.C().boundaryField()[patchi]);
1804 }
1805 }
1806 else
1807 {
1808 const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
1809 const surfaceVectorField::Boundary& Sfb =
1810 mesh_.Sf().boundaryField();
1811 tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
1812 const volSymmTensorField::Boundary& devRhoReffb
1813 = tdevRhoReff().boundaryField();
1814 // Scale pRef by density for incompressible simulations
1815 scalar pRef = pRef_ / rho(p);
1816 forAllConstIter(labelHashSet, patchSet_, iter)
1817 {
1818 label patchi = iter.key();
1819 vectorField Md
1820 (
1821 mesh_.C().boundaryField()[patchi] - coordSys_.origin()
1822 );
1823 vectorField fN
1824 (
1825 rho(p) * Sfb[patchi] * (p.boundaryField()[patchi] - pRef)
1826 );
1827 vectorField fT(Sfb[patchi] & devRhoReffb[patchi]);
1828 vectorField fP(Md.size(), Zero);
1829 applyBins(Md, fN, fT, fP, mesh_.C().boundaryField()[patchi]);
1830 }
1831 }
1832
1833 if (porosity_)
1834 {
1835 const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
1836 const volScalarField rho(this->rho());
1837 const volScalarField mu(this->mu());
1838 const HashTable<const porosityModel*> models =
1839 obr_.lookupClass<porosityModel>();
1840
1841 if (models.empty())
1842 {
1843 WarningInFunction
1844 << "Porosity effects requested, but no porosity models found "
1845 << "in the database"
1846 << endl;
1847 }
1848
1849 forAllConstIter(HashTable<const porosityModel*>, models, iter)
1850 {
1851 // non-const access required if mesh is changing
1852 porosityModel& pm = const_cast<porosityModel&>(* iter());
1853 vectorField fPTot(pm.force(U, rho, mu));
1854 const labelList& cellZoneIDs = pm.cellZoneIDs();
1855 forAll(cellZoneIDs, i)
1856 {
1857 label zoneI = cellZoneIDs[i];
1858 const cellZone& cZone = mesh_.cellZones()[zoneI];
1859 const vectorField d(mesh_.C(), cZone);
1860 const vectorField fP(fPTot, cZone);
1861 const vectorField Md(d - coordSys_.origin());
1862 const vectorField fDummy(Md.size(), Zero);
1863 applyBins(Md, fDummy, fDummy, fP, d);
1864 }
1865 }
1866 }
1867
1868 Pstream::listCombineGather(force_, plusEqOp<vectorField>());
1869 Pstream::listCombineGather(moment_, plusEqOp<vectorField>());
1870 Pstream::listCombineScatter(force_);
1871 Pstream::listCombineScatter(moment_);
1872}
1873
1874
1875Foam::vector Foam::functionObjects::ITHACAforces::forceEff() const
1876{
1877 return sum(force_[0]) + sum(force_[1]) + sum(force_[2]);
1878}
1879
1880
1881Foam::vector Foam::functionObjects::ITHACAforces::momentEff() const
1882{
1883 return sum(moment_[0]) + sum(moment_[1]) + sum(moment_[2]);
1884}
1885
1886
1887bool Foam::functionObjects::ITHACAforces::execute()
1888{
1889 return true;
1890}
1891
1892Foam::vector Foam::functionObjects::ITHACAforces::forcePorous() const
1893{
1894 return sum(force_[2]);
1895}
1896
1897Foam::vector Foam::functionObjects::ITHACAforces::forcePressure() const
1898{
1899 return sum(force_[0]);
1900}
1901
1902Foam::vector Foam::functionObjects::ITHACAforces::forceTau() const
1903{
1904 return sum(force_[1]);
1905}
1906
1907
1908bool Foam::functionObjects::ITHACAforces::write()
1909{
1910 calcForcesMoment();
1911
1912 if (Pstream::master())
1913 {
1914 logFiles::write();
1915 writeForces();
1916 writeBins();
1917 Log << endl;
1918 }
1919
1920 return true;
1921}
1922
1923#endif
1924// ************************************************************************* //
forAll
forAll(example_CG.gList, solutionI)
Definition CGtest.H:21
runTime
Foam::Time & runTime
Definition createTime.H:33
read
turbulence read()
Foam::functionObjects::ITHACAforces::moment_
List< Field< vector > > moment_
Definition ITHACAforces.H:67
Foam::functionObjects::ITHACAforces::forcePressure
virtual vector forcePressure() const
Definition ITHACAforces.C:1897
Foam::functionObjects::ITHACAforces::writeBinnedForceMoment
void writeBinnedForceMoment(const List< Field< vector > > &fm, autoPtr< OFstream > &osPtr) const
Foam::functionObjects::ITHACAforces::writeIntegratedForceMoment
void writeIntegratedForceMoment(const string &descriptor, const vectorField &fm0, const vectorField &fm1, const vectorField &fm2, autoPtr< OFstream > &osPtr) const
Foam::functionObjects::ITHACAforces::read
virtual bool read(const dictionary &)
Definition ITHACAforces.C:1642
Foam::functionObjects::ITHACAforces::writeFileHeader
virtual void writeFileHeader(const label i)
Definition ITHACAforces.C:1132
Foam::functionObjects::ITHACAforces::forcePorous
virtual vector forcePorous() const
Definition ITHACAforces.C:1892
Foam::functionObjects::ITHACAforces::rho
tmp< volScalarField > rho() const
Definition ITHACAforces.C:1363
Foam::functionObjects::ITHACAforces::mu
tmp< volScalarField > mu() const
Definition ITHACAforces.C:1324
Foam::functionObjects::ITHACAforces::applyBins
void applyBins(const vectorField &Md, const vectorField &fN, const vectorField &fT, const vectorField &fP, const vectorField &d)
Definition ITHACAforces.C:1411
Foam::functionObjects::ITHACAforces::createFileNames
wordList createFileNames(const dictionary &dict) const
Definition ITHACAforces.C:1107
Foam::functionObjects::ITHACAforces::addToFields
void addToFields(const label patchi, const vectorField &Md, const vectorField &fN, const vectorField &fT, const vectorField &fP)
Foam::functionObjects::ITHACAforces::fieldName
word fieldName(const word &name) const
Foam::functionObjects::ITHACAforces::writeBinHeader
void writeBinHeader(const word &header, Ostream &os) const
Foam::functionObjects::ITHACAforces::force_
List< Field< vector > > force_
Definition ITHACAforces.H:64
Foam::functionObjects::ITHACAforces::writeIntegratedHeader
void writeIntegratedHeader(const word &header, Ostream &os) const
Foam::functionObjects::ITHACAforces::devRhoReff
tmp< volSymmTensorField > devRhoReff() const
Definition ITHACAforces.C:1272
Foam::functionObjects::ITHACAforces::ITHACAforces
ITHACAforces(const ITHACAforces &)
nu
dimensionedScalar & nu
Definition readConstants.H:30
EigenFunctions::max
T max(Eigen::SparseMatrix< T > &mat, label &ind_row, label &ind_col)
Find the maximum of a sparse Matrix (Useful for DEIM)
Definition EigenFunctions.H:277
EigenFunctions::min
T min(Eigen::SparseMatrix< T > &mat, label &ind_row, label &ind_col)
Find the minimum of a sparse Matrix (Useful for DEIM)
Definition EigenFunctions.H:301
Foam::functionObjects::defineTypeNameAndDebug
defineTypeNameAndDebug(ITHACAforces, 0)
Foam::functionObjects::addToRunTimeSelectionTable
addToRunTimeSelectionTable(functionObject, ITHACAforces, dictionary)
Foam::names
const char * NamedEnum< Foam::incompressible::turbulentHeatFluxTemperatureFvPatchScalarField::heatSourceType, 2 >::names[]
Definition turbulentHeatFluxTemperatureFvPatchScalarField.C:41
U
volVectorField & U
Definition createFields.H:42
p
volScalarField & p
Definition createFields.H:42
rho
volScalarField & rho
Definition createFields.H:67
thermo
fluidThermo & thermo
Definition createFields.H:33
transportProperties
IOdictionary transportProperties(IOobject("transportProperties", runTime.constant(), mesh, IOobject::MUST_READ_IF_MODIFIED, IOobject::NO_WRITE))
i
label i
Definition pEqn.H:46
Generated by   doxygen 1.13.2