zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

solvers::shockFluid: Now solves for rho, U and e while conserving rho*U and rho*E

Browse Source 
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.
This commit is contained in:
Henry Weller
2023-01-24 18:17:46 +00:00
parent 4e827263ff
commit 9fadb7fccf
16 changed files with 84 additions and 145 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/files
View File

@ -1,6 +0,0 @@
mixedFixedValueSlip/mixedFixedValueSlipFvPatchFields.C
U/maxwellSlipUFvPatchVectorField.C
T/smoluchowskiJumpTFvPatchScalarField.C
rho/fixedRhoFvPatchScalarField.C
LIB = $(FOAM_LIBBIN)/librhoCentralFoam

11
applications/solvers/modules/shockFluid/derivedFvPatchFields/Make/options
View File

@ -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

26
applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.C
View File

@ -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<word>("rho", "rho")),
psiName_(dict.lookupOrDefault<word>("psi", "psi")),
muName_(dict.lookupOrDefault<word>("mu", "mu")),
tauMCName_(dict.lookupOrDefault<word>("tauMC", "tauMC")),
accommodationCoeff_(dict.lookup<scalar>("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<scalar>& ppsi =
patch().lookupPatchField<volScalarField, scalar>(psiName_);
Field<scalar> C1
const Field<scalar> C1
(
sqrt(ppsi*constant::mathematical::piByTwo)
* (2.0 - accommodationCoeff_)/accommodationCoeff_
);
Field<scalar> pnu(pmu/prho);
const Field<scalar> 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<volScalarField>(TName_);
label patchi = this->patch().index();
const label patchi = this->patch().index();
const fvPatchScalarField& pT = vsfT.boundaryField()[patchi];
Field<vector> gradpT(fvc::grad(vsfT)().boundaryField()[patchi]);
vectorField n(patch().nf());
const Field<vector> 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<volTensorField, tensor>(tauMCName_);
vectorField n(patch().nf());
const fvsPatchVectorField& divDevTauFlux =
patch().lookupPatchField<surfaceVectorField, vector>
(
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<word>(os, "rho", "rho", rhoName_);
writeEntryIfDifferent<word>(os, "psi", "psi", psiName_);
writeEntryIfDifferent<word>(os, "mu", "mu", muName_);
writeEntryIfDifferent<word>(os, "tauMC", "tauMC", tauMCName_);
writeEntry(os, "accommodationCoeff", accommodationCoeff_);
writeEntry(os, "Uwall", Uwall_);

3
applications/solvers/modules/shockFluid/derivedFvPatchFields/U/maxwellSlipUFvPatchVectorField.H
View File

@ -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_;

1
applications/solvers/modules/shockFluid/fluxPredictor.C
View File

@ -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());

56
applications/solvers/modules/shockFluid/momentumPredictor.C
View File

@ -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<fvVectorMatrix> 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();
}
}

32
applications/solvers/modules/shockFluid/shockFluid.C
View File

@ -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();

19
applications/solvers/modules/shockFluid/shockFluid.H
View File

@ -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<surfaceScalarField> aphiv_pos;
tmp<surfaceScalarField> aphiv_neg;
tmp<volScalarField> muEff;
tmp<volTensorField> tauMC;
tmp<surfaceVectorField> devTau;
//- Optional LTS reciprocal time-step field
tmp<volScalarField> trDeltaT;

51
applications/solvers/modules/shockFluid/thermophysicalPredictor.C
View File

@ -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();
}

2
tutorials/modules/shockFluid/LadenburgJet60psi/system/fvSchemes
View File

@ -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;

2
tutorials/modules/shockFluid/biconic25-55Run35/system/fvSchemes
View File

@ -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;

2
tutorials/modules/shockFluid/forwardStep/system/fvSchemes
View File

@ -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;

2
tutorials/modules/shockFluid/movingCone/system/fvSchemes
View File

@ -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;

2
tutorials/modules/shockFluid/obliqueShock/system/fvSchemes
View File

@ -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;

2
tutorials/modules/shockFluid/wedge15Ma5/system/fvSchemes
View File

@ -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;

12
tutorials/modules/shockFluid/wedge15Ma5/system/fvSolution
View File

@ -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;
}
}