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-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 "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
179  (
180  UIndirectList<vector>(mesh_.C(), cells)()
181  );
182 
183  alpha_[zoneI] = R.transformTensor(alpha_[zoneI]);
184  beta_[zoneI] = R.transformTensor(beta_[zoneI]);
185  }
186  }
187 }
188 
189 
191 (
192  const volVectorField& U,
193  const volScalarField& rho,
194  const volScalarField& mu,
195  vectorField& force
196 ) const
197 {
198  scalarField Udiag(U.size(), 0.0);
199  vectorField Usource(U.size(), Zero);
200  const scalarField& V = mesh_.V();
201  scalar rhoRef = coeffs_.lookup<scalar>("rhoRef");
202 
203  apply(Udiag, Usource, V, U, rhoRef);
204 
205  force = Udiag*U - Usource;
206 }
207 
208 
210 (
211  fvVectorMatrix& UEqn
212 ) const
213 {
214  const vectorField& U = UEqn.psi();
215  const scalarField& V = mesh_.V();
216  scalarField& Udiag = UEqn.diag();
217  vectorField& Usource = UEqn.source();
218 
219  scalar rho = 1.0;
220  if (UEqn.dimensions() == dimForce)
221  {
222  coeffs_.lookup("rhoRef") >> rho;
223  }
224 
225  apply(Udiag, Usource, V, U, rho);
226 }
227 
228 
230 (
231  const fvVectorMatrix& UEqn,
232  volTensorField& AU
233 ) const
234 {
235  const vectorField& U = UEqn.psi();
236 
237  scalar rho = 1.0;
238  if (UEqn.dimensions() == dimForce)
239  {
240  coeffs_.lookup("rhoRef") >> rho;
241  }
242 
243  apply(AU, U, rho);
244 }
245 
246 
248 {
249  os << indent << name_ << endl;
250  dict_.write(os);
251 
252  return true;
253 }
254 
255 
256 // ************************************************************************* //
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:156
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
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual ~fixedCoeff()
Destructor.
Definition: fixedCoeff.C:128
const dimensionSet dimless
Macros for easy insertion into run-time selection tables.
const dimensionSet dimLength
const dimensionSet dimTime
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)
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
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
virtual void correct(fvVectorMatrix &UEqn) const
Add resistance.
Definition: fixedCoeff.C:210
#define R(A, B, C, D, E, F, K, M)
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
A List with indirect addressing.
Definition: fvMatrix.H:106
bool writeData(Ostream &os) const
Write.
Definition: fixedCoeff.C:247
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:191
virtual tensor transformTensor(const tensor &st) const =0
Transform tensor using transformation tensorField.
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
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:54
Namespace for OpenFOAM.
const dimensionSet & dimensions() const
Definition: fvMatrix.H:287