applications/solvers/multiphase/cavitatingFoam/pEqn.H

   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)
  37           - fvm::laplacian(rUAf, p)
  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: ************************ //