fixedCoeff.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-2019 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 "fixedCoeff.H"
28 #include "fvMatrices.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34  namespace porosityModels
35  {
36  defineTypeNameAndDebug(fixedCoeff, 0);
37  addToRunTimeSelectionTable(porosityModel, fixedCoeff, mesh);
38  }
39 }
40 
41 
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
43 
44 void Foam::porosityModels::fixedCoeff::apply
45 (
46  scalarField& Udiag,
47  vectorField& Usource,
48  const scalarField& V,
49  const vectorField& U,
50  const scalar rho
51 ) const
52 {
53  forAll(cellZoneIDs_, zoneI)
54  {
55  const tensorField& alphaZones = alpha_[zoneI];
56  const tensorField& betaZones = beta_[zoneI];
57 
58  const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
59 
60  forAll(cells, i)
61  {
62  const label celli = cells[i];
63  const label j = fieldIndex(i);
64  const tensor Cd = rho*(alphaZones[j] + betaZones[j]*mag(U[celli]));
65  const scalar isoCd = tr(Cd);
66 
67  Udiag[celli] += V[celli]*isoCd;
68  Usource[celli] -= V[celli]*((Cd - I*isoCd) & U[celli]);
69  }
70  }
71 }
72 
73 
74 void Foam::porosityModels::fixedCoeff::apply
75 (
76  tensorField& AU,
77  const vectorField& U,
78  const scalar rho
79 ) const
80 {
81 
82  forAll(cellZoneIDs_, zoneI)
83  {
84  const tensorField& alphaZones = alpha_[zoneI];
85  const tensorField& betaZones = beta_[zoneI];
86 
87  const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
88 
89  forAll(cells, i)
90  {
91  const label celli = cells[i];
92  const label j = fieldIndex(i);
93  const tensor alpha = alphaZones[j];
94  const tensor beta = betaZones[j];
95 
96  AU[celli] += rho*(alpha + beta*mag(U[celli]));
97  }
98  }
99 }
100 
101 
102 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
103 
105 (
106  const word& name,
107  const word& modelType,
108  const fvMesh& mesh,
109  const dictionary& dict,
110  const word& cellZoneName
111 )
112 :
113  porosityModel(name, modelType, mesh, dict, cellZoneName),
114  alphaXYZ_("alpha", dimless/dimTime, coeffs_),
115  betaXYZ_("beta", dimless/dimLength, coeffs_),
116  alpha_(cellZoneIDs_.size()),
117  beta_(cellZoneIDs_.size())
118 {
119  adjustNegativeResistance(alphaXYZ_);
120  adjustNegativeResistance(betaXYZ_);
121 
122  calcTransformModelData();
123 }
124 
125 
126 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
127 
129 {}
130 
131 
132 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
133 
135 {
136  if (coordSys_.R().uniform())
137  {
138  forAll(cellZoneIDs_, zoneI)
139  {
140  alpha_[zoneI].setSize(1);
141  beta_[zoneI].setSize(1);
142 
143  alpha_[zoneI][0] = Zero;
144  alpha_[zoneI][0].xx() = alphaXYZ_.value().x();
145  alpha_[zoneI][0].yy() = alphaXYZ_.value().y();
146  alpha_[zoneI][0].zz() = alphaXYZ_.value().z();
147  alpha_[zoneI][0] = coordSys_.R().transformTensor(alpha_[zoneI][0]);
148 
149  beta_[zoneI][0] = Zero;
150  beta_[zoneI][0].xx() = betaXYZ_.value().x();
151  beta_[zoneI][0].yy() = betaXYZ_.value().y();
152  beta_[zoneI][0].zz() = betaXYZ_.value().z();
153  beta_[zoneI][0] = coordSys_.R().transformTensor(beta_[zoneI][0]);
154  }
155  }
156  else
157  {
158  forAll(cellZoneIDs_, zoneI)
159  {
160  const labelList& cells = mesh_.cellZones()[cellZoneIDs_[zoneI]];
161 
162  alpha_[zoneI].setSize(cells.size());
163  beta_[zoneI].setSize(cells.size());
164 
165  forAll(cells, i)
166  {
167  alpha_[zoneI][i] = Zero;
168  alpha_[zoneI][i].xx() = alphaXYZ_.value().x();
169  alpha_[zoneI][i].yy() = alphaXYZ_.value().y();
170  alpha_[zoneI][i].zz() = alphaXYZ_.value().z();
171 
172  beta_[zoneI][i] = Zero;
173  beta_[zoneI][i].xx() = betaXYZ_.value().x();
174  beta_[zoneI][i].yy() = betaXYZ_.value().y();
175  beta_[zoneI][i].zz() = betaXYZ_.value().z();
176  }
177 
178  const coordinateRotation& R = coordSys_.R(mesh_, cells);
179 
180  alpha_[zoneI] = R.transformTensor(alpha_[zoneI], cells);
181  beta_[zoneI] = R.transformTensor(beta_[zoneI], cells);
182  }
183  }
184 }
185 
186 
188 (
189  const volVectorField& U,
190  const volScalarField& rho,
191  const volScalarField& mu,
192  vectorField& force
193 ) const
194 {
195  scalarField Udiag(U.size(), 0.0);
196  vectorField Usource(U.size(), Zero);
197  const scalarField& V = mesh_.V();
198  scalar rhoRef = coeffs_.lookup<scalar>("rhoRef");
199 
200  apply(Udiag, Usource, V, U, rhoRef);
201 
202  force = Udiag*U - Usource;
203 }
204 
205 
207 (
208  fvVectorMatrix& UEqn
209 ) const
210 {
211  const vectorField& U = UEqn.psi();
212  const scalarField& V = mesh_.V();
213  scalarField& Udiag = UEqn.diag();
214  vectorField& Usource = UEqn.source();
215 
216  scalar rho = 1.0;
217  if (UEqn.dimensions() == dimForce)
218  {
219  coeffs_.lookup("rhoRef") >> rho;
220  }
221 
222  apply(Udiag, Usource, V, U, rho);
223 }
224 
225 
227 (
228  fvVectorMatrix& UEqn,
229  const volScalarField&,
230  const volScalarField&
231 ) const
232 {
233  const vectorField& U = UEqn.psi();
234  const scalarField& V = mesh_.V();
235  scalarField& Udiag = UEqn.diag();
236  vectorField& Usource = UEqn.source();
237 
238  scalar rho = 1.0;
239  if (UEqn.dimensions() == dimForce)
240  {
241  coeffs_.lookup("rhoRef") >> rho;
242  }
243 
244  apply(Udiag, Usource, V, U, rho);
245 }
246 
247 
249 (
250  const fvVectorMatrix& UEqn,
251  volTensorField& AU
252 ) const
253 {
254  const vectorField& U = UEqn.psi();
255 
256  scalar rho = 1.0;
257  if (UEqn.dimensions() == dimForce)
258  {
259  coeffs_.lookup("rhoRef") >> rho;
260  }
261 
262  apply(AU, U, rho);
263 }
264 
265 
267 {
268  os << indent << name_ << endl;
269  dict_.write(os);
270 
271  return true;
272 }
273 
274 
275 // ************************************************************************* //
fixedCoeff(const word &name, const word &modelType, const fvMesh &mesh, const dictionary &dict, const word &cellZoneName)
Definition: fixedCoeff.C:105
Abstract base class for coordinate rotation.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual Ostream & write(const char)=0
Write character.
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
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:158
addToRunTimeSelectionTable(porosityModel, powerLawLopesdaCosta, mesh)
const GeometricField< Type, fvPatchField, volMesh > & psi() const
Definition: fvMatrix.H:282
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
virtual ~fixedCoeff()
Destructor.
Definition: fixedCoeff.C:128
Macros for easy insertion into run-time selection tables.
const cellShapeList & cells
static const Identity< scalar > I
Definition: Identity.H:93
A class for handling words, derived from string.
Definition: word.H:59
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual void calcTransformModelData()
Transform the model data wrt mesh changes.
Definition: fixedCoeff.C:134
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
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
dimensionedScalar tr(const dimensionedSphericalTensor &dt)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
const dimensionSet dimForce
Field< Type > & source()
Definition: fvMatrix.H:292
virtual tmp< tensorField > transformTensor(const tensorField &st) const =0
Transform tensor field using transformation tensorField.
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
virtual void correct(fvVectorMatrix &UEqn) const
Add resistance.
Definition: fixedCoeff.C:207
#define R(A, B, C, D, E, F, K, M)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
bool writeData(Ostream &os) const
Write.
Definition: fixedCoeff.C:266
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
A special matrix type and solver, designed for finite volume solutions of scalar equations.
virtual void calcForce(const volVectorField &U, const volScalarField &rho, const volScalarField &mu, vectorField &force) const
Calculate the porosity force.
Definition: fixedCoeff.C:188
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
scalarField & diag()
Definition: lduMatrix.C:186
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
Top level model for porosity models.
Definition: porosityModel.H:55
Namespace for OpenFOAM.
const dimensionSet & dimensions() const
Definition: fvMatrix.H:287