multiphaseMangrovesTurbulenceModel.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) 2017 IH-Cantabria
9 Copyright (C) 2017-2021 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 "fvMesh.H"
32#include "fvMatrices.H"
33#include "fvmSup.H"
35
36// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37
38namespace Foam
39{
40namespace fv
41{
44 (
45 option,
48 );
49}
50}
51
52
53// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
54
57(
58 const volVectorField& U
59) const
60{
61 auto tkCoeff = tmp<volScalarField>::New
62 (
64 (
65 typeName + ":kCoeff",
67 mesh_.time(),
70 ),
71 mesh_,
73 );
74
75 volScalarField& kCoeff = tkCoeff.ref();
76
78 {
79 const scalar a = aZone_[i];
80 const scalar N = NZone_[i];
81 const scalar Ckp = CkpZone_[i];
82 const scalar Cd = CdZone_[i];
83
84 for (label zonei : zoneIDs_[i])
85 {
86 const cellZone& cz = mesh_.cellZones()[zonei];
87
88 for (label celli : cz)
89 {
90 kCoeff[celli] = Ckp*Cd*a*N*mag(U[celli]);
91 }
92 }
93 }
94
95 kCoeff.correctBoundaryConditions();
96
97 return tkCoeff;
98}
99
100
103(
104 const volVectorField& U
105) const
106{
107 auto tepsilonCoeff = tmp<volScalarField>::New
108 (
110 (
111 typeName + ":epsilonCoeff",
112 mesh_.time().timeName(),
113 mesh_.time(),
116 ),
117 mesh_,
119 );
120
121 volScalarField& epsilonCoeff = tepsilonCoeff.ref();
122
123 forAll(zoneIDs_, i)
124 {
125 const scalar a = aZone_[i];
126 const scalar N = NZone_[i];
127 const scalar Cep = CepZone_[i];
128 const scalar Cd = CdZone_[i];
129
130 for (label zonei : zoneIDs_[i])
131 {
132 const cellZone& cz = mesh_.cellZones()[zonei];
133
134 for (label celli : cz)
135 {
136 epsilonCoeff[celli] = Cep*Cd*a*N*mag(U[celli]);
137 }
138 }
139 }
140
141 epsilonCoeff.correctBoundaryConditions();
142
143 return tepsilonCoeff;
144}
145
146
147// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
148
150(
151 const word& name,
152 const word& modelType,
153 const dictionary& dict,
154 const fvMesh& mesh
155)
156:
157 fv::option(name, modelType, dict, mesh),
158 aZone_(),
159 NZone_(),
160 CkpZone_(),
161 CepZone_(),
162 CdZone_(),
163 UName_("U"),
164 kName_("k"),
165 epsilonName_("epsilon")
166{
167 read(dict);
168}
169
170
171// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
172
174(
175 fvMatrix<scalar>& eqn,
176 const label fieldi
177)
178{
179 const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
180
181 if (eqn.psi().name() == epsilonName_)
182 {
183 fvMatrix<scalar> epsilonEqn
184 (
185 - fvm::Sp(epsilonCoeff(U), eqn.psi())
186 );
187 eqn += epsilonEqn;
188 }
189 else if (eqn.psi().name() == kName_)
190 {
192 (
193 - fvm::Sp(kCoeff(U), eqn.psi())
194 );
195 eqn += kEqn;
196 }
197}
198
199
201(
202 const volScalarField& rho,
203 fvMatrix<scalar>& eqn,
204 const label fieldi
205)
206{
207 const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
208
209 if (eqn.psi().name() == epsilonName_)
210 {
211 fvMatrix<scalar> epsilonEqn
212 (
213 - fvm::Sp(rho*epsilonCoeff(U), eqn.psi())
214 );
215 eqn += epsilonEqn;
216 }
217 else if (eqn.psi().name() == kName_)
218 {
220 (
221 - fvm::Sp(rho*kCoeff(U), eqn.psi())
222 );
223 eqn += kEqn;
224 }
225}
226
227
229{
231 {
232 if (!coeffs_.readIfPresent("epsilonNames", fieldNames_))
233 {
234 if (coeffs_.found("epsilon"))
235 {
236 fieldNames_.resize(1);
237 coeffs_.readEntry("epsilon", fieldNames_.first());
238 }
239 else
240 {
241 fieldNames_.resize(2);
242 fieldNames_[0] = "epsilon";
243 fieldNames_[1] = "k";
244 }
245 }
247
248 // Create the Mangroves models - 1 per region
249 const dictionary& regionsDict(coeffs_.subDict("regions"));
250 const wordList regionNames(regionsDict.toc());
251 aZone_.setSize(regionNames.size(), 1);
252 NZone_.setSize(regionNames.size(), 1);
253 CkpZone_.setSize(regionNames.size(), 1);
254 CepZone_.setSize(regionNames.size(), 1);
255 CdZone_.setSize(regionNames.size(), 1);
256 zoneIDs_.setSize(regionNames.size());
257
258 forAll(zoneIDs_, i)
259 {
260 const word& regionName = regionNames[i];
261 const dictionary& modelDict = regionsDict.subDict(regionName);
262
263 const word zoneName(modelDict.get<word>("cellZone"));
264
265 zoneIDs_[i] = mesh_.cellZones().indices(zoneName);
266 if (zoneIDs_[i].empty())
267 {
269 << "Unable to find cellZone " << zoneName << nl
270 << "Valid cellZones are:" << mesh_.cellZones().names()
271 << exit(FatalError);
272 }
273
274 modelDict.readEntry("a", aZone_[i]);
275 modelDict.readEntry("N", NZone_[i]);
276 modelDict.readEntry("Ckp", CkpZone_[i]);
277 modelDict.readEntry("Cep", CepZone_[i]);
278 modelDict.readEntry("Cd", CdZone_[i]);
279 }
280
281 return true;
282 }
283
284 return false;
285}
286
287
288// ************************************************************************* //
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
Internal & ref(const bool updateAccessTime=true)
Return a reference to the dimensioned internal field.
void correctBoundaryConditions()
Correct boundary field.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
virtual bool read()
Re-read model coefficients if they have changed.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
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 & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Definition: dictionary.C:460
wordList toc() const
Return the table of contents.
Definition: dictionary.C:602
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvMatrix.H:121
const GeometricField< Type, fvPatchField, volMesh > & psi(const label i=0) const
Return psi.
Definition: fvMatrix.H:412
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
tmp< volScalarField > kCoeff(const volVectorField &U) const
Return the k coefficient.
virtual void addSup(fvMatrix< scalar > &eqn, const label fieldi)
Add implicit contribution to momentum equation.
virtual bool read(const dictionary &dict)
Read dictionary.
scalarList aZone_
Width of the vegetation element.
scalarList NZone_
Number of plants per unit of area.
tmp< volScalarField > epsilonCoeff(const volVectorField &U) const
Return the epsilon coefficient.
Base abstract class for handling finite volume options (i.e. fvOption).
Definition: fvOption.H:127
const fvMesh & mesh_
Reference to the mesh database.
Definition: fvOption.H:139
void resetApplied()
Resize/reset applied flag list for all fieldNames_ entries.
Definition: fvOption.C:48
const cellZoneMesh & cellZones() const noexcept
Return cell zone mesh.
Definition: polyMesh.H:504
A class for managing temporary objects.
Definition: tmp.H:65
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
dynamicFvMesh & mesh
Foam::word regionName(Foam::polyMesh::defaultRegion)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Calculate the finiteVolume matrix for implicit and explicit sources.
wordList regionNames
zeroField Sp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
Namespace for OpenFOAM.
const dimensionSet dimless
Dimensionless.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList fv(nPoints)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
const Vector< label > N(dict.get< Vector< label > >("N"))