FreeFOAM The Cross-Platform CFD Toolkit
pEqn.H
Go to the documentation of this file.
1 {
2  if (nOuterCorr == 1)
3  {
4  p =
5  (
6  rho
7  - (1.0 - gamma)*rhol0
8  - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
9  )/psi;
10  }
11 
12  surfaceScalarField rhof = fvc::interpolate(rho, "rhof");
13 
14  volScalarField rUA = 1.0/UEqn.A();
15  surfaceScalarField rUAf("rUAf", rhof*fvc::interpolate(rUA));
16  volVectorField HbyA = rUA*UEqn.H();
17 
18  phiv = (fvc::interpolate(HbyA) & mesh.Sf())
19  + fvc::ddtPhiCorr(rUA, rho, U, phiv);
20 
21  p.boundaryField().updateCoeffs();
22 
23  surfaceScalarField phiGradp = rUAf*mesh.magSf()*fvc::snGrad(p);
24 
25  phiv -= phiGradp/rhof;
26 
27  #include "resetPhivPatches.H"
28 
29  for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
30  {
31  fvScalarMatrix pEqn
32  (
33  fvm::ddt(psi, p)
34  - (rhol0 + (psil - psiv)*pSat)*fvc::ddt(gamma) - pSat*fvc::ddt(psi)
35  + fvc::div(phiv, rho)
36  + fvc::div(phiGradp)
38  );
39 
40  if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
41  {
42  pEqn.solve(mesh.solver(p.name() + "Final"));
43  }
44  else
45  {
46  pEqn.solve(mesh.solver(p.name()));
47  }
48 
49  if (nonOrth == nNonOrthCorr)
50  {
51  phiv += (phiGradp + pEqn.flux())/rhof;
52  }
53  }
54 
55  Info<< "Predicted p max-min : " << max(p).value()
56  << " " << min(p).value() << endl;
57 
58  rho == max
59  (
60  psi*p
61  + (1.0 - gamma)*rhol0
62  + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
63  rhoMin
64  );
65 
66  #include "gammaPsi.H"
67 
68  p =
69  (
70  rho
71  - (1.0 - gamma)*rhol0
72  - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
73  )/psi;
74 
75  p.correctBoundaryConditions();
76 
77  Info<< "Phase-change corrected p max-min : " << max(p).value()
78  << " " << min(p).value() << endl;
79 
80  // Correct velocity
81 
82  U = HbyA - rUA*fvc::grad(p);
83 
84  // Remove the swirl component of velocity for "wedge" cases
85  if (piso.found("removeSwirl"))
86  {
87  label swirlCmpt(readLabel(piso.lookup("removeSwirl")));
88 
89  Info<< "Removing swirl component-" << swirlCmpt << " of U" << endl;
90  U.field().replace(swirlCmpt, 0.0);
91  }
92 
93  U.correctBoundaryConditions();
94 
95  Info<< "max(U) " << max(mag(U)).value() << endl;
96 }
97 
98 // ************************ vim: set sw=4 sts=4 et: ************************ //