diff --git a/applications/solvers/compressible/rhoCentralFoam/Make/files b/applications/solvers/compressible/rhoCentralFoam/Make/files deleted file mode 100644 index 7c0c3704ee..0000000000 --- a/applications/solvers/compressible/rhoCentralFoam/Make/files +++ /dev/null @@ -1,3 +0,0 @@ -rhoCentralFoam.C - -EXE = $(FOAM_APPBIN)/rhoCentralFoam diff --git a/applications/solvers/compressible/rhoCentralFoam/createFieldRefs.H b/applications/solvers/compressible/rhoCentralFoam/createFieldRefs.H deleted file mode 100644 index 5522ccd6aa..0000000000 --- a/applications/solvers/compressible/rhoCentralFoam/createFieldRefs.H +++ /dev/null @@ -1,10 +0,0 @@ -volScalarField& p = thermo.p(); -const volScalarField& T = thermo.T(); -const volScalarField& psi = thermo.psi(); -const volScalarField& mu = thermo.mu(); - -bool inviscid(true); -if (max(mu.primitiveField()) > 0.0) -{ - inviscid = false; -} diff --git a/applications/solvers/compressible/rhoCentralFoam/createFields.H b/applications/solvers/compressible/rhoCentralFoam/createFields.H deleted file mode 100644 index b1a44f1e98..0000000000 --- a/applications/solvers/compressible/rhoCentralFoam/createFields.H +++ /dev/null @@ -1,108 +0,0 @@ -#include "createRDeltaT.H" - -Info<< "Reading thermophysical properties\n" << endl; - -autoPtr pThermo -( - psiThermo::New(mesh) -); -psiThermo& thermo = pThermo(); - -volScalarField& e = thermo.he(); - -Info<< "Reading field U\n" << endl; -volVectorField U -( - IOobject - ( - "U", - runTime.name(), - mesh, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - mesh -); - -volScalarField rho -( - IOobject - ( - "rho", - runTime.name(), - mesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - thermo.renameRho() -); - -volVectorField rhoU -( - IOobject - ( - "rhoU", - runTime.name(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - rho*U -); - -volScalarField rhoE -( - IOobject - ( - "rhoE", - runTime.name(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - rho*(e + 0.5*magSqr(U)) -); - -surfaceScalarField pos -( - IOobject - ( - "pos", - runTime.name(), - mesh - ), - mesh, - dimensionedScalar(dimless, 1.0) -); - -surfaceScalarField neg -( - IOobject - ( - "neg", - runTime.name(), - mesh - ), - mesh, - dimensionedScalar(dimless, -1.0) -); - -surfaceScalarField phi("phi", fvc::flux(rhoU)); - -Info<< "Creating turbulence model\n" << endl; -autoPtr turbulence -( - compressible::momentumTransportModel::New - ( - rho, - U, - phi, - thermo - ) -); - -Info<< "Creating thermophysical transport model\n" << endl; -autoPtr thermophysicalTransport -( - fluidThermoThermophysicalTransportModel::New(turbulence(), thermo) -); diff --git a/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H b/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H deleted file mode 100644 index 058258ba95..0000000000 --- a/applications/solvers/compressible/rhoCentralFoam/directionInterpolate.H +++ /dev/null @@ -1,32 +0,0 @@ -namespace Foam -{ - -//- Interpolate field vf according to direction dir -template -tmp> interpolate -( - const VolField& vf, - const surfaceScalarField& dir, - const word& reconFieldName = word::null -) -{ - tmp> tsf - ( - fvc::interpolate - ( - vf, - dir, - "reconstruct(" - + (reconFieldName != word::null ? reconFieldName : vf.name()) - + ')' - ) - ); - - SurfaceField& sf = tsf.ref(); - - sf.rename(vf.name() + '_' + dir.name()); - - return tsf; -} - -} diff --git a/applications/solvers/compressible/rhoCentralFoam/readFluxScheme.H b/applications/solvers/compressible/rhoCentralFoam/readFluxScheme.H deleted file mode 100644 index 4cba5629f7..0000000000 --- a/applications/solvers/compressible/rhoCentralFoam/readFluxScheme.H +++ /dev/null @@ -1,16 +0,0 @@ -word fluxScheme("Kurganov"); -if (mesh.schemes().dict().readIfPresent("fluxScheme", fluxScheme)) -{ - if ((fluxScheme == "Tadmor") || (fluxScheme == "Kurganov")) - { - Info<< "fluxScheme: " << fluxScheme << endl; - } - else - { - FatalErrorInFunction - << "fluxScheme: " << fluxScheme - << " is not a valid choice. " - << "Options are: Tadmor, Kurganov" - << abort(FatalError); - } -} diff --git a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C b/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C deleted file mode 100644 index f531536a70..0000000000 --- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C +++ /dev/null @@ -1,295 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2023 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 . - -Application - rhoCentralFoam - -Description - Density-based compressible flow solver based on central-upwind schemes of - Kurganov and Tadmor with support for mesh-motion and topology changes. - -\*---------------------------------------------------------------------------*/ - -#include "fvCFD.H" -#include "psiThermo.H" -#include "compressibleMomentumTransportModels.H" -#include "fluidThermoThermophysicalTransportModel.H" -#include "fixedRhoFvPatchScalarField.H" -#include "directionInterpolate.H" -#include "localEulerDdtScheme.H" -#include "fvcSmooth.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -int main(int argc, char *argv[]) -{ - #define NO_CONTROL - #include "postProcess.H" - - #include "setRootCaseLists.H" - #include "createTime.H" - #include "createMesh.H" - #include "createFields.H" - #include "createFieldRefs.H" - #include "createTimeControls.H" - - turbulence->validate(); - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - - #include "readFluxScheme.H" - - dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0); - - // Courant numbers used to adjust the time-step - scalar CoNum = 0.0; - - Info<< "\nStarting time loop\n" << endl; - - while (runTime.run()) - { - #include "readTimeControls.H" - - if (!LTS) - { - #include "setDeltaT.H" - - // Update the mesh for topology change, mesh to mesh mapping - mesh.update(); - - runTime++; - - // Move the mesh - mesh.move(); - } - - // --- Directed interpolation of primitive fields onto faces - - const surfaceScalarField rho_pos(interpolate(rho, pos)); - const surfaceScalarField rho_neg(interpolate(rho, neg)); - - const surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name())); - const surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name())); - - const surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos); - const surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg); - - const volScalarField rPsi("rPsi", 1.0/psi); - const surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name())); - const surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name())); - - const surfaceScalarField p_pos("p_pos", rho_pos*rPsi_pos); - const 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(); - } - - const volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi)); - const surfaceScalarField cSf_pos - ( - "cSf_pos", - interpolate(c, pos, T.name())*mesh.magSf() - ); - const surfaceScalarField cSf_neg - ( - "cSf_neg", - interpolate(c, neg, T.name())*mesh.magSf() - ); - - const surfaceScalarField ap - ( - "ap", - max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero) - ); - const surfaceScalarField am - ( - "am", - min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero) - ); - - const surfaceScalarField a_pos - ( - "a_pos", - fluxScheme == "Tadmor" - ? surfaceScalarField::New("a_pos", mesh, 0.5) - : ap/(ap - am) - ); - - const surfaceScalarField a_neg("a_neg", 1.0 - a_pos); - - phiv_pos *= a_pos; - phiv_neg *= a_neg; - - const surfaceScalarField aSf - ( - "aSf", - fluxScheme == "Tadmor" - ? -0.5*max(mag(am), mag(ap)) - : am*a_pos - ); - - const surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf); - const surfaceScalarField aphiv_neg("aphiv_neg", phiv_neg + aSf); - - { - const surfaceScalarField amaxSf - ( - max(mag(aphiv_pos), mag(aphiv_neg)) - ); - - #include "centralCourantNo.H" - - if (LTS) - { - #include "setRDeltaT.H" - runTime++; - } - } - - Info<< "Time = " << runTime.userTimeName() << nl << endl; - - phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg; - - // --- Solve density - solve(fvm::ddt(rho) + fvc::div(phi)); - - turbulence->predict(); - thermophysicalTransport->predict(); - - const volScalarField muEff("muEff", rho*turbulence->nuEff()); - const volTensorField tauMC("tauMC", muEff*dev2(Foam::T(fvc::grad(U)))); - - // --- Solve momentum - { - const surfaceVectorField phiUp - ( - (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg) - + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf() - ); - - solve(fvm::ddt(rhoU) + fvc::div(phiUp)); - - U.ref() = - rhoU() - /rho(); - U.correctBoundaryConditions(); - rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); - - if (!inviscid) - { - solve - ( - fvm::ddt(rho, U) - fvc::ddt(rho, U) - - fvm::laplacian(muEff, U) - - fvc::div(tauMC) - ); - rhoU = rho*U; - } - } - - // --- Solve energy - { - const surfaceScalarField e_pos(interpolate(e, pos, T.name())); - const surfaceScalarField e_neg(interpolate(e, neg, T.name())); - - 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) - + aSf*p_pos - aSf*p_neg - ); - - // Make flux for pressure-work absolute - if (mesh.moving()) - { - phiEp += mesh.phi()*(a_pos*p_pos + a_neg*p_neg); - } - - solve - ( - fvm::ddt(rhoE) - + fvc::div(phiEp) - ); - - e = rhoE/rho - 0.5*magSqr(U); - e.correctBoundaryConditions(); - thermo.correct(); - rhoE.boundaryFieldRef() == - rho.boundaryField()* - ( - e.boundaryField() + 0.5*magSqr(U.boundaryField()) - ); - - if (!inviscid) - { - const surfaceScalarField sigmaDotU - ( - "sigmaDotU", - ( - fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U) - + fvc::dotInterpolate(mesh.Sf(), tauMC) - ) - & (a_pos*U_pos + a_neg*U_neg) - ); - - solve - ( - fvm::ddt(rho, e) - fvc::ddt(rho, e) - + thermophysicalTransport->divq(e) - - fvc::div(sigmaDotU) - ); - thermo.correct(); - rhoE = rho*(e + 0.5*magSqr(U)); - } - } - - { - p.ref() = rho()/psi(); - p.correctBoundaryConditions(); - rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField(); - } - - turbulence->correct(); - thermophysicalTransport->correct(); - - runTime.write(); - - Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" - << " ClockTime = " << runTime.elapsedClockTime() << " s" - << nl << endl; - } - - Info<< "End\n" << endl; - - return 0; -} - -// ************************************************************************* // diff --git a/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H b/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H deleted file mode 100644 index b126cb82a0..0000000000 --- a/applications/solvers/compressible/rhoCentralFoam/setRDeltaT.H +++ /dev/null @@ -1,29 +0,0 @@ -{ - volScalarField& rDeltaT = trDeltaT.ref(); - - scalar rDeltaTSmoothingCoeff - ( - runTime.controlDict().lookupOrDefault - ( - "rDeltaTSmoothingCoeff", - 0.02 - ) - ); - - // Set the reciprocal time-step from the local Courant number - rDeltaT.ref() = max - ( - 1/dimensionedScalar(dimTime, maxDeltaT), - fvc::surfaceSum(amaxSf)()() - /((2*maxCo)*mesh.V()) - ); - - // Update the boundary values of the reciprocal time-step - rDeltaT.correctBoundaryConditions(); - - fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff); - - Info<< "Flow time scale min/max = " - << gMin(1/rDeltaT.primitiveField()) - << ", " << gMax(1/rDeltaT.primitiveField()) << endl; -} diff --git a/applications/solvers/modules/shockFluid/Make/files b/applications/solvers/modules/shockFluid/Make/files new file mode 100644 index 0000000000..ab60cc9c60 --- /dev/null +++ b/applications/solvers/modules/shockFluid/Make/files @@ -0,0 +1,16 @@ +setRDeltaT.C +moveMesh.C +prePredictor.C +fluxPredictor.C +correctDensity.C +momentumPredictor.C +thermophysicalPredictor.C +pressureCorrector.C +shockFluid.C + +derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C +derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.C +derivedFvPatchFields/T/smoluchowskiJumpTFvPatchScalarField.C +derivedFvPatchFields/rho/fixedRhoFvPatchScalarField.C + +LIB = $(FOAM_LIBBIN)/libshockFluid diff --git a/applications/solvers/compressible/rhoCentralFoam/Make/options b/applications/solvers/modules/shockFluid/Make/options similarity index 81% rename from applications/solvers/compressible/rhoCentralFoam/Make/options rename to applications/solvers/modules/shockFluid/Make/options index f7b1cf9d1e..abed2f88fb 100644 --- a/applications/solvers/compressible/rhoCentralFoam/Make/options +++ b/applications/solvers/modules/shockFluid/Make/options @@ -1,24 +1,23 @@ EXE_INC = \ - -IBCs/lnInclude \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \ -I$(LIB_SRC)/physicalProperties/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \ -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \ -I$(LIB_SRC)/ThermophysicalTransportModels/thermophysicalTransportModel/lnInclude \ -I$(LIB_SRC)/ThermophysicalTransportModels/fluid/lnInclude \ -I$(LIB_SRC)/ThermophysicalTransportModels/fluidThermo/lnInclude \ + -I$(LIB_SRC)/finiteVolume/cfdTools \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude -EXE_LIBS = \ - -lfiniteVolume \ - -lfvModels \ - -lfvConstraints \ +LIB_LIBS = \ + -lfluidSolver \ -lfluidThermophysicalModels \ - -lspecie \ - -lrhoCentralFoam \ -lmomentumTransportModels \ -lcompressibleMomentumTransportModels \ -lcoupledThermophysicalTransportModels \ - -lmeshTools + -lfiniteVolume \ + -lmeshTools \ + -lfvModels \ + -lfvConstraints diff --git a/applications/solvers/compressible/rhoCentralFoam/centralCourantNo.H b/applications/solvers/modules/shockFluid/correctDensity.C similarity index 69% rename from applications/solvers/compressible/rhoCentralFoam/centralCourantNo.H rename to applications/solvers/modules/shockFluid/correctDensity.C index 64973905e6..3862cf654e 100644 --- a/applications/solvers/compressible/rhoCentralFoam/centralCourantNo.H +++ b/applications/solvers/modules/shockFluid/correctDensity.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -21,24 +21,27 @@ License 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. - \*---------------------------------------------------------------------------*/ +#include "shockFluid.H" + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::solvers::shockFluid::correctDensity() { - const scalarField sumAmaxSf(fvc::surfaceSum(amaxSf)().primitiveField()); + fvScalarMatrix rhoEqn + ( + fvm::ddt(rho) + fvc::div(phi) + == + fvModels().source(rho) + ); - CoNum = 0.5*gMax(sumAmaxSf/mesh.V().field())*runTime.deltaTValue(); + fvConstraints().constrain(rhoEqn); - scalar meanCoNum = - 0.5*(gSum(sumAmaxSf)/gSum(mesh.V().field()))*runTime.deltaTValue(); + rhoEqn.solve(); - Info<< "Mean and max Courant Numbers = " - << meanCoNum << " " << CoNum << endl; + fvConstraints().constrain(rho); } + // ************************************************************************* // diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/Make/files b/applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/files similarity index 100% rename from applications/solvers/compressible/rhoCentralFoam/BCs/Make/files rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/files diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/Make/options b/applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/options similarity index 100% rename from applications/solvers/compressible/rhoCentralFoam/BCs/Make/options rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/options diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C b/applications/solvers/modules/shockFluid/derivedFvPatchFields/T/smoluchowskiJumpTFvPatchScalarField.C similarity index 99% rename from applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/T/smoluchowskiJumpTFvPatchScalarField.C index 84dba6d080..32f506e703 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/T/smoluchowskiJumpTFvPatchScalarField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.H b/applications/solvers/modules/shockFluid/derivedFvPatchFields/T/smoluchowskiJumpTFvPatchScalarField.H similarity index 98% rename from applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.H rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/T/smoluchowskiJumpTFvPatchScalarField.H index c866d92412..30bf9f946e 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.H +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/T/smoluchowskiJumpTFvPatchScalarField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C b/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.C similarity index 99% rename from applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.C index 9215cb78e3..60dbaa21cc 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.H b/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.H similarity index 98% rename from applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.H rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.H index 8adc3c2ac8..61a88ec642 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.H +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.C b/applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.C similarity index 98% rename from applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.C rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.C index ed425e357d..825c640f2f 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.C +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.H b/applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.H similarity index 99% rename from applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.H rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.H index b3c98733e9..aa4646761f 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.H +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C b/applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C similarity index 96% rename from applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C index 9ba51165bc..dd77350f8f 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.H b/applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.H similarity index 96% rename from applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.H rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.H index 022da6e257..3426e3a5bb 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.H +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFieldsFwd.H b/applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFieldsFwd.H similarity index 96% rename from applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFieldsFwd.H rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFieldsFwd.H index c3e8e547cf..433e7c1a89 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFieldsFwd.H +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/mixedFixedValueSlip/mixedFixedValueSlipFvPatchFieldsFwd.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.C b/applications/solvers/modules/shockFluid/derivedFvPatchFields/rho/fixedRhoFvPatchScalarField.C similarity index 98% rename from applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.C rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/rho/fixedRhoFvPatchScalarField.C index 626520f093..0cba45e24b 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.C +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/rho/fixedRhoFvPatchScalarField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.H b/applications/solvers/modules/shockFluid/derivedFvPatchFields/rho/fixedRhoFvPatchScalarField.H similarity index 98% rename from applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.H rename to applications/solvers/modules/shockFluid/derivedFvPatchFields/rho/fixedRhoFvPatchScalarField.H index ab74d62202..3451375a0a 100644 --- a/applications/solvers/compressible/rhoCentralFoam/BCs/rho/fixedRhoFvPatchScalarField.H +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/rho/fixedRhoFvPatchScalarField.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License diff --git a/applications/solvers/modules/shockFluid/fluxPredictor.C b/applications/solvers/modules/shockFluid/fluxPredictor.C new file mode 100644 index 0000000000..bdf95b9424 --- /dev/null +++ b/applications/solvers/modules/shockFluid/fluxPredictor.C @@ -0,0 +1,130 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2023 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 . + +\*---------------------------------------------------------------------------*/ + +#include "shockFluid.H" + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::solvers::shockFluid::fluxPredictor() +{ + if (!pos.valid()) + { + pos = surfaceScalarField::New + ( + "pos", + mesh, + dimensionedScalar(dimless, 1) + ); + + neg = surfaceScalarField::New + ( + "neg", + mesh, + dimensionedScalar(dimless, -1.0) + ); + } + + rho_pos = interpolate(rho, pos()); + rho_neg = interpolate(rho, neg()); + + rhoU_pos = interpolate(rhoU, pos(), U.name()); + rhoU_neg = interpolate(rhoU, neg(), U.name()); + + U_pos = surfaceVectorField::New("U_pos", rhoU_pos()/rho_pos()); + U_neg = surfaceVectorField::New("U_neg", rhoU_neg()/rho_neg()); + + const volScalarField& T = thermo.T(); + + const volScalarField rPsi("rPsi", 1.0/thermo.psi()); + const surfaceScalarField rPsi_pos(interpolate(rPsi, pos(), T.name())); + const surfaceScalarField rPsi_neg(interpolate(rPsi, neg(), T.name())); + + p_pos = surfaceScalarField::New("p_pos", rho_pos()*rPsi_pos); + p_neg = surfaceScalarField::New("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(); + } + + const volScalarField c("c", sqrt(thermo.Cp()/thermo.Cv()*rPsi)); + const surfaceScalarField cSf_pos + ( + "cSf_pos", + interpolate(c, pos(), T.name())*mesh.magSf() + ); + const surfaceScalarField cSf_neg + ( + "cSf_neg", + interpolate(c, neg(), T.name())*mesh.magSf() + ); + + const dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0); + + const surfaceScalarField ap + ( + "ap", + max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero) + ); + const surfaceScalarField am + ( + "am", + min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero) + ); + + a_pos = surfaceScalarField::New + ( + "a_pos", + fluxScheme == "Tadmor" + ? surfaceScalarField::New("a_pos", mesh, 0.5) + : ap/(ap - am) + ); + + a_neg = surfaceScalarField::New("a_neg", 1.0 - a_pos()); + + phiv_pos *= a_pos(); + phiv_neg *= a_neg(); + + aSf = surfaceScalarField::New + ( + "aSf", + fluxScheme == "Tadmor" + ? -0.5*max(mag(am), mag(ap)) + : am*a_pos() + ); + + aphiv_pos = surfaceScalarField::New("aphiv_pos", phiv_pos - aSf()); + aphiv_neg = surfaceScalarField::New("aphiv_neg", phiv_neg + aSf()); + + phi = aphiv_pos()*rho_pos() + aphiv_neg()*rho_neg(); +} + + +// ************************************************************************* // diff --git a/applications/solvers/modules/shockFluid/momentumPredictor.C b/applications/solvers/modules/shockFluid/momentumPredictor.C new file mode 100644 index 0000000000..4b5c920947 --- /dev/null +++ b/applications/solvers/modules/shockFluid/momentumPredictor.C @@ -0,0 +1,81 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2023 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 . + +\*---------------------------------------------------------------------------*/ + +#include "shockFluid.H" +#include "fvmDdt.H" +#include "fvmDiv.H" +#include "fvcDdt.H" + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::solvers::shockFluid::momentumPredictor() +{ + if (!inviscid) + { + muEff.clear(); + tauMC.clear(); + + muEff = volScalarField::New("muEff", rho*momentumTransport->nuEff()); + tauMC = new volTensorField + ( + "tauMC", + muEff()*dev2(Foam::T(fvc::grad(U))) + ); + } + + const surfaceVectorField phiUp + ( + (aphiv_pos()*rhoU_pos() + aphiv_neg()*rhoU_neg()) + + (a_pos()*p_pos() + a_neg()*p_neg())*mesh.Sf() + ); + + solve(fvm::ddt(rhoU) + fvc::div(phiUp)); + + U.ref() = rhoU()/rho(); + U.correctBoundaryConditions(); + + if (!inviscid) + { + solve + ( + fvm::ddt(rho, U) - fvc::ddt(rho, U) + - fvm::laplacian(muEff(), U) + - fvc::div(tauMC()) + == + fvModels().source(rho, U) + ); + + rhoU == rho*U; + } + else + { + rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); + } + + fvConstraints().constrain(U); +} + + +// ************************************************************************* // diff --git a/applications/solvers/compressible/rhoCentralFoam/setDeltaT.H b/applications/solvers/modules/shockFluid/moveMesh.C similarity index 68% rename from applications/solvers/compressible/rhoCentralFoam/setDeltaT.H rename to applications/solvers/modules/shockFluid/moveMesh.C index 1f5fd4221a..54ffb27f45 100644 --- a/applications/solvers/compressible/rhoCentralFoam/setDeltaT.H +++ b/applications/solvers/modules/shockFluid/moveMesh.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -21,30 +21,33 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . -Global - setDeltaT - -Description - Reset the timestep to maintain a constant maximum courant Number. - Reduction of time-step is immediate, but increase is damped to avoid - unstable oscillations. - \*---------------------------------------------------------------------------*/ -if (adjustTimeStep) -{ - scalar deltaT = 1.2*runTime.deltaTValue(); +#include "shockFluid.H" - if (CoNum > small) +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +bool Foam::solvers::shockFluid::moveMesh() +{ + if (pimple.firstIter() || pimple.moveMeshOuterCorrectors()) { - deltaT = min(deltaT, maxCo/CoNum*runTime.deltaTValue()); + // Move the mesh + mesh.move(); + + if (mesh.changing()) + { + if (mesh.topoChanged()) + { + } + + meshCourantNo(); + + return true; + } } - deltaT = min(deltaT, maxDeltaT); - - runTime.setDeltaT(deltaT); - - Info<< "deltaT = " << runTime.deltaTValue() << endl; + return false; } + // ************************************************************************* // diff --git a/applications/solvers/modules/shockFluid/prePredictor.C b/applications/solvers/modules/shockFluid/prePredictor.C new file mode 100644 index 0000000000..764b09dd21 --- /dev/null +++ b/applications/solvers/modules/shockFluid/prePredictor.C @@ -0,0 +1,46 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2023 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 . + +\*---------------------------------------------------------------------------*/ + +#include "shockFluid.H" + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::solvers::shockFluid::prePredictor() +{ + fluxPredictor(); + + correctDensity(); + + fvModels().correct(); + + if (!inviscid && pimple.predictTransport()) + { + momentumTransport->predict(); + thermophysicalTransport->predict(); + } +} + + +// ************************************************************************* // diff --git a/applications/solvers/modules/shockFluid/pressureCorrector.C b/applications/solvers/modules/shockFluid/pressureCorrector.C new file mode 100644 index 0000000000..397e8805f5 --- /dev/null +++ b/applications/solvers/modules/shockFluid/pressureCorrector.C @@ -0,0 +1,39 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2023 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 . + +\*---------------------------------------------------------------------------*/ + +#include "shockFluid.H" + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::solvers::shockFluid::pressureCorrector() +{ + const volScalarField& psi = thermo.psi(); + p.ref() = rho()/psi(); + p.correctBoundaryConditions(); + rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField(); +} + + +// ************************************************************************* // diff --git a/applications/solvers/modules/shockFluid/setRDeltaT.C b/applications/solvers/modules/shockFluid/setRDeltaT.C new file mode 100644 index 0000000000..5e1902f807 --- /dev/null +++ b/applications/solvers/modules/shockFluid/setRDeltaT.C @@ -0,0 +1,75 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2023 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 . + +\*---------------------------------------------------------------------------*/ + +#include "shockFluid.H" +#include "fvcSmooth.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::solvers::shockFluid::setRDeltaT(const surfaceScalarField& amaxSf) +{ + volScalarField& rDeltaT = trDeltaT.ref(); + + const dictionary& pimpleDict = pimple.dict(); + + const scalar maxCo + ( + pimpleDict.lookupOrDefault("maxCo", 0.8) + ); + + const scalar rDeltaTSmoothingCoeff + ( + pimpleDict.lookupOrDefault + ( + "rDeltaTSmoothingCoeff", + 0.02 + ) + ); + + const scalar maxDeltaT + ( + pimpleDict.lookupOrDefault("maxDeltaT", great) + ); + + // Set the reciprocal time-step from the local Courant number + rDeltaT.ref() = max + ( + 1/dimensionedScalar(dimTime, maxDeltaT), + fvc::surfaceSum(amaxSf)()() + /((2*maxCo)*mesh.V()) + ); + + // Update the boundary values of the reciprocal time-step + rDeltaT.correctBoundaryConditions(); + + fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff); + + Info<< "Flow time scale min/max = " + << gMin(1/rDeltaT.primitiveField()) + << ", " << gMax(1/rDeltaT.primitiveField()) << endl; +} + + +// ************************************************************************* // diff --git a/applications/solvers/modules/shockFluid/shockFluid.C b/applications/solvers/modules/shockFluid/shockFluid.C new file mode 100644 index 0000000000..1eae1bd1b5 --- /dev/null +++ b/applications/solvers/modules/shockFluid/shockFluid.C @@ -0,0 +1,309 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2023 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 . + +\*---------------------------------------------------------------------------*/ + +#include "shockFluid.H" +#include "localEulerDdtScheme.H" +#include "hydrostaticInitialisation.H" +#include "fvcMeshPhi.H" +#include "fvcVolumeIntegrate.H" +#include "fvcReconstruct.H" +#include "fvcSnGrad.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solvers +{ + defineTypeNameAndDebug(shockFluid, 0); + addToRunTimeSelectionTable(solver, shockFluid, fvMesh); +} +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::solvers::shockFluid::correctCoNum(const surfaceScalarField& amaxSf) +{ + const scalarField sumAmaxSf(fvc::surfaceSum(amaxSf)().primitiveField()); + + CoNum = 0.5*gMax(sumAmaxSf/mesh.V().field())*runTime.deltaTValue(); + + const scalar meanCoNum = + 0.5*(gSum(sumAmaxSf)/gSum(mesh.V().field()))*runTime.deltaTValue(); + + Info<< "Courant Number mean: " << meanCoNum + << " max: " << CoNum << endl; +} + + +// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // + +void Foam::solvers::shockFluid::clearTemporaryFields() +{ + rho_pos.clear(); + rho_neg.clear(); + + rhoU_pos.clear(); + rhoU_neg.clear(); + + U_pos.clear(); + U_neg.clear(); + + p_pos.clear(); + p_neg.clear(); + + a_pos.clear(); + a_neg.clear(); + + aSf.clear(); + + aphiv_pos.clear(); + aphiv_neg.clear(); + + muEff.clear(); + tauMC.clear(); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::solvers::shockFluid::shockFluid(fvMesh& mesh) +: + fluidSolver(mesh), + + thermo_(psiThermo::New(mesh)), + + thermo(thermo_()), + + p(thermo.p()), + + rho + ( + IOobject + ( + "rho", + runTime.name(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + thermo.renameRho() + ), + + U + ( + IOobject + ( + "U", + runTime.name(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ), + + rhoU + ( + IOobject + ( + "rhoU", + runTime.name(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + rho*U + ), + + rhoE + ( + IOobject + ( + "rhoE", + runTime.name(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + rho*(thermo.he() + 0.5*magSqr(U)) + ), + + phi + ( + IOobject + ( + "phi", + runTime.name(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + linearInterpolate(rho*U) & mesh.Sf() + ), + + inviscid + ( + max(thermo.mu()().primitiveField()) > 0 + ? false + : true + ), + + momentumTransport + ( + inviscid + ? autoPtr(nullptr) + : compressible::momentumTransportModel::New + ( + rho, + U, + phi, + thermo + ) + ), + + thermophysicalTransport + ( + inviscid + ? autoPtr(nullptr) + : fluidThermoThermophysicalTransportModel::New + ( + momentumTransport(), + thermo + ) + ), + + fluxScheme + ( + mesh.schemes().dict().lookupOrDefault("fluxScheme", "Kurganov") + ) +{ + // Read the controls + read(); + + thermo.validate(type(), "e"); + + if (momentumTransport.valid()) + { + momentumTransport->validate(); + } + + fluxPredictor(); + + if (transient()) + { + const surfaceScalarField amaxSf + ( + max(mag(aphiv_pos()), mag(aphiv_neg())) + ); + + correctCoNum(amaxSf); + } + else if (LTS) + { + Info<< "Using LTS" << endl; + + trDeltaT = tmp + ( + new volScalarField + ( + IOobject + ( + fv::localEulerDdt::rDeltaTName, + runTime.name(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar(dimless/dimTime, 1), + extrapolatedCalculatedFvPatchScalarField::typeName + ) + ); + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::solvers::shockFluid::~shockFluid() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::solvers::shockFluid::preSolve() +{ + // Read the controls + read(); + + { + const surfaceScalarField amaxSf + ( + max(mag(aphiv_pos()), mag(aphiv_neg())) + ); + + if (transient()) + { + correctCoNum(amaxSf); + } + else if (LTS) + { + setRDeltaT(amaxSf); + } + } + + fvModels().preUpdateMesh(); + + if (mesh.topoChanging()) + { + pos.clear(); + neg.clear(); + + clearTemporaryFields(); + } + + // Update the mesh for topology change, mesh to mesh mapping + mesh.update(); +} + + +void Foam::solvers::shockFluid::postCorrector() +{ + if (!inviscid && pimple.correctTransport()) + { + momentumTransport->correct(); + thermophysicalTransport->correct(); + } +} + + +void Foam::solvers::shockFluid::postSolve() +{} + + +// ************************************************************************* // diff --git a/applications/solvers/modules/shockFluid/shockFluid.H b/applications/solvers/modules/shockFluid/shockFluid.H new file mode 100644 index 0000000000..42331c8395 --- /dev/null +++ b/applications/solvers/modules/shockFluid/shockFluid.H @@ -0,0 +1,269 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2023 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 . + +Class + Foam::solvers::shockFluid + +Description + Solver module for density-based solution of compressible flow + + Based on central-upwind schemes of Kurganov and Tadmor with support for + mesh-motion and topology change. + + Reference: + \verbatim + Greenshields, C. J., Weller, H. G., Gasparini, L., + & Reese, J. M. (2010). + Implementation of semi‐discrete, non‐staggered central schemes + in a colocated, polyhedral, finite volume framework, + for high‐speed viscous flows. + International journal for numerical methods in fluids, 63(1), 1-21. + \endverbatim + +SourceFiles + shockFluid.C + +See also + Foam::solvers::fluidSolver + Foam::solvers::incompressibleFluid + +\*---------------------------------------------------------------------------*/ + +#ifndef shockFluid_H +#define shockFluid_H + +#include "fluidSolver.H" +#include "psiThermo.H" +#include "compressibleMomentumTransportModels.H" +#include "fluidThermoThermophysicalTransportModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solvers +{ + +/*---------------------------------------------------------------------------*\ + Class shockFluid Declaration +\*---------------------------------------------------------------------------*/ + +class shockFluid +: + public fluidSolver +{ + +protected: + + // Thermophysical properties + + //- Pointer to the fluid thermophysical properties + autoPtr thermo_; + + //- Reference to the fluid thermophysical properties + psiThermo& thermo; + + //- Reference to the pressure field + volScalarField& p; + + //- The continuity density field + volScalarField rho; + + //- Velocity field + volVectorField U; + + //- The momentum field + volVectorField rhoU; + + //- The energy field + volScalarField rhoE; + + + // Kinematic properties + + //- Mass-flux field + surfaceScalarField phi; + + + // Momentum transport + + bool inviscid; + + //- Pointer to the momentum transport model + autoPtr momentumTransport; + + + // Thermophysical transport + + autoPtr + thermophysicalTransport; + + + // Controls + + word fluxScheme; + + + // Cached temporary fields + + tmp pos; + tmp neg; + + tmp rho_pos; + tmp rho_neg; + + tmp rhoU_pos; + tmp rhoU_neg; + + tmp U_pos; + tmp U_neg; + + tmp p_pos; + tmp p_neg; + + tmp a_pos; + tmp a_neg; + + tmp aSf; + + tmp aphiv_pos; + tmp aphiv_neg; + + tmp muEff; + tmp tauMC; + + //- Optional LTS reciprocal time-step field + tmp trDeltaT; + + +private: + + // Private Member Functions + + //- Set rDeltaT for LTS + void setRDeltaT(const surfaceScalarField& amaxSf); + + //- Correct the cached Courant numbers + void correctCoNum(const surfaceScalarField& amaxSf); + + //- Construct the continuity equation and correct the density + void correctDensity(); + + //- Interpolate field vf according to direction dir + template + inline tmp> interpolate + ( + const VolField& vf, + const surfaceScalarField& dir, + const word& reconFieldName = word::null + ) + { + tmp> tsf + ( + fvc::interpolate + ( + vf, + dir, + "reconstruct(" + + (reconFieldName != word::null ? reconFieldName : vf.name()) + + ')' + ) + ); + + SurfaceField& sf = tsf.ref(); + + sf.rename(vf.name() + '_' + dir.name()); + + return tsf; + } + + void fluxPredictor(); + + void clearTemporaryFields(); + + +public: + + //- Runtime type information + TypeName("shockFluid"); + + + // Constructors + + //- Construct from region mesh + shockFluid(fvMesh& mesh); + + //- Disallow default bitwise copy construction + shockFluid(const shockFluid&) = delete; + + + //- Destructor + virtual ~shockFluid(); + + + // Member Functions + + //- Called at the start of the time-step, before the PIMPLE loop + virtual void preSolve(); + + //- Called at the start of the PIMPLE loop to move the mesh + virtual bool moveMesh(); + + //- Called at the start of the PIMPLE loop + virtual void prePredictor(); + + //- Construct and optionally solve the momentum equation + virtual void momentumPredictor(); + + //- Construct and solve the energy equation, + // convert to temperature + // and update thermophysical and transport properties + virtual void thermophysicalPredictor(); + + //- Construct and solve the pressure equation in the PISO loop + virtual void pressureCorrector(); + + //- Correct the momentum and thermophysical transport modelling + virtual void postCorrector(); + + //- Called after the PIMPLE loop at the end of the time-step + virtual void postSolve(); + + + // Member Operators + + //- Disallow default bitwise assignment + void operator=(const shockFluid&) = delete; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solvers +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/modules/shockFluid/thermophysicalPredictor.C b/applications/solvers/modules/shockFluid/thermophysicalPredictor.C new file mode 100644 index 0000000000..e775ae8830 --- /dev/null +++ b/applications/solvers/modules/shockFluid/thermophysicalPredictor.C @@ -0,0 +1,93 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2023 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 . + +\*---------------------------------------------------------------------------*/ + +#include "shockFluid.H" +#include "fvcDdt.H" +#include "fvmDiv.H" + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::solvers::shockFluid::thermophysicalPredictor() +{ + volScalarField& e = thermo.he(); + + const surfaceScalarField e_pos(interpolate(e, pos, thermo.T().name())); + const surfaceScalarField e_neg(interpolate(e, neg, thermo.T().name())); + + 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()) + + aSf()*(p_pos() - p_neg()) + ); + + // Make flux for pressure-work absolute + if (mesh.moving()) + { + phiEp += mesh.phi()*(a_pos()*p_pos() + a_neg()*p_neg()); + } + + solve(fvm::ddt(rhoE) + fvc::div(phiEp)); + + e = rhoE/rho - 0.5*magSqr(U); + e.correctBoundaryConditions(); + thermo.correct(); + rhoE.boundaryFieldRef() == + rho.boundaryField()* + ( + e.boundaryField() + 0.5*magSqr(U.boundaryField()) + ); + + if (!inviscid) + { + const surfaceScalarField sigmaDotU + ( + "sigmaDotU", + ( + fvc::interpolate(muEff())*mesh.magSf()*fvc::snGrad(U) + + fvc::dotInterpolate(mesh.Sf(), tauMC()) + ) + & (a_pos()*U_pos() + a_neg()*U_neg()) + ); + + solve + ( + fvm::ddt(rho, e) - fvc::ddt(rho, e) + + thermophysicalTransport->divq(e) + - fvc::div(sigmaDotU) + == + fvModels().source(rho, e) + ); + + fvConstraints().constrain(e); + + thermo.correct(); + rhoE = rho*(e + 0.5*magSqr(U)); + } +} + + +// ************************************************************************* // diff --git a/bin/rhoCentralFoam b/bin/rhoCentralFoam new file mode 100755 index 0000000000..7adb60a1d0 --- /dev/null +++ b/bin/rhoCentralFoam @@ -0,0 +1,46 @@ +#!/bin/sh +#------------------------------------------------------------------------------ +# ========= | +# \\ / F ield | OpenFOAM: The Open Source CFD Toolbox +# \\ / O peration | Website: https://openfoam.org +# \\ / A nd | Copyright (C) 2023 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 . +# +# Script +# rhoCentralFoam +# +# Description +# Script to inform the user that rhoCentralFoam has been superseded +# and replaced by the more general shockFluid solver module +# executed by the foamRun application. +# +#------------------------------------------------------------------------------ + +cat < /dev/null 2>&1 +then + echo 'gnuplot not found - skipping graph creation' >&2 + exit 1 +fi + +time=$(foamListTimes -latestTime) + +graphFile=postProcessing/sample/$time/data.xy + +gnuplot< /dev/null 2>&1 +then + echo 'gnuplot not found - skipping graph creation' >&2 + exit 1 +fi + +time=$(foamListTimes -latestTime) + +graphFile=postProcessing/sample/$time/data.xy + +gnuplot<