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 // Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes
35 tmp<surfaceScalarField> phiF1;
36 tmp<surfaceScalarField> phiF2;
37 {
38  // Turbulent-dispersion diffusivity
39  volScalarField D(fluid.D());
40 
41  // Phase-1 turbulent dispersion and particle-pressure flux
43  (
45  (
46  rAU1*(D + phase1.turbulence().pPrime())
47  )
48  );
49 
50  // Phase-2 turbulent dispersion and particle-pressure flux
52  (
54  (
55  rAU2*(D + phase2.turbulence().pPrime())
56  )
57  );
58 
59  // Cache the net diffusivity for implicit diffusion treatment in the
60  // phase-fraction equation
61  if (implicitPhasePressure)
62  {
63  fluid.pPrimeByA() = Df1 + Df2;
64  }
65 
66  // Lift and wall-lubrication forces
67  volVectorField F(fluid.F());
68 
69  // Phase-fraction face-gradient
71 
72  // Phase-1 dispersion, lift and wall-lubrication flux
73  phiF1 =
74  (
76  + (fvc::interpolate(rAU1*F) & mesh.Sf())
77  );
78 
79  // Phase-1 dispersion, lift and wall-lubrication flux
80  phiF2 =
81  (
83  - (fvc::interpolate(rAU2*F) & mesh.Sf())
84  );
85 }
86 
87 
88 // --- Pressure corrector loop
89 while (pimple.correct())
90 {
91  // Update continuity errors due to temperature changes
92  #include "correctContErrs.H"
93 
94  // Correct fixed-flux BCs to be consistent with the velocity BCs
95  MRF.correctBoundaryFlux(U1, phi1);
96  MRF.correctBoundaryFlux(U2, phi2);
97 
98  volVectorField HbyA1
99  (
100  IOobject::groupName("HbyA", phase1.name()),
101  U1
102  );
103  HbyA1 =
104  rAU1
105  *(
106  U1Eqn.H()
107  + max(phase1.residualAlpha() - alpha1, scalar(0))
108  *rho1*U1.oldTime()/runTime.deltaT()
109  );
110 
111  volVectorField HbyA2
112  (
113  IOobject::groupName("HbyA", phase2.name()),
114  U2
115  );
116  HbyA2 =
117  rAU2
118  *(
119  U2Eqn.H()
120  + max(phase2.residualAlpha() - alpha2, scalar(0))
121  *rho2*U2.oldTime()/runTime.deltaT()
122  );
123 
124  // Mean density for buoyancy force and p_rgh -> p
125  volScalarField rho("rho", fluid.rho());
126 
128  (
129  "ghSnGradRho",
130  ghf*fvc::snGrad(rho)*mesh.magSf()
131  );
132 
133  surfaceScalarField phig1
134  (
135  alpharAUf1
136  *(
138  - alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
139  )
140  );
141 
142  surfaceScalarField phig2
143  (
144  alpharAUf2
145  *(
147  - alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
148  )
149  );
150 
151 
152  // ddtPhiCorr filter -- only apply in pure(ish) phases
154  surfaceScalarField phiCorrCoeff1(pos(alphaf1Bar - 0.99));
155  surfaceScalarField phiCorrCoeff2(pos(0.01 - alphaf1Bar));
156 
157  forAll(mesh.boundary(), patchi)
158  {
159  // Set ddtPhiCorr to 0 on non-coupled boundaries
160  if
161  (
162  !mesh.boundary()[patchi].coupled()
163  || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
164  )
165  {
166  phiCorrCoeff1.boundaryField()[patchi] = 0;
167  phiCorrCoeff2.boundaryField()[patchi] = 0;
168  }
169  }
170 
171  // Phase-1 predicted flux
172  surfaceScalarField phiHbyA1
173  (
174  IOobject::groupName("phiHbyA", phase1.name()),
175  (fvc::interpolate(HbyA1) & mesh.Sf())
176  + phiCorrCoeff1*fvc::interpolate(alpha1.oldTime()*rho1.oldTime()*rAU1)
177  *(
178  MRF.absolute(phi1.oldTime())
179  - (fvc::interpolate(U1.oldTime()) & mesh.Sf())
180  )/runTime.deltaT()
181  - phiF1()
182  - phig1
183  );
184 
185  // Phase-2 predicted flux
186  surfaceScalarField phiHbyA2
187  (
188  IOobject::groupName("phiHbyA", phase2.name()),
189  (fvc::interpolate(HbyA2) & mesh.Sf())
190  + phiCorrCoeff2*fvc::interpolate(alpha2.oldTime()*rho2.oldTime()*rAU2)
191  *(
192  MRF.absolute(phi2.oldTime())
193  - (fvc::interpolate(U2.oldTime()) & mesh.Sf())
194  )/runTime.deltaT()
195  - phiF2()
196  - phig2
197  );
198 
199  // Face-drag coefficients
202 
203  // Construct the mean predicted flux
204  // including explicit drag contributions based on absolute fluxes
206  (
207  "phiHbyA",
208  alphaf1*(phiHbyA1 + rAUKd1*MRF.absolute(phi2))
209  + alphaf2*(phiHbyA2 + rAUKd2*MRF.absolute(phi1))
210  );
211  MRF.makeRelative(phiHbyA);
212 
213  // Construct pressure "diffusivity"
215  (
216  "rAUf",
218  );
219 
220  // Update the fixedFluxPressure BCs to ensure flux consistency
222  (
223  p_rgh.boundaryField(),
224  (
225  phiHbyA.boundaryField()
226  - (
227  alphaf1.boundaryField()*phi1.boundaryField()
228  + alphaf2.boundaryField()*phi2.boundaryField()
229  )
230  )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
231  );
232 
233  tmp<fvScalarMatrix> pEqnComp1;
234  tmp<fvScalarMatrix> pEqnComp2;
235 
236  // Construct the compressibility parts of the pressure equation
237  if (pimple.transonic())
238  {
239  surfaceScalarField phid1
240  (
241  IOobject::groupName("phid", phase1.name()),
243  );
244  surfaceScalarField phid2
245  (
246  IOobject::groupName("phid", phase2.name()),
248  );
249 
250  pEqnComp1 =
251  (
252  contErr1
254  )/rho1
256  (
258  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
259  );
260  deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
261  pEqnComp1().relax();
262 
263  pEqnComp2 =
264  (
265  contErr2
267  )/rho2
269  (
271  + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
272  );
273  deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
274  pEqnComp2().relax();
275  }
276  else
277  {
278  pEqnComp1 =
279  (
280  contErr1
282  )/rho1
284 
285  pEqnComp2 =
286  (
287  contErr2
289  )/rho2
291  }
292 
293  // Cache p prior to solve for density update
295 
296  // Iterate over the pressure equation to correct for non-orthogonality
297  while (pimple.correctNonOrthogonal())
298  {
299  // Construct the transport part of the pressure equation
300  fvScalarMatrix pEqnIncomp
301  (
302  fvc::div(phiHbyA)
303  - fvm::laplacian(rAUf, p_rgh)
304  );
305 
306  solve
307  (
308  pEqnComp1() + pEqnComp2() + pEqnIncomp,
309  mesh.solver(p_rgh.select(pimple.finalInnerIter()))
310  );
311 
312  // Correct fluxes and velocities on last non-orthogonal iteration
313  if (pimple.finalNonOrthogonalIter())
314  {
315  phi = phiHbyA + pEqnIncomp.flux();
316 
317  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
318 
319  // Partial-elimination phase-flux corrector
320  {
321  surfaceScalarField phi1s
322  (
323  phiHbyA1 + alpharAUf1*mSfGradp
324  );
325 
326  surfaceScalarField phi2s
327  (
328  phiHbyA2 + alpharAUf2*mSfGradp
329  );
330 
332  (
333  ((phi1s + rAUKd1*phi2s) - (phi2s + rAUKd2*phi1s))
334  /(1 - rAUKd1*rAUKd2)
335  );
336 
337  phi1 = phi + alphaf2*phir;
338  phi2 = phi - alphaf1*phir;
339  }
340 
341  // Compressibility correction for phase-fraction equations
342  fluid.dgdt() =
343  (
344  alpha1*(pEqnComp2 & p_rgh)
345  - alpha2*(pEqnComp1 & p_rgh)
346  );
347 
348  // Optionally relax pressure for velocity correction
349  p_rgh.relax();
350 
351  mSfGradp = pEqnIncomp.flux()/rAUf;
352 
353  // Partial-elimination phase-velocity corrector
354  {
355  volVectorField Us1
356  (
357  HbyA1
358  + fvc::reconstruct(alpharAUf1*mSfGradp - phiF1() - phig1)
359  );
360 
361  volVectorField Us2
362  (
363  HbyA2
364  + fvc::reconstruct(alpharAUf2*mSfGradp - phiF2() - phig2)
365  );
366 
367  volScalarField D1(rAU1*Kd);
368  volScalarField D2(rAU2*Kd);
369 
370  U = alpha1*(Us1 + D1*U2) + alpha2*(Us2 + D2*U1);
371  volVectorField Ur(((1 - D2)*Us1 - (1 - D1)*Us2)/(1 - D1*D2));
372 
373  U1 = U + alpha2*Ur;
374  U1.correctBoundaryConditions();
375  fvOptions.correct(U1);
376 
377  U2 = U - alpha1*Ur;
378  U2.correctBoundaryConditions();
379  fvOptions.correct(U2);
380 
381  U = fluid.U();
382  }
383  }
384  }
385 
386  // Update and limit the static pressure
387  p = max(p_rgh + rho*gh, pMin);
388 
389  // Limit p_rgh
390  p_rgh = p - rho*gh;
391 
392  // Update densities from change in p_rgh
393  rho1 += psi1*(p_rgh - p_rgh_0);
394  rho2 += psi2*(p_rgh - p_rgh_0);
395 
396  // Correct p_rgh for consistency with p and the updated densities
397  rho = fluid.rho();
398  p_rgh = p - rho*gh;
399  p_rgh.correctBoundaryConditions();
400 }
401 
402 // Update the phase kinetic energies
403 K1 = 0.5*magSqr(U1);
404 K2 = 0.5*magSqr(U2);
405 
406 // Update the pressure time-derivative if required
407 if (thermo1.dpdt() || thermo2.dpdt())
408 {
409  dpdt = fvc::ddt(p);
410 }
phaseModel & phase1
Definition: createFields.H:12
volScalarField & dpdt
phaseModel & phase2
Definition: createFields.H:13
K1
Definition: pEqn.H:403
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
contErr1
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
dimensioned< scalar > magSqr(const dimensioned< Type > &)
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
contErr2
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() ))
rhoThermo & thermo1
Definition: createFields.H:39
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
surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1))
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
K2
Definition: pEqn.H:404
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
surfaceScalarField Df1(fvc::interpolate(D+phase1.turbulence().pPrime()))
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
rhoThermo & thermo2
Definition: createFields.H:40
dimensionedScalar pMin("pMin", dimPressure, fluid)
surfaceScalarField Df2(fvc::interpolate(D+phase2.turbulence().pPrime()))
volVectorField F(fluid.F())
U
Definition: pEqn.H:82
const surfaceScalarField & ghf
const volScalarField & gh