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-2015 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"
33 
34 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
35 
36 template<class CloudType>
38 (
39  const dictionary& dict,
40  CloudType& owner
41 )
42 :
43  PackingModel<CloudType>(dict, owner, typeName),
44  alpha_
45  (
46  this->owner().name() + ":alpha",
47  this->owner().theta()
48  ),
49  phiCorrect_(NULL),
50  uCorrect_(NULL),
51  applyGravity_(this->coeffDict().lookup("applyGravity")),
52  alphaMin_(readScalar(this->coeffDict().lookup("alphaMin"))),
53  rhoMin_(readScalar(this->coeffDict().lookup("rhoMin")))
54 {
55  alpha_.oldTime();
56 }
57 
58 
59 template<class CloudType>
61 (
62  const Implicit<CloudType>& cm
63 )
64 :
66  alpha_(cm.alpha_),
67  phiCorrect_(cm.phiCorrect_()),
68  uCorrect_(cm.uCorrect_()),
69  applyGravity_(cm.applyGravity_),
70  alphaMin_(cm.alphaMin_),
71  rhoMin_(cm.rhoMin_)
72 {
73  alpha_.oldTime();
74 }
75 
76 
77 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
78 
79 template<class CloudType>
81 {}
82 
83 
84 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
85 
86 template<class CloudType>
88 {
90 
91  if (store)
92  {
93  const fvMesh& mesh = this->owner().mesh();
94  const dimensionedScalar deltaT = this->owner().db().time().deltaT();
95  const word& cloudName = this->owner().name();
96 
97  const dimensionedVector& g = this->owner().g();
98  const volScalarField& rhoc = this->owner().rho();
99 
100  const AveragingMethod<scalar>& rhoAverage =
102  (
103  cloudName + ":rhoAverage"
104  );
105  const AveragingMethod<scalar>& uSqrAverage =
107  (
108  cloudName + ":uSqrAverage"
109  );
110 
111  mesh.setFluxRequired(alpha_.name());
112 
113  // Property fields
114  // ~~~~~~~~~~~~~~~
115 
116  // volume fraction field
117  alpha_ = max(this->owner().theta(), alphaMin_);
118  alpha_.correctBoundaryConditions();
119 
120  // average density
122  (
123  IOobject
124  (
125  cloudName + ":rho",
126  this->owner().db().time().timeName(),
127  mesh,
128  IOobject::NO_READ,
129  IOobject::NO_WRITE
130  ),
131  mesh,
132  dimensionedScalar("zero", dimDensity, 0),
134  );
135  rho.internalField() = max(rhoAverage.internalField(), rhoMin_);
137 
138  //Info << " x: " << mesh.C().internalField().component(2) << endl;
139  //Info << " alpha: " << alpha_.internalField() << endl;
140  //Info << "alphaOld: " << alpha_.oldTime().internalField() << endl;
141  //Info << " rho: " << rho.internalField() << endl;
142  //Info << endl;
143 
144 
145  // Stress field
146  // ~~~~~~~~~~~~
147 
148  // stress derivative wrt volume fraction
149  volScalarField tauPrime
150  (
151  IOobject
152  (
153  cloudName + ":tauPrime",
154  this->owner().db().time().timeName(),
155  mesh,
156  IOobject::NO_READ,
157  IOobject::NO_WRITE
158  ),
159  mesh,
160  dimensionedScalar("zero", dimPressure, 0),
162  );
163 
164  tauPrime.internalField() =
165  this->particleStressModel_->dTaudTheta
166  (
167  alpha_.internalField(),
168  rho.internalField(),
169  uSqrAverage.internalField()
170  )();
171 
172  tauPrime.correctBoundaryConditions();
173 
174 
175  // Implicit solution for the volume fraction
176  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
177 
179  tauPrimeByRhoAf
180  (
181  "tauPrimeByRhoAf",
182  fvc::interpolate(deltaT*tauPrime/rho)
183  );
184 
185  fvScalarMatrix alphaEqn
186  (
187  fvm::ddt(alpha_)
188  - fvc::ddt(alpha_)
189  - fvm::laplacian(tauPrimeByRhoAf, alpha_)
190  );
191 
192  if (applyGravity_)
193  {
195  phiGByA
196  (
197  "phiGByA",
198  deltaT*(g & mesh.Sf())*fvc::interpolate(1.0 - rhoc/rho)
199  );
200 
201  alphaEqn += fvm::div(phiGByA, alpha_);
202  }
203 
204  alphaEqn.solve();
205 
206 
207  // Generate correction fields
208  // ~~~~~~~~~~~~~~~~~
209 
210  // correction volumetric flux
211  phiCorrect_ = tmp<surfaceScalarField>
212  (
214  (
215  cloudName + ":phiCorrect",
216  alphaEqn.flux()/fvc::interpolate(alpha_)
217  )
218  );
219 
220  // correction velocity
221  uCorrect_ = tmp<volVectorField>
222  (
223  new volVectorField
224  (
225  cloudName + ":uCorrect",
226  fvc::reconstruct(phiCorrect_())
227  )
228 
229  );
230  uCorrect_->correctBoundaryConditions();
231 
232  //Info << endl;
233  //Info << " alpha: " << alpha_.internalField() << endl;
234  //Info << "phiCorrect: " << phiCorrect_->internalField() << endl;
235  //Info << " uCorrect: " << uCorrect_->internalField() << endl;
236  //Info << endl;
237  }
238  else
239  {
240  alpha_.oldTime();
241  phiCorrect_.clear();
242  uCorrect_.clear();
243  }
244 }
245 
246 
247 template<class CloudType>
249 (
250  typename CloudType::parcelType& p,
251  const scalar deltaT
252 ) const
253 {
254  const fvMesh& mesh = this->owner().mesh();
255 
256  // containing tetrahedron and parcel coordinates within
257  const label cellI = p.cell();
258  const label faceI = p.tetFace();
259  const tetIndices tetIs(cellI, faceI, p.tetPt(), mesh);
260  List<scalar> tetCoordinates(4);
261  tetIs.tet(mesh).barycentric(p.position(), tetCoordinates);
262 
263  // cell velocity
264  const vector U = uCorrect_()[cellI];
265 
266  // face geometry
267  vector nHat = mesh.faces()[faceI].normal(mesh.points());
268  const scalar nMag = mag(nHat);
269  nHat /= nMag;
270 
271  // get face flux
272  scalar phi;
273  const label patchI = mesh.boundaryMesh().whichPatch(faceI);
274  if (patchI == -1)
275  {
276  phi = phiCorrect_()[faceI];
277  }
278  else
279  {
280  phi =
281  phiCorrect_().boundaryField()[patchI]
282  [
283  mesh.boundaryMesh()[patchI].whichFace(faceI)
284  ];
285  }
286 
287  // interpolant equal to 1 at the cell centre and 0 at the face
288  const scalar t = tetCoordinates[0];
289 
290  // the normal component of the velocity correction is interpolated linearly
291  // the tangential component is equal to that at the cell centre
292  return U + (1.0 - t)*nHat*(phi/nMag - (U & nHat));
293 }
294 
295 
296 // ************************************************************************* //
const dimensionSet dimPressure
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
const word cloudName(propsDict.lookup("cloudName"))
#define readScalar
Definition: doubleScalar.C:38
dimensioned< scalar > mag(const dimensioned< Type > &)
const surfaceVectorField & Sf() const
Return cell face area vectors.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > flux() const
Return the face-flux field from the matrix.
Definition: fvMatrix.C:883
void clear() const
If object pointer points to valid object:
Definition: tmpI.H:174
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
void setFluxRequired(const word &name) const
Definition: fvSchemes.C:497
A class for handling words, derived from string.
Definition: word.H:59
virtual vector velocityCorrection(typename CloudType::parcelType &p, const scalar deltaT) const
Calculate the velocity correction.
Definition: Implicit.C:249
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
Calculate the matrix for the divergence of the given field and flux.
Calulate the matrix for the first temporal derivative.
InternalField & internalField()
Return internal field.
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
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:68
tmp< surfaceScalarField > interpolate(const RhoType &rho)
Base class for packing models.
Definition: MPPICCloud.H:53
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
dynamicFvMesh & mesh
const dimensionSet dimDensity
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
ParcelType parcelType
Type of parcel the cloud was instantiated for.
Definition: DSMCCloud.H:217
This boundary condition applies a zero-gradient condition from the patch internal field onto the patc...
dictionary dict
stressControl lookup("compactNormalStress") >> compactNormalStress
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Reconstruct volField from a face flux field.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
const Type & lookupObject(const word &name) const
Lookup and return the object of the given Type.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
solverPerformance solve(const dictionary &)
Solve segregated or coupled returning the solution statistics.
Definition: fvMatrixSolve.C:57
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
surfaceScalarField & phi
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
Calculate the matrix for the laplacian of the field.
virtual ~Implicit()
Destructor.
Definition: Implicit.C:80
const dimensionedVector & g
virtual tmp< Field< Type > > internalField() const =0
Return an internal field of the average.
void correctBoundaryConditions()
Correct boundary field.
virtual void cacheFields(const bool store)
Calculate the inter particles stresses.
Definition: Implicit.C:87
Implicit(const dictionary &dict, CloudType &owner)
Construct from components.
Definition: Implicit.C:38
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1060
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:73
Implicit model for applying an inter-particle stress to the particles.
Definition: Implicit.H:56
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
U
Definition: pEqn.H:82
word timeName
Definition: getTimeIndex.H:3