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
74 
75  // Phase-1 dispersion, lift and wall-lubrication flux
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  (
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  (
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
290 
291  pEqnComp2 =
292  (
293  contErr2
295  )/rho2
297  }
298 
299  // Cache p prior to solve for density update
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  (
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 }
phiF2
const surfaceScalarField & phiF2
Definition: pEqn.H:51
runTime
engineTime & runTime
Definition: createEngineTime.H:13
Foam::fvc::reconstruct
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcReconstruct.C:56
Foam::fvc::snGrad
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:47
U1
volVectorField & U1
Definition: setRegionFluidFields.H:11
contErr2
contErr2
Definition: correctContErrs.H:5
Foam::fvc::flux
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
ghf
const surfaceScalarField & ghf
Definition: setRegionFluidFields.H:18
pMin
const dimensionedScalar & pMin
Definition: setRegionFluidFields.H:58
alpha2
const volScalarField & alpha2
Definition: setRegionFluidFields.H:9
alphaf2
const surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1)
contErr1
contErr1
Definition: correctContErrs.H:1
correctContErrs.H
phi1
surfaceScalarField & phi1
Definition: setRegionFluidFields.H:12
phiHbyA
phiHbyA
Definition: pEqn.H:20
Sp
zeroField Sp
Definition: alphaSuSp.H:2
rAU1
const volScalarField & rAU1
Definition: pEqn.H:29
rAU2
const volScalarField & rAU2
Definition: pEqn.H:30
phiF1
const surfaceScalarField & phiF1
Definition: pEqn.H:50
F
volVectorField F(fluid.F())
Foam::correction
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Foam::pos0
dimensionedScalar pos0(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:188
alpha1
const volScalarField & alpha1
Definition: setRegionFluidFields.H:8
fluid
twoPhaseSystem & fluid
Definition: setRegionFluidFields.H:3
Kd
const volScalarField Kd(fluid.Kd())
gh
const volScalarField & gh
Definition: setRegionFluidFields.H:17
Foam::fac::div
tmp< GeometricField< Type, faPatchField, areaMesh > > div(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Definition: facDiv.C:50
psi2
const volScalarField & psi2
Definition: setRegionFluidFields.H:31
phi2
surfaceScalarField & phi2
Definition: setRegionFluidFields.H:16
solve
CEqn solve()
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
forAll
forAll(phases, phasei)
Definition: pEqn.H:3
MRF
IOMRFZoneList & MRF
Definition: setRegionFluidFields.H:22
pimple
pimpleControl & pimple
Definition: setRegionFluidFields.H:56
rAUKd1
const volScalarField rAUKd1(rAU1 *Kd)
phir
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
K1
K1
Definition: pEqn.H:409
Foam::fvScalarMatrix
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:44
fvOptions
fv::options & fvOptions
Definition: setRegionFluidFields.H:23
Foam::fac::average
tmp< GeometricField< Type, faPatchField, areaMesh > > average(const GeometricField< Type, faePatchField, edgeMesh > &ssf)
Area-weighted average a edgeField creating a areaField.
Definition: facAverage.C:47
alphaf1
const surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1))
rho2
volScalarField & rho2
Definition: setRegionFluidFields.H:30
Foam::volScalarField
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
thermo1
mixture thermo1().correctRho(psi1 *(p_rgh - p_rgh_0))
rho1
volScalarField & rho1
Definition: setRegionFluidFields.H:27
U2
volVectorField & U2
Definition: setRegionFluidFields.H:15
phase2
phaseModel & phase2
Definition: setRegionFluidFields.H:6
p
p
Definition: pEqn.H:50
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
g
const uniformDimensionedVectorField & g
Definition: createFluidFields.H:24
alphaPhi1
const surfaceScalarField & alphaPhi1
Definition: setRegionFluidFields.H:13
psi1
const volScalarField & psi1
Definition: setRegionFluidFields.H:28
snGradAlpha1
surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1) *mesh.magSf())
K2
K2
Definition: pEqn.H:410
Foam::volVectorField
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:62
U
U
Definition: pEqn.H:72
p_rgh_0
volScalarField p_rgh_0(p_rgh)
Foam::surfaceScalarField
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Definition: surfaceFieldsFwd.H:54
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::fac::ddt
tmp< GeometricField< Type, faPatchField, areaMesh > > ddt(const dimensioned< Type > dt, const faMesh &mesh)
Definition: facDdt.C:47
ghSnGradRho
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho) *mesh.magSf())
phi
phi
Definition: pEqn.H:18
setSnGrad< fixedFluxPressureFvPatchScalarField >
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryFieldRef(),(phiHbyA.boundaryField() - MRF.relative(phib))/(mesh.magSf().boundaryField() *rAUf.boundaryField()))
Df1
surfaceScalarField Df1(fvc::interpolate(rAU1 *(D+phase1.turbulence().pPrime())))
D
const dimensionedScalar & D
Definition: solveBulkSurfactant.H:4
thermo2
mixture thermo2().correctRho(psi2 *(p_rgh - p_rgh_0))
alpharAUf2
const surfaceScalarField alpharAUf2(fvc::interpolate(max(alpha2, phase2.residualAlpha()) *rAU2))
Df2
surfaceScalarField Df2(fvc::interpolate(rAU2 *(D+phase2.turbulence().pPrime())))
dpdt
volScalarField & dpdt
Definition: setRegionFluidFields.H:32
phase1
phaseModel & phase1
Definition: setRegionFluidFields.H:5
p_rgh
p_rgh
Definition: pEqn.H:139
Foam::fac::laplacian
tmp< GeometricField< Type, faPatchField, areaMesh > > laplacian(const GeometricField< Type, faPatchField, areaMesh > &vf, const word &name)
Definition: facLaplacian.C:47
rAUKd2
const volScalarField rAUKd2(rAU2 *Kd)
alpharAUf1
const surfaceScalarField alpharAUf1(fvc::interpolate(max(alpha1, phase1.residualAlpha()) *rAU1))
rho
rho
Definition: pEqn.H:1
interpolate
mesh interpolate(rAU)
alphaPhi2
const surfaceScalarField & alphaPhi2
Definition: setRegionFluidFields.H:17
rAUf
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))