ThermalPhaseChangePhaseSystem.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) 2015-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 
28 #include "fvcVolumeIntegrate.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 template<class BasePhaseSystem>
35 (
36  const fvMesh& mesh
37 )
38 :
39  HeatAndMassTransferPhaseSystem<BasePhaseSystem>(mesh),
40  volatile_(this->lookup("volatile")),
41  saturationModel_(saturationModel::New(this->subDict("saturationModel"))),
42  massTransfer_(this->lookup("massTransfer"))
43 {
44 
46  (
48  this->phasePairs_,
49  phasePairIter
50  )
51  {
52  const phasePair& pair(phasePairIter());
53 
54  if (pair.ordered())
55  {
56  continue;
57  }
58 
59  // Initially assume no mass transfer
60  iDmdt_.insert
61  (
62  pair,
63  new volScalarField
64  (
65  IOobject
66  (
67  IOobject::groupName("iDmdt", pair.name()),
68  this->mesh().time().timeName(),
69  this->mesh(),
72  ),
73  this->mesh(),
75  )
76  );
77  }
78 }
79 
80 
81 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
82 
83 template<class BasePhaseSystem>
86 {}
87 
88 
89 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
90 
91 template<class BasePhaseSystem>
94 {
95  return saturationModel_();
96 }
97 
98 
99 template<class BasePhaseSystem>
102 {
103  typedef compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
104  alphatPhaseChangeWallFunction;
105 
106  autoPtr<phaseSystem::heatTransferTable> eqnsPtr =
108 
109  phaseSystem::heatTransferTable& eqns = eqnsPtr();
110 
111  // Accumulate mDotL contributions from boundaries
113  (
115  this->phasePairs_,
116  phasePairIter
117  )
118  {
119  const phasePair& pair(phasePairIter());
120 
121  if (pair.ordered())
122  {
123  continue;
124  }
125 
126  const phaseModel& phase = pair.phase1();
127  const phaseModel& otherPhase = pair.phase2();
128 
129  volScalarField mDotL
130  (
131  IOobject
132  (
133  "mDotL",
134  phase.mesh().time().timeName(),
135  phase.mesh(),
138  false
139  ),
140  phase.mesh(),
141  dimensionedScalar("",dimensionSet(1,-1,-3,0,0),0.0)
142  );
143 
144  if
145  (
146  otherPhase.mesh().foundObject<volScalarField>
147  (
148  "alphat." + otherPhase.name()
149  )
150  )
151  {
152  const volScalarField& alphat =
153  otherPhase.mesh().lookupObject<volScalarField>
154  (
155  "alphat." + otherPhase.name()
156  );
157 
158  const fvPatchList& patches = this->mesh().boundary();
159  forAll(patches, patchi)
160  {
161  const fvPatch& currPatch = patches[patchi];
162 
163  if
164  (
165  isA<alphatPhaseChangeWallFunction>
166  (
167  alphat.boundaryField()[patchi]
168  )
169  )
170  {
171  const scalarField& patchMDotL =
172  refCast<const alphatPhaseChangeWallFunction>
173  (
174  alphat.boundaryField()[patchi]
175  ).mDotL();
176 
177  forAll(patchMDotL,facei)
178  {
179  label faceCelli = currPatch.faceCells()[facei];
180  mDotL[faceCelli] = patchMDotL[facei];
181  }
182  }
183  }
184  }
185 
186  *eqns[otherPhase.name()] -= mDotL;
187 
188  }
189 
190  return eqnsPtr;
191 }
192 
193 
194 template<class BasePhaseSystem>
197 {
198  // Create a mass transfer matrix for each species of each phase
199  autoPtr<phaseSystem::massTransferTable> eqnsPtr
200  (
202  );
203 
204  phaseSystem::massTransferTable& eqns = eqnsPtr();
205 
206  forAll(this->phaseModels_, phasei)
207  {
208  const phaseModel& phase = this->phaseModels_[phasei];
209 
210  const PtrList<volScalarField>& Yi = phase.Y();
211 
212  forAll(Yi, i)
213  {
214  eqns.insert
215  (
216  Yi[i].name(),
217  new fvScalarMatrix(Yi[i], dimMass/dimTime)
218  );
219  }
220  }
221 
223  (
225  this->phasePairs_,
226  phasePairIter
227  )
228  {
229  const phasePair& pair(phasePairIter());
230 
231  if (pair.ordered())
232  {
233  continue;
234  }
235  const phaseModel& phase = pair.phase1();
236  const phaseModel& otherPhase = pair.phase2();
237 
238  const word name
239  (
240  IOobject::groupName(volatile_, phase.name())
241  );
242 
243  const word otherName
244  (
245  IOobject::groupName(volatile_, otherPhase.name())
246  );
247 
248  const volScalarField dmdt(this->dmdt(pair));
249  const volScalarField dmdt12(posPart(dmdt));
250  const volScalarField dmdt21(negPart(dmdt));
251 
252  *eqns[name] += fvm::Sp(dmdt21, eqns[name]->psi()) - dmdt21;
253  *eqns[otherName] += dmdt12 - fvm::Sp(dmdt12, eqns[otherName]->psi());
254  }
255 
256  return eqnsPtr;
257 }
258 
259 
260 template<class BasePhaseSystem>
263 (
264  const phasePairKey& key
265 ) const
266 {
267  const scalar dmdtSign(Pair<word>::compare(iDmdt_.find(key).key(), key));
268 
269  return dmdtSign**iDmdt_[key];
270 }
271 
272 
273 template<class BasePhaseSystem>
276 (
277  const Foam::phaseModel& phase
278 ) const
279 {
280  tmp<volScalarField> tiDmdt
281  (
282  new volScalarField
283  (
284  IOobject
285  (
286  IOobject::groupName("iDmdt", phase.name()),
287  this->mesh_.time().timeName(),
288  this->mesh_
289  ),
290  this->mesh_,
292  )
293  );
294 
296  (
298  this->phasePairs_,
299  phasePairIter
300  )
301  {
302  const phasePair& pair(phasePairIter());
303 
304  if (pair.ordered())
305  {
306  continue;
307  }
308 
309  const phaseModel* phase1 = &pair.phase1();
310  const phaseModel* phase2 = &pair.phase2();
311 
312  forAllConstIter(phasePair, pair, iter)
313  {
314  if (phase1 == &phase)
315  {
316  tiDmdt.ref() += this->iDmdt(pair);
317  }
318 
319  Swap(phase1, phase2);
320  }
321  }
322 
323  return tiDmdt;
324 }
325 
326 
327 template<class BasePhaseSystem>
329 {
330  typedef compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
331  alphatPhaseChangeWallFunction;
332 
334 
336  (
338  this->phasePairs_,
339  phasePairIter
340  )
341  {
342  const phasePair& pair(phasePairIter());
343 
344  if (pair.ordered())
345  {
346  continue;
347  }
348 
349  const phaseModel& phase1 = pair.phase1();
350  const phaseModel& phase2 = pair.phase2();
351 
352  Info<< phase1.name() << " min/max T "
353  << min(phase1.thermo().T()).value()
354  << " - "
355  << max(phase1.thermo().T()).value()
356  << endl;
357 
358  Info<< phase2.name() << " min/max T "
359  << min(phase2.thermo().T()).value()
360  << " - "
361  << max(phase2.thermo().T()).value()
362  << endl;
363 
364  const volScalarField& T1(phase1.thermo().T());
365  const volScalarField& T2(phase2.thermo().T());
366 
367  const volScalarField& he1(phase1.thermo().he());
368  const volScalarField& he2(phase2.thermo().he());
369 
370  volScalarField& dmdt(*this->dmdt_[pair]);
371  volScalarField& iDmdt(*this->iDmdt_[pair]);
372 
373  volScalarField& Tf = *this->Tf_[pair];
374 
375  volScalarField hef1(phase1.thermo().he(phase1.thermo().p(), Tf));
376  volScalarField hef2(phase2.thermo().he(phase2.thermo().p(), Tf));
377 
379  (
380  min
381  (
382  (pos0(iDmdt)*he2 + neg(iDmdt)*hef2)
383  - (neg(iDmdt)*he1 + pos0(iDmdt)*hef1),
384  0.3*mag(hef2 - hef1)
385  )
386  );
387 
388  volScalarField iDmdtNew(iDmdt);
389 
390  if (massTransfer_ )
391  {
392  volScalarField H1
393  (
394  this->heatTransferModels_[pair][pair.first()]->K(0)
395  );
396 
397  volScalarField H2
398  (
399  this->heatTransferModels_[pair][pair.second()]->K(0)
400  );
401 
402  Tf = saturationModel_->Tsat(phase1.thermo().p());
403 
404  iDmdtNew =
405  (H1*(Tf - T1) + H2*(Tf - T2))/L;
406 
407  }
408  else
409  {
410  iDmdtNew == dimensionedScalar("0",dmdt.dimensions(), 0);
411  }
412 
413  volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
414  volScalarField H2(this->heatTransferModels_[pair][pair.second()]->K());
415 
416  // Limit the H[12] boundary field to avoid /0
417  const scalar HLimit = 1e-4;
418  H1.boundaryFieldRef() =
419  max(H1.boundaryField(), phase1.boundaryField()*HLimit);
420  H2.boundaryFieldRef() =
421  max(H2.boundaryField(), phase2.boundaryField()*HLimit);
422 
423  Tf = (H1*T1 + H2*T2 + iDmdtNew*L)/(H1 + H2);
424 
425  Info<< "Tf." << pair.name()
426  << ": min = " << min(Tf.primitiveField())
427  << ", mean = " << average(Tf.primitiveField())
428  << ", max = " << max(Tf.primitiveField())
429  << endl;
430 
431  scalar iDmdtRelax(this->mesh().fieldRelaxationFactor("iDmdt"));
432  iDmdt = (1 - iDmdtRelax)*iDmdt + iDmdtRelax*iDmdtNew;
433 
434  if (massTransfer_ )
435  {
436  Info<< "iDmdt." << pair.name()
437  << ": min = " << min(iDmdt.primitiveField())
438  << ", mean = " << average(iDmdt.primitiveField())
439  << ", max = " << max(iDmdt.primitiveField())
440  << ", integral = " << fvc::domainIntegrate(iDmdt).value()
441  << endl;
442  }
443 
444  // Accumulate dmdt contributions from boundaries
445  volScalarField wDmdt
446  (
447  IOobject
448  (
449  IOobject::groupName("wDmdt", pair.name()),
450  this->mesh().time().timeName(),
451  this->mesh(),
454  false
455  ),
456  this->mesh(),
458  );
459 
460  if
461  (
462  phase2.mesh().foundObject<volScalarField>
463  (
464  "alphat." + phase2.name()
465  )
466  )
467  {
468  const volScalarField& alphat =
469  phase2.mesh().lookupObject<volScalarField>
470  (
471  "alphat." + phase2.name()
472  );
473 
474  const fvPatchList& patches = this->mesh().boundary();
475  forAll(patches, patchi)
476  {
477  const fvPatch& currPatch = patches[patchi];
478 
479  if
480  (
481  isA<alphatPhaseChangeWallFunction>
482  (
483  alphat.boundaryField()[patchi]
484  )
485  )
486  {
487  const scalarField& patchDmdt =
488  refCast<const alphatPhaseChangeWallFunction>
489  (
490  alphat.boundaryField()[patchi]
491  ).dmdt();
492 
493  forAll(patchDmdt,facei)
494  {
495  label faceCelli = currPatch.faceCells()[facei];
496  wDmdt[faceCelli] += patchDmdt[facei];
497  }
498  }
499  }
500 
501  Info<< "wDmdt." << pair.name()
502  << ": min = " << min(wDmdt.primitiveField())
503  << ", mean = " << average(wDmdt.primitiveField())
504  << ", max = " << max(wDmdt.primitiveField())
505  << ", integral = " << fvc::domainIntegrate(wDmdt).value()
506  << endl;
507  }
508 
509  dmdt = wDmdt + iDmdt;
510 
511  Info<< "dmdt." << pair.name()
512  << ": min = " << min(dmdt.primitiveField())
513  << ", mean = " << average(dmdt.primitiveField())
514  << ", max = " << max(dmdt.primitiveField())
515  << ", integral = " << fvc::domainIntegrate(dmdt).value()
516  << endl;
517  }
518 }
519 
520 
521 template<class BasePhaseSystem>
523 {
524  if (BasePhaseSystem::read())
525  {
526  bool readOK = true;
527 
528  // Models ...
529 
530  return readOK;
531  }
532  else
533  {
534  return false;
535  }
536 }
537 
538 
539 // ************************************************************************* //
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
label phasei
Definition: pEqn.H:27
virtual autoPtr< phaseSystem::heatTransferTable > heatTransfer() const
Return the heat transfer matrices.
const double e
Elementary charge.
Definition: doubleFloat.H:78
HashTable< autoPtr< phasePair >, phasePairKey, phasePairKey::hash > phasePairTable
Definition: phaseSystem.H:133
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
ThermalPhaseChangePhaseSystem(const fvMesh &)
Construct from fvMesh.
static int compare(const Pair< word > &a, const Pair< word > &b)
Compare Pairs.
Definition: Pair.H:141
virtual void correctThermo()
Correct the thermodynamics.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
virtual autoPtr< phaseSystem::massTransferTable > massTransfer() const
Return the mass transfer matrices.
const word & name() const
Definition: phaseModel.H:109
dimensioned< Type > domainIntegrate(const GeometricField< Type, fvPatchField, volMesh > &vf)
void read(Istream &, label &, const dictionary &)
In-place read with dictionary lookup.
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
dimensionedScalar posPart(const dimensionedScalar &ds)
dimensionedScalar neg(const dimensionedScalar &ds)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
HashPtrTable< fvScalarMatrix, word, string::hash > heatTransferTable
Definition: phaseSystem.H:109
stressControl lookup("compactNormalStress") >> compactNormalStress
dynamicFvMesh & mesh
void Swap(T &a, T &b)
Definition: Swap.H:43
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static word groupName(Name name, const word &group)
virtual autoPtr< phaseSystem::heatTransferTable > heatTransfer() const
Return the heat transfer matrices.
Volume integrate volField creating a volField.
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
dimensionedScalar pos0(const dimensionedScalar &ds)
volScalarField & he2
Definition: EEqns.H:3
const Mesh & mesh() const
Return mesh.
virtual bool read()
Read base phaseProperties dictionary.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
const dimensionSet dimDensity
const saturationModel & saturation() const
Return the saturationModel.
PtrList< fvPatch > fvPatchList
container classes for fvPatch
Definition: fvPatchList.H:45
label patchi
virtual ~ThermalPhaseChangePhaseSystem()
Destructor.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
const volScalarField & psi
HashPtrTable< fvScalarMatrix, word, string::hash > massTransferTable
Definition: phaseSystem.H:118
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
virtual tmp< volScalarField > iDmdt(const phasePairKey &key) const
Return the interfacial mass flow rate.
A class for managing temporary objects.
Definition: PtrList.H:53
Single incompressible phase derived from the phase-fraction. Used as part of the multiPhaseMixture fo...
Definition: phaseModel.H:53
mixture correctThermo()
dimensionedScalar negPart(const dimensionedScalar &ds)