pEqn.H
Go to the documentation of this file.
2 surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1);
3 
5 (
6  IOobject::groupName("rAU", phase1.name()),
7  1.0
8  /(
9  U1Eqn.A()
10  + max(phase1.residualAlpha() - alpha1, scalar(0))
11  *rho1/runTime.deltaT()
12  )
13 );
15 (
16  IOobject::groupName("rAU", phase2.name()),
17  1.0
18  /(
19  U2Eqn.A()
20  + max(phase2.residualAlpha() - alpha2, scalar(0))
21  *rho2/runTime.deltaT()
22  )
23 );
24 
26 (
27  fvc::interpolate(max(alpha1, phase1.residualAlpha())*rAU1)
28 );
30 (
31  fvc::interpolate(max(alpha2, phase2.residualAlpha())*rAU2)
32 );
33 
34 volScalarField Kd(fluid.Kd());
35 
36 // Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes
37 tmp<surfaceScalarField> phiF1;
38 tmp<surfaceScalarField> phiF2;
39 {
40  // Turbulent-dispersion diffusivity
41  volScalarField D(fluid.D());
42 
43  // Phase-1 turbulent dispersion and particle-pressure flux
44  tmp<surfaceScalarField> DbyA1
45  (
47  (
48  rAU1*(D + phase1.turbulence().pPrime())
49  )
50  );
51 
52  // Phase-2 turbulent dispersion and particle-pressure flux
53  tmp<surfaceScalarField> DbyA2
54  (
56  (
57  rAU2*(D + phase2.turbulence().pPrime())
58  )
59  );
60 
61  // Lift and wall-lubrication forces
62  volVectorField F(fluid.F());
63 
64  // Phase-fraction face-gradient
66 
67  // Phase-1 dispersion, lift and wall-lubrication flux
68  phiF1 =
69  (
71  + (fvc::interpolate(rAU1*F) & mesh.Sf())
72  );
73 
74  // Phase-2 dispersion, lift and wall-lubrication flux
75  phiF2 =
76  (
78  - (fvc::interpolate(rAU2*F) & mesh.Sf())
79  );
80 
81  // Cache the phase diffusivities for implicit treatment in the
82  // phase-fraction equation
83  if (implicitPhasePressure)
84  {
85  phase1.DbyA(DbyA1);
86  phase2.DbyA(DbyA2);
87  }
88 }
89 
90 
91 // --- Pressure corrector loop
92 while (pimple.correct())
93 {
94  // Update continuity errors due to temperature changes
95  fluid.correct();
96 
97  volScalarField rho("rho", fluid.rho());
98 
99  // Correct p_rgh for consistency with p and the updated densities
100  p_rgh = p - rho*gh;
101 
102  // Correct fixed-flux BCs to be consistent with the velocity BCs
103  MRF.correctBoundaryFlux(U1, phi1);
104  MRF.correctBoundaryFlux(U2, phi2);
105 
106  volVectorField HbyA1
107  (
108  IOobject::groupName("HbyA", phase1.name()),
109  U1
110  );
111  HbyA1 =
112  rAU1
113  *(
114  U1Eqn.H()
115  + max(phase1.residualAlpha() - alpha1, scalar(0))
116  *rho1*U1.oldTime()/runTime.deltaT()
117  );
118 
119  volVectorField HbyA2
120  (
121  IOobject::groupName("HbyA", phase2.name()),
122  U2
123  );
124  HbyA2 =
125  rAU2
126  *(
127  U2Eqn.H()
128  + max(phase2.residualAlpha() - alpha2, scalar(0))
129  *rho2*U2.oldTime()/runTime.deltaT()
130  );
131 
133  (
134  "ghSnGradRho",
135  ghf*fvc::snGrad(rho)*mesh.magSf()
136  );
137 
138  surfaceScalarField phig1
139  (
140  alpharAUf1
141  *(
143  - alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
144  )
145  );
146 
147  surfaceScalarField phig2
148  (
149  alpharAUf2
150  *(
152  - alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
153  )
154  );
155 
156 
157  // ddtPhiCorr filter -- only apply in pure(ish) phases
159  surfaceScalarField phiCorrCoeff1(pos(alphaf1Bar - 0.99));
160  surfaceScalarField phiCorrCoeff2(pos(0.01 - alphaf1Bar));
161 
162  forAll(mesh.boundary(), patchi)
163  {
164  // Set ddtPhiCorr to 0 on non-coupled boundaries
165  if
166  (
167  !mesh.boundary()[patchi].coupled()
168  || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
169  )
170  {
171  phiCorrCoeff1.boundaryField()[patchi] = 0;
172  phiCorrCoeff2.boundaryField()[patchi] = 0;
173  }
174  }
175 
176  // Phase-1 predicted flux
177  surfaceScalarField phiHbyA1
178  (
179  IOobject::groupName("phiHbyA", phase1.name()),
180  (fvc::interpolate(HbyA1) & mesh.Sf())
181  + phiCorrCoeff1*fvc::interpolate(alpha1.oldTime()*rho1.oldTime()*rAU1)
182  *(
183  MRF.absolute(phi1.oldTime())
184  - (fvc::interpolate(U1.oldTime()) & mesh.Sf())
185  )/runTime.deltaT()
186  - phiF1()
187  - phig1
188  );
189 
190  // Phase-2 predicted flux
191  surfaceScalarField phiHbyA2
192  (
193  IOobject::groupName("phiHbyA", phase2.name()),
194  (fvc::interpolate(HbyA2) & mesh.Sf())
195  + phiCorrCoeff2*fvc::interpolate(alpha2.oldTime()*rho2.oldTime()*rAU2)
196  *(
197  MRF.absolute(phi2.oldTime())
198  - (fvc::interpolate(U2.oldTime()) & mesh.Sf())
199  )/runTime.deltaT()
200  - phiF2()
201  - phig2
202  );
203 
204  // Face-drag coefficients
207 
208  // Construct the mean predicted flux
209  // including explicit drag contributions based on absolute fluxes
211  (
212  "phiHbyA",
213  alphaf1*(phiHbyA1 + rAUKd1*MRF.absolute(phi2))
214  + alphaf2*(phiHbyA2 + rAUKd2*MRF.absolute(phi1))
215  );
216  MRF.makeRelative(phiHbyA);
217 
218  // Construct pressure "diffusivity"
220  (
221  "rAUf",
223  );
224 
225  // Update the fixedFluxPressure BCs to ensure flux consistency
227  (
228  p_rgh.boundaryField(),
229  (
230  phiHbyA.boundaryField()
231  - (
232  alphaf1.boundaryField()*phi1.boundaryField()
233  + alphaf2.boundaryField()*phi2.boundaryField()
234  )
235  )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
236  );
237 
238  tmp<fvScalarMatrix> pEqnComp1;
239  tmp<fvScalarMatrix> pEqnComp2;
240 
241  // Construct the compressibility parts of the pressure equation
242  if (pimple.transonic())
243  {
244  if (phase1.compressible())
245  {
246  surfaceScalarField phid1
247  (
248  IOobject::groupName("phid", phase1.name()),
250  );
251 
252  pEqnComp1 =
253  (
254  phase1.continuityError()
256  )/rho1
258  (
259  psi1*fvm::ddt(p_rgh)
260  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
261  );
262 
263  deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
264  pEqnComp1().relax();
265  }
266 
267  if (phase2.compressible())
268  {
269  surfaceScalarField phid2
270  (
271  IOobject::groupName("phid", phase2.name()),
273  );
274 
275  pEqnComp2 =
276  (
277  phase2.continuityError()
279  )/rho2
281  (
282  psi2*fvm::ddt(p_rgh)
283  + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
284  );
285  deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
286  pEqnComp2().relax();
287  }
288  }
289  else
290  {
291  if (phase1.compressible())
292  {
293  pEqnComp1 =
294  (
295  phase1.continuityError()
297  )/rho1
298  + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
299  }
300 
301  if (phase2.compressible())
302  {
303  pEqnComp2 =
304  (
305  phase2.continuityError()
307  )/rho2
308  + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
309  }
310  }
311 
312  if (fluid.transfersMass())
313  {
314  if (pEqnComp1.valid())
315  {
316  pEqnComp1() -= fluid.dmdt()/rho1;
317  }
318  else
319  {
320  pEqnComp1 = fvm::Su(-fluid.dmdt()/rho1, p_rgh);
321  }
322 
323  if (pEqnComp2.valid())
324  {
325  pEqnComp2() += fluid.dmdt()/rho2;
326  }
327  else
328  {
329  pEqnComp2 = fvm::Su(fluid.dmdt()/rho2, p_rgh);
330  }
331  }
332 
333  // Cache p prior to solve for density update
334  volScalarField p_rgh_0(p_rgh);
335 
336  // Iterate over the pressure equation to correct for non-orthogonality
337  while (pimple.correctNonOrthogonal())
338  {
339  // Construct the transport part of the pressure equation
340  fvScalarMatrix pEqnIncomp
341  (
342  fvc::div(phiHbyA)
343  - fvm::laplacian(rAUf, p_rgh)
344  );
345 
346  {
347  fvScalarMatrix pEqn(pEqnIncomp);
348 
349  if (pEqnComp1.valid())
350  {
351  pEqn += pEqnComp1();
352  }
353 
354  if (pEqnComp2.valid())
355  {
356  pEqn += pEqnComp2();
357  }
358 
359  solve
360  (
361  pEqn,
362  mesh.solver(p_rgh.select(pimple.finalInnerIter()))
363  );
364  }
365 
366  // Correct fluxes and velocities on last non-orthogonal iteration
367  if (pimple.finalNonOrthogonalIter())
368  {
369  phi = phiHbyA + pEqnIncomp.flux();
370 
371  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
372 
373  // Partial-elimination phase-flux corrector
374  {
375  surfaceScalarField phi1s
376  (
377  phiHbyA1 + alpharAUf1*mSfGradp
378  );
379 
380  surfaceScalarField phi2s
381  (
382  phiHbyA2 + alpharAUf2*mSfGradp
383  );
384 
386  (
387  ((phi1s + rAUKd1*phi2s) - (phi2s + rAUKd2*phi1s))
388  /(1 - rAUKd1*rAUKd2)
389  );
390 
391  phi1 = phi + alphaf2*phir;
392  phi2 = phi - alphaf1*phir;
393  }
394 
395  // Set the phase dilatation rates
396  if (phase1.compressible())
397  {
398  phase1.divU(-pEqnComp1 & p_rgh);
399  }
400  if (phase2.compressible())
401  {
402  phase2.divU(-pEqnComp2 & p_rgh);
403  }
404 
405  // Optionally relax pressure for velocity correction
406  p_rgh.relax();
407 
408  mSfGradp = pEqnIncomp.flux()/rAUf;
409 
410  // Partial-elimination phase-velocity corrector
411  {
412  volVectorField Us1
413  (
414  HbyA1
415  + fvc::reconstruct(alpharAUf1*mSfGradp - phiF1() - phig1)
416  );
417 
418  volVectorField Us2
419  (
420  HbyA2
421  + fvc::reconstruct(alpharAUf2*mSfGradp - phiF2() - phig2)
422  );
423 
424  volScalarField D1(rAU1*Kd);
425  volScalarField D2(rAU2*Kd);
426 
427  volVectorField U(alpha1*(Us1 + D1*U2) + alpha2*(Us2 + D2*U1));
428  volVectorField Ur(((1 - D2)*Us1 - (1 - D1)*Us2)/(1 - D1*D2));
429 
430  U1 = U + alpha2*Ur;
431  U1.correctBoundaryConditions();
432  fvOptions.correct(U1);
433 
434  U2 = U - alpha1*Ur;
435  U2.correctBoundaryConditions();
436  fvOptions.correct(U2);
437  }
438  }
439  }
440 
441  // Update and limit the static pressure
442  p = max(p_rgh + rho*gh, pMin);
443 
444  // Limit p_rgh
445  p_rgh = p - rho*gh;
446 
447  // Update densities from change in p_rgh
448  rho1 += psi1*(p_rgh - p_rgh_0);
449  rho2 += psi2*(p_rgh - p_rgh_0);
450 
451  // Correct p_rgh for consistency with p and the updated densities
452  rho = fluid.rho();
453  p_rgh = p - rho*gh;
454  p_rgh.correctBoundaryConditions();
455 }
phaseModel & phase1
Definition: createFields.H:12
tmp< surfaceScalarField > DbyA1(fvc::interpolate( rAU1 *(D+phase1.turbulence().pPrime()) ))
phaseModel & phase2
Definition: createFields.H:13
tmp< surfaceScalarField > phiF1
Definition: pEqn.H:37
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix.
volScalarField Kd(fluid.Kd())
volVectorField & U1
Definition: createFields.H:18
multiphaseSystem & fluid
Definition: createFields.H:10
phi
Definition: pEqn.H:20
dimensioned< scalar > mag(const dimensioned< Type > &)
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryField(),( phiHbyA.boundaryField() -MRF.relative(mesh.Sf().boundaryField()&U.boundaryField()) *rho.boundaryField() )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField()))
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
surfaceScalarField & alphaPhi1
Definition: createFields.H:20
volVectorField & U2
Definition: createFields.H:23
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
surfaceScalarField phir(phic *interface.nHatf())
surfaceScalarField alphaf2("alphaf2", scalar(1)-alphaf1)
void deleteDemandDrivenData(DataPtr &dataPtr)
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
psi2
Definition: TEqns.H:35
volScalarField & p_rgh
tmp< surfaceScalarField > interpolate(const RhoType &rho)
surfaceScalarField & phi2
Definition: createFields.H:24
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
dynamicFvMesh & mesh
const dictionary & pimple
p
Definition: pEqn.H:59
volScalarField rAU1(IOobject::groupName("rAU", phase1.name()), 1.0/( U1Eqn.A() +max(phase1.residualAlpha()-alpha1, scalar(0)) *rho1/runTime.deltaT() ))
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:44
rhoEqn solve()
fv::IOoptionList & fvOptions
surfaceScalarField & alphaPhi2
Definition: createFields.H:25
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf())
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho)*mesh.magSf())
IOMRFZoneList & MRF
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
forAll(phases, phasei)
Definition: pEqn.H:5
tmp< surfaceScalarField > DbyA2(fvc::interpolate( rAU2 *(D+phase2.turbulence().pPrime()) ))
surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1))
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
rho1
Definition: pEqn.H:124
label patchi
volScalarField p_rgh_0(p_rgh)
phiHbyA
Definition: pEqn.H:21
rho2
Definition: pEqn.H:125
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
surfaceScalarField & phi1
Definition: createFields.H:19
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
alpha2
Definition: alphaEqn.H:112
const dimensionedVector & g
tmp< surfaceScalarField > phiF2
Definition: pEqn.H:38
tmp< GeometricField< Type, fvPatchField, volMesh > > average(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Area-weighted average a surfaceField creating a volField.
Definition: fvcAverage.C:46
surfaceScalarField alpharAUf2(fvc::interpolate(max(alpha2, phase2.residualAlpha())*rAU2))
psi1
Definition: TEqns.H:34
surfaceScalarField alpharAUf1(fvc::interpolate(max(alpha1, phase1.residualAlpha())*rAU1))
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
dimensionedScalar pos(const dimensionedScalar &ds)
volScalarField rAU2(IOobject::groupName("rAU", phase2.name()), 1.0/( U2Eqn.A() +max(phase2.residualAlpha()-alpha2, scalar(0)) *rho2/runTime.deltaT() ))
rho
Definition: pEqn.H:1
dimensionedScalar pMin("pMin", dimPressure, fluid)
volVectorField F(fluid.F())
U
Definition: pEqn.H:82
const surfaceScalarField & ghf
const volScalarField & gh