From 0be1e89204598e5161dce56e5b7861b9921bff58 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Wed, 8 Apr 2015 12:19:23 +0100
Subject: [PATCH] twoPhaseEulerFoam: Interpolate lift, wall-lubrication and
turbulent dispersion forces Reduces or eliminates staggering patterns due to
cell-force imbalances Resolves bug-report
http://www.openfoam.org/mantisbt/view.php?id=1363
---
.../multiphase/twoPhaseEulerFoam/UEqns.H | 9 -
.../turbulentDispersionModels/Burns/Burns.C | 4 +-
.../turbulentDispersionModels/Burns/Burns.H | 6 +-
.../turbulentDispersionModels/Gosman/Gosman.C | 4 +-
.../turbulentDispersionModels/Gosman/Gosman.H | 6 +-
.../LopezDeBertodano/LopezDeBertodano.C | 5 +-
.../LopezDeBertodano/LopezDeBertodano.H | 6 +-
.../constantTurbulentDispersionCoefficient.C | 4 +-
.../constantTurbulentDispersionCoefficient.H | 6 +-
.../noTurbulentDispersion.C | 54 ++--
.../noTurbulentDispersion.H | 10 +-
.../turbulentDispersionModel.C | 6 +-
.../turbulentDispersionModel.H | 13 +-
.../multiphase/twoPhaseEulerFoam/pEqn.H | 267 ------------------
.../multiphase/twoPhaseEulerFoam/pEqnElim.H | 98 ++++---
.../BlendedInterfacialModel.C | 55 ++++
.../BlendedInterfacialModel.H | 7 +-
.../twoPhaseSystem/twoPhaseSystem.C | 7 +
.../twoPhaseSystem/twoPhaseSystem.H | 6 +-
19 files changed, 190 insertions(+), 383 deletions(-)
delete mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
index bfd5faf8011..1cee6d04813 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H
@@ -9,9 +9,6 @@ volScalarField dragCoeff(fluid.dragCoeff());
{
volScalarField virtualMassCoeff(fluid.virtualMassCoeff());
- volVectorField liftForce(fluid.liftForce());
- volVectorField wallLubricationForce(fluid.wallLubricationForce());
- volVectorField turbulentDispersionForce(fluid.turbulentDispersionForce());
{
U1Eqn =
@@ -21,9 +18,6 @@ volScalarField dragCoeff(fluid.dragCoeff());
+ mrfZones(alpha1*rho1 + virtualMassCoeff, U1)
+ phase1.turbulence().divDevRhoReff(U1)
==
- - liftForce
- - wallLubricationForce
- - turbulentDispersionForce
- virtualMassCoeff
*(
fvm::ddt(U1)
@@ -46,9 +40,6 @@ volScalarField dragCoeff(fluid.dragCoeff());
+ mrfZones(alpha2*rho2 + virtualMassCoeff, U2)
+ phase2.turbulence().divDevRhoReff(U2)
==
- liftForce
- + wallLubricationForce
- + turbulentDispersionForce
- virtualMassCoeff
*(
fvm::ddt(U2)
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C
index 1de7c87de59..178326dd37c 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -70,7 +70,7 @@ Foam::turbulentDispersionModels::Burns::~Burns()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
-Foam::turbulentDispersionModels::Burns::Fprime() const
+Foam::turbulentDispersionModels::Burns::D() const
{
const fvMesh& mesh(pair_.phase1().mesh());
const dragModel&
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.H
index 0d414bd5f3d..9a0491ebe9e 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -103,9 +103,9 @@ public:
// Member Functions
- //- Turbulent dispersion force coefficient
+ //- Turbulent diffusivity
// multiplying the gradient of the phase-fraction
- virtual tmp<volScalarField> Fprime() const;
+ virtual tmp<volScalarField> D() const;
};
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.C
index 005ef7c1a59..a40dbf07703 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -69,7 +69,7 @@ Foam::turbulentDispersionModels::Gosman::~Gosman()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
-Foam::turbulentDispersionModels::Gosman::Fprime() const
+Foam::turbulentDispersionModels::Gosman::D() const
{
const fvMesh& mesh(pair_.phase1().mesh());
const dragModel&
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.H
index 7de9411f56c..49c48e6b122 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -92,9 +92,9 @@ public:
// Member Functions
- //- Turbulent dispersion force coefficient
+ //- Turbulent diffusivity
// multiplying the gradient of the phase-fraction
- virtual tmp<volScalarField> Fprime() const;
+ virtual tmp<volScalarField> D() const;
};
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.C
index 419c25ca0de..52adccb9b44 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -67,7 +67,7 @@ Foam::turbulentDispersionModels::LopezDeBertodano::~LopezDeBertodano()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::volScalarField>
-Foam::turbulentDispersionModels::LopezDeBertodano::Fprime() const
+Foam::turbulentDispersionModels::LopezDeBertodano::D() const
{
return
Ctd_
@@ -75,4 +75,5 @@ Foam::turbulentDispersionModels::LopezDeBertodano::Fprime() const
*pair_.continuous().turbulence().k();
}
+
// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.H
index 6f83bf0389d..ce3c9996a06 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -98,9 +98,9 @@ public:
// Member Functions
- //- Turbulent dispersion force coefficient
+ //- Turbulent diffusivity
// multiplying the gradient of the phase-fraction
- virtual tmp<volScalarField> Fprime() const;
+ virtual tmp<volScalarField> D() const;
};
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.C
index e1754a72ce0..67566b9f885 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -70,7 +70,7 @@ Foam::turbulentDispersionModels::constantTurbulentDispersionCoefficient::
Foam::tmp<Foam::volScalarField>
Foam::turbulentDispersionModels::constantTurbulentDispersionCoefficient::
-Fprime() const
+D() const
{
return
Ctd_
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.H
index 0a4146fc15f..eb6385670cd 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -83,9 +83,9 @@ public:
// Member Functions
- //- Turbulent dispersion force coefficient
+ //- Turbulent diffusivity
// multiplying the gradient of the phase-fraction
- virtual tmp<volScalarField> Fprime() const;
+ virtual tmp<volScalarField> D() const;
};
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.C
index 44a4c3504d0..446abe98149 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -65,43 +65,45 @@ Foam::turbulentDispersionModels::noTurbulentDispersion::
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-Foam::tmp<Foam::volVectorField>
-Foam::turbulentDispersionModels::noTurbulentDispersion::F() const
+Foam::tmp<Foam::volScalarField>
+Foam::turbulentDispersionModels::noTurbulentDispersion::D() const
{
const fvMesh& mesh(this->pair_.phase1().mesh());
- return
- tmp<volVectorField>
+ return tmp<volScalarField>
+ (
+ new volScalarField
(
- new volVectorField
+ IOobject
(
- IOobject
- (
- "zero",
- mesh.time().timeName(),
- mesh
- ),
+ "zero",
+ mesh.time().timeName(),
mesh,
- dimensionedVector
- (
- "zero",
- dimensionSet(1, -2, -2, 0, 0),
- vector::zero
- )
+ IOobject::NO_READ,
+ IOobject::NO_WRITE,
+ false
+ ),
+ mesh,
+ dimensionedScalar
+ (
+ "zero",
+ dimensionSet(1, -2, 1, 0, 0),
+ 0
)
- );
+ )
+ );
}
-Foam::tmp<Foam::volScalarField>
-Foam::turbulentDispersionModels::noTurbulentDispersion::Fprime() const
+Foam::tmp<Foam::volVectorField>
+Foam::turbulentDispersionModels::noTurbulentDispersion::F() const
{
const fvMesh& mesh(this->pair_.phase1().mesh());
return
- tmp<volScalarField>
+ tmp<volVectorField>
(
- new volScalarField
+ new volVectorField
(
IOobject
(
@@ -110,11 +112,11 @@ Foam::turbulentDispersionModels::noTurbulentDispersion::Fprime() const
mesh
),
mesh,
- dimensionedScalar
+ dimensionedVector
(
"zero",
- dimensionSet(1, -2, 1, 0, 0),
- 0
+ dimensionSet(1, -2, -2, 0, 0),
+ vector::zero
)
)
);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.H
index 2153ff6866f..8bca7fd04bb 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -76,12 +76,12 @@ public:
// Member Functions
+ //- Turbulent diffusivity
+ // multiplying the gradient of the phase-fraction
+ virtual tmp<volScalarField> D() const;
+
//- Turbulent dispersion force
virtual tmp<volVectorField> F() const;
-
- //- Turbulent dispersion force coefficient
- // multiplying the gradient of the phase-fraction
- virtual tmp<volScalarField> Fprime() const;
};
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C
index d4199b440f7..d5d805c0851 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -35,6 +35,7 @@ namespace Foam
defineRunTimeSelectionTable(turbulentDispersionModel, dictionary);
}
+const Foam::dimensionSet Foam::turbulentDispersionModel::dimD(1, -1, -2, 0, 0);
const Foam::dimensionSet Foam::turbulentDispersionModel::dimF(1, -2, -2, 0, 0);
@@ -61,8 +62,7 @@ Foam::turbulentDispersionModel::~turbulentDispersionModel()
Foam::tmp<Foam::volVectorField>
Foam::turbulentDispersionModel::F() const
{
- return Fprime()*fvc::grad(pair_.dispersed());
-
+ return D()*fvc::grad(pair_.dispersed());
}
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H
index 5862452544a..932fdb3801c 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2014-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -82,6 +82,9 @@ public:
// Static data members
+ //- Diffusivity dimensions
+ static const dimensionSet dimD;
+
//- Force dimensions
static const dimensionSet dimF;
@@ -111,12 +114,12 @@ public:
// Member Functions
+ //- Turbulent diffusivity
+ // multiplying the gradient of the phase-fraction
+ virtual tmp<volScalarField> D() const = 0;
+
//- Turbulent dispersion force
virtual tmp<volVectorField> F() const;
-
- //- Turbulent dispersion force coefficient
- // multiplying the gradient of the phase-fraction
- virtual tmp<volScalarField> Fprime() const = 0;
};
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
deleted file mode 100644
index 9e015406d18..00000000000
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H
+++ /dev/null
@@ -1,267 +0,0 @@
-// --- Pressure corrector loop
-while (pimple.correct())
-{
- surfaceScalarField alpha1f("alpha1f", fvc::interpolate(alpha1));
- surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f);
-
- rAU1 = 1.0/U1Eqn.A();
- rAU2 = 1.0/U2Eqn.A();
-
- surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rho1*rAU1));
- surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rho2*rAU2));
-
- volVectorField HbyA1
- (
- IOobject::groupName("HbyA", phase1.name()),
- U1
- );
- HbyA1 = rAU1*U1Eqn.H();
-
- volVectorField HbyA2
- (
- IOobject::groupName("HbyA", phase2.name()),
- U2
- );
- HbyA2 = rAU2*U2Eqn.H();
-
- mrfZones.makeAbsolute(phi1.oldTime());
- mrfZones.makeAbsolute(phi1);
- mrfZones.makeAbsolute(phi2.oldTime());
- mrfZones.makeAbsolute(phi2);
-
- // Phase-1 pressure flux (e.g. due to particle-particle pressure)
- surfaceScalarField phiF1
- (
- "phiF1",
- fvc::interpolate(rAU1*phase1.turbulence().pPrime())
- *fvc::snGrad(alpha1)*mesh.magSf()
- );
-
- // Phase-2 pressure flux (e.g. due to particle-particle pressure)
- surfaceScalarField phiF2
- (
- "phiF2",
- fvc::interpolate(rAU2*phase2.turbulence().pPrime())
- *fvc::snGrad(alpha2)*mesh.magSf()
- );
-
- volScalarField rho("rho", fluid.rho());
- surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
-
- // Add the phase-1 buoyancy force
- phiF1 +=
- rAlphaAU1f
- *(
- ghSnGradRho
- - alpha2f*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
- )/fvc::interpolate(rho1);
-
- // Add the phase-2 buoyancy force
- phiF2 +=
- rAlphaAU2f
- *(
- ghSnGradRho
- - alpha1f*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
- )/fvc::interpolate(rho2);
-
- // Phase-1 predicted flux
- surfaceScalarField phiHbyA1
- (
- IOobject::groupName("phiHbyA", phase1.name()),
- (fvc::interpolate(HbyA1) & mesh.Sf())
- + rAlphaAU1f*fvc::ddtCorr(U1, phi1)
- + fvc::interpolate(rAU1*dragCoeff)*phi2
- - phiF1
- );
-
- // Phase-2 predicted flux
- surfaceScalarField phiHbyA2
- (
- IOobject::groupName("phiHbyA", phase2.name()),
- (fvc::interpolate(HbyA2) & mesh.Sf())
- + rAlphaAU2f*fvc::ddtCorr(U2, phi2)
- + fvc::interpolate(rAU2*dragCoeff)*phi1
- - phiF2
- );
-
- mrfZones.makeRelative(phiHbyA1);
- mrfZones.makeRelative(phiHbyA2);
- mrfZones.makeRelative(phi1.oldTime());
- mrfZones.makeRelative(phi1);
- mrfZones.makeRelative(phi2.oldTime());
- mrfZones.makeRelative(phi2);
-
- surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
-
- HbyA1 += rAU1*dragCoeff*U2;
- HbyA2 += rAU2*dragCoeff*U1;
-
- surfaceScalarField rAUf
- (
- "rAUf",
- mag
- (
- alpha1f*rAlphaAU1f/fvc::interpolate(rho1)
- + alpha2f*rAlphaAU2f/fvc::interpolate(rho2)
- )
- );
-
- // Update the fixedFluxPressure BCs to ensure flux consistency
- setSnGrad<fixedFluxPressureFvPatchScalarField>
- (
- p_rgh.boundaryField(),
- (
- phiHbyA.boundaryField()
- - mrfZones.relative
- (
- alpha1f.boundaryField()
- *(mesh.Sf().boundaryField() & U1.boundaryField())
- + alpha2f.boundaryField()
- *(mesh.Sf().boundaryField() & U2.boundaryField())
- )
- )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
- );
-
- tmp<fvScalarMatrix> pEqnComp1;
- tmp<fvScalarMatrix> pEqnComp2;
-
- if (pimple.transonic())
- {
- surfaceScalarField phid1
- (
- IOobject::groupName("phid", phase1.name()),
- fvc::interpolate(psi1)*phi1
- );
- surfaceScalarField phid2
- (
- IOobject::groupName("phid", phase2.name()),
- fvc::interpolate(psi2)*phi2
- );
-
- pEqnComp1 =
- (
- contErr1
- - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
- )/rho1
- + (alpha1/rho1)*correction
- (
- psi1*fvm::ddt(p_rgh)
- + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
- );
- deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
- pEqnComp1().relax();
-
- pEqnComp2 =
- (
- contErr2
- - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
- )/rho2
- + (alpha2/rho2)*correction
- (
- psi2*fvm::ddt(p_rgh)
- + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
- );
- deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
- pEqnComp2().relax();
- }
- else
- {
- pEqnComp1 =
- (
- contErr1
- - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
- )/rho1
- + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
-
- pEqnComp2 =
- (
- contErr2
- - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
- )/rho2
- + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
- }
-
- // Cache p prior to solve for density update
- volScalarField p_rgh_0(p_rgh);
-
- while (pimple.correctNonOrthogonal())
- {
- fvScalarMatrix pEqnIncomp
- (
- fvc::div(phiHbyA)
- - fvm::laplacian(rAUf, p_rgh)
- );
-
- solve
- (
- pEqnComp1() + pEqnComp2() + pEqnIncomp,
- mesh.solver(p_rgh.select(pimple.finalInnerIter()))
- );
-
- if (pimple.finalNonOrthogonalIter())
- {
- surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
-
- phi1.boundaryField() ==
- mrfZones.relative
- (
- mesh.Sf().boundaryField() & U1.boundaryField()
- );
- phi1 = phiHbyA1 + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1);
-
- phi2.boundaryField() ==
- mrfZones.relative
- (
- mesh.Sf().boundaryField() & U2.boundaryField()
- );
- phi2 = phiHbyA2 + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2);
-
- phi = alpha1f*phi1 + alpha2f*phi2;
-
- fluid.dgdt() =
- (
- alpha1*(pEqnComp2 & p_rgh)
- - alpha2*(pEqnComp1 & p_rgh)
- );
-
- p_rgh.relax();
-
- p = max(p_rgh + rho*gh, pMin);
- p_rgh = p - rho*gh;
-
- mSfGradp = pEqnIncomp.flux()/rAUf;
-
- U1 = HbyA1
- + fvc::reconstruct
- (
- rAlphaAU1f*mSfGradp/fvc::interpolate(rho1)
- - phiF1
- );
- U1.correctBoundaryConditions();
- fvOptions.correct(U1);
-
- U2 = HbyA2
- + fvc::reconstruct
- (
- rAlphaAU2f*mSfGradp/fvc::interpolate(rho2)
- - phiF2
- );
- U2.correctBoundaryConditions();
- fvOptions.correct(U2);
-
- U = fluid.U();
- }
- }
-
- // Update densities from change in p
- rho1 += psi1*(p_rgh - p_rgh_0);
- rho2 += psi2*(p_rgh - p_rgh_0);
-
- K1 = 0.5*magSqr(U1);
- K2 = 0.5*magSqr(U2);
-
- if (thermo1.dpdt() || thermo2.dpdt())
- {
- dpdt = fvc::ddt(p);
- }
-}
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqnElim.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqnElim.H
index 638dc11b868..edac7acd49f 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqnElim.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pEqnElim.H
@@ -4,8 +4,8 @@ surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f);
rAU1 = 1.0/U1Eqn.A();
rAU2 = 1.0/U2Eqn.A();
-surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rho1*rAU1));
-surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rho2*rAU2));
+surfaceScalarField alpharAU1f(fvc::interpolate(alpha1*rAU1));
+surfaceScalarField alpharAU2f(fvc::interpolate(alpha2*rAU2));
// --- Pressure corrector loop
while (pimple.correct())
@@ -24,21 +24,43 @@ while (pimple.correct())
);
HbyA2 = rAU2*U2Eqn.H();
- // Phase pressure flux (e.g. due to particle-particle pressure)
- surfaceScalarField phiF1
- (
- "phiF1",
- fvc::interpolate(rAU1*phase1.turbulence().pPrime())
- *fvc::snGrad(alpha1)*mesh.magSf()
- );
+ // Face force fluxes
+ tmp<surfaceScalarField> phiF1;
+ tmp<surfaceScalarField> phiF2;
- // Phase-2 pressure flux (e.g. due to particle-particle pressure)
- surfaceScalarField phiF2
- (
- "phiF2",
- fvc::interpolate(rAU2*phase2.turbulence().pPrime())
- *fvc::snGrad(alpha2)*mesh.magSf()
- );
+ // Turbulent diffusion, particle-pressure lift and wall-lubrication fluxes
+ {
+ volScalarField turbulentDiffusivity(fluid.turbulentDiffusivity());
+ volVectorField liftForce(fluid.liftForce());
+ volVectorField wallLubricationForce(fluid.wallLubricationForce());
+ surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
+
+ phiF1 =
+ (
+ fvc::interpolate
+ (
+ rAU1*(turbulentDiffusivity + phase1.turbulence().pPrime())
+ )*snGradAlpha1
+
+ + (
+ fvc::interpolate(rAU1*(wallLubricationForce + liftForce))
+ & mesh.Sf()
+ )
+ );
+
+ phiF2 =
+ (
+ - fvc::interpolate
+ (
+ rAU2*(turbulentDiffusivity + phase2.turbulence().pPrime())
+ )*snGradAlpha1
+
+ - (
+ fvc::interpolate(rAU2*(wallLubricationForce + liftForce))
+ & mesh.Sf()
+ )
+ );
+ }
// Mean density for buoyancy force and p_rgh -> p
volScalarField rho("rho", fluid.rho());
@@ -47,19 +69,19 @@ while (pimple.correct())
{
surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
- phiF1 +=
- rAlphaAU1f
+ phiF1() +=
+ alpharAU1f
*(
ghSnGradRho
- alpha2f*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
- )/fvc::interpolate(rho1);
+ );
- phiF2 +=
- rAlphaAU2f
+ phiF2() +=
+ alpharAU2f
*(
ghSnGradRho
- alpha1f*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
- )/fvc::interpolate(rho2);
+ );
}
// ddtPhiCorr filter -- only apply in pure(ish) phases
@@ -86,12 +108,12 @@ while (pimple.correct())
(
IOobject::groupName("phiHbyA", phase1.name()),
(fvc::interpolate(HbyA1) & mesh.Sf())
- + phiCorrCoeff1*rAlphaAU1f
+ + phiCorrCoeff1*fvc::interpolate(alpha1.oldTime()*rho1.oldTime()*rAU1)
*(
mrfZones.absolute(phi1.oldTime())
- (fvc::interpolate(U1.oldTime()) & mesh.Sf())
)/runTime.deltaT()
- - phiF1
+ - phiF1()
);
// Phase-2 predicted flux
@@ -99,12 +121,12 @@ while (pimple.correct())
(
IOobject::groupName("phiHbyA", phase2.name()),
(fvc::interpolate(HbyA2) & mesh.Sf())
- + phiCorrCoeff2*rAlphaAU2f
+ + phiCorrCoeff2*fvc::interpolate(alpha2.oldTime()*rho2.oldTime()*rAU2)
*(
mrfZones.absolute(phi2.oldTime())
- (fvc::interpolate(U2.oldTime()) & mesh.Sf())
)/runTime.deltaT()
- - phiF2
+ - phiF2()
);
// Face-drag coefficients
@@ -125,11 +147,7 @@ while (pimple.correct())
surfaceScalarField rAUf
(
"rAUf",
- mag
- (
- alpha1f*rAlphaAU1f/fvc::interpolate(rho1)
- + alpha2f*rAlphaAU2f/fvc::interpolate(rho2)
- )
+ mag(alpha1f*alpharAU1f + alpha2f*alpharAU2f)
);
// Update the fixedFluxPressure BCs to ensure flux consistency
@@ -238,12 +256,12 @@ while (pimple.correct())
{
surfaceScalarField phi1s
(
- phiHbyA1 + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1)
+ phiHbyA1 + alpharAU1f*mSfGradp
);
surfaceScalarField phi2s
(
- phiHbyA2 + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2)
+ phiHbyA2 + alpharAU2f*mSfGradp
);
surfaceScalarField phir
@@ -286,22 +304,12 @@ while (pimple.correct())
{
volVectorField U1s
(
- HbyA1
- + fvc::reconstruct
- (
- rAlphaAU1f*mSfGradp/fvc::interpolate(rho1)
- - phiF1
- )
+ HbyA1 + fvc::reconstruct(alpharAU1f*mSfGradp - phiF1)
);
volVectorField U2s
(
- HbyA2
- + fvc::reconstruct
- (
- rAlphaAU2f*mSfGradp/fvc::interpolate(rho2)
- - phiF2
- )
+ HbyA2 + fvc::reconstruct(alpharAU2f*mSfGradp - phiF2)
);
volScalarField D1(rAU1*dragCoeff);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
index f011791d259..c82e4db247f 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C
@@ -226,6 +226,61 @@ Foam::BlendedInterfacialModel<modelType>::F() const
}
+template<class modelType>
+Foam::tmp<Foam::volScalarField>
+Foam::BlendedInterfacialModel<modelType>::D() const
+{
+ tmp<volScalarField> f1, f2;
+
+ if (model_.valid() || model1In2_.valid())
+ {
+ f1 = blending_.f1(pair1In2_.dispersed(), pair2In1_.dispersed());
+ }
+
+ if (model_.valid() || model2In1_.valid())
+ {
+ f2 = blending_.f2(pair1In2_.dispersed(), pair2In1_.dispersed());
+ }
+
+ tmp<volScalarField> x
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ modelType::typeName + "Coeff",
+ pair_.phase1().mesh().time().timeName(),
+ pair_.phase1().mesh()
+ ),
+ pair_.phase1().mesh(),
+ dimensionedScalar("zero", modelType::dimD, 0)
+ )
+ );
+
+ if (model_.valid())
+ {
+ x() += model_->D()*(f1() - f2());
+ }
+
+ if (model1In2_.valid())
+ {
+ x() += model1In2_->D()*(1 - f1);
+ }
+
+ if (model2In1_.valid())
+ {
+ x() += model2In1_->D()*f2;
+ }
+
+ if (model_.valid() || model1In2_.valid() || model2In1_.valid())
+ {
+ correctFixedFluxBCs(x());
+ }
+
+ return x;
+}
+
+
template<class modelType>
const modelType& Foam::BlendedInterfacialModel<modelType>::phaseModel
(
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H
index de51082f865..aa3bb3df2e0 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H
@@ -113,13 +113,16 @@ public:
// Member Functions
- //- Return the implicit coefficient
+ //- Return the blended force coefficient
tmp<volScalarField> K() const;
- //- Return the explicit value
+ //- Return the blended force
template<class Type>
tmp<GeometricField<Type, fvPatchField, volMesh> > F() const;
+ //- Return the blended diffusivity
+ tmp<volScalarField> D() const;
+
//- Return the model for the supplied phase
const modelType& phaseModel(const phaseModel& phase) const;
};
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index f536e93a608..d7a58a38096 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -336,6 +336,13 @@ Foam::twoPhaseSystem::turbulentDispersionForce() const
}
+Foam::tmp<Foam::volScalarField>
+Foam::twoPhaseSystem::turbulentDiffusivity() const
+{
+ return turbulentDispersion_->D();
+}
+
+
void Foam::twoPhaseSystem::solve()
{
const Time& runTime = mesh_.time();
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
index 8b8ab19f41b..094670ca68f 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
@@ -155,7 +155,11 @@ public:
//- Return the wall lubrication force
tmp<volVectorField> wallLubricationForce() const;
- //- Return the wall lubrication force
+ //- Return the turbulent diffusivity
+ // Multiplies the phase-fraction gradient
+ tmp<volScalarField> turbulentDiffusivity() const;
+
+ //- Return the turbulent dispersion force
tmp<volVectorField> turbulentDispersionForce() const;
//- Solve for the two-phase-fractions
--
GitLab