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