if (!pimple.SIMPLErho()) { rho = thermo.rho(); } // Thermodynamic density needs to be updated by psi*d(p) after the // pressure solution const volScalarField psip0(psi*p); volScalarField rAU(1.0/UEqn.A()); surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); surfaceScalarField phiHbyA ( "phiHbyA", ( fvc::flux(rho*HbyA) + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi)) ) + phig ); fvc::makeRelative(phiHbyA, rho, U); MRF.makeRelative(fvc::interpolate(rho), phiHbyA); // Update the pressure BCs to ensure flux consistency constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF); fvScalarMatrix p_rghDDtEqn ( fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) + fvc::div(phiHbyA) == parcels.Srho() + surfaceFilm.Srho() + fvOptions(psi, p_rgh, rho.name()) ); while (pimple.correctNonOrthogonal()) { fvScalarMatrix p_rghEqn ( p_rghDDtEqn - fvm::laplacian(rhorAUf, p_rgh) ); p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter()))); if (pimple.finalNonOrthogonalIter()) { phi = phiHbyA + p_rghEqn.flux(); // Explicitly relax pressure for momentum corrector p_rgh.relax(); U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf); U.correctBoundaryConditions(); fvOptions.correct(U); K = 0.5*magSqr(U); } } p = p_rgh + rho*gh; // Thermodynamic density update thermo.correctRho(psi*p - psip0); #include "rhoEqn.H" #include "compressibleContinuityErrs.H" if (pressureControl.limit(p)) { p.correctBoundaryConditions(); rho = thermo.rho(); p_rgh = p - rho*gh; } else if (pimple.SIMPLErho()) { rho = thermo.rho(); } // Correct rhoUf if the mesh is moving fvc::correctRhoUf(rhoUf, rho, U, phi); if (thermo.dpdt()) { dpdt = fvc::ddt(p); if (mesh.moving()) { dpdt -= fvc::div(fvc::meshPhi(rho, U), p); } }