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             false   // Use slices for the couples
  39         );
  40
  41
  42         // Create the complete convection flux for alpha1
  43         surfaceScalarField phiAlpha1
  44         (
  45             fvc::flux
  46             (
  47                 phi,
  48                 alpha1,
  49                 alphaScheme
  50             )
  51           + fvc::flux
  52             (
  53                 -fvc::flux(-phir, alpha2, alpharScheme),
  54                 alpha1,
  55                 alpharScheme
  56             )
  57           + fvc::flux
  58             (
  59                 -fvc::flux(-phir, alpha3, alpharScheme),
  60                 alpha1,
  61                 alpharScheme
  62             )
  63         );
  64
  65         // Create the bounded (upwind) flux for alpha1
  66         surfaceScalarField phiAlpha1BD
  67         (
  68             upwind<scalar>(mesh, phi).flux(alpha1)
  69         );
  70
  71         // Calculate the flux correction for alpha1
  72         phiAlpha1 -= phiAlpha1BD;
  73
  74         // Calculate the limiter for alpha1
  75         MULES::limiter
  76         (
  77             allLambda,
  78             geometricOneField(),
  79             alpha1,
  80             phiAlpha1BD,
  81             phiAlpha1,
  82             zeroField(),
  83             zeroField(),
  84             1,
  85             0,
  86             3
  87         );
  88
  89         // Create the complete flux for alpha2
  90         surfaceScalarField phiAlpha2
  91         (
  92             fvc::flux
  93             (
  94                 phi,
  95                 alpha2,
  96                 alphaScheme
  97             )
  98           + fvc::flux
  99             (
 100                 -fvc::flux(phir, alpha1, alpharScheme),
 101                 alpha2,
 102                 alpharScheme
 103             )
 104         );
 105
 106         // Create the bounded (upwind) flux for alpha2
 107         surfaceScalarField phiAlpha2BD
 108         (
 109             upwind<scalar>(mesh, phi).flux(alpha2)
 110         );
 111
 112         // Calculate the flux correction for alpha2
 113         phiAlpha2 -= phiAlpha2BD;
 114
 115         // Further limit the limiter for alpha2
 116         MULES::limiter
 117         (
 118             allLambda,
 119             geometricOneField(),
 120             alpha2,
 121             phiAlpha2BD,
 122             phiAlpha2,
 123             zeroField(),
 124             zeroField(),
 125             1,
 126             0,
 127             3
 128         );
 129
 130         // Construct the limited fluxes
 131         phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
 132         phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
 133
 134         // Solve for alpha1
 135         solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
 136
 137         // Create the diffusion coefficients for alpha2<->alpha3
 138         volScalarField Dc23(D23*max(alpha3, scalar(0))*pos(alpha2));
 139         volScalarField Dc32(D23*max(alpha2, scalar(0))*pos(alpha3));
 140
 141         // Add the diffusive flux for alpha3->alpha2
 142         phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
 143
 144         // Solve for alpha2
 145         fvScalarMatrix alpha2Eqn
 146         (
 147             fvm::ddt(alpha2)
 148           + fvc::div(phiAlpha2)
 149           - fvm::laplacian(Dc23 + Dc32, alpha2)
 150         );
 151         alpha2Eqn.solve();
 152
 153         // Construct the complete mass flux
 154         rhoPhi =
 155               phiAlpha1*(rho1 - rho3)
 156             + (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
 157             + phi*rho3;
 158
 159         alpha3 = 1.0 - alpha1 - alpha2;
 160     }
 161
 162     Info<< "Air phase volume fraction = "
 163         << alpha1.weightedAverage(mesh.V()).value()
 164         << "  Min(alpha1) = " << min(alpha1).value()
 165         << "  Max(alpha1) = " << max(alpha1).value()
 166         << endl;
 167
 168     Info<< "Liquid phase volume fraction = "
 169         << alpha2.weightedAverage(mesh.V()).value()
 170         << "  Min(alpha2) = " << min(alpha2).value()
 171         << "  Max(alpha2) = " << max(alpha2).value()
 172         << endl;
 173 }