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         (
  11             fvc::flux
  12             (
  13                 phi,
  14                 alpha1,
  15                 alphaScheme
  16             )
  17           + fvc::flux
  18             (
  19                 -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
  20                 alpha1,
  21                 alpharScheme
  22             )
  23         );
  24
  25         Pair<tmp<volScalarField> > vDotAlphal =
  26             twoPhaseProperties->vDotAlphal();
  27         const volScalarField& vDotcAlphal = vDotAlphal[0]();
  28         const volScalarField& vDotvAlphal = vDotAlphal[1]();
  29
  30         volScalarField Sp
  31         (
  32             IOobject
  33             (
  34                 "Sp",
  35                 runTime.timeName(),
  36                 mesh
  37             ),
  38             vDotvAlphal - vDotcAlphal
  39         );
  40
  41         volScalarField Su
  42         (
  43             IOobject
  44             (
  45                 "Su",
  46                 runTime.timeName(),
  47                 mesh
  48             ),
  49             // Divergence term is handled explicitly to be
  50             // consistent with the explicit transport solution
  51             divU*alpha1
  52           + vDotcAlphal
  53         );
  54
  55         //MULES::explicitSolve
  56         //(
  57         //    geometricOneField(),
  58         //    alpha1,
  59         //    phi,
  60         //    phiAlpha,
  61         //    Sp,
  62         //    Su,
  63         //    1,
  64         //    0
  65         //);
  66
  67         MULES::implicitSolve
  68         (
  69             geometricOneField(),
  70             alpha1,
  71             phi,
  72             phiAlpha,
  73             Sp,
  74             Su,
  75             1,
  76             0
  77         );
  78
  79         rhoPhi +=
  80             (runTime.deltaT()/totalDeltaT)
  81            *(phiAlpha*(rho1 - rho2) + phi*rho2);
  82     }
  83
  84     Info<< "Liquid phase volume fraction = "
  85         << alpha1.weightedAverage(mesh.V()).value()
  86         << "  Min(alpha1) = " << min(alpha1).value()
  87         << "  Max(alpha1) = " << max(alpha1).value()
  88         << endl;
  89 }