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