applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

   1 {
   2     surfaceScalarField alphaf(fvc::interpolate(alpha));
   3     surfaceScalarField betaf(scalar(1) - alphaf);
   4
   5     volScalarField rUaA(1.0/UaEqn.A());
   6     volScalarField rUbA(1.0/UbEqn.A());
   7
   8     phia == (fvc::interpolate(Ua) & mesh.Sf());
   9     phib == (fvc::interpolate(Ub) & mesh.Sf());
  10
  11     rUaAf = fvc::interpolate(rUaA);
  12     surfaceScalarField rUbAf(fvc::interpolate(rUbA));
  13
  14     Ua = rUaA*UaEqn.H();
  15     Ub = rUbA*UbEqn.H();
  16
  17     surfaceScalarField phiDraga
  18     (
  19         fvc::interpolate(beta/rhoa*K*rUaA)*phib + rUaAf*(g & mesh.Sf())
  20     );
  21
  22     if (g0.value() > 0.0)
  23     {
  24         phiDraga -= ppMagf*fvc::snGrad(alpha)*mesh.magSf();
  25     }
  26
  27     if (kineticTheory.on())
  28     {
  29         phiDraga -= rUaAf*fvc::snGrad(kineticTheory.pa()/rhoa)*mesh.magSf();
  30     }
  31
  32     surfaceScalarField phiDragb
  33     (
  34         fvc::interpolate(alpha/rhob*K*rUbA)*phia + rUbAf*(g & mesh.Sf())
  35     );
  36
  37     // Fix for gravity on outlet boundary.
  38     forAll(p.boundaryField(), patchi)
  39     {
  40         if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
  41         {
  42             phiDraga.boundaryField()[patchi] = 0.0;
  43             phiDragb.boundaryField()[patchi] = 0.0;
  44         }
  45     }
  46
  47     phia = (fvc::interpolate(Ua) & mesh.Sf()) + fvc::ddtPhiCorr(rUaA, Ua, phia)
  48          + phiDraga;
  49
  50     phib = (fvc::interpolate(Ub) & mesh.Sf()) + fvc::ddtPhiCorr(rUbA, Ub, phib)
  51          + phiDragb;
  52
  53     phi = alphaf*phia + betaf*phib;
  54
  55     surfaceScalarField Dp
  56     (
  57         "(rho*(1|A(U)))",
  58         alphaf*rUaAf/rhoa + betaf*rUbAf/rhob
  59     );
  60
  61     for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
  62     {
  63         fvScalarMatrix pEqn
  64         (
  65             fvm::laplacian(Dp, p) == fvc::div(phi)
  66         );
  67
  68         pEqn.setReference(pRefCell, pRefValue);
  69
  70         pEqn.solve
  71         (
  72             mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
  73         );
  74
  75         if (nonOrth == pimple.nNonOrthCorr())
  76         {
  77             surfaceScalarField SfGradp(pEqn.flux()/Dp);
  78
  79             phia -= rUaAf*SfGradp/rhoa;
  80             phib -= rUbAf*SfGradp/rhob;
  81             phi = alphaf*phia + betaf*phib;
  82
  83             p.relax();
  84             SfGradp = pEqn.flux()/Dp;
  85
  86             Ua += fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa);
  87             Ua.correctBoundaryConditions();
  88
  89             Ub += fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob);
  90             Ub.correctBoundaryConditions();
  91
  92             U = alpha*Ua + beta*Ub;
  93         }
  94     }
  95 }
  96
  97 #include "continuityErrs.H"