pEqn.H
Go to the documentation of this file.
1 PtrList<surfaceScalarField> alphafs(phases.size());
2 PtrList<volScalarField> rAUs(phases.size());
3 PtrList<surfaceScalarField> alpharAUfs(phases.size());
4 
6 {
7  phaseModel& phase = phases[phasei];
8  const volScalarField& alpha = phase;
9 
10  alphafs.set(phasei, fvc::interpolate(alpha).ptr());
11  alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
12 
13  rAUs.set
14  (
15  phasei,
16  new volScalarField
17  (
18  IOobject::groupName("rAU", phase.name()),
19  1.0
20  /(
21  UEqns[phasei].A()
22  + max(phase.residualAlpha() - alpha, scalar(0))
23  *phase.rho()/runTime.deltaT()
24  )
25  )
26  );
27 
28  alpharAUfs.set
29  (
30  phasei,
31  (
32  fvc::interpolate(max(alpha, phase.residualAlpha())*rAUs[phasei])
33  ).ptr()
34  );
35 }
36 
37 // Lift, wall-lubrication and turbulent diffusion fluxes
38 PtrList<surfaceScalarField> phiFs(phases.size());
39 {
40  autoPtr<PtrList<volVectorField>> Fs = fluid.Fs();
41 
43  {
44  phaseModel& phase = phases[phasei];
45 
46  if (Fs().set(phasei))
47  {
48  phiFs.set
49  (
50  phasei,
52  (
53  IOobject::groupName("phiF", phase.name()),
54  fvc::flux(rAUs[phasei]*Fs()[phasei])
55  )
56  );
57  }
58  }
59 }
60 {
61  autoPtr<PtrList<surfaceScalarField>> phiDs = fluid.phiDs(rAUs);
62 
64  {
65  phaseModel& phase = phases[phasei];
66 
67  if (phiDs().set(phasei))
68  {
69  if (phiFs.set(phasei))
70  {
71  phiFs[phasei] += phiDs()[phasei];
72  }
73  else
74  {
75  phiFs.set
76  (
77  phasei,
79  (
80  IOobject::groupName("phiF", phase.name()),
81  phiDs()[phasei]
82  )
83  );
84  }
85  }
86  }
87 }
88 
89 // --- Pressure corrector loop
90 while (pimple.correct())
91 {
92  // Update continuity errors due to temperature changes
93  fluid.correct();
94 
95  volScalarField rho("rho", fluid.rho());
96 
97  // Correct p_rgh for consistency with p and the updated densities
98  p_rgh = p - rho*gh;
99 
100  PtrList<volVectorField> HbyAs(phases.size());
101 
103  {
104  phaseModel& phase = phases[phasei];
105  const volScalarField& alpha = phase;
106 
107  // Correct fixed-flux BCs to be consistent with the velocity BCs
108  MRF.correctBoundaryFlux(phase.U(), phase.phi());
109 
110  HbyAs.set
111  (
112  phasei,
113  new volVectorField
114  (
115  IOobject::groupName("HbyA", phase.name()),
116  phase.U()
117  )
118  );
119 
120  HbyAs[phasei] =
121  rAUs[phasei]
122  *(
123  UEqns[phasei].H()
124  + max(phase.residualAlpha() - alpha, scalar(0))
125  *phase.rho()*phase.U().oldTime()/runTime.deltaT()
126  );
127  }
128 
130  (
131  "ghSnGradRho",
132  ghf*fvc::snGrad(rho)*mesh.magSf()
133  );
134 
135  PtrList<surfaceScalarField> phigFs(phases.size());
137  {
138  phaseModel& phase = phases[phasei];
139 
140  phigFs.set
141  (
142  phasei,
143  (
145  *(
147  - (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
148  - fluid.surfaceTension(phase)*mesh.magSf()
149  )
150  ).ptr()
151  );
152 
153  if (phiFs.set(phasei))
154  {
155  phigFs[phasei] += phiFs[phasei];
156  }
157  }
158 
159  PtrList<surfaceScalarField> phiHbyAs(phases.size());
160 
162  (
163  IOobject
164  (
165  "phiHbyA",
166  runTime.timeName(),
167  mesh,
168  IOobject::NO_READ,
169  IOobject::AUTO_WRITE
170  ),
171  mesh,
172  dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
173  );
174 
176  {
177  phaseModel& phase = phases[phasei];
178  const volScalarField& alpha = phase;
179 
180  // ddtPhiCorr filter -- only apply in pure(ish) phases
181  surfaceScalarField alphafBar
182  (
184  );
185  surfaceScalarField phiCorrCoeff(pos(alphafBar - 0.99));
186 
187  surfaceScalarField::Boundary& phiCorrCoeffBf =
188  phiCorrCoeff.boundaryFieldRef();
189 
190  forAll(mesh.boundary(), patchi)
191  {
192  // Set ddtPhiCorr to 0 on non-coupled boundaries
193  if
194  (
195  !mesh.boundary()[patchi].coupled()
196  || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
197  )
198  {
199  phiCorrCoeffBf[patchi] = 0;
200  }
201  }
202 
203  phiHbyAs.set
204  (
205  phasei,
207  (
208  IOobject::groupName("phiHbyA", phase.name()),
209  fvc::flux(HbyAs[phasei])
210  + phiCorrCoeff
212  (
213  alpha.oldTime()*phase.rho()().oldTime()*rAUs[phasei]
214  )
215  *(
216  MRF.absolute(phase.phi().oldTime())
217  - fvc::flux(phase.U().oldTime())
218  )/runTime.deltaT()
219  - phigFs[phasei]
220  )
221  );
222 
224  (
225  phaseSystem::KdTable,
226  fluid.Kds(),
227  KdIter
228  )
229  {
230  const volScalarField& K(*KdIter());
231 
232  const phasePair& pair(fluid.phasePairs()[KdIter.key()]);
233 
234  const phaseModel* phase1 = &pair.phase1();
235  const phaseModel* phase2 = &pair.phase2();
236 
237  forAllConstIter(phasePair, pair, iter)
238  {
239  if (phase1 == &phase)
240  {
241  phiHbyAs[phasei] +=
242  fvc::interpolate(rAUs[phasei]*K)
243  *MRF.absolute(phase2->phi());
244 
245  HbyAs[phasei] += rAUs[phasei]*K*phase2->U();
246  }
247 
248  Swap(phase1, phase2);
249  }
250  }
251 
253  }
254 
255  MRF.makeRelative(phiHbyA);
256 
257  // Construct pressure "diffusivity"
259  (
260  IOobject
261  (
262  "rAUf",
263  runTime.timeName(),
264  mesh
265  ),
266  mesh,
267  dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
268  );
269 
271  {
272  rAUf += alphafs[phasei]*alpharAUfs[phasei];
273  }
274  rAUf = mag(rAUf);
275 
276 
277  // Update the fixedFluxPressure BCs to ensure flux consistency
278  {
279  surfaceScalarField::Boundary phib(phi.boundaryField());
280  phib = 0;
282  {
283  phaseModel& phase = phases[phasei];
284  phib += alphafs[phasei].boundaryField()*phase.phi().boundaryField();
285  }
286 
288  (
289  p_rgh.boundaryFieldRef(),
290  (
291  phiHbyA.boundaryField() - phib
292  )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
293  );
294  }
295 
296  PtrList<fvScalarMatrix> pEqnComps(phases.size());
298  {
299  phaseModel& phase = phases[phasei];
300 
301  if (phase.compressible())
302  {
303  const volScalarField& alpha = phase;
304  const volScalarField& rho = phase.rho();
305 
306  if (pimple.transonic())
307  {
309  (
310  IOobject::groupName("phid", phase.name()),
311  fvc::interpolate(phase.thermo().psi())*phase.phi()
312  );
313 
314  pEqnComps.set
315  (
316  phasei,
317  (
318  (
319  phase.continuityError()
320  - fvc::Sp
321  (
322  fvc::ddt(alpha) + fvc::div(phase.alphaPhi()),
323  rho
324  )
325  )/rho
326  + correction
327  (
328  (alpha/rho)*
329  (
330  phase.thermo().psi()*fvm::ddt(p_rgh)
331  + fvm::div(phid, p_rgh)
332  - fvm::Sp(fvc::div(phid), p_rgh)
333  )
334  )
335  ).ptr()
336  );
337 
339  (
340  pEqnComps[phasei].faceFluxCorrectionPtr()
341  );
342  pEqnComps[phasei].relax();
343  }
344  else
345  {
346  pEqnComps.set
347  (
348  phasei,
349  (
350  (
351  phase.continuityError()
352  - fvc::Sp
353  (
354  (fvc::ddt(alpha) + fvc::div(phase.alphaPhi())),
355  rho
356  )
357  )/rho
358  + (alpha*phase.thermo().psi()/rho)
359  *correction(fvm::ddt(p_rgh))
360  ).ptr()
361  );
362  }
363 
364  if (fluid.transfersMass(phase))
365  {
366  if (pEqnComps.set(phasei))
367  {
368  pEqnComps[phasei] -= fluid.dmdt(phase)/rho;
369  }
370  else
371  {
372  pEqnComps.set
373  (
374  phasei,
375  fvm::Su(-fluid.dmdt(phase)/rho, p_rgh)
376  );
377  }
378  }
379  }
380  }
381 
382  // Cache p prior to solve for density update
383  volScalarField p_rgh_0(p_rgh);
384 
385  // Iterate over the pressure equation to correct for non-orthogonality
386  while (pimple.correctNonOrthogonal())
387  {
388  // Construct the transport part of the pressure equation
389  fvScalarMatrix pEqnIncomp
390  (
392  - fvm::laplacian(rAUf, p_rgh)
393  );
394 
395  {
396  fvScalarMatrix pEqn(pEqnIncomp);
397 
399  {
400  if (pEqnComps.set(phasei))
401  {
402  pEqn += pEqnComps[phasei];
403  }
404  }
405 
406  solve
407  (
408  pEqn,
409  mesh.solver(p_rgh.select(pimple.finalInnerIter()))
410  );
411  }
412 
413  // Correct fluxes and velocities on last non-orthogonal iteration
414  if (pimple.finalNonOrthogonalIter())
415  {
416  phi = phiHbyA + pEqnIncomp.flux();
417 
418  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
419 
421  {
422  phaseModel& phase = phases[phasei];
423 
424  phase.phi() = phiHbyAs[phasei] + alpharAUfs[phasei]*mSfGradp;
425 
426  // Set the phase dilatation rates
427  if (phase.compressible())
428  {
429  phase.divU(-pEqnComps[phasei] & p_rgh);
430  }
431  }
432 
433  // Optionally relax pressure for velocity correction
434  p_rgh.relax();
435 
436  mSfGradp = pEqnIncomp.flux()/rAUf;
437 
439  {
440  phaseModel& phase = phases[phasei];
441 
442  phase.U() =
443  HbyAs[phasei]
445  (
446  alpharAUfs[phasei]*mSfGradp
447  - phigFs[phasei]
448  );
449  phase.U().correctBoundaryConditions();
450  fvOptions.correct(phase.U());
451  }
452  }
453  }
454 
455  // Update and limit the static pressure
456  p = max(p_rgh + rho*gh, pMin);
457 
458  // Limit p_rgh
459  p_rgh = p - rho*gh;
460 
461  // Update densities from change in p_rgh
463  {
464  phaseModel& phase = phases[phasei];
465  phase.thermo().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
466  }
467 
468  // Correct p_rgh for consistency with p and the updated densities
469  rho = fluid.rho();
470  p_rgh = p - rho*gh;
471  p_rgh.correctBoundaryConditions();
472 }
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix.
K
Definition: pEqn.H:86
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
rho oldTime()
multiphaseSystem & fluid
Definition: createFields.H:10
p
Definition: pEqn.H:50
label phasei
Definition: pEqn.H:27
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
phiHbyA
Definition: pEqn.H:20
multiphaseSystem::phaseModelList & phases
Definition: createFields.H:11
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
const dictionary & pimple
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
tmp< surfaceScalarField > interpolate(const RhoType &rho)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
forAll(phases, phasei)
Definition: pEqn.H:5
fv::options & fvOptions
const surfaceScalarField & ghf
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
dimensionedScalar pos(const dimensionedScalar &ds)
dynamicFvMesh & mesh
void Swap(T &a, T &b)
Definition: Swap.H:43
rhoEqn solve()
IOMRFZoneList & MRF
PtrList< volScalarField > rAUs(fluid.phases().size())
PtrList< surfaceScalarField > alpharAUfs(phases.size())
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho)*mesh.magSf())
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
surfaceScalarField phid("phid", fvc::interpolate(psi)*(fvc::flux(HbyA)+rhorAUf *fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)))
phaseModel & phase1
p_rgh
Definition: pEqn.H:120
const dimensionedVector & g
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryFieldRef(),(phiHbyA.boundaryField()-MRF.relative(phib))/(mesh.magSf().boundaryField()*rAUf.boundaryField()))
dimensionedScalar pMin("pMin", dimPressure, fluid)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
PtrList< volVectorField > HbyAs(fluid.phases().size())
PtrList< surfaceScalarField > alphafs(phases.size())
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
label patchi
PtrList< fvVectorMatrix > UEqns(fluid.phases().size())
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const volScalarField & gh
volScalarField p_rgh_0(p_rgh)
PtrList< surfaceScalarField > phiHbyAs(fluid.phases().size())
tmp< GeometricField< Type, fvPatchField, volMesh > > average(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Area-weighted average a surfaceField creating a volField.
Definition: fvcAverage.C:46
dimensioned< scalar > mag(const dimensioned< Type > &)
phi
Definition: pEqn.H:18
phib
Definition: pEqn.H:191
phaseModel & phase2
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
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:44
rho
Definition: pEqn.H:1
void deleteDemandDrivenData(DataPtr &dataPtr)
PtrList< surfaceScalarField > phiFs(phases.size())
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
const dimensionSet dimArea(sqr(dimLength))
Definition: dimensionSets.H:57
const dimensionSet dimVelocity