pEqn.H
Go to the documentation of this file.
2 const surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1);
3 
4 PtrList<volScalarField> rAUs;
5 rAUs.append
6 (
8  (
9  IOobject::groupName("rAU", phase1.name()),
10  1.0
11  /(
12  U1Eqn.A()
13  + byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1)
14  )
15  )
16 );
17 rAUs.append
18 (
19  new volScalarField
20  (
21  IOobject::groupName("rAU", phase2.name()),
22  1.0
23  /(
24  U2Eqn.A()
25  + byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2)
26  )
27  )
28 );
29 const volScalarField& rAU1 = rAUs[0];
30 const volScalarField& rAU2 = rAUs[1];
31 
33 (
34  fvc::interpolate(max(alpha1, phase1.residualAlpha())*rAU1)
35 );
37 (
38  fvc::interpolate(max(alpha2, phase2.residualAlpha())*rAU2)
39 );
40 
41 // Drag coefficients
42 const volScalarField Kd(fluid.Kd());
47 
48 // Explicit force fluxes
49 PtrList<surfaceScalarField> phiFs(fluid.phiFs(rAUs));
52 
53 // --- Pressure corrector loop
54 while (pimple.correct())
55 {
56  volScalarField rho("rho", fluid.rho());
57 
58  // Correct p_rgh for consistency with p and the updated densities
59  p_rgh = p - rho*gh;
60 
61  // Correct fixed-flux BCs to be consistent with the velocity BCs
62  MRF.correctBoundaryFlux(U1, phi1);
63  MRF.correctBoundaryFlux(U2, phi2);
64 
65  // Combined buoyancy and force fluxes
67  (
68  "ghSnGradRho",
69  ghf*fvc::snGrad(rho)*mesh.magSf()
70  );
71 
72  const surfaceScalarField phigF1
73  (
75  *(
77  - alphaf2*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
78  )
79  + phiF1
80  );
81 
82  const surfaceScalarField phigF2
83  (
85  *(
87  - alphaf1*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
88  )
89  + phiF2
90  );
91 
92  // Predicted velocities
93  volVectorField HbyA1
94  (
95  IOobject::groupName("HbyA", phase1.name()),
96  U1
97  );
98  HbyA1 =
99  rAU1
100  *(
101  U1Eqn.H()
102  + byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1)
103  *U1.oldTime()
104  );
105 
106  volVectorField HbyA2
107  (
108  IOobject::groupName("HbyA", phase2.name()),
109  U2
110  );
111  HbyA2 =
112  rAU2
113  *(
114  U2Eqn.H()
115  + byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2)
116  *U2.oldTime()
117  );
118 
119  // Correction force fluxes
120  PtrList<surfaceScalarField> ddtCorrByAs(fluid.ddtCorrByAs(rAUs));
121 
122  // Predicted fluxes
123  const surfaceScalarField phiHbyA1
124  (
125  IOobject::groupName("phiHbyA", phase1.name()),
126  fvc::flux(HbyA1) - phigF1 - ddtCorrByAs[0]
127  );
128 
129  const surfaceScalarField phiHbyA2
130  (
131  IOobject::groupName("phiHbyA", phase2.name()),
132  fvc::flux(HbyA2) - phigF2 - ddtCorrByAs[1]
133  );
134 
135  ddtCorrByAs.clear();
136 
137  // Drag fluxes
138  PtrList<surfaceScalarField> phiKdPhis(fluid.phiKdPhis(rAUs));
139  const surfaceScalarField& phiKdPhi1 = phiKdPhis[0];
140  const surfaceScalarField& phiKdPhi2 = phiKdPhis[1];
141 
142  // Total predicted flux
144  (
145  "phiHbyA",
146  alphaf1*(phiHbyA1 - phiKdPhi1) + alphaf2*(phiHbyA2 - phiKdPhi2)
147  );
148 
149  MRF.makeRelative(phiHbyA);
150 
151  phiKdPhis.clear();
152 
153  // Construct pressure "diffusivity"
155  (
156  "rAUf",
158  );
159 
160  // Update the fixedFluxPressure BCs to ensure flux consistency
162  (
163  p_rgh.boundaryFieldRef(),
164  (
165  phiHbyA.boundaryField()
166  - (
167  alphaf1.boundaryField()*phi1.boundaryField()
168  + alphaf2.boundaryField()*phi2.boundaryField()
169  )
170  )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
171  );
172 
173  // Construct the compressibility parts of the pressure equation
174  tmp<fvScalarMatrix> pEqnComp1, pEqnComp2;
175  if (phase1.compressible())
176  {
177  if (pimple.transonic())
178  {
179  const surfaceScalarField phid1
180  (
181  IOobject::groupName("phid", phase1.name()),
183  );
184 
185  pEqnComp1 =
186  (
187  fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi())
189  )/rho1
190  + correction
191  (
192  (alpha1/rho1)*
193  (
195  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
196  )
197  );
198 
199  pEqnComp1.ref().relax();
200  }
201  else
202  {
203  pEqnComp1 =
204  (
205  fvc::ddt(alpha1, rho1) + fvc::div(phase1.alphaRhoPhi())
207  )/rho1
209  }
210  }
211  if (phase2.compressible())
212  {
213  if (pimple.transonic())
214  {
215  const surfaceScalarField phid2
216  (
217  IOobject::groupName("phid", phase2.name()),
219  );
220 
221  pEqnComp2 =
222  (
223  fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi())
225  )/rho2
226  + correction
227  (
228  (alpha2/rho2)*
229  (
231  + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
232  )
233  );
234 
235  pEqnComp2.ref().relax();
236  }
237  else
238  {
239  pEqnComp2 =
240  (
241  fvc::ddt(alpha2, rho2) + fvc::div(phase2.alphaRhoPhi())
243  )/rho2
245  }
246  }
247 
248  // Add option sources
249  {
250  if (fvOptions.appliesToField(rho1.name()))
251  {
252  tmp<fvScalarMatrix> optEqn1 = fvOptions(alpha1, rho1);
253  if (pEqnComp1.valid())
254  {
255  pEqnComp1.ref() -= (optEqn1 & rho1)/rho1;
256  }
257  else
258  {
259  pEqnComp1 = fvm::Su(- (optEqn1 & rho1)/rho1, p_rgh);
260  }
261  }
262  if (fvOptions.appliesToField(rho2.name()))
263  {
264  tmp<fvScalarMatrix> optEqn2 = fvOptions(alpha2, rho2);
265  if (pEqnComp2.valid())
266  {
267  pEqnComp2.ref() -= (optEqn2 & rho2)/rho2;
268  }
269  else
270  {
271  pEqnComp2 = fvm::Su(- (optEqn2 & rho2)/rho2, p_rgh);
272  }
273  }
274  }
275 
276  // Add mass transfer
277  {
278  PtrList<volScalarField> dmdts(fluid.dmdts());
279  if (dmdts.set(0))
280  {
281  if (pEqnComp1.valid())
282  {
283  pEqnComp1.ref() -= dmdts[0]/rho1;
284  }
285  else
286  {
287  pEqnComp1 = fvm::Su(- dmdts[0]/rho1, p_rgh);
288  }
289  }
290  if (dmdts.set(1))
291  {
292  if (pEqnComp2.valid())
293  {
294  pEqnComp2.ref() -= dmdts[1]/rho2;
295  }
296  else
297  {
298  pEqnComp2 = fvm::Su(- dmdts[1]/rho2, p_rgh);
299  }
300  }
301  }
302 
303  // Cache p prior to solve for density update
305 
306  // Iterate over the pressure equation to correct for non-orthogonality
307  while (pimple.correctNonOrthogonal())
308  {
309  // Construct the transport part of the pressure equation
310  fvScalarMatrix pEqnIncomp
311  (
314  );
315 
316  {
317  fvScalarMatrix pEqn(pEqnIncomp);
318 
319  if (pEqnComp1.valid())
320  {
321  pEqn += pEqnComp1();
322  }
323 
324  if (pEqnComp2.valid())
325  {
326  pEqn += pEqnComp2();
327  }
328 
329  pEqn.solve();
330  }
331 
332  // Correct fluxes and velocities on last non-orthogonal iteration
333  if (pimple.finalNonOrthogonalIter())
334  {
335  phi = phiHbyA + pEqnIncomp.flux();
336 
337  surfaceScalarField mSfGradp
338  (
339  "mSfGradp",
340  pEqnIncomp.flux()/rAUf
341  );
342 
343  // Partial-elimination phase-flux corrector
344  {
345  const surfaceScalarField phi1s
346  (
347  phiHbyA1 + alpharAUf1*mSfGradp
348  );
349 
350  const surfaceScalarField phi2s
351  (
352  phiHbyA2 + alpharAUf2*mSfGradp
353  );
354 
356  (
357  ((phi1s + phiKd1*phi2s) - (phi2s + phiKd2*phi1s))
358  /(1 - phiKd1*phiKd2)
359  );
360 
361  phi1 = phi + alphaf2*phir;
362  phi2 = phi - alphaf1*phir;
363  }
364 
365  // Set the phase dilatation rates
366  if (pEqnComp1.valid())
367  {
368  phase1.divU(-pEqnComp1 & p_rgh);
369  }
370  if (pEqnComp2.valid())
371  {
372  phase2.divU(-pEqnComp2 & p_rgh);
373  }
374 
375  // Optionally relax pressure for velocity correction
376  p_rgh.relax();
377 
378  mSfGradp = pEqnIncomp.flux()/rAUf;
379 
380  // Partial-elimination phase-velocity corrector
381  {
382  const volVectorField Us1
383  (
384  HbyA1
385  + fvc::reconstruct(alpharAUf1*mSfGradp - phigF1)
386  );
387 
388  const volVectorField Us2
389  (
390  HbyA2
391  + fvc::reconstruct(alpharAUf2*mSfGradp - phigF2)
392  );
393 
394  const volVectorField U
395  (
396  alpha1*(Us1 + rAUKd1*U2) + alpha2*(Us2 + rAUKd2*U1)
397  );
398 
399  const volVectorField Ur
400  (
401  ((1 - rAUKd2)*Us1 - (1 - rAUKd1)*Us2)/(1 - rAUKd1*rAUKd2)
402  );
403 
404  U1 = U + alpha2*Ur;
405  U1.correctBoundaryConditions();
406  fvOptions.correct(U1);
407 
408  U2 = U - alpha1*Ur;
409  U2.correctBoundaryConditions();
410  fvOptions.correct(U2);
411  }
412  }
413  }
414 
415  // Update and limit the static pressure
416  p = max(p_rgh + rho*gh, pMin);
417 
418  // Limit p_rgh
419  p_rgh = p - rho*gh;
420 
421  // Update densities from change in p_rgh
422  rho1 += psi1*(p_rgh - p_rgh_0);
423  rho2 += psi2*(p_rgh - p_rgh_0);
424 
425  // Correct p_rgh for consistency with p and the updated densities
426  rho = fluid.rho();
427  p_rgh = p - rho*gh;
428  p_rgh.correctBoundaryConditions();
429 }
phiF2
const surfaceScalarField & phiF2
Definition: pEqn.H:51
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
rAUs
PtrList< volScalarField > rAUs
Definition: pEqn.H:4
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)
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
Foam::correction
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
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
phiKd2
const surfaceScalarField phiKd2(fvc::interpolate(rAUKd2))
psi2
const volScalarField & psi2
Definition: setRegionFluidFields.H:31
phi2
surfaceScalarField & phi2
Definition: setRegionFluidFields.H:16
Su
zeroField Su
Definition: alphaSuSp.H:1
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()))
phiKd1
const surfaceScalarField phiKd1(fvc::interpolate(rAUKd1))
Foam::fvScalarMatrix
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:44
fvOptions
fv::options & fvOptions
Definition: setRegionFluidFields.H:23
phiFs
PtrList< surfaceScalarField > phiFs(fluid.phiFs(rAUs))
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
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:26
alphaPhi1
const surfaceScalarField & alphaPhi1
Definition: setRegionFluidFields.H:13
psi1
const volScalarField & psi1
Definition: setRegionFluidFields.H:28
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()))
alpharAUf2
const surfaceScalarField alpharAUf2(fvc::interpolate(max(alpha2, phase2.residualAlpha()) *rAU2))
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
Foam::byDt
tmp< volScalarField > byDt(const volScalarField &vf)
Definition: phaseSystem.C:438
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))