pEqn.H
Go to the documentation of this file.
2 surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1);
3 
5 (
6  "alphaRhof10",
8  (
9  max(alpha1.oldTime(), phase1.residualAlpha())
10  *rho1.oldTime()
11  )
12 );
13 
15 (
16  "alphaRhof20",
18  (
19  max(alpha2.oldTime(), phase2.residualAlpha())
20  *rho2.oldTime()
21  )
22 );
23 
24 // Drag coefficient
25 surfaceScalarField Kdf("Kdf", fluid.Kdf());
26 
27 // Virtual-mass coefficient
28 surfaceScalarField Vmf("Vmf", fluid.Vmf());
29 
31 (
32  IOobject::groupName("rAUf", phase1.name()),
33  1.0
34  /(
35  byDt(alphaRhof10 + Vmf)
36  + fvc::interpolate(U1Eqn.A())
37  + Kdf
38  )
39 );
40 
42 (
43  IOobject::groupName("rAUf", phase2.name()),
44  1.0
45  /(
46  byDt(alphaRhof20 + Vmf)
47  + fvc::interpolate(U2Eqn.A())
48  + Kdf
49  )
50 );
51 
52 
53 // Turbulent dispersion, particle-pressure, lift and wall-lubrication forces
54 tmp<surfaceScalarField> Ff1;
55 tmp<surfaceScalarField> Ff2;
56 {
57  // Turbulent-dispersion diffusivity
58  volScalarField D(fluid.D());
59 
60  // Phase-1 turbulent dispersion and particle-pressure diffusivity
62  (
63  fvc::interpolate(D + phase1.turbulence().pPrime())
64  );
65 
66  // Phase-2 turbulent dispersion and particle-pressure diffusivity
68  (
69  fvc::interpolate(D + phase2.turbulence().pPrime())
70  );
71 
72  // Cache the phase diffusivities for implicit treatment in the
73  // phase-fraction equation
74  if (implicitPhasePressure)
75  {
76  phase1.DbyA(rAUf1*Df1);
77  phase2.DbyA(rAUf2*Df2);
78  }
79 
80  // Lift and wall-lubrication forces
82 
83  // Phase-fraction face-gradient
85 
86  // Phase-1 dispersion, lift and wall-lubrication force
87  Ff1 = Df1*snGradAlpha1 + Ff;
88 
89  // Phase-2 dispersion, lift and wall-lubrication force
90  Ff2 = -Df2*snGradAlpha1 - Ff;
91 }
92 
93 
94 while (pimple.correct())
95 {
96  volScalarField rho("rho", fluid.rho());
97 
98  // Correct p_rgh for consistency with p and the updated densities
99  p_rgh = p - rho*gh;
100 
103 
104  // Correct fixed-flux BCs to be consistent with the velocity BCs
105  MRF.correctBoundaryFlux(U1, phi1);
106  MRF.correctBoundaryFlux(U2, phi2);
107 
109  (
110  IOobject::groupName("alpharAUf", phase1.name()),
111  max(alphaf1, phase1.residualAlpha())*rAUf1
112  );
113 
115  (
116  IOobject::groupName("alpharAUf", phase2.name()),
117  max(alphaf2, phase2.residualAlpha())*rAUf2
118  );
119 
121  (
122  "ghSnGradRho",
123  ghf*fvc::snGrad(rho)*mesh.magSf()
124  );
125 
126  // Phase-1 buoyancy flux
127  surfaceScalarField phig1
128  (
129  IOobject::groupName("phig", phase1.name()),
130  alpharAUf1
131  *(
133  - alphaf2*(rhof1 - rhof2)*(g & mesh.Sf())
134  )
135  );
136 
137  // Phase-2 buoyancy flux
138  surfaceScalarField phig2
139  (
140  IOobject::groupName("phig", phase2.name()),
141  alpharAUf2
142  *(
144  - alphaf1*(rhof2 - rhof1)*(g & mesh.Sf())
145  )
146  );
147 
148 
149  // Phase-1 predicted flux
150  surfaceScalarField phiHbyA1
151  (
152  IOobject::groupName("phiHbyA", phase1.name()),
153  phi1
154  );
155 
156  phiHbyA1 =
157  rAUf1
158  *(
159  (alphaRhof10 + Vmf)
160  *byDt(MRF.absolute(phi1.oldTime()))
161  + fvc::flux(U1Eqn.H())
162  + Vmf*tddtPhi2()
163  + Kdf*MRF.absolute(phi2)
164  - Ff1()
165  );
166 
167  // Phase-2 predicted flux
168  surfaceScalarField phiHbyA2
169  (
170  IOobject::groupName("phiHbyA", phase2.name()),
171  phi2
172  );
173 
174  phiHbyA2 =
175  rAUf2
176  *(
177  (alphaRhof20 + Vmf)
178  *byDt(MRF.absolute(phi2.oldTime()))
179  + fvc::flux(U2Eqn.H())
180  + Vmf*tddtPhi1()
181  + Kdf*MRF.absolute(phi1)
182  - Ff2()
183  );
184 
185 
187  (
188  "phiHbyA",
189  alphaf1*(phiHbyA1 - phig1) + alphaf2*(phiHbyA2 - phig2)
190  );
191  MRF.makeRelative(phiHbyA);
192 
193  phiHbyA1 -= phig1;
194  phiHbyA2 -= phig2;
195 
197  (
198  "rAUf",
199  mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2)
200  );
201 
202  // Update the fixedFluxPressure BCs to ensure flux consistency
204  (
205  p_rgh.boundaryFieldRef(),
206  (
207  phiHbyA.boundaryField()
208  - (
209  alphaf1.boundaryField()*phi1.boundaryField()
210  + alphaf2.boundaryField()*phi2.boundaryField()
211  )
212  )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
213  );
214 
215  tmp<fvScalarMatrix> pEqnComp1;
216  tmp<fvScalarMatrix> pEqnComp2;
217 
218  // Construct the compressibility parts of the pressure equation
219 
220  if (phase1.compressible())
221  {
222  if (pimple.transonic())
223  {
224  surfaceScalarField phid1
225  (
226  IOobject::groupName("phid", phase1.name()),
228  );
229 
230  pEqnComp1 =
231  (
232  phase1.continuityError()
234  )/rho1
235  + correction
236  (
237  (alpha1/rho1)*
238  (
239  psi1*fvm::ddt(p_rgh)
240  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
241  )
242  );
243 
244  deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
245  pEqnComp1.ref().relax();
246  }
247  else
248  {
249  pEqnComp1 =
250  (
251  phase1.continuityError()
253  )/rho1
254  + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
255  }
256  }
257  else
258  {
259  pEqnComp1 = fvm::Su(-(fvOptions(alpha1, rho1)&rho1)/rho1, p_rgh);
260  }
261 
262  if (phase2.compressible())
263  {
264  if (pimple.transonic())
265  {
266  surfaceScalarField phid2
267  (
268  IOobject::groupName("phid", phase2.name()),
270  );
271 
272  pEqnComp2 =
273  (
274  phase2.continuityError()
276  )/rho2
277  + correction
278  (
279  (alpha2/rho2)*
280  (
281  psi2*fvm::ddt(p_rgh)
282  + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
283  )
284  );
285 
286  deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
287  pEqnComp2.ref().relax();
288  }
289  else
290  {
291  pEqnComp2 =
292  (
293  phase2.continuityError()
295  )/rho2
296  + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
297  }
298  }
299  else
300  {
301  pEqnComp2 = fvm::Su(-(fvOptions(alpha2, rho2)&rho2)/rho2, p_rgh);
302  }
303 
304  if (fluid.transfersMass())
305  {
306  if (pEqnComp1.valid())
307  {
308  pEqnComp1.ref() -= fluid.dmdt()/rho1;
309  }
310  else
311  {
312  pEqnComp1 = fvm::Su(-fluid.dmdt()/rho1, p_rgh);
313  }
314 
315  if (pEqnComp2.valid())
316  {
317  pEqnComp2.ref() += fluid.dmdt()/rho2;
318  }
319  else
320  {
321  pEqnComp2 = fvm::Su(fluid.dmdt()/rho2, p_rgh);
322  }
323  }
324 
325  // Cache p prior to solve for density update
326  volScalarField p_rgh_0("p_rgh_0", p_rgh);
327 
328  while (pimple.correctNonOrthogonal())
329  {
330  fvScalarMatrix pEqnIncomp
331  (
332  fvc::div(phiHbyA)
333  - fvm::laplacian(rAUf, p_rgh)
334  );
335 
336  {
337  fvScalarMatrix pEqn(pEqnIncomp);
338 
339  if (pEqnComp1.valid())
340  {
341  pEqn += pEqnComp1();
342  }
343 
344  if (pEqnComp2.valid())
345  {
346  pEqn += pEqnComp2();
347  }
348 
349  solve
350  (
351  pEqn,
352  mesh.solver(p_rgh.select(pimple.finalInnerIter()))
353  );
354  }
355 
356  if (pimple.finalNonOrthogonalIter())
357  {
358  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
359 
360  phi = phiHbyA + pEqnIncomp.flux();
361 
362  surfaceScalarField phi1s
363  (
364  phiHbyA1
365  + alpharAUf1*mSfGradp
366  - rAUf1*Kdf*MRF.absolute(phi2)
367  );
368 
369  surfaceScalarField phi2s
370  (
371  phiHbyA2
372  + alpharAUf2*mSfGradp
373  - rAUf2*Kdf*MRF.absolute(phi1)
374  );
375 
377  (
378  ((phi2s + rAUf2*Kdf*phi1s) - (phi1s + rAUf1*Kdf*phi2s))
379  /(1.0 - rAUf1*rAUf2*sqr(Kdf))
380  );
381 
382  phi1 = phi - alphaf2*phir;
383  phi2 = phi + alphaf1*phir;
384 
385  U1 = fvc::reconstruct(MRF.absolute(phi1));
386  U1.correctBoundaryConditions();
387  fvOptions.correct(U1);
388 
389  U2 = fvc::reconstruct(MRF.absolute(phi2));
390  U2.correctBoundaryConditions();
391  fvOptions.correct(U2);
392 
393  // Set the phase dilatation rates
394  if (pEqnComp1.valid())
395  {
396  phase1.divU(-pEqnComp1 & p_rgh);
397  }
398  if (pEqnComp2.valid())
399  {
400  phase2.divU(-pEqnComp2 & p_rgh);
401  }
402  }
403  }
404 
405  // Update and limit the static pressure
406  p = max(p_rgh + rho*gh, pMin);
407 
408  // Limit p_rgh
409  p_rgh = p - rho*gh;
410 
411  // Update densities from change in p_rgh
412  rho1 += psi1*(p_rgh - p_rgh_0);
413  rho2 += psi2*(p_rgh - p_rgh_0);
414 
415  // Correct p_rgh for consistency with p and the updated densities
416  rho = fluid.rho();
417  p_rgh = p - rho*gh;
418  p_rgh.correctBoundaryConditions();
419 }
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix.
psi2
Definition: TEqns.H:35
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1))
surfaceScalarField alpharAUf2(fvc::interpolate(max(alpha2, phase2.residualAlpha()) *rAU2))
zeroField Su
Definition: alphaSuSp.H:1
surfaceScalarField & phi2
multiphaseSystem & fluid
Definition: createFields.H:11
p
Definition: pEqn.H:50
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
surfaceScalarField Df2(fvc::interpolate(D+phase2.turbulence().pPrime()))
phiHbyA
Definition: pEqn.H:20
dimensionedSymmTensor sqr(const dimensionedVector &dv)
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
surfaceScalarField & alphaPhi2
const dictionary & pimple
surfaceScalarField Df1(fvc::interpolate(D+phase1.turbulence().pPrime()))
const surfaceScalarField & alphaPhi1
surfaceScalarField & phi1
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
fv::options & fvOptions
const surfaceScalarField & ghf
surfaceScalarField Ff(fluid.Ff())
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
surfaceScalarField Kdf("Kdf", fluid.Kdf())
dynamicFvMesh & mesh
IOMRFZoneList & MRF
surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1) *mesh.magSf())
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
rhoEqn solve()
volVectorField & U1
surfaceScalarField alphaRhof10("alphaRhof10", fvc::interpolate(max(alpha1.oldTime(), phase1.residualAlpha()) *rho1.oldTime()))
alpha2
Definition: alphaEqn.H:112
phaseModel & phase1
p_rgh
Definition: pEqn.H:134
tmp< surfaceScalarField > tddtPhi2
Definition: createDDtU.H:2
const dimensionedVector & g
volScalarField & alpha1
dimensionedScalar pMin("pMin", dimPressure, fluid)
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
surfaceScalarField alphaRhof20("alphaRhof20", fvc::interpolate(max(alpha2.oldTime(), phase2.residualAlpha()) *rho2.oldTime()))
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
tmp< surfaceScalarField > tddtPhi1
Definition: createDDtU.H:1
psi1
Definition: TEqns.H:34
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho) *mesh.magSf())
surfaceScalarField alpharAUf1(fvc::interpolate(max(alpha1, phase1.residualAlpha()) *rAU1))
tmp< surfaceScalarField > Ff2
Definition: pEqn.H:55
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryFieldRef(),(phiHbyA.boundaryField() - MRF.relative(phib))/(mesh.magSf().boundaryField() *rAUf.boundaryField()))
volVectorField & U2
const volScalarField & gh
volScalarField p_rgh_0(p_rgh)
dimensioned< scalar > mag(const dimensioned< Type > &)
surfaceScalarField rAUf1(IOobject::groupName("rAUf", phase1.name()), 1.0/(byDt(alphaRhof10+Vmf)+fvc::interpolate(U1Eqn.A())+Kdf))
const dimensionedScalar & rho2
Definition: createFields.H:40
phi
Definition: pEqn.H:18
surfaceScalarField Vmf("Vmf", fluid.Vmf())
phaseModel & phase2
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:32
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
surfaceScalarField alphaf2("alphaf2", scalar(1) - alphaf1)
const dimensionedScalar & rho1
Definition: createFields.H:39
rho
Definition: pEqn.H:1
void deleteDemandDrivenData(DataPtr &dataPtr)
surfaceScalarField rAUf2(IOobject::groupName("rAUf", phase2.name()), 1.0/(byDt(alphaRhof20+Vmf)+fvc::interpolate(U2Eqn.A())+Kdf))
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
zeroField Sp
Definition: alphaSuSp.H:2
tmp< surfaceScalarField > Ff1
Definition: pEqn.H:54