pEqn.H
Go to the documentation of this file.
1surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1));
2surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1);
3
4surfaceScalarField alphaRhof10
5(
6 "alphaRhof10",
7 fvc::interpolate
8 (
9 max(alpha1.oldTime(), phase1.residualAlpha())
10 *rho1.oldTime()
11 )
12);
13
14surfaceScalarField alphaRhof20
15(
16 "alphaRhof20",
17 fvc::interpolate
18 (
19 max(alpha2.oldTime(), phase2.residualAlpha())
20 *rho2.oldTime()
21 )
22);
23
24// Drag coefficient
25surfaceScalarField Kdf("Kdf", fluid.Kdf());
26
27// Virtual-mass coefficient
28surfaceScalarField Vmf("Vmf", fluid.Vmf());
29
30surfaceScalarField rAUf1
31(
32 IOobject::groupName("rAUf", phase1.name()),
33 1.0
34 /(
35 (alphaRhof10 + Vmf)/runTime.deltaT()
36 + fvc::interpolate(U1Eqn.A())
37 + Kdf
38 )
39);
40
41surfaceScalarField rAUf2
42(
43 IOobject::groupName("rAUf", phase2.name()),
44 1.0
45 /(
46 (alphaRhof20 + Vmf)/runTime.deltaT()
47 + fvc::interpolate(U2Eqn.A())
48 + Kdf
49 )
50);
51
52
53// Turbulent dispersion, particle-pressure, lift and wall-lubrication forces
54tmp<surfaceScalarField> Ff1;
55tmp<surfaceScalarField> Ff2;
56{
57 // Turbulent-dispersion diffusivity
58 volScalarField D(fluid.D());
59
60 // Phase-1 turbulent dispersion and particle-pressure diffusivity
61 surfaceScalarField Df1
62 (
63 fvc::interpolate(D + phase1.turbulence().pPrime())
64 );
65
66 // Phase-2 turbulent dispersion and particle-pressure diffusivity
67 surfaceScalarField Df2
68 (
69 fvc::interpolate(D + phase2.turbulence().pPrime())
70 );
71
72 // Cache the net diffusivity for implicit diffusion treatment in the
73 // phase-fraction equation
74 if (implicitPhasePressure)
75 {
76 fluid.pPrimeByA() = rAUf1*Df1 + rAUf2*Df2;
77 }
78
79 // Lift and wall-lubrication forces
80 surfaceScalarField Ff(fluid.Ff());
81
82 // Phase-fraction face-gradient
83 surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
84
85 // Phase-1 dispersion, lift and wall-lubrication force
87
88 // Phase-2 dispersion, lift and wall-lubrication force
90}
91
92
93while (pimple.correct())
94{
95 // Update continuity errors due to temperature changes
96 #include "correctContErrs.H"
97
98 volScalarField rho("rho", fluid.rho());
99
100 // Correct p_rgh for consistency with p and the updated densities
101 p_rgh = p - rho*gh;
102
103 surfaceScalarField rhof1(fvc::interpolate(rho1));
104 surfaceScalarField rhof2(fvc::interpolate(rho2));
105
106 // Correct fixed-flux BCs to be consistent with the velocity BCs
107 MRF.correctBoundaryFlux(U1, phi1);
108 MRF.correctBoundaryFlux(U2, phi2);
109
110 surfaceScalarField alpharAUf1
111 (
112 IOobject::groupName("alpharAUf", phase1.name()),
113 max(alphaf1, phase1.residualAlpha())*rAUf1
114 );
115
116 surfaceScalarField alpharAUf2
117 (
118 IOobject::groupName("alpharAUf", phase2.name()),
119 max(alphaf2, phase2.residualAlpha())*rAUf2
120 );
121
122 surfaceScalarField ghSnGradRho
123 (
124 "ghSnGradRho",
125 ghf*fvc::snGrad(rho)*mesh.magSf()
126 );
127
128 // Phase-1 buoyancy flux
129 surfaceScalarField phig1
130 (
131 IOobject::groupName("phig", phase1.name()),
133 *(
135 - alphaf2*(rhof1 - rhof2)*(g & mesh.Sf())
136 )
137 );
138
139 // Phase-2 buoyancy flux
140 surfaceScalarField phig2
141 (
142 IOobject::groupName("phig", phase2.name()),
144 *(
146 - alphaf1*(rhof2 - rhof1)*(g & mesh.Sf())
147 )
148 );
149
150
151 // Phase-1 predicted flux
152 surfaceScalarField phiHbyA1
153 (
154 IOobject::groupName("phiHbyA", phase1.name()),
155 phi1
156 );
157
158 phiHbyA1 =
159 rAUf1
160 *(
161 (alphaRhof10 + Vmf)
162 *MRF.absolute(phi1.oldTime())/runTime.deltaT()
163 + fvc::flux(U1Eqn.H())
164 + Vmf*ddtPhi2
165 + Kdf*MRF.absolute(phi2)
166 - Ff1()
167 );
168
169 // Phase-2 predicted flux
170 surfaceScalarField phiHbyA2
171 (
172 IOobject::groupName("phiHbyA", phase2.name()),
173 phi2
174 );
175
176 phiHbyA2 =
177 rAUf2
178 *(
179 (alphaRhof20 + Vmf)
180 *MRF.absolute(phi2.oldTime())/runTime.deltaT()
181 + fvc::flux(U2Eqn.H())
182 + Vmf*ddtPhi1
183 + Kdf*MRF.absolute(phi1)
184 - Ff2()
185 );
186
187
188 surfaceScalarField phiHbyA
189 (
190 "phiHbyA",
191 alphaf1*(phiHbyA1 - phig1) + alphaf2*(phiHbyA2 - phig2)
192 );
193 MRF.makeRelative(phiHbyA);
194
195 phiHbyA1 -= phig1;
196 phiHbyA2 -= phig2;
197
198 surfaceScalarField rAUf
199 (
200 "rAUf",
202 );
203
204 // Update the fixedFluxPressure BCs to ensure flux consistency
206 (
207 p_rgh.boundaryFieldRef(),
208 (
209 phiHbyA.boundaryField()
210 - (
211 alphaf1.boundaryField()*phi1.boundaryField()
212 + alphaf2.boundaryField()*phi2.boundaryField()
213 )
214 )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
215 );
216
217 tmp<fvScalarMatrix> pEqnComp1;
218 tmp<fvScalarMatrix> pEqnComp2;
219
220 if (pimple.transonic())
221 {
222 surfaceScalarField phid1
223 (
224 IOobject::groupName("phid", phase1.name()),
225 fvc::interpolate(psi1)*phi1
226 );
227 surfaceScalarField phid2
228 (
229 IOobject::groupName("phid", phase2.name()),
230 fvc::interpolate(psi2)*phi2
231 );
232
233 pEqnComp1 =
234 (
236 - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
237 )/rho1
238 + correction
239 (
240 (alpha1/rho1)*
241 (
242 psi1*fvm::ddt(p_rgh)
243 + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
244 )
245 );
246 pEqnComp1.ref().relax();
247
248 pEqnComp2 =
249 (
251 - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
252 )/rho2
253 + correction
254 (
255 (alpha2/rho2)*
256 (
257 psi2*fvm::ddt(p_rgh)
258 + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
259 )
260 );
261 pEqnComp2.ref().relax();
262 }
263 else
264 {
265 pEqnComp1 =
266 (
268 - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
269 )/rho1
270 + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
271
272 pEqnComp2 =
273 (
275 - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
276 )/rho2
277 + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
278 }
279
280 // Cache p prior to solve for density update
281 volScalarField p_rgh_0("p_rgh_0", p_rgh);
282
283 while (pimple.correctNonOrthogonal())
284 {
285 fvScalarMatrix pEqnIncomp
286 (
287 fvc::div(phiHbyA)
288 - fvm::laplacian(rAUf, p_rgh)
289 );
290
291 solve
292 (
293 pEqnComp1() + pEqnComp2() + pEqnIncomp,
294 mesh.solver(p_rgh.select(pimple.finalInnerIter()))
295 );
296
297 if (pimple.finalNonOrthogonalIter())
298 {
299 surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
300
301 phi = phiHbyA + pEqnIncomp.flux();
302
303 surfaceScalarField phi1s
304 (
305 phiHbyA1
306 + alpharAUf1*mSfGradp
307 - rAUf1*Kdf*MRF.absolute(phi2)
308 );
309
310 surfaceScalarField phi2s
311 (
312 phiHbyA2
313 + alpharAUf2*mSfGradp
314 - rAUf2*Kdf*MRF.absolute(phi1)
315 );
316
317 surfaceScalarField phir
318 (
319 ((phi2s + rAUf2*Kdf*phi1s) - (phi1s + rAUf1*Kdf*phi2s))
320 /(1.0 - rAUf1*rAUf2*sqr(Kdf))
321 );
322
323 phi1 = phi - alphaf2*phir;
324 phi2 = phi + alphaf1*phir;
325
326 U1 = fvc::reconstruct(MRF.absolute(phi1));
327 U1.correctBoundaryConditions();
328 fvOptions.correct(U1);
329
330 U2 = fvc::reconstruct(MRF.absolute(phi2));
331 U2.correctBoundaryConditions();
332 fvOptions.correct(U2);
333
334 U = fluid.U();
335
336 fluid.dgdt() =
337 (
338 alpha1*(pEqnComp2 & p_rgh)
339 - alpha2*(pEqnComp1 & p_rgh)
340 );
341 }
342 }
343
344 // Update and limit the static pressure
345 p = max(p_rgh + rho*gh, pMin);
346
347 // Limit p_rgh
348 p_rgh = p - rho*gh;
349
350 // Update densities from change in p_rgh
351 rho1 += psi1*(p_rgh - p_rgh_0);
352 rho2 += psi2*(p_rgh - p_rgh_0);
353
354 // Correct p_rgh for consistency with p and the updated densities
355 rho = fluid.rho();
356 p_rgh = p - rho*gh;
357 p_rgh.correctBoundaryConditions();
358}
359
360// Update the phase kinetic energies
361K1 = 0.5*magSqr(U1);
362K2 = 0.5*magSqr(U2);
363
364// Update the pressure time-derivative if required
365if (thermo1.dpdt() || thermo2.dpdt())
366{
367 dpdt = fvc::ddt(p);
368}
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 alphaf2("alphaf2", scalar(1) - alphaf1)
const surfaceScalarField alpharAUf1(fvc::interpolate(max(alpha1, phase1.residualAlpha()) *rAU1))
const surfaceScalarField alpharAUf2(fvc::interpolate(max(alpha2, phase2.residualAlpha()) *rAU2))
const surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1))
const surfaceScalarField & rAUf1
Definition: pEqn.H:57
const surfaceScalarField & rAUf2
Definition: pEqn.H:58
surfaceScalarField alphaRhof20("alphaRhof20", fvc::interpolate(max(alpha2.oldTime(), phase2.residualAlpha()) *rho2.oldTime()))
surfaceScalarField alphaRhof10("alphaRhof10", fvc::interpolate(max(alpha1.oldTime(), phase1.residualAlpha()) *rho1.oldTime()))
const surfaceScalarField Kdf("Kdf", fluid.Kdf())
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())))
tmp< surfaceScalarField > Ff2
Definition: pEqn.H:55
surfaceScalarField Vmf("Vmf", fluid.Vmf())
tmp< surfaceScalarField > Ff1
Definition: pEqn.H:54
surfaceScalarField Ff(fluid.Ff())
const dimensionedScalar & D
CEqn solve()
ddtPhi1
Definition: DDtU.H:1
ddtPhi2
Definition: DDtU.H:6