powerLawLopesdaCosta.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) 2018-2023 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 
27 #include "powerLawLopesdaCosta.H"
28 #include "geometricOneField.H"
29 #include "fvMatrices.H"
30 #include "triSurfaceMesh.H"
31 #include "triSurfaceTools.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  namespace porosityModels
38  {
41  }
42 }
43 
44 
45 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
46 
48 (
49  const word& name,
50  const word& modelType,
51  const fvMesh& mesh,
52  const dictionary& dict
53 )
54 :
55  zoneName_(name + ":porousZone")
56 {
57  dictionary coeffs(dict.optionalSubDict(modelType + "Coeffs"));
58 
59  // Vertical direction within porous region
60  vector zDir(coeffs.lookup("zDir"));
61 
62  // Span of the search for the cells in the porous region
63  scalar searchSpan(coeffs.lookup<scalar>("searchSpan"));
64 
65  // Top surface file name defining the extent of the porous region
66  word topSurfaceFileName(coeffs.lookup("topSurface"));
67 
68  // List of the ground patches defining the lower surface
69  // of the porous region
70  List<wordRe> groundPatches(coeffs.lookup("groundPatches"));
71 
72  // Functional form of the porosity surface area per unit volume as a
73  // function of the normalised vertical position
75  (
77  (
78  dict.lookupEntryBackwardsCompatible
79  (
80  {"Av", "Sigma"},
81  false,
82  true
83  ).keyword(),
84  dict
85  )
86  );
87 
88  // Searchable triSurface for the top of the porous region
89  triSurfaceMesh searchSurf
90  (
91  IOobject
92  (
93  topSurfaceFileName,
94  mesh.time().constant(),
96  mesh.time()
97  )
98  );
99 
100  // Limit the porous cell search to those within the lateral and vertical
101  // extent of the top surface
102 
103  boundBox surfBounds(searchSurf.points());
104  boundBox searchBounds
105  (
106  surfBounds.min() - searchSpan*zDir, surfBounds.max()
107  );
108 
109  const pointField& C = mesh.cellCentres();
110 
111  // List of cells within the porous region
112  labelList porousCells(C.size());
113  label porousCelli = 0;
114 
115  forAll(C, celli)
116  {
117  if (searchBounds.contains(C[celli]))
118  {
119  porousCells[porousCelli++] = celli;
120  }
121  }
122 
123  porousCells.setSize(porousCelli);
124 
125  pointField start(porousCelli);
126  pointField end(porousCelli);
127 
128  forAll(porousCells, porousCelli)
129  {
130  start[porousCelli] = C[porousCells[porousCelli]];
131  end[porousCelli] = start[porousCelli] + searchSpan*zDir;
132  }
133 
134  // Field of distance between the cell centre and the corresponding point
135  // on the porous region top surface
136  scalarField zTop(porousCelli);
137 
138  // Search the porous cells for those with a corresponding point on the
139  // porous region top surface
140  List<pointIndexHit> hitInfo;
141  searchSurf.findLine(start, end, hitInfo);
142 
143  porousCelli = 0;
144 
145  forAll(porousCells, celli)
146  {
147  const pointIndexHit& hit = hitInfo[celli];
148 
149  if (hit.hit())
150  {
151  porousCells[porousCelli] = porousCells[celli];
152 
153  zTop[porousCelli] =
154  (hit.hitPoint() - C[porousCells[porousCelli]]) & zDir;
155 
156  porousCelli++;
157  }
158  }
159 
160  // Resize the porous cells list and height field
161  porousCells.setSize(porousCelli);
162  zTop.setSize(porousCelli);
163 
164  // Create a triSurface for the ground patch(s)
165  triSurface groundSurface
166  (
168  (
169  mesh.boundaryMesh(),
170  mesh.boundaryMesh().patchSet(groundPatches),
171  searchBounds
172  )
173  );
174 
175  // Combine the ground triSurfaces across all processors
176  if (Pstream::parRun())
177  {
178  List<List<labelledTri>> groundSurfaceProcTris(Pstream::nProcs());
179  List<pointField> groundSurfaceProcPoints(Pstream::nProcs());
180 
181  groundSurfaceProcTris[Pstream::myProcNo()] = groundSurface;
182  groundSurfaceProcPoints[Pstream::myProcNo()] = groundSurface.points();
183 
184  Pstream::gatherList(groundSurfaceProcTris);
185  Pstream::scatterList(groundSurfaceProcTris);
186  Pstream::gatherList(groundSurfaceProcPoints);
187  Pstream::scatterList(groundSurfaceProcPoints);
188 
189  label nTris = 0;
190  forAll(groundSurfaceProcTris, i)
191  {
192  nTris += groundSurfaceProcTris[i].size();
193  }
194 
195  List<labelledTri> groundSurfaceTris(nTris);
196  label trii = 0;
197  label offset = 0;
198  forAll(groundSurfaceProcTris, i)
199  {
200  forAll(groundSurfaceProcTris[i], j)
201  {
202  groundSurfaceTris[trii] = groundSurfaceProcTris[i][j];
203  groundSurfaceTris[trii][0] += offset;
204  groundSurfaceTris[trii][1] += offset;
205  groundSurfaceTris[trii][2] += offset;
206  trii++;
207  }
208  offset += groundSurfaceProcPoints[i].size();
209  }
210 
211  label nPoints = 0;
212  forAll(groundSurfaceProcPoints, i)
213  {
214  nPoints += groundSurfaceProcPoints[i].size();
215  }
216 
217  pointField groundSurfacePoints(nPoints);
218 
219  label pointi = 0;
220  forAll(groundSurfaceProcPoints, i)
221  {
222  forAll(groundSurfaceProcPoints[i], j)
223  {
224  groundSurfacePoints[pointi++] = groundSurfaceProcPoints[i][j];
225  }
226  }
227 
228  groundSurface = triSurface(groundSurfaceTris, groundSurfacePoints);
229  }
230 
231  // Create a searchable triSurface for the ground
232  triSurfaceSearch groundSearch(groundSurface);
233 
234  // Search the porous cells for the corresponding point on the ground
235 
236  start.setSize(porousCelli);
237  end.setSize(porousCelli);
238 
239  forAll(porousCells, porousCelli)
240  {
241  start[porousCelli] = C[porousCells[porousCelli]];
242  end[porousCelli] = start[porousCelli] - searchSpan*zDir;
243  }
244 
245  groundSearch.findLine(start, end, hitInfo);
246 
247  scalarField zBottom(porousCelli);
248 
249  forAll(porousCells, porousCelli)
250  {
251  const pointIndexHit& hit = hitInfo[porousCelli];
252 
253  if (hit.hit())
254  {
255  zBottom[porousCelli] =
256  (C[porousCells[porousCelli]] - hit.hitPoint()) & zDir;
257  }
258  }
259 
260  // Create the normalised height field
261  scalarField zNorm(zBottom/(zBottom + zTop));
262 
263  // Create the porosity surface area per unit volume zone field
264  Av_ = AvFunc->value(zNorm);
265 
266  // Create the porous region cellZone and add to the mesh cellZones
267 
268  meshCellZones& cellZones = const_cast<meshCellZones&>(mesh.cellZones());
269 
270  label zoneID = cellZones.findZoneID(zoneName_);
271 
272  if (zoneID == -1)
273  {
274  zoneID = cellZones.size();
275  cellZones.setSize(zoneID + 1);
276 
277  cellZones.set
278  (
279  zoneID,
280  new cellZone
281  (
282  zoneName_,
283  porousCells,
284  zoneID,
285  cellZones
286  )
287  );
288  }
289  else
290  {
292  << "Unable to create porous cellZone " << zoneName_
293  << ": zone already exists"
294  << abort(FatalError);
295  }
296 }
297 
298 
300 (
301  const word& name,
302  const word& modelType,
303  const fvMesh& mesh,
304  const dictionary& dict,
305  const word& dummyCellZoneName
306 )
307 :
308  powerLawLopesdaCostaZone(name, modelType, mesh, dict),
310  (
311  name,
312  modelType,
313  mesh,
314  dict,
315  powerLawLopesdaCostaZone::zoneName_
316  ),
317  Cd_(coeffs_.lookup<scalar>("Cd")),
318  C1_(coeffs_.lookup<scalar>("C1")),
319  rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho"))
320 {}
321 
322 
323 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
324 
326 {}
327 
328 
329 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
330 
331 const Foam::scalarField&
333 {
334  return Av_;
335 }
336 
337 
339 {}
340 
341 
343 (
344  const volVectorField& U,
345  const volScalarField& rho,
346  const volScalarField& mu,
347  vectorField& force
348 ) const
349 {
350  scalarField Udiag(U.size(), 0.0);
351  const scalarField& V = mesh_.V();
352 
353  apply(Udiag, V, rho, U);
354 
355  force = Udiag*U;
356 }
357 
358 
360 (
362 ) const
363 {
364  const vectorField& U = UEqn.psi();
365  const scalarField& V = mesh_.V();
366  scalarField& Udiag = UEqn.diag();
367 
368  if (UEqn.dimensions() == dimForce)
369  {
370  const volScalarField& rho =
371  mesh_.lookupObject<volScalarField>(rhoName_);
372 
373  apply(Udiag, V, rho, U);
374  }
375  else
376  {
377  apply(Udiag, V, geometricOneField(), U);
378  }
379 }
380 
381 
383 (
385  const volScalarField& rho,
386  const volScalarField& mu
387 ) const
388 {
389  const vectorField& U = UEqn.psi();
390  const scalarField& V = mesh_.V();
391  scalarField& Udiag = UEqn.diag();
392 
393  apply(Udiag, V, rho, U);
394 }
395 
396 
398 (
399  const fvVectorMatrix& UEqn,
400  volTensorField& AU
401 ) const
402 {
403  const vectorField& U = UEqn.psi();
404 
405  if (UEqn.dimensions() == dimForce)
406  {
407  const volScalarField& rho =
408  mesh_.lookupObject<volScalarField>(rhoName_);
409 
410  apply(AU, rho, U);
411  }
412  else
413  {
414  apply(AU, geometricOneField(), U);
415  }
416 }
417 
418 
420 {
421  os << indent << name_ << endl;
422  dict_.write(os);
423 
424  return true;
425 }
426 
427 
428 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
Graphite solid properties.
Definition: C.H:51
Run-time selectable general function of one variable.
Definition: Function1.H:64
Generic GeometricField class.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:99
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void setSize(const label)
Reset size of List.
Definition: List.C:281
label findZoneID(const word &zoneName) const
Find zone index given a name.
Definition: MeshZones.C:341
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:57
virtual Ostream & write(const char)=0
Write character.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:54
const Point & hitPoint() const
Return hit point.
bool hit() const
Is there a hit.
const Field< PointType > & points() const
Return reference to global points.
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
static const word & constant()
Return constant name.
Definition: TimePaths.H:122
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: autoPtr.H:51
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:59
const point & min() const
Minimum point defining the bounding box.
Definition: boundBoxI.H:54
const point & max() const
Maximum point defining the bounding box.
Definition: boundBoxI.H:60
bool contains(const point &) const
Contains point? (inside or on edge)
Definition: boundBoxI.H:170
A subset of mesh cells.
Definition: cellZone.H:64
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:160
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:860
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:118
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:101
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:406
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
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.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:403
const meshCellZones & cellZones() const
Return cell zones.
Definition: polyMesh.H:451
Top level model for porosity models.
Definition: porosityModel.H:57
const word zoneName_
Automatically generated zone name for this porous zone.
const scalarField & Av() const
Return the porosity surface area per unit volume zone field.
scalarField Av_
Porosity surface area per unit volume zone field.
powerLawLopesdaCostaZone(const word &name, const word &modelType, const fvMesh &mesh, const dictionary &dict)
Constructor.
Variant of the power law porosity model with spatially varying drag coefficient.
powerLawLopesdaCosta(const word &name, const word &modelType, const fvMesh &mesh, const dictionary &dict, const word &cellZoneName)
virtual void calcForce(const volVectorField &U, const volScalarField &rho, const volScalarField &mu, vectorField &force) const
Calculate the porosity force.
virtual void correct(fvVectorMatrix &UEqn) const
Add resistance.
virtual void calcTransformModelData()
Transform the model data wrt mesh changes.
const vectorField & cellCentres() const
static const word & geometryDir()
Return the geometry directory name.
A surface geometry formed of discrete facets, e.g. triangles and/or quadrilaterals,...
virtual void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &) const
Find first intersection on segment from start to end.
virtual tmp< pointField > points() const
Get the points that define the surface.
Helper class to search on triSurface.
void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &info) const
static triSurface triangulate(const polyBoundaryMesh &mBesh, const labelHashSet &includePatches, const bool verbose=false)
Simple triangulation of (selected patches of) boundaryMesh. Needs.
Triangulated surface description with patch information.
Definition: triSurface.H:69
A class for handling words, derived from string.
Definition: word.H:62
fvVectorMatrix & UEqn
Definition: UEqn.H:11
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
A special matrix type and solver, designed for finite volume solutions of scalar equations.
label nPoints
U
Definition: pEqn.H:72
const dimensionedScalar mu
Atomic mass unit.
addToRunTimeSelectionTable(porosityModel, powerLawLopesdaCosta, mesh)
defineTypeNameAndDebug(powerLawLopesdaCosta, 0)
Namespace for OpenFOAM.
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const dimensionSet dimForce
error FatalError
void offset(label &lst, const label o)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:221
dictionary dict