Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- {
- // surfaceScalarField betaf = scalar(1) - fvc::interpolate(sumAlpha);
- surfaceScalarField betaf("betaf", scalar(1) - fvc::interpolate(sumAlpha));
- volScalarField rUbA = 1.0/UbEqn.A();
- // surfaceScalarField rUbAf = fvc::interpolate(rUbA);
- surfaceScalarField rUbAf(fvc::interpolate(beta*rhob*rUbA));
- Ub = rUbA*UbEqn.H();
- surfaceScalarField phiDragb = rUbAf*(g & mesh.Sf());
- forAll(phasea, i)
- {
- volScalarField rUaA = 1.0/UaEqn[i].A();
- rUaAf[i] = fvc::interpolate(rUaA);
- Ua[i] = rUaA*UaEqn[i].H();
- phiDraga[i] = fvc::interpolate(beta/phasea[i].rho()*K[i]*rUaA)*phib
- + rUaAf[i]*(g & mesh.Sf());
- // phiDragb += fvc::interpolate(alpha[i]/rhob*K[i]*rUbA)*phia[i];
- phiDragb += fvc::interpolate((alpha[i]/rhob)*K[i]*rUbA)*phia[i];
- phia[i] = (fvc::interpolate(Ua[i]) & mesh.Sf()) + rUaAf[i]*fvc::ddtCorr(Ua[i], phia[i]);
- }
- // Fix for gravity on outlet boundary.
- forAll(p.boundaryField(), patchi)
- {
- if (isA<zeroGradientFvPatchScalarField>(p.boundaryField()[patchi]))
- {
- forAll(phasea, i)
- {
- phiDraga[i].boundaryField()[patchi] = 0.0;
- }
- phiDragb.boundaryField()[patchi] = 0.0;
- }
- }
- phib = (fvc::interpolate(Ub) & mesh.Sf()) + rUbAf*fvc::ddtCorr(Ub, phib)
- + phiDragb;
- phi = betaf*phib;
- forAll(phasea, i)
- {
- phia[i] += phiDraga[i];
- surfaceScalarField alphaf = fvc::interpolate(alpha[i]);
- phi += alphaf*phia[i];
- }
- surfaceScalarField Dp
- (
- "(rho*(1|A(U)))",
- betaf*rUbAf/rhob
- );
- forAll(phasea, i)
- {
- surfaceScalarField alphaf = fvc::interpolate(alpha[i]);
- Dp == Dp + alphaf*rUaAf[i]/phasea[i].rho();
- }
- for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
- {
- fvScalarMatrix pEqn
- (
- fvm::laplacian(Dp, p) == fvc::div(phi)
- );
- pEqn.setReference(pRefCell, pRefValue);
- pEqn.solve();
- if (nonOrth == nNonOrthCorr)
- {
- surfaceScalarField SfGradp = pEqn.flux()/Dp;
- phib -= rUbAf*SfGradp/rhob;
- phi = betaf*phib;
- forAll(phia, i)
- {
- surfaceScalarField alphaf = fvc::interpolate(alpha[i]);
- phia[i] -= rUaAf[i]*SfGradp/phasea[i].rho();
- phi += alphaf*phia[i];
- }
- p.relax();
- SfGradp = pEqn.flux()/Dp;
- Ub += fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob);
- Ub.correctBoundaryConditions();
- U = beta*Ub;
- forAll(Ua, i)
- {
- Ua[i] += fvc::reconstruct(phiDraga[i] - rUaAf[i]*SfGradp/phasea[i].rho());
- Ua[i].correctBoundaryConditions();
- U += alpha[i]*Ua[i];
- }
- }
- }
- }
- #include "continuityErrs.H"
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement