applications/solvers/multiphase/twoPhaseEulerFoam/alphaEqn.H

   1 {
   2     word scheme("div(phi,alpha)");
   3     word schemer("div(phir,alpha)");
   4
   5     surfaceScalarField phic = phi;
   6     surfaceScalarField phir = phia - phib;
   7
   8     if (g0.value() > 0.0)
   9     {
  10         surfaceScalarField alphaf = fvc::interpolate(alpha);
  11         surfaceScalarField phipp = ppMagf*fvc::snGrad(alpha)*mesh.magSf();
  12         phir += phipp;
  13         phic += fvc::interpolate(alpha)*phipp;
  14     }
  15
  16     for (int acorr=0; acorr<nAlphaCorr; acorr++)
  17     {
  18         fvScalarMatrix alphaEqn
  19         (
  20              fvm::ddt(alpha)
  21            + fvm::div(phic, alpha, scheme)
  22            + fvm::div(-fvc::flux(-phir, beta, schemer), alpha, schemer)
  23         );
  24
  25         if (g0.value() > 0.0)
  26         {
  27             ppMagf = rUaAf*fvc::interpolate
  28             (
  29                 (1.0/(rhoa*(alpha + scalar(0.0001))))
  30                *g0*min(exp(preAlphaExp*(alpha - alphaMax)), expMax)
  31             );
  32
  33             alphaEqn -= fvm::laplacian
  34             (
  35                 (fvc::interpolate(alpha) + scalar(0.0001))*ppMagf,
  36                 alpha,
  37                 "laplacian(alphaPpMag,alpha)"
  38             );
  39         }
  40
  41         alphaEqn.relax();
  42         alphaEqn.solve();
  43
  44 #       include "packingLimiter.H"
  45
  46         beta = scalar(1) - alpha;
  47
  48         Info<< "Dispersed phase volume fraction = "
  49             << alpha.weightedAverage(mesh.V()).value()
  50             << "  Min(alpha) = " << min(alpha).value()
  51             << "  Max(alpha) = " << max(alpha).value()
  52             << endl;
  53     }
  54 }
  55
  56 rho = alpha*rhoa + beta*rhob;