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  (alphaRhof10 + Vmf)/runTime.deltaT()
36  + fvc::interpolate(U1Eqn.A())
37  + Kdf
38  )
39 );
40 
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
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  // Update continuity errors due to temperature changes
97  fluid.correct();
98 
99  volScalarField rho("rho", fluid.rho());
100 
101  // Correct p_rgh for consistency with p and the updated densities
102  p_rgh = p - rho*gh;
103 
106 
107  // Correct fixed-flux BCs to be consistent with the velocity BCs
108  MRF.correctBoundaryFlux(U1, phi1);
109  MRF.correctBoundaryFlux(U2, phi2);
110 
112  (
113  IOobject::groupName("alpharAUf", phase1.name()),
114  max(alphaf1, phase1.residualAlpha())*rAUf1
115  );
116 
118  (
119  IOobject::groupName("alpharAUf", phase2.name()),
120  max(alphaf2, phase2.residualAlpha())*rAUf2
121  );
122 
124  (
125  "ghSnGradRho",
126  ghf*fvc::snGrad(rho)*mesh.magSf()
127  );
128 
129  // Phase-1 buoyancy flux
130  surfaceScalarField phig1
131  (
132  IOobject::groupName("phig", phase1.name()),
133  alpharAUf1
134  *(
136  - alphaf2*(rhof1 - rhof2)*(g & mesh.Sf())
137  )
138  );
139 
140  // Phase-2 buoyancy flux
141  surfaceScalarField phig2
142  (
143  IOobject::groupName("phig", phase2.name()),
144  alpharAUf2
145  *(
147  - alphaf1*(rhof2 - rhof1)*(g & mesh.Sf())
148  )
149  );
150 
151 
152  // Phase-1 predicted flux
153  surfaceScalarField phiHbyA1
154  (
155  IOobject::groupName("phiHbyA", phase1.name()),
156  phi1
157  );
158 
159  phiHbyA1 =
160  rAUf1
161  *(
162  (alphaRhof10 + Vmf)
163  *MRF.absolute(phi1.oldTime())/runTime.deltaT()
164  + fvc::flux(U1Eqn.H())
165  + Vmf*ddtPhi2
166  + Kdf*MRF.absolute(phi2)
167  - Ff1()
168  );
169 
170  // Phase-2 predicted flux
171  surfaceScalarField phiHbyA2
172  (
173  IOobject::groupName("phiHbyA", phase2.name()),
174  phi2
175  );
176 
177  phiHbyA2 =
178  rAUf2
179  *(
180  (alphaRhof20 + Vmf)
181  *MRF.absolute(phi2.oldTime())/runTime.deltaT()
182  + fvc::flux(U2Eqn.H())
183  + Vmf*ddtPhi1
184  + Kdf*MRF.absolute(phi1)
185  - Ff2()
186  );
187 
188 
190  (
191  "phiHbyA",
192  alphaf1*(phiHbyA1 - phig1) + alphaf2*(phiHbyA2 - phig2)
193  );
194  MRF.makeRelative(phiHbyA);
195 
196  phiHbyA1 -= phig1;
197  phiHbyA2 -= phig2;
198 
200  (
201  "rAUf",
202  mag(alphaf1*alpharAUf1 + alphaf2*alpharAUf2)
203  );
204 
205  // Update the fixedFluxPressure BCs to ensure flux consistency
207  (
208  p_rgh.boundaryFieldRef(),
209  (
210  phiHbyA.boundaryField()
211  - (
212  alphaf1.boundaryField()*phi1.boundaryField()
213  + alphaf2.boundaryField()*phi2.boundaryField()
214  )
215  )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
216  );
217 
218  tmp<fvScalarMatrix> pEqnComp1;
219  tmp<fvScalarMatrix> pEqnComp2;
220 
221  if (pimple.transonic())
222  {
223  surfaceScalarField phid1
224  (
225  IOobject::groupName("phid", phase1.name()),
227  );
228  surfaceScalarField phid2
229  (
230  IOobject::groupName("phid", phase2.name()),
232  );
233 
234  if (phase1.compressible())
235  {
236  pEqnComp1 =
237  (
238  phase1.continuityError()
240  )/rho1
241  + correction
242  (
243  (alpha1/rho1)*
244  (
245  psi1*fvm::ddt(p_rgh)
246  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
247  )
248  );
249  deleteDemandDrivenData(pEqnComp1.ref().faceFluxCorrectionPtr());
250  pEqnComp1.ref().relax();
251  }
252 
253  if (phase2.compressible())
254  {
255  pEqnComp2 =
256  (
257  phase2.continuityError()
259  )/rho2
260  + correction
261  (
262  (alpha2/rho2)*
263  (
264  psi2*fvm::ddt(p_rgh)
265  + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
266  )
267  );
268  deleteDemandDrivenData(pEqnComp2.ref().faceFluxCorrectionPtr());
269  pEqnComp2.ref().relax();
270  }
271  }
272  else
273  {
274  if (phase1.compressible())
275  {
276  pEqnComp1 =
277  (
278  phase1.continuityError()
280  )/rho1
281  + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
282  }
283 
284  if (phase2.compressible())
285  {
286  pEqnComp2 =
287  (
288  phase2.continuityError()
290  )/rho2
291  + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
292  }
293  }
294 
295  if (fluid.transfersMass())
296  {
297  if (pEqnComp1.valid())
298  {
299  pEqnComp1.ref() -= fluid.dmdt()/rho1;
300  }
301  else
302  {
303  pEqnComp1 = fvm::Su(-fluid.dmdt()/rho1, p_rgh);
304  }
305 
306  if (pEqnComp2.valid())
307  {
308  pEqnComp2.ref() += fluid.dmdt()/rho2;
309  }
310  else
311  {
312  pEqnComp2 = fvm::Su(fluid.dmdt()/rho2, p_rgh);
313  }
314  }
315 
316  // Cache p prior to solve for density update
317  volScalarField p_rgh_0("p_rgh_0", p_rgh);
318 
319  while (pimple.correctNonOrthogonal())
320  {
321  fvScalarMatrix pEqnIncomp
322  (
323  fvc::div(phiHbyA)
324  - fvm::laplacian(rAUf, p_rgh)
325  );
326 
327  {
328  fvScalarMatrix pEqn(pEqnIncomp);
329 
330  if (pEqnComp1.valid())
331  {
332  pEqn += pEqnComp1();
333  }
334 
335  if (pEqnComp2.valid())
336  {
337  pEqn += pEqnComp2();
338  }
339 
340  solve
341  (
342  pEqn,
343  mesh.solver(p_rgh.select(pimple.finalInnerIter()))
344  );
345  }
346 
347  if (pimple.finalNonOrthogonalIter())
348  {
349  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
350 
351  phi = phiHbyA + pEqnIncomp.flux();
352 
353  surfaceScalarField phi1s
354  (
355  phiHbyA1
356  + alpharAUf1*mSfGradp
357  - rAUf1*Kdf*MRF.absolute(phi2)
358  );
359 
360  surfaceScalarField phi2s
361  (
362  phiHbyA2
363  + alpharAUf2*mSfGradp
364  - rAUf2*Kdf*MRF.absolute(phi1)
365  );
366 
368  (
369  ((phi2s + rAUf2*Kdf*phi1s) - (phi1s + rAUf1*Kdf*phi2s))
370  /(1.0 - rAUf1*rAUf2*sqr(Kdf))
371  );
372 
373  phi1 = phi - alphaf2*phir;
374  phi2 = phi + alphaf1*phir;
375 
376  U1 = fvc::reconstruct(MRF.absolute(phi1));
377  U1.correctBoundaryConditions();
378  fvOptions.correct(U1);
379 
380  U2 = fvc::reconstruct(MRF.absolute(phi2));
381  U2.correctBoundaryConditions();
382  fvOptions.correct(U2);
383 
384  // Set the phase dilatation rates
385  if (pEqnComp1.valid())
386  {
387  phase1.divU(-pEqnComp1 & p_rgh);
388  }
389  if (pEqnComp2.valid())
390  {
391  phase2.divU(-pEqnComp2 & p_rgh);
392  }
393  }
394  }
395 
396  // Update and limit the static pressure
397  p = max(p_rgh + rho*gh, pMin);
398 
399  // Limit p_rgh
400  p_rgh = p - rho*gh;
401 
402  // Update densities from change in p_rgh
403  rho1 += psi1*(p_rgh - p_rgh_0);
404  rho2 += psi2*(p_rgh - p_rgh_0);
405 
406  // Correct p_rgh for consistency with p and the updated densities
407  rho = fluid.rho();
408  p_rgh = p - rho*gh;
409  p_rgh.correctBoundaryConditions();
410 }
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 & phi2
multiphaseSystem & fluid
Definition: createFields.H:10
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
surfaceScalarField alphaRhof10("alphaRhof10", fvc::interpolate(max(alpha1.oldTime(), phase1.residualAlpha())*rho1.oldTime()))
const dictionary & pimple
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
surfaceScalarField Df1(fvc::interpolate(D+phase1.turbulence().pPrime()))
surfaceScalarField & phi1
ddtPhi2
Definition: DDtU.H:6
tmp< surfaceScalarField > interpolate(const RhoType &rho)
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
fv::options & fvOptions
surfaceScalarField phir(phic *interface.nHatf())
rho1
Definition: pEqn.H:114
const surfaceScalarField & ghf
surfaceScalarField Ff(fluid.Ff())
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf())
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
surfaceScalarField Kdf("Kdf", fluid.Kdf())
dynamicFvMesh & mesh
rhoEqn solve()
IOMRFZoneList & MRF
rho2
Definition: pEqn.H:115
ddtPhi1
Definition: DDtU.H:1
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho)*mesh.magSf())
surfaceScalarField alpharAUf1(fvc::interpolate(max(alpha1, phase1.residualAlpha())*rAU1))
volVectorField & U1
surfaceScalarField alpharAUf2(fvc::interpolate(max(alpha2, phase2.residualAlpha())*rAU2))
alpha2
Definition: alphaEqn.H:112
phaseModel & phase1
surfaceScalarField rAUf1(IOobject::groupName("rAUf", phase1.name()), 1.0/((alphaRhof10+Vmf)/runTime.deltaT()+fvc::interpolate(U1Eqn.A())+Kdf))
p_rgh
Definition: pEqn.H:120
const dimensionedVector & g
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryFieldRef(),(phiHbyA.boundaryField()-MRF.relative(phib))/(mesh.magSf().boundaryField()*rAUf.boundaryField()))
volScalarField & alpha1
dimensionedScalar pMin("pMin", dimPressure, fluid)
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
psi1
Definition: TEqns.H:34
surfaceScalarField rAUf2(IOobject::groupName("rAUf", phase2.name()), 1.0/((alphaRhof20+Vmf)/runTime.deltaT()+fvc::interpolate(U2Eqn.A())+Kdf))
tmp< surfaceScalarField > Ff2
Definition: pEqn.H:55
volVectorField & U2
const volScalarField & gh
volScalarField p_rgh_0(p_rgh)
dimensioned< scalar > mag(const dimensioned< Type > &)
phi
Definition: pEqn.H:18
surfaceScalarField Vmf("Vmf", fluid.Vmf())
surfaceScalarField alphaf2("alphaf2", scalar(1)-alphaf1)
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 alphaRhof20("alphaRhof20", fvc::interpolate(max(alpha2.oldTime(), phase2.residualAlpha())*rho2.oldTime()))
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:44
rho
Definition: pEqn.H:1
void deleteDemandDrivenData(DataPtr &dataPtr)
surfaceScalarField & alphaPhi1
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
tmp< surfaceScalarField > Ff1
Definition: pEqn.H:54