Advertisement
Guest User

pEqn.H

a guest
Sep 18th, 2014
389
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 3.10 KB | None | 0 0
  1. {
  2. //    surfaceScalarField betaf = scalar(1) - fvc::interpolate(sumAlpha);
  3.     surfaceScalarField betaf("betaf", scalar(1) - fvc::interpolate(sumAlpha));
  4.  
  5.     volScalarField rUbA = 1.0/UbEqn.A();
  6.    
  7. //    surfaceScalarField rUbAf = fvc::interpolate(rUbA);
  8.     surfaceScalarField rUbAf(fvc::interpolate(beta*rhob*rUbA));
  9.  
  10.     Ub = rUbA*UbEqn.H();
  11.  
  12.         surfaceScalarField phiDragb  = rUbAf*(g & mesh.Sf());      
  13.        
  14.     forAll(phasea, i)
  15.     {
  16.         volScalarField rUaA = 1.0/UaEqn[i].A();
  17.         rUaAf[i] = fvc::interpolate(rUaA);
  18.        
  19.         Ua[i] = rUaA*UaEqn[i].H();
  20.        
  21.         phiDraga[i] = fvc::interpolate(beta/phasea[i].rho()*K[i]*rUaA)*phib
  22.                     + rUaAf[i]*(g & mesh.Sf());
  23.    
  24. //        phiDragb   += fvc::interpolate(alpha[i]/rhob*K[i]*rUbA)*phia[i];
  25.         phiDragb   += fvc::interpolate((alpha[i]/rhob)*K[i]*rUbA)*phia[i];
  26.  
  27.         phia[i] = (fvc::interpolate(Ua[i]) & mesh.Sf()) + rUaAf[i]*fvc::ddtCorr(Ua[i], phia[i]);
  28.  
  29.     }
  30.  
  31.     // Fix for gravity on outlet boundary.
  32.     forAll(p.boundaryField(), patchi)
  33.     {
  34.         if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
  35.         {
  36.             forAll(phasea, i)
  37.             {
  38.                 phiDraga[i].boundaryField()[patchi] = 0.0;
  39.             }
  40.             phiDragb.boundaryField()[patchi] = 0.0;
  41.         }
  42.     }
  43.  
  44.     phib = (fvc::interpolate(Ub) & mesh.Sf()) + rUbAf*fvc::ddtCorr(Ub, phib)
  45.          + phiDragb;
  46.     phi  = betaf*phib;
  47.  
  48.     forAll(phasea, i)
  49.     {
  50.         phia[i] += phiDraga[i];
  51.        
  52.         surfaceScalarField alphaf = fvc::interpolate(alpha[i]);
  53.         phi     += alphaf*phia[i];
  54.     }
  55.  
  56.     surfaceScalarField Dp
  57.                       (
  58.                        "(rho*(1|A(U)))",
  59.                        betaf*rUbAf/rhob
  60.                       );
  61.  
  62.     forAll(phasea, i)
  63.     {
  64.         surfaceScalarField alphaf = fvc::interpolate(alpha[i]);
  65.         Dp == Dp + alphaf*rUaAf[i]/phasea[i].rho();
  66.     }
  67.  
  68.     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  69.     {
  70.         fvScalarMatrix pEqn
  71.         (
  72.             fvm::laplacian(Dp, p) == fvc::div(phi)
  73.         );
  74.  
  75.         pEqn.setReference(pRefCell, pRefValue);
  76.         pEqn.solve();
  77.  
  78.         if (nonOrth == nNonOrthCorr)
  79.         {
  80.             surfaceScalarField SfGradp = pEqn.flux()/Dp;
  81.  
  82.             phib  -= rUbAf*SfGradp/rhob;
  83.             phi    = betaf*phib;
  84.  
  85.             forAll(phia, i)
  86.             {
  87.                 surfaceScalarField alphaf = fvc::interpolate(alpha[i]);
  88.                 phia[i] -= rUaAf[i]*SfGradp/phasea[i].rho();
  89.                 phi     += alphaf*phia[i];
  90.             }
  91.  
  92.             p.relax();
  93.             SfGradp = pEqn.flux()/Dp;
  94.  
  95.             Ub += fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob);
  96.             Ub.correctBoundaryConditions();
  97.             U   = beta*Ub;
  98.  
  99.             forAll(Ua, i)
  100.             {
  101.                 Ua[i] += fvc::reconstruct(phiDraga[i] - rUaAf[i]*SfGradp/phasea[i].rho());
  102.                 Ua[i].correctBoundaryConditions();
  103.                 U    += alpha[i]*Ua[i];
  104.             }
  105.  
  106.         }
  107.     }
  108.  
  109. }
  110.  
  111. #include "continuityErrs.H"
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement