From 98f7f0975ad997f803eb8136fef7438242330849 Mon Sep 17 00:00:00 2001 From: Henrik Rusche Date: Wed, 11 Nov 2015 16:14:33 +0100 Subject: [PATCH] Updated steadyCompressible*Foam solvers & steadyCompressibleFoam tutorials --- .../steadyCompressibleFoam/createFields.H | 3 +- .../steadyCompressibleFoam/hEqn.H | 1 + .../steadyCompressibleFoam/pEqn.H | 29 ++++++++--------- .../steadyCompressibleMRFFoam/pEqn.H | 16 ++++++---- .../steadyCompressibleSRFFoam/createFields.H | 3 +- .../steadyCompressibleSRFFoam/iEqn.H | 3 +- .../steadyCompressibleSRFFoam/pEqn.H | 32 +++++++++---------- .../2bump/system/fvSolution | 8 ++--- .../transonicBump/system/fvSchemes | 2 +- .../transonicBump/system/fvSolution | 2 -- 10 files changed, 50 insertions(+), 49 deletions(-) diff --git a/applications/solvers/compressible/steadyCompressibleFoam/createFields.H b/applications/solvers/compressible/steadyCompressibleFoam/createFields.H index e1d363a75..9bf40b3a6 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/createFields.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/createFields.H @@ -9,7 +9,8 @@ volScalarField& p = thermo.p(); volScalarField& h = thermo.h(); const volScalarField& T = thermo.T(); - const volScalarField& psi = thermo.psi(); + volScalarField psis("psi", thermo.psi()/thermo.Cp()*thermo.Cv()); + psis.oldTime(); volScalarField rho ( diff --git a/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H index 780c9401d..91a8e1dc3 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H @@ -34,4 +34,5 @@ h.correctBoundaryConditions(); thermo.correct(); + psis = thermo.psi()/thermo.Cp()*thermo.Cv(); } diff --git a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H index 5a214ca34..a85bad7be 100644 --- a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H @@ -1,24 +1,22 @@ { volScalarField rUA = 1.0/UEqn.A(); + surfaceScalarField psisf = fvc::interpolate(psis); + surfaceScalarField rhof = fvc::interpolate(rho); + + // Needs to be outside of loop since p is changing, but psi and rho are not. + surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p); + for (int corr = 0; corr < nCorr; corr++) { U = rUA*UEqn.H(); - // Execute ddtPhiCorr before recalculating flux - // HJ, 27/Apr/2010 - surfaceScalarField phid - ( - "phid", - fvc::interpolate(psi)* - ( - (fvc::interpolate(U) & mesh.Sf()) - + fvc::ddtPhiCorr(rUA, rho, U, phi) - ) - ); - // Calculate phi for boundary conditions - phi = fvc::interpolate(rho*U) & mesh.Sf(); + phi = rhof*fvc::interpolate(U) & mesh.Sf(); + + surfaceScalarField phid2 = rhoReff/rhof*phi; + + surfaceScalarField phid("phid", psisf/rhof*phi); p.storePrevIter(); @@ -26,8 +24,9 @@ { fvScalarMatrix pEqn ( - fvm::ddt(psi, p) + fvm::ddt(psis, p) + fvm::div(phid, p) + + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); @@ -42,7 +41,7 @@ // Calculate the flux if (nonOrth == nNonOrthCorr) { - phi = pEqn.flux(); + phi = phid2 + pEqn.flux(); } } diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H index 1cd8b9efb..b45e1573d 100644 --- a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H @@ -11,13 +11,17 @@ { U = rUA*UEqn.H(); - phi = rhoReff*(fvc::interpolate(U) & mesh.Sf()); + // Calculate phi for boundary conditions + phi = rhof*fvc::interpolate(U) & mesh.Sf(); - surfaceScalarField phid("phid", psisf/rhoReff*phi); + surfaceScalarField phid2 = rhoReff/rhof*phi; + + surfaceScalarField phid("phid", psisf/rhof*phi); // Make fluxes relative within the MRF zone - mrfZones.relativeFlux(psisf, phid); mrfZones.relativeFlux(rhoReff, phi); + mrfZones.relativeFlux(psisf, phid); + mrfZones.relativeFlux(rhoReff, phid2); p.storePrevIter(); @@ -27,7 +31,7 @@ ( fvm::ddt(psis, p) + fvm::div(phid, p) - + fvc::div(phi) + + fvc::div(phid2) - fvm::laplacian(rho*rUA, p) ); @@ -42,13 +46,13 @@ // Calculate the flux if (nonOrth == nNonOrthCorr) { - phi += pEqn.flux(); + phi = phid2 + pEqn.flux(); } } # include "compressibleContinuityErrs.H" - // Explicitly relax the pressure for momentum corrector + // Relax the pressure p.relax(); U -= rUA*fvc::grad(p); diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/createFields.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/createFields.H index 49b4b0b97..6b84b46ba 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/createFields.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/createFields.H @@ -9,7 +9,8 @@ volScalarField& p = thermo.p(); volScalarField& h = thermo.h(); const volScalarField& T = thermo.T(); - const volScalarField& psi = thermo.psi(); + volScalarField psis("psi", thermo.psi()/thermo.Cp()*thermo.Cv()); + psis.oldTime(); volScalarField rho ( diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H index c80d715ca..d91a1ffb3 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H @@ -1,5 +1,5 @@ { - // Solve the enthalpy equation + // Solve the rothalpy equation T.storePrevIter(); // Calculate face velocity from flux @@ -44,4 +44,5 @@ i = h - 0.5*magSqr(Urot); thermo.correct(); + psis = thermo.psi()/thermo.Cp()*thermo.Cv(); } diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H index ed734fc7d..f161a3c4d 100644 --- a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H +++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H @@ -1,24 +1,22 @@ { volScalarField rUrelA = 1.0/UrelEqn.A(); + surfaceScalarField psisf = fvc::interpolate(psis); + surfaceScalarField rhof = fvc::interpolate(rho); + + // Needs to be outside of loop since p is changing, but psi and rho are not. + surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p); + for (int corr = 0; corr < nCorr; corr++) { Urel = rUrelA*UrelEqn.H(); - // Execute ddtPhiCorr before recalculating flux - // HJ, 27/Apr/2010 - surfaceScalarField phid - ( - "phid", - fvc::interpolate(thermo.psi())* - ( - (fvc::interpolate(Urel) & mesh.Sf()) - + fvc::ddtPhiCorr(rUrelA, rho, Urel, phi) - ) - ); - // Calculate phi for boundary conditions - phi = fvc::interpolate(rho*Urel) & mesh.Sf(); + phi = rhof*fvc::interpolate(Urel) & mesh.Sf(); + + surfaceScalarField phid2 = rhoReff/rhof*phi; + + surfaceScalarField phid("phid", psisf/rhof*phi); p.storePrevIter(); @@ -26,8 +24,9 @@ { fvScalarMatrix pEqn ( - fvm::ddt(psi, p) + fvm::ddt(psis, p) + fvm::div(phid, p) + + fvc::div(phid2) - fvm::laplacian(rho*rUrelA, p) ); @@ -42,14 +41,13 @@ // Calculate the flux if (nonOrth == nNonOrthCorr) { - phi = pEqn.flux(); + phi = phid2 + pEqn.flux(); } } # include "compressibleContinuityErrs.H" - - // Explicitly relax the pressure for momentum corrector + // Relax the pressure p.relax(); Urel -= rUrelA*fvc::grad(p); diff --git a/tutorials/compressible/steadyCompressibleFoam/2bump/system/fvSolution b/tutorials/compressible/steadyCompressibleFoam/2bump/system/fvSolution index c567ebbed..cf8cf1e2c 100644 --- a/tutorials/compressible/steadyCompressibleFoam/2bump/system/fvSolution +++ b/tutorials/compressible/steadyCompressibleFoam/2bump/system/fvSolution @@ -58,11 +58,9 @@ PIMPLE relaxationFactors { // Note: under-relaxation factors used in wave-transmissive schemes - U 0.5; - p 0.2; - h 0.5; - rho 0.5; - T 1; + U 0.7; + p 0.7; + h 0.7; } fieldBounds diff --git a/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSchemes b/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSchemes index 1bb526f98..8c0c582e2 100644 --- a/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSchemes +++ b/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSchemes @@ -20,7 +20,7 @@ ddtSchemes ddt(rho,U) steadyState; ddt(rho,h) steadyState; - ddt(psi,p) steadyInertial phi rho 0.5; + ddt(psi,p) steadyInertial phi rho 0.95; U steadyState; T steadyState; diff --git a/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution b/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution index 5a67d8053..e499e4b91 100644 --- a/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution +++ b/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution @@ -76,8 +76,6 @@ relaxationFactors U 0.1; p 0.25; h 0.1; - rho 0.5; - T 0.5; } fieldBounds