DarcyForchheimer.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 | Copyright (C) 2012-2014 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  {
39  }
40 }
41 
42 
43 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
44 
45 Foam::porosityModels::DarcyForchheimer::DarcyForchheimer
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_(coeffs_.lookup("d")),
56  fXYZ_(coeffs_.lookup("f")),
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  calcTranformModelData();
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] = tensor::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] = tensor::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] = tensor::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] = tensor::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(mesh_, cells);
127 
128  D_[zoneI] = R.transformTensor(D_[zoneI], cells);
129  F_[zoneI] = R.transformTensor(F_[zoneI], cells);
130  }
131  }
132 
133  if (debug && mesh_.time().outputTime())
134  {
135  volTensorField Dout
136  (
137  IOobject
138  (
139  typeName + ":D",
140  mesh_.time().timeName(),
141  mesh_,
144  ),
145  mesh_,
146  dimensionedTensor("0", dXYZ_.dimensions(), tensor::zero)
147  );
148  volTensorField Fout
149  (
150  IOobject
151  (
152  typeName + ":F",
153  mesh_.time().timeName(),
154  mesh_,
157  ),
158  mesh_,
159  dimensionedTensor("0", fXYZ_.dimensions(), tensor::zero)
160  );
161 
162  UIndirectList<tensor>(Dout, mesh_.cellZones()[cellZoneIDs_[0]]) = D_[0];
163  UIndirectList<tensor>(Fout, mesh_.cellZones()[cellZoneIDs_[0]]) = F_[0];
164 
165  Dout.write();
166  Fout.write();
167  }
168 }
169 
170 
172 (
173  const volVectorField& U,
174  const volScalarField& rho,
175  const volScalarField& mu,
176  vectorField& force
177 ) const
178 {
179  scalarField Udiag(U.size(), 0.0);
180  vectorField Usource(U.size(), vector::zero);
181  const scalarField& V = mesh_.V();
182 
183  apply(Udiag, Usource, V, rho, mu, U);
184 
185  force = Udiag*U - Usource;
186 }
187 
188 
190 (
191  fvVectorMatrix& UEqn
192 ) const
193 {
194  const volVectorField& U = UEqn.psi();
195  const scalarField& V = mesh_.V();
196  scalarField& Udiag = UEqn.diag();
197  vectorField& Usource = UEqn.source();
198 
199  word rhoName(IOobject::groupName(rhoName_, U.group()));
200  word muName(IOobject::groupName(muName_, U.group()));
201  word nuName(IOobject::groupName(nuName_, U.group()));
202 
203  if (UEqn.dimensions() == dimForce)
204  {
205  const volScalarField& rho = mesh_.lookupObject<volScalarField>(rhoName);
206 
207  if (mesh_.foundObject<volScalarField>(muName))
208  {
209  const volScalarField& mu =
210  mesh_.lookupObject<volScalarField>(muName);
211 
212  apply(Udiag, Usource, V, rho, mu, U);
213  }
214  else
215  {
216  const volScalarField& nu =
217  mesh_.lookupObject<volScalarField>(nuName);
218 
219  apply(Udiag, Usource, V, rho, rho*nu, U);
220  }
221  }
222  else
223  {
224  if (mesh_.foundObject<volScalarField>(nuName))
225  {
226  const volScalarField& nu =
227  mesh_.lookupObject<volScalarField>(nuName);
228 
229  apply(Udiag, Usource, V, geometricOneField(), nu, U);
230  }
231  else
232  {
233  const volScalarField& rho =
234  mesh_.lookupObject<volScalarField>(rhoName);
235  const volScalarField& mu =
236  mesh_.lookupObject<volScalarField>(muName);
237 
238  apply(Udiag, Usource, V, geometricOneField(), mu/rho, U);
239  }
240  }
241 }
242 
243 
245 (
246  fvVectorMatrix& UEqn,
247  const volScalarField& rho,
248  const volScalarField& mu
249 ) const
250 {
251  const vectorField& U = UEqn.psi();
252  const scalarField& V = mesh_.V();
253  scalarField& Udiag = UEqn.diag();
254  vectorField& Usource = UEqn.source();
255 
256  apply(Udiag, Usource, V, rho, mu, U);
257 }
258 
259 
261 (
262  const fvVectorMatrix& UEqn,
263  volTensorField& AU
264 ) const
265 {
266  const volVectorField& U = UEqn.psi();
267 
268  word rhoName(IOobject::groupName(rhoName_, U.group()));
269  word muName(IOobject::groupName(muName_, U.group()));
270  word nuName(IOobject::groupName(nuName_, U.group()));
271 
272  if (UEqn.dimensions() == dimForce)
273  {
274  const volScalarField& rho = mesh_.lookupObject<volScalarField>(rhoName);
275  const volScalarField& mu = mesh_.lookupObject<volScalarField>(muName);
276 
277  apply(AU, rho, mu, U);
278  }
279  else
280  {
281  if (mesh_.foundObject<volScalarField>(nuName))
282  {
283  const volScalarField& nu =
284  mesh_.lookupObject<volScalarField>(nuName);
285 
286  apply(AU, geometricOneField(), nu, U);
287  }
288  else
289  {
290  const volScalarField& rho =
291  mesh_.lookupObject<volScalarField>(rhoName);
292  const volScalarField& mu =
293  mesh_.lookupObject<volScalarField>(muName);
294 
295  apply(AU, geometricOneField(), mu/rho, U);
296  }
297  }
298 }
299 
300 
302 {
303  os << indent << name_ << endl;
304  dict_.write(os);
305 
306  return true;
307 }
308 
309 
310 // ************************************************************************* //
virtual Ostream & write(const token &)=0
Write next token to stream.
static word groupName(Name name, const word &group)
const dimensionSet & dimensions() const
Definition: fvMatrix.H:286
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static const Tensor zero
Definition: Tensor.H:80
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
#define forAll(list, i)
Definition: UList.H:421
virtual void calcTranformModelData()
Transform the model data wrt mesh changes.
const GeometricField< Type, fvPatchField, volMesh > & psi() const
Definition: fvMatrix.H:281
volScalarField & nu
A special matrix type and solver, designed for finite volume solutions of scalar equations.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
addToRunTimeSelectionTable(porosityModel, DarcyForchheimer, mesh)
Abstract base class for coordinate rotation.
virtual void correct(fvVectorMatrix &UEqn) const
Add resistance.
dimensioned< tensor > dimensionedTensor
Dimensioned tensor obtained from generic dimensioned type.
A class for handling words, derived from string.
Definition: word.H:59
scalarField & diag()
Definition: lduMatrix.C:183
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.
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:221
bool writeData(Ostream &os) const
Write.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
Top level model for porosity models.
Definition: porosityModel.H:55
Field< Type > & source()
Definition: fvMatrix.H:291
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:68
#define R(A, B, C, D, E, F, K, M)
A List with indirect addressing.
Definition: fvMatrix.H:106
const cellShapeList & cells
Macros for easy insertion into run-time selection tables.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
word group() const
Return group (extension part of name)
Definition: IOobject.C:263
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
const Cmpt & xx() const
Definition: TensorI.H:160
virtual tmp< tensorField > transformTensor(const tensorField &st) const =0
Transform tensor field using transformation tensorField.
const dimensionSet dimForce
defineTypeNameAndDebug(DarcyForchheimer, 0)
static const Vector zero
Definition: Vector.H:80