pEqn.H
Go to the documentation of this file.
1 // Face volume fractions
2 PtrList<surfaceScalarField> alphafs(phases.size());
3 PtrList<surfaceScalarField> alphaRho0fs(phases.size());
5 {
6  phaseModel& phase = phases[phasei];
7  const volScalarField& alpha = phase;
8 
9  alphafs.set(phasei, fvc::interpolate(alpha).ptr());
10  alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
11 
12  alphaRho0fs.set
13  (
14  phasei,
15  (
17  (
18  max(alpha.oldTime(), phase.residualAlpha())
19  *phase.rho()().oldTime()
20  )
21  ).ptr()
22  );
23 }
24 
25 // Diagonal coefficients
26 PtrList<surfaceScalarField> rAUfs(phases.size());
27 {
28  PtrList<surfaceScalarField> AFfs(fluid.AFfs());
29 
30  forAll(fluid.movingPhases(), movingPhasei)
31  {
32  phaseModel& phase = fluid.movingPhases()[movingPhasei];
33 
34  rAUfs.set
35  (
36  phase.index(),
38  (
39  IOobject::groupName("rAUf", phase.name()),
40  1.0
41  /(
42  byDt(alphaRho0fs[phase.index()])
43  + fvc::interpolate(UEqns[phase.index()].A())
44  + AFfs[phase.index()]
45  )
46  )
47  );
48  }
49 }
50 fluid.fillFields("rAUf", dimTime/dimDensity, rAUfs);
51 
52 // Phase diagonal coefficients
53 PtrList<surfaceScalarField> alpharAUfs(phases.size());
55 {
56  phaseModel& phase = phases[phasei];
57  alpharAUfs.set
58  (
59  phase.index(),
60  (
61  max(alphafs[phase.index()], phase.residualAlpha())
62  *rAUfs[phase.index()]
63  ).ptr()
64  );
65 }
66 
67 // Explicit force fluxes
68 PtrList<surfaceScalarField> phiFfs(fluid.phiFfs(rAUfs));
69 
70 // --- Pressure corrector loop
71 while (pimple.correct())
72 {
73  volScalarField rho("rho", fluid.rho());
74 
75  // Correct p_rgh for consistency with p and the updated densities
76  p_rgh = p - rho*gh;
77 
78  // Correct fixed-flux BCs to be consistent with the velocity BCs
79  forAll(fluid.movingPhases(), movingPhasei)
80  {
81  phaseModel& phase = fluid.movingPhases()[movingPhasei];
82  MRF.correctBoundaryFlux(phase.U(), phase.phiRef());
83  }
84 
85  // Combined buoyancy and force fluxes
86  PtrList<surfaceScalarField> phigFs(phases.size());
87  {
89  (
90  "ghSnGradRho",
91  ghf*fvc::snGrad(rho)*mesh.magSf()
92  );
93 
95  {
96  phaseModel& phase = phases[phasei];
97 
98  phigFs.set
99  (
100  phasei,
101  (
103  *(
105  - (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
106  - fluid.surfaceTension(phase)*mesh.magSf()
107  )
108  ).ptr()
109  );
110 
111  if (phiFfs.set(phasei))
112  {
113  phigFs[phasei] += phiFfs[phasei];
114  }
115  }
116  }
117 
118  // Predicted fluxes for each phase
119  PtrList<surfaceScalarField> phiHbyAs(phases.size());
120  forAll(fluid.movingPhases(), movingPhasei)
121  {
122  phaseModel& phase = fluid.movingPhases()[movingPhasei];
123 
124  phiHbyAs.set
125  (
126  phase.index(),
128  (
129  IOobject::groupName("phiHbyA", phase.name()),
130  rAUfs[phase.index()]
131  *(
132  fvc::flux(UEqns[phase.index()].H())
133  + alphaRho0fs[phase.index()]
134  *byDt(MRF.absolute(phase.phi()().oldTime()))
135  )
136  - phigFs[phase.index()]
137  )
138  );
139  }
140  fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs);
141 
142  // Add explicit drag forces and fluxes if not doing partial elimination
143  if (!partialElimination)
144  {
145  PtrList<surfaceScalarField> phiKdPhifs(fluid.phiKdPhifs(rAUfs));
146 
148  {
149  if (phiKdPhifs.set(phasei))
150  {
151  phiHbyAs[phasei] -= phiKdPhifs[phasei];
152  }
153  }
154  }
155 
156  // Total predicted flux
158  (
159  IOobject
160  (
161  "phiHbyA",
162  runTime.timeName(),
163  mesh,
164  IOobject::NO_READ,
165  IOobject::AUTO_WRITE
166  ),
167  mesh,
168  dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
169  );
170 
172  {
174  }
175 
176  // Add explicit drag fluxes if doing partial elimination
177  if (partialElimination)
178  {
179  PtrList<surfaceScalarField> phiKdPhifs(fluid.phiKdPhifs(rAUfs));
180 
182  {
183  if (phiKdPhifs.set(phasei))
184  {
185  phiHbyA -= alphafs[phasei]*phiKdPhifs[phasei];
186  }
187  }
188  }
189 
190  MRF.makeRelative(phiHbyA);
191 
192  // Construct pressure "diffusivity"
194  (
195  IOobject
196  (
197  "rAUf",
198  runTime.timeName(),
199  mesh
200  ),
201  mesh,
202  dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
203  );
204 
206  {
207  rAUf += alphafs[phasei]*alpharAUfs[phasei];
208  }
209 
210  rAUf = mag(rAUf);
211 
212  // Update the fixedFluxPressure BCs to ensure flux consistency
213  {
214  surfaceScalarField::Boundary phib(phi.boundaryField());
215  phib = 0;
217  {
218  phaseModel& phase = phases[phasei];
219  phib +=
220  alphafs[phasei].boundaryField()*phase.phi()().boundaryField();
221  }
222 
224  (
225  p_rgh.boundaryFieldRef(),
226  (
227  phiHbyA.boundaryField() - phib
228  )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
229  );
230  }
231 
232  // Compressible pressure equations
233  PtrList<fvScalarMatrix> pEqnComps(phases.size());
234  PtrList<volScalarField> dmdts(fluid.dmdts());
236  {
237  phaseModel& phase = phases[phasei];
238  const volScalarField& alpha = phase;
239  volScalarField& rho = phase.thermoRef().rho();
240 
241  if (phase.compressible())
242  {
243  if (pimple.transonic())
244  {
246  (
247  IOobject::groupName("phid", phase.name()),
248  fvc::interpolate(phase.thermo().psi())*phase.phi()
249  );
250 
251  pEqnComps.set
252  (
253  phasei,
254  (
255  (
256  fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
257  - fvc::Sp
258  (
259  fvc::ddt(alpha) + fvc::div(phase.alphaPhi()),
260  rho
261  )
262  )/rho
263  + correction
264  (
265  (alpha/rho)*
266  (
267  phase.thermo().psi()*fvm::ddt(p_rgh)
268  + fvm::div(phid, p_rgh)
269  - fvm::Sp(fvc::div(phid), p_rgh)
270  )
271  )
272  ).ptr()
273  );
274 
276  (
277  pEqnComps[phasei].faceFluxCorrectionPtr()
278  );
279 
280  pEqnComps[phasei].relax();
281  }
282  else
283  {
284  pEqnComps.set
285  (
286  phasei,
287  (
288  (
289  fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
290  - fvc::Sp
291  (
292  (fvc::ddt(alpha) + fvc::div(phase.alphaPhi())),
293  rho
294  )
295  )/rho
296  + (alpha*phase.thermo().psi()/rho)
297  *correction(fvm::ddt(p_rgh))
298  ).ptr()
299  );
300  }
301  }
302 
303  if (fvOptions.appliesToField(rho.name()))
304  {
305  tmp<fvScalarMatrix> optEqn = fvOptions(alpha, rho);
306  if (pEqnComps.set(phasei))
307  {
308  pEqnComps[phasei] -= (optEqn&rho)/rho;
309  }
310  else
311  {
312  pEqnComps.set
313  (
314  phasei,
315  fvm::Su(- (optEqn&rho)/rho, p_rgh).ptr()
316  );
317  }
318  }
319 
320  if (dmdts.set(phasei))
321  {
322  if (pEqnComps.set(phasei))
323  {
324  pEqnComps[phasei] -= dmdts[phasei]/rho;
325  }
326  else
327  {
328  pEqnComps.set
329  (
330  phasei,
331  fvm::Su(- dmdts[phasei]/rho, p_rgh)
332  );
333  }
334  }
335  }
336 
337  // Cache p prior to solve for density update
338  volScalarField p_rgh_0(p_rgh);
339 
340  // Iterate over the pressure equation to correct for non-orthogonality
341  while (pimple.correctNonOrthogonal())
342  {
343  // Construct the transport part of the pressure equation
344  fvScalarMatrix pEqnIncomp
345  (
347  - fvm::laplacian(rAUf, p_rgh)
348  );
349 
350  {
351  fvScalarMatrix pEqn(pEqnIncomp);
352 
354  {
355  if (pEqnComps.set(phasei))
356  {
357  pEqn += pEqnComps[phasei];
358  }
359  }
360 
361  solve
362  (
363  pEqn,
364  mesh.solver(p_rgh.select(pimple.finalInnerIter()))
365  );
366  }
367 
368  // Correct fluxes and velocities on last non-orthogonal iteration
369  if (pimple.finalNonOrthogonalIter())
370  {
371  phi = phiHbyA + pEqnIncomp.flux();
372 
373  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
374 
375  forAll(fluid.movingPhases(), movingPhasei)
376  {
377  phaseModel& phase = fluid.movingPhases()[movingPhasei];
378 
379  phase.phiRef() =
380  phiHbyAs[phase.index()]
381  + alpharAUfs[phase.index()]*mSfGradp;
382 
383  // Set the phase dilatation rates
384  if (pEqnComps.set(phase.index()))
385  {
386  phase.divU(-pEqnComps[phase.index()] & p_rgh);
387  }
388  }
389 
390  if (partialElimination)
391  {
392  fluid.partialEliminationf(rAUfs);
393  }
394 
395  // Optionally relax pressure for velocity correction
396  p_rgh.relax();
397 
398  mSfGradp = pEqnIncomp.flux()/rAUf;
399 
400  forAll(fluid.movingPhases(), movingPhasei)
401  {
402  phaseModel& phase = fluid.movingPhases()[movingPhasei];
403 
404  phase.URef() = fvc::reconstruct(MRF.absolute(phase.phi()));
405  phase.URef().correctBoundaryConditions();
406  fvOptions.correct(phase.URef());
407  }
408  }
409  }
410 
411  // Update and limit the static pressure
412  p = max(p_rgh + rho*gh, pMin);
413 
414  // Limit p_rgh
415  p_rgh = p - rho*gh;
416 
417  // Update densities from change in p_rgh
419  {
420  phaseModel& phase = phases[phasei];
421  phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
422  }
423 
424  // Correct p_rgh for consistency with p and the updated densities
425  rho = fluid.rho();
426  p_rgh = p - rho*gh;
427  p_rgh.correctBoundaryConditions();
428 }
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix.
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
rho oldTime()
fv::options & fvOptions
rho
Definition: pEqn.H:1
pimpleNoLoopControl & pimple
zeroField Su
Definition: alphaSuSp.H:1
multiphaseSystem & fluid
Definition: createFields.H:11
p
Definition: pEqn.H:50
label phasei
Definition: pEqn.H:27
IOMRFZoneList & MRF
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
phiHbyA
Definition: pEqn.H:20
multiphaseSystem::phaseModelList & phases
Definition: createFields.H:12
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
engineTime & runTime
PtrList< surfaceScalarField > AFfs(fluid.AFfs())
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:3
surfaceScalarField phid("phid", fvc::interpolate(psi) *(fvc::flux(HbyA)+MRF.zeroFilter(rhorAUf *fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho))))
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
dynamicFvMesh & mesh
PtrList< surfaceScalarField > alpharAUfs(phases.size())
rhoEqn solve()
PtrList< surfaceScalarField > rAUfs
Definition: pEqn.H:30
p_rgh
Definition: pEqn.H:152
const dimensionSet dimForce
const surfaceScalarField & ghf
dimensionedScalar pMin("pMin", dimPressure, fluid)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
tmp< volScalarField > byDt(const volScalarField &vf)
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.
const dimensionSet dimDensity
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)
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho) *mesh.magSf())
PtrList< fvVectorMatrix > UEqns(fluid.phases().size())
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
PtrList< surfaceScalarField > alphaRho0fs(phases.size())
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryFieldRef(),(phiHbyA.boundaryField() - MRF.relative(phib))/(mesh.magSf().boundaryField() *rAUf.boundaryField()))
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
volScalarField p_rgh_0(p_rgh)
PtrList< surfaceScalarField > phiHbyAs(fluid.phases().size())
const volScalarField & gh
dimensioned< scalar > mag(const dimensioned< Type > &)
phi
Definition: pEqn.H:18
phib
Definition: pEqn.H:194
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 dimensionedVector & g
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
void deleteDemandDrivenData(DataPtr &dataPtr)
PtrList< surfaceScalarField > phiFfs(fluid.phiFfs(rAUfs))
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
zeroField Sp
Definition: alphaSuSp.H:2
const dimensionSet dimVelocity