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 | Copyright (C) 2013-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 
26 #include "twoPhaseSystem.H"
28 #include "BlendedInterfacialModel.H"
29 #include "virtualMassModel.H"
30 #include "heatTransferModel.H"
31 #include "liftModel.H"
32 #include "wallLubricationModel.H"
33 #include "turbulentDispersionModel.H"
34 #include "fvMatrix.H"
35 #include "surfaceInterpolate.H"
36 #include "MULES.H"
37 #include "subCycle.H"
38 #include "fvcDdt.H"
39 #include "fvcDiv.H"
40 #include "fvcSnGrad.H"
41 #include "fvcFlux.H"
42 #include "fvcCurl.H"
43 #include "fvmDdt.H"
44 #include "fvmLaplacian.H"
46 #include "blendingMethod.H"
47 #include "HashPtrTable.H"
48 
49 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
50 
52 (
53  const fvMesh& mesh,
54  const dimensionedVector& g
55 )
56 :
57  IOdictionary
58  (
59  IOobject
60  (
61  "phaseProperties",
62  mesh.time().constant(),
63  mesh,
64  IOobject::MUST_READ_IF_MODIFIED,
65  IOobject::NO_WRITE
66  )
67  ),
68 
69  mesh_(mesh),
70 
71  phase1_
72  (
73  *this,
74  *this,
75  wordList(lookup("phases"))[0]
76  ),
77 
78  phase2_
79  (
80  *this,
81  *this,
82  wordList(lookup("phases"))[1]
83  ),
84 
85  phi_
86  (
87  IOobject
88  (
89  "phi",
90  mesh.time().timeName(),
91  mesh,
92  IOobject::NO_READ,
93  IOobject::AUTO_WRITE
94  ),
95  this->calcPhi()
96  ),
97 
98  dgdt_
99  (
100  IOobject
101  (
102  "dgdt",
103  mesh.time().timeName(),
104  mesh,
105  IOobject::READ_IF_PRESENT,
106  IOobject::AUTO_WRITE
107  ),
108  mesh,
109  dimensionedScalar("dgdt", dimless/dimTime, 0)
110  )
111 {
112  phase2_.volScalarField::operator=(scalar(1) - phase1_);
113 
114 
115  // Blending
116  forAllConstIter(dictionary, subDict("blending"), iter)
117  {
118  blendingMethods_.insert
119  (
120  iter().dict().dictName(),
122  (
123  iter().dict(),
124  wordList(lookup("phases"))
125  )
126  );
127  }
128 
129 
130  // Pairs
131 
132  phasePair::scalarTable sigmaTable(lookup("sigma"));
133  phasePair::dictTable aspectRatioTable(lookup("aspectRatio"));
134 
135  pair_.set
136  (
137  new phasePair
138  (
139  phase1_,
140  phase2_,
141  g,
142  sigmaTable
143  )
144  );
145 
146  pair1In2_.set
147  (
148  new orderedPhasePair
149  (
150  phase1_,
151  phase2_,
152  g,
153  sigmaTable,
154  aspectRatioTable
155  )
156  );
157 
158  pair2In1_.set
159  (
160  new orderedPhasePair
161  (
162  phase2_,
163  phase1_,
164  g,
165  sigmaTable,
166  aspectRatioTable
167  )
168  );
169 
170 
171  // Models
172 
173  drag_.set
174  (
175  new BlendedInterfacialModel<dragModel>
176  (
177  lookup("drag"),
178  (
179  blendingMethods_.found("drag")
180  ? blendingMethods_["drag"]
181  : blendingMethods_["default"]
182  ),
183  pair_,
184  pair1In2_,
185  pair2In1_,
186  false // Do not zero drag coefficent at fixed-flux BCs
187  )
188  );
189 
190  virtualMass_.set
191  (
192  new BlendedInterfacialModel<virtualMassModel>
193  (
194  lookup("virtualMass"),
195  (
196  blendingMethods_.found("virtualMass")
197  ? blendingMethods_["virtualMass"]
198  : blendingMethods_["default"]
199  ),
200  pair_,
201  pair1In2_,
202  pair2In1_
203  )
204  );
205 
206  heatTransfer_.set
207  (
208  new BlendedInterfacialModel<heatTransferModel>
209  (
210  lookup("heatTransfer"),
211  (
212  blendingMethods_.found("heatTransfer")
213  ? blendingMethods_["heatTransfer"]
214  : blendingMethods_["default"]
215  ),
216  pair_,
217  pair1In2_,
218  pair2In1_
219  )
220  );
221 
222  lift_.set
223  (
224  new BlendedInterfacialModel<liftModel>
225  (
226  lookup("lift"),
227  (
228  blendingMethods_.found("lift")
229  ? blendingMethods_["lift"]
230  : blendingMethods_["default"]
231  ),
232  pair_,
233  pair1In2_,
234  pair2In1_
235  )
236  );
237 
238  wallLubrication_.set
239  (
240  new BlendedInterfacialModel<wallLubricationModel>
241  (
242  lookup("wallLubrication"),
243  (
244  blendingMethods_.found("wallLubrication")
245  ? blendingMethods_["wallLubrication"]
246  : blendingMethods_["default"]
247  ),
248  pair_,
249  pair1In2_,
250  pair2In1_
251  )
252  );
253 
254  turbulentDispersion_.set
255  (
256  new BlendedInterfacialModel<turbulentDispersionModel>
257  (
258  lookup("turbulentDispersion"),
259  (
260  blendingMethods_.found("turbulentDispersion")
261  ? blendingMethods_["turbulentDispersion"]
262  : blendingMethods_["default"]
263  ),
264  pair_,
265  pair1In2_,
266  pair2In1_
267  )
268  );
269 }
270 
271 
272 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
273 
275 {}
276 
277 
278 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
279 
281 {
282  return phase1_*phase1_.thermo().rho() + phase2_*phase2_.thermo().rho();
283 }
284 
285 
287 {
288  return phase1_*phase1_.U() + phase2_*phase2_.U();
289 }
290 
291 
292 Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::calcPhi() const
293 {
294  return
295  fvc::interpolate(phase1_)*phase1_.phi()
296  + fvc::interpolate(phase2_)*phase2_.phi();
297 }
298 
299 
301 {
302  return drag_->K();
303 }
304 
305 
307 {
308  return drag_->Kf();
309 }
310 
311 
313 {
314  return virtualMass_->K();
315 }
316 
317 
319 {
320  return virtualMass_->Kf();
321 }
322 
323 
325 {
326  return heatTransfer_->K();
327 }
328 
329 
331 {
332  return lift_->F<vector>() + wallLubrication_->F<vector>();
333 }
334 
335 
337 {
338  return lift_->Ff() + wallLubrication_->Ff();
339 }
340 
341 
343 {
344  return turbulentDispersion_->D();
345 }
346 
347 
349 {
350  const Time& runTime = mesh_.time();
351 
352  volScalarField& alpha1 = phase1_;
353  volScalarField& alpha2 = phase2_;
354 
355  const surfaceScalarField& phi1 = phase1_.phi();
356  const surfaceScalarField& phi2 = phase2_.phi();
357 
358  const dictionary& alphaControls = mesh_.solverDict
359  (
360  alpha1.name()
361  );
362 
363  label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
364  label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
365 
366  word alphaScheme("div(phi," + alpha1.name() + ')');
367  word alpharScheme("div(phir," + alpha1.name() + ')');
368 
369  alpha1.correctBoundaryConditions();
370 
371 
372  surfaceScalarField phic("phic", phi_);
373  surfaceScalarField phir("phir", phi1 - phi2);
374 
375  tmp<surfaceScalarField> alpha1alpha2f;
376 
377  if (pPrimeByA_.valid())
378  {
379  alpha1alpha2f =
380  fvc::interpolate(max(alpha1, scalar(0)))
381  *fvc::interpolate(max(alpha2, scalar(0)));
382 
383  surfaceScalarField phiP
384  (
385  pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
386  );
387 
388  phir += phiP;
389  }
390 
391  for (int acorr=0; acorr<nAlphaCorr; acorr++)
392  {
394  (
395  IOobject
396  (
397  "Sp",
398  runTime.timeName(),
399  mesh_
400  ),
401  mesh_,
402  dimensionedScalar("Sp", dgdt_.dimensions(), 0.0)
403  );
404 
406  (
407  IOobject
408  (
409  "Su",
410  runTime.timeName(),
411  mesh_
412  ),
413  // Divergence term is handled explicitly to be
414  // consistent with the explicit transport solution
415  fvc::div(phi_)*min(alpha1, scalar(1))
416  );
417 
418  forAll(dgdt_, celli)
419  {
420  if (dgdt_[celli] > 0.0)
421  {
422  Sp[celli] -= dgdt_[celli]/max(1.0 - alpha1[celli], 1e-4);
423  Su[celli] += dgdt_[celli]/max(1.0 - alpha1[celli], 1e-4);
424  }
425  else if (dgdt_[celli] < 0.0)
426  {
427  Sp[celli] += dgdt_[celli]/max(alpha1[celli], 1e-4);
428  }
429  }
430 
431  surfaceScalarField alphaPhic1
432  (
433  fvc::flux
434  (
435  phic,
436  alpha1,
437  alphaScheme
438  )
439  + fvc::flux
440  (
441  -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
442  alpha1,
444  )
445  );
446 
447  phase1_.correctInflowOutflow(alphaPhic1);
448 
449  if (nAlphaSubCycles > 1)
450  {
451  for
452  (
453  subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
454  !(++alphaSubCycle).end();
455  )
456  {
457  surfaceScalarField alphaPhic10(alphaPhic1);
458 
460  (
461  geometricOneField(),
462  alpha1,
463  phi_,
464  alphaPhic10,
465  (alphaSubCycle.index()*Sp)(),
466  (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
467  phase1_.alphaMax(),
468  0
469  );
470 
471  if (alphaSubCycle.index() == 1)
472  {
473  phase1_.alphaPhi() = alphaPhic10;
474  }
475  else
476  {
477  phase1_.alphaPhi() += alphaPhic10;
478  }
479  }
480 
481  phase1_.alphaPhi() /= nAlphaSubCycles;
482  }
483  else
484  {
486  (
487  geometricOneField(),
488  alpha1,
489  phi_,
490  alphaPhic1,
491  Sp,
492  Su,
493  phase1_.alphaMax(),
494  0
495  );
496 
497  phase1_.alphaPhi() = alphaPhic1;
498  }
499 
500  if (pPrimeByA_.valid())
501  {
502  fvScalarMatrix alpha1Eqn
503  (
504  fvm::ddt(alpha1) - fvc::ddt(alpha1)
505  - fvm::laplacian(alpha1alpha2f()*pPrimeByA_(), alpha1, "bounded")
506  );
507 
508  alpha1Eqn.relax();
509  alpha1Eqn.solve();
510 
511  phase1_.alphaPhi() += alpha1Eqn.flux();
512  }
513 
514  phase1_.alphaRhoPhi() =
515  fvc::interpolate(phase1_.rho())*phase1_.alphaPhi();
516 
517  phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
518  phase2_.correctInflowOutflow(phase2_.alphaPhi());
519  phase2_.alphaRhoPhi() =
520  fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();
521 
522  Info<< alpha1.name() << " volume fraction = "
523  << alpha1.weightedAverage(mesh_.V()).value()
524  << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
525  << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
526  << endl;
527 
528  // Ensure the phase-fractions are bounded
529  alpha1.max(0);
530  alpha1.min(1);
531 
532  alpha2 = scalar(1) - alpha1;
533  }
534 }
535 
536 
538 {
539  phase1_.correct();
540  phase2_.correct();
541 }
542 
543 
545 {
546  phase1_.turbulence().correct();
547  phase2_.turbulence().correct();
548 }
549 
550 
552 {
553  if (regIOobject::read())
554  {
555  bool readOK = true;
556 
557  readOK &= phase1_.read(*this);
558  readOK &= phase2_.read(*this);
559 
560  // models ...
561 
562  return readOK;
563  }
564  else
565  {
566  return false;
567  }
568 }
569 
570 
571 const Foam::dragModel& Foam::twoPhaseSystem::drag(const phaseModel& phase) const
572 {
573  return drag_->phaseModel(phase);
574 }
575 
576 
578 Foam::twoPhaseSystem::virtualMass(const phaseModel& phase) const
579 {
580  return virtualMass_->phaseModel(phase);
581 }
582 
583 
585 {
586  return pair_->sigma();
587 }
588 
589 
590 // ************************************************************************* //
tmp< volVectorField > F() const
Return the combined force (lift + wall-lubrication)
HashTable< dictionary, phasePairKey, phasePairKey::hash > dictTable
Dictionary hash table.
Definition: phasePair.H:59
dictionary dict
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
tmp< volScalarField > Kh() const
Return the heat transfer coefficient.
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
zeroField Su
Definition: alphaSuSp.H:1
void correct()
Correct two-phase properties other than turbulence.
const double e
Elementary charge.
Definition: doubleFloat.H:78
virtual bool read()
Read object.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
word alpharScheme("div(phirb,alpha)")
tmp< volScalarField > sigma() const
Return the surface tension coefficient.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
Calculate the matrix for the laplacian of the field.
Calculate the snGrad of the given volField.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Generic dimensioned Type class.
DimensionedField< scalar, volMesh > Internal
Type of the internal field from which this GeometricField is derived.
label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")))
tmp< volScalarField > Kd() const
Return the drag coefficient.
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
Calculate the first temporal derivative.
stressControl lookup("compactNormalStress") >> compactNormalStress
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
Calculate the curl of the given volField by constructing the Hodge-dual of the symmetric part of the ...
Calculate the face-flux of the given field.
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
Calulate the matrix for the first temporal derivative.
word timeName
Definition: getTimeIndex.H:3
twoPhaseSystem(const fvMesh &)
Construct from fvMesh.
const word dictName("particleTrackDict")
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
label readLabel(Istream &is)
Definition: label.H:64
static autoPtr< blendingMethod > New(const dictionary &dict, const wordList &phaseNames)
tmp< surfaceScalarField > Kdf() const
Return the face drag coefficient.
tmp< volScalarField > D() const
Return the turbulent diffusivity.
bool read()
Read base phaseProperties dictionary.
Calculate the divergence of the given field.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
volScalarField & alpha1
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< volVectorField > U() const
Return the mixture velocity.
void correctTurbulence()
Correct two-phase turbulence.
List< word > wordList
A List of words.
Definition: fileName.H:54
virtual void solve()
Solve for the phase fractions.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
surfaceScalarField phic(mixture.cAlpha() *mag(phi/mesh.magSf()))
virtual ~twoPhaseSystem()
Destructor.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:46
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
tmp< surfaceScalarField > Vmf() const
Return the face virtual mass coefficient.
tmp< surfaceScalarField > Ff() const
Return the combined face-force (lift + wall-lubrication)
messageStream Info
const dragModel & drag(const phaseModel &phase) const
Return the drag model for the given phase.
const virtualMassModel & virtualMass(const phaseModel &phase) const
Return the virtual mass model for the given phase.
MULES: Multidimensional universal limiter for explicit solution.
HashTable< scalar, phasePairKey, phasePairKey::hash > scalarTable
Scalar hash table.
Definition: phasePair.H:63
A class for managing temporary objects.
Definition: PtrList.H:53
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:32
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
tmp< volScalarField > rho() const
Return the mixture density.
tmp< volScalarField > Vm() const
Return the virtual mass coefficient.
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")))
zeroField Sp
Definition: alphaSuSp.H:2