pEqn.H
Go to the documentation of this file.
1// Face volume fractions
2PtrList<surfaceScalarField> alphafs(phases.size());
3PtrList<surfaceScalarField> alphaRho0fs(phases.size());
5{
6 phaseModel& phase = phases[phasei];
7 const volScalarField& alpha = phase;
8
9 alphafs.set(phasei, fvc::interpolate(alpha).ptr());
10 alphafs[phasei].rename("pEqn" + alphafs[phasei].name());
11
12 alphaRho0fs.set
13 (
14 phasei,
15 (
16 fvc::interpolate
17 (
18 max(alpha.oldTime(), phase.residualAlpha())
19 *phase.rho()().oldTime()
20 )
21 ).ptr()
22 );
23}
24
25// Diagonal coefficients
26PtrList<surfaceScalarField> rAUfs(phases.size());
27{
28 PtrList<surfaceScalarField> AFfs(fluid.AFfs());
29
30 forAll(fluid.movingPhases(), movingPhasei)
31 {
32 phaseModel& phase = fluid.movingPhases()[movingPhasei];
33
34 rAUfs.set
35 (
36 phase.index(),
37 new surfaceScalarField
38 (
39 IOobject::groupName("rAUf", phase.name()),
40 1.0
41 /(
42 byDt(alphaRho0fs[phase.index()])
43 + fvc::interpolate(UEqns[phase.index()].A())
44 + AFfs[phase.index()]
45 )
46 )
47 );
48 }
49}
50fluid.fillFields("rAUf", dimTime/dimDensity, rAUfs);
51
52// Phase diagonal coefficients
53PtrList<surfaceScalarField> alpharAUfs(phases.size());
55{
56 phaseModel& phase = phases[phasei];
57 alpharAUfs.set
58 (
59 phase.index(),
60 (
61 max(alphafs[phase.index()], phase.residualAlpha())
62 *rAUfs[phase.index()]
63 ).ptr()
64 );
65}
66
67// Explicit force fluxes
68PtrList<surfaceScalarField> phiFfs(fluid.phiFfs(rAUfs));
69
70// --- Pressure corrector loop
71while (pimple.correct())
72{
73 volScalarField rho("rho", fluid.rho());
74
75 // Correct p_rgh for consistency with p and the updated densities
76 p_rgh = p - rho*gh;
77
78 // Correct fixed-flux BCs to be consistent with the velocity BCs
79 forAll(fluid.movingPhases(), movingPhasei)
80 {
81 phaseModel& phase = fluid.movingPhases()[movingPhasei];
82 MRF.correctBoundaryFlux(phase.U(), phase.phiRef());
83 }
84
85 // Combined buoyancy and force fluxes
86 PtrList<surfaceScalarField> phigFs(phases.size());
87 {
88 const surfaceScalarField ghSnGradRho
89 (
90 "ghSnGradRho",
91 ghf*fvc::snGrad(rho)*mesh.magSf()
92 );
93
95 {
96 phaseModel& phase = phases[phasei];
97
98 phigFs.set
99 (
100 phasei,
101 (
103 *(
105 - (fvc::interpolate(phase.rho() - rho))*(g & mesh.Sf())
106 - fluid.surfaceTension(phase)*mesh.magSf()
107 )
108 ).ptr()
109 );
110
111 if (phiFfs.set(phasei))
112 {
113 phigFs[phasei] += phiFfs[phasei];
114 }
115 }
116 }
117
118 // Predicted fluxes for each phase
119 PtrList<surfaceScalarField> phiHbyAs(phases.size());
120 forAll(fluid.movingPhases(), movingPhasei)
121 {
122 phaseModel& phase = fluid.movingPhases()[movingPhasei];
123
124 phiHbyAs.set
125 (
126 phase.index(),
127 new surfaceScalarField
128 (
129 IOobject::groupName("phiHbyA", phase.name()),
130 rAUfs[phase.index()]
131 *(
132 fvc::flux(UEqns[phase.index()].H())
133 + alphaRho0fs[phase.index()]
134 *byDt(MRF.absolute(phase.phi()().oldTime()))
135 )
136 - phigFs[phase.index()]
137 )
138 );
139 }
140 fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs);
141
142 // Add explicit drag forces and fluxes if not doing partial elimination
143 if (!partialElimination)
144 {
145 PtrList<surfaceScalarField> phiKdPhifs(fluid.phiKdPhifs(rAUfs));
146
148 {
149 if (phiKdPhifs.set(phasei))
150 {
151 phiHbyAs[phasei] -= phiKdPhifs[phasei];
152 }
153 }
154 }
155
156 // Total predicted flux
157 surfaceScalarField phiHbyA
158 (
159 IOobject
160 (
161 "phiHbyA",
162 runTime.timeName(),
163 mesh,
164 IOobject::NO_READ,
165 IOobject::AUTO_WRITE
166 ),
167 mesh,
168 dimensionedScalar("phiHbyA", dimArea*dimVelocity, 0)
169 );
170
172 {
174 }
175
176 // Add explicit drag fluxes if doing partial elimination
177 if (partialElimination)
178 {
179 PtrList<surfaceScalarField> phiKdPhifs(fluid.phiKdPhifs(rAUfs));
180
182 {
183 if (phiKdPhifs.set(phasei))
184 {
185 phiHbyA -= alphafs[phasei]*phiKdPhifs[phasei];
186 }
187 }
188 }
189
190 MRF.makeRelative(phiHbyA);
191
192 // Construct pressure "diffusivity"
193 surfaceScalarField rAUf
194 (
195 IOobject
196 (
197 "rAUf",
198 runTime.timeName(),
199 mesh
200 ),
201 mesh,
202 dimensionedScalar("rAUf", dimensionSet(-1, 3, 1, 0, 0), 0)
203 );
204
206 {
208 }
209
210 rAUf = mag(rAUf);
211
212 // Update the fixedFluxPressure BCs to ensure flux consistency
213 {
214 surfaceScalarField::Boundary phib(phi.boundaryField());
215 phib = 0;
217 {
218 phaseModel& phase = phases[phasei];
219 phib +=
220 alphafs[phasei].boundaryField()*phase.phi()().boundaryField();
221 }
222
224 (
225 p_rgh.boundaryFieldRef(),
226 (
227 phiHbyA.boundaryField() - phib
228 )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
229 );
230 }
231
232 // Compressible pressure equations
233 PtrList<fvScalarMatrix> pEqnComps(phases.size());
234 PtrList<volScalarField> dmdts(fluid.dmdts());
236 {
237 phaseModel& phase = phases[phasei];
238 const volScalarField& alpha = phase;
239 volScalarField& rho = phase.thermoRef().rho();
240
241 if (phase.compressible())
242 {
243 if (pimple.transonic())
244 {
245 surfaceScalarField phid
246 (
247 IOobject::groupName("phid", phase.name()),
248 fvc::interpolate(phase.thermo().psi())*phase.phi()
249 );
250
251 pEqnComps.set
252 (
253 phasei,
254 (
255 (
256 fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
257 - fvc::Sp
258 (
259 fvc::ddt(alpha) + fvc::div(phase.alphaPhi()),
260 rho
261 )
262 )/rho
263 + correction
264 (
265 (alpha/rho)*
266 (
267 phase.thermo().psi()*fvm::ddt(p_rgh)
268 + fvm::div(phid, p_rgh)
269 - fvm::Sp(fvc::div(phid), p_rgh)
270 )
271 )
272 ).ptr()
273 );
274
275 deleteDemandDrivenData
276 (
277 pEqnComps[phasei].faceFluxCorrectionPtr()
278 );
279
280 pEqnComps[phasei].relax();
281 }
282 else
283 {
284 pEqnComps.set
285 (
286 phasei,
287 (
288 (
289 fvc::ddt(alpha, rho) + fvc::div(phase.alphaRhoPhi())
290 - fvc::Sp
291 (
292 (fvc::ddt(alpha) + fvc::div(phase.alphaPhi())),
293 rho
294 )
295 )/rho
296 + (alpha*phase.thermo().psi()/rho)
297 *correction(fvm::ddt(p_rgh))
298 ).ptr()
299 );
300 }
301 }
302
303 if (fvOptions.appliesToField(rho.name()))
304 {
305 tmp<fvScalarMatrix> optEqn = fvOptions(alpha, rho);
306 if (pEqnComps.set(phasei))
307 {
308 pEqnComps[phasei] -= (optEqn&rho)/rho;
309 }
310 else
311 {
312 pEqnComps.set
313 (
314 phasei,
315 fvm::Su(- (optEqn&rho)/rho, p_rgh).ptr()
316 );
317 }
318 }
319
320 if (dmdts.set(phasei))
321 {
322 if (pEqnComps.set(phasei))
323 {
324 pEqnComps[phasei] -= dmdts[phasei]/rho;
325 }
326 else
327 {
328 pEqnComps.set
329 (
330 phasei,
331 fvm::Su(- dmdts[phasei]/rho, p_rgh)
332 );
333 }
334 }
335 }
336
337 // Cache p prior to solve for density update
338 volScalarField p_rgh_0(p_rgh);
339
340 // Iterate over the pressure equation to correct for non-orthogonality
341 while (pimple.correctNonOrthogonal())
342 {
343 // Construct the transport part of the pressure equation
344 fvScalarMatrix pEqnIncomp
345 (
346 fvc::div(phiHbyA)
347 - fvm::laplacian(rAUf, p_rgh)
348 );
349
350 {
351 fvScalarMatrix pEqn(pEqnIncomp);
352
354 {
355 if (pEqnComps.set(phasei))
356 {
357 pEqn += pEqnComps[phasei];
358 }
359 }
360
361 solve
362 (
363 pEqn,
364 mesh.solver(p_rgh.select(pimple.finalInnerIter()))
365 );
366 }
367
368 // Correct fluxes and velocities on last non-orthogonal iteration
369 if (pimple.finalNonOrthogonalIter())
370 {
371 phi = phiHbyA + pEqnIncomp.flux();
372
373 surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
374
375 forAll(fluid.movingPhases(), movingPhasei)
376 {
377 phaseModel& phase = fluid.movingPhases()[movingPhasei];
378
379 phase.phiRef() =
380 phiHbyAs[phase.index()]
381 + alpharAUfs[phase.index()]*mSfGradp;
382
383 // Set the phase dilatation rates
384 if (pEqnComps.set(phase.index()))
385 {
386 phase.divU(-pEqnComps[phase.index()] & p_rgh);
387 }
388 }
389
390 if (partialElimination)
391 {
392 fluid.partialEliminationf(rAUfs);
393 }
394
395 // Optionally relax pressure for velocity correction
396 p_rgh.relax();
397
398 mSfGradp = pEqnIncomp.flux()/rAUf;
399
400 forAll(fluid.movingPhases(), movingPhasei)
401 {
402 phaseModel& phase = fluid.movingPhases()[movingPhasei];
403
404 phase.URef() = fvc::reconstruct(MRF.absolute(phase.phi()));
405 phase.URef().correctBoundaryConditions();
406 fvOptions.correct(phase.URef());
407 }
408 }
409 }
410
411 // Update and limit the static pressure
412 p = max(p_rgh + rho*gh, pMin);
413
414 // Limit p_rgh
415 p_rgh = p - rho*gh;
416
417 // Update densities from change in p_rgh
419 {
420 phaseModel& phase = phases[phasei];
421 phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
422 }
423
424 // Correct p_rgh for consistency with p and the updated densities
425 rho = fluid.rho();
426 p_rgh = p - rho*gh;
427 p_rgh.correctBoundaryConditions();
428}
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
twoPhaseSystem & fluid
const dimensionedScalar & pMin
volScalarField & p
phiHbyA
Definition: pcEqn.H:73
surfaceScalarField phid("phid", fvc::interpolate(psi) *(fvc::flux(HbyA)+MRF.zeroFilter(rhorAUf *fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho))))
dynamicFvMesh & mesh
engineTime & runTime
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
PtrList< fvVectorMatrix > UEqns(fluid.phases().size())
volScalarField p_rgh_0(p_rgh)
label phasei
Definition: pEqn.H:27
phib
Definition: pEqn.H:189
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho) *mesh.magSf())
PtrList< surfaceScalarField > phiHbyAs(fluid.phases().size())
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryFieldRef(),(phiHbyA.boundaryField() - MRF.relative(phib))/(mesh.magSf().boundaryField() *rAUf.boundaryField()))
PtrList< surfaceScalarField > alphafs(phases.size())
PtrList< surfaceScalarField > alpharAUfs(phases.size())
PtrList< surfaceScalarField > phiFfs(fluid.phiFfs(rAUfs))
PtrList< surfaceScalarField > alphaRho0fs(phases.size())
PtrList< surfaceScalarField > rAUfs
Definition: pEqn.H:30
PtrList< surfaceScalarField > AFfs(fluid.AFfs())
volScalarField & alpha
CEqn solve()
multiphaseSystem::phaseModelList & phases
Definition: createFields.H:12
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333