Improved formulation of steady compressible solver suite

This commit is contained in:
Hrvoje Jasak 2017-04-04 10:36:16 +01:00
parent 7f020db73e
commit ef86f503f7
17 changed files with 39 additions and 130 deletions

View file

@ -4,7 +4,6 @@
fvVectorMatrix UEqn
(
fvm::div(phi, U)
+ fvm::SuSp(-fvc::div(phi), U)
+ turbulence->divDevRhoReff()
);

View file

@ -21,12 +21,6 @@
p.storePrevIter();
volScalarField divPhid
(
"divPhid",
fvc::div(phid)
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
@ -45,7 +39,8 @@
}
}
# include "compressibleContinuityErrs.H"
// Use incompressible continuity error check: div(rho U) = 0
# include "continuityErrs.H"
// Relax the pressure
p.relax();

View file

@ -1,7 +1,6 @@
fvVectorMatrix UEqn
(
fvm::div(phi, U)
+ fvm::SuSp(-fvc::div(phi), U)
+ turbulence->divDevRhoReff()
);

View file

@ -24,6 +24,7 @@
);
hEqn.relax();
hEqn.solve();
// Bounding of enthalpy taken out
thermo.correct();

View file

@ -47,7 +47,8 @@
}
}
# include "compressibleContinuityErrs.H"
// Use incompressible continuity error check: div(rho U) = 0
# include "continuityErrs.H"
// Relax the pressure
p.relax();

View file

@ -4,7 +4,6 @@
fvVectorMatrix UrelEqn
(
fvm::div(phi, Urel)
+ fvm::SuSp(-fvc::div(phi), Urel)
+ turbulence->divDevRhoReff()
+ rho*SRF->Su()
);

View file

@ -25,8 +25,8 @@
fvScalarMatrix pEqn
(
fvm::div(phid, p)
- fvm::laplacian(rho*rUrelA, p)
+ fvc::div(phid2)
- fvm::laplacian(rho*rUrelA, p)
);
pEqn.solve();
@ -38,7 +38,8 @@
}
}
# include "compressibleContinuityErrs.H"
// Use incompressible continuity error check: div(rho U) = 0
# include "continuityErrs.H"
// Relax the pressure
p.relax();

View file

@ -3,8 +3,7 @@
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
fvm::div(phi, U)
+ turbulence->divDevRhoReff()
);

View file

@ -10,8 +10,8 @@
fvScalarMatrix hEqn
(
fvm::ddt(rho, h)
+ fvm::div(phi, h)
fvm::div(phi, h)
+ fvm::SuSp(-fvc::div(phi), h)
- fvm::laplacian(turbulence->alphaEff(), h)
==
fvc::div(faceU, p, "div(U,p)")
@ -22,7 +22,6 @@
);
hEqn.relax();
hEqn.solve();
// Bounding of enthalpy taken out

View file

@ -24,8 +24,7 @@
{
fvScalarMatrix pEqn
(
fvm::ddt(psis, p)
+ fvm::div(phid, p)
fvm::div(phid, p)
+ fvc::div(phid2)
- fvm::laplacian(rho*rUA, p)
);
@ -39,8 +38,8 @@
}
}
// Use custom continuity error check
# include "universalContinuityErrs.H"
// Use incompressible continuity error check: div(rho U) = 0
# include "continuityErrs.H"
// Relax the pressure
p.relax();

View file

@ -12,4 +12,5 @@
rho = Foam::min(rho, rhoMax);
rho = Foam::max(rho, rhoMin);
rho.relax();
rho.correctBoundaryConditions();
}

View file

@ -1,49 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Global
continuityErrs
Description
Calculates and prints the continuity errors.
\*---------------------------------------------------------------------------*/
{
volScalarField contErr = fvc::ddt(rho) + fvc::div(phi);
sumLocalContErr = runTime.deltaT().value()*
mag(contErr)().weightedAverage(mesh.V()).value();
globalContErr = runTime.deltaT().value()*
contErr.weightedAverage(mesh.V()).value();
cumulativeContErr += globalContErr;
Info<< "time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr
<< endl;
}
// ************************************************************************* //

View file

@ -3,8 +3,7 @@
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
fvm::div(phi, U)
+ turbulence->divDevRhoReff()
);

View file

@ -3,14 +3,29 @@
Urel == U;
mrfZones.relativeVelocity(Urel);
// Bound the relative velocity to preserve the i to h conversion bound
// HJ, 22/Mar/2017
volScalarField magUrel = mag(Urel);
if (max(magUrel) > UMax)
{
volScalarField UrelLimiter = pos(magUrel - UMax)*UMax/(magUrel + smallU)
+ neg(magUrel - UMax);
UrelLimiter.max(scalar(0));
UrelLimiter.min(scalar(1));
Urel *= UrelLimiter;
Urel.correctBoundaryConditions();
}
// Create rotational velocity (= omega x r)
Urot == U - Urel;
fvScalarMatrix iEqn
(
fvm::ddt(rho, i)
+ fvm::div(phi, i)
- fvm::laplacian(turbulence->alphaEff(), i)
fvm::div(phi, i)
+ fvm::SuSp(-fvc::div(phi), i)
- fvm::laplacian(turbulence->alphaEff(), i)
==
// Viscous heating: note sign (devRhoReff has a minus in it)
- (turbulence->devRhoReff() && fvc::grad(U))
@ -23,8 +38,8 @@
// From rothalpy, calculate enthalpy after solution of rothalpy equation
h = i + 0.5*(magSqr(Urot) - magSqr(Urel));
h.correctBoundaryConditions();
// Update thermo for new h
thermo.correct();
psis = thermo.psi()/thermo.Cp()*thermo.Cv();
}

View file

@ -33,8 +33,7 @@
{
fvScalarMatrix pEqn
(
fvm::ddt(psis, p)
+ fvm::div(phid, p)
fvm::div(phid, p)
+ fvc::div(phid2)
- fvm::laplacian(rho*rUA, p)
);
@ -48,8 +47,8 @@
}
}
// Use custom continuity error check
# include "universalContinuityErrs.H"
// Use incompressible continuity error check: div(rho U) = 0
# include "continuityErrs.H"
// Relax the pressure
p.relax();

View file

@ -12,4 +12,5 @@
rho = Foam::min(rho, rhoMax);
rho = Foam::max(rho, rhoMin);
rho.relax();
rho.correctBoundaryConditions();
}

View file

@ -1,49 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation, either version 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Global
continuityErrs
Description
Calculates and prints the continuity errors.
\*---------------------------------------------------------------------------*/
{
volScalarField contErr = fvc::ddt(rho) + fvc::div(phi);
sumLocalContErr = runTime.deltaT().value()*
mag(contErr)().weightedAverage(mesh.V()).value();
globalContErr = runTime.deltaT().value()*
contErr.weightedAverage(mesh.V()).value();
cumulativeContErr += globalContErr;
Info<< "time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr
<< endl;
}
// ************************************************************************* //