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 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 
27 #include "fvCFD.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  // Initially assume no mass transfer
79  wDmdt_.insert
80  (
81  pair,
82  new volScalarField
83  (
84  IOobject
85  (
86  IOobject::groupName("wDmdt", pair.name()),
87  this->mesh().time().timeName(),
88  this->mesh(),
91  ),
92  this->mesh(),
94  )
95  );
96  }
97 }
98 
99 
100 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
101 
102 template<class BasePhaseSystem>
105 {}
106 
107 
108 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
109 
110 template<class BasePhaseSystem>
113 {
114  return saturationModel_();
115 }
116 
117 
118 template<class BasePhaseSystem>
121 {
122  // Create a mass transfer matrix for each species of each phase
123  autoPtr<phaseSystem::massTransferTable> eqnsPtr
124  (
126  );
127 
128  phaseSystem::massTransferTable& eqns = eqnsPtr();
129 
130  forAll(this->phaseModels_, phasei)
131  {
132  const phaseModel& phase = this->phaseModels_[phasei];
133 
134  const PtrList<volScalarField>& Yi = phase.Y();
135 
136  forAll(Yi, i)
137  {
138  eqns.insert
139  (
140  Yi[i].name(),
141  new fvScalarMatrix(Yi[i], dimMass/dimTime)
142  );
143  }
144  }
145 
147  (
149  this->phasePairs_,
150  phasePairIter
151  )
152  {
153  const phasePair& pair(phasePairIter());
154 
155  if (pair.ordered())
156  {
157  continue;
158  }
159  const phaseModel& phase = pair.phase1();
160  const phaseModel& otherPhase = pair.phase2();
161 
162  const word name
163  (
164  IOobject::groupName(volatile_, phase.name())
165  );
166 
167  const word otherName
168  (
169  IOobject::groupName(volatile_, otherPhase.name())
170  );
171 
172  const volScalarField dmdt(this->dmdt(pair));
173  const volScalarField dmdt12(posPart(dmdt));
174  const volScalarField dmdt21(negPart(dmdt));
175 
176  *eqns[name] += fvm::Sp(dmdt21, eqns[name]->psi()) - dmdt21;
177  *eqns[otherName] += dmdt12 - fvm::Sp(dmdt12, eqns[otherName]->psi());
178  }
179 
180  return eqnsPtr;
181 }
182 
183 
184 template<class BasePhaseSystem>
186 {
187  typedef compressible::alphatPhaseChangeWallFunctionFvPatchScalarField
188  alphatPhaseChangeWallFunction;
189 
191 
193  (
195  this->phasePairs_,
196  phasePairIter
197  )
198  {
199  const phasePair& pair(phasePairIter());
200 
201  if (pair.ordered())
202  {
203  continue;
204  }
205 
206  const phaseModel& phase1 = pair.phase1();
207  const phaseModel& phase2 = pair.phase2();
208 
209  const volScalarField& T1(phase1.thermo().T());
210  const volScalarField& T2(phase2.thermo().T());
211 
212  const volScalarField& he1(phase1.thermo().he());
213  const volScalarField& he2(phase2.thermo().he());
214 
215  volScalarField& dmdt(*this->dmdt_[pair]);
216  volScalarField& iDmdt(*this->iDmdt_[pair]);
217  volScalarField& wDmdt(*this->wDmdt_[pair]);
218 
219  volScalarField& Tf = *this->Tf_[pair];
220 
221  volScalarField hef1(phase1.thermo().he(phase1.thermo().p(), Tf));
222  volScalarField hef2(phase2.thermo().he(phase2.thermo().p(), Tf));
223 
224  if (massTransfer_ )
225  {
226  volScalarField H1
227  (
228  this->heatTransferModels_[pair][pair.first()]->K(0)
229  );
230 
231  volScalarField H2
232  (
233  this->heatTransferModels_[pair][pair.second()]->K(0)
234  );
235 
236  Tf = saturationModel_->Tsat(phase1.thermo().p());
237 
238  scalar iDmdtRelax(this->mesh().fieldRelaxationFactor("iDmdt"));
239 
240  iDmdt =
241  (1 - iDmdtRelax)*iDmdt
242  + iDmdtRelax*(H1*(Tf - T1) + H2*(Tf - T2))
243  /min
244  (
245  (pos(iDmdt)*he2 + neg(iDmdt)*hef2)
246  - (neg(iDmdt)*he1 + pos(iDmdt)*hef1),
247  0.3*mag(hef2 - hef1)
248  );
249 
250  Info<< "iDmdt." << pair.name()
251  << ": min = " << min(iDmdt.internalField())
252  << ", mean = " << average(iDmdt.internalField())
253  << ", max = " << max(iDmdt.internalField())
254  << ", integral = " << fvc::domainIntegrate(iDmdt).value()
255  << endl;
256  }
257  else
258  {
259  iDmdt == dimensionedScalar("0", dmdt.dimensions(), 0);
260  }
261 
262  volScalarField H1(this->heatTransferModels_[pair][pair.first()]->K());
263  volScalarField H2(this->heatTransferModels_[pair][pair.second()]->K());
264 
265  // Limit the H[12] boundary field to avoid /0
266  const scalar HLimit = 1e-4;
267  H1.boundaryField() =
268  max(H1.boundaryField(), phase1.boundaryField()*HLimit);
269  H2.boundaryField() =
270  max(H2.boundaryField(), phase2.boundaryField()*HLimit);
271 
272  volScalarField mDotL
273  (
274  iDmdt*
275  (
276  (pos(iDmdt)*he2 + neg(iDmdt)*hef2)
277  - (neg(iDmdt)*he1 + pos(iDmdt)*hef1)
278  )
279  );
280 
281  Tf = (H1*T1 + H2*T2 + mDotL)/(H1 + H2);
282 
283  Info<< "Tf." << pair.name()
284  << ": min = " << min(Tf.internalField())
285  << ", mean = " << average(Tf.internalField())
286  << ", max = " << max(Tf.internalField())
287  << endl;
288 
289  // Accumulate dmdt contributions from boundaries
290  if
291  (
292  phase2.mesh().foundObject<volScalarField>
293  (
294  "alphat." + phase2.name()
295  )
296  )
297  {
298  scalar wDmdtRelax(this->mesh().fieldRelaxationFactor("wDmdt"));
299  wDmdt *= (1 - wDmdtRelax);
300 
301  const volScalarField& alphat =
302  phase2.mesh().lookupObject<volScalarField>
303  (
304  "alphat." + phase2.name()
305  );
306 
307  const fvPatchList& patches = this->mesh().boundary();
308  forAll(patches, patchi)
309  {
310  const fvPatch& currPatch = patches[patchi];
311 
312  if
313  (
314  isA<alphatPhaseChangeWallFunction>
315  (
316  alphat.boundaryField()[patchi]
317  )
318  )
319  {
320  const scalarField& patchDmdt =
321  refCast<const alphatPhaseChangeWallFunction>
322  (
323  alphat.boundaryField()[patchi]
324  ).dmdt();
325 
326  forAll(patchDmdt,facei)
327  {
328  label faceCelli = currPatch.faceCells()[facei];
329  wDmdt[faceCelli] += wDmdtRelax*patchDmdt[facei];
330  }
331  }
332  }
333 
334  Info<< "wDmdt." << pair.name()
335  << ": min = " << min(wDmdt.internalField())
336  << ", mean = " << average(wDmdt.internalField())
337  << ", max = " << max(wDmdt.internalField())
338  << ", integral = " << fvc::domainIntegrate(wDmdt).value()
339  << endl;
340  }
341 
342  dmdt = wDmdt + iDmdt;
343 
344  Info<< "dmdt." << pair.name()
345  << ": min = " << min(dmdt.internalField())
346  << ", mean = " << average(dmdt.internalField())
347  << ", max = " << max(dmdt.internalField())
348  << ", integral = " << fvc::domainIntegrate(dmdt).value()
349  << endl;
350  }
351 }
352 
353 
354 template<class BasePhaseSystem>
356 {
357  if (BasePhaseSystem::read())
358  {
359  bool readOK = true;
360 
361  // Models ...
362 
363  return readOK;
364  }
365  else
366  {
367  return false;
368  }
369 }
370 
371 
372 // ************************************************************************* //
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
virtual autoPtr< phaseSystem::massTransferTable > massTransfer() const
Return the mass transfer matrices.
const saturationModel & saturation() const
Return the saturationModel.
dimensionedScalar neg(const dimensionedScalar &ds)
HashTable< autoPtr< phasePair >, phasePairKey, phasePairKey::hash > phasePairTable
Definition: phaseSystem.H:133
dimensioned< scalar > mag(const dimensioned< Type > &)
label phasei
Definition: pEqn.H:37
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
dimensionedScalar posPart(const dimensionedScalar &ds)
virtual ~ThermalPhaseChangePhaseSystem()
Destructor.
fluid correctThermo()
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
messageStream Info
const Mesh & mesh() const
Return mesh.
dynamicFvMesh & mesh
const dimensionSet dimDensity
tmp< fvMatrix< Type > > Sp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
autoPtr< BasicCompressibleTurbulenceModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleTurbulenceModel::transportModel &transport, const word &propertiesName)
virtual void correctThermo()
Correct the thermodynamics.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static word groupName(Name name, const word &group)
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
const double e
Elementary charge.
Definition: doubleFloat.H:78
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
stressControl lookup("compactNormalStress") >> compactNormalStress
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
dimensionedScalar negPart(const dimensionedScalar &ds)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
HashPtrTable< fvScalarMatrix, word, string::hash > massTransferTable
Definition: phaseSystem.H:118
#define forAll(list, i)
Definition: UList.H:421
label patchi
ThermalPhaseChangePhaseSystem(const fvMesh &)
Construct from fvMesh.
const volScalarField & psi
Definition: createFields.H:24
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
virtual bool read()
Read base phaseProperties dictionary.
dimensioned< Type > domainIntegrate(const GeometricField< Type, fvPatchField, volMesh > &vf)
bool read(const char *, int32_t &)
Definition: int32IO.C:87
dimensionedScalar pos(const dimensionedScalar &ds)
volScalarField & he2
Definition: EEqns.H:3
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
PtrList< fvPatch > fvPatchList
container classes for fvPatch
Definition: fvPatchList.H:45
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:117