powerLawLopesdaCosta.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2018 OpenFOAM Foundation
9 Copyright (C) 2020-2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
31#include "geometricOneField.H"
32#include "fvMatrices.H"
33#include "triSurfaceMesh.H"
34#include "triSurfaceTools.H"
35
36// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37
38namespace Foam
39{
40 namespace porosityModels
41 {
44 }
45}
46
47
48// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
49
51(
52 const word& name,
53 const word& modelType,
54 const fvMesh& mesh,
55 const dictionary& dict
56)
57:
58 zoneName_(name + ":porousZone")
59{
60 const dictionary& coeffs(dict.optionalSubDict(modelType + "Coeffs"));
61
62 // Vertical direction within porous region
63 vector zDir(coeffs.get<vector>("zDir"));
64
65 // Span of the search for the cells in the porous region
66 scalar searchSpan(coeffs.get<scalar>("searchSpan"));
67
68 // Top surface file name defining the extent of the porous region
69 const word topSurfaceFileName(coeffs.get<word>("topSurface"));
70
71 // List of the ground patches defining the lower surface
72 // of the porous region
73 wordRes groundPatches(coeffs.get<wordRes>("groundPatches"));
74
75 // Functional form of the porosity surface area per unit volume as a
76 // function of the normalized vertical position
78 (
80 );
81
82 // Searchable triSurface for the top of the porous region
83 triSurfaceMesh searchSurf
84 (
86 (
87 topSurfaceFileName,
88 mesh.time().constant(),
89 "triSurface",
90 mesh.time()
91 )
92 );
93
94 // Limit the porous cell search to those within the lateral and vertical
95 // extent of the top surface
96
97 boundBox surfBounds(searchSurf.points());
98 boundBox searchBounds
99 (
100 surfBounds.min() - searchSpan*zDir, surfBounds.max()
101 );
102
103 const pointField& C = mesh.cellCentres();
104
105 // List of cells within the porous region
106 labelList porousCells(C.size());
107 label porousCelli = 0;
108
109 forAll(C, celli)
110 {
111 if (searchBounds.contains(C[celli]))
112 {
113 porousCells[porousCelli++] = celli;
114 }
115 }
116
117 porousCells.setSize(porousCelli);
118
119 pointField start(porousCelli);
120 pointField end(porousCelli);
121
122 forAll(porousCells, porousCelli)
123 {
124 start[porousCelli] = C[porousCells[porousCelli]];
125 end[porousCelli] = start[porousCelli] + searchSpan*zDir;
126 }
127
128 // Field of distance between the cell centre and the corresponding point
129 // on the porous region top surface
130 scalarField zTop(porousCelli);
131
132 // Search the porous cells for those with a corresponding point on the
133 // porous region top surface
134 List<pointIndexHit> hitInfo;
135 searchSurf.findLine(start, end, hitInfo);
136
137 porousCelli = 0;
138
139 forAll(porousCells, celli)
140 {
141 const pointIndexHit& hit = hitInfo[celli];
142
143 if (hit.hit())
144 {
145 porousCells[porousCelli] = porousCells[celli];
146
147 zTop[porousCelli] =
148 (hit.hitPoint() - C[porousCells[porousCelli]]) & zDir;
149
150 porousCelli++;
151 }
152 }
153
154 // Resize the porous cells list and height field
155 porousCells.setSize(porousCelli);
156 zTop.setSize(porousCelli);
157
158 // Create a triSurface for the ground patch(s)
159 triSurface groundSurface
160 (
162 (
164 mesh.boundaryMesh().patchSet(groundPatches),
165 searchBounds
166 )
167 );
168
169 // Combine the ground triSurfaces across all processors
170 if (Pstream::parRun())
171 {
172 List<List<labelledTri>> groundSurfaceProcTris(Pstream::nProcs());
173 List<pointField> groundSurfaceProcPoints(Pstream::nProcs());
174
175 groundSurfaceProcTris[Pstream::myProcNo()] = groundSurface;
176 groundSurfaceProcPoints[Pstream::myProcNo()] = groundSurface.points();
177
178 Pstream::allGatherList(groundSurfaceProcTris);
179 Pstream::allGatherList(groundSurfaceProcPoints);
180
181 label nTris = 0;
182 forAll(groundSurfaceProcTris, i)
183 {
184 nTris += groundSurfaceProcTris[i].size();
185 }
186
187 List<labelledTri> groundSurfaceTris(nTris);
188 label trii = 0;
189 label offset = 0;
190 forAll(groundSurfaceProcTris, i)
191 {
192 forAll(groundSurfaceProcTris[i], j)
193 {
194 groundSurfaceTris[trii] = groundSurfaceProcTris[i][j];
195 groundSurfaceTris[trii][0] += offset;
196 groundSurfaceTris[trii][1] += offset;
197 groundSurfaceTris[trii][2] += offset;
198 trii++;
199 }
200 offset += groundSurfaceProcPoints[i].size();
201 }
202
203 label nPoints = 0;
204 forAll(groundSurfaceProcPoints, i)
205 {
206 nPoints += groundSurfaceProcPoints[i].size();
207 }
208
209 pointField groundSurfacePoints(nPoints);
210
211 label pointi = 0;
212 forAll(groundSurfaceProcPoints, i)
213 {
214 forAll(groundSurfaceProcPoints[i], j)
215 {
216 groundSurfacePoints[pointi++] = groundSurfaceProcPoints[i][j];
217 }
218 }
219
220 groundSurface = triSurface(groundSurfaceTris, groundSurfacePoints);
221 }
222
223 // Create a searchable triSurface for the ground
224 triSurfaceSearch groundSearch(groundSurface);
225
226 // Search the porous cells for the corresponding point on the ground
227
228 start.setSize(porousCelli);
229 end.setSize(porousCelli);
230
231 forAll(porousCells, porousCelli)
232 {
233 start[porousCelli] = C[porousCells[porousCelli]];
234 end[porousCelli] = start[porousCelli] - searchSpan*zDir;
235 }
236
237 groundSearch.findLine(start, end, hitInfo);
238
239 scalarField zBottom(porousCelli);
240
241 forAll(porousCells, porousCelli)
242 {
243 const pointIndexHit& hit = hitInfo[porousCelli];
244
245 if (hit.hit())
246 {
247 zBottom[porousCelli] =
248 (C[porousCells[porousCelli]] - hit.hitPoint()) & zDir;
249 }
250 }
251
252 // Create the normalized height field
253 scalarField zNorm(zBottom/(zBottom + zTop));
254
255 // Create the porosity surface area per unit volume zone field
256 Sigma_ = SigmaFunc->value(zNorm);
257
258 // Create the porous region cellZone and add to the mesh cellZones
259
260 cellZoneMesh& cellZones = const_cast<cellZoneMesh&>(mesh.cellZones());
261
262 label zoneID = cellZones.findZoneID(zoneName_);
263
264 if (zoneID == -1)
265 {
266 zoneID = cellZones.size();
267 cellZones.setSize(zoneID + 1);
268
269 cellZones.set
270 (
271 zoneID,
272 new cellZone
273 (
274 zoneName_,
275 porousCells,
276 zoneID,
277 cellZones
278 )
279 );
280 }
281 else
282 {
284 << "Unable to create porous cellZone " << zoneName_
285 << ": zone already exists"
286 << abort(FatalError);
287 }
288}
289
290
292(
293 const word& name,
294 const word& modelType,
295 const fvMesh& mesh,
296 const dictionary& dict,
297 const word& dummyCellZoneName
298)
299:
302 (
303 name,
304 modelType,
305 mesh,
306 dict,
307 powerLawLopesdaCostaZone::zoneName_
308 ),
309 Cd_(coeffs_.get<scalar>("Cd")),
310 C1_(coeffs_.get<scalar>("C1")),
311 rhoName_(coeffs_.getOrDefault<word>("rho", "rho"))
312{}
313
314
315// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
316
318{}
319
320
321// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
322
325{
326 return Sigma_;
327}
328
329
331{}
332
333
335(
336 const volVectorField& U,
337 const volScalarField& rho,
338 const volScalarField& mu,
339 vectorField& force
340) const
341{
342 scalarField Udiag(U.size(), Zero);
343 const scalarField& V = mesh_.V();
344
345 apply(Udiag, V, rho, U);
346
347 force = Udiag*U;
348}
349
350
352(
354) const
355{
356 const vectorField& U = UEqn.psi();
357 const scalarField& V = mesh_.V();
358 scalarField& Udiag = UEqn.diag();
359
360 if (UEqn.dimensions() == dimForce)
361 {
362 const volScalarField& rho =
363 mesh_.lookupObject<volScalarField>(rhoName_);
364
365 apply(Udiag, V, rho, U);
366 }
367 else
368 {
369 apply(Udiag, V, geometricOneField(), U);
370 }
371}
372
373
375(
377 const volScalarField& rho,
378 const volScalarField& mu
379) const
380{
381 const vectorField& U = UEqn.psi();
382 const scalarField& V = mesh_.V();
383 scalarField& Udiag = UEqn.diag();
384
385 apply(Udiag, V, rho, U);
386}
387
388
390(
391 const fvVectorMatrix& UEqn,
393) const
394{
395 const vectorField& U = UEqn.psi();
396
397 if (UEqn.dimensions() == dimForce)
398 {
399 const volScalarField& rho =
400 mesh_.lookupObject<volScalarField>(rhoName_);
401
402 apply(AU, rho, U);
403 }
404 else
405 {
406 apply(AU, geometricOneField(), U);
407 }
408}
409
410
412{
413 dict_.writeEntry(name_, os);
414
415 return true;
416}
417
418
419// ************************************************************************* //
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
Graphite solid properties.
Definition: C.H:53
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
Definition: Function1.H:96
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
void setSize(const label n)
Alias for resize()
Definition: List.H:218
virtual label triangulate()
Triangulate in-place, returning the number of triangles added.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
bool hit() const noexcept
Is there a hit?
const point_type & hitPoint() const
Return hit point. Fatal if not hit.
const Field< point_type > & points() const noexcept
Return reference to global points.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
static void allGatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
const T * set(const label i) const
Definition: PtrList.H:138
void setSize(const label newLen)
Same as resize()
Definition: PtrList.H:151
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
const word & constant() const
Return constant name.
Definition: TimePathsI.H:96
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
label findZoneID(const word &zoneName) const
Find zone index by name, return -1 if not found.
Definition: ZoneMesh.C:525
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
bool contains(const point &pt) const
Contains point? (inside or on edge)
Definition: boundBoxI.H:271
A subset of mesh cells.
Definition: cellZone.H:65
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const dictionary & optionalSubDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary, otherwise return this dictionary.
Definition: dictionary.C:577
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:121
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
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 useGroups=true) const
Return the set of patch IDs corresponding to the given names.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
const cellZoneMesh & cellZones() const noexcept
Return cell zone mesh.
Definition: polyMesh.H:504
Top level model for porosity models.
Definition: porosityModel.H:60
const word zoneName_
Automatically generated zone name for this porous zone.
scalarField Sigma_
Porosity surface area per unit volume zone field.
const scalarField & Sigma() const
Return the porosity surface area per unit volume zone field.
Variant of the power law porosity model with spatially varying drag coefficient.
virtual void calcForce(const volVectorField &U, const volScalarField &rho, const volScalarField &mu, vectorField &force) const
Calculate the porosity force.
virtual void calcTransformModelData()
Transform the model data wrt mesh changes.
const vectorField & cellCentres() const
int myProcNo() const noexcept
Return processor number.
IOoject and searching on triSurface.
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
Triangulated surface description with patch information.
Definition: triSurface.H:79
A List of wordRe with additional matching capabilities.
Definition: wordRes.H:54
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
U
Definition: pEqn.H:72
fvVectorMatrix & UEqn
Definition: UEqn.H:13
const volScalarField & mu
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
A special matrix type and solver, designed for finite volume solutions of scalar equations.
label nPoints
const labelIOList & zoneID
Namespace for OpenFOAM.
static void apply(bitSet &selection, const Detail::parcelSelection::actionType action, const Predicate &accept, const UList< Type > &list, const AccessOp &aop)
errorManip< error > abort(error &err)
Definition: errorManip.H:144
const dimensionSet dimForce
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333