forces.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 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 "fvcGrad.H"
28 #include "porosityModel.H"
33 #include "fluidThermo.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40 namespace functionObjects
41 {
42  defineTypeNameAndDebug(forces, 0);
43 
44  addToRunTimeSelectionTable(functionObject, forces, dictionary);
45 }
46 }
47 
48 
49 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
50 
52 (
53  const dictionary& dict
54 ) const
55 {
56  DynamicList<word> names(1);
57 
58  const word forceType(dict.lookup("type"));
59 
60  // Name for file(fileID::mainFile=0)
61  names.append(forceType);
62 
63  if (dict.found("binData"))
64  {
65  const dictionary& binDict(dict.subDict("binData"));
66  const label nb = binDict.lookup<label>("nBin");
67  if (nb > 0)
68  {
69  // Name for file(fileID::binsFile=1)
70  names.append(forceType + "_bins");
71  }
72  }
73 
74  return move(names);
75 }
76 
77 
79 {
80  const word forceTypes
81  (
82  porosity_
83  ? "(pressure viscous porous)"
84  : "(pressure viscous)"
85  );
86 
87  switch (fileID(i))
88  {
89  case fileID::mainFile:
90  {
91  // force data
92 
93  writeHeader(file(i), "Forces");
94  writeHeaderValue(file(i), "CofR", coordSys_.origin());
95  writeCommented(file(i), "Time");
96 
97  file(i)
98  << "forces" << forceTypes << tab
99  << "moments" << forceTypes;
100 
101  if (localSystem_)
102  {
103  file(i)
104  << tab
105  << "localForces" << forceTypes << tab
106  << "localMoments" << forceTypes;
107  }
108 
109  break;
110  }
111  case fileID::binsFile:
112  {
113  // bin data
114 
115  writeHeader(file(i), "Force bins");
116  writeHeaderValue(file(i), "bins", nBin_);
117  writeHeaderValue(file(i), "start", binMin_);
118  writeHeaderValue(file(i), "delta", binDx_);
119  writeHeaderValue(file(i), "direction", binDir_);
120 
121  vectorField binPoints(nBin_);
122  writeCommented(file(i), "x co-ords :");
123  forAll(binPoints, pointi)
124  {
125  binPoints[pointi] = (binMin_ + (pointi + 1)*binDx_)*binDir_;
126  file(i) << tab << binPoints[pointi].x();
127  }
128  file(i) << nl;
129 
130  writeCommented(file(i), "y co-ords :");
131  forAll(binPoints, pointi)
132  {
133  file(i) << tab << binPoints[pointi].y();
134  }
135  file(i) << nl;
136 
137  writeCommented(file(i), "z co-ords :");
138  forAll(binPoints, pointi)
139  {
140  file(i) << tab << binPoints[pointi].z();
141  }
142  file(i) << nl;
143 
144  writeCommented(file(i), "Time");
145 
146  for (label j = 0; j < nBin_; j++)
147  {
148  const word jn('(' + Foam::name(j) + ')');
149  const word f("forces" + jn + forceTypes);
150  const word m("moments" + jn + forceTypes);
151 
152  file(i)<< tab << f << tab << m;
153  }
154  if (localSystem_)
155  {
156  for (label j = 0; j < nBin_; j++)
157  {
158  const word jn('(' + Foam::name(j) + ')');
159  const word f("localForces" + jn + forceTypes);
160  const word m("localMoments" + jn + forceTypes);
161 
162  file(i)<< tab << f << tab << m;
163  }
164  }
165 
166  break;
167  }
168  default:
169  {
171  << "Unhandled file index: " << i
172  << abort(FatalError);
173  }
174  }
175 
176  file(i)<< endl;
177 }
178 
179 
181 {
182  if (initialised_)
183  {
184  return;
185  }
186 
187  if (directForceDensity_)
188  {
189  if (!obr_.foundObject<volVectorField>(fDName_))
190  {
192  << "Could not find " << fDName_ << " in database."
193  << exit(FatalError);
194  }
195  }
196  else
197  {
198  if
199  (
200  !obr_.foundObject<volVectorField>(UName_)
201  || !obr_.foundObject<volScalarField>(pName_)
202 
203  )
204  {
206  << "Could not find " << UName_ << ", " << pName_
207  << exit(FatalError);
208  }
209 
210  if
211  (
212  rhoName_ != "rhoInf"
213  && !obr_.foundObject<volScalarField>(rhoName_)
214  )
215  {
217  << "Could not find " << rhoName_
218  << exit(FatalError);
219  }
220  }
221 
222  initialised_ = true;
223 }
224 
225 
228 {
230  typedef compressible::momentumTransportModel cmpModel;
231  typedef phaseIncompressible::momentumTransportModel phaseIcoModel;
232  typedef phaseCompressible::momentumTransportModel phaseCmpModel;
233 
234  const word& modelName = momentumTransportModel::typeName;
235  const word phaseModelName =
236  phaseName_ == word::null
237  ? word::null
238  : IOobject::groupName(momentumTransportModel::typeName, phaseName_);
239 
240  if (obr_.foundObject<icoModel>(modelName))
241  {
243  obr_.lookupObject<icoModel>(modelName);
244 
245  return alpha()*rho()*model.devSigma();
246  }
247  else if (obr_.foundObject<cmpModel>(modelName))
248  {
249  const cmpModel& model =
250  obr_.lookupObject<cmpModel>(modelName);
251 
252  return alpha()*model.devTau();
253  }
254  else if (obr_.foundObject<phaseIcoModel>(phaseModelName))
255  {
256  const phaseIcoModel& model =
257  obr_.lookupObject<phaseIcoModel>(phaseModelName);
258 
259  return rho()*model.devSigma();
260  }
261  else if (obr_.foundObject<phaseCmpModel>(phaseModelName))
262  {
263  const phaseCmpModel& model =
264  obr_.lookupObject<phaseCmpModel>(phaseModelName);
265 
266  return model.devTau();
267  }
268  else if (obr_.foundObject<dictionary>("transportProperties"))
269  {
270  // Legacy support for icoFoam
271 
272  const dictionary& transportProperties =
273  obr_.lookupObject<dictionary>("transportProperties");
274 
275  const dimensionedScalar nu
276  (
277  "nu",
278  dimViscosity,
279  transportProperties.lookup("nu")
280  );
281 
282  const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
283 
284  return -rho()*nu*dev(twoSymm(fvc::grad(U)));
285  }
286  else
287  {
289  << "No valid model for viscous stress calculation"
290  << exit(FatalError);
291 
292  return volSymmTensorField::null();
293  }
294 }
295 
296 
298 {
300  typedef compressible::momentumTransportModel cmpModel;
301  typedef phaseIncompressible::momentumTransportModel phaseIcoModel;
302  typedef phaseCompressible::momentumTransportModel phaseCmpModel;
303 
304  const word& modelName = momentumTransportModel::typeName;
305  const word phaseModelName =
306  phaseName_ == word::null
307  ? word::null
308  : IOobject::groupName(momentumTransportModel::typeName, phaseName_);
309 
310  if (obr_.foundObject<icoModel>(modelName))
311  {
313  obr_.lookupObject<icoModel>(modelName);
314 
315  return rho()*model.transport().nu();
316  }
317  else if (obr_.foundObject<cmpModel>(modelName))
318  {
319  const cmpModel& model =
320  obr_.lookupObject<cmpModel>(modelName);
321 
322  return model.transport().mu();
323  }
324  else if (obr_.foundObject<phaseIcoModel>(phaseModelName))
325  {
326  const phaseIcoModel& model =
327  obr_.lookupObject<phaseIcoModel>(phaseModelName);
328 
329  return rho()*model.transport().nu();
330  }
331  else if (obr_.foundObject<phaseCmpModel>(phaseModelName))
332  {
333  const phaseCmpModel& model =
334  obr_.lookupObject<phaseCmpModel>(phaseModelName);
335 
336  return model.transport().mu();
337  }
338  else if (obr_.foundObject<dictionary>("transportProperties"))
339  {
340  // Legacy support for icoFoam
341 
342  const dictionary& transportProperties =
343  obr_.lookupObject<dictionary>("transportProperties");
344 
345  const dimensionedScalar nu
346  (
347  "nu",
348  dimViscosity,
349  transportProperties.lookup("nu")
350  );
351 
352  return rho()*nu;
353  }
354  else
355  {
357  << "No valid model for dynamic viscosity calculation"
358  << exit(FatalError);
359 
360  return volScalarField::null();
361  }
362 }
363 
364 
366 {
367  if (rhoName_ == "rhoInf")
368  {
369  return volScalarField::New
370  (
371  "rho",
372  mesh_,
373  dimensionedScalar(dimDensity, rhoRef_)
374  );
375  }
376  else
377  {
378  return(obr_.lookupObject<volScalarField>(rhoName_));
379  }
380 }
381 
382 
384 {
385  if (p.dimensions() == dimPressure)
386  {
387  return 1.0;
388  }
389  else
390  {
391  if (rhoName_ != "rhoInf")
392  {
394  << "Dynamic pressure is expected but kinematic is provided."
395  << exit(FatalError);
396  }
397 
398  return rhoRef_;
399  }
400 }
401 
402 
404 {
405  if (phaseName_ == word::null)
406  {
407  return volScalarField::New
408  (
409  "alpha",
410  mesh_,
412  );
413  }
414  else
415  {
416  return obr_.lookupObject<volScalarField>
417  (
418  IOobject::groupName("alpha", phaseName_)
419  );
420  }
421 }
422 
423 
425 (
426  const label patchi
427 ) const
428 {
429  if (phaseName_ == word::null)
430  {
431  return tmp<scalarField>
432  (
433  new scalarField(mesh_.boundary()[patchi].size(), 1)
434  );
435  }
436  else
437  {
438  return obr_.lookupObject<volScalarField>
439  (
440  IOobject::groupName("alpha", phaseName_)
441  ).boundaryField()[patchi];
442  }
443 }
444 
445 
447 (
448  const vectorField& Md,
449  const vectorField& fN,
450  const vectorField& fT,
451  const vectorField& fP,
452  const vectorField& d
453 )
454 {
455  if (nBin_ == 1)
456  {
457  force_[0][0] += sum(fN);
458  force_[1][0] += sum(fT);
459  force_[2][0] += sum(fP);
460  moment_[0][0] += sum(Md^fN);
461  moment_[1][0] += sum(Md^fT);
462  moment_[2][0] += sum(Md^fP);
463  }
464  else
465  {
466  scalarField dd((d & binDir_) - binMin_);
467 
468  forAll(dd, i)
469  {
470  label bini = min(max(floor(dd[i]/binDx_), 0), force_[0].size() - 1);
471 
472  force_[0][bini] += fN[i];
473  force_[1][bini] += fT[i];
474  force_[2][bini] += fP[i];
475  moment_[0][bini] += Md[i]^fN[i];
476  moment_[1][bini] += Md[i]^fT[i];
477  moment_[2][bini] += Md[i]^fP[i];
478  }
479  }
480 }
481 
482 
484 {
485  Log << type() << " " << name() << " write:" << nl
486  << " sum of forces:" << nl
487  << " pressure : " << sum(force_[0]) << nl
488  << " viscous : " << sum(force_[1]) << nl
489  << " porous : " << sum(force_[2]) << nl
490  << " sum of moments:" << nl
491  << " pressure : " << sum(moment_[0]) << nl
492  << " viscous : " << sum(moment_[1]) << nl
493  << " porous : " << sum(moment_[2])
494  << endl;
495 
496  writeTime(file(fileID::mainFile));
497 
498  if (porosity_)
499  {
500  file(fileID::mainFile) << tab << setw(1) << '('
501  << sum(force_[0]) << setw(1) << ' '
502  << sum(force_[1]) << setw(1) << ' '
503  << sum(force_[2]) << setw(3) << ") ("
504  << sum(moment_[0]) << setw(1) << ' '
505  << sum(moment_[1]) << setw(1) << ' '
506  << sum(moment_[2]) << setw(1) << ')';
507  }
508  else
509  {
510  file(fileID::mainFile) << tab << setw(1) << '('
511  << sum(force_[0]) << setw(1) << ' '
512  << sum(force_[1]) << setw(3) << ") ("
513  << sum(moment_[0]) << setw(1) << ' '
514  << sum(moment_[1]) << setw(1) << ')';
515  }
516 
517  if (localSystem_)
518  {
519  if (porosity_)
520  {
521  file(fileID::mainFile) << tab << setw(1) << '('
522  << sum(coordSys_.localVector(force_[0])) << setw(1) << ' '
523  << sum(coordSys_.localVector(force_[1])) << setw(1) << ' '
524  << sum(coordSys_.localVector(force_[2])) << setw(3) << ") ("
525  << sum(coordSys_.localVector(moment_[0])) << setw(1) << ' '
526  << sum(coordSys_.localVector(moment_[1])) << setw(1) << ' '
527  << sum(coordSys_.localVector(moment_[2])) << setw(1) << ')';
528  }
529  else
530  {
531  file(fileID::mainFile) << tab << setw(1) << '('
532  << sum(coordSys_.localVector(force_[0])) << setw(1) << ' '
533  << sum(coordSys_.localVector(force_[1])) << setw(3) << ") ("
534  << sum(coordSys_.localVector(moment_[0])) << setw(1) << ' '
535  << sum(coordSys_.localVector(moment_[1])) << setw(1) << ')';
536  }
537  }
538 
539  file(fileID::mainFile) << endl;
540 }
541 
542 
544 {
545  if (nBin_ == 1)
546  {
547  return;
548  }
549 
550  List<Field<vector>> f(force_);
551  List<Field<vector>> m(moment_);
552 
553  if (binCumulative_)
554  {
555  for (label i = 1; i < f[0].size(); i++)
556  {
557  f[0][i] += f[0][i-1];
558  f[1][i] += f[1][i-1];
559  f[2][i] += f[2][i-1];
560 
561  m[0][i] += m[0][i-1];
562  m[1][i] += m[1][i-1];
563  m[2][i] += m[2][i-1];
564  }
565  }
566 
567  writeTime(file(fileID::binsFile));
568 
569 
570  forAll(f[0], i)
571  {
572  if (porosity_)
573  {
574  file(fileID::binsFile)
575  << tab << setw(1) << '('
576  << f[0][i] << setw(1) << ' '
577  << f[1][i] << setw(1) << ' '
578  << f[2][i] << setw(3) << ") ("
579  << m[0][i] << setw(1) << ' '
580  << m[1][i] << setw(1) << ' '
581  << m[2][i] << setw(1) << ')';
582  }
583  else
584  {
585  file(fileID::binsFile)
586  << tab << setw(1) << '('
587  << f[0][i] << setw(1) << ' '
588  << f[1][i] << setw(3) << ") ("
589  << m[0][i] << setw(1) << ' '
590  << m[1][i] << setw(1) << ')';
591  }
592  }
593 
594  if (localSystem_)
595  {
596  List<Field<vector>> lf(3);
597  List<Field<vector>> lm(3);
598  lf[0] = coordSys_.localVector(force_[0]);
599  lf[1] = coordSys_.localVector(force_[1]);
600  lf[2] = coordSys_.localVector(force_[2]);
601  lm[0] = coordSys_.localVector(moment_[0]);
602  lm[1] = coordSys_.localVector(moment_[1]);
603  lm[2] = coordSys_.localVector(moment_[2]);
604 
605  if (binCumulative_)
606  {
607  for (label i = 1; i < lf[0].size(); i++)
608  {
609  lf[0][i] += lf[0][i-1];
610  lf[1][i] += lf[1][i-1];
611  lf[2][i] += lf[2][i-1];
612  lm[0][i] += lm[0][i-1];
613  lm[1][i] += lm[1][i-1];
614  lm[2][i] += lm[2][i-1];
615  }
616  }
617 
618  forAll(lf[0], i)
619  {
620  if (porosity_)
621  {
622  file(fileID::binsFile)
623  << tab << setw(1) << '('
624  << lf[0][i] << setw(1) << ' '
625  << lf[1][i] << setw(1) << ' '
626  << lf[2][i] << setw(3) << ") ("
627  << lm[0][i] << setw(1) << ' '
628  << lm[1][i] << setw(1) << ' '
629  << lm[2][i] << setw(1) << ')';
630  }
631  else
632  {
633  file(fileID::binsFile)
634  << tab << setw(1) << '('
635  << lf[0][i] << setw(1) << ' '
636  << lf[1][i] << setw(3) << ") ("
637  << lm[0][i] << setw(1) << ' '
638  << lm[1][i] << setw(1) << ')';
639  }
640  }
641  }
642 
643  file(fileID::binsFile) << endl;
644 }
645 
646 
647 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
648 
650 (
651  const word& name,
652  const Time& runTime,
653  const dictionary& dict
654 )
655 :
656  fvMeshFunctionObject(name, runTime, dict),
657  logFiles(obr_, name),
658  force_(3),
659  moment_(3),
660  patchSet_(),
661  pName_(word::null),
662  UName_(word::null),
663  rhoName_(word::null),
664  phaseName_(word::null),
665  directForceDensity_(false),
666  fDName_(""),
667  rhoRef_(vGreat),
668  pRef_(0),
669  coordSys_("coordSys", vector::zero),
670  localSystem_(false),
671  porosity_(false),
672  nBin_(1),
673  binDir_(Zero),
674  binDx_(0.0),
675  binMin_(great),
676  binPoints_(),
677  binCumulative_(true),
678  initialised_(false)
679 {
680  read(dict);
681 }
682 
683 
685 (
686  const word& name,
687  const objectRegistry& obr,
688  const dictionary& dict
689 )
690 :
691  fvMeshFunctionObject(name, obr, dict),
692  logFiles(obr_, name),
693  force_(3),
694  moment_(3),
695  patchSet_(),
696  pName_(word::null),
697  UName_(word::null),
698  rhoName_(word::null),
699  phaseName_(word::null),
700  directForceDensity_(false),
701  fDName_(""),
702  rhoRef_(vGreat),
703  pRef_(0),
704  coordSys_("coordSys", vector::zero),
705  localSystem_(false),
706  porosity_(false),
707  nBin_(1),
708  binDir_(Zero),
709  binDx_(0.0),
710  binMin_(great),
711  binPoints_(),
712  binCumulative_(true),
713  initialised_(false)
714 {
715  read(dict);
716 }
717 
718 
719 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
720 
722 {}
723 
724 
725 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
726 
728 {
730 
731  initialised_ = false;
732 
733  Log << type() << " " << name() << ":" << nl;
734 
735  directForceDensity_ = dict.lookupOrDefault("directForceDensity", false);
736 
737  const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
738 
739  patchSet_ = pbm.patchSet(wordReList(dict.lookup("patches")));
740 
741  if (directForceDensity_)
742  {
743  // Optional entry for fDName
744  fDName_ = dict.lookupOrDefault<word>("fD", "fD");
745  }
746  else
747  {
748  // Optional phase entry
749  phaseName_ = dict.lookupOrDefault<word>("phase", word::null);
750 
751  // Optional U, p and rho entries
752  pName_ =
753  dict.lookupOrDefault<word>
754  (
755  "p",
756  IOobject::groupName("p", phaseName_)
757  );
758  UName_ =
759  dict.lookupOrDefault<word>
760  (
761  "U",
762  IOobject::groupName("U", phaseName_)
763  );
764  rhoName_ =
765  dict.lookupOrDefault<word>
766  (
767  "rho",
768  IOobject::groupName("rho", phaseName_)
769  );
770 
771  // Reference density needed for incompressible calculations
772  if (rhoName_ == "rhoInf")
773  {
774  dict.lookup("rhoInf") >> rhoRef_;
775  }
776 
777  // Reference pressure, 0 by default
778  pRef_ = dict.lookupOrDefault<scalar>("pRef", 0.0);
779  }
780 
781  // Centre of rotation for moment calculations
782  // specified directly, from coordinate system, or implicitly (0 0 0)
783  if (dict.found("CofR"))
784  {
785  coordSys_ = coordinateSystem("coordSys", vector(dict.lookup("CofR")));
786  localSystem_ = false;
787  }
788  else
789  {
790  coordSys_ = coordinateSystem("coordSys", dict);
791  localSystem_ = true;
792  }
793 
794  dict.readIfPresent("porosity", porosity_);
795  if (porosity_)
796  {
797  Log << " Including porosity effects" << endl;
798  }
799  else
800  {
801  Log << " Not including porosity effects" << endl;
802  }
803 
804  if (dict.found("binData"))
805  {
806  const dictionary& binDict(dict.subDict("binData"));
807  binDict.lookup("nBin") >> nBin_;
808 
809  if (nBin_ < 0)
810  {
812  << "Number of bins (nBin) must be zero or greater"
813  << exit(FatalIOError);
814  }
815  else if ((nBin_ == 0) || (nBin_ == 1))
816  {
817  nBin_ = 1;
818  forAll(force_, i)
819  {
820  force_[i].setSize(1);
821  moment_[i].setSize(1);
822  }
823  }
824 
825  if (nBin_ > 1)
826  {
827  binDict.lookup("direction") >> binDir_;
828  binDir_ /= mag(binDir_);
829 
830  binMin_ = great;
831  scalar binMax = -great;
832  forAllConstIter(labelHashSet, patchSet_, iter)
833  {
834  const label patchi = iter.key();
835  const polyPatch& pp = pbm[patchi];
836  const scalarField d(pp.faceCentres() & binDir_);
837  binMin_ = min(min(d), binMin_);
838  binMax = max(max(d), binMax);
839  }
840  reduce(binMin_, minOp<scalar>());
841  reduce(binMax, maxOp<scalar>());
842 
843  // slightly boost binMax so that region of interest is fully
844  // within bounds
845  binMax = 1.0001*(binMax - binMin_) + binMin_;
846 
847  binDx_ = (binMax - binMin_)/scalar(nBin_);
848 
849  // create the bin points used for writing
850  binPoints_.setSize(nBin_);
851  forAll(binPoints_, i)
852  {
853  binPoints_[i] = (i + 0.5)*binDir_*binDx_;
854  }
855 
856  binDict.lookup("cumulative") >> binCumulative_;
857 
858  // allocate storage for forces and moments
859  forAll(force_, i)
860  {
861  force_[i].setSize(nBin_);
862  moment_[i].setSize(nBin_);
863  }
864  }
865  }
866 
867  if (nBin_ == 1)
868  {
869  // allocate storage for forces and moments
870  force_[0].setSize(1);
871  force_[1].setSize(1);
872  force_[2].setSize(1);
873  moment_[0].setSize(1);
874  moment_[1].setSize(1);
875  moment_[2].setSize(1);
876  }
877 
878  resetNames(createFileNames(dict));
879 
880  return true;
881 }
882 
883 
885 {
886  initialise();
887 
888  force_[0] = Zero;
889  force_[1] = Zero;
890  force_[2] = Zero;
891 
892  moment_[0] = Zero;
893  moment_[1] = Zero;
894  moment_[2] = Zero;
895 
896  if (directForceDensity_)
897  {
898  const volVectorField& fD = obr_.lookupObject<volVectorField>(fDName_);
899 
900  const surfaceVectorField::Boundary& Sfb =
901  mesh_.Sf().boundaryField();
902 
903  forAllConstIter(labelHashSet, patchSet_, iter)
904  {
905  const label patchi = iter.key();
906 
907  const vectorField Md
908  (
909  mesh_.C().boundaryField()[patchi] - coordSys_.origin()
910  );
911 
912  const scalarField sA(mag(Sfb[patchi]));
913 
914  // Normal force = surfaceUnitNormal*(surfaceNormal & forceDensity)
915  const vectorField fN
916  (
917  Sfb[patchi]/sA
918  *(
919  Sfb[patchi] & fD.boundaryField()[patchi]
920  )
921  );
922 
923  // Tangential force (total force minus normal fN)
924  const vectorField fT(sA*fD.boundaryField()[patchi] - fN);
925 
926  //- Porous force
927  const vectorField fP(Md.size(), Zero);
928 
929  applyBins(Md, fN, fT, fP, mesh_.C().boundaryField()[patchi]);
930  }
931  }
932  else
933  {
934  const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);
935 
936  const surfaceVectorField::Boundary& Sfb =
937  mesh_.Sf().boundaryField();
938 
939  tmp<volSymmTensorField> tdevTau = devTau();
940  const volSymmTensorField::Boundary& devTaub =
941  tdevTau().boundaryField();
942 
943  // Scale pRef by density for incompressible simulations
944  const scalar pRef = pRef_/rho(p);
945 
946  forAllConstIter(labelHashSet, patchSet_, iter)
947  {
948  const label patchi = iter.key();
949 
950  const vectorField Md
951  (
952  mesh_.C().boundaryField()[patchi] - coordSys_.origin()
953  );
954 
955  const vectorField fN
956  (
957  alpha(patchi)
958  *rho(p)
959  *Sfb[patchi]
960  *(p.boundaryField()[patchi] - pRef)
961  );
962 
963  const vectorField fT(Sfb[patchi] & devTaub[patchi]);
964 
965  const vectorField fP(Md.size(), Zero);
966 
967  applyBins(Md, fN, fT, fP, mesh_.C().boundaryField()[patchi]);
968  }
969  }
970 
971  if (porosity_)
972  {
973  const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
974  const volScalarField rho(this->rho());
975  const volScalarField mu(this->mu());
976 
977  const HashTable<const porosityModel*> models =
978  obr_.lookupClass<porosityModel>();
979 
980  if (models.empty())
981  {
983  << "Porosity effects requested, but no porosity models found "
984  << "in the database"
985  << endl;
986  }
987 
989  {
990  const porosityModel& pm = *iter();
991 
992  const vectorField fPTot(pm.force(U, rho, mu));
993 
994  const labelList& cellZoneIDs = pm.cellZoneIDs();
995 
996  forAll(cellZoneIDs, i)
997  {
998  const label zoneI = cellZoneIDs[i];
999  const cellZone& cZone = mesh_.cellZones()[zoneI];
1000 
1001  const vectorField d(mesh_.C(), cZone);
1002  const vectorField fP(fPTot, cZone);
1003  const vectorField Md(d - coordSys_.origin());
1004 
1005  const vectorField fDummy(Md.size(), Zero);
1006 
1007  applyBins(Md, fDummy, fDummy, fP, d);
1008  }
1009  }
1010  }
1011 
1015  Pstream::listCombineScatter(moment_);
1016 }
1017 
1018 
1020 {
1021  return sum(force_[0]) + sum(force_[1]) + sum(force_[2]);
1022 }
1023 
1024 
1026 {
1027  return sum(moment_[0]) + sum(moment_[1]) + sum(moment_[2]);
1028 }
1029 
1030 
1032 {
1033  return true;
1034 }
1035 
1036 
1038 {
1039  calcForcesMoment();
1040 
1041  if (Pstream::master())
1042  {
1043  logFiles::write();
1044 
1045  writeForces();
1046 
1047  writeBins();
1048 
1049  Log << endl;
1050  }
1051 
1052  return true;
1053 }
1054 
1055 
1056 // ************************************************************************* //
virtual vector momentEff() const
Return the total moment.
Definition: forces.C:1025
Base class for other coordinate system specifications.
Templated abstract base class for multiphase compressible turbulence models.
const dimensionSet dimViscosity
virtual void writeFileHeader(const label i)
Output file header information.
Definition: forces.C:78
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:643
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)
Loop across all elements in list.
Definition: UList.H:434
const transportModel & transport() const
Access function to incompressible transport model.
virtual bool write()
Write function.
Definition: logFiles.C:167
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static const char tab
Definition: Ostream.H:259
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
Templated abstract base class for single-phase compressible turbulence models.
forces(const word &name, const Time &runTime, const dictionary &dict)
Construct from Time and dictionary.
Definition: forces.C:650
const Boundary & boundaryField() const
Return const-reference to the boundary field.
const dimensionSet dimPressure
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
addToRunTimeSelectionTable(functionObject, Qdot, dictionary)
bool empty() const
Return true if the hash table is empty.
Definition: HashTableI.H:72
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
virtual vector forceEff() const
Return the total force.
Definition: forces.C:1019
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const dimensionSet dimless
const labelList & cellZoneIDs() const
Return const access to the cell zone IDs.
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.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
virtual bool read(const dictionary &)
Read the forces data.
Definition: forces.C:727
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
Macros for easy insertion into run-time selection tables.
void writeBins()
Helper function to write bin data.
Definition: forces.C:543
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Templated abstract base class for multiphase incompressible turbulence models.
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:982
virtual tmp< vectorField > force(const volVectorField &U, const volScalarField &rho, const volScalarField &mu) const
Return the force over the cell zone(s)
virtual bool execute()
Execute, currently does nothing.
Definition: forces.C:1031
const dimensionSet & dimensions() const
Return dimensions.
virtual bool read(const dictionary &)
Read optional controls.
bool read(const char *, int32_t &)
Definition: int32IO.C:85
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
void writeHeader(std::ostream &, const bool isBinary, const std::string &title)
Write header.
Calculate the gradient of the given field.
A class for handling words, derived from string.
Definition: word.H:59
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
tmp< volScalarField > mu() const
Dynamic viscosity field.
Definition: forces.C:297
static word groupName(Name name, const word &group)
dimensionedScalar jn(const int n, const dimensionedScalar &ds)
const dimensionSet dimDensity
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
static const word null
An empty word.
Definition: word.H:77
virtual void calcForcesMoment()
Calculate the forces and moments.
Definition: forces.C:884
static const zero Zero
Definition: zero.H:97
tmp< volScalarField > rho() const
Return rho if specified otherwise rhoRef.
Definition: forces.C:365
void applyBins(const vectorField &Md, const vectorField &fN, const vectorField &fT, const vectorField &fP, const vectorField &d)
Accumulate bin data.
Definition: forces.C:447
An STL-conforming hash table.
Definition: HashTable.H:61
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Foam::polyBoundaryMesh.
static const char nl
Definition: Ostream.H:260
const dimensionedScalar mu
Atomic mass unit.
void initialise()
Initialise the fields.
Definition: forces.C:180
scalar pRef
Definition: createFields.H:6
labelList f(nPoints)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
A subset of mesh cells.
Definition: cellZone.H:61
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
defineTypeNameAndDebug(Qdot, 0)
static const GeometricField< symmTensor, fvPatchField, volMesh > & null()
Return a null geometric field.
void writeForces()
Helper function to write force data.
Definition: forces.C:483
Templated abstract base class for single-phase incompressible turbulence models.
label patchi
U
Definition: pEqn.H:72
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
#define WarningInFunction
Report a warning using Foam::Warning.
tmp< volSymmTensorField > devTau() const
Return the effective viscous stress (laminar + turbulent).
Definition: forces.C:227
#define FatalIOErrorInFunction(ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:335
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
tmp< volScalarField > alpha() const
Get the volume fraction field.
Definition: forces.C:403
#define Log
Report write to Foam::Info if the local log switch is true.
dimensioned< scalar > mag(const dimensioned< Type > &)
Specialisation of Foam::functionObject for an Foam::fvMesh, providing a reference to the Foam::fvMesh...
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
List< wordRe > wordReList
A List of wordRe (word or regular expression)
Definition: wordReList.H:50
volScalarField & p
A class for managing temporary objects.
Definition: PtrList.H:53
Registry of regIOobjects.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
virtual ~forces()
Destructor.
Definition: forces.C:721
fileID
Enumeration for ensuring the right file is accessed.
Definition: forces.H:215
virtual tmp< volSymmTensorField > devSigma() const
Return the effective stress tensor.
Top level model for porosity models.
Definition: porosityModel.H:54
virtual bool write()
Write the forces.
Definition: forces.C:1037
wordList createFileNames(const dictionary &dict) const
Create file names for forces and bins.
Definition: forces.C:52
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844
functionObject base class for creating, maintaining and writing log files e.g. integrated of averaged...
Definition: logFiles.H:57
IOerror FatalIOError
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:297