twoPhaseSystem.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2013-2018 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 #include "UniformField.H"
49 
50 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
51 
53 (
54  const fvMesh& mesh,
55  const dimensionedVector& g
56 )
57 :
58  IOdictionary
59  (
60  IOobject
61  (
62  "phaseProperties",
63  mesh.time().constant(),
64  mesh,
65  IOobject::MUST_READ_IF_MODIFIED,
66  IOobject::NO_WRITE
67  )
68  ),
69 
70  mesh_(mesh),
71 
72  phase1_
73  (
74  *this,
75  *this,
76  wordList(lookup("phases"))[0]
77  ),
78 
79  phase2_
80  (
81  *this,
82  *this,
83  wordList(lookup("phases"))[1]
84  ),
85 
86  phi_
87  (
88  IOobject
89  (
90  "phi",
91  mesh.time().timeName(),
92  mesh,
93  IOobject::NO_READ,
94  IOobject::AUTO_WRITE
95  ),
96  this->calcPhi()
97  ),
98 
99  dgdt_
100  (
101  IOobject
102  (
103  "dgdt",
104  mesh.time().timeName(),
105  mesh,
106  IOobject::READ_IF_PRESENT,
107  IOobject::AUTO_WRITE
108  ),
109  mesh,
111  )
112 {
113  phase2_.volScalarField::operator=(scalar(1) - phase1_);
114 
115 
116  // Blending
117  forAllConstIter(dictionary, subDict("blending"), iter)
118  {
119  blendingMethods_.insert
120  (
121  iter().dict().dictName(),
123  (
124  iter().dict(),
125  wordList(lookup("phases"))
126  )
127  );
128  }
129 
130 
131  // Pairs
132 
133  phasePair::scalarTable sigmaTable(lookup("sigma"));
134  phasePair::dictTable aspectRatioTable(lookup("aspectRatio"));
135 
136  pair_.set
137  (
138  new phasePair
139  (
140  phase1_,
141  phase2_,
142  g,
143  sigmaTable
144  )
145  );
146 
147  pair1In2_.set
148  (
149  new orderedPhasePair
150  (
151  phase1_,
152  phase2_,
153  g,
154  sigmaTable,
155  aspectRatioTable
156  )
157  );
158 
159  pair2In1_.set
160  (
161  new orderedPhasePair
162  (
163  phase2_,
164  phase1_,
165  g,
166  sigmaTable,
167  aspectRatioTable
168  )
169  );
170 
171 
172  // Models
173 
174  drag_.set
175  (
176  new BlendedInterfacialModel<dragModel>
177  (
178  lookup("drag"),
179  (
180  blendingMethods_.found("drag")
181  ? blendingMethods_["drag"]
182  : blendingMethods_["default"]
183  ),
184  pair_,
185  pair1In2_,
186  pair2In1_,
187  false // Do not zero drag coefficient at fixed-flux BCs
188  )
189  );
190 
191  virtualMass_.set
192  (
193  new BlendedInterfacialModel<virtualMassModel>
194  (
195  lookup("virtualMass"),
196  (
197  blendingMethods_.found("virtualMass")
198  ? blendingMethods_["virtualMass"]
199  : blendingMethods_["default"]
200  ),
201  pair_,
202  pair1In2_,
203  pair2In1_
204  )
205  );
206 
207  heatTransfer_.set
208  (
209  new BlendedInterfacialModel<heatTransferModel>
210  (
211  lookup("heatTransfer"),
212  (
213  blendingMethods_.found("heatTransfer")
214  ? blendingMethods_["heatTransfer"]
215  : blendingMethods_["default"]
216  ),
217  pair_,
218  pair1In2_,
219  pair2In1_
220  )
221  );
222 
223  lift_.set
224  (
225  new BlendedInterfacialModel<liftModel>
226  (
227  lookup("lift"),
228  (
229  blendingMethods_.found("lift")
230  ? blendingMethods_["lift"]
231  : blendingMethods_["default"]
232  ),
233  pair_,
234  pair1In2_,
235  pair2In1_
236  )
237  );
238 
239  wallLubrication_.set
240  (
241  new BlendedInterfacialModel<wallLubricationModel>
242  (
243  lookup("wallLubrication"),
244  (
245  blendingMethods_.found("wallLubrication")
246  ? blendingMethods_["wallLubrication"]
247  : blendingMethods_["default"]
248  ),
249  pair_,
250  pair1In2_,
251  pair2In1_
252  )
253  );
254 
255  turbulentDispersion_.set
256  (
257  new BlendedInterfacialModel<turbulentDispersionModel>
258  (
259  lookup("turbulentDispersion"),
260  (
261  blendingMethods_.found("turbulentDispersion")
262  ? blendingMethods_["turbulentDispersion"]
263  : blendingMethods_["default"]
264  ),
265  pair_,
266  pair1In2_,
267  pair2In1_
268  )
269  );
270 }
271 
272 
273 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
274 
276 {}
277 
278 
279 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
280 
282 {
283  return phase1_*phase1_.thermo().rho() + phase2_*phase2_.thermo().rho();
284 }
285 
286 
288 {
289  return phase1_*phase1_.U() + phase2_*phase2_.U();
290 }
291 
292 
293 Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::calcPhi() const
294 {
295  return
296  fvc::interpolate(phase1_)*phase1_.phi()
297  + fvc::interpolate(phase2_)*phase2_.phi();
298 }
299 
300 
302 {
303  return drag_->K();
304 }
305 
306 
308 {
309  return drag_->Kf();
310 }
311 
312 
314 {
315  return virtualMass_->K();
316 }
317 
318 
320 {
321  return virtualMass_->Kf();
322 }
323 
324 
326 {
327  return heatTransfer_->K();
328 }
329 
330 
332 {
333  return lift_->F<vector>() + wallLubrication_->F<vector>();
334 }
335 
336 
338 {
339  return lift_->Ff() + wallLubrication_->Ff();
340 }
341 
342 
344 {
345  return turbulentDispersion_->D();
346 }
347 
348 
350 {
351  const Time& runTime = mesh_.time();
352 
353  volScalarField& alpha1 = phase1_;
354  volScalarField& alpha2 = phase2_;
355 
356  const surfaceScalarField& phi1 = phase1_.phi();
357  const surfaceScalarField& phi2 = phase2_.phi();
358 
359  const dictionary& alphaControls = mesh_.solverDict
360  (
361  alpha1.name()
362  );
363 
364  label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
365  label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
366 
367  word alphaScheme("div(phi," + alpha1.name() + ')');
368  word alpharScheme("div(phir," + alpha1.name() + ')');
369 
370  alpha1.correctBoundaryConditions();
371 
372 
373  surfaceScalarField phic("phic", phi_);
374  surfaceScalarField phir("phir", phi1 - phi2);
375 
376  tmp<surfaceScalarField> alpha1alpha2f;
377 
378  if (pPrimeByA_.valid())
379  {
380  alpha1alpha2f =
381  fvc::interpolate(max(alpha1, scalar(0)))
382  *fvc::interpolate(max(alpha2, scalar(0)));
383 
384  surfaceScalarField phiP
385  (
386  pPrimeByA_()*fvc::snGrad(alpha1, "bounded")*mesh_.magSf()
387  );
388 
389  phir += phiP;
390  }
391 
392  for (int acorr=0; acorr<nAlphaCorr; acorr++)
393  {
395  (
396  IOobject
397  (
398  "Sp",
399  runTime.timeName(),
400  mesh_
401  ),
402  mesh_,
403  dimensionedScalar(dgdt_.dimensions(), 0)
404  );
405 
407  (
408  IOobject
409  (
410  "Su",
411  runTime.timeName(),
412  mesh_
413  ),
414  // Divergence term is handled explicitly to be
415  // consistent with the explicit transport solution
416  fvc::div(phi_)*min(alpha1, scalar(1))
417  );
418 
419  forAll(dgdt_, celli)
420  {
421  if (dgdt_[celli] > 0.0)
422  {
423  Sp[celli] -= dgdt_[celli]/max(1.0 - alpha1[celli], 1e-4);
424  Su[celli] += dgdt_[celli]/max(1.0 - alpha1[celli], 1e-4);
425  }
426  else if (dgdt_[celli] < 0.0)
427  {
428  Sp[celli] += dgdt_[celli]/max(alpha1[celli], 1e-4);
429  }
430  }
431 
432  surfaceScalarField alphaPhic1
433  (
434  fvc::flux
435  (
436  phic,
437  alpha1,
438  alphaScheme
439  )
440  + fvc::flux
441  (
442  -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
443  alpha1,
445  )
446  );
447 
448  phase1_.correctInflowOutflow(alphaPhic1);
449 
450  if (nAlphaSubCycles > 1)
451  {
452  for
453  (
454  subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
455  !(++alphaSubCycle).end();
456  )
457  {
458  surfaceScalarField alphaPhic10(alphaPhic1);
459 
461  (
462  geometricOneField(),
463  alpha1,
464  phi_,
465  alphaPhic10,
466  (alphaSubCycle.index()*Sp)(),
467  (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
468  UniformField<scalar>(phase1_.alphaMax()),
469  zeroField()
470  );
471 
472  if (alphaSubCycle.index() == 1)
473  {
474  phase1_.alphaPhi() = alphaPhic10;
475  }
476  else
477  {
478  phase1_.alphaPhi() += alphaPhic10;
479  }
480  }
481 
482  phase1_.alphaPhi() /= nAlphaSubCycles;
483  }
484  else
485  {
487  (
488  geometricOneField(),
489  alpha1,
490  phi_,
491  alphaPhic1,
492  Sp,
493  Su,
494  UniformField<scalar>(phase1_.alphaMax()),
495  zeroField()
496  );
497 
498  phase1_.alphaPhi() = alphaPhic1;
499  }
500 
501  if (pPrimeByA_.valid())
502  {
503  fvScalarMatrix alpha1Eqn
504  (
505  fvm::ddt(alpha1) - fvc::ddt(alpha1)
506  - fvm::laplacian(alpha1alpha2f()*pPrimeByA_(), alpha1, "bounded")
507  );
508 
509  alpha1Eqn.relax();
510  alpha1Eqn.solve();
511 
512  phase1_.alphaPhi() += alpha1Eqn.flux();
513  }
514 
515  phase1_.alphaRhoPhi() =
516  fvc::interpolate(phase1_.rho())*phase1_.alphaPhi();
517 
518  phase2_.alphaPhi() = phi_ - phase1_.alphaPhi();
519  phase2_.correctInflowOutflow(phase2_.alphaPhi());
520  phase2_.alphaRhoPhi() =
521  fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();
522 
523  Info<< alpha1.name() << " volume fraction = "
524  << alpha1.weightedAverage(mesh_.V()).value()
525  << " Min(" << alpha1.name() << ") = " << min(alpha1).value()
526  << " Max(" << alpha1.name() << ") = " << max(alpha1).value()
527  << endl;
528 
529  // Ensure the phase-fractions are bounded
530  alpha1.max(0);
531  alpha1.min(1);
532 
533  alpha2 = scalar(1) - alpha1;
534  }
535 }
536 
537 
539 {
540  phase1_.correct();
541  phase2_.correct();
542 }
543 
544 
546 {
547  phase1_.turbulence().correct();
548  phase2_.turbulence().correct();
549 }
550 
551 
553 {
554  if (regIOobject::read())
555  {
556  bool readOK = true;
557 
558  readOK &= phase1_.read(*this);
559  readOK &= phase2_.read(*this);
560 
561  // models ...
562 
563  return readOK;
564  }
565  else
566  {
567  return false;
568  }
569 }
570 
571 
573 {
574  return pair_->sigma();
575 }
576 
577 
578 // ************************************************************************* //
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:434
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.
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:256
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.
static autoPtr< blendingMethod > New(const word &modelName, const dictionary &dict, const wordList &phaseNames)
const volScalarField & alpha1
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
Calculate 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
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 > &)
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
MULES: Multidimensional universal limiter for explicit solution.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
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