research/OpenFOAM-2.4.x
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

rhoCentralFoam: Corrected the BCs for positive (outgoing) fluxes

Browse source 
Resolves bug-report http://openfoam.org/mantisbt/view.php?id=1548

Also upgraded rhoCentralDyMFoam to the version in OpenFOAM-dev which
support mesh-motion and topology change
This commit is contained in:
Henry 2015-02-26 17:08:09 +00:00
parent f219e26d11
commit 11ddd07109
5 changed files with 218 additions and 154 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

49
applications/solvers/compressible/rhoCentralFoam/centralCourantNo.H Normal file
View file

@ -0,0 +1,49 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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.
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Global
centralCourantNo
Description
Calculates the mean and maximum wave speed based Courant Numbers.
\*---------------------------------------------------------------------------*/
if (mesh.nInternalFaces())
{
surfaceScalarField amaxSfbyDelta
(
mesh.surfaceInterpolation::deltaCoeffs()*amaxSf
);
CoNum = max(amaxSfbyDelta/mesh.magSf()).value()*runTime.deltaTValue();
meanCoNum =
(sum(amaxSfbyDelta)/sum(mesh.magSf())).value()
*runTime.deltaTValue();
}
Info<< "Mean and max Courant Numbers = "
<< meanCoNum << " " << CoNum << endl;
// ************************************************************************* //

45
applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H Normal file
View file

@ -0,0 +1,45 @@
namespace Foam
{
//- Interpolate field vf according to direction dir
template<class Type>
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate
(
const GeometricField<Type, fvPatchField, volMesh>& vf,
const surfaceScalarField& dir,
const word& reconFieldName = word::null
)
{
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > tsf
(
fvc::interpolate
(
vf,
dir,
"reconstruct("
+ (reconFieldName != word::null ? reconFieldName : vf.name())
+ ')'
)
);
GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf();
sf.rename(vf.name() + '_' + dir.name());
// Correct BCs of the positive (outgoing) fluxes
if (dir[0] > 0)
{
forAll(sf.boundaryField(), patchi)
{
if (!sf.boundaryField()[patchi].coupled())
{
sf.boundaryField()[patchi] =
vf.boundaryField()[patchi].patchInternalField();
}
}
}
return tsf;
}
}

3
applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/options
View file

@ -6,6 +6,7 @@ EXE_INC = \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
@ -17,4 +18,6 @@ EXE_LIBS = \
-lcompressibleRASModels \
-lcompressibleLESModels \
-ldynamicMesh \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-lmeshTools

204
applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -26,15 +26,17 @@ Application
Description
Density-based compressible flow solver based on central-upwind schemes of
Kurganov and Tadmor
Kurganov and Tadmor with support for mesh-motion and topology changes
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "psiThermo.H"
#include "turbulenceModel.H"
#include "zeroGradientFvPatchFields.H"
#include "fixedRhoFvPatchScalarField.H"
#include "directionInterpolate.H"
#include "motionSolver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -42,9 +44,8 @@ Description
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createDynamicFvMesh.H"
#include "createFields.H"
#include "readTimeControls.H"
@ -54,107 +55,14 @@ int main(int argc, char *argv[])
dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0);
Info<< "\nStarting time loop\n" << endl;
// Courant numbers used to adjust the time-step
scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
autoPtr<Foam::motionSolver> motionPtr = motionSolver::New(mesh);
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
// --- upwind interpolation of primitive fields on faces
surfaceScalarField rho_pos
(
fvc::interpolate(rho, pos, "reconstruct(rho)")
);
surfaceScalarField rho_neg
(
fvc::interpolate(rho, neg, "reconstruct(rho)")
);
surfaceVectorField rhoU_pos
(
fvc::interpolate(rhoU, pos, "reconstruct(U)")
);
surfaceVectorField rhoU_neg
(
fvc::interpolate(rhoU, neg, "reconstruct(U)")
);
volScalarField rPsi(1.0/psi);
surfaceScalarField rPsi_pos
(
fvc::interpolate(rPsi, pos, "reconstruct(T)")
);
surfaceScalarField rPsi_neg
(
fvc::interpolate(rPsi, neg, "reconstruct(T)")
);
surfaceScalarField e_pos
(
fvc::interpolate(e, pos, "reconstruct(T)")
);
surfaceScalarField e_neg
(
fvc::interpolate(e, neg, "reconstruct(T)")
);
surfaceVectorField U_pos(rhoU_pos/rho_pos);
surfaceVectorField U_neg(rhoU_neg/rho_neg);
surfaceScalarField p_pos(rho_pos*rPsi_pos);
surfaceScalarField p_neg(rho_neg*rPsi_neg);
surfaceScalarField phiv_pos(U_pos & mesh.Sf());
surfaceScalarField phiv_neg(U_neg & mesh.Sf());
fvc::makeRelative(phiv_pos, U);
fvc::makeRelative(phiv_neg, U);
volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi));
surfaceScalarField cSf_pos
(
fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf()
);
surfaceScalarField cSf_neg
(
fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf()
);
surfaceScalarField ap
(
max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
);
surfaceScalarField am
(
min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
);
surfaceScalarField a_pos(ap/(ap - am));
surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
surfaceScalarField aSf(am*a_pos);
if (fluxScheme == "Tadmor")
{
aSf = -0.5*amaxSf;
a_pos = 0.5;
}
surfaceScalarField a_neg(1.0 - a_pos);
phiv_pos *= a_pos;
phiv_neg *= a_neg;
surfaceScalarField aphiv_pos(phiv_pos - aSf);
surfaceScalarField aphiv_neg(phiv_neg + aSf);
// Reuse amaxSf for the maximum positive and negative fluxes
// estimated by the central scheme
amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
#include "compressibleCourantNo.H"
#include "readTimeControls.H"
#include "setDeltaT.H"
@ -162,7 +70,88 @@ int main(int argc, char *argv[])
Info<< "Time = " << runTime.timeName() << nl << endl;
mesh.movePoints(motionPtr->newPoints());
// Do any mesh changes
mesh.update();
// --- Directed interpolation of primitive fields onto faces
surfaceScalarField rho_pos(interpolate(rho, pos));
surfaceScalarField rho_neg(interpolate(rho, neg));
surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
volScalarField rPsi("rPsi", 1.0/psi);
surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
surfaceScalarField e_pos(interpolate(e, pos, T.name()));
surfaceScalarField e_neg(interpolate(e, neg, T.name()));
surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos);
surfaceScalarField p_neg("p_neg", rho_neg*rPsi_neg);
surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
// Make fluxes relative to mesh-motion
if (mesh.moving())
{
phiv_pos -= mesh.phi();
phiv_neg -= mesh.phi();
}
volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
surfaceScalarField cSf_pos
(
"cSf_pos",
interpolate(c, pos, T.name())*mesh.magSf()
);
surfaceScalarField cSf_neg
(
"cSf_neg",
interpolate(c, neg, T.name())*mesh.magSf()
);
surfaceScalarField ap
(
"ap",
max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero)
);
surfaceScalarField am
(
"am",
min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero)
);
surfaceScalarField a_pos("a_pos", ap/(ap - am));
surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
surfaceScalarField aSf("aSf", am*a_pos);
if (fluxScheme == "Tadmor")
{
aSf = -0.5*amaxSf;
a_pos = 0.5;
}
surfaceScalarField a_neg("a_neg", 1.0 - a_pos);
phiv_pos *= a_pos;
phiv_neg *= a_neg;
surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);
surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf);
// Reuse amaxSf for the maximum positive and negative fluxes
// estimated by the central scheme
amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
#include "centralCourantNo.H"
phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
@ -174,13 +163,19 @@ int main(int argc, char *argv[])
surfaceScalarField phiEp
(
"phiEp",
aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
+ aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
+ mesh.phi()*(a_pos*p_pos + a_neg*p_neg)
+ aSf*p_pos - aSf*p_neg
);
volScalarField muEff(turbulence->muEff());
// Make flux for pressure-work absolute
if (mesh.moving())
{
phiEp += mesh.phi()*(a_pos*p_pos + a_neg*p_neg);
}
volScalarField muEff("muEff", turbulence->muEff());
volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
// --- Solve density
@ -209,6 +204,7 @@ int main(int argc, char *argv[])
// --- Solve energy
surfaceScalarField sigmaDotU
(
"sigmaDotU",
(
fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
+ (mesh.Sf() & fvc::interpolate(tauMC))

71
applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C
View file

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -35,6 +35,7 @@ Description
#include "turbulenceModel.H"
#include "zeroGradientFvPatchFields.H"
#include "fixedRhoFvPatchScalarField.H"
#include "directionInterpolate.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -53,56 +54,28 @@ int main(int argc, char *argv[])
dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0);
// Courant numbers used to adjust the time-step
scalar CoNum = 0.0;
scalar meanCoNum = 0.0;
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
// --- upwind interpolation of primitive fields on faces
// --- Directed interpolation of primitive fields onto faces
surfaceScalarField rho_pos
(
"rho_pos",
fvc::interpolate(rho, pos, "reconstruct(rho)")
);
surfaceScalarField rho_neg
(
"rho_neg",
fvc::interpolate(rho, neg, "reconstruct(rho)")
);
surfaceScalarField rho_pos(interpolate(rho, pos));
surfaceScalarField rho_neg(interpolate(rho, neg));
surfaceVectorField rhoU_pos
(
"rhoU_pos",
fvc::interpolate(rhoU, pos, "reconstruct(U)")
);
surfaceVectorField rhoU_neg
(
"rhoU_neg",
fvc::interpolate(rhoU, neg, "reconstruct(U)")
);
surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
volScalarField rPsi(1.0/psi);
surfaceScalarField rPsi_pos
(
"rPsi_pos",
fvc::interpolate(rPsi, pos, "reconstruct(T)")
);
surfaceScalarField rPsi_neg
(
"rPsi_neg",
fvc::interpolate(rPsi, neg, "reconstruct(T)")
);
volScalarField rPsi("rPsi", 1.0/psi);
surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
surfaceScalarField e_pos
(
"e_pos",
fvc::interpolate(e, pos, "reconstruct(T)")
);
surfaceScalarField e_neg
(
"e_neg",
fvc::interpolate(e, neg, "reconstruct(T)")
);
surfaceScalarField e_pos(interpolate(e, pos, T.name()));
surfaceScalarField e_neg(interpolate(e, neg, T.name()));
surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
@ -113,16 +86,16 @@ int main(int argc, char *argv[])
surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());
surfaceScalarField phiv_neg("phiv_neg", U_neg & mesh.Sf());
volScalarField c(sqrt(thermo.Cp()/thermo.Cv()*rPsi));
volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi));
surfaceScalarField cSf_pos
(
"cSf_pos",
fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf()
interpolate(c, pos, T.name())*mesh.magSf()
);
surfaceScalarField cSf_neg
(
"cSf_neg",
fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf()
interpolate(c, neg, T.name())*mesh.magSf()
);
surfaceScalarField ap
@ -160,7 +133,7 @@ int main(int argc, char *argv[])
// estimated by the central scheme
amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
#include "compressibleCourantNo.H"
#include "centralCourantNo.H"
#include "readTimeControls.H"
#include "setDeltaT.H"
@ -184,7 +157,7 @@ int main(int argc, char *argv[])
+ aSf*p_pos - aSf*p_neg
);
volScalarField muEff(turbulence->muEff());
volScalarField muEff("muEff", turbulence->muEff());
volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U))));
// --- Solve density
@ -199,8 +172,6 @@ int main(int argc, char *argv[])
U.correctBoundaryConditions();
rhoU.boundaryField() = rho.boundaryField()*U.boundaryField();
volScalarField rhoBydt(rho/runTime.deltaT());
if (!inviscid)
{
solve