MassTransferPhaseSystem.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) 2017-2022 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
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
29#include "HashPtrTable.H"
30#include "fvcDiv.H"
31#include "fvmSup.H"
32#include "fvMatrix.H"
33#include "volFields.H"
35
36using namespace Foam::multiphaseInter;
37
38// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
39
40template<class BasePhaseSystem>
42(
43 const fvMesh& mesh
44)
45:
46 BasePhaseSystem(mesh)
47{
48 this->generatePairsAndSubModels("massTransferModel", massTransferModels_);
49
51 {
52 if (!dmdt_.found(iterModel()->pair()))
53 {
55 (
56 iterModel()->pair(),
58 (
60 (
61 IOobject::groupName("dmdt",iterModel()->pair().name()),
62 this->mesh().time().timeName(),
63 this->mesh(),
66 ),
67 this->mesh(),
69 )
70 );
71 }
72 }
73}
74
75
76// * * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * //
77
78template<class BasePhaseSystem>
81(
82 const volScalarField& dmdtNetki,
83 const phasePairKey& keyik,
84 const phasePairKey& keyki,
85 const volScalarField& T
86) const
87{
89 (
91 (
92 "tL",
93 this->mesh().time().timeName(),
94 this->mesh(),
97 ),
98 this->mesh(),
100 );
101 auto& L = tL.ref();
102
103 if (massTransferModels_.found(keyik))
104 {
105 const autoPtr<interfaceCompositionModel>& interfacePtr =
106 massTransferModels_[keyik];
107
108 word speciesName = interfacePtr->transferSpecie();
109
110 const word species(speciesName.substr(0, speciesName.find('.')));
111
112 L -= neg(dmdtNetki)*interfacePtr->L(species, T);
113 }
114
115 if (massTransferModels_.found(keyki))
116 {
117 const autoPtr<interfaceCompositionModel>& interfacePtr =
118 massTransferModels_[keyki];
119
120 word speciesName = interfacePtr->transferSpecie();
121
122 const word species(speciesName.substr(0, speciesName.find('.')));
123
124 L -= pos(dmdtNetki)*interfacePtr->L(species, T);
125 }
126
127 return tL;
128}
129
130
131// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
132
133template<class BasePhaseSystem>
136(
137 const phasePairKey& key
138) const
139{
140 auto tdmdt = tmp<volScalarField>::New
141 (
143 (
144 "dmdt",
145 this->mesh().time().timeName(),
146 this->mesh()
147 ),
148 this->mesh(),
150 );
151 auto& dmdt = tdmdt.ref();
152
153 if (dmdt_.found(key))
154 {
155 dmdt = *dmdt_[key];
156 }
157
158 return tdmdt;
159}
160
161
162template<class BasePhaseSystem>
165(
166 const volScalarField& T
167)
168{
170 auto& eqn = teqn.ref();
171
172 forAllConstIters(this->phaseModels_, iteri)
173 {
174 const phaseModel& phasei = iteri()();
175
176 auto iterk = iteri;
177
178 for (++iterk; iterk != this->phaseModels_.end(); ++iterk)
179 {
180 if (iteri()().name() != iterk()().name())
181 {
182 const phaseModel& phasek = iterk()();
183
184 // Phase i to phase k
185 const phasePairKey keyik(phasei.name(), phasek.name(), true);
186
187 // Phase k to phase i
188 const phasePairKey keyki(phasek.name(), phasei.name(), true);
189
190 // Net mass transfer from k to i phase
191 auto tdmdtNetki = tmp<volScalarField>::New
192 (
194 (
195 "tdmdtYki",
196 this->mesh().time().timeName(),
197 this->mesh()
198 ),
199 this->mesh(),
201 );
202 auto& dmdtNetki = tdmdtNetki.ref();
203
204 auto tSp = tmp<volScalarField>::New
205 (
207 (
208 "Sp",
209 this->mesh().time().timeName(),
210 this->mesh()
211 ),
212 this->mesh(),
214 );
215 auto& Sp = tSp.ref();
216
217 auto tSu = tmp<volScalarField>::New
218 (
220 (
221 "Su",
222 this->mesh().time().timeName(),
223 this->mesh()
224 ),
225 this->mesh(),
227 );
228 auto& Su = tSu.ref();
229
230
231 if (massTransferModels_.found(keyik))
232 {
234 massTransferModels_[keyik];
235
236 dmdtNetki -= *dmdt_[keyik];
237
239 interfacePtr->KSp(interfaceCompositionModel::T, T);
240
241 if (KSp.valid())
242 {
243 Sp += KSp.ref();
244 }
245
247 interfacePtr->KSu(interfaceCompositionModel::T, T);
248
249 if (KSu.valid())
250 {
251 Su += KSu.ref();
252 }
253
254 // If linearization is not provided used full explicit
255 if (!KSp.valid() && !KSu.valid())
256 {
257 Su += *dmdt_[keyik];
258 }
259 }
260
261 // Looking for mass transfer in the other direction (k to i)
262 if (massTransferModels_.found(keyki))
263 {
265 massTransferModels_[keyki];
266
267 dmdtNetki += *dmdt_[keyki];
268
269
271 interfacePtr->KSp(interfaceCompositionModel::T, T);
272
273 if (KSp.valid())
274 {
275 Sp += KSp.ref();
276 }
277
279 interfacePtr->KSu(interfaceCompositionModel::T, T);
280
281 if (KSu.valid())
282 {
283 Su += KSu.ref();
284 }
285
286 // If linearization is not provided used full explicit
287 if (!KSp.valid() && !KSu.valid())
288 {
289 Su += *dmdt_[keyki];
290 }
291 }
292
293 tmp<volScalarField> L = calculateL(dmdtNetki, keyik, keyki, T);
294
295 //eqn -= dmdtNetki*L;
296 eqn -= fvm::Sp(Sp*L.ref(), T) + Su*L.ref();
297 }
298 }
299 }
300 return teqn;
301}
302
303
304template<class BasePhaseSystem>
307(
308 const volScalarField& p
309)
310{
312 auto& eqn = teqn.ref();
313
314 auto tSp = tmp<volScalarField>::New
315 (
317 (
318 "Sp",
319 this->mesh().time().timeName(),
320 this->mesh()
321 ),
322 this->mesh(),
324 );
325 auto& Sp = tSp.ref();
326
327 auto tSu = tmp<volScalarField>::New
328 (
330 (
331 "Su",
332 this->mesh().time().timeName(),
333 this->mesh()
334 ),
335 this->mesh(),
337 );
338 auto& Su = tSu.ref();
339
340 forAllConstIters(this->totalPhasePairs(), iter)
341 {
342 const phasePair& pair = iter()();
343
344 const phaseModel& phase1 = pair.phase1();
345 const phaseModel& phase2 = pair.phase2();
346
347 const phasePairKey key12
348 (
349 phase1.name(),
350 phase2.name(),
351 true
352 );
353
354 if (massTransferModels_.found(key12))
355 {
357 massTransferModels_[key12];
358
360 interfacePtr->KSp(interfaceCompositionModel::P, p);
361
362 if (KSp.valid())
363 {
364 Sp -=
365 KSp.ref()
366 *(
367 - this->coeffs(phase1.name())
368 + this->coeffs(phase2.name())
369 );
370 }
371
373 interfacePtr->KSu(interfaceCompositionModel::P, p);
374
375 if (KSu.valid())
376 {
377 Su -=
378 KSu.ref()
379 *(
380 - this->coeffs(phase1.name())
381 + this->coeffs(phase2.name())
382 );
383 }
384
385 // If linearization is not provided used full explicit
386 if (!KSp.valid() && !KSu.valid())
387 {
388 Su -=
389 *dmdt_[key12]
390 *(
391 - this->coeffs(phase1.name())
392 + this->coeffs(phase2.name())
393 );
394 }
395 }
396
397 const phasePairKey key21
398 (
399 phase2.name(),
400 phase1.name(),
401 true
402 );
403
404 if (massTransferModels_.found(key21))
405 {
407 massTransferModels_[key21];
408
410 interfacePtr->KSp(interfaceCompositionModel::P, p);
411
412 if (KSp.valid())
413 {
414 Sp +=
415 KSp.ref()
416 *(
417 - this->coeffs(phase1.name())
418 + this->coeffs(phase2.name())
419 );
420 }
421
423 interfacePtr->KSu(interfaceCompositionModel::P, p);
424
425 if (KSu.valid())
426 {
427 Su +=
428 KSu.ref()
429 *(
430 - this->coeffs(phase1.name())
431 + this->coeffs(phase2.name())
432 );
433 }
434
435 // If linearization is not provided used full explicit
436 if (!KSp.valid() && !KSu.valid())
437 {
438 Su +=
439 *dmdt_[key21]
440 *(
441 - this->coeffs(phase1.name())
442 + this->coeffs(phase2.name())
443 );
444 }
445 }
446
447 }
448
449 eqn += fvm::Sp(Sp, p) + Su;
450 return teqn;
451}
452
453
454template<class BasePhaseSystem>
456(
457 const volScalarField& T
458)
459{
460 forAllConstIters(this->phaseModels_, iteri)
461 {
462 const phaseModel& phasei = iteri()();
463
464 auto iterk = iteri;
465
466 for (++iterk; iterk != this->phaseModels_.end(); ++iterk)
467 {
468 if (iteri()().name() != iterk()().name())
469 {
470 const phaseModel& phasek = iterk()();
471
472 // Phase i to phase k
473 const phasePairKey keyik(phasei.name(), phasek.name(), true);
474
475 // Phase k to phase i
476 const phasePairKey keyki(phasek.name(), phasei.name(), true);
477
478 if (massTransferModels_.found(keyik))
479 {
481 massTransferModels_[keyik];
482
483 tmp<volScalarField> Kexp = interfacePtr->Kexp(T);
484
485 *dmdt_[keyik] = Kexp.ref();
486
487 }
488
489 if (massTransferModels_.found(keyki))
490 {
492 massTransferModels_[keyki];
493
494 // Explicit temperature mass transfer rate
495 const tmp<volScalarField> Kexp = interfacePtr->Kexp(T);
496
497 *dmdt_[keyki] = Kexp.ref();
498 }
499 }
500 }
501 }
502}
503
504
505template<class BasePhaseSystem>
507(
508 SuSpTable& Su,
510)
511{
512 // This term adds/subtracts alpha*div(U) as a source term
513 // for alpha, substituting div(U) = mDot(1/rho1 - 1/rho2)
514 bool includeDivU(true);
515
516 forAllConstIters(this->totalPhasePairs(), iter)
517 {
518 const phasePair& pair = iter()();
519
520 const phaseModel& phase1 = pair.phase1();
521 const phaseModel& phase2 = pair.phase2();
522
523 const volScalarField& alpha1 = pair.phase1();
524 const volScalarField& alpha2 = pair.phase2();
525
526 tmp<volScalarField> tCoeffs1 = this->coeffs(phase1.name());
527 const volScalarField& coeffs1 = tCoeffs1();
528
529 tmp<volScalarField> tCoeffs2 = this->coeffs(phase2.name());
530 const volScalarField& coeffs2 = tCoeffs2();
531
532 // Phase 1 to phase 2
533 const phasePairKey key12
534 (
535 phase1.name(),
536 phase2.name(),
537 true
538 );
539
540 tmp<volScalarField> tdmdt12(this->dmdt(key12));
541 volScalarField& dmdt12 = tdmdt12.ref();
542
543 if (massTransferModels_.found(key12))
544 {
546 massTransferModels_[key12];
547
549 interfacePtr->KSu(interfaceCompositionModel::alpha, phase1);
550
551 if (KSu.valid())
552 {
553 dmdt12 = KSu.ref();
554 }
555
556 includeDivU = interfacePtr->includeDivU();
557 }
558
559
560 // Phase 2 to phase 1
561 const phasePairKey key21
562 (
563 phase2.name(),
564 phase1.name(),
565 true
566 );
567
568 tmp<volScalarField> tdmdt21(this->dmdt(key21));
569 volScalarField& dmdt21 = tdmdt21.ref();
570
571 if (massTransferModels_.found(key21))
572 {
574 massTransferModels_[key21];
575
577 interfacePtr->KSu(interfaceCompositionModel::alpha, phase2);
578
579 if (KSu.valid())
580 {
581 dmdt21 = KSu.ref();
582 }
583
584 includeDivU = interfacePtr->includeDivU();
585 }
586
587 volScalarField::Internal& SpPhase1 = Sp[phase1.name()];
588
589 volScalarField::Internal& SuPhase1 = Su[phase1.name()];
590
591 volScalarField::Internal& SpPhase2 = Sp[phase2.name()];
592
593 volScalarField::Internal& SuPhase2 = Su[phase2.name()];
594
595 const volScalarField dmdtNet(dmdt21 - dmdt12);
596
597 const volScalarField coeffs12(coeffs1 - coeffs2);
598
599 const surfaceScalarField& phi = this->phi();
600
601 if (includeDivU)
602 {
603 SuPhase1 +=
604 fvc::div(phi)*min(max(alpha1, scalar(0)), scalar(1));
605
606 SuPhase2 +=
607 fvc::div(phi)*min(max(alpha2, scalar(0)), scalar(1));
608 }
609
610 // NOTE: dmdtNet is distributed in terms =
611 // Source for phase 1 =
612 // dmdtNet/rho1
613 // - alpha1*dmdtNet(1/rho1 - 1/rho2)
614
615 forAll(dmdtNet, celli)
616 {
617 scalar dmdt21 = dmdtNet[celli];
618 scalar coeffs12Cell = coeffs12[celli];
619
620 scalar alpha1Limited = max(min(alpha1[celli], 1.0), 0.0);
621
622 // exp.
623 SuPhase1[celli] += coeffs1[celli]*dmdt21;
624
625 if (includeDivU)
626 {
627 if (dmdt21 > 0)
628 {
629 if (coeffs12Cell > 0)
630 {
631 // imp
632 SpPhase1[celli] -= dmdt21*coeffs12Cell;
633 }
634 else if (coeffs12Cell < 0)
635 {
636 // exp
637 SuPhase1[celli] -=
638 dmdt21*coeffs12Cell*alpha1Limited;
639 }
640 }
641 else if (dmdt21 < 0)
642 {
643 if (coeffs12Cell > 0)
644 {
645 // exp
646 SuPhase1[celli] -=
647 dmdt21*coeffs12Cell*alpha1Limited;
648 }
649 else if (coeffs12Cell < 0)
650 {
651 // imp
652 SpPhase1[celli] -= dmdt21*coeffs12Cell;
653 }
654 }
655 }
656 }
657
658 forAll(dmdtNet, celli)
659 {
660 scalar dmdt12 = -dmdtNet[celli];
661 scalar coeffs21Cell = -coeffs12[celli];
662
663 scalar alpha2Limited = max(min(alpha2[celli], 1.0), 0.0);
664
665 // exp
666 SuPhase2[celli] += coeffs2[celli]*dmdt12;
667
668 if (includeDivU)
669 {
670 if (dmdt12 > 0)
671 {
672 if (coeffs21Cell > 0)
673 {
674 // imp
675 SpPhase2[celli] -= dmdt12*coeffs21Cell;
676 }
677 else if (coeffs21Cell < 0)
678 {
679 // exp
680 SuPhase2[celli] -=
681 dmdt12*coeffs21Cell*alpha2Limited;
682 }
683 }
684 else if (dmdt12 < 0)
685 {
686 if (coeffs21Cell > 0)
687 {
688 // exp
689 SuPhase2[celli] -=
690 coeffs21Cell*dmdt12*alpha2Limited;
691 }
692 else if (coeffs21Cell < 0)
693 {
694 // imp
695 SpPhase2[celli] -= dmdt12*coeffs21Cell;
696 }
697 }
698 }
699
700 }
701
702 // Update ddtAlphaMax
703 this->ddtAlphaMax_ =
704 max(gMax((dmdt21*coeffs1)()), gMax((dmdt12*coeffs2)()));
705 }
706}
707
708
709template<class BasePhaseSystem>
711(
712 const phaseModel& phase,
715 const word speciesName
716)
717{
718 // Fill the volumetric mass transfer for species
719 forAllConstIters(massTransferModels_, iter)
720 {
721 if (iter()->transferSpecie() == speciesName)
722 {
723 // Explicit source
724 Su +=
725 this->Su()[phase.name()]
726 + this->Sp()[phase.name()]*phase.oldTime();
727 }
728 }
729}
730
731
732template<class BasePhaseSystem>
734{
735 bool includeVolChange(true);
736 forAllIters(massTransferModels_, iter)
737 {
738 if (!iter()->includeVolChange())
739 {
740 includeVolChange = false;
741 }
742 }
743 return includeVolChange;
744}
745
746// ************************************************************************* //
surfaceScalarField & phi
const volScalarField & alpha1
phaseModel & phase1
const volScalarField & alpha2
phaseModel & phase2
tmp< volScalarField > dmdt() const
Return the blended mass transfer rate.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
const GeometricField< Type, PatchField, GeoMesh > & oldTime() const
Return old time field.
bool set(const Key &key, T *ptr)
Assign a new entry, overwriting existing entries.
bool found(const Key &key) const
Return true if hashed entry is found in table.
Definition: HashTableI.H:100
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
static word groupName(StringType base, const word &group)
Create dot-delimited name.group string.
Class for mass transfer between phases.
virtual void correctMassSources(const volScalarField &T)
Correct/calculates mass sources dmdt for phases.
virtual void massSpeciesTransfer(const Foam::phaseModel &phase, volScalarField::Internal &Su, volScalarField::Internal &Sp, const word speciesName)
Calculate mass transfer for species.
dmdtTable dmdt_
Overall inter-phase mass transfer rates [Kg/s].
virtual tmp< fvScalarMatrix > volTransfer(const volScalarField &p)
Return the volumetric rate transfer matrix.
virtual bool includeVolChange()
Add volume change in pEq.
massTransferModelTable massTransferModels_
Mass transfer models.
virtual void alphaTransfer(SuSpTable &Su, SuSpTable &Sp)
Calculate mass transfer for alpha's.
tmp< volScalarField > calculateL(const volScalarField &dmdtNetki, const phasePairKey &keyik, const phasePairKey &keyki, const volScalarField &T) const
Calculate L between phases.
const HeatTransferModel< ThermoCloud< CloudType > > & heatTransfer() const
Return reference to heat transfer model.
Definition: ThermoCloudI.H:82
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
Single incompressible phase derived from the phase-fraction. Used as part of the multiPhaseMixture fo...
Definition: phaseModel.H:61
const word & name() const
Definition: phaseModel.H:144
An ordered or unorder pair of phase names. Typically specified as follows.
Definition: phasePairKey.H:68
Description for mass transfer between a pair of phases. The direction of the mass transfer is from th...
Definition: phasePair.H:56
const multiphaseInter::phaseModel & phase1() const
Definition: phasePairI.H:30
const multiphaseInter::phaseModel & phase2() const
Definition: phasePairI.H:36
Single incompressible phase derived from the phase-fraction. Used as part of the multiPhaseMixture fo...
Definition: phase.H:57
A class for managing temporary objects.
Definition: tmp.H:65
bool valid() const noexcept
Identical to good(), or bool operator.
Definition: tmp.H:292
T & ref() const
Definition: tmpI.H:227
A class for handling words, derived from Foam::string.
Definition: word.H:68
volScalarField & p
const volScalarField & T
dynamicFvMesh & mesh
Fundamental dimensioned constants.
Calculate the divergence of the given field.
Calculate the finiteVolume matrix for implicit and explicit sources.
zeroField Su
Definition: alphaSuSp.H:1
zeroField Sp
Definition: alphaSuSp.H:2
word timeName
Definition: getTimeIndex.H:3
label phasei
Definition: pEqn.H:27
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:49
zeroField Sp(const Foam::zero, const GeometricField< Type, fvPatchField, volMesh > &)
A no-op source.
const dimensionSet dimPressure
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
dimensionedScalar pos(const dimensionedScalar &ds)
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.
const dimensionSet dimEnergy
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
const dimensionSet dimTemperature(0, 0, 0, 1, 0, 0, 0)
Definition: dimensionSets.H:54
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
dimensionedScalar neg(const dimensionedScalar &ds)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
const dimensionSet dimVolume(pow3(dimLength))
Definition: dimensionSets.H:60
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
const dimensionSet dimDensity
Type gMax(const FieldField< Field, Type > &f)
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:51
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllIters(container, iter)
Iterate across all elements in the container object.
Definition: stdFoam.H:260
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278
const vector L(dict.get< vector >("L"))