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::interpolate(U1Eqn.H()) & mesh.Sf())
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::interpolate(U2Eqn.H()) & mesh.Sf())
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.boundaryField(),
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
242  (
243  psi1*fvm::ddt(p_rgh)
244  + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
245  );
246  deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
247  pEqnComp1().relax();
248  }
249 
250  if (phase2.compressible())
251  {
252  pEqnComp2 =
253  (
254  phase2.continuityError()
256  )/rho2
258  (
259  psi2*fvm::ddt(p_rgh)
260  + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
261  );
262  deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
263  pEqnComp2().relax();
264  }
265  }
266  else
267  {
268  if (phase1.compressible())
269  {
270  pEqnComp1 =
271  (
272  phase1.continuityError()
274  )/rho1
275  + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
276  }
277 
278  if (phase2.compressible())
279  {
280  pEqnComp2 =
281  (
282  phase2.continuityError()
284  )/rho2
285  + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
286  }
287  }
288 
289  if (fluid.transfersMass())
290  {
291  if (pEqnComp1.valid())
292  {
293  pEqnComp1() -= fluid.dmdt()/rho1;
294  }
295  else
296  {
297  pEqnComp1 = fvm::Su(-fluid.dmdt()/rho1, p_rgh);
298  }
299 
300  if (pEqnComp2.valid())
301  {
302  pEqnComp2() += fluid.dmdt()/rho2;
303  }
304  else
305  {
306  pEqnComp2 = fvm::Su(fluid.dmdt()/rho2, p_rgh);
307  }
308  }
309 
310  // Cache p prior to solve for density update
311  volScalarField p_rgh_0("p_rgh_0", p_rgh);
312 
313  while (pimple.correctNonOrthogonal())
314  {
315  fvScalarMatrix pEqnIncomp
316  (
317  fvc::div(phiHbyA)
318  - fvm::laplacian(rAUf, p_rgh)
319  );
320 
321  {
322  fvScalarMatrix pEqn(pEqnIncomp);
323 
324  if (pEqnComp1.valid())
325  {
326  pEqn += pEqnComp1();
327  }
328 
329  if (pEqnComp2.valid())
330  {
331  pEqn += pEqnComp2();
332  }
333 
334  solve
335  (
336  pEqn,
337  mesh.solver(p_rgh.select(pimple.finalInnerIter()))
338  );
339  }
340 
341  if (pimple.finalNonOrthogonalIter())
342  {
343  surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
344 
345  phi = phiHbyA + pEqnIncomp.flux();
346 
347  surfaceScalarField phi1s
348  (
349  phiHbyA1
350  + alpharAUf1*mSfGradp
351  - rAUf1*Kdf*MRF.absolute(phi2)
352  );
353 
354  surfaceScalarField phi2s
355  (
356  phiHbyA2
357  + alpharAUf2*mSfGradp
358  - rAUf2*Kdf*MRF.absolute(phi1)
359  );
360 
362  (
363  ((phi2s + rAUf2*Kdf*phi1s) - (phi1s + rAUf1*Kdf*phi2s))
364  /(1.0 - rAUf1*rAUf2*sqr(Kdf))
365  );
366 
367  phi1 = phi - alphaf2*phir;
368  phi2 = phi + alphaf1*phir;
369 
370  U1 = fvc::reconstruct(MRF.absolute(phi1));
371  U1.correctBoundaryConditions();
372  fvOptions.correct(U1);
373 
374  U2 = fvc::reconstruct(MRF.absolute(phi2));
375  U2.correctBoundaryConditions();
376  fvOptions.correct(U2);
377 
378  // Set the phase dilatation rates
379  if (pEqnComp1.valid())
380  {
381  phase1.divU(-pEqnComp1 & p_rgh);
382  }
383  if (pEqnComp2.valid())
384  {
385  phase2.divU(-pEqnComp2 & p_rgh);
386  }
387  }
388  }
389 
390  Info<< "min(p) = " << min(p_rgh + rho*gh).value() << endl;
391  if (min(p_rgh + rho*gh) < pMin)
392  {
393  Info<< "Clipping p" << endl;
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 }
ddtPhi2
Definition: DDtU.H:6
phaseModel & phase1
Definition: createFields.H:12
phaseModel & phase2
Definition: createFields.H:13
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix.
volVectorField & U1
Definition: createFields.H:18
multiphaseSystem & fluid
Definition: createFields.H:10
phi
Definition: pEqn.H:20
dimensioned< scalar > mag(const dimensioned< Type > &)
setSnGrad< fixedFluxPressureFvPatchScalarField >(p_rgh.boundaryField(),( phiHbyA.boundaryField() -MRF.relative(mesh.Sf().boundaryField()&U.boundaryField()) *rho.boundaryField() )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField()))
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU))
surfaceScalarField & alphaPhi1
Definition: createFields.H:20
volVectorField & U2
Definition: createFields.H:23
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> reconstruct(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
surfaceScalarField phir(phic *interface.nHatf())
surfaceScalarField alphaf2("alphaf2", scalar(1)-alphaf1)
void deleteDemandDrivenData(DataPtr &dataPtr)
surfaceScalarField rAUf2(IOobject::groupName("rAUf", phase2.name()), 1.0/( (alphaRhof20+Vmf)/runTime.deltaT() +fvc::interpolate(U2Eqn.A()) +Kdf ))
tmp< surfaceScalarField > Ff1
Definition: pEqn.H:54
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
tmp< GeometricField< Type, fvPatchField, volMesh > > laplacian(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcLaplacian.C:45
psi2
Definition: TEqns.H:35
volScalarField & p_rgh
ddtPhi1
Definition: DDtU.H:1
messageStream Info
tmp< surfaceScalarField > interpolate(const RhoType &rho)
surfaceScalarField & phi2
Definition: createFields.H:24
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
dynamicFvMesh & mesh
const dictionary & pimple
p
Definition: pEqn.H:59
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:44
rhoEqn solve()
fv::IOoptionList & fvOptions
surfaceScalarField Kdf("Kdf", fluid.Kdf())
surfaceScalarField & alphaPhi2
Definition: createFields.H:25
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf())
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
surfaceScalarField ghSnGradRho(ghf *fvc::snGrad(rho)*mesh.magSf())
IOMRFZoneList & MRF
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
surfaceScalarField alphaf1("alphaf1", fvc::interpolate(alpha1))
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
surfaceScalarField alphaRhof10("alphaRhof10", fvc::interpolate( max(alpha1.oldTime(), phase1.residualAlpha()) *rho1.oldTime() ))
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
rho1
Definition: pEqn.H:124
volScalarField p_rgh_0(p_rgh)
phiHbyA
Definition: pEqn.H:21
rho2
Definition: pEqn.H:125
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
volScalarField & alpha1
Definition: createFields.H:15
surfaceScalarField & phi1
Definition: createFields.H:19
surfaceScalarField alphaRhof20("alphaRhof20", fvc::interpolate( max(alpha2.oldTime(), phase2.residualAlpha()) *rho2.oldTime() ))
surfaceScalarField Df1(fvc::interpolate(D+phase1.turbulence().pPrime()))
alpha2
Definition: alphaEqn.H:112
tmp< surfaceScalarField > Ff2
Definition: pEqn.H:55
surfaceScalarField rAUf1(IOobject::groupName("rAUf", phase1.name()), 1.0/( (alphaRhof10+Vmf)/runTime.deltaT() +fvc::interpolate(U1Eqn.A()) +Kdf ))
const dimensionedVector & g
surfaceScalarField alpharAUf2(fvc::interpolate(max(alpha2, phase2.residualAlpha())*rAU2))
psi1
Definition: TEqns.H:34
surfaceScalarField alpharAUf1(fvc::interpolate(max(alpha1, phase1.residualAlpha())*rAU1))
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
rho
Definition: pEqn.H:1
surfaceScalarField Ff(fluid.Ff())
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedScalar pMin("pMin", dimPressure, fluid)
surfaceScalarField Vmf("Vmf", fluid.Vmf())
surfaceScalarField Df2(fvc::interpolate(D+phase2.turbulence().pPrime()))
const surfaceScalarField & ghf
const volScalarField & gh