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 = Df1*snGradAlpha1 + fvc::flux(rAU1*F);
74 
75  // Phase-1 dispersion, lift and wall-lubrication flux
76  phiF2 = - Df2*snGradAlpha1 - fvc::flux(rAU2*F);
77 }
78 
79 
80 // --- Pressure corrector loop
81 while (pimple.correct())
82 {
83  // Update continuity errors due to temperature changes
84  #include "correctContErrs.H"
85 
86  volScalarField rho("rho", fluid.rho());
87 
88  // Correct p_rgh for consistency with p and the updated densities
89  p_rgh = p - rho*gh;
90 
91  // Correct fixed-flux BCs to be consistent with the velocity BCs
92  MRF.correctBoundaryFlux(U1, phi1);
93  MRF.correctBoundaryFlux(U2, phi2);
94 
95  volVectorField HbyA1
96  (
97  IOobject::groupName("HbyA", phase1.name()),
98  U1
99  );
100  HbyA1 =
101  rAU1
102  *(
103  U1Eqn.H()
104  + max(phase1.residualAlpha() - alpha1, scalar(0))
105  *rho1*U1.oldTime()/runTime.deltaT()
106  );
107 
108  volVectorField HbyA2
109  (
110  IOobject::groupName("HbyA", phase2.name()),
111  U2
112  );
113  HbyA2 =
114  rAU2
115  *(
116  U2Eqn.H()
117  + max(phase2.residualAlpha() - alpha2, scalar(0))
118  *rho2*U2.oldTime()/runTime.deltaT()
119  );
120 
122  (
123  "ghSnGradRho",
124  ghf*fvc::snGrad(rho)*mesh.magSf()
125  );
126 
127  surfaceScalarField phig1
128  (
129  alpharAUf1
130  *(
132  - alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
133  )
134  );
135 
136  surfaceScalarField phig2
137  (
138  alpharAUf2
139  *(
141  - alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
142  )
143  );
144 
145 
146  // ddtPhiCorr filter -- only apply in pure(ish) phases
148  surfaceScalarField phiCorrCoeff1(pos0(alphaf1Bar - 0.99));
149  surfaceScalarField phiCorrCoeff2(pos0(0.01 - alphaf1Bar));
150 
151  {
152  surfaceScalarField::Boundary& phiCorrCoeff1Bf =
153  phiCorrCoeff1.boundaryFieldRef();
154 
155  surfaceScalarField::Boundary& phiCorrCoeff2Bf =
156  phiCorrCoeff2.boundaryFieldRef();
157 
158  forAll(mesh.boundary(), patchi)
159  {
160  // Set ddtPhiCorr to 0 on non-coupled boundaries
161  if
162  (
163  !mesh.boundary()[patchi].coupled()
164  || isA<cyclicAMIFvPatch>(mesh.boundary()[patchi])
165  )
166  {
167  phiCorrCoeff1Bf[patchi] = 0;
168  phiCorrCoeff2Bf[patchi] = 0;
169  }
170  }
171  }
172 
173  // Phase-1 predicted flux
174  surfaceScalarField phiHbyA1
175  (
176  IOobject::groupName("phiHbyA", phase1.name()),
177  fvc::flux(HbyA1)
178  + phiCorrCoeff1*fvc::interpolate(alpha1.oldTime()*rho1.oldTime()*rAU1)
179  *(
180  MRF.absolute(phi1.oldTime())
181  - fvc::flux(U1.oldTime())
182  )/runTime.deltaT()
183  - phiF1()
184  - phig1
185  );
186 
187  // Phase-2 predicted flux
188  surfaceScalarField phiHbyA2
189  (
190  IOobject::groupName("phiHbyA", phase2.name()),
191  fvc::flux(HbyA2)
192  + phiCorrCoeff2*fvc::interpolate(alpha2.oldTime()*rho2.oldTime()*rAU2)
193  *(
194  MRF.absolute(phi2.oldTime())
195  - fvc::flux(U2.oldTime())
196  )/runTime.deltaT()
197  - phiF2()
198  - phig2
199  );
200 
201  // Face-drag coefficients
204 
205  // Construct the mean predicted flux
206  // including explicit drag contributions based on absolute fluxes
208  (
209  "phiHbyA",
210  alphaf1*(phiHbyA1 + rAUKd1*MRF.absolute(phi2))
211  + alphaf2*(phiHbyA2 + rAUKd2*MRF.absolute(phi1))
212  );
213  MRF.makeRelative(phiHbyA);
214 
215  // Construct pressure "diffusivity"
217  (
218  "rAUf",
220  );
221 
222  // Update the fixedFluxPressure BCs to ensure flux consistency
224  (
225  p_rgh.boundaryFieldRef(),
226  (
227  phiHbyA.boundaryField()
228  - (
229  alphaf1.boundaryField()*phi1.boundaryField()
230  + alphaf2.boundaryField()*phi2.boundaryField()
231  )
232  )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
233  );
234 
235  tmp<fvScalarMatrix> pEqnComp1;
236  tmp<fvScalarMatrix> pEqnComp2;
237 
238  // Construct the compressibility parts of the pressure equation
239  if (pimple.transonic())
240  {
241  surfaceScalarField phid1
242  (
243  IOobject::groupName("phid", phase1.name()),
245  );
246  surfaceScalarField phid2
247  (
248  IOobject::groupName("phid", phase2.name()),
250  );
251 
252  pEqnComp1 =
253  (
254  contErr1
256  )/rho1
257  + correction
258  (
259  (alpha1/rho1)*
260  (
261  psi1*fvm::ddt(p_rgh)
262  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
263  )
264  );
265  pEqnComp1.ref().relax();
266 
267  pEqnComp2 =
268  (
269  contErr2
271  )/rho2
272  + correction
273  (
274  (alpha2/rho2)*
275  (
276  psi2*fvm::ddt(p_rgh)
277  + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
278  )
279  );
280  pEqnComp2.ref().relax();
281  }
282  else
283  {
284  pEqnComp1 =
285  (
286  contErr1
288  )/rho1
289  + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
290 
291  pEqnComp2 =
292  (
293  contErr2
295  )/rho2
296  + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
297  }
298 
299  // Cache p prior to solve for density update
300  volScalarField p_rgh_0(p_rgh);
301 
302  // Iterate over the pressure equation to correct for non-orthogonality
303  while (pimple.correctNonOrthogonal())
304  {
305  // Construct the transport part of the pressure equation
306  fvScalarMatrix pEqnIncomp
307  (
308  fvc::div(phiHbyA)
309  - fvm::laplacian(rAUf, p_rgh)
310  );
311 
312  solve
313  (
314  pEqnComp1() + pEqnComp2() + pEqnIncomp,
315  mesh.solver(p_rgh.select(pimple.finalInnerIter()))
316  );
317 
318  // Correct fluxes and velocities on last non-orthogonal iteration
319  if (pimple.finalNonOrthogonalIter())
320  {
321  phi = phiHbyA + pEqnIncomp.flux();
322 
323  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
324 
325  // Partial-elimination phase-flux corrector
326  {
327  surfaceScalarField phi1s
328  (
329  phiHbyA1 + alpharAUf1*mSfGradp
330  );
331 
332  surfaceScalarField phi2s
333  (
334  phiHbyA2 + alpharAUf2*mSfGradp
335  );
336 
338  (
339  ((phi1s + rAUKd1*phi2s) - (phi2s + rAUKd2*phi1s))
340  /(1 - rAUKd1*rAUKd2)
341  );
342 
343  phi1 = phi + alphaf2*phir;
344  phi2 = phi - alphaf1*phir;
345  }
346 
347  // Compressibility correction for phase-fraction equations
348  fluid.dgdt() =
349  (
350  alpha1*(pEqnComp2 & p_rgh)
351  - alpha2*(pEqnComp1 & p_rgh)
352  );
353 
354  // Optionally relax pressure for velocity correction
355  p_rgh.relax();
356 
357  mSfGradp = pEqnIncomp.flux()/rAUf;
358 
359  // Partial-elimination phase-velocity corrector
360  {
361  volVectorField Us1
362  (
363  HbyA1
364  + fvc::reconstruct(alpharAUf1*mSfGradp - phiF1() - phig1)
365  );
366 
367  volVectorField Us2
368  (
369  HbyA2
370  + fvc::reconstruct(alpharAUf2*mSfGradp - phiF2() - phig2)
371  );
372 
373  volScalarField D1(rAU1*Kd);
374  volScalarField D2(rAU2*Kd);
375 
376  U = alpha1*(Us1 + D1*U2) + alpha2*(Us2 + D2*U1);
377  volVectorField Ur(((1 - D2)*Us1 - (1 - D1)*Us2)/(1 - D1*D2));
378 
379  U1 = U + alpha2*Ur;
380  U1.correctBoundaryConditions();
381  fvOptions.correct(U1);
382 
383  U2 = U - alpha1*Ur;
384  U2.correctBoundaryConditions();
385  fvOptions.correct(U2);
386 
387  U = fluid.U();
388  }
389  }
390  }
391 
392  // Update and limit the static pressure
393  p = max(p_rgh + rho*gh, pMin);
394 
395  // Limit p_rgh
396  p_rgh = p - rho*gh;
397 
398  // Update densities from change in p_rgh
399  rho1 += psi1*(p_rgh - p_rgh_0);
400  rho2 += psi2*(p_rgh - p_rgh_0);
401 
402  // Correct p_rgh for consistency with p and the updated densities
403  rho = fluid.rho();
404  p_rgh = p - rho*gh;
405  p_rgh.correctBoundaryConditions();
406 }
407 
408 // Update the phase kinetic energies
409 K1 = 0.5*magSqr(U1);
410 K2 = 0.5*magSqr(U2);
411 
412 // Update the pressure time-derivative if required
413 if (thermo1.dpdt() || thermo2.dpdt())
414 {
415  dpdt = fvc::ddt(p);
416 }
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix.
contErr1
psi2
Definition: TEqns.H:35
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
const surfaceScalarField & alphaPhi2
fv::options & fvOptions
rho
Definition: pEqn.H:1
pimpleNoLoopControl & pimple
surfaceScalarField & phi2
const volScalarField rAUKd2(rAU2 *Kd)
multiphaseSystem & fluid
Definition: createFields.H:11
p
Definition: pEqn.H:50
const surfaceScalarField alpharAUf2(fvc::interpolate(max(alpha2, phase2.residualAlpha()) *rAU2))
surfaceScalarField Df1(fvc::interpolate(rAU1 *(D+phase1.turbulence().pPrime())))
IOMRFZoneList & MRF
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
contErr2
phiHbyA
Definition: pEqn.H:20
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
engineTime & runTime
const surfaceScalarField & alphaPhi1
surfaceScalarField & phi1
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
const volScalarField Kd(fluid.Kd())
forAll(phases, phasei)
Definition: pEqn.H:3
const surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1)
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
volScalarField & dpdt
dynamicFvMesh & mesh
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
const volScalarField & alpha1
rhoEqn solve()
volVectorField F(fluid.F())
K2
Definition: pEqn.H:410
volVectorField & U1
alpha2
Definition: alphaEqn.H:115
phaseModel & phase1
p_rgh
Definition: pEqn.H:152
dimensioned< scalar > magSqr(const dimensioned< Type > &)
dimensionedScalar pos0(const dimensionedScalar &ds)
const surfaceScalarField & ghf
dimensionedScalar pMin("pMin", dimPressure, fluid)
const volScalarField & rAU2
Definition: pEqn.H:30
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.
mixture thermo2().correctRho(psi2 *(p_rgh - p_rgh_0))
K1
Definition: pEqn.H:409
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
surfaceScalarField Df2(fvc::interpolate(rAU2 *(D+phase2.turbulence().pPrime())))
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
label patchi
psi1
Definition: TEqns.H:34
U
Definition: pEqn.H:72
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho) *mesh.magSf())
surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1) *mesh.magSf())
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryFieldRef(),(phiHbyA.boundaryField() - MRF.relative(phib))/(mesh.magSf().boundaryField() *rAUf.boundaryField()))
const volScalarField & rAU1
Definition: pEqn.H:29
volVectorField & U2
volScalarField p_rgh_0(p_rgh)
const surfaceScalarField & phiF1
Definition: pEqn.H:50
const volScalarField & gh
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 > &)
const dimensionedScalar & rho2
Definition: createFields.H:40
phi
Definition: pEqn.H:18
const volScalarField rAUKd1(rAU1 *Kd)
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 dimensionedVector & g
const dimensionedScalar & rho1
Definition: createFields.H:39
const surfaceScalarField & phiF2
Definition: pEqn.H:51
const surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1))
const surfaceScalarField alpharAUf1(fvc::interpolate(max(alpha1, phase1.residualAlpha()) *rAU1))
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
zeroField Sp
Definition: alphaSuSp.H:2
mixture thermo1().correctRho(psi1 *(p_rgh - p_rgh_0))