Implicit.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) 2013-2017 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(),
56  dimensionedScalar("zero", dimless, 0.0),
57  zeroGradientFvPatchScalarField::typeName
58  ),
59  phiCorrect_(nullptr),
60  uCorrect_(nullptr),
61  applyLimiting_(this->coeffDict().lookup("applyLimiting")),
62  applyGravity_(this->coeffDict().lookup("applyGravity")),
63  alphaMin_(readScalar(this->coeffDict().lookup("alphaMin"))),
64  rhoMin_(readScalar(this->coeffDict().lookup("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.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,
150  dimensionedScalar("zero", dimDensity, 0),
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,
171  dimensionedScalar("zero", dimPressure, 0),
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 = tmp<surfaceScalarField>
194  (
196  (
197  "phiGByA",
198  deltaT*(g & mesh.Sf())*fvc::interpolate(1.0 - rhoc/rho)
199  )
200  )
201  );
202 
203 
204  // Implicit solution for the volume fraction
205  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
206 
208  tauPrimeByRhoAf
209  (
210  "tauPrimeByRhoAf",
211  fvc::interpolate(deltaT*tauPrime/rho)
212  );
213 
214  fvScalarMatrix alphaEqn
215  (
216  fvm::ddt(alpha_)
217  - fvc::ddt(alpha_)
218  - fvm::laplacian(tauPrimeByRhoAf, alpha_)
219  );
220 
221  if (applyGravity_)
222  {
223  alphaEqn += fvm::div(phiGByA(), alpha_);
224  }
225 
226  alphaEqn.solve();
227 
228 
229  // Generate correction fields
230  // ~~~~~~~~~~~~~~~~~
231 
232  // correction volumetric flux
233  phiCorrect_ = tmp<surfaceScalarField>
234  (
236  (
237  cloudName + ":phiCorrect",
238  alphaEqn.flux()/fvc::interpolate(alpha_)
239  )
240  );
241 
242  // limit the correction flux
243  if (applyLimiting_)
244  {
246  (
247  IOobject
248  (
249  cloudName + ":U",
250  this->owner().db().time().timeName(),
251  mesh,
252  IOobject::NO_READ,
253  IOobject::NO_WRITE
254  ),
255  mesh,
258  );
259  U.primitiveFieldRef() = uAverage.primitiveField();
261 
263  (
264  cloudName + ":phi",
265  linearInterpolate(U) & mesh.Sf()
266  );
267 
268  if (applyGravity_)
269  {
270  phiCorrect_.ref() -= phiGByA();
271  }
272 
273  forAll(phiCorrect_(), facei)
274  {
275  // Current and correction fluxes
276  const scalar phiCurr = phi[facei];
277  scalar& phiCorr = phiCorrect_.ref()[facei];
278 
279  // Don't limit if the correction is in the opposite direction to
280  // the flux. We need all the help we can get in this state.
281  if (phiCurr*phiCorr < 0)
282  {}
283 
284  // If the correction and the flux are in the same direction then
285  // don't apply any more correction than is already present in
286  // the flux.
287  else if (phiCorr > 0)
288  {
289  phiCorr = max(phiCorr - phiCurr, 0);
290  }
291  else
292  {
293  phiCorr = min(phiCorr - phiCurr, 0);
294  }
295  }
296 
297  if (applyGravity_)
298  {
299  phiCorrect_.ref() += phiGByA();
300  }
301  }
302 
303  // correction velocity
304  uCorrect_ = tmp<volVectorField>
305  (
306  new volVectorField
307  (
308  cloudName + ":uCorrect",
309  fvc::reconstruct(phiCorrect_())
310  )
311 
312  );
313  uCorrect_->correctBoundaryConditions();
314 
315  //Info << endl;
316  //Info << " alpha: " << alpha_.primitiveField() << endl;
317  //Info << "phiCorrect: " << phiCorrect_->primitiveField() << endl;
318  //Info << " uCorrect: " << uCorrect_->primitiveField() << endl;
319  //Info << endl;
320  }
321  else
322  {
323  alpha_.oldTime();
324  phiCorrect_.clear();
325  uCorrect_.clear();
326  }
327 }
328 
329 
330 template<class CloudType>
332 (
333  typename CloudType::parcelType& p,
334  const scalar deltaT
335 ) const
336 {
337  const fvMesh& mesh = this->owner().mesh();
338 
339  // containing tetrahedron and parcel coordinates within
340  const label celli = p.cell();
341  const label facei = p.tetFace();
342 
343  // cell velocity
344  const vector U = uCorrect_()[celli];
345 
346  // face geometry
347  vector nHat = mesh.faces()[facei].normal(mesh.points());
348  const scalar nMag = mag(nHat);
349  nHat /= nMag;
350 
351  // get face flux
352  scalar phi;
353  const label patchi = mesh.boundaryMesh().whichPatch(facei);
354  if (patchi == -1)
355  {
356  phi = phiCorrect_()[facei];
357  }
358  else
359  {
360  phi =
361  phiCorrect_().boundaryField()[patchi]
362  [
363  mesh.boundaryMesh()[patchi].whichFace(facei)
364  ];
365  }
366 
367  // interpolant equal to 1 at the cell centre and 0 at the face
368  const scalar t = p.coordinates()[0];
369 
370  // the normal component of the velocity correction is interpolated linearly
371  // the tangential component is equal to that at the cell centre
372  return U + (1.0 - t)*nHat*(phi/nMag - (U & nHat));
373 }
374 
375 
376 // ************************************************************************* //
#define readScalar
Definition: doubleScalar.C:38
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:424
surfaceScalarField & phi
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:428
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.
U
Definition: pEqn.H:83
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
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:230
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > linearInterpolate(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: linear.H:107
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
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:497
const Internal::FieldType & primitiveField() const
Return a const-reference to the internal field.
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
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:1011
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
dynamicFvMesh & mesh
A class for handling words, derived from string.
Definition: word.H:59
Calulate the matrix for the first temporal derivative.
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:91
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1036
const dimensionSet dimPressure
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:218
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:57
const dimensionedVector & g
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
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.
const dimensionSet dimDensity
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
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
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:876
const word cloudName(propsDict.lookup("cloudName"))
virtual vector velocityCorrection(typename CloudType::parcelType &p, const scalar deltaT) const
Calculate the velocity correction.
Definition: Implicit.C:332
A class for managing temporary objects.
Definition: PtrList.H:53
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
const dimensionSet dimVelocity