DarcyForchheimer.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) 2012-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 
27 #include "DarcyForchheimer.H"
28 #include "geometricOneField.H"
29 #include "fvMatrices.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  namespace porosityModels
36  {
37  defineTypeNameAndDebug(DarcyForchheimer, 0);
38  addToRunTimeSelectionTable(porosityModel, DarcyForchheimer, mesh);
39  }
40 }
41 
42 
43 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
44 
46 (
47  const word& name,
48  const word& modelType,
49  const fvMesh& mesh,
50  const dictionary& dict,
51  const word& cellZoneName
52 )
53 :
54  porosityModel(name, modelType, mesh, dict, cellZoneName),
55  dXYZ_("d", dimless/sqr(dimLength), coeffs_),
56  fXYZ_("f", dimless/dimLength, coeffs_),
57  D_(cellZoneIDs_.size()),
58  F_(cellZoneIDs_.size()),
59  rhoName_(coeffs_.lookupOrDefault<word>("rho", "rho")),
60  muName_(coeffs_.lookupOrDefault<word>("mu", "thermo:mu")),
61  nuName_(coeffs_.lookupOrDefault<word>("nu", "nu"))
62 {
63  adjustNegativeResistance(dXYZ_);
64  adjustNegativeResistance(fXYZ_);
65 
66  calcTransformModelData();
67 }
68 
69 
70 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
71 
73 {}
74 
75 
76 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
77 
79 {
80  if (coordSys_.R().uniform())
81  {
82  forAll(cellZoneIDs_, zoneI)
83  {
84  D_[zoneI].setSize(1);
85  F_[zoneI].setSize(1);
86 
87  D_[zoneI][0] = Zero;
88  D_[zoneI][0].xx() = dXYZ_.value().x();
89  D_[zoneI][0].yy() = dXYZ_.value().y();
90  D_[zoneI][0].zz() = dXYZ_.value().z();
91 
92  D_[zoneI][0] = coordSys_.R().transformTensor(D_[zoneI][0]);
93 
94  // leading 0.5 is from 1/2*rho
95  F_[zoneI][0] = Zero;
96  F_[zoneI][0].xx() = 0.5*fXYZ_.value().x();
97  F_[zoneI][0].yy() = 0.5*fXYZ_.value().y();
98  F_[zoneI][0].zz() = 0.5*fXYZ_.value().z();
99 
100  F_[zoneI][0] = coordSys_.R().transformTensor(F_[zoneI][0]);
101  }
102  }
103  else
104  {
105  forAll(cellZoneIDs_, zoneI)
106  {
107  const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
108 
109  D_[zoneI].setSize(cells.size());
110  F_[zoneI].setSize(cells.size());
111 
112  forAll(cells, i)
113  {
114  D_[zoneI][i] = Zero;
115  D_[zoneI][i].xx() = dXYZ_.value().x();
116  D_[zoneI][i].yy() = dXYZ_.value().y();
117  D_[zoneI][i].zz() = dXYZ_.value().z();
118 
119  // leading 0.5 is from 1/2*rho
120  F_[zoneI][i] = Zero;
121  F_[zoneI][i].xx() = 0.5*fXYZ_.value().x();
122  F_[zoneI][i].yy() = 0.5*fXYZ_.value().y();
123  F_[zoneI][i].zz() = 0.5*fXYZ_.value().z();
124  }
125 
126  const coordinateRotation& R = coordSys_.R
127  (
128  UIndirectList<vector>(mesh_.C(), cells)()
129  );
130 
131  D_[zoneI] = R.transformTensor(D_[zoneI]);
132  F_[zoneI] = R.transformTensor(F_[zoneI]);
133  }
134  }
135 
136  if (debug && (mesh_.time().writeTime() || mesh_.time().timeIndex() == 0))
137  {
138  volTensorField Dout
139  (
140  IOobject
141  (
142  typeName + ":D",
143  mesh_.time().timeName(),
144  mesh_,
147  ),
148  mesh_,
149  dimensionedTensor(dXYZ_.dimensions(), Zero)
150  );
151  volTensorField Fout
152  (
153  IOobject
154  (
155  typeName + ":F",
156  mesh_.time().timeName(),
157  mesh_,
160  ),
161  mesh_,
162  dimensionedTensor(fXYZ_.dimensions(), Zero)
163  );
164 
165  UIndirectList<tensor>(Dout, mesh_.cellZones()[cellZoneIDs_[0]]) = D_[0];
166  UIndirectList<tensor>(Fout, mesh_.cellZones()[cellZoneIDs_[0]]) = F_[0];
167 
168  Dout.write();
169  Fout.write();
170  }
171 }
172 
173 
175 (
176  const volVectorField& U,
177  const volScalarField& rho,
178  const volScalarField& mu,
179  vectorField& force
180 ) const
181 {
182  scalarField Udiag(U.size(), 0.0);
183  vectorField Usource(U.size(), Zero);
184  const scalarField& V = mesh_.V();
185 
186  apply(Udiag, Usource, V, rho, mu, U);
187 
188  force = Udiag*U - Usource;
189 }
190 
191 
193 (
194  fvVectorMatrix& UEqn
195 ) const
196 {
197  const volVectorField& U = UEqn.psi();
198  const scalarField& V = mesh_.V();
199  scalarField& Udiag = UEqn.diag();
200  vectorField& Usource = UEqn.source();
201 
202  word rhoName(IOobject::groupName(rhoName_, U.group()));
203  word muName(IOobject::groupName(muName_, U.group()));
204  word nuName(IOobject::groupName(nuName_, U.group()));
205 
206  if (UEqn.dimensions() == dimForce)
207  {
208  const volScalarField& rho = mesh_.lookupObject<volScalarField>(rhoName);
209 
210  if (mesh_.foundObject<volScalarField>(muName))
211  {
212  const volScalarField& mu =
213  mesh_.lookupObject<volScalarField>(muName);
214 
215  apply(Udiag, Usource, V, rho, mu, U);
216  }
217  else
218  {
219  const volScalarField& nu =
220  mesh_.lookupObject<volScalarField>(nuName);
221 
222  apply(Udiag, Usource, V, rho, rho*nu, U);
223  }
224  }
225  else
226  {
227  if (mesh_.foundObject<volScalarField>(nuName))
228  {
229  const volScalarField& nu =
230  mesh_.lookupObject<volScalarField>(nuName);
231 
232  apply(Udiag, Usource, V, geometricOneField(), nu, U);
233  }
234  else
235  {
236  const volScalarField& rho =
237  mesh_.lookupObject<volScalarField>(rhoName);
238  const volScalarField& mu =
239  mesh_.lookupObject<volScalarField>(muName);
240 
241  apply(Udiag, Usource, V, geometricOneField(), mu/rho, U);
242  }
243  }
244 }
245 
246 
248 (
249  const fvVectorMatrix& UEqn,
250  volTensorField& AU
251 ) const
252 {
253  const volVectorField& U = UEqn.psi();
254 
255  word rhoName(IOobject::groupName(rhoName_, U.group()));
256  word muName(IOobject::groupName(muName_, U.group()));
257  word nuName(IOobject::groupName(nuName_, U.group()));
258 
259  if (UEqn.dimensions() == dimForce)
260  {
261  const volScalarField& rho = mesh_.lookupObject<volScalarField>(rhoName);
262  const volScalarField& mu = mesh_.lookupObject<volScalarField>(muName);
263 
264  apply(AU, rho, mu, U);
265  }
266  else
267  {
268  if (mesh_.foundObject<volScalarField>(nuName))
269  {
270  const volScalarField& nu =
271  mesh_.lookupObject<volScalarField>(nuName);
272 
273  apply(AU, geometricOneField(), nu, U);
274  }
275  else
276  {
277  const volScalarField& rho =
278  mesh_.lookupObject<volScalarField>(rhoName);
279  const volScalarField& mu =
280  mesh_.lookupObject<volScalarField>(muName);
281 
282  apply(AU, geometricOneField(), mu/rho, U);
283  }
284  }
285 }
286 
287 
289 {
290  os << indent << name_ << endl;
291  dict_.write(os);
292 
293  return true;
294 }
295 
296 
297 // ************************************************************************* //
static word group(const word &name)
Return group (extension part of name)
Definition: IOobject.C:134
Abstract base class for coordinate rotation.
bool writeData(Ostream &os) const
Write.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual Ostream & write(const char)=0
Write character.
defineTypeNameAndDebug(powerLawLopesdaCosta, 0)
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:221
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
addToRunTimeSelectionTable(porosityModel, powerLawLopesdaCosta, mesh)
const GeometricField< Type, fvPatchField, volMesh > & psi() const
Definition: fvMatrix.H:282
dimensionedSymmTensor sqr(const dimensionedVector &dv)
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const dimensionSet dimless
Macros for easy insertion into run-time selection tables.
const dimensionSet dimLength
DarcyForchheimer(const word &name, const word &modelType, const fvMesh &mesh, const dictionary &dict, const word &cellZoneName)
virtual void calcTransformModelData()
Transform the model data wrt mesh changes.
A class representing the concept of a GeometricField of 1 used to avoid unnecessary manipulations for...
virtual void calcForce(const volVectorField &U, const volScalarField &rho, const volScalarField &mu, vectorField &force) const
Calculate the porosity force.
const cellShapeList & cells
A class for handling words, derived from string.
Definition: word.H:59
static word groupName(Name name, const word &group)
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:72
static const zero Zero
Definition: zero.H:97
dimensioned< tensor > dimensionedTensor
Dimensioned tensor obtained from generic dimensioned type.
const dimensionSet dimForce
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
Field< Type > & source()
Definition: fvMatrix.H:292
#define R(A, B, C, D, E, F, K, M)
virtual void correct(fvVectorMatrix &UEqn) const
Add resistance.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
A List with indirect addressing.
Definition: fvMatrix.H:106
A special matrix type and solver, designed for finite volume solutions of scalar equations.
virtual tensor transformTensor(const tensor &st) const =0
Transform tensor using transformation tensorField.
scalarField & diag()
Definition: lduMatrix.C:186
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Top level model for porosity models.
Definition: porosityModel.H:54
Namespace for OpenFOAM.
const dimensionSet & dimensions() const
Definition: fvMatrix.H:287