pEqn.H
Go to the documentation of this file.
1surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1));
2surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1);
3
4volScalarField rAU1
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);
14volScalarField rAU2
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
25surfaceScalarField alpharAUf1
26(
27 fvc::interpolate(max(alpha1, phase1.residualAlpha())*rAU1)
28);
29surfaceScalarField alpharAUf2
30(
31 fvc::interpolate(max(alpha2, phase2.residualAlpha())*rAU2)
32);
33
34// Turbulent diffusion, particle-pressure, lift and wall-lubrication fluxes
35tmp<surfaceScalarField> phiF1;
36tmp<surfaceScalarField> phiF2;
37{
38 // Turbulent-dispersion diffusivity
39 volScalarField D(fluid.D());
40
41 // Phase-1 turbulent dispersion and particle-pressure flux
42 surfaceScalarField Df1
43 (
44 fvc::interpolate
45 (
46 rAU1*(D + phase1.turbulence().pPrime())
47 )
48 );
49
50 // Phase-2 turbulent dispersion and particle-pressure flux
51 surfaceScalarField Df2
52 (
53 fvc::interpolate
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
70 surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
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
81while (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
121 surfaceScalarField ghSnGradRho
122 (
123 "ghSnGradRho",
124 ghf*fvc::snGrad(rho)*mesh.magSf()
125 );
126
127 surfaceScalarField phig1
128 (
130 *(
132 - alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
133 )
134 );
135
136 surfaceScalarField phig2
137 (
139 *(
141 - alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
142 )
143 );
144
145
146 // ddtPhiCorr filter -- only apply in pure(ish) phases
147 surfaceScalarField alphaf1Bar(fvc::interpolate(fvc::average(alphaf1)));
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
202 surfaceScalarField rAUKd1(fvc::interpolate(rAU1*Kd));
203 surfaceScalarField rAUKd2(fvc::interpolate(rAU2*Kd));
204
205 // Construct the mean predicted flux
206 // including explicit drag contributions based on absolute fluxes
207 surfaceScalarField phiHbyA
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"
216 surfaceScalarField rAUf
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()),
244 fvc::interpolate(psi1)*phi1
245 );
246 surfaceScalarField phid2
247 (
248 IOobject::groupName("phid", phase2.name()),
249 fvc::interpolate(psi2)*phi2
250 );
251
252 pEqnComp1 =
253 (
255 - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
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 (
270 - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
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 (
287 - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
288 )/rho1
289 + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
290
291 pEqnComp2 =
292 (
294 - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
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
337 surfaceScalarField phir
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
409K1 = 0.5*magSqr(U1);
410K2 = 0.5*magSqr(U2);
411
412// Update the pressure time-derivative if required
413if (thermo1.dpdt() || thermo2.dpdt())
414{
415 dpdt = fvc::ddt(p);
416}
Y[inertIndex] max(0.0)
const uniformDimensionedVectorField & g
volScalarField & p_rgh
fv::options & fvOptions
surfaceScalarField & phi
const surfaceScalarField & ghf
IOMRFZoneList & MRF
const volScalarField & gh
pimpleControl & pimple
rhoThermo & thermo2
const volScalarField & alpha1
rhoThermo & thermo1
const volScalarField & psi2
surfaceScalarField & phi2
twoPhaseSystem & fluid
phaseModel & phase1
volScalarField & rho2
const volScalarField & alpha2
const surfaceScalarField & alphaPhi1
surfaceScalarField & phi1
phaseModel & phase2
const volScalarField & psi1
volVectorField & U1
const dimensionedScalar & pMin
volScalarField & rho1
volVectorField & U2
const surfaceScalarField & alphaPhi2
U
Definition: pEqn.H:72
volScalarField & p
phiHbyA
Definition: pcEqn.H:73
contErr2
contErr1
dynamicFvMesh & mesh
engineTime & runTime
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
volScalarField & dpdt
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
volScalarField p_rgh_0(p_rgh)
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho) *mesh.magSf())
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryFieldRef(),(phiHbyA.boundaryField() - MRF.relative(phib))/(mesh.magSf().boundaryField() *rAUf.boundaryField()))
const surfaceScalarField & phiF1
Definition: pEqn.H:50
const surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1)
const surfaceScalarField alpharAUf1(fvc::interpolate(max(alpha1, phase1.residualAlpha()) *rAU1))
const volScalarField Kd(fluid.Kd())
const surfaceScalarField & phiF2
Definition: pEqn.H:51
const volScalarField & rAU2
Definition: pEqn.H:30
const volScalarField rAUKd2(rAU2 *Kd)
const surfaceScalarField alpharAUf2(fvc::interpolate(max(alpha2, phase2.residualAlpha()) *rAU2))
const volScalarField rAUKd1(rAU1 *Kd)
const surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1))
const volScalarField & rAU1
Definition: pEqn.H:29
volVectorField F(fluid.F())
K2
Definition: pEqn.H:410
surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1) *mesh.magSf())
surfaceScalarField Df1(fvc::interpolate(rAU1 *(D+phase1.turbulence().pPrime())))
K1
Definition: pEqn.H:409
surfaceScalarField Df2(fvc::interpolate(rAU2 *(D+phase2.turbulence().pPrime())))
const dimensionedScalar & D
CEqn solve()
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333