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