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-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 
26 #include "twoPhaseSystem.H"
27 #include "dragModel.H"
28 #include "virtualMassModel.H"
29 
30 #include "MULES.H"
31 #include "subCycle.H"
32 
33 #include "fvcDdt.H"
34 #include "fvcDiv.H"
35 #include "fvcSnGrad.H"
36 #include "fvcFlux.H"
37 #include "fvcSup.H"
38 
39 #include "fvmDdt.H"
40 #include "fvmLaplacian.H"
41 #include "fvmSup.H"
42 
43 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44 
45 namespace Foam
46 {
47  defineTypeNameAndDebug(twoPhaseSystem, 0);
48  defineRunTimeSelectionTable(twoPhaseSystem, dictionary);
49 }
50 
51 
52 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
53 
55 (
56  const fvMesh& mesh
57 )
58 :
59  phaseSystem(mesh),
60  phase1_(phaseModels_[0]),
61  phase2_(phaseModels_[1])
62 {
63  phase2_.volScalarField::operator=(scalar(1) - phase1_);
64 
65  volScalarField& alpha1 = phase1_;
66  mesh.setFluxRequired(alpha1.name());
67 }
68 
69 
70 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
71 
73 {}
74 
75 
76 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
77 
80 {
81  return sigma
82  (
83  phasePairKey(phase1().name(), phase2().name())
84  );
85 }
86 
87 
88 const Foam::dragModel& Foam::twoPhaseSystem::drag(const phaseModel& phase) const
89 {
90  return lookupSubModel<dragModel>(phase, otherPhase(phase));
91 }
92 
93 
96 {
97  return Kd
98  (
99  phasePairKey(phase1().name(), phase2().name())
100  );
101 }
102 
103 
106 {
107  return Kdf
108  (
109  phasePairKey(phase1().name(), phase2().name())
110  );
111 }
112 
113 
115 Foam::twoPhaseSystem::virtualMass(const phaseModel& phase) const
116 {
117  return lookupSubModel<virtualMassModel>(phase, otherPhase(phase));
118 }
119 
120 
123 {
124  return Vm
125  (
126  phasePairKey(phase1().name(), phase2().name())
127  );
128 }
129 
130 
133 {
134  return Vmf
135  (
136  phasePairKey(phase1().name(), phase2().name())
137  );
138 }
139 
140 
143 {
144  return F
145  (
146  phasePairKey(phase1().name(), phase2().name())
147  );
148 }
149 
150 
153 {
154  return Ff
155  (
156  phasePairKey(phase1().name(), phase2().name())
157  );
158 }
159 
160 
163 {
164  return D
165  (
166  phasePairKey(phase1().name(), phase2().name())
167  );
168 }
169 
170 
172 {
173  return transfersMass(phase1());
174 }
175 
176 
179 {
180  return dmdt
181  (
182  phasePairKey(phase1().name(), phase2().name())
183  );
184 }
185 
186 
188 {
189  const Time& runTime = mesh_.time();
190 
191  volScalarField& alpha1 = phase1_;
192  volScalarField& alpha2 = phase2_;
193 
194  const dictionary& alphaControls = mesh_.solverDict(alpha1.name());
195 
196  label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
197  label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")));
198 
199  bool LTS = fv::localEulerDdt::enabled(mesh_);
200 
201  word alphaScheme("div(phi," + alpha1.name() + ')');
202  word alpharScheme("div(phir," + alpha1.name() + ')');
203 
204  const surfaceScalarField& phi = this->phi();
205  const surfaceScalarField& phi1 = phase1_.phi();
206  const surfaceScalarField& phi2 = phase2_.phi();
207 
208  // Construct the dilatation rate source term
209  tmp<volScalarField::DimensionedInternalField> tdgdt;
210 
211  if (phase1_.divU().valid() && phase2_.divU().valid())
212  {
213  tdgdt =
214  (
215  alpha2.dimensionedInternalField()
216  *phase1_.divU()().dimensionedInternalField()
217  - alpha1.dimensionedInternalField()
218  *phase2_.divU()().dimensionedInternalField()
219  );
220  }
221  else if (phase1_.divU().valid())
222  {
223  tdgdt =
224  (
225  alpha2.dimensionedInternalField()
226  *phase1_.divU()().dimensionedInternalField()
227  );
228  }
229  else if (phase2_.divU().valid())
230  {
231  tdgdt =
232  (
233  - alpha1.dimensionedInternalField()
234  *phase2_.divU()().dimensionedInternalField()
235  );
236  }
237 
238  alpha1.correctBoundaryConditions();
239  surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, scalar(0))));
240 
241  surfaceScalarField phic("phic", phi);
242  surfaceScalarField phir("phir", phi1 - phi2);
243 
244  tmp<surfaceScalarField> alphaDbyA;
245 
246  if (notNull(phase1_.DbyA()) && notNull(phase2_.DbyA()))
247  {
248  surfaceScalarField DbyA(phase1_.DbyA() + phase2_.DbyA());
249 
250  alphaDbyA =
251  fvc::interpolate(max(alpha1, scalar(0)))
252  *fvc::interpolate(max(alpha2, scalar(0)))
253  *DbyA;
254 
255  phir += DbyA*fvc::snGrad(alpha1, "bounded")*mesh_.magSf();
256  }
257 
258  for (int acorr=0; acorr<nAlphaCorr; acorr++)
259  {
261  (
262  IOobject
263  (
264  "Sp",
265  runTime.timeName(),
266  mesh_
267  ),
268  mesh_,
269  dimensionedScalar("Sp", dimless/dimTime, 0.0)
270  );
271 
273  (
274  IOobject
275  (
276  "Su",
277  runTime.timeName(),
278  mesh_
279  ),
280  // Divergence term is handled explicitly to be
281  // consistent with the explicit transport solution
282  fvc::div(phi)*min(alpha1, scalar(1))
283  );
284 
285  if (tdgdt.valid())
286  {
287  scalarField& dgdt = tdgdt();
288 
289  forAll(dgdt, celli)
290  {
291  if (dgdt[celli] > 0.0)
292  {
293  Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
294  Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
295  }
296  else if (dgdt[celli] < 0.0)
297  {
298  Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
299  }
300  }
301  }
302 
303  surfaceScalarField alphaPhic1
304  (
305  fvc::flux
306  (
307  phic,
308  alpha1,
309  alphaScheme
310  )
311  + fvc::flux
312  (
313  -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
314  alpha1,
316  )
317  );
318 
319  // Ensure that the flux at inflow BCs is preserved
320  forAll(alphaPhic1.boundaryField(), patchi)
321  {
322  fvsPatchScalarField& alphaPhic1p =
323  alphaPhic1.boundaryField()[patchi];
324 
325  if (!alphaPhic1p.coupled())
326  {
327  const scalarField& phi1p = phi1.boundaryField()[patchi];
328  const scalarField& alpha1p = alpha1.boundaryField()[patchi];
329 
330  forAll(alphaPhic1p, facei)
331  {
332  if (phi1p[facei] < 0)
333  {
334  alphaPhic1p[facei] = alpha1p[facei]*phi1p[facei];
335  }
336  }
337  }
338  }
339 
340  if (nAlphaSubCycles > 1)
341  {
342  tmp<volScalarField> trSubDeltaT;
343 
344  if (LTS)
345  {
346  trSubDeltaT =
348  }
349 
350  for
351  (
352  subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
353  !(++alphaSubCycle).end();
354  )
355  {
356  surfaceScalarField alphaPhic10(alphaPhic1);
357 
359  (
360  geometricOneField(),
361  alpha1,
362  phi,
363  alphaPhic10,
364  (alphaSubCycle.index()*Sp)(),
365  (Su - (alphaSubCycle.index() - 1)*Sp*alpha1)(),
366  phase1_.alphaMax(),
367  0
368  );
369 
370  if (alphaSubCycle.index() == 1)
371  {
372  phase1_.alphaPhi() = alphaPhic10;
373  }
374  else
375  {
376  phase1_.alphaPhi() += alphaPhic10;
377  }
378  }
379 
380  phase1_.alphaPhi() /= nAlphaSubCycles;
381  }
382  else
383  {
385  (
386  geometricOneField(),
387  alpha1,
388  phi,
389  alphaPhic1,
390  Sp,
391  Su,
392  phase1_.alphaMax(),
393  0
394  );
395 
396  phase1_.alphaPhi() = alphaPhic1;
397  }
398 
399  if (alphaDbyA.valid())
400  {
401  fvScalarMatrix alpha1Eqn
402  (
403  fvm::ddt(alpha1) - fvc::ddt(alpha1)
404  - fvm::laplacian(alphaDbyA, alpha1, "bounded")
405  );
406 
407  alpha1Eqn.relax();
408  alpha1Eqn.solve();
409 
410  phase1_.alphaPhi() += alpha1Eqn.flux();
411  }
412 
413  phase1_.alphaRhoPhi() =
414  fvc::interpolate(phase1_.rho())*phase1_.alphaPhi();
415 
416  phase2_.alphaPhi() = phi - phase1_.alphaPhi();
417  alpha2 = scalar(1) - alpha1;
418  phase2_.alphaRhoPhi() =
419  fvc::interpolate(phase2_.rho())*phase2_.alphaPhi();
420 
421  Info<< alpha1.name() << " volume fraction = "
422  << alpha1.weightedAverage(mesh_.V()).value()
423  << " Min(alpha1) = " << min(alpha1).value()
424  << " Max(alpha1) = " << max(alpha1).value()
425  << endl;
426  }
427 }
428 
429 
430 // ************************************************************************* //
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &vf, const surfaceScalarField &faceFlux, Istream &schemeData)
phaseModel & phase1
Definition: createFields.H:12
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
const virtualMassModel & virtualMass(const phaseModel &phase) const
Return the virtual mass model for the given phase.
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
phaseModel & phase2
Definition: createFields.H:13
Calculate the snGrad of the given volField.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > flux(const surfaceScalarField &phi, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcFlux.C:45
tmp< surfaceScalarField > Kdf() const
Return the face drag coefficient.
volScalarField Kd(fluid.Kd())
Calculate the matrix for implicit and explicit sources.
const phaseModel & phase() const
Definition: IATEsource.H:138
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
surfaceScalarField phir(phic *interface.nHatf())
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
tmp< fvMatrix< Type > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmLaplacian.C:46
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
twoPhaseSystem(const fvMesh &)
Construct from fvMesh.
Calulate the matrix for the first temporal derivative.
bool notNull(const T &t)
Return true if t is not a reference to the nullObject of type T.
Definition: nullObjectI.H:46
messageStream Info
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
MULES: Multidimensional universal limiter for explicit solution.
Namespace for OpenFOAM.
fvsPatchField< scalar > fvsPatchScalarField
label readLabel(Istream &is)
Definition: label.H:64
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
tmp< volScalarField > Kd() const
Return the drag coefficient.
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:44
virtual void solve()
Solve for the phase fractions.
surfaceScalarField Kdf("Kdf", fluid.Kdf())
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
Calculate the first temporal derivative.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
label nAlphaCorr(readLabel(alphaControls.lookup("nAlphaCorr")))
Calculate the field for explicit evaluation of implicit and explicit sources.
tmp< volVectorField > F() const
Return the combined force (lift + wall-lubrication)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
static tmp< volScalarField > localRSubDeltaT(const fvMesh &mesh, const label nAlphaSubCycles)
Calculate and return the reciprocal of the local sub-cycling.
Definition: localEulerDdt.C:57
#define forAll(list, i)
Definition: UList.H:421
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
tmp< surfaceScalarField > Vmf() const
Return the face virtual mass coefficient.
tmp< volScalarField > Vm() const
Return the virtual mass coefficient.
label patchi
word alpharScheme("div(phirb,alpha)")
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
volScalarField & alpha1
Definition: createFields.H:15
const phaseModel & otherPhase() const
Definition: IATEsource.H:148
rDeltaT dimensionedInternalField()
Calculate the matrix for the laplacian of the field.
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")))
Calculate the divergence of the given field.
void correctBoundaryConditions()
Correct boundary field.
virtual ~twoPhaseSystem()
Destructor.
DimensionedField< scalar, volMesh > DimensionedInternalField
tmp< volScalarField > D() const
Return the turbulent diffusivity.
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
tmp< surfaceScalarField > Ff() const
Return the combined face-force (lift + wall-lubrication)
bool transfersMass() const
Return true if there is mass transfer.
surfaceScalarField Ff(fluid.Ff())
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const dimensionedScalar F
Faraday constant: default SI units: [C/mol].
const dragModel & drag(const phaseModel &phase) const
Return the drag model for the given phase.
surfaceScalarField Vmf("Vmf", fluid.Vmf())
tmp< volScalarField > dmdt() const
Return the interfacial mass flow rate.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A class for managing temporary objects.
Definition: PtrList.H:118
tmp< volScalarField > sigma() const
Return the surface tension coefficient.
static bool enabled(const fvMesh &mesh)
Return true if LTS is enabled.
Definition: localEulerDdt.C:36
defineTypeNameAndDebug(combustionModel, 0)
Calculate the face-flux of the given field.