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         (
  50             fvc::flux
  51             (
  52                 phi,
  53                 alpha1,
  54                 alphaScheme
  55             )
  56           + fvc::flux
  57             (
  58                 -fvc::flux(-phir, alpha2, alpharScheme),
  59                 alpha1,
  60                 alpharScheme
  61             )
  62         );
  63
  64         MULES::explicitSolve
  65         (
  66             geometricOneField(),
  67             alpha1,
  68             phi,
  69             phiAlpha1,
  70             Sp,
  71             Su,
  72             1,
  73             0
  74         );
  75
  76         surfaceScalarField rho1f(fvc::interpolate(rho1));
  77         surfaceScalarField rho2f(fvc::interpolate(rho2));
  78         rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
  79
  80         alpha2 = scalar(1) - alpha1;
  81     }
  82
  83     Info<< "Liquid phase volume fraction = "
  84         << alpha1.weightedAverage(mesh.V()).value()
  85         << "  Min(alpha1) = " << min(alpha1).value()
  86         << "  Min(alpha2) = " << min(alpha2).value()
  87         << endl;
  88 }