forces.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "forces.H"
27 #include "volFields.H"
28 #include "dictionary.H"
29 #include "Time.H"
30 #include "wordReList.H"
31 #include "fvcGrad.H"
32 #include "porosityModel.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40  defineTypeNameAndDebug(forces, 0);
41 }
42 
43 
44 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
45 
47 {
48  DynamicList<word> names(1);
49 
50  const word forceType(dict.lookup("type"));
51 
52  if (dict.found("binData"))
53  {
54  const dictionary& binDict(dict.subDict("binData"));
55  label nb = readLabel(binDict.lookup("nBin"));
56  if (nb > 0)
57  {
58  names.append(forceType + "_bins");
59  }
60  }
61 
62  names.append(forceType);
63 
64  return names;
65 }
66 
67 
69 {
70  if (i == 0)
71  {
72  // force data
73 
74  writeHeader(file(i), "Forces");
75  writeHeaderValue(file(i), "CofR", coordSys_.origin());
76  writeCommented(file(i), "Time");
77 
78  file(i)
79  << "forces(pressure viscous porous) "
80  << "moment(pressure viscous porous)";
81 
82  if (localSystem_)
83  {
84  file(i)
85  << tab
86  << "localForces(pressure,viscous,porous) "
87  << "localMoments(pressure,viscous,porous)";
88  }
89  }
90  else if (i == 1)
91  {
92  // bin data
93 
94  writeHeader(file(i), "Force bins");
95  writeHeaderValue(file(i), "bins", nBin_);
96  writeHeaderValue(file(i), "start", binMin_);
97  writeHeaderValue(file(i), "delta", binDx_);
98  writeHeaderValue(file(i), "direction", binDir_);
99 
100  vectorField binPoints(nBin_);
101  writeCommented(file(i), "x co-ords :");
102  forAll(binPoints, pointI)
103  {
104  binPoints[pointI] = (binMin_ + (pointI + 1)*binDx_)*binDir_;
105  file(i) << tab << binPoints[pointI].x();
106  }
107  file(i) << nl;
108 
109  writeCommented(file(i), "y co-ords :");
110  forAll(binPoints, pointI)
111  {
112  file(i) << tab << binPoints[pointI].y();
113  }
114  file(i) << nl;
115 
116  writeCommented(file(i), "z co-ords :");
117  forAll(binPoints, pointI)
118  {
119  file(i) << tab << binPoints[pointI].z();
120  }
121  file(i) << nl;
122 
123  writeCommented(file(i), "Time");
124 
125  for (label j = 0; j < nBin_; j++)
126  {
127  const word jn('(' + Foam::name(j) + ')');
128  const word f("forces" + jn + "[pressure,viscous,porous]");
129  const word m("moments" + jn + "[pressure,viscous,porous]");
130 
131  file(i)<< tab << f << tab << m;
132  }
133  if (localSystem_)
134  {
135  for (label j = 0; j < nBin_; j++)
136  {
137  const word jn('(' + Foam::name(j) + ')');
138  const word f("localForces" + jn + "[pressure,viscous,porous]");
139  const word m("localMoments" + jn + "[pressure,viscous,porous]");
140 
141  file(i)<< tab << f << tab << m;
142  }
143  }
144  }
145  else
146  {
147  FatalErrorIn("void Foam::forces::writeFileHeader(const label)")
148  << "Unhandled file index: " << i
149  << abort(FatalError);
150  }
151 
152  file(i)<< endl;
153 }
154 
155 
157 {
158  if (initialised_ || !active_)
159  {
160  return;
161  }
162 
163  if (directForceDensity_)
164  {
165  if (!obr_.foundObject<volVectorField>(fDName_))
166  {
167  active_ = false;
168  WarningIn("void Foam::forces::initialise()")
169  << "Could not find " << fDName_ << " in database." << nl
170  << " De-activating forces."
171  << endl;
172  }
173  }
174  else
175  {
176  if
177  (
178  !obr_.foundObject<volVectorField>(UName_)
179  || !obr_.foundObject<volScalarField>(pName_)
180  || (
181  rhoName_ != "rhoInf"
182  && !obr_.foundObject<volScalarField>(rhoName_)
183  )
184  )
185  {
186  active_ = false;
187 
188  WarningIn("void Foam::forces::initialise()")
189  << "Could not find " << UName_ << ", " << pName_;
190 
191  if (rhoName_ != "rhoInf")
192  {
193  Info<< " or " << rhoName_;
194  }
195 
196  Info<< " in database." << nl
197  << " De-activating forces." << endl;
198  }
199  }
200 
201  initialised_ = true;
202 }
203 
204 
206 {
207  typedef compressible::turbulenceModel cmpTurbModel;
208  typedef incompressible::turbulenceModel icoTurbModel;
209 
210  if (obr_.foundObject<cmpTurbModel>(cmpTurbModel::propertiesName))
211  {
212  const cmpTurbModel& turb =
213  obr_.lookupObject<cmpTurbModel>(cmpTurbModel::propertiesName);
214 
215  return turb.devRhoReff();
216  }
217  else if (obr_.foundObject<icoTurbModel>(icoTurbModel::propertiesName))
218  {
220  obr_.lookupObject<icoTurbModel>(icoTurbModel::propertiesName);
221 
222  return rho()*turb.devReff();
223  }
224  else if (obr_.foundObject<fluidThermo>(fluidThermo::typeName))
225  {
226  const fluidThermo& thermo =
227  obr_.lookupObject<fluidThermo>(fluidThermo::typeName);
228 
229  const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
230 
231  return -thermo.mu()*dev(twoSymm(fvc::grad(U)));
232  }
233  else if
234  (
235  obr_.foundObject<transportModel>("transportProperties")
236  )
237  {
238  const transportModel& laminarT =
239  obr_.lookupObject<transportModel>("transportProperties");
240 
241  const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
242 
243  return -rho()*laminarT.nu()*dev(twoSymm(fvc::grad(U)));
244  }
245  else if (obr_.foundObject<dictionary>("transportProperties"))
246  {
247  const dictionary& transportProperties =
248  obr_.lookupObject<dictionary>("transportProperties");
249 
250  dimensionedScalar nu(transportProperties.lookup("nu"));
251 
252  const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
253 
254  return -rho()*nu*dev(twoSymm(fvc::grad(U)));
255  }
256  else
257  {
258  FatalErrorIn("forces::devRhoReff()")
259  << "No valid model for viscous stress calculation"
260  << exit(FatalError);
261 
262  return volSymmTensorField::null();
263  }
264 }
265 
266 
268 {
269  if (obr_.foundObject<fluidThermo>(basicThermo::dictName))
270  {
271  const fluidThermo& thermo =
272  obr_.lookupObject<fluidThermo>(basicThermo::dictName);
273 
274  return thermo.mu();
275  }
276  else if
277  (
278  obr_.foundObject<transportModel>("transportProperties")
279  )
280  {
281  const transportModel& laminarT =
282  obr_.lookupObject<transportModel>("transportProperties");
283 
284  return rho()*laminarT.nu();
285  }
286  else if (obr_.foundObject<dictionary>("transportProperties"))
287  {
288  const dictionary& transportProperties =
289  obr_.lookupObject<dictionary>("transportProperties");
290 
291  dimensionedScalar nu(transportProperties.lookup("nu"));
292 
293  return rho()*nu;
294  }
295  else
296  {
297  FatalErrorIn("forces::mu()")
298  << "No valid model for dynamic viscosity calculation"
299  << exit(FatalError);
300 
301  return volScalarField::null();
302  }
303 }
304 
305 
307 {
308  if (rhoName_ == "rhoInf")
309  {
310  const fvMesh& mesh = refCast<const fvMesh>(obr_);
311 
312  return tmp<volScalarField>
313  (
314  new volScalarField
315  (
316  IOobject
317  (
318  "rho",
319  mesh.time().timeName(),
320  mesh
321  ),
322  mesh,
323  dimensionedScalar("rho", dimDensity, rhoRef_)
324  )
325  );
326  }
327  else
328  {
329  return(obr_.lookupObject<volScalarField>(rhoName_));
330  }
331 }
332 
333 
334 Foam::scalar Foam::forces::rho(const volScalarField& p) const
335 {
336  if (p.dimensions() == dimPressure)
337  {
338  return 1.0;
339  }
340  else
341  {
342  if (rhoName_ != "rhoInf")
343  {
344  FatalErrorIn("forces::rho(const volScalarField& p)")
345  << "Dynamic pressure is expected but kinematic is provided."
346  << exit(FatalError);
347  }
348 
349  return rhoRef_;
350  }
351 }
352 
353 
355 (
356  const vectorField& Md,
357  const vectorField& fN,
358  const vectorField& fT,
359  const vectorField& fP,
360  const vectorField& d
361 )
362 {
363  if (nBin_ == 1)
364  {
365  force_[0][0] += sum(fN);
366  force_[1][0] += sum(fT);
367  force_[2][0] += sum(fP);
368  moment_[0][0] += sum(Md^fN);
369  moment_[1][0] += sum(Md^fT);
370  moment_[2][0] += sum(Md^fP);
371  }
372  else
373  {
374  scalarField dd((d & binDir_) - binMin_);
375 
376  forAll(dd, i)
377  {
378  label bini = min(max(floor(dd[i]/binDx_), 0), force_[0].size() - 1);
379 
380  force_[0][bini] += fN[i];
381  force_[1][bini] += fT[i];
382  force_[2][bini] += fP[i];
383  moment_[0][bini] += Md[i]^fN[i];
384  moment_[1][bini] += Md[i]^fT[i];
385  moment_[2][bini] += Md[i]^fP[i];
386  }
387  }
388 }
389 
390 
392 {
393  if (log_) Info
394  << type() << " " << name_ << " output:" << nl
395  << " sum of forces:" << nl
396  << " pressure : " << sum(force_[0]) << nl
397  << " viscous : " << sum(force_[1]) << nl
398  << " porous : " << sum(force_[2]) << nl
399  << " sum of moments:" << nl
400  << " pressure : " << sum(moment_[0]) << nl
401  << " viscous : " << sum(moment_[1]) << nl
402  << " porous : " << sum(moment_[2])
403  << endl;
404 
405  file(0) << obr_.time().value() << tab << setw(1) << '('
406  << sum(force_[0]) << setw(1) << ' '
407  << sum(force_[1]) << setw(1) << ' '
408  << sum(force_[2]) << setw(3) << ") ("
409  << sum(moment_[0]) << setw(1) << ' '
410  << sum(moment_[1]) << setw(1) << ' '
411  << sum(moment_[2]) << setw(1) << ')'
412  << endl;
413 
414  if (localSystem_)
415  {
416  vectorField localForceN(coordSys_.localVector(force_[0]));
417  vectorField localForceT(coordSys_.localVector(force_[1]));
418  vectorField localForceP(coordSys_.localVector(force_[2]));
419  vectorField localMomentN(coordSys_.localVector(moment_[0]));
420  vectorField localMomentT(coordSys_.localVector(moment_[1]));
421  vectorField localMomentP(coordSys_.localVector(moment_[2]));
422 
423  file(0) << obr_.time().value() << tab << setw(1) << '('
424  << sum(localForceN) << setw(1) << ' '
425  << sum(localForceT) << setw(1) << ' '
426  << sum(localForceP) << setw(3) << ") ("
427  << sum(localMomentN) << setw(1) << ' '
428  << sum(localMomentT) << setw(1) << ' '
429  << sum(localMomentP) << setw(1) << ')'
430  << endl;
431  }
432 }
433 
434 
436 {
437  if (nBin_ == 1)
438  {
439  return;
440  }
441 
442  List<Field<vector> > f(force_);
443  List<Field<vector> > m(moment_);
444 
445  if (binCumulative_)
446  {
447  for (label i = 1; i < f[0].size(); i++)
448  {
449  f[0][i] += f[0][i-1];
450  f[1][i] += f[1][i-1];
451  f[2][i] += f[2][i-1];
452 
453  m[0][i] += m[0][i-1];
454  m[1][i] += m[1][i-1];
455  m[2][i] += m[2][i-1];
456  }
457  }
458 
459  file(1) << obr_.time().value();
460 
461  forAll(f[0], i)
462  {
463  file(1)
464  << tab << setw(1) << '('
465  << f[0][i] << setw(1) << ' '
466  << f[1][i] << setw(1) << ' '
467  << f[2][i] << setw(3) << ") ("
468  << m[0][i] << setw(1) << ' '
469  << m[1][i] << setw(1) << ' '
470  << m[2][i] << setw(1) << ')';
471  }
472 
473  if (localSystem_)
474  {
475  List<Field<vector> > lf(3);
476  List<Field<vector> > lm(3);
477  lf[0] = coordSys_.localVector(force_[0]);
478  lf[1] = coordSys_.localVector(force_[1]);
479  lf[2] = coordSys_.localVector(force_[2]);
480  lm[0] = coordSys_.localVector(moment_[0]);
481  lm[1] = coordSys_.localVector(moment_[1]);
482  lm[2] = coordSys_.localVector(moment_[2]);
483 
484  if (binCumulative_)
485  {
486  for (label i = 1; i < lf[0].size(); i++)
487  {
488  lf[0][i] += lf[0][i-1];
489  lf[1][i] += lf[1][i-1];
490  lf[2][i] += lf[2][i-1];
491  lm[0][i] += lm[0][i-1];
492  lm[1][i] += lm[1][i-1];
493  lm[2][i] += lm[2][i-1];
494  }
495  }
496 
497  forAll(lf[0], i)
498  {
499  file(1)
500  << tab << setw(1) << '('
501  << lf[0][i] << setw(1) << ' '
502  << lf[1][i] << setw(1) << ' '
503  << lf[2][i] << setw(3) << ") ("
504  << lm[0][i] << setw(1) << ' '
505  << lm[1][i] << setw(1) << ' '
506  << lm[2][i] << setw(1) << ')';
507  }
508  }
509 
510  file(1) << endl;
511 }
512 
513 
514 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
515 
517 (
518  const word& name,
519  const objectRegistry& obr,
520  const dictionary& dict,
521  const bool loadFromFiles,
522  const bool readFields
523 )
524 :
525  functionObjectFile(obr, name, createFileNames(dict)),
526  name_(name),
527  obr_(obr),
528  active_(true),
529  log_(true),
530  force_(3),
531  moment_(3),
532  patchSet_(),
533  pName_(word::null),
534  UName_(word::null),
535  rhoName_(word::null),
536  directForceDensity_(false),
537  fDName_(""),
538  rhoRef_(VGREAT),
539  pRef_(0),
540  coordSys_(),
541  localSystem_(false),
542  porosity_(false),
543  nBin_(1),
544  binDir_(vector::zero),
545  binDx_(0.0),
546  binMin_(GREAT),
547  binPoints_(),
548  binCumulative_(true),
549  initialised_(false)
550 {
551  // Check if the available mesh is an fvMesh otherise deactivate
552  if (isA<fvMesh>(obr_))
553  {
554  if (readFields)
555  {
556  read(dict);
557  Info<< endl;
558  }
559  }
560  else
561  {
562  active_ = false;
563  WarningIn
564  (
565  "Foam::forces::forces"
566  "("
567  "const word&, "
568  "const objectRegistry&, "
569  "const dictionary&, "
570  "const bool"
571  ")"
572  ) << "No fvMesh available, deactivating " << name_
573  << endl;
574  }
575 
576 }
577 
578 
580 (
581  const word& name,
582  const objectRegistry& obr,
583  const labelHashSet& patchSet,
584  const word& pName,
585  const word& UName,
586  const word& rhoName,
587  const scalar rhoInf,
588  const scalar pRef,
589  const coordinateSystem& coordSys
590 )
591 :
592  functionObjectFile(obr, name, typeName),
593  name_(name),
594  obr_(obr),
595  active_(true),
596  log_(true),
597  force_(3),
598  moment_(3),
599  patchSet_(patchSet),
600  pName_(pName),
601  UName_(UName),
602  rhoName_(rhoName),
603  directForceDensity_(false),
604  fDName_(""),
605  rhoRef_(rhoInf),
606  pRef_(pRef),
607  coordSys_(coordSys),
608  localSystem_(false),
609  porosity_(false),
610  nBin_(1),
611  binDir_(vector::zero),
612  binDx_(0.0),
613  binMin_(GREAT),
614  binPoints_(),
615  binCumulative_(true),
616  initialised_(false)
617 {
618  forAll(force_, i)
619  {
620  force_[i].setSize(nBin_);
621  moment_[i].setSize(nBin_);
622  }
623 }
624 
625 
626 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
627 
629 {}
630 
631 
632 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
633 
635 {
636  if (active_)
637  {
638  initialised_ = false;
639 
640  log_ = dict.lookupOrDefault<Switch>("log", false);
641 
642  if (log_) Info<< type() << " " << name_ << ":" << nl;
643 
644  directForceDensity_ = dict.lookupOrDefault("directForceDensity", false);
645 
646  const fvMesh& mesh = refCast<const fvMesh>(obr_);
647  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
648 
649  patchSet_ = pbm.patchSet(wordReList(dict.lookup("patches")));
650 
651  if (directForceDensity_)
652  {
653  // Optional entry for fDName
654  fDName_ = dict.lookupOrDefault<word>("fDName", "fD");
655  }
656  else
657  {
658  // Optional entries U and p
659  pName_ = dict.lookupOrDefault<word>("pName", "p");
660  UName_ = dict.lookupOrDefault<word>("UName", "U");
661  rhoName_ = dict.lookupOrDefault<word>("rhoName", "rho");
662 
663  // Reference density needed for incompressible calculations
664  rhoRef_ = readScalar(dict.lookup("rhoInf"));
665 
666  // Reference pressure, 0 by default
667  pRef_ = dict.lookupOrDefault<scalar>("pRef", 0.0);
668  }
669 
670  coordSys_.clear();
671 
672  // Centre of rotation for moment calculations
673  // specified directly, from coordinate system, or implicitly (0 0 0)
674  if (!dict.readIfPresent<point>("CofR", coordSys_.origin()))
675  {
676  coordSys_ = coordinateSystem(obr_, dict);
677  localSystem_ = true;
678  }
679 
680  dict.readIfPresent("porosity", porosity_);
681  if (porosity_)
682  {
683  if (log_) Info<< " Including porosity effects" << endl;
684  }
685  else
686  {
687  if (log_) Info<< " Not including porosity effects" << endl;
688  }
689 
690  if (dict.found("binData"))
691  {
692  const dictionary& binDict(dict.subDict("binData"));
693  binDict.lookup("nBin") >> nBin_;
694 
695  if (nBin_ < 0)
696  {
698  (
699  "void Foam::forces::read(const dictionary&)", dict
700  ) << "Number of bins (nBin) must be zero or greater"
701  << exit(FatalIOError);
702  }
703  else if ((nBin_ == 0) || (nBin_ == 1))
704  {
705  nBin_ = 1;
706  forAll(force_, i)
707  {
708  force_[i].setSize(1);
709  moment_[i].setSize(1);
710  }
711  }
712 
713  if (nBin_ > 1)
714  {
715  binDict.lookup("direction") >> binDir_;
716  binDir_ /= mag(binDir_);
717 
718  binMin_ = GREAT;
719  scalar binMax = -GREAT;
720  forAllConstIter(labelHashSet, patchSet_, iter)
721  {
722  label patchI = iter.key();
723  const polyPatch& pp = pbm[patchI];
724  scalarField d(pp.faceCentres() & binDir_);
725  binMin_ = min(min(d), binMin_);
726  binMax = max(max(d), binMax);
727  }
728  reduce(binMin_, minOp<scalar>());
729  reduce(binMax, maxOp<scalar>());
730 
731  // slightly boost binMax so that region of interest is fully
732  // within bounds
733  binMax = 1.0001*(binMax - binMin_) + binMin_;
734 
735  binDx_ = (binMax - binMin_)/scalar(nBin_);
736 
737  // create the bin points used for writing
738  binPoints_.setSize(nBin_);
739  forAll(binPoints_, i)
740  {
741  binPoints_[i] = (i + 0.5)*binDir_*binDx_;
742  }
743 
744  binDict.lookup("cumulative") >> binCumulative_;
745 
746  // allocate storage for forces and moments
747  forAll(force_, i)
748  {
749  force_[i].setSize(nBin_);
750  moment_[i].setSize(nBin_);
751  }
752  }
753  }
754 
755  if (nBin_ == 1)
756  {
757  // allocate storage for forces and moments
758  force_[0].setSize(1);
759  force_[1].setSize(1);
760  force_[2].setSize(1);
761  moment_[0].setSize(1);
762  moment_[1].setSize(1);
763  moment_[2].setSize(1);
764  }
765  }
766 }
767 
768 
770 {
771  // Do nothing - only valid on write
772 }
773 
774 
776 {
777  // Do nothing - only valid on write
778 }
779 
780 
782 {
783  // Do nothing - only valid on write
784 }
785 
786 
788 {
789  calcForcesMoment();
790 
791  if (!active_)
792  {
793  return;
794  }
795 
796  if (Pstream::master())
797  {
799 
800  writeForces();
801 
802  writeBins();
803 
804  if (log_) Info<< endl;
805  }
806 }
807 
808 
810 {
811  initialise();
812 
813  if (!active_)
814  {
815  return;
816  }
817 
818  force_[0] = vector::zero;
819  force_[1] = vector::zero;
820  force_[2] = vector::zero;
821 
822  moment_[0] = vector::zero;
823  moment_[1] = vector::zero;
824  moment_[2] = vector::zero;
825 
826  if (directForceDensity_)
827  {
828  const volVectorField& fD = obr_.lookupObject<volVectorField>(fDName_);
829 
830  const fvMesh& mesh = fD.mesh();
831 
832  const surfaceVectorField::GeometricBoundaryField& Sfb =
833  mesh.Sf().boundaryField();
834 
835  forAllConstIter(labelHashSet, patchSet_, iter)
836  {
837  label patchI = iter.key();
838 
839  vectorField Md
840  (
841  mesh.C().boundaryField()[patchI] - coordSys_.origin()
842  );
843 
844  scalarField sA(mag(Sfb[patchI]));
845 
846  // Normal force = surfaceUnitNormal*(surfaceNormal & forceDensity)
847  vectorField fN
848  (
849  Sfb[patchI]/sA
850  *(
851  Sfb[patchI] & fD.boundaryField()[patchI]
852  )
853  );
854 
855  // Tangential force (total force minus normal fN)
856  vectorField fT(sA*fD.boundaryField()[patchI] - fN);
857 
858  //- Porous force
859  vectorField fP(Md.size(), vector::zero);
860 
861  applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]);
862  }
863  }
864  else
865  {
866  const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
867  const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
868 
869  const fvMesh& mesh = U.mesh();
870 
871  const surfaceVectorField::GeometricBoundaryField& Sfb =
872  mesh.Sf().boundaryField();
873 
874  tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
875  const volSymmTensorField::GeometricBoundaryField& devRhoReffb
876  = tdevRhoReff().boundaryField();
877 
878  // Scale pRef by density for incompressible simulations
879  scalar pRef = pRef_/rho(p);
880 
881  forAllConstIter(labelHashSet, patchSet_, iter)
882  {
883  label patchI = iter.key();
884 
885  vectorField Md
886  (
887  mesh.C().boundaryField()[patchI] - coordSys_.origin()
888  );
889 
890  vectorField fN
891  (
892  rho(p)*Sfb[patchI]*(p.boundaryField()[patchI] - pRef)
893  );
894 
895  vectorField fT(Sfb[patchI] & devRhoReffb[patchI]);
896 
897  vectorField fP(Md.size(), vector::zero);
898 
899  applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchI]);
900  }
901  }
902 
903  if (porosity_)
904  {
905  const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
906  const volScalarField rho(this->rho());
907  const volScalarField mu(this->mu());
908 
909  const fvMesh& mesh = U.mesh();
910 
911  const HashTable<const porosityModel*> models =
912  obr_.lookupClass<porosityModel>();
913 
914  if (models.empty())
915  {
916  WarningIn("void Foam::forces::calcForcesMoment()")
917  << "Porosity effects requested, but no porosity models found "
918  << "in the database"
919  << endl;
920  }
921 
923  {
924  // non-const access required if mesh is changing
925  porosityModel& pm = const_cast<porosityModel&>(*iter());
926 
927  vectorField fPTot(pm.force(U, rho, mu));
928 
929  const labelList& cellZoneIDs = pm.cellZoneIDs();
930 
931  forAll(cellZoneIDs, i)
932  {
933  label zoneI = cellZoneIDs[i];
934  const cellZone& cZone = mesh.cellZones()[zoneI];
935 
936  const vectorField d(mesh.C(), cZone);
937  const vectorField fP(fPTot, cZone);
938  const vectorField Md(d - coordSys_.origin());
939 
940  const vectorField fDummy(Md.size(), vector::zero);
941 
942  applyBins(Md, fDummy, fDummy, fP, d);
943  }
944  }
945  }
946 
951 }
952 
953 
955 {
956  return sum(force_[0]) + sum(force_[1]) + sum(force_[2]);
957 }
958 
959 
961 {
962  return sum(moment_[0]) + sum(moment_[1]) + sum(moment_[2]);
963 }
964 
965 
966 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Definition: polyMesh.H:469
const dimensionSet dimPressure
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:306
virtual void timeSet()
Called when time was set at the end of the Time::operator++.
Definition: forces.C:781
const dimensionedScalar mu
Atomic mass unit.
Definition: createFields.H:13
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
void writeBins()
Helper function to write bin data.
Definition: forces.C:435
virtual void calcForcesMoment()
Calculate the forces and moments.
Definition: forces.C:809
dimensioned< scalar > mag(const dimensioned< Type > &)
GeometricBoundaryField & boundaryField()
Return reference to GeometricBoundaryField.
A subset of mesh cells.
Definition: cellZone.H:61
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
labelList f(nPoints)
tmp< volScalarField > rho() const
Return rho if rhoName is specified otherwise rhoRef.
Definition: forces.C:306
tmp< volSymmTensorField > devRhoReff() const
Return the effective viscous stress (laminar + turbulent).
Definition: forces.C:205
Base class for other coordinate system specifications.
dimensionedScalar jn(const int n, const dimensionedScalar &ds)
Templated abstract base class for single-phase incompressible turbulence models.
An STL-conforming hash table.
Definition: HashTable.H:61
List< wordRe > wordReList
A List of wordRe (word or regular expression)
Definition: wordReList.H:50
Top level model for porosity models.
Definition: porosityModel.H:55
A class for handling words, derived from string.
Definition: word.H:59
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:741
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
Definition: Switch.H:60
virtual void write()
Write the forces.
Definition: forces.C:787
messageStream Info
virtual tmp< volScalarField > mu() const =0
Dynamic viscosity of mixture [kg/m/s].
const Mesh & mesh() const
Return mesh.
dynamicFvMesh & mesh
const dimensionSet dimDensity
forces(const forces &)
Disallow default bitwise copy construct.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
virtual tmp< vectorField > force(const volVectorField &U, const volScalarField &rho, const volScalarField &mu)
Return the force over the cell zone(s)
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:638
Namespace for OpenFOAM.
compressible::turbulenceModel & turb
label readLabel(Istream &is)
Definition: label.H:64
void applyBins(const vectorField &Md, const vectorField &fN, const vectorField &fT, const vectorField &fP, const vectorField &d)
Accumulate bin data.
Definition: forces.C:355
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Calculate the gradient of the given field.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:284
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
tmp< volScalarField > mu() const
Dynamic viscosity field.
Definition: forces.C:267
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOerror FatalIOError
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
#define WarningIn(functionName)
Report a warning using Foam::Warning.
const volVectorField & C() const
Return cell centres as volVectorField.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
Foam::polyBoundaryMesh.
virtual void end()
Execute at the final time-loop, currently does nothing.
Definition: forces.C:775
volScalarField & p
Definition: createFields.H:51
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
#define forAll(list, i)
Definition: UList.H:421
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
const dimensionSet & dimensions() const
Return dimensions.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Fundamental fluid thermodynamic properties.
Definition: fluidThermo.H:49
fileName::Type type(const fileName &)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:589
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
virtual void write()
Write function.
labelHashSet patchSet(const UList< wordRe > &patchNames, const bool warnNotFound=true, const bool usePatchGroups=true) const
Return the set of patch IDs corresponding to the given names.
virtual vector momentEff() const
Return the total moment.
Definition: forces.C:960
static const GeometricField< symmTensor, fvPatchField, volMesh > & null()
Return a null geometric field.
static const char tab
Definition: Ostream.H:259
virtual tmp< volSymmTensorField > devReff() const
Return the effective stress tensor.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
void writeForces()
Helper function to write force data.
Definition: forces.C:391
error FatalError
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
virtual void read(const dictionary &)
Read the forces data.
Definition: forces.C:634
Registry of regIOobjects.
static const Vector zero
Definition: Vector.H:80
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
virtual void execute()
Execute, currently does nothing.
Definition: forces.C:769
virtual void writeFileHeader(const label i)
Output file header information.
Definition: forces.C:68
Base class for output file data handling.
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:398
bool read(const char *, int32_t &)
Definition: int32IO.C:87
volScalarField & nu
Base-class for all transport models used by the incompressible turbulence models. ...
void initialise()
Initialise the fields.
Definition: forces.C:156
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
virtual ~forces()
Destructor.
Definition: forces.C:628
const word dictName() const
Return the local dictionary name (final part of scoped name)
Definition: dictionary.H:115
scalar pRef
Definition: createFields.H:19
bool empty() const
Return true if the hash table is empty.
Definition: HashTableI.H:72
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
#define FatalIOErrorIn(functionName, ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:325
virtual vector forceEff() const
Return the total force.
Definition: forces.C:954
Templated wrapper class to provide compressible turbulence models thermal diffusivity based thermal t...
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
wordList createFileNames(const dictionary &dict) const
Create file names for forces and bins.
Definition: forces.C:46
U
Definition: pEqn.H:82
static const word null
An empty word.
Definition: word.H:77
defineTypeNameAndDebug(combustionModel, 0)
virtual tmp< volScalarField > nu() const =0
Return the laminar viscosity.
const labelList & cellZoneIDs() const
Return const access to the cell zone IDs.