From 9fadb7fccf519c8a86fcb134965b4dffa15f6dc2 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Tue, 24 Jan 2023 18:17:46 +0000 Subject: [PATCH] solvers::shockFluid: Now solves for rho, U and e while conserving rho*U and rho*E By solving for U and e rather than rhoU and rhoE the convection and stress matrices can be combined and solved together avoiding the need for Strang splitting. Conservation of rho*U and rho*E is ensured by constructing and solving the three equations in sequence, constructing each using the results of the solution of the previous equations. --- .../derivedFvPatchFields/Make/files | 6 -- .../derivedFvPatchFields/Make/options | 11 ---- .../U/maxwellSlipUFvPatchVectorField.C | 26 ++++----- .../U/maxwellSlipUFvPatchVectorField.H | 3 - .../modules/shockFluid/fluxPredictor.C | 1 + .../modules/shockFluid/momentumPredictor.C | 56 +++++++++---------- .../solvers/modules/shockFluid/shockFluid.C | 32 ++--------- .../solvers/modules/shockFluid/shockFluid.H | 19 +++---- .../shockFluid/thermophysicalPredictor.C | 51 +++++++---------- .../LadenburgJet60psi/system/fvSchemes | 2 +- .../biconic25-55Run35/system/fvSchemes | 2 +- .../shockFluid/forwardStep/system/fvSchemes | 2 +- .../shockFluid/movingCone/system/fvSchemes | 2 +- .../shockFluid/obliqueShock/system/fvSchemes | 2 +- .../shockFluid/wedge15Ma5/system/fvSchemes | 2 +- .../shockFluid/wedge15Ma5/system/fvSolution | 12 ++-- 16 files changed, 84 insertions(+), 145 deletions(-) delete mode 100644 applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/files delete mode 100644 applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/options diff --git a/applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/files b/applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/files deleted file mode 100644 index c666b09751..0000000000 --- a/applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/files +++ /dev/null @@ -1,6 +0,0 @@ -mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C -U/maxwellSlipUFvPatchVectorField.C -T/smoluchowskiJumpTFvPatchScalarField.C -rho/fixedRhoFvPatchScalarField.C - -LIB = $(FOAM_LIBBIN)/librhoCentralFoam diff --git a/applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/options b/applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/options deleted file mode 100644 index 3c16cbd9d5..0000000000 --- a/applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/options +++ /dev/null @@ -1,11 +0,0 @@ -EXE_INC = \ - -I$(LIB_SRC)/finiteVolume/lnInclude \ - -I$(LIB_SRC)/physicalProperties/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ - -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude - -LIB_LIBS = \ - -lfiniteVolume \ - -lfluidThermophysicalModels \ - -lspecie - diff --git a/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.C b/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.C index 60dbaa21cc..4a7a2a302f 100644 --- a/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.C +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.C @@ -43,7 +43,6 @@ Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField rhoName_("rho"), psiName_("psi"), muName_("mu"), - tauMCName_("tauMC"), accommodationCoeff_(1.0), Uwall_(p.size(), vector(0.0, 0.0, 0.0)), thermalCreep_(true), @@ -63,7 +62,6 @@ Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField rhoName_(dict.lookupOrDefault("rho", "rho")), psiName_(dict.lookupOrDefault("psi", "psi")), muName_(dict.lookupOrDefault("mu", "mu")), - tauMCName_(dict.lookupOrDefault("tauMC", "tauMC")), accommodationCoeff_(dict.lookup("accommodationCoeff")), Uwall_("Uwall", dict, p.size()), thermalCreep_(dict.lookupOrDefault("thermalCreep", true)), @@ -118,7 +116,6 @@ Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField rhoName_(mspvf.rhoName_), psiName_(mspvf.psiName_), muName_(mspvf.muName_), - tauMCName_(mspvf.tauMCName_), accommodationCoeff_(mspvf.accommodationCoeff_), Uwall_(mapper(mspvf.Uwall_)), thermalCreep_(mspvf.thermalCreep_), @@ -137,7 +134,6 @@ Foam::maxwellSlipUFvPatchVectorField::maxwellSlipUFvPatchVectorField rhoName_(mspvf.rhoName_), psiName_(mspvf.psiName_), muName_(mspvf.muName_), - tauMCName_(mspvf.tauMCName_), accommodationCoeff_(mspvf.accommodationCoeff_), Uwall_(mspvf.Uwall_), thermalCreep_(mspvf.thermalCreep_), @@ -200,36 +196,39 @@ void Foam::maxwellSlipUFvPatchVectorField::updateCoeffs() const fvPatchField& ppsi = patch().lookupPatchField(psiName_); - Field C1 + const Field C1 ( sqrt(ppsi*constant::mathematical::piByTwo) * (2.0 - accommodationCoeff_)/accommodationCoeff_ ); - Field pnu(pmu/prho); + const Field pnu(pmu/prho); valueFraction() = (1.0/(1.0 + patch().deltaCoeffs()*C1*pnu)); refValue() = Uwall_; + const vectorField n(patch().nf()); + if (thermalCreep_) { const volScalarField& vsfT = this->db().objectRegistry::lookupObject(TName_); - label patchi = this->patch().index(); + const label patchi = this->patch().index(); const fvPatchScalarField& pT = vsfT.boundaryField()[patchi]; - Field gradpT(fvc::grad(vsfT)().boundaryField()[patchi]); - vectorField n(patch().nf()); + const Field gradpT(fvc::grad(vsfT)().boundaryField()[patchi]); refValue() -= 3.0*pnu/(4.0*pT)*transform(I - n*n, gradpT); } if (curvature_) { - const fvPatchTensorField& ptauMC = - patch().lookupPatchField(tauMCName_); - vectorField n(patch().nf()); + const fvsPatchVectorField& divDevTauFlux = + patch().lookupPatchField + ( + IOobject::groupName("devTauCorrFlux", internalField().group()) + ); - refValue() -= C1/prho*transform(I - n*n, (n & ptauMC)); + refValue() += C1/prho*transform(I - n*n, divDevTauFlux/patch().magSf()); } mixedFixedValueSlipFvPatchVectorField::updateCoeffs(); @@ -243,7 +242,6 @@ void Foam::maxwellSlipUFvPatchVectorField::write(Ostream& os) const writeEntryIfDifferent(os, "rho", "rho", rhoName_); writeEntryIfDifferent(os, "psi", "psi", psiName_); writeEntryIfDifferent(os, "mu", "mu", muName_); - writeEntryIfDifferent(os, "tauMC", "tauMC", tauMCName_); writeEntry(os, "accommodationCoeff", accommodationCoeff_); writeEntry(os, "Uwall", Uwall_); diff --git a/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.H b/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.H index 61a88ec642..69e6b842b4 100644 --- a/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.H +++ b/applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.H @@ -66,9 +66,6 @@ class maxwellSlipUFvPatchVectorField //- Dynamic viscosity field name, default = "mu" word muName_; - //- tauMC field name, default = "tauMC" - word tauMCName_; - // Accommodation coefficient scalar accommodationCoeff_; diff --git a/applications/solvers/modules/shockFluid/fluxPredictor.C b/applications/solvers/modules/shockFluid/fluxPredictor.C index bdf95b9424..2940a044ce 100644 --- a/applications/solvers/modules/shockFluid/fluxPredictor.C +++ b/applications/solvers/modules/shockFluid/fluxPredictor.C @@ -49,6 +49,7 @@ void Foam::solvers::shockFluid::fluxPredictor() rho_pos = interpolate(rho, pos()); rho_neg = interpolate(rho, neg()); + const volVectorField rhoU(rho*U); rhoU_pos = interpolate(rhoU, pos(), U.name()); rhoU_neg = interpolate(rhoU, neg(), U.name()); diff --git a/applications/solvers/modules/shockFluid/momentumPredictor.C b/applications/solvers/modules/shockFluid/momentumPredictor.C index 4b5c920947..86e72d0f5a 100644 --- a/applications/solvers/modules/shockFluid/momentumPredictor.C +++ b/applications/solvers/modules/shockFluid/momentumPredictor.C @@ -32,49 +32,45 @@ License 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)); + // Construct the divDevTau matrix first + // so that the maxwellSlipU BC can access the explicit part + tmp divDevTau; + if (!inviscid) + { + divDevTau = momentumTransport->divDevTau(U); + } - U.ref() = rhoU()/rho(); - U.correctBoundaryConditions(); + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + fvc::div(phiUp) + == + fvModels().source(rho, U) + ); if (!inviscid) { - solve - ( - fvm::ddt(rho, U) - fvc::ddt(rho, U) - - fvm::laplacian(muEff(), U) - - fvc::div(tauMC()) - == - fvModels().source(rho, U) - ); + UEqn += divDevTau(); + } - rhoU == rho*U; - } - else - { - rhoU.boundaryFieldRef() == rho.boundaryField()*U.boundaryField(); - } + UEqn.relax(); + + fvConstraints().constrain(UEqn); + + solve(UEqn); fvConstraints().constrain(U); + K = 0.5*magSqr(U); + + if (!inviscid) + { + devTau = divDevTau->flux(); + } } diff --git a/applications/solvers/modules/shockFluid/shockFluid.C b/applications/solvers/modules/shockFluid/shockFluid.C index 1eae1bd1b5..64a76a78de 100644 --- a/applications/solvers/modules/shockFluid/shockFluid.C +++ b/applications/solvers/modules/shockFluid/shockFluid.C @@ -84,8 +84,7 @@ void Foam::solvers::shockFluid::clearTemporaryFields() aphiv_pos.clear(); aphiv_neg.clear(); - muEff.clear(); - tauMC.clear(); + devTau.clear(); } @@ -127,32 +126,6 @@ Foam::solvers::shockFluid::shockFluid(fvMesh& mesh) 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 @@ -166,6 +139,8 @@ Foam::solvers::shockFluid::shockFluid(fvMesh& mesh) linearInterpolate(rho*U) & mesh.Sf() ), + K("K", 0.5*magSqr(U)), + inviscid ( max(thermo.mu()().primitiveField()) > 0 @@ -210,6 +185,7 @@ Foam::solvers::shockFluid::shockFluid(fvMesh& mesh) if (momentumTransport.valid()) { momentumTransport->validate(); + mesh.schemes().setFluxRequired(U.name()); } fluxPredictor(); diff --git a/applications/solvers/modules/shockFluid/shockFluid.H b/applications/solvers/modules/shockFluid/shockFluid.H index 42331c8395..9bc2422b32 100644 --- a/applications/solvers/modules/shockFluid/shockFluid.H +++ b/applications/solvers/modules/shockFluid/shockFluid.H @@ -89,21 +89,19 @@ protected: //- The continuity density field volScalarField rho; - //- Velocity field - volVectorField U; - - //- The momentum field - volVectorField rhoU; - - //- The energy field - volScalarField rhoE; - // Kinematic properties + //- Velocity field + volVectorField U; + //- Mass-flux field surfaceScalarField phi; + //- Kinetic energy field + // Used in the energy equation + volScalarField K; + // Momentum transport @@ -149,8 +147,7 @@ protected: tmp aphiv_pos; tmp aphiv_neg; - tmp muEff; - tmp tauMC; + tmp devTau; //- Optional LTS reciprocal time-step field tmp trDeltaT; diff --git a/applications/solvers/modules/shockFluid/thermophysicalPredictor.C b/applications/solvers/modules/shockFluid/thermophysicalPredictor.C index e775ae8830..74e39ba182 100644 --- a/applications/solvers/modules/shockFluid/thermophysicalPredictor.C +++ b/applications/solvers/modules/shockFluid/thermophysicalPredictor.C @@ -50,43 +50,34 @@ void Foam::solvers::shockFluid::thermophysicalPredictor() 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()) - ); + fvScalarMatrix EEqn + ( + fvm::ddt(rho, e) + fvc::div(phiEp) + + fvc::ddt(rho, K) + == + fvModels().source(rho, e) + ); if (!inviscid) { - const surfaceScalarField sigmaDotU + const surfaceScalarField devTauDotU ( - "sigmaDotU", - ( - fvc::interpolate(muEff())*mesh.magSf()*fvc::snGrad(U) - + fvc::dotInterpolate(mesh.Sf(), tauMC()) - ) - & (a_pos()*U_pos() + a_neg()*U_neg()) + "devTauDotU", + devTau() & (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)); + EEqn += thermophysicalTransport->divq(e) + fvc::div(devTauDotU); } + + EEqn.relax(); + + fvConstraints().constrain(EEqn); + + EEqn.solve(); + + fvConstraints().constrain(e); + + thermo.correct(); } diff --git a/tutorials/modules/shockFluid/LadenburgJet60psi/system/fvSchemes b/tutorials/modules/shockFluid/LadenburgJet60psi/system/fvSchemes index 3f43397e5d..e9831ac15d 100644 --- a/tutorials/modules/shockFluid/LadenburgJet60psi/system/fvSchemes +++ b/tutorials/modules/shockFluid/LadenburgJet60psi/system/fvSchemes @@ -29,7 +29,6 @@ gradSchemes divSchemes { default none; - div(tauMC) Gauss linear; } laplacianSchemes @@ -40,6 +39,7 @@ laplacianSchemes interpolationSchemes { default linear; + reconstruct(rho) vanLeer; reconstruct(U) vanLeerV; reconstruct(T) vanLeer; diff --git a/tutorials/modules/shockFluid/biconic25-55Run35/system/fvSchemes b/tutorials/modules/shockFluid/biconic25-55Run35/system/fvSchemes index 01e7c3ec9e..513ddbbe5d 100644 --- a/tutorials/modules/shockFluid/biconic25-55Run35/system/fvSchemes +++ b/tutorials/modules/shockFluid/biconic25-55Run35/system/fvSchemes @@ -29,7 +29,6 @@ gradSchemes divSchemes { default none; - div(tauMC) Gauss linear; } laplacianSchemes @@ -40,6 +39,7 @@ laplacianSchemes interpolationSchemes { default linear; + reconstruct(rho) vanLeer; reconstruct(U) vanLeerV; reconstruct(T) vanLeer; diff --git a/tutorials/modules/shockFluid/forwardStep/system/fvSchemes b/tutorials/modules/shockFluid/forwardStep/system/fvSchemes index 3f43397e5d..e9831ac15d 100644 --- a/tutorials/modules/shockFluid/forwardStep/system/fvSchemes +++ b/tutorials/modules/shockFluid/forwardStep/system/fvSchemes @@ -29,7 +29,6 @@ gradSchemes divSchemes { default none; - div(tauMC) Gauss linear; } laplacianSchemes @@ -40,6 +39,7 @@ laplacianSchemes interpolationSchemes { default linear; + reconstruct(rho) vanLeer; reconstruct(U) vanLeerV; reconstruct(T) vanLeer; diff --git a/tutorials/modules/shockFluid/movingCone/system/fvSchemes b/tutorials/modules/shockFluid/movingCone/system/fvSchemes index 73becce720..9349e4351f 100644 --- a/tutorials/modules/shockFluid/movingCone/system/fvSchemes +++ b/tutorials/modules/shockFluid/movingCone/system/fvSchemes @@ -29,7 +29,6 @@ gradSchemes divSchemes { default none; - div(tauMC) Gauss linear; } laplacianSchemes @@ -42,6 +41,7 @@ laplacianSchemes interpolationSchemes { default linear; + reconstruct(rho) vanLeer; reconstruct(U) vanLeerV; reconstruct(T) vanLeer; diff --git a/tutorials/modules/shockFluid/obliqueShock/system/fvSchemes b/tutorials/modules/shockFluid/obliqueShock/system/fvSchemes index 3f43397e5d..e9831ac15d 100644 --- a/tutorials/modules/shockFluid/obliqueShock/system/fvSchemes +++ b/tutorials/modules/shockFluid/obliqueShock/system/fvSchemes @@ -29,7 +29,6 @@ gradSchemes divSchemes { default none; - div(tauMC) Gauss linear; } laplacianSchemes @@ -40,6 +39,7 @@ laplacianSchemes interpolationSchemes { default linear; + reconstruct(rho) vanLeer; reconstruct(U) vanLeerV; reconstruct(T) vanLeer; diff --git a/tutorials/modules/shockFluid/wedge15Ma5/system/fvSchemes b/tutorials/modules/shockFluid/wedge15Ma5/system/fvSchemes index 3f43397e5d..e9831ac15d 100644 --- a/tutorials/modules/shockFluid/wedge15Ma5/system/fvSchemes +++ b/tutorials/modules/shockFluid/wedge15Ma5/system/fvSchemes @@ -29,7 +29,6 @@ gradSchemes divSchemes { default none; - div(tauMC) Gauss linear; } laplacianSchemes @@ -40,6 +39,7 @@ laplacianSchemes interpolationSchemes { default linear; + reconstruct(rho) vanLeer; reconstruct(U) vanLeerV; reconstruct(T) vanLeer; diff --git a/tutorials/modules/shockFluid/wedge15Ma5/system/fvSolution b/tutorials/modules/shockFluid/wedge15Ma5/system/fvSolution index 98cba61940..2ffe60c678 100644 --- a/tutorials/modules/shockFluid/wedge15Ma5/system/fvSolution +++ b/tutorials/modules/shockFluid/wedge15Ma5/system/fvSolution @@ -16,7 +16,7 @@ FoamFile solvers { - "(rho|rhoU|rhoE).*" + "rho.*" { solver diagonal; } @@ -25,16 +25,16 @@ solvers { solver smoothSolver; smoother GaussSeidel; - nSweeps 2; - tolerance 1e-09; - relTol 0.01; + + tolerance 1e-9; + relTol 0; } - "h.*" + "e.*" { $U; + tolerance 1e-10; - relTol 0; } }