applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

   1 {
   2     word alphaScheme("div(phi,alpha)");
   3     word alpharScheme("div(phirb,alpha)");
   4
   5     surfaceScalarField phir("phir", phic*interface.nHatf());
   6
   7     for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
   8     {
   9         surfaceScalarField phiAlpha =
  10             fvc::flux
  11             (
  12                 phi,
  13                 alpha1,
  14                 alphaScheme
  15             )
  16           + fvc::flux
  17             (
  18                 -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
  19                 alpha1,
  20                 alpharScheme
  21             );
  22
  23         Pair<tmp<volScalarField> > vDotAlphal =
  24             twoPhaseProperties->vDotAlphal();
  25         const volScalarField& vDotcAlphal = vDotAlphal[0]();
  26         const volScalarField& vDotvAlphal = vDotAlphal[1]();
  27
  28         volScalarField Sp
  29         (
  30             IOobject
  31             (
  32                 "Sp",
  33                 runTime.timeName(),
  34                 mesh
  35             ),
  36             vDotvAlphal - vDotcAlphal
  37         );
  38
  39         volScalarField Su
  40         (
  41             IOobject
  42             (
  43                 "Su",
  44                 runTime.timeName(),
  45                 mesh
  46             ),
  47             // Divergence term is handled explicitly to be
  48             // consistent with the explicit transport solution
  49             divU*alpha1
  50           + vDotcAlphal
  51         );
  52
  53         //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
  54         //MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
  55         MULES::implicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
  56
  57         rhoPhi +=
  58             (runTime.deltaT()/totalDeltaT)
  59            *(phiAlpha*(rho1 - rho2) + phi*rho2);
  60     }
  61
  62     Info<< "Liquid phase volume fraction = "
  63         << alpha1.weightedAverage(mesh.V()).value()
  64         << "  Min(alpha1) = " << min(alpha1).value()
  65         << "  Max(alpha1) = " << max(alpha1).value()
  66         << endl;
  67 }