applications/solvers/multiphase/porousInterFoam/alphaEqn.H

   1 {
   2     // Creates the porosity field for MULES
   3     volScalarField porosity
   4     (
   5         IOobject
   6         (
   7             "porosity",
   8             runTime.timeName(),
   9             mesh,
  10             IOobject::NO_READ,
  11             IOobject::NO_WRITE
  12         ),
  13         mesh,
  14         dimensionedScalar("NULL", dimless, 1.0),
  15         "zeroGradient"
  16     );
  17
  18     forAll (pZones, zoneI)
  19     {
  20         const label zoneId = pZones[zoneI].zoneId();
  21
  22         const labelList& cells= mesh.cellZones()[zoneId];
  23
  24         const scalar& zonePorosity = pZones[zoneI].porosity();
  25
  26         forAll (cells, cellI)
  27         {
  28             porosity[cells[cellI]] = zonePorosity;
  29         }
  30     }
  31
  32     // Ordinary MULES except for the argument list to explicitSolve
  33     word alphaScheme("div(phi,alpha)");
  34     word alpharScheme("div(phirb,alpha)");
  35
  36     surfaceScalarField phic = mag(phi/mesh.magSf());
  37     phic = min(interface.cAlpha()*phic, max(phic));
  38     surfaceScalarField phir = phic*interface.nHatf();
  39
  40     for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
  41     {
  42         surfaceScalarField phiAlpha =
  43             fvc::flux
  44             (
  45                 phi,
  46                 alpha1,
  47                 alphaScheme
  48             )
  49           + fvc::flux
  50             (
  51                 -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
  52                 alpha1,
  53                 alpharScheme
  54             );
  55
  56         MULES::explicitSolve
  57         (
  58             porosity,
  59             alpha1,
  60             phi,
  61             phiAlpha,
  62             zeroField(),
  63             zeroField(),
  64             1,
  65             0
  66         );
  67
  68         rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
  69     }
  70
  71     // The weightedAverage will change even in a closed box,
  72     // however, the volume of water should
  73     // remain constant - hence both values are written in the log
  74     Info<< "Liquid phase volume fraction = "
  75         << alpha1.weightedAverage(mesh.V()).value()
  76         << "  Volume of water = "
  77         << gSum( alpha1.internalField()*porosity.internalField()*mesh.V()
  78         << "  Min(alpha1) = " << min(alpha1).value()
  79         << "  Max(alpha1) = " << max(alpha1).value()
  80         << endl;
  81 }