applications/solvers/multiphase/interMixingFoam/alphaEqns.H

   1 {
   2     word alphaScheme("div(phi,alpha)");
   3     word alpharScheme("div(phirb,alpha)");
   4
   5     surfaceScalarField phir
   6     (
   7         IOobject
   8         (
   9             "phir",
  10             runTime.timeName(),
  11             mesh,
  12             IOobject::NO_READ,
  13             IOobject::NO_WRITE
  14         ),
  15         interface.cAlpha()*mag(phi/mesh.magSf())*interface.nHatf()
  16     );
  17
  18     for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
  19     {
  20         // Create the limiter to be used for all phase-fractions
  21         scalarField allLambda(mesh.nFaces(), 1.0);
  22
  23         // Split the limiter into a surfaceScalarField
  24         slicedSurfaceScalarField lambda
  25         (
  26             IOobject
  27             (
  28                 "lambda",
  29                 mesh.time().timeName(),
  30                 mesh,
  31                 IOobject::NO_READ,
  32                 IOobject::NO_WRITE,
  33                 false
  34             ),
  35             mesh,
  36             dimless,
  37             allLambda
  38         );
  39
  40
  41         // Create the complete convection flux for alpha1
  42         surfaceScalarField phiAlpha1 =
  43             fvc::flux
  44             (
  45                 phi,
  46                 alpha1,
  47                 alphaScheme
  48             )
  49           + fvc::flux
  50             (
  51                 -fvc::flux(-phir, alpha2, alpharScheme),
  52                 alpha1,
  53                 alpharScheme
  54             )
  55           + fvc::flux
  56             (
  57                 -fvc::flux(-phir, alpha3, alpharScheme),
  58                 alpha1,
  59                 alpharScheme
  60             );
  61
  62         // Create the bounded (upwind) flux for alpha1
  63         surfaceScalarField phiAlpha1BD =
  64             upwind<scalar>(mesh, phi).flux(alpha1);
  65
  66         // Calculate the flux correction for alpha1
  67         phiAlpha1 -= phiAlpha1BD;
  68
  69         // Calculate the limiter for alpha1
  70         MULES::limiter
  71         (
  72             allLambda,
  73             geometricOneField(),
  74             alpha1,
  75             phiAlpha1BD,
  76             phiAlpha1,
  77             zeroField(),
  78             zeroField(),
  79             1,
  80             0,
  81             3
  82         );
  83
  84         // Create the complete flux for alpha2
  85         surfaceScalarField phiAlpha2 =
  86             fvc::flux
  87             (
  88                 phi,
  89                 alpha2,
  90                 alphaScheme
  91             )
  92           + fvc::flux
  93             (
  94                 -fvc::flux(phir, alpha1, alpharScheme),
  95                 alpha2,
  96                 alpharScheme
  97             );
  98
  99         // Create the bounded (upwind) flux for alpha2
 100         surfaceScalarField phiAlpha2BD =
 101             upwind<scalar>(mesh, phi).flux(alpha2);
 102
 103         // Calculate the flux correction for alpha2
 104         phiAlpha2 -= phiAlpha2BD;
 105
 106         // Further limit the limiter for alpha2
 107         MULES::limiter
 108         (
 109             allLambda,
 110             geometricOneField(),
 111             alpha2,
 112             phiAlpha2BD,
 113             phiAlpha2,
 114             zeroField(),
 115             zeroField(),
 116             1,
 117             0,
 118             3
 119         );
 120
 121         // Construct the limited fluxes
 122         phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
 123         phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
 124
 125         // Solve for alpha1
 126         solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
 127
 128         // Create the diffusion coefficients for alpha2<->alpha3
 129         volScalarField Dc23 = D23*max(alpha3, scalar(0))*pos(alpha2);
 130         volScalarField Dc32 = D23*max(alpha2, scalar(0))*pos(alpha3);
 131
 132         // Add the diffusive flux for alpha3->alpha2
 133         phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
 134
 135         // Solve for alpha2
 136         fvScalarMatrix alpha2Eqn
 137         (
 138             fvm::ddt(alpha2)
 139           + fvc::div(phiAlpha2)
 140           - fvm::laplacian(Dc23 + Dc32, alpha2)
 141         );
 142         alpha2Eqn.solve();
 143
 144         // Construct the complete mass flux
 145         rhoPhi =
 146               phiAlpha1*(rho1 - rho3)
 147             + (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
 148             + phi*rho3;
 149
 150         alpha3 = 1.0 - alpha1 - alpha2;
 151     }
 152
 153     Info<< "Air phase volume fraction = "
 154         << alpha1.weightedAverage(mesh.V()).value()
 155         << "  Min(alpha1) = " << min(alpha1).value()
 156         << "  Max(alpha1) = " << max(alpha1).value()
 157         << endl;
 158
 159     Info<< "Liquid phase volume fraction = "
 160         << alpha2.weightedAverage(mesh.V()).value()
 161         << "  Min(alpha2) = " << min(alpha2).value()
 162         << "  Max(alpha2) = " << max(alpha2).value()
 163         << endl;
 164 }