Bugfix: use total enthalpy

This commit is contained in:
Hrvoje Jasak 2016-03-15 12:56:56 +00:00
parent c1ad46b28f
commit 5ad3e2dbd8
34 changed files with 132 additions and 173 deletions

View file

@ -40,6 +40,13 @@
mesh
);
// Store kinetic energy for total enthalpy equation
volScalarField K
(
"K",
0.5*(magSqr(U))
);
# include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;

View file

@ -1,23 +1,16 @@
{
// Solve the enthalpy equation
T.storePrevIter();
// Solve the enthalpy equation in total enthalpy formulation (see K)
// Calculate face velocity from flux
surfaceScalarField faceU
(
"faceU",
phi/fvc::interpolate(rho)
);
K = 0.5*(magSqr(U));
fvScalarMatrix hEqn
(
fvm::ddt(rho, h)
+ fvc::ddt(rho, K)
+ fvm::div(phi, h)
+ fvc::div(phi, K)
- fvm::laplacian(turbulence->alphaEff(), h)
==
// ddt(p) term removed: steady-state. HJ, 27/Apr/2010
fvc::div(faceU, p, "div(U,p)")
- p*fvc::div(faceU)
// Viscous heating: note sign (devRhoReff has a minus in it)
- (turbulence->devRhoReff() && fvc::grad(U))
);

View file

@ -40,6 +40,13 @@
mesh
);
// Store kinetic energy for total enthalpy equation
volScalarField K
(
"K",
0.5*(magSqr(U))
);
# include "compressibleCreatePhi.H"
Info<< "Creating turbulence model\n" << endl;

View file

@ -1,25 +1,12 @@
{
// Solve the rothalpy equation
T.storePrevIter();
// Create relative velocity
Urel == U;
mrfZones.relativeVelocity(Urel);
// Create rotational velocity (= omega x r)
Urot = U - Urel;
// Calculate face velocity from absolute flux
surfaceScalarField rhof = fvc::interpolate(rho);
surfaceScalarField phiAbs
(
"phiAbs",
phi
);
mrfZones.absoluteFlux(rhof, phiAbs);
surfaceScalarField faceU("faceU", phiAbs/rhof);
Urot == U - Urel;
fvScalarMatrix iEqn
(
@ -27,9 +14,6 @@
+ fvm::div(phi, i)
- fvm::laplacian(turbulence->alphaEff(), i)
==
// ddt(p) term removed: steady-state. HJ, 27/Apr/2010
fvc::div(faceU, p, "div(U,p)")
- p*fvc::div(faceU)
// Viscous heating: note sign (devRhoReff has a minus in it)
- (turbulence->devRhoReff() && fvc::grad(Urel))
);
@ -40,7 +24,7 @@
maxResidual = max(eqnResidual, maxResidual);
// Calculate enthalpy out of rothalpy
h = i + 0.5*magSqr(Urot);
h = i + 0.5*(magSqr(Urot) - magSqr(Urel));
h.correctBoundaryConditions();
thermo.correct();

View file

@ -73,6 +73,9 @@
SRF::SRFModel::New(Urel)
);
// Create rotational velocity
volVectorField Urot("Urot", SRF->U());
// Create absolute velocity field
volVectorField Uabs
(
@ -84,7 +87,7 @@
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
Urel + SRF->U()
Urel + Urot
);
// Create rothalpy, in two steps to preserve boundary conditions
@ -100,4 +103,4 @@
),
h
);
i -= 0.5*magSqr(SRF->U());
i == h - 0.5*(magSqr(Urot) - magSqr(Urel));

View file

@ -1,13 +1,5 @@
{
// Solve the rothalpy equation
T.storePrevIter();
// Calculate face velocity from flux
surfaceScalarField faceU
(
"faceU",
phi/fvc::interpolate(rho) + (SRF->faceU() & mesh.Sf())
);
fvScalarMatrix iEqn
(
@ -15,9 +7,6 @@
+ fvm::div(phi, i)
- fvm::laplacian(turbulence->alphaEff(), i)
==
// ddt(p) term removed: steady-state. HJ, 27/Apr/2010
fvc::div(faceU, p, "div(U,p)")
- p*fvc::div(faceU)
// Viscous heating: note sign (devRhoReff has a minus in it)
- (turbulence->devRhoReff() && fvc::grad(Urel))
);
@ -28,15 +17,9 @@
maxResidual = max(eqnResidual, maxResidual);
// Calculate enthalpy out of rothalpy
volVectorField Urot("Urot", SRF->U());
h = i + 0.5*magSqr(Urot);
h = i + 0.5*(magSqr(Urot) - magSqr(Urel));
h.correctBoundaryConditions();
// Re-initialise rothalpy based on limited enthalpy
i = h - 0.5*magSqr(Urot);
// Bounding of enthalpy taken out
thermo.correct();
psis = thermo.psi()/thermo.Cp()*thermo.Cv();
}

View file

@ -4,9 +4,10 @@
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.
// Needs to be outside of loop since p is changing, but psi and rho are not
surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p);
// --- PISO loop
for (int corr = 0; corr < nCorr; corr++)
{
Urel = rUrelA*UrelEqn.H();

View file

@ -65,15 +65,22 @@ int main(int argc, char *argv[])
# include "initConvergenceCheck.H"
# include "UEqn.H"
# include "iEqn.H"
# include "pEqn.H"
// # include "hEqn.H"
// Solving for rothalpy
# include "iEqn.H"
# include "rhoFromP.H"
// Correct turbulence
turbulence->correct();
Uabs = Urel + SRF->U();
// Update rotational velocity
Urot = SRF->U();
// Update absolute velocity
Uabs = Urel + Urot;
runTime.write();

View file

@ -26,7 +26,7 @@ deltaT 1;
writeControl timeStep;
writeInterval 50;
writeInterval 100;
purgeWrite 0;
@ -34,7 +34,7 @@ writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
writeCompression compressed;
timeFormat general;

View file

@ -20,8 +20,8 @@ ddtSchemes
ddt(rho,U) steadyState;
ddt(rho,h) steadyState;
// ddt(rho,h) steadyInertial phi rho 1;
ddt(psi,p) steadyInertial phi rho 1;
ddt(psi,p) steadyState;
ddt(rho,K) steadyState;
U steadyState;
T steadyState;
@ -36,15 +36,11 @@ gradSchemes
divSchemes
{
default none;
// div(phi,U) Gauss upwind;
// div(phi,h) Gauss upwind;
// div(phid,p) Gauss upwind;
div(phi,U) Gauss vanLeerDC;
div(phi,h) Gauss vanLeerDC;
div(phid,p) Gauss vanLeer;
div(phid,p) Gauss upwind;
div(U,p) Gauss linear;
div(phi,K) Gauss linear;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;

View file

@ -24,28 +24,28 @@ solvers
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.01;
};
relTol 0.0;
}
U
{
solver BiCGStab;
preconditioner DILU;
minIter 0;
maxIter 1000;
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.01;
};
relTol 0.0;
}
h
{
solver BiCGStab;
preconditioner DILU;
minIter 0;
maxIter 1000;
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.01;
};
relTol 0.0;
}
}
PIMPLE
@ -58,18 +58,13 @@ PIMPLE
relaxationFactors
{
// Note: under-relaxation factors used in wave-transmissive schemes
U 0.7;
p 0.7;
h 0.7;
U 0.5;
p 0.3;
h 0.5;
}
fieldBounds
{
// No bounding
// p 0 1e7;
// T 0 10000;
// U 1e6;
// With bounding
p 50 1e6;
T 20 3000;

View file

@ -22,14 +22,14 @@ boundaryField
{
INLE1
{
type fixedValue;
type pressureInletOutletVelocity;
value uniform (234.636777 0 0);
}
PRES2
{
type inletOutlet;
inletValue uniform (0 0 0);
type pressureInletOutletVelocity;
value uniform (234.636777 0 0);
}
WALL3

View file

@ -3,7 +3,7 @@
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.2 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{

View file

@ -20,8 +20,9 @@ ddtSchemes
ddt(rho,U) steadyState;
ddt(rho,h) steadyState;
ddt(psi,p) steadyInertial phi rho 0.95;
ddt(psi,p) steadyState;
ddt(rho,K) steadyState;
U steadyState;
T steadyState;
p steadyState;
@ -29,7 +30,7 @@ ddtSchemes
gradSchemes
{
default Gauss linear;
default cellLimited Gauss linear 1;
}
divSchemes
@ -37,22 +38,13 @@ divSchemes
default none;
div(phi,U) Gauss upwind;
// div(phi,U) Gauss linearUpwind faceLimited Gauss linear 1.0;
// div(phi,U) Gauss vanLeerDC;
div(phi,h) Gauss upwind;
// div(phi,h) Gauss vanLeerDC;
div(phid,p) Gauss upwind;
// div(phid,p) Gauss vanLeerDC;
div(phi,K) Gauss linear;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div((muEff*dev2(grad(U).T()))) Gauss linear;
div(U,p) Gauss linear;
div(U) Gauss linear;
}
laplacianSchemes

View file

@ -32,19 +32,18 @@ solvers
solver BiCGStab;
preconditioner DILU;
minIter 0;
maxIter 1000;
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.0;
}
h
{
solver BiCGStab;
preconditioner DILU;
minIter 0;
maxIter 1000;
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.0;
}
@ -53,16 +52,16 @@ solvers
PIMPLE
{
nOuterCorrectors 1;
nCorrectors 3;
nCorrectors 2;
nNonOrthogonalCorrectors 0;
}
relaxationFactors
{
// Note: under-relaxation factors used in wave-transmissive schemes
U 0.1;
p 0.25;
h 0.1;
U 0.4;
p 0.3;
h 0.5;
}
fieldBounds

View file

@ -3,7 +3,7 @@
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.2 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{

View file

@ -20,7 +20,7 @@ ddtSchemes
ddt(rho,U) steadyState;
ddt(rho,i) steadyState;
ddt(psi,p) steadyInertial phi rho 0.25;
ddt(psi,p) steadyState;
ddt(rho,k) steadyState;
ddt(rho,epsilon) steadyState;

View file

@ -31,8 +31,8 @@ solvers
solver BiCGStab;
preconditioner DILU;
minIter 0;
maxIter 1000;
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.01;
}
@ -78,10 +78,9 @@ PIMPLE
relaxationFactors
{
// Note: under-relaxation factors used in wave-transmissive schemes
U 0.7;
p 1;
i 0.1;
rho 0.25;
U 0.5;
p 0.3;
i 0.5;
k 0.5;
epsilon 0.5;
@ -89,11 +88,6 @@ relaxationFactors
fieldBounds
{
// No bounding
// p 0 1e7;
// T 0 10000;
// U 1e6;
// With bounding
p 2e4 1e6;
T 200 500;

View file

@ -20,7 +20,7 @@ ddtSchemes
ddt(rho,U) steadyState;
ddt(rho,i) steadyState;
ddt(psi,p) steadyInertial phi rho 0.25;
ddt(psi,p) steadyState;
ddt(rho,k) steadyState;
ddt(rho,epsilon) steadyState;

View file

@ -31,8 +31,8 @@ solvers
solver BiCGStab;
preconditioner DILU;
minIter 0;
maxIter 1000;
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.01;
}
@ -79,13 +79,11 @@ relaxationFactors
{
// Note: under-relaxation factors used in wave-transmissive schemes
U 0.5;
p 0.2;
i 0.1;
h 0.5;
rho 0.25;
p 0.3;
i 0.5;
k 0.2;
epsilon 0.2;
k 0.5;
epsilon 0.5;
}
fieldBounds

View file

@ -22,7 +22,7 @@ boundaryField
{
inlet
{
type totalTemperature;
type isentropicTotalTemperature;
phi phi;
rho none;
psi psi;
@ -35,7 +35,6 @@ boundaryField
outlet
{
type zeroGradient;
value $internalField;
}
blade

View file

@ -28,8 +28,9 @@ boundaryField
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
type inletOutlet;
inletValue uniform (0 0 0);
value uniform (100 0 0);
}
blade

View file

@ -22,7 +22,7 @@ boundaryField
{
inlet
{
type totalPressure;
type isentropicTotalPressure;
phi phi;
rho none;
psi psi;

View file

@ -3,7 +3,7 @@
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.2 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{

View file

@ -34,7 +34,7 @@ writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
writeCompression compressed;
timeFormat general;

View file

@ -20,7 +20,7 @@ ddtSchemes
ddt(rho,Urel) steadyState;
ddt(rho,i) steadyState;
ddt(psi,p) steadyInertial phi rho 0.25;
ddt(psi,p) steadyState;
ddt(rho,k) steadyState;
ddt(rho,epsilon) steadyState;
@ -38,21 +38,21 @@ gradSchemes
divSchemes
{
default none;
div(phi,Urel) Gauss upwind;
div(phi,i) Gauss upwind;
div(phid,p) Gauss upwind;
div(phi,Urel) Gauss upwind;
div(phi,i) Gauss upwind;
div(phid,p) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div((muEff*dev2(grad(Urel).T()))) Gauss linear;
div(U,p) Gauss linear;
div(U,p) Gauss linear;
}
laplacianSchemes
{
default Gauss linear corrected;
default Gauss linear limited 0.5;
}
interpolationSchemes
@ -62,7 +62,7 @@ interpolationSchemes
snGradSchemes
{
default corrected;
default limited 0.5;
}
fluxRequired

View file

@ -31,10 +31,10 @@ solvers
solver BiCGStab;
preconditioner DILU;
minIter 0;
maxIter 1000;
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.01;
relTol 0.0;
}
i
{
@ -44,7 +44,7 @@ solvers
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.01;
relTol 0.0;
}
k
{
@ -79,12 +79,12 @@ relaxationFactors
{
// Note: under-relaxation factors used in wave-transmissive schemes
Urel 0.5;
p 0.2;
p 0.3;
i 0.1;
rho 0.5;
k 0.2;
epsilon 0.2;
k 0.5;
epsilon 0.5;
}
fieldBounds

View file

@ -22,7 +22,7 @@ boundaryField
{
inlet
{
type totalTemperature;
type isentropicTotalTemperature;
phi phi;
rho none;
psi psi;
@ -35,7 +35,6 @@ boundaryField
outlet
{
type zeroGradient;
value $internalField;
}
blade

View file

@ -28,8 +28,9 @@ boundaryField
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
type inletOutlet;
inletValue uniform (0 0 0);
value uniform (100 0 0);
}
blade

View file

@ -22,7 +22,7 @@ boundaryField
{
inlet
{
type totalPressure;
type isentropicTotalPressure;
phi phi;
rho none;
psi psi;

View file

@ -3,7 +3,7 @@
| \\ / F ield | foam-extend: Open Source CFD |
| \\ / O peration | Version: 3.2 |
| \\ / A nd | Web: http://www.foam-extend.org |
| \\/ M anipulation | |
| \\/ M anipulation | For copyright notice see file Copyright |
\*---------------------------------------------------------------------------*/
FoamFile
{

View file

@ -34,7 +34,7 @@ writeFormat ascii;
writePrecision 6;
writeCompression uncompressed;
writeCompression compressed;
timeFormat general;

View file

@ -20,7 +20,7 @@ ddtSchemes
ddt(rho,Urel) steadyState;
ddt(rho,i) steadyState;
ddt(psi,p) steadyInertial phi rho 0.25;
ddt(psi,p) steadyState;
ddt(rho,k) steadyState;
ddt(rho,epsilon) steadyState;
@ -38,21 +38,21 @@ gradSchemes
divSchemes
{
default none;
div(phi,Urel) Gauss upwind;
div(phi,i) Gauss upwind;
div(phid,p) Gauss upwind;
div(phi,Urel) Gauss upwind;
div(phi,i) Gauss upwind;
div(phid,p) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div((muEff*dev2(grad(Urel).T()))) Gauss linear;
div(U,p) Gauss linear;
div(U,p) Gauss linear;
}
laplacianSchemes
{
default Gauss linear corrected;
default Gauss linear limited 0.5;
}
interpolationSchemes
@ -62,7 +62,7 @@ interpolationSchemes
snGradSchemes
{
default corrected;
default limited 0.5;
}
fluxRequired

View file

@ -31,10 +31,10 @@ solvers
solver BiCGStab;
preconditioner DILU;
minIter 0;
maxIter 1000;
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.01;
relTol 0.0;
}
i
{
@ -44,7 +44,7 @@ solvers
minIter 0;
maxIter 1000;
tolerance 1e-8;
relTol 0.01;
relTol 0.0;
}
k
{
@ -79,12 +79,12 @@ relaxationFactors
{
// Note: under-relaxation factors used in wave-transmissive schemes
Urel 0.5;
p 0.2;
p 0.3;
i 0.1;
rho 0.5;
k 0.2;
epsilon 0.2;
k 0.5;
epsilon 0.5;
}
fieldBounds