Maxwell.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) 2016-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 
26 #include "Maxwell.H"
27 #include "fvModels.H"
28 #include "fvConstraints.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 namespace laminarModels
36 {
37 
38 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
39 
40 template<class BasicMomentumTransportModel>
41 PtrList<dimensionedScalar>
43 (
44  const word& name,
45  const dimensionSet& dims
46 ) const
47 {
48  PtrList<dimensionedScalar> modeCoeffs(nModes_);
49 
50  if (modeCoefficients_.size())
51  {
52  if (this->coeffDict().found(name))
53  {
54  IOWarningInFunction(this->coeffDict())
55  << "Using 'modes' list, '" << name << "' entry will be ignored."
56  << endl;
57  }
58 
59  forAll(modeCoefficients_, modei)
60  {
61  modeCoeffs.set
62  (
63  modei,
65  (
66  name,
67  dims,
68  modeCoefficients_[modei].lookup(name)
69  )
70  );
71  }
72  }
73  else
74  {
75  modeCoeffs.set
76  (
77  0,
79  (
80  name,
81  dims,
82  this->coeffDict_.lookup(name)
83  )
84  );
85  }
86 
87  return modeCoeffs;
88 }
89 
90 
91 template<class BasicMomentumTransportModel>
93 (
94  const label modei,
95  volSymmTensorField& sigma
96 ) const
97 {
98  return -fvm::Sp(this->alpha_*this->rho_/lambdas_[modei], sigma);
99 }
100 
101 
102 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
103 
104 template<class BasicMomentumTransportModel>
106 (
107  const alphaField& alpha,
108  const rhoField& rho,
109  const volVectorField& U,
110  const surfaceScalarField& alphaRhoPhi,
111  const surfaceScalarField& phi,
112  const transportModel& transport,
113  const word& type
114 )
115 :
117  (
118  type,
119  alpha,
120  rho,
121  U,
122  alphaRhoPhi,
123  phi,
124  transport
125  ),
126 
127  modeCoefficients_
128  (
129  this->coeffDict().found("modes")
131  (
132  this->coeffDict().lookup("modes")
133  )
135  ),
136 
137  nModes_(modeCoefficients_.size() ? modeCoefficients_.size() : 1),
138 
139  nuM_("nuM", dimViscosity, this->coeffDict_.lookup("nuM")),
140 
141  lambdas_(readModeCoefficients("lambda", dimTime)),
142 
143  sigma_
144  (
145  IOobject
146  (
147  IOobject::groupName("sigma", alphaRhoPhi.group()),
148  this->runTime_.timeName(),
149  this->mesh_,
152  ),
153  this->mesh_
154  )
155 {
156  if (nModes_ > 1)
157  {
158  sigmas_.setSize(nModes_);
159 
160  forAll(sigmas_, modei)
161  {
162  IOobject header
163  (
164  IOobject::groupName("sigma" + name(modei), alphaRhoPhi.group()),
165  this->runTime_.timeName(),
166  this->mesh_,
168  );
169 
170  // Check if mode field exists and can be read
171  if (header.typeHeaderOk<volSymmTensorField>(true))
172  {
173  Info<< " Reading mode stress field "
174  << header.name() << endl;
175 
176  sigmas_.set
177  (
178  modei,
180  (
181  IOobject
182  (
183  header.name(),
184  this->runTime_.timeName(),
185  this->mesh_,
188  ),
189  this->mesh_
190  )
191  );
192  }
193  else
194  {
195  sigmas_.set
196  (
197  modei,
199  (
200  IOobject
201  (
202  header.name(),
203  this->runTime_.timeName(),
204  this->mesh_,
207  ),
208  sigma_
209  )
210  );
211  }
212  }
213  }
214 
215  if (type == typeName)
216  {
217  this->printCoeffs(type);
218  }
219 }
220 
221 
222 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
223 
224 template<class BasicMomentumTransportModel>
226 {
228  {
229  if (modeCoefficients_.size())
230  {
231  this->coeffDict().lookup("modes") >> modeCoefficients_;
232  }
233 
234  nuM_.read(this->coeffDict());
235 
236  lambdas_ = readModeCoefficients("lambda", dimTime);
237 
238  return true;
239  }
240  else
241  {
242  return false;
243  }
244 }
245 
246 
247 template<class BasicMomentumTransportModel>
249 {
250  return volScalarField::New
251  (
252  IOobject::groupName("nuEff", this->alphaRhoPhi_.group()),
253  this->nu()
254  );
255 }
256 
257 
258 template<class BasicMomentumTransportModel>
260 (
261  const label patchi
262 ) const
263 {
264  return this->nu(patchi);
265 }
266 
267 
268 template<class BasicMomentumTransportModel>
270 {
271  return sigma_;
272 }
273 
274 
275 template<class BasicMomentumTransportModel>
277 {
279  (
280  IOobject::groupName("devTau", this->alphaRhoPhi_.group()),
281  this->alpha_*this->rho_*sigma_
282  - (this->alpha_*this->rho_*this->nu())
283  *dev(twoSymm(fvc::grad(this->U_)))
284  );
285 }
286 
287 
288 template<class BasicMomentumTransportModel>
290 (
291  volVectorField& U
292 ) const
293 {
294  return
295  (
296  fvc::div
297  (
298  this->alpha_*this->rho_*this->nuM_*fvc::grad(U)
299  )
300  + fvc::div(this->alpha_*this->rho_*sigma_)
301  - fvc::div(this->alpha_*this->rho_*this->nu()*dev2(T(fvc::grad(U))))
302  - fvm::laplacian(this->alpha_*this->rho_*nu0(), U)
303  );
304 }
305 
306 
307 template<class BasicMomentumTransportModel>
310 (
311  const volScalarField& rho,
312  volVectorField& U
313 ) const
314 {
315  return
316  (
317  fvc::div
318  (
319  this->alpha_*rho*this->nuM_*fvc::grad(U)
320  )
321  + fvc::div(this->alpha_*rho*sigma_)
322  - fvc::div(this->alpha_*rho*this->nu()*dev2(T(fvc::grad(U))))
323  - fvm::laplacian(this->alpha_*rho*nu0(), U)
324  );
325 }
326 
327 
328 template<class BasicMomentumTransportModel>
330 {
331  // Local references
332  const alphaField& alpha = this->alpha_;
333  const rhoField& rho = this->rho_;
334  const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
335  const volVectorField& U = this->U_;
336  const Foam::fvModels& fvModels(Foam::fvModels::New(this->mesh_));
338  (
339  Foam::fvConstraints::New(this->mesh_)
340  );
341 
343 
344  tmp<volTensorField> tgradU(fvc::grad(U));
345  const volTensorField& gradU = tgradU();
346 
347  forAll(lambdas_, modei)
348  {
349  volSymmTensorField& sigma = nModes_ == 1 ? sigma_ : sigmas_[modei];
350 
352  (
353  IOobject
354  (
356  (
357  "rLambda"
358  + (nModes_ == 1 ? word::null : name(modei)),
359  this->alphaRhoPhi_.group()
360  ),
361  this->runTime_.constant(),
362  this->mesh_
363  ),
364  1/lambdas_[modei]
365  );
366 
367  // Note sigma is positive on lhs of momentum eqn
368  const volSymmTensorField P
369  (
370  twoSymm(sigma & gradU)
371  - nuM_*rLambda*twoSymm(gradU)
372  );
373 
374  // Viscoelastic stress equation
375  fvSymmTensorMatrix sigmaEqn
376  (
377  fvm::ddt(alpha, rho, sigma)
378  + fvm::div
379  (
380  alphaRhoPhi,
381  sigma,
382  "div(" + alphaRhoPhi.name() + ',' + sigma_.name() + ')'
383  )
384  ==
385  alpha*rho*P
386  + sigmaSource(modei, sigma)
387  + fvModels.source(alpha, rho, sigma)
388  );
389 
390  sigmaEqn.relax();
391  fvConstraints.constrain(sigmaEqn);
392  sigmaEqn.solve("sigma");
393  fvConstraints.constrain(sigma);
394  }
395 
396  if (sigmas_.size())
397  {
398  volSymmTensorField sigmaSum("sigmaSum", sigmas_[0]);
399 
400  for (label modei = 1; modei<sigmas_.size(); modei++)
401  {
402  sigmaSum += sigmas_[modei];
403  }
404 
405  sigma_ == sigmaSum;
406  }
407 }
408 
409 
410 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
411 
412 } // End namespace laminarModels
413 } // End namespace Foam
414 
415 // ************************************************************************* //
const dimensionSet dimViscosity
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
Templated abstract base class for laminar transport models.
Definition: laminarModel.H:49
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
const word & name() const
Return name.
Definition: IOobject.H:303
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
virtual tmp< volSymmTensorField > devTau() const
Return the effective stress tensor.
Definition: Maxwell.C:276
Maxwell(const alphaField &alpha, const rhoField &rho, const volVectorField &U, const surfaceScalarField &alphaRhoPhi, const surfaceScalarField &phi, const transportModel &transport, const word &type=typeName)
Construct from components.
Definition: Maxwell.C:106
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
BasicMomentumTransportModel::rhoField rhoField
Definition: laminarModel.H:77
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
BasicMomentumTransportModel::alphaField alphaField
Definition: laminarModel.H:76
PtrList< dimensionedScalar > readModeCoefficients(const word &name, const dimensionSet &dims) const
Definition: Maxwell.C:43
BasicMomentumTransportModel::transportModel transportModel
Definition: laminarModel.H:78
virtual tmp< volScalarField > nuEff() const
Return the effective viscosity, i.e. the laminar viscosity.
Definition: Maxwell.C:248
dimensionedSymmTensor twoSymm(const dimensionedSymmTensor &dt)
const dimensionSet dimTime
Dimension set for the base types.
Definition: dimensionSet.H:120
stressControl lookup("compactNormalStress") >> compactNormalStress
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
Foam::fvConstraints & fvConstraints
virtual tmp< fvSymmTensorMatrix > sigmaSource(const label modei, volSymmTensorField &sigma) const
Definition: Maxwell.C:93
A class for handling words, derived from string.
Definition: word.H:59
static word groupName(Name name, const word &group)
Finite volume models.
Definition: fvModels.H:60
virtual tmp< fvVectorMatrix > divDevTau(volVectorField &U) const
Return the source term for the momentum equation.
Definition: Maxwell.C:290
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
static const word null
An empty word.
Definition: word.H:77
phi
Definition: correctPhi.H:3
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
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
static autoPtr< dictionary > New(Istream &)
Construct top-level dictionary on freestore from Istream.
Definition: dictionaryIO.C:96
virtual tmp< volSymmTensorField > sigma() const
Return the stress tensor [m^2/s^2].
Definition: Maxwell.C:269
virtual bool read()
Read model coefficients if they have changed.
Definition: Maxwell.C:225
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:521
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
bool constrain(fvMatrix< Type > &eqn) const
Apply constraints to an equation.
Foam::fvModels & fvModels
virtual void correct()
Solve the turbulence equations and correct eddy-Viscosity and.
Definition: Maxwell.C:329
dimensionedSymmTensor dev2(const dimensionedSymmTensor &dt)
U
Definition: pEqn.H:72
Finite volume constraints.
Definition: fvConstraints.H:57
virtual void correct()
Correct the laminar transport.
Definition: laminarModel.C:247
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:46
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
messageStream Info
#define IOWarningInFunction(ios)
Report an IO warning using Foam::Warning.
A class for managing temporary objects.
Definition: PtrList.H:53
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Namespace for OpenFOAM.