applications/solvers/multiphase/cavitatingFoam/pEqn.H

   1 {
   2     if (pimple.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 rAU(1.0/UEqn.A());
  15     surfaceScalarField rAUf("rAUf", rhof*fvc::interpolate(rAU));
  16     volVectorField HbyA(rAU*UEqn.H());
  17
  18     phiv = (fvc::interpolate(HbyA) & mesh.Sf())
  19          + fvc::ddtPhiCorr(rAU, rho, U, phiv);
  20
  21     p.boundaryField().updateCoeffs();
  22
  23     surfaceScalarField phiGradp(rAUf*mesh.magSf()*fvc::snGrad(p));
  24
  25     phiv -= phiGradp/rhof;
  26
  27     #include "resetPhivPatches.H"
  28
  29     for (int nonOrth=0; nonOrth<=pimple.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)
  37           - fvm::laplacian(rAUf, p)
  38         );
  39
  40         pEqn.solve
  41         (
  42             mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
  43         );
  44
  45         if (nonOrth == pimple.nNonOrthCorr())
  46         {
  47             phiv += (phiGradp + pEqn.flux())/rhof;
  48         }
  49     }
  50
  51     Info<< "Predicted p max-min : " << max(p).value()
  52         << " " << min(p).value() << endl;
  53
  54     rho == max
  55     (
  56         psi*p
  57       + (1.0 - gamma)*rhol0
  58       + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
  59         rhoMin
  60     );
  61
  62     #include "gammaPsi.H"
  63
  64     p =
  65     (
  66         rho
  67       - (1.0 - gamma)*rhol0
  68       - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
  69     )/psi;
  70
  71     p.correctBoundaryConditions();
  72
  73     Info<< "Phase-change corrected p max-min : " << max(p).value()
  74         << " " << min(p).value() << endl;
  75
  76     // Correct velocity
  77
  78     U = HbyA - rAU*fvc::grad(p);
  79
  80     // Remove the swirl component of velocity for "wedge" cases
  81     if (pimple.dict().found("removeSwirl"))
  82     {
  83         label swirlCmpt(readLabel(pimple.dict().lookup("removeSwirl")));
  84
  85         Info<< "Removing swirl component-" << swirlCmpt << " of U" << endl;
  86         U.field().replace(swirlCmpt, 0.0);
  87     }
  88
  89     U.correctBoundaryConditions();
  90
  91     Info<< "max(U) " << max(mag(U)).value() << endl;
  92 }