applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H

   1 {
   2     word alphaScheme("div(phi,alpha)");
   3     word alpharScheme("div(phirb,alpha)");
   4
   5     surfaceScalarField phir = phic*interface.nHatf();
   6
   7     for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
   8     {
   9         volScalarField::DimensionedInternalField Sp
  10         (
  11             IOobject
  12             (
  13                 "Sp",
  14                 runTime.timeName(),
  15                 mesh
  16             ),
  17             mesh,
  18             dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
  19         );
  20
  21         volScalarField::DimensionedInternalField Su
  22         (
  23             IOobject
  24             (
  25                 "Su",
  26                 runTime.timeName(),
  27                 mesh
  28             ),
  29             // Divergence term is handled explicitly to be
  30             // consistent with the explicit transport solution
  31             divU*min(alpha1, scalar(1))
  32         );
  33
  34         forAll(dgdt, celli)
  35         {
  36             if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
  37             {
  38                 Sp[celli] -= dgdt[celli]*alpha1[celli];
  39                 Su[celli] += dgdt[celli]*alpha1[celli];
  40             }
  41             else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
  42             {
  43                 Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
  44             }
  45         }
  46
  47
  48         surfaceScalarField phiAlpha1 =
  49             fvc::flux
  50             (
  51                 phi,
  52                 alpha1,
  53                 alphaScheme
  54             )
  55           + fvc::flux
  56             (
  57                 -fvc::flux(-phir, alpha2, alpharScheme),
  58                 alpha1,
  59                 alpharScheme
  60             );
  61
  62         MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
  63
  64         surfaceScalarField rho1f = fvc::interpolate(rho1);
  65         surfaceScalarField rho2f = fvc::interpolate(rho2);
  66         rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
  67
  68         alpha2 = scalar(1) - alpha1;
  69     }
  70
  71     Info<< "Liquid phase volume fraction = "
  72         << alpha1.weightedAverage(mesh.V()).value()
  73         << "  Min(alpha1) = " << min(alpha1).value()
  74         << "  Min(alpha2) = " << min(alpha2).value()
  75         << endl;
  76 }