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

solvers::multiphase: Improved CorrectPhi handling for compressible multiphase flows

Browse Source 
The mixture compressibility/density is now included in CorrectPhi for
compressible mixtures, consistent with the compressibility handling in the
pressure equation.  This improves consistency, robustness and convergence of the
pcorr equation.
This commit is contained in:
Henry Weller
2023-03-29 15:59:13 +01:00
parent f67445212d
commit 113d07862c
20 changed files with 249 additions and 53 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
applications/solvers/modules/VoFSolver/VoFSolver.C
View File

@ -25,8 +25,9 @@ License
#include "VoFSolver.H"
#include "localEulerDdtScheme.H"
#include "CorrectPhi.H"
#include "geometricZeroField.H"
#include "linear.H"
#include "fvcDiv.H"
#include "fvcMeshPhi.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

9
applications/solvers/modules/VoFSolver/VoFSolver.H
View File

@ -169,9 +169,16 @@ protected:
//- Return the pressure reference
virtual const Foam::pressureReference& pressureReference() const = 0;
//- Is the flow incompressible?
virtual bool incompressible() const = 0;
//- Is the flow divergent?
// i.e. compressible or include phase-fraction sources
virtual bool divergent() = 0;
virtual bool divergent() const = 0;
//- Return the mixture compressibility/density
// Used by CorrectPhi for compressible mixtures following mesh change
virtual tmp<volScalarField> psiByRho() const = 0;
//- Return the momentum equation stress term
virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U) = 0;

43
applications/solvers/modules/VoFSolver/moveMesh.C
View File

@ -65,29 +65,44 @@ void Foam::solvers::VoFSolver::moveMesh()
correctUphiBCs(U_, phi, true);
if (divU.valid())
if (incompressible())
{
CorrectPhi
(
phi,
U,
p_rgh,
surfaceScalarField("rAUf", fvc::interpolate(rAU())),
divU(),
pressureReference(),
pimple
);
if (divU.valid())
{
CorrectPhi
(
phi,
U,
p_rgh,
surfaceScalarField("rAUf", fvc::interpolate(rAU())),
divU(),
pressureReference(),
pimple
);
}
else
{
CorrectPhi
(
phi,
U,
p_rgh,
surfaceScalarField("rAUf", fvc::interpolate(rAU())),
geometricZeroField(),
pressureReference(),
pimple
);
}
}
else
{
CorrectPhi
(
phi,
U,
p_rgh,
psiByRho(),
surfaceScalarField("rAUf", fvc::interpolate(rAU())),
geometricZeroField(),
pressureReference(),
divU(),
pimple
);
}

2
applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.C
View File

@ -24,9 +24,9 @@ License
\*---------------------------------------------------------------------------*/
#include "compressibleMultiphaseVoF.H"
#include "CorrectPhi.H"
#include "geometricZeroField.H"
#include "fvcDdt.H"
#include "fvcDiv.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

20
applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoF.H
View File

@ -128,10 +128,24 @@ protected:
return pressureReference_;
}
//- Compressible flow is divergent
virtual bool divergent()
//- The flow is incompressible if all phases are incompressible
virtual bool incompressible() const
{
return true;
return mixture.incompressible();
}
//- The flow is divergent if it is not incompressible
// Mass sources are not currently supported
virtual bool divergent() const
{
return !incompressible();
}
//- Return the mixture compressibility/density
// Used by CorrectPhi for compressible mixtures following mesh change
virtual tmp<volScalarField> psiByRho() const
{
return mixture.psiByRho();
}
//- Return the momentum equation stress term

33
applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.C
View File

@ -67,6 +67,20 @@ Foam::compressibleMultiphaseVoFMixture::compressibleMultiphaseVoFMixture
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::compressibleMultiphaseVoFMixture::incompressible() const
{
bool incompressible = true;
forAll(phases_, phasei)
{
incompressible =
incompressible && phases_[phasei].thermo().incompressible();
}
return incompressible;
}
Foam::tmp<Foam::volScalarField>
Foam::compressibleMultiphaseVoFMixture::nu() const
{
@ -103,6 +117,25 @@ Foam::tmp<Foam::scalarField> Foam::compressibleMultiphaseVoFMixture::nu
}
Foam::tmp<Foam::volScalarField>
Foam::compressibleMultiphaseVoFMixture::psiByRho() const
{
tmp<volScalarField> tpsiByRho
(
phases_[0]*phases_[0].thermo().psi()/phases_[0].thermo().rho()
);
for (label phasei=1; phasei<phases_.size(); phasei++)
{
tpsiByRho.ref() +=
phases_[phasei]*phases_[phasei].thermo().psi()
/phases_[phasei].thermo().rho();
}
return tpsiByRho;
}
Foam::tmp<Foam::volScalarField> Foam::compressibleMultiphaseVoFMixture::alphaEff
(
const volScalarField& nut

6
applications/solvers/modules/compressibleMultiphaseVoF/compressibleMultiphaseVoFMixture/compressibleMultiphaseVoFMixture.H
View File

@ -104,6 +104,9 @@ public:
return rho_;
}
//- Return true if all phases are incompressible
bool incompressible() const;
//- Return the kinematic laminar viscosity
virtual tmp<volScalarField> nu() const;
@ -113,6 +116,9 @@ public:
//- Return the face-interpolated dynamic laminar viscosity
tmp<surfaceScalarField> nuf() const;
//- Return the mixture compressibility/density
virtual tmp<volScalarField> psiByRho() const;
//- Correct the thermodynamics of each phase
virtual void correctThermo();

11
applications/solvers/modules/compressibleVoF/compressibleTwoPhaseVoFMixture/compressibleTwoPhaseVoFMixture.C
View File

@ -203,6 +203,17 @@ Foam::tmp<Foam::scalarField> Foam::compressibleTwoPhaseVoFMixture::nu
}
Foam::tmp<Foam::volScalarField>
Foam::compressibleTwoPhaseVoFMixture::psiByRho() const
{
return
(
alpha1()*thermo1_->psi()/thermo1_->rho()
+ alpha2()*thermo2_->psi()/thermo2_->rho()
);
}
bool Foam::compressibleTwoPhaseVoFMixture::read()
{
if (twoPhaseVoFMixture::read())

11
applications/solvers/modules/compressibleVoF/compressibleTwoPhaseVoFMixture/compressibleTwoPhaseVoFMixture.H
View File

@ -161,12 +161,23 @@ public:
return thermo2_->rho();
}
//- The fluid is incompressible if both phases are incompressible
bool incompressible() const
{
return
thermo1_->incompressible()
&& thermo2_->incompressible();
}
//- Return mixture density [kg/m^3]
virtual const volScalarField& rho() const
{
return rho_;
}
//- Return the mixture compressibility/density
virtual tmp<volScalarField> psiByRho() const;
//- Correct the thermodynamics of each phase
virtual void correctThermo();

21
applications/solvers/modules/compressibleVoF/compressibleVoF.H
View File

@ -148,10 +148,25 @@ protected:
return pressureReference_;
}
//- Compressible flow is divergent
virtual bool divergent()
//- The fluid is incompressible if both phases are incompressible
virtual bool incompressible() const
{
return true;
return mixture.incompressible();
}
//- Compressible flow is divergent
virtual bool divergent() const
{
return
!incompressible()
|| fvModels().addsSupToField(alpha1.name());
}
//- Return the mixture compressibility/density
// Used by CorrectPhi for compressible mixtures following mesh change
virtual tmp<volScalarField> psiByRho() const
{
return mixture.psiByRho();
}
//- Calculate the alpha equation sources

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

@ -111,13 +111,26 @@ protected:
return pressureReference_;
}
//- Is the flow divergent?
// i.e. includes phase-fraction sources
virtual bool divergent()
//- The flow is incompressible
virtual bool incompressible() const
{
return true;
}
//- The flow is not incompressible and hence not divergent
// Mass sources are not currently supported
virtual bool divergent() const
{
return false;
}
//- Return the mixture compressibility/density
// Not required for incompressible fluids
virtual tmp<volScalarField> psiByRho() const
{
return tmp<volScalarField>(nullptr);
}
//- Return the momentum equation stress term
virtual tmp<fvVectorMatrix> divDevTau(volVectorField& U)
{

15
applications/solvers/modules/incompressibleVoF/incompressibleVoF.H
View File

@ -112,13 +112,26 @@ protected:
return pressureReference_;
}
//- The flow is incompressible
virtual bool incompressible() const
{
return true;
}
//- Is the flow divergent?
// i.e. includes phase-fraction sources
virtual bool divergent()
virtual bool divergent() const
{
return fvModels().addsSupToField(alpha1.name());
}
//- Return the mixture compressibility/density
// Not required for incompressible fluids
virtual tmp<volScalarField> psiByRho() const
{
return tmp<volScalarField>(nullptr);
}
//- Calculate the alpha equation sources
virtual void alphaSuSp
(

1
applications/solvers/modules/isothermalFluid/moveMesh.C
View File

@ -61,7 +61,6 @@ void Foam::solvers::isothermalFluid::moveMesh()
(
phi,
buoyancy.valid() ? p_rgh : p,
rho,
thermo.psi(),
dimensionedScalar("rAUf", dimTime, 1),
divrhoU(),

1
applications/solvers/modules/multiphaseEuler/multiphaseEuler/moveMesh.C
View File

@ -40,6 +40,7 @@ void Foam::solvers::multiphaseEuler::moveMesh()
if
(
(correctPhi || mesh.topoChanged())
// && divergent()
&& !divU.valid()
)
{

70
applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C
View File

@ -774,17 +774,65 @@ void Foam::phaseSystem::correctPhi
phi_ += alphafs[phasei]*(mesh_.Sf() & phase.UfRef());
}
CorrectPhi
(
phi_,
movingPhases()[0].U(),
p_rgh,
// surfaceScalarField("rAUf", fvc::interpolate(rAU())),
dimensionedScalar(dimTime/dimDensity, 1),
divU(),
pressureReference,
pimple
);
if (incompressible())
{
if (divU.valid())
{
CorrectPhi
(
phi_,
movingPhases()[0].U(),
p_rgh,
dimensionedScalar(dimTime/dimDensity, 1),
divU(),
pressureReference,
pimple
);
}
else
{
CorrectPhi
(
phi_,
movingPhases()[0].U(),
p_rgh,
dimensionedScalar(dimTime/dimDensity, 1),
geometricZeroField(),
pressureReference,
pimple
);
}
}
else
{
volScalarField psi
(
volScalarField::New
(
"psi",
mesh_,
dimensionedScalar(dimless/dimPressure, 0)
)
);
forAll(phases(), phasei)
{
phaseModel& phase = phases()[phasei];
const volScalarField& alpha = phase;
psi += alpha*phase.thermo().psi()/phase.thermo().rho();
}
CorrectPhi
(
phi_,
p_rgh,
psi,
dimensionedScalar(dimTime/dimDensity, 1),
divU(),
pimple
);
}
// Make the flux relative to the mesh motion
fvc::makeRelative(phi_, movingPhases()[0].U());

3
applications/solvers/modules/multiphaseVoFSolver/multiphaseVoFSolver.C
View File

@ -25,8 +25,7 @@ License
#include "multiphaseVoFSolver.H"
#include "localEulerDdtScheme.H"
#include "CorrectPhi.H"
#include "geometricZeroField.H"
#include "fvcAverage.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

3
applications/solvers/modules/twoPhaseVoFSolver/twoPhaseVoFSolver.C
View File

@ -25,8 +25,7 @@ License
#include "twoPhaseVoFSolver.H"
#include "localEulerDdtScheme.H"
#include "CorrectPhi.H"
#include "geometricZeroField.H"
#include "fvcAverage.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

3
src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2015-2022 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -118,7 +118,6 @@ void Foam::CorrectPhi
(
surfaceScalarField& phi,
const volScalarField& p,
const volScalarField& rho,
const volScalarField& psi,
const RAUfType& rAUf,
const DivRhoUType& divRhoU,

3
src/finiteVolume/cfdTools/general/CorrectPhi/CorrectPhi.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2015-2021 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -87,7 +87,6 @@ namespace Foam
(
surfaceScalarField& phi,
const volScalarField& p,
const volScalarField& rho,
const volScalarField& psi,
const RAUfType& rAUf,
const DivRhoUType& divRhoU,

23
tutorials/modules/multiphaseEuler/mixerVessel2D/system/fvSolution
View File

@ -25,9 +25,10 @@ solvers
p_rgh
{
solver GAMG;
smoother DIC;
tolerance 1e-7;
relTol 0.01;
smoother GaussSeidel;
}
p_rghFinal
@ -36,20 +37,30 @@ solvers
preconditioner
{
preconditioner GAMG;
smoother GaussSeidel;
tolerance 1e-7;
relTol 0;
nVcycles 2;
smoother GaussSeidel;
}
tolerance 1e-7;
relTol 0;
maxIter 30;
maxIter 50;
}
"pcorr.*"
{
$p_rghFinal;
tolerance 1e-5;
tolerance 0.001;
relTol 0;
}
MeshPhi
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-2;
relTol 0;
}
@ -57,6 +68,7 @@ solvers
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-5;
relTol 0;
minIter 1;
@ -66,6 +78,7 @@ solvers
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-8;
relTol 0;
minIter 1;
@ -79,7 +92,7 @@ PIMPLE
nCorrectors 3;
nNonOrthogonalCorrectors 0;
correctPhi no;
correctPhi yes;
correctMeshPhi no;
}