diff --git a/applications/solvers/compressible/rhoCentralFoam/centralCourantNo.H b/applications/solvers/compressible/rhoCentralFoam/centralCourantNo.H new file mode 100644 index 00000000..09c13ff2 --- /dev/null +++ b/applications/solvers/compressible/rhoCentralFoam/centralCourantNo.H @@ -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 . + +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; + +// ************************************************************************* // diff --git a/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H b/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H new file mode 100644 index 00000000..443bab5f --- /dev/null +++ b/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H @@ -0,0 +1,45 @@ +namespace Foam +{ + +//- Interpolate field vf according to direction dir +template +tmp > interpolate +( + const GeometricField& vf, + const surfaceScalarField& dir, + const word& reconFieldName = word::null +) +{ + tmp > tsf + ( + fvc::interpolate + ( + vf, + dir, + "reconstruct(" + + (reconFieldName != word::null ? reconFieldName : vf.name()) + + ')' + ) + ); + + GeometricField& 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; +} + +} diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/options b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/options index 33707728..55856a5b 100644 --- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/options +++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/options @@ -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 diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C index a4282534..b9fca440 100644 --- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C +++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C @@ -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 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)) diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C index 3ae35694..1c4cc517 100644 --- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C +++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C @@ -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