twoPhaseSystem.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-2018 OpenFOAM Foundation
9 Copyright (C) 2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "twoPhaseSystem.H"
31#include "BlendedInterfacialModel.H"
32#include "virtualMassModel.H"
33#include "heatTransferModel.H"
34#include "liftModel.H"
35#include "wallLubricationModel.H"
36#include "turbulentDispersionModel.H"
37#include "fvMatrix.H"
38#include "surfaceInterpolate.H"
39#include "MULES.H"
40#include "subCycle.H"
41#include "fvcDdt.H"
42#include "fvcDiv.H"
43#include "fvcSnGrad.H"
44#include "fvcFlux.H"
45#include "fvcCurl.H"
46#include "fvmDdt.H"
47#include "fvmLaplacian.H"
49#include "blendingMethod.H"
50#include "HashPtrTable.H"
51#include "UniformField.H"
52
53// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
54
56(
57 const fvMesh& mesh,
58 const dimensionedVector& g
59)
60:
62 (
64 (
65 "phaseProperties",
66 mesh.time().constant(),
67 mesh,
68 IOobject::MUST_READ_IF_MODIFIED,
69 IOobject::NO_WRITE
70 )
71 ),
72
73 mesh_(mesh),
74
75 phase1_
76 (
77 *this,
78 *this,
79 this->get<wordList>("phases")[0]
80 ),
81
82 phase2_
83 (
84 *this,
85 *this,
86 this->get<wordList>("phases")[1]
87 ),
88
89 phi_
90 (
92 (
93 "phi",
94 mesh.time().timeName(),
95 mesh,
96 IOobject::NO_READ,
97 IOobject::AUTO_WRITE
98 ),
99 this->calcPhi()
100 ),
101
102 dgdt_
103 (
105 (
106 "dgdt",
107 mesh.time().timeName(),
108 mesh,
109 IOobject::READ_IF_PRESENT,
110 IOobject::AUTO_WRITE
111 ),
112 mesh,
113 dimensionedScalar("dgdt", dimless/dimTime, Zero)
114 )
115{
116 phase2_.volScalarField::operator=(scalar(1) - phase1_);
117
118
119 // Blending
120 forAllConstIter(dictionary, subDict("blending"), iter)
121 {
122 blendingMethods_.insert
123 (
124 iter().dict().dictName(),
126 (
127 iter().dict(),
128 this->get<wordList>("phases")
129 )
130 );
131 }
132
133
134 // Pairs
135
136 phasePair::scalarTable sigmaTable(lookup("sigma"));
137 phasePair::dictTable aspectRatioTable(lookup("aspectRatio"));
138
139 pair_.reset
140 (
141 new phasePair
142 (
143 phase1_,
144 phase2_,
145 g,
146 sigmaTable
147 )
148 );
149
150 pair1In2_.reset
151 (
153 (
154 phase1_,
155 phase2_,
156 g,
157 sigmaTable,
158 aspectRatioTable
159 )
160 );
161
162 pair2In1_.reset
163 (
165 (
166 phase2_,
167 phase1_,
168 g,
169 sigmaTable,
170 aspectRatioTable
171 )
172 );
173
174
175 // Models
176
177 drag_.reset
178 (
180 (
181 lookup("drag"),
182 (
183 blendingMethods_.found("drag")
184 ? blendingMethods_["drag"]
185 : blendingMethods_["default"]
186 ),
187 pair_,
188 pair1In2_,
189 pair2In1_,
190 false // Do not zero drag coefficient at fixed-flux BCs
191 )
192 );
193
194 virtualMass_.reset
195 (
197 (
198 lookup("virtualMass"),
199 (
200 blendingMethods_.found("virtualMass")
201 ? blendingMethods_["virtualMass"]
202 : blendingMethods_["default"]
203 ),
204 pair_,
205 pair1In2_,
206 pair2In1_
207 )
208 );
209
210 heatTransfer_.reset
211 (
213 (
214 lookup("heatTransfer"),
215 (
216 blendingMethods_.found("heatTransfer")
217 ? blendingMethods_["heatTransfer"]
218 : blendingMethods_["default"]
219 ),
220 pair_,
221 pair1In2_,
222 pair2In1_
223 )
224 );
225
226 lift_.reset
227 (
229 (
230 lookup("lift"),
231 (
232 blendingMethods_.found("lift")
233 ? blendingMethods_["lift"]
234 : blendingMethods_["default"]
235 ),
236 pair_,
237 pair1In2_,
238 pair2In1_
239 )
240 );
241
242 wallLubrication_.reset
243 (
245 (
246 lookup("wallLubrication"),
247 (
248 blendingMethods_.found("wallLubrication")
249 ? blendingMethods_["wallLubrication"]
250 : blendingMethods_["default"]
251 ),
252 pair_,
253 pair1In2_,
254 pair2In1_
255 )
256 );
257
258 turbulentDispersion_.reset
259 (
261 (
262 lookup("turbulentDispersion"),
263 (
264 blendingMethods_.found("turbulentDispersion")
265 ? blendingMethods_["turbulentDispersion"]
266 : blendingMethods_["default"]
267 ),
268 pair_,
269 pair1In2_,
270 pair2In1_
271 )
272 );
273}
274
275
276// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
277
279{}
280
281
282// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
283
285{
286 return phase1_*phase1_.thermo().rho() + phase2_*phase2_.thermo().rho();
287}
288
289
291{
292 return phase1_*phase1_.U() + phase2_*phase2_.U();
293}
294
295
296Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::calcPhi() const
297{
298 return
299 fvc::interpolate(phase1_)*phase1_.phi()
300 + fvc::interpolate(phase2_)*phase2_.phi();
301}
302
303
305{
306 return drag_->K();
307}
308
309
311{
312 return drag_->Kf();
313}
314
315
317{
318 return virtualMass_->K();
319}
320
321
323{
324 return virtualMass_->Kf();
325}
326
327
329{
330 return heatTransfer_->K();
331}
332
333
335{
336 return lift_->F<vector>() + wallLubrication_->F<vector>();
337}
338
339
341{
342 return lift_->Ff() + wallLubrication_->Ff();
343}
344
345
347{
348 return turbulentDispersion_->D();
349}
350
351
353{
354 const Time& runTime = mesh_.time();
355
356 volScalarField& alpha1 = phase1_;
357 volScalarField& alpha2 = phase2_;
358
359 const surfaceScalarField& phi1 = phase1_.phi();
360 const surfaceScalarField& phi2 = phase2_.phi();
361
362 const dictionary& alphaControls = mesh_.solverDict
363 (
364 alpha1.name()
365 );
366
367 label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
368 label nAlphaCorr(alphaControls.get<label>("nAlphaCorr"));
369
370 word alphaScheme("div(phi," + alpha1.name() + ')');
371 word alpharScheme("div(phir," + alpha1.name() + ')');
372
374
375
376 surfaceScalarField phic("phic", phi_);
377 surfaceScalarField phir("phir", phi1 - phi2);
378
379 tmp<surfaceScalarField> alpha1alpha2f;
380
381 if (pPrimeByA_.valid())
382 {
383 alpha1alpha2f =
384 fvc::interpolate(max(alpha1, scalar(0)))
385 *fvc::interpolate(max(alpha2, scalar(0)));
386
388 (
389 pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
390 );
391
392 phir += phiP;
393 }
394
395 for (int acorr=0; acorr<nAlphaCorr; acorr++)
396 {
398 (
399 IOobject
400 (
401 "Sp",
403 mesh_
404 ),
405 mesh_,
406 dimensionedScalar("Sp", dgdt_.dimensions(), 0.0)
407 );
408
410 (
411 IOobject
412 (
413 "Su",
415 mesh_
416 ),
417 // Divergence term is handled explicitly to be
418 // consistent with the explicit transport solution
419 fvc::div(phi_)*min(alpha1, scalar(1))
420 );
421
422 forAll(dgdt_, celli)
423 {
424 if (dgdt_[celli] > 0.0)
425 {
426 Sp[celli] -= dgdt_[celli]/max(1.0 - alpha1[celli], 1e-4);
427 Su[celli] += dgdt_[celli]/max(1.0 - alpha1[celli], 1e-4);
428 }
429 else if (dgdt_[celli] < 0.0)
430 {
431 Sp[celli] += dgdt_[celli]/max(alpha1[celli], 1e-4);
432 }
433 }
434
435 surfaceScalarField alphaPhic1
436 (
438 (
439 phic,
440 alpha1,
441 alphaScheme
442 )
443 + fvc::flux
444 (
445 -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
446 alpha1,
448 )
449 );
450
451 phase1_.correctInflowOutflow(alphaPhic1);
452
453 if (nAlphaSubCycles > 1)
454 {
455 for
456 (
457 subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
458 !(++alphaSubCycle).end();
459 )
460 {
461 surfaceScalarField alphaPhic10(alphaPhic1);
462
464 (
465 geometricOneField(),
466 alpha1,
467 phi_,
468 alphaPhic10,
469 (alphaSubCycle.index()*Sp)(),
470 (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
471 UniformField<scalar>(phase1_.alphaMax()),
472 zeroField()
473 );
474
475 if (alphaSubCycle.index() == 1)
476 {
477 phase1_.alphaPhi() = alphaPhic10;
478 }
479 else
480 {
481 phase1_.alphaPhi() += alphaPhic10;
482 }
483 }
484
485 phase1_.alphaPhi() /= nAlphaSubCycles;
486 }
487 else
488 {
490 (
491 geometricOneField(),
492 alpha1,
493 phi_,
494 alphaPhic1,
495 Sp,
496 Su,
497 UniformField<scalar>(phase1_.alphaMax()),
498 zeroField()
499 );
500
501 phase1_.alphaPhi() = alphaPhic1;
502 }
503
504 if (pPrimeByA_.valid())
505 {
506 fvScalarMatrix alpha1Eqn
507 (
509 - fvm::laplacian(alpha1alpha2f()*pPrimeByA_(), alpha1, "bounded")
510 );
511
512 alpha1Eqn.relax();
513 alpha1Eqn.solve();
514
515 phase1_.alphaPhi() += alpha1Eqn.flux();
516 }
517
518 phase1_.alphaRhoPhi() =
519 fvc::interpolate(phase1_.rho())*phase1_.alphaPhi();
520
521 phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
522 phase2_.correctInflowOutflow(phase2_.alphaPhi());
523 phase2_.alphaRhoPhi() =
524 fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();
525
526 Info<< alpha1.name() << " volume fraction = "
527 << alpha1.weightedAverage(mesh_.V()).value()
528 << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
529 << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
530 << endl;
531
532 // Ensure the phase-fractions are bounded
533 alpha1.max(0);
534 alpha1.min(1);
535
536 alpha2 = scalar(1) - alpha1;
537 }
538}
539
540
542{
543 phase1_.correct();
544 phase2_.correct();
545}
546
547
549{
550 phase1_.turbulence().correct();
551 phase2_.turbulence().correct();
552}
553
554
556{
557 if (regIOobject::read())
558 {
559 bool readOK = true;
560
561 readOK &= phase1_.read(*this);
562 readOK &= phase2_.read(*this);
563
564 // models ...
565
566 return readOK;
567 }
568
569 return false;
570}
571
572
574{
575 return pair_->sigma();
576}
577
578
579// ************************************************************************* //
surfaceScalarField phic(mixture.cAlpha() *mag(alphaPhic/mesh.magSf()))
MULES: Multidimensional universal limiter for explicit solution.
const uniformDimensionedVectorField & g
const volScalarField & alpha1
surfaceScalarField & phi2
const volScalarField & alpha2
surfaceScalarField & phi1
dimensioned< Type > weightedAverage(const DimensionedField< scalar, GeoMesh > &weightField) const
Calculate and return weighted average.
DimensionedField< scalar, volMesh > Internal
The internal field type from which this GeometricField is derived.
void max(const dimensioned< Type > &dt)
Use the maximum of the field and specified value.
void min(const dimensioned< Type > &dt)
Use the minimum of the field and specified value.
void correctBoundaryConditions()
Correct boundary field.
bool found(const Key &key) const
Return true if hashed entry is found in table.
Definition: HashTableI.H:100
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:180
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:57
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
static word timeName(const scalar t, const int precision=precision_)
Definition: Time.C:780
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Definition: dictionary.C:460
word dictName() const
The local dictionary name (final part of scoped name)
Definition: dictionaryI.H:60
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
const Time & time() const noexcept
Return time registry.
constant condensation/saturation model.
Description for mass transfer between a pair of phases. The direction of the mass transfer is from th...
Definition: phasePair.H:56
Lookup type of boundary radiation properties.
Definition: lookup.H:66
virtual bool read()
Read object.
A class for managing temporary objects.
Definition: tmp.H:65
Class which solves the volume fraction equations for two phases.
tmp< surfaceScalarField > Ff() const
Return the combined face-force (lift + wall-lubrication)
tmp< volVectorField > F() const
Return the combined force (lift + wall-lubrication)
tmp< volScalarField > D() const
Return the turbulent diffusivity.
tmp< volScalarField > Kd() const
Return the drag coefficient.
void correct()
Correct two-phase properties other than turbulence.
virtual ~twoPhaseSystem()
Destructor.
phaseModel & phase1_
Phase model 1.
tmp< surfaceScalarField > Kdf() const
Return the face drag coefficient.
tmp< surfaceScalarField > Vmf() const
Return the face virtual mass coefficient.
tmp< volVectorField > U() const
Return the mixture velocity.
phaseModel & phase2_
Phase model 2.
void correctTurbulence()
Correct two-phase turbulence.
tmp< volScalarField > sigma() const
Return the surface tension coefficient.
tmp< volScalarField > rho() const
Return the mixture density.
tmp< volScalarField > Kh() const
Return the heat transfer coefficient.
virtual void solve()
Solve for the phase fractions.
bool read()
Read base phaseProperties dictionary.
tmp< volScalarField > Vm() const
Return the virtual mass coefficient.
A class for handling words, derived from Foam::string.
Definition: word.H:68
dynamicFvMesh & mesh
engineTime & runTime
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
word alpharScheme("div(phirb,alpha)")
Calculate the curl of the given volField by constructing the Hodge-dual of the symmetric part of the ...
Calculate the first temporal derivative.
Calculate the divergence of the given field.
Calculate the face-flux of the given field.
Calculate the snGrad of the given volField.
Calculate the matrix for the first temporal derivative.
Calculate the matrix for the laplacian of the field.
zeroField Su
Definition: alphaSuSp.H:1
zeroField Sp
Definition: alphaSuSp.H:2
word timeName
Definition: getTimeIndex.H:3
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
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.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:47
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:49
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:47
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:34
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:48
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:48
vtkSmartPointer< vtkFloatArray > zeroField(const word &name, const label size)
Create named field initialized to zero.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
const dimensionSet dimless
Dimensionless.
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:45
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
messageStream Info
Information stream (stdout output on master, null elsewhere)
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:158
dictionary dict
volScalarField & e
Definition: createFields.H:11
const dictionary & alphaControls
Definition: alphaControls.H:1
label nAlphaSubCycles(alphaControls.get< label >("nAlphaSubCycles"))
label nAlphaCorr(alphaControls.get< label >("nAlphaCorr"))
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object.
Definition: stdFoam.H:381