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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2013-2017 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "Implicit.H"
30 #include "fvmDdt.H"
31 #include "fvmDiv.H"
32 #include "fvmLaplacian.H"
33 #include "fvcReconstruct.H"
34 #include "volPointInterpolation.H"
36 
37 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
38 
39 template<class CloudType>
41 (
42  const dictionary& dict,
43  CloudType& owner
44 )
45 :
46  PackingModel<CloudType>(dict, owner, typeName),
47  alpha_
48  (
49  IOobject
50  (
51  this->owner().name() + ":alpha",
52  this->owner().db().time().timeName(),
53  this->owner().mesh(),
54  IOobject::NO_READ,
55  IOobject::NO_WRITE
56  ),
57  this->owner().mesh(),
59  zeroGradientFvPatchScalarField::typeName
60  ),
61  phiCorrect_(nullptr),
62  uCorrect_(nullptr),
63  applyLimiting_(this->coeffDict().lookup("applyLimiting")),
64  applyGravity_(this->coeffDict().lookup("applyGravity")),
65  alphaMin_(this->coeffDict().getScalar("alphaMin")),
66  rhoMin_(this->coeffDict().getScalar("rhoMin"))
67 {
68  alpha_ = this->owner().theta();
69  alpha_.oldTime();
70 }
71 
72 
73 template<class CloudType>
75 (
76  const Implicit<CloudType>& cm
77 )
78 :
80  alpha_(cm.alpha_),
81  phiCorrect_(cm.phiCorrect_()),
82  uCorrect_(cm.uCorrect_()),
83  applyLimiting_(cm.applyLimiting_),
84  applyGravity_(cm.applyGravity_),
85  alphaMin_(cm.alphaMin_),
86  rhoMin_(cm.rhoMin_)
87 {
88  alpha_.oldTime();
89 }
90 
91 
92 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
93 
94 template<class CloudType>
96 {}
97 
98 
99 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
100 
101 template<class CloudType>
103 {
105 
106  if (store)
107  {
108  const fvMesh& mesh = this->owner().mesh();
109  const dimensionedScalar deltaT = this->owner().db().time().deltaT();
110  const word& cloudName = this->owner().name();
111 
112  const dimensionedVector& g = this->owner().g();
113  const volScalarField& rhoc = this->owner().rho();
114 
115  const AveragingMethod<scalar>& rhoAverage =
116  mesh.lookupObject<AveragingMethod<scalar>>
117  (
118  cloudName + ":rhoAverage"
119  );
120  const AveragingMethod<vector>& uAverage =
121  mesh.lookupObject<AveragingMethod<vector>>
122  (
123  cloudName + ":uAverage"
124  );
125  const AveragingMethod<scalar>& uSqrAverage =
126  mesh.lookupObject<AveragingMethod<scalar>>
127  (
128  cloudName + ":uSqrAverage"
129  );
130 
131  mesh.setFluxRequired(alpha_.name());
132 
133  // Property fields
134  // ~~~~~~~~~~~~~~~
135 
136  // volume fraction field
137  alpha_ = max(this->owner().theta(), alphaMin_);
138  alpha_.correctBoundaryConditions();
139 
140  // average density
142  (
143  IOobject
144  (
145  cloudName + ":rho",
146  this->owner().db().time().timeName(),
147  mesh,
148  IOobject::NO_READ,
149  IOobject::NO_WRITE
150  ),
151  mesh,
154  );
155  rho.primitiveFieldRef() = max(rhoAverage.primitiveField(), rhoMin_);
156  rho.correctBoundaryConditions();
157 
158  // Stress field
159  // ~~~~~~~~~~~~
160 
161  // stress derivative wrt volume fraction
162  volScalarField tauPrime
163  (
164  IOobject
165  (
166  cloudName + ":tauPrime",
167  this->owner().db().time().timeName(),
168  mesh,
169  IOobject::NO_READ,
170  IOobject::NO_WRITE
171  ),
172  mesh,
175  );
176 
177  tauPrime.primitiveFieldRef() =
178  this->particleStressModel_->dTaudTheta
179  (
180  alpha_.primitiveField(),
181  rho.primitiveField(),
182  uSqrAverage.primitiveField()
183  )();
184 
185  tauPrime.correctBoundaryConditions();
186 
187 
188  // Gravity flux
189  // ~~~~~~~~~~~~
190 
191  tmp<surfaceScalarField> phiGByA;
192 
193  if (applyGravity_)
194  (
195  phiGByA = tmp<surfaceScalarField>
196  (
198  (
199  "phiGByA",
200  deltaT*(g & mesh.Sf())*fvc::interpolate(1.0 - rhoc/rho)
201  )
202  )
203  );
204 
205 
206  // Implicit solution for the volume fraction
207  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
208 
210  tauPrimeByRhoAf
211  (
212  "tauPrimeByRhoAf",
213  fvc::interpolate(deltaT*tauPrime/rho)
214  );
215 
216  fvScalarMatrix alphaEqn
217  (
218  fvm::ddt(alpha_)
219  - fvc::ddt(alpha_)
220  - fvm::laplacian(tauPrimeByRhoAf, alpha_)
221  );
222 
223  if (applyGravity_)
224  {
225  alphaEqn += fvm::div(phiGByA(), alpha_);
226  }
227 
228  alphaEqn.solve();
229 
230 
231  // Generate correction fields
232  // ~~~~~~~~~~~~~~~~~
233 
234  // correction volumetric flux
235  phiCorrect_ = tmp<surfaceScalarField>
236  (
238  (
239  cloudName + ":phiCorrect",
240  alphaEqn.flux()/fvc::interpolate(alpha_)
241  )
242  );
243 
244  // limit the correction flux
245  if (applyLimiting_)
246  {
248  (
249  IOobject
250  (
251  cloudName + ":U",
252  this->owner().db().time().timeName(),
253  mesh,
254  IOobject::NO_READ,
255  IOobject::NO_WRITE
256  ),
257  mesh,
260  );
261  U.primitiveFieldRef() = uAverage.primitiveField();
262  U.correctBoundaryConditions();
263 
265  (
266  cloudName + ":phi",
267  linearInterpolate(U) & mesh.Sf()
268  );
269 
270  if (applyGravity_)
271  {
272  phiCorrect_.ref() -= phiGByA();
273  }
274 
275  forAll(phiCorrect_(), facei)
276  {
277  // Current and correction fluxes
278  const scalar phiCurr = phi[facei];
279  scalar& phiCorr = phiCorrect_.ref()[facei];
280 
281  // Don't limit if the correction is in the opposite direction to
282  // the flux. We need all the help we can get in this state.
283  if (phiCurr*phiCorr < 0)
284  {}
285 
286  // If the correction and the flux are in the same direction then
287  // don't apply any more correction than is already present in
288  // the flux.
289  else if (phiCorr > 0)
290  {
291  phiCorr = max(phiCorr - phiCurr, 0);
292  }
293  else
294  {
295  phiCorr = min(phiCorr - phiCurr, 0);
296  }
297  }
298 
299  if (applyGravity_)
300  {
301  phiCorrect_.ref() += phiGByA();
302  }
303  }
304 
305  // correction velocity
306  uCorrect_ = tmp<volVectorField>
307  (
308  new volVectorField
309  (
310  cloudName + ":uCorrect",
311  fvc::reconstruct(phiCorrect_())
312  )
313 
314  );
315  uCorrect_->correctBoundaryConditions();
316 
317  //Info << endl;
318  //Info << " alpha: " << alpha_.primitiveField() << endl;
319  //Info << "phiCorrect: " << phiCorrect_->primitiveField() << endl;
320  //Info << " uCorrect: " << uCorrect_->primitiveField() << endl;
321  //Info << endl;
322  }
323  else
324  {
325  alpha_.oldTime();
326  phiCorrect_.clear();
327  uCorrect_.clear();
328  }
329 }
330 
331 
332 template<class CloudType>
334 (
335  typename CloudType::parcelType& p,
336  const scalar deltaT
337 ) const
338 {
339  const fvMesh& mesh = this->owner().mesh();
340 
341  // containing tetrahedron and parcel coordinates within
342  const label celli = p.cell();
343  const label facei = p.tetFace();
344 
345  // cell velocity
346  const vector U = uCorrect_()[celli];
347 
348  // face geometry
349  vector nHat = mesh.faces()[facei].areaNormal(mesh.points());
350  const scalar nMag = mag(nHat);
351  nHat /= nMag;
352 
353  // get face flux
354  scalar phi;
355  const label patchi = mesh.boundaryMesh().whichPatch(facei);
356  if (patchi == -1)
357  {
358  phi = phiCorrect_()[facei];
359  }
360  else
361  {
362  phi =
363  phiCorrect_().boundaryField()[patchi]
364  [
365  mesh.boundaryMesh()[patchi].whichFace(facei)
366  ];
367  }
368 
369  // interpolant equal to 1 at the cell centre and 0 at the face
370  const scalar t = p.coordinates()[0];
371 
372  // the normal component of the velocity correction is interpolated linearly
373  // the tangential component is equal to that at the cell centre
374  return U + (1.0 - t)*nHat*(phi/nMag - (U & nHat));
375 }
376 
377 
378 // ************************************************************************* //
Foam::AveragingMethod< scalar >
Foam::fvc::reconstruct
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcReconstruct.C:56
Foam::dimPressure
const dimensionSet dimPressure
Foam::IOobject
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:169
p
volScalarField & p
Definition: createFieldRefs.H:8
cloudName
const word cloudName(propsDict.get< word >("cloud"))
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::dimVelocity
const dimensionSet dimVelocity
Foam::dimDensity
const dimensionSet dimDensity
Foam::AveragingMethod::primitiveField
virtual tmp< Field< Type > > primitiveField() const =0
Return an internal field of the average.
Foam::PackingModels::Implicit::~Implicit
virtual ~Implicit()
Destructor.
Definition: Implicit.C:95
Foam::PackingModels::Implicit::cacheFields
virtual void cacheFields(const bool store)
Calculate the inter particles stresses.
Definition: Implicit.C:102
Foam::fac::div
tmp< GeometricField< Type, faPatchField, areaMesh > > div(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Definition: facDiv.C:50
rho
rho
Definition: readInitialConditions.H:88
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::zeroGradientFvPatchField
This boundary condition applies a zero-gradient condition from the patch internal field onto the patc...
Definition: zeroGradientFvPatchField.H:64
Foam::fixedValueFvPatchField
This boundary condition supplies a fixed value constraint, and is the base class for a number of othe...
Definition: fixedValueFvPatchField.H:80
fixedValueFvsPatchField.H
fvmDiv.H
Calculate the matrix for the divergence of the given field and flux.
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::dimensionedVector
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Definition: dimensionedVector.H:50
Foam::PackingModels::Implicit::velocityCorrection
virtual vector velocityCorrection(typename CloudType::parcelType &p, const scalar deltaT) const
Calculate the velocity correction.
Definition: Implicit.C:334
Foam::DSMCCloud::mesh
const fvMesh & mesh() const
Return reference to the mesh.
Definition: DSMCCloudI.H:44
fvmLaplacian.H
Calculate the matrix for the laplacian of the field.
Foam::dimensionedScalar
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Definition: dimensionedScalarFwd.H:42
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::radiation::lookup
Lookup type of boundary radiation properties.
Definition: lookup.H:63
Foam::DSMCCloud
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:71
timeName
word timeName
Definition: getTimeIndex.H:3
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::PackingModel
Base class for packing models.
Definition: KinematicCloud.H:98
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::dimensioned< scalar >
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:85
g
const uniformDimensionedVectorField & g
Definition: createFluidFields.H:26
Foam::PackingModels::Implicit::Implicit
Implicit(const dictionary &dict, CloudType &owner)
Construct from components.
Definition: Implicit.C:41
Foam::GeometricField::primitiveFieldRef
Internal::FieldType & primitiveFieldRef(const bool updateAccessTime=true)
Return a reference to the internal field.
Definition: GeometricField.C:766
Foam::GeometricField::correctBoundaryConditions
void correctBoundaryConditions()
Correct boundary field.
Definition: GeometricField.C:940
volPointInterpolation.H
U
U
Definition: pEqn.H:72
Foam::linearInterpolate
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > linearInterpolate(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: linear.H:112
Foam::Vector< scalar >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::fac::ddt
tmp< GeometricField< Type, faPatchField, areaMesh > > ddt(const dimensioned< Type > dt, const faMesh &mesh)
Definition: facDdt.C:47
Foam::PackingModels::Implicit
Implicit model for applying an inter-particle stress to the particles.
Definition: Implicit.H:61
Foam::fvMatrix
A special matrix type and solver, designed for finite volume solutions of scalar equations....
Definition: fvPatchField.H:68
Foam::fvMatrix::solve
SolverPerformance< Type > solve(const dictionary &)
Solve returning the solution statistics.
Definition: fvMatrixSolve.C:319
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
fvcReconstruct.H
Reconstruct volField from a face flux field.
Foam::DSMCCloud::parcelType
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:220
Foam::fac::laplacian
tmp< GeometricField< Type, faPatchField, areaMesh > > laplacian(const GeometricField< Type, faPatchField, areaMesh > &vf, const word &name)
Definition: facLaplacian.C:47
Foam::GeometricField< scalar, fvPatchField, volMesh >
zeroGradientFvPatchFields.H
Foam::fvMatrix::flux
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > flux() const
Return the face-flux field from the matrix.
Definition: fvMatrix.C:1534
Foam::fac::interpolate
static tmp< GeometricField< Type, faePatchField, edgeMesh > > interpolate(const GeometricField< Type, faPatchField, areaMesh > &tvf, const edgeScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
Implicit.H
fvmDdt.H
Calculate the matrix for the first temporal derivative.
Foam::dimless
const dimensionSet dimless
Dimensionless.
Definition: dimensionSets.C:189