CFD Online Discussion Forums - solver for sub-sonic flow

Hi Foamers

I am trying to develop a subsonic solver based on original PISO for low mach number. I have started with the sonicFoam solver and changed that solver accordingly.
Changes are:

1. The energy equation is coupled with the pressure equation to solve it iteratively. (as mentioned in the original PISO)
2. The pressure equation is changed and reference is taken from other sonicfoam solver. (like XiFoam, as suggested by Henery)
3. solution of continuity equation at the start of iteration is abandoned to preserve the pressure fluctuation at low ach number which would die otherwise.

the source file is:

#include "fvCFD.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { # include "setRootCase.H" # include "createTime.H" # include "createMesh.H" # include "readThermodynamicProperties.H" # include "readTransportProperties.H" # include "createFields.H" # include "initContinuityErrs.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; for (runTime++; !runTime.end(); runTime++) { Info<< "Time = " << runTime.timeName() << nl << endl; # include "readPISOControls.H" # include "compressibleCourantNo.H" // removed solving of the continuity equation to make sure low pressure fluctuations do not die out //# include "rhoEqn.H" fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(phi, U) - fvm::laplacian(mu, U) ); solve(UEqn == -fvc::grad(p)); // --- PISO loop for (int corr=0; corr<nCorr; corr++) { // coupled energy and pressure equations !! solve ( fvm::ddt(rho, e) + fvm::div(phi, e) - fvm::laplacian(mu, e) == - p*fvc::div(phi/fvc::interpolate(rho)) + mu*magSqr(symm(fvc::grad(U))) ); T = e/Cv; psi = 1.0/(R*T); volScalarField rUA = 1.0/UEqn.A(); U = rUA*UEqn.H(); // pressure equation has been copied from other solvers solving for sub sonic flow cases for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::ddt(psi, p) + fvc::div(phi) - fvm::laplacian(rho*rUA, p) ); pEqn.solve(); phi = pEqn.flux(); } # include "compressibleContinuityErrs.H" U -= rUA*fvc::grad(p); U.correctBoundaryConditions(); } rho = psi*p; runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return(0); }

I have tested this solver with plane channel flow and result are good. I have checked it with forward step problem also.
Let me know if this solver is correct or not? or if there is any improvement you can suggest.

with regards,

Nishant Singh