Implicit.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) 2013-2022 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 "Implicit.H"
28 #include "fvmDdt.H"
29 #include "fvmDiv.H"
30 #include "fvmLaplacian.H"
31 #include "fvcReconstruct.H"
32 #include "volPointInterpolation.H"
34 
35 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
36 
37 template<class CloudType>
39 (
40  const dictionary& dict,
41  CloudType& owner
42 )
43 :
44  PackingModel<CloudType>(dict, owner, typeName),
45  alpha_
46  (
47  IOobject
48  (
49  this->owner().name() + ":alpha",
50  this->owner().db().time().timeName(),
51  this->owner().mesh(),
52  IOobject::NO_READ,
53  IOobject::NO_WRITE
54  ),
55  this->owner().mesh(),
57  zeroGradientFvPatchScalarField::typeName
58  ),
59  phiCorrect_(nullptr),
60  uCorrect_(nullptr),
61  applyLimiting_(this->coeffDict().lookup("applyLimiting")),
62  applyGravity_(this->coeffDict().lookup("applyGravity")),
63  alphaMin_(this->coeffDict().template lookup<scalar>("alphaMin")),
64  rhoMin_(this->coeffDict().template lookup<scalar>("rhoMin"))
65 {
66  alpha_ = this->owner().theta();
67  alpha_.oldTime();
68 }
69 
70 
71 template<class CloudType>
73 (
74  const Implicit<CloudType>& cm
75 )
76 :
78  alpha_(cm.alpha_),
79  phiCorrect_(cm.phiCorrect_()),
80  uCorrect_(cm.uCorrect_()),
81  applyLimiting_(cm.applyLimiting_),
82  applyGravity_(cm.applyGravity_),
83  alphaMin_(cm.alphaMin_),
84  rhoMin_(cm.rhoMin_)
85 {
86  alpha_.oldTime();
87 }
88 
89 
90 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
91 
92 template<class CloudType>
94 {}
95 
96 
97 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
98 
99 template<class CloudType>
101 {
103 
104  if (store)
105  {
106  const fvMesh& mesh = this->owner().mesh();
107  const dimensionedScalar deltaT = this->owner().db().time().deltaT();
108  const word& cloudName = this->owner().name();
109 
110  const dimensionedVector& g = this->owner().g();
111  const volScalarField& rhoc = this->owner().rho();
112 
113  const AveragingMethod<scalar>& rhoAverage =
115  (
116  cloudName + ":rhoAverage"
117  );
118  const AveragingMethod<vector>& uAverage =
120  (
121  cloudName + ":uAverage"
122  );
123  const AveragingMethod<scalar>& uSqrAverage =
125  (
126  cloudName + ":uSqrAverage"
127  );
128 
129  mesh.schemes().setFluxRequired(alpha_.name());
130 
131  // Property fields
132  // ~~~~~~~~~~~~~~~
133 
134  // volume fraction field
135  alpha_ = max(this->owner().theta(), alphaMin_);
136  alpha_.correctBoundaryConditions();
137 
138  // average density
140  (
141  IOobject
142  (
143  cloudName + ":rho",
144  this->owner().db().time().timeName(),
145  mesh,
146  IOobject::NO_READ,
147  IOobject::NO_WRITE
148  ),
149  mesh,
152  );
153  rho.primitiveFieldRef() = max(rhoAverage.primitiveField(), rhoMin_);
155 
156  // Stress field
157  // ~~~~~~~~~~~~
158 
159  // stress derivative wrt volume fraction
160  volScalarField tauPrime
161  (
162  IOobject
163  (
164  cloudName + ":tauPrime",
165  this->owner().db().time().timeName(),
166  mesh,
167  IOobject::NO_READ,
168  IOobject::NO_WRITE
169  ),
170  mesh,
173  );
174 
175  tauPrime.primitiveFieldRef() =
176  this->particleStressModel_->dTaudTheta
177  (
178  alpha_.primitiveField(),
179  rho.primitiveField(),
180  uSqrAverage.primitiveField()
181  )();
182 
183  tauPrime.correctBoundaryConditions();
184 
185 
186  // Gravity flux
187  // ~~~~~~~~~~~~
188 
189  tmp<surfaceScalarField> phiGByA;
190 
191  if (applyGravity_)
192  (
193  phiGByA = surfaceScalarField::New
194  (
195  "phiGByA",
196  deltaT*(g & mesh.Sf())*fvc::interpolate(1.0 - rhoc/rho)
197  )
198  );
199 
200 
201  // Implicit solution for the volume fraction
202  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
203 
205  tauPrimeByRhoAf
206  (
207  "tauPrimeByRhoAf",
208  fvc::interpolate(deltaT*tauPrime/rho)
209  );
210 
211  fvScalarMatrix alphaEqn
212  (
213  fvm::ddt(alpha_)
214  - fvc::ddt(alpha_)
215  - fvm::laplacian(tauPrimeByRhoAf, alpha_)
216  );
217 
218  if (applyGravity_)
219  {
220  alphaEqn += fvm::div(phiGByA(), alpha_);
221  }
222 
223  alphaEqn.solve();
224 
225 
226  // Generate correction fields
227  // ~~~~~~~~~~~~~~~~~
228 
229  // correction volumetric flux
230  phiCorrect_ = surfaceScalarField::New
231  (
232  cloudName + ":phiCorrect",
233  alphaEqn.flux()/fvc::interpolate(alpha_)
234  );
235 
236  // limit the correction flux
237  if (applyLimiting_)
238  {
240  (
241  IOobject
242  (
243  cloudName + ":U",
244  this->owner().db().time().timeName(),
245  mesh,
246  IOobject::NO_READ,
247  IOobject::NO_WRITE
248  ),
249  mesh,
252  );
253  U.primitiveFieldRef() = uAverage.primitiveField();
255 
257  (
258  cloudName + ":phi",
259  linearInterpolate(U) & mesh.Sf()
260  );
261 
262  if (applyGravity_)
263  {
264  phiCorrect_.ref() -= phiGByA();
265  }
266 
267  forAll(phiCorrect_(), facei)
268  {
269  // Current and correction fluxes
270  const scalar phiCurr = phi[facei];
271  scalar& phiCorr = phiCorrect_.ref()[facei];
272 
273  // Don't limit if the correction is in the opposite direction to
274  // the flux. We need all the help we can get in this state.
275  if (phiCurr*phiCorr < 0)
276  {}
277 
278  // If the correction and the flux are in the same direction then
279  // don't apply any more correction than is already present in
280  // the flux.
281  else if (phiCorr > 0)
282  {
283  phiCorr = max(phiCorr - phiCurr, 0);
284  }
285  else
286  {
287  phiCorr = min(phiCorr - phiCurr, 0);
288  }
289  }
290 
291  if (applyGravity_)
292  {
293  phiCorrect_.ref() += phiGByA();
294  }
295  }
296 
297  // correction velocity
298  uCorrect_ = volVectorField::New
299  (
300  cloudName + ":uCorrect",
301  fvc::reconstruct(phiCorrect_())
302  );
303  uCorrect_->correctBoundaryConditions();
304 
305  // Info << endl;
306  // Info << " alpha: " << alpha_.primitiveField() << endl;
307  // Info << "phiCorrect: " << phiCorrect_->primitiveField() << endl;
308  // Info << " uCorrect: " << uCorrect_->primitiveField() << endl;
309  // Info << endl;
310  }
311  else
312  {
313  alpha_.oldTime();
314  phiCorrect_.clear();
315  uCorrect_.clear();
316  }
317 }
318 
319 
320 template<class CloudType>
322 (
323  typename CloudType::parcelType& p,
324  const scalar deltaT
325 ) const
326 {
327  const fvMesh& mesh = this->owner().mesh();
328 
329  // containing tetrahedron and parcel coordinates within
330  const label celli = p.cell();
331  const label facei = p.tetFace();
332 
333  // cell velocity
334  const vector U = uCorrect_()[celli];
335 
336  // face geometry
337  vector nHat = mesh.faces()[facei].area(mesh.points());
338  const scalar nMag = mag(nHat);
339  nHat /= nMag;
340 
341  // get face flux
342  scalar phi;
343  const label patchi = mesh.boundaryMesh().whichPatch(facei);
344  if (patchi == -1)
345  {
346  phi = phiCorrect_()[facei];
347  }
348  else
349  {
350  phi =
351  phiCorrect_().boundaryField()[patchi]
352  [
353  mesh.boundaryMesh()[patchi].whichFace(facei)
354  ];
355  }
356 
357  // interpolant equal to 1 at the cell centre and 0 at the face
358  const scalar t = p.coordinates()[0];
359 
360  // the normal component of the velocity correction is interpolated linearly
361  // the tangential component is equal to that at the cell centre
362  return U + (1.0 - t)*nHat*(phi/nMag - (U & nHat));
363 }
364 
365 
366 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
Implicit model for applying an inter-particle stress to the particles.
Definition: Implicit.H:56
This boundary condition supplies a fixed value constraint, and is the base class for a number of othe...
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
virtual tmp< Field< Type > > primitiveField() const =0
Return an internal field of the average.
virtual void cacheFields(const bool store)
Calculate the inter particles stresses.
Definition: Implicit.C:100
const surfaceVectorField & Sf() const
Return cell face area vectors.
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
autoPtr< CompressibleMomentumTransportModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const viscosity &viscosity)
const fvSchemes & schemes() const
Return the fvSchemes.
Definition: fvMesh.C:1672
U
Definition: pEqn.H:72
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > linearInterpolate(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: linear.H:108
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
const dimensionSet dimPressure
Reconstruct volField from a face flux field.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Base class for packing models.
Definition: MPPICCloud.H:53
virtual ~Implicit()
Destructor.
Definition: Implicit.C:93
Calculate the matrix for the laplacian of the field.
void setFluxRequired(const word &name) const
Definition: fvSchemes.C:487
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
const dimensionSet dimless
fvMesh & mesh
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1211
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
stressControl lookup("compactNormalStress") >> compactNormalStress
Implicit(const dictionary &dict, CloudType &owner)
Construct from components.
Definition: Implicit.C:39
A class for handling words, derived from string.
Definition: word.H:59
const dimensionSet dimDensity
Calculate the matrix for the first temporal derivative.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
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
word timeName
Definition: getTimeIndex.H:3
static const zero Zero
Definition: zero.H:97
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1224
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:221
volScalarField rhoc(IOobject(rhocValue.name(), runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE), mesh, rhocValue)
Internal::FieldType & primitiveFieldRef()
Return a reference to the internal field.
SolverPerformance< Type > solve(const dictionary &)
Solve segregated or coupled returning the solution statistics.
Definition: fvMatrixSolve.C:58
const dimensionSet dimVelocity
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
This boundary condition applies a zero-gradient condition from the patch internal field onto the patc...
Calculate the matrix for the divergence of the given field and flux.
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
label patchi
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:95
void correctBoundaryConditions()
Correct boundary field.
dimensioned< scalar > mag(const dimensioned< Type > &)
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > flux() const
Return the face-flux field from the matrix.
Definition: fvMatrix.C:927
const word cloudName(propsDict.lookup("cloudName"))
virtual vector velocityCorrection(typename CloudType::parcelType &p, const scalar deltaT) const
Calculate the velocity correction.
Definition: Implicit.C:322
A class for managing temporary objects.
Definition: PtrList.H:53
const dimensionedVector & g
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:98
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
label whichPatch(const label faceIndex) const
Return patch index for a given face label.