diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H index 2a5f6c8bd0..9d4af292b2 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H @@ -5,7 +5,7 @@ volScalarField Cpv1("Cpv1", thermo1.Cpv()); volScalarField Cpv2("Cpv2", thermo2.Cpv()); - volScalarField heatTransferCoeff(fluid.heatTransferCoeff()); + volScalarField Kh(fluid.Kh()); fvScalarMatrix he1Eqn ( @@ -32,9 +32,9 @@ he1Eqn -= ( - heatTransferCoeff*(thermo2.T() - thermo1.T()) - + heatTransferCoeff*he1/Cpv1 - - fvm::Sp(heatTransferCoeff/Cpv1, he1) + Kh*(thermo2.T() - thermo1.T()) + + Kh*he1/Cpv1 + - fvm::Sp(Kh/Cpv1, he1) + fvOptions(alpha1, rho1, he1) ); @@ -63,9 +63,9 @@ he2Eqn -= ( - heatTransferCoeff*(thermo1.T() - thermo2.T()) - + heatTransferCoeff*he2/Cpv2 - - fvm::Sp(heatTransferCoeff/Cpv2, he2) + Kh*(thermo1.T() - thermo2.T()) + + Kh*he2/Cpv2 + - fvm::Sp(Kh/Cpv2, he2) + fvOptions(alpha2, rho2, he2) ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H index efb2310399..08282e78b3 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H @@ -70,52 +70,6 @@ fluid.U() ); - Info<< "Calculating field DDtU1 and DDtU2\n" << endl; - - volVectorField DDtU1 - ( - "DDtU1", - fvc::ddt(U1) - + fvc::div(phi1, U1) - - fvc::div(phi1)*U1 - ); - - volVectorField DDtU2 - ( - "DDtU2", - fvc::ddt(U2) - + fvc::div(phi2, U2) - - fvc::div(phi2)*U2 - ); - - volScalarField rAU1 - ( - IOobject - ( - IOobject::groupName("rAU", phase1.name()), - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0.0) - ); - - volScalarField rAU2 - ( - IOobject - ( - IOobject::groupName("rAU", phase2.name()), - runTime.timeName(), - mesh, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - mesh, - dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0.0) - ); - label pRefCell = 0; scalar pRefValue = 0.0; setRefCell @@ -144,3 +98,60 @@ Info<< "Creating field kinetic energy K\n" << endl; volScalarField K1(IOobject::groupName("K", phase1.name()), 0.5*magSqr(U1)); volScalarField K2(IOobject::groupName("K", phase2.name()), 0.5*magSqr(U2)); + + + volScalarField rAU1 + ( + IOobject + ( + IOobject::groupName("rAU", phase1.name()), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0) + ); + + volScalarField rAU2 + ( + IOobject + ( + IOobject::groupName("rAU", phase2.name()), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0) + ); + + surfaceScalarField rAU1f + ( + IOobject + ( + IOobject::groupName("rAUf", phase1.name()), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0) + ); + + surfaceScalarField rAU2f + ( + IOobject + ( + IOobject::groupName("rAUf", phase2.name()), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0) + ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C index e9cb5c3341..e4b9a8b276 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -26,6 +26,7 @@ License #include "dragModel.H" #include "phasePair.H" #include "swarmCorrection.H" +#include "surfaceInterpolate.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -102,12 +103,11 @@ Foam::dragModel::~dragModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::tmp Foam::dragModel::K() const +Foam::tmp Foam::dragModel::Ki() const { return 0.75 *CdRe() - *max(pair_.dispersed(), residualAlpha_) *swarmCorrection_->Cs() *pair_.continuous().rho() *pair_.continuous().nu() @@ -115,6 +115,20 @@ Foam::tmp Foam::dragModel::K() const } +Foam::tmp Foam::dragModel::K() const +{ + return max(pair_.dispersed(), residualAlpha_)*Ki(); +} + + +Foam::tmp Foam::dragModel::Kf() const +{ + return + max(fvc::interpolate(pair_.dispersed()), residualAlpha_) + *fvc::interpolate(Ki()); +} + + bool Foam::dragModel::writeData(Ostream& os) const { return os.good(); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H index 54c784412e..4d4c5693ba 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -130,14 +130,32 @@ public: // Member Functions + //- Return the residual phase-fraction + // used to stabilize the phase momentum as the phase-fraction -> 0 + const dimensionedScalar& residualAlpha() const + { + return residualAlpha_; + } + //- Drag coefficient virtual tmp CdRe() const = 0; - //- The drag function K used in the momentum equation + //- Return the phase-intensive drag coefficient Ki + // used in the momentum equations + // ddt(alpha1*rho1*U1) + ... = ... alphad*K*(U1-U2) + // ddt(alpha2*rho2*U2) + ... = ... alphad*K*(U2-U1) + virtual tmp Ki() const; + + //- Return the drag coefficient K + // used in the momentum equations // ddt(alpha1*rho1*U1) + ... = ... K*(U1-U2) // ddt(alpha2*rho2*U2) + ... = ... K*(U2-U1) virtual tmp K() const; + //- Return the drag coefficient Kf + // used in the face-momentum equations + virtual tmp Kf() const; + //- Dummy write for regIOobject bool writeData(Ostream& os) const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.C index 9f3df40829..cb6a9fc34f 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.C @@ -26,6 +26,7 @@ License #include "segregated.H" #include "phasePair.H" #include "fvcGrad.H" +#include "surfaceInterpolate.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -147,4 +148,10 @@ Foam::tmp Foam::dragModels::segregated::K() const } +Foam::tmp Foam::dragModels::segregated::Kf() const +{ + return fvc::interpolate(K()); +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.H index dac7c04d20..3379c02103 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.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 @@ -99,6 +99,9 @@ public: //- The drag function used in the momentum equation virtual tmp K() const; + + //- The drag function Kf used in the face-momentum equations + virtual tmp Kf() const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C index c48a05f646..79e52a2856 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C @@ -26,7 +26,7 @@ License #include "liftModel.H" #include "phasePair.H" #include "fvcCurl.H" -#include "addToRunTimeSelectionTable.H" +#include "surfaceInterpolate.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -59,11 +59,10 @@ Foam::liftModel::~liftModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::tmp Foam::liftModel::F() const +Foam::tmp Foam::liftModel::Fi() const { return Cl() - *pair_.dispersed() *pair_.continuous().rho() *( pair_.Ur() ^ fvc::curl(pair_.continuous().U()) @@ -71,4 +70,20 @@ Foam::tmp Foam::liftModel::F() const } +Foam::tmp Foam::liftModel::F() const +{ + return pair_.dispersed()*Fi(); +} + + +Foam::tmp Foam::liftModel::Ff() const +{ + const fvMesh& mesh(this->pair_.phase1().mesh()); + + return + fvc::interpolate(pair_.dispersed()) + *(fvc::interpolate(Fi()) & mesh.Sf()); +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.H index d20d43216b..665c95dc38 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.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 @@ -112,11 +112,17 @@ public: // Member Functions - //- Lift coefficient + //- Return lift coefficient virtual tmp Cl() const = 0; - //- Lift force + //- Return phase-intensive lift force + virtual tmp Fi() const; + + //- Return lift force virtual tmp F() const; + + //- Return face lift force + virtual tmp Ff() const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/noLift/noLift.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/noLift/noLift.C index d6ad1c89f8..d4e13d5c7f 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/noLift/noLift.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/noLift/noLift.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 @@ -63,34 +63,71 @@ Foam::tmp Foam::liftModels::noLift::Cl() const { const fvMesh& mesh(this->pair_.phase1().mesh()); - return - tmp + return tmp + ( + new volScalarField ( - new volScalarField + IOobject ( - IOobject - ( - "Cl", - mesh.time().timeName(), - mesh - ), + "Cl", + mesh.time().timeName(), mesh, - dimensionedScalar("Cl", dimless, 0) - ) - ); + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimensionedScalar("Cl", dimless, 0) + ) + ); } Foam::tmp Foam::liftModels::noLift::F() const { - return - Cl() - *dimensionedVector + const fvMesh& mesh(this->pair_.phase1().mesh()); + + return tmp + ( + new volVectorField ( - "zero", - dimensionSet(1, -2, -2, 0, 0), - vector::zero - ); + IOobject + ( + "noLift:F", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimensionedVector("zero", dimF, vector::zero) + ) + ); +} + + +Foam::tmp Foam::liftModels::noLift::Ff() const +{ + const fvMesh& mesh(this->pair_.phase1().mesh()); + + return tmp + ( + new surfaceScalarField + ( + IOobject + ( + "noLift:Ff", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimensionedScalar("zero", dimF*dimArea, 0) + ) + ); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/noLift/noLift.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/noLift/noLift.H index 0b0760a67d..db3a90cb13 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/noLift/noLift.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/noLift/noLift.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 @@ -81,6 +81,9 @@ public: //- Lift force virtual tmp F() const; + + //- Lift force on faces + virtual tmp Ff() const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/constantVirtualMassCoefficient/constantVirtualMassCoefficient.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/constantVirtualMassCoefficient/constantVirtualMassCoefficient.C index bdf1ab08b4..2fa5f49305 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/constantVirtualMassCoefficient/constantVirtualMassCoefficient.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/constantVirtualMassCoefficient/constantVirtualMassCoefficient.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 @@ -73,21 +73,23 @@ Foam::virtualMassModels::constantVirtualMassCoefficient::Cvm() const { const fvMesh& mesh(this->pair_.phase1().mesh()); - return - tmp + return tmp + ( + new volScalarField ( - new volScalarField + IOobject ( - IOobject - ( - "zero", - mesh.time().timeName(), - mesh - ), + "Cvm", + mesh.time().timeName(), mesh, - Cvm_ - ) - ); + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + Cvm_ + ) + ); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.C index 34473a6c31..a2ece4192d 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.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 @@ -25,6 +25,7 @@ License #include "virtualMassModel.H" #include "phasePair.H" +#include "surfaceInterpolate.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -70,9 +71,22 @@ Foam::virtualMassModel::~virtualMassModel() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +Foam::tmp Foam::virtualMassModel::Ki() const +{ + return Cvm()*pair_.continuous().rho(); +} + + Foam::tmp Foam::virtualMassModel::K() const { - return Cvm()*pair_.dispersed()*pair_.continuous().rho(); + return pair_.dispersed()*Ki(); +} + + +Foam::tmp Foam::virtualMassModel::Kf() const +{ + return + fvc::interpolate(pair_.dispersed())*fvc::interpolate(Ki()); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.H index 65ea9fd547..0af7097cf9 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.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 @@ -116,14 +116,25 @@ public: // Member Functions - //- Virtual mass coefficient + //- Return the virtual mass coefficient virtual tmp Cvm() const = 0; - //- The virtual mass function K used in the momentum equation + //- Return the phase-intensive virtual mass coefficient Ki + // used in the momentum equation + // ddt(alpha1*rho1*U1) + ... = ... alphad*K*(DU1_Dt - DU2_Dt) + // ddt(alpha2*rho2*U2) + ... = ... alphad*K*(DU1_Dt - DU2_Dt) + virtual tmp Ki() const; + + //- Return the virtual mass coefficient K + // used in the momentum equation // ddt(alpha1*rho1*U1) + ... = ... K*(DU1_Dt - DU2_Dt) // ddt(alpha2*rho2*U2) + ... = ... K*(DU1_Dt - DU2_Dt) virtual tmp K() const; + //- Return the virtual mass coefficient Kf + // used in the face-momentum equations + virtual tmp Kf() const; + // Dummy write for regIOobject bool writeData(Ostream& os) const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.C index 208d6b2d10..fa9b2749f0 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.C @@ -66,7 +66,7 @@ Foam::wallLubricationModels::Antal::~Antal() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::tmp Foam::wallLubricationModels::Antal::F() const +Foam::tmp Foam::wallLubricationModels::Antal::Fi() const { volVectorField Ur(pair_.Ur()); @@ -78,7 +78,6 @@ Foam::tmp Foam::wallLubricationModels::Antal::F() const dimensionedScalar("zero", dimless/dimLength, 0), Cw1_/pair_.dispersed().d() + Cw2_/yWall() ) - *pair_.dispersed() *pair_.continuous().rho() *magSqr(Ur - (Ur & n)*n) *n; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.H index 3be657da18..7639712417 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.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 @@ -95,8 +95,8 @@ public: // Member Functions - //- Wall lubrication force - tmp F() const; + //- Return phase-intensive wall lubrication force + tmp Fi() const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C index 53ab61493c..f1b11a32cc 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C @@ -67,7 +67,7 @@ Foam::wallLubricationModels::Frank::~Frank() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::tmp Foam::wallLubricationModels::Frank::F() const +Foam::tmp Foam::wallLubricationModels::Frank::Fi() const { volVectorField Ur(pair_.Ur()); @@ -88,7 +88,6 @@ Foam::tmp Foam::wallLubricationModels::Frank::F() const dimensionedScalar("zero", dimless/dimLength, 0.0), (1.0 - yTilde)/(Cwd_*y*pow(yTilde, p_ - 1.0)) ) - *pair_.dispersed() *pair_.continuous().rho() *magSqr(Ur - (Ur & n)*n) *n; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.H index 1b67b16d33..51527fd80b 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.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 @@ -105,8 +105,8 @@ public: // Member Functions - //- Wall lubrication force - tmp F() const; + //- Return phase-intensive wall lubrication force + tmp Fi() const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C index e28cb9b29b..c3db3597af 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C @@ -66,7 +66,7 @@ Foam::wallLubricationModels::TomiyamaWallLubrication::~TomiyamaWallLubrication() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // Foam::tmp -Foam::wallLubricationModels::TomiyamaWallLubrication::F() const +Foam::wallLubricationModels::TomiyamaWallLubrication::Fi() const { volVectorField Ur(pair_.Ur()); @@ -87,7 +87,6 @@ Foam::wallLubricationModels::TomiyamaWallLubrication::F() const 1/sqr(y) - 1/sqr(D_ - y) ) - *pair_.dispersed() *pair_.continuous().rho() *magSqr(Ur - (Ur & n)*n) *n; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.H index 426b0d3ec4..9d37bbbd17 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.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,8 +98,8 @@ public: // Member Functions - //- Wall lubrication force - tmp F() const; + //- Return phase-intensive wall lubrication force + tmp Fi() const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/noWallLubrication/noWallLubrication.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/noWallLubrication/noWallLubrication.C index 4fbdffc085..8de7e30260 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/noWallLubrication/noWallLubrication.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/noWallLubrication/noWallLubrication.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 @@ -64,31 +64,53 @@ Foam::wallLubricationModels::noWallLubrication::~noWallLubrication() // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +Foam::tmp +Foam::wallLubricationModels::noWallLubrication::Fi() const +{ + const fvMesh& mesh(this->pair_.phase1().mesh()); + + return tmp + ( + new volVectorField + ( + IOobject + ( + "noWallLubrication:Fi", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimensionedVector("zero", dimF, vector::zero) + ) + ); +} + + Foam::tmp Foam::wallLubricationModels::noWallLubrication::F() const { const fvMesh& mesh(this->pair_.phase1().mesh()); - return - tmp + return tmp + ( + new volVectorField ( - new volVectorField + IOobject ( - IOobject - ( - "zero", - mesh.time().timeName(), - mesh - ), + "noWallLubrication:F", + mesh.time().timeName(), mesh, - dimensionedVector - ( - "zero", - dimensionSet(1, -2, -2, 0, 0), - vector::zero - ) - ) - ); + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimensionedVector("zero", dimF, vector::zero) + ) + ); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/noWallLubrication/noWallLubrication.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/noWallLubrication/noWallLubrication.H index 4861c25675..81a2c0327e 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/noWallLubrication/noWallLubrication.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/noWallLubrication/noWallLubrication.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,6 +76,9 @@ public: // Member Functions + //- Return phase-intensive wall lubrication force + tmp Fi() const; + //- Wall lubrication force tmp F() const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C index b7a7edc385..2f3cb1aaa0 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.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 @@ -25,8 +25,7 @@ License #include "wallLubricationModel.H" #include "phasePair.H" - -const Foam::dimensionSet Foam::wallLubricationModel::dimF(1, -2, -2, 0, 0); +#include "surfaceInterpolate.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -36,6 +35,8 @@ namespace Foam defineRunTimeSelectionTable(wallLubricationModel, dictionary); } +const Foam::dimensionSet Foam::wallLubricationModel::dimF(1, -2, -2, 0, 0); + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -56,4 +57,22 @@ Foam::wallLubricationModel::~wallLubricationModel() {} +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::tmp Foam::wallLubricationModel::F() const +{ + return pair_.dispersed()*Fi(); +} + + +Foam::tmp Foam::wallLubricationModel::Ff() const +{ + const fvMesh& mesh(this->pair_.phase1().mesh()); + + return + fvc::interpolate(pair_.dispersed()) + *(fvc::interpolate(Fi()) & mesh.Sf()); +} + + // ************************************************************************* // diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.H index 8f77e21352..b55415aeb6 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.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 @@ -115,8 +115,14 @@ public: // Member Functions - //- Wall lubrication force - virtual tmp F() const = 0; + //- Return phase-intensive wall lubrication force + virtual tmp Fi() const = 0; + + //- Return wall lubrication force + virtual tmp F() const; + + //- Return face wall lubrication force + virtual tmp Ff() const; }; diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/DDtU.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/DDtU.H similarity index 100% rename from applications/solvers/multiphase/twoPhaseEulerFoam/DDtU.H rename to applications/solvers/multiphase/twoPhaseEulerFoam/pU/DDtU.H diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/UEqns.H similarity index 78% rename from applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H rename to applications/solvers/multiphase/twoPhaseEulerFoam/pU/UEqns.H index abd3c8b646..0814a162db 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/UEqns.H @@ -5,20 +5,20 @@ mrfZones.correctBoundaryVelocity(U); fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime); fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime); -volScalarField dragCoeff(fluid.dragCoeff()); +volScalarField Kd(fluid.Kd()); { - volScalarField virtualMassCoeff(fluid.virtualMassCoeff()); + volScalarField Vm(fluid.Vm()); { U1Eqn = ( fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1) - fvm::Sp(contErr1, U1) - + mrfZones(alpha1*rho1 + virtualMassCoeff, U1) + + mrfZones(alpha1*rho1 + Vm, U1) + phase1.turbulence().divDevRhoReff(U1) == - - virtualMassCoeff + - Vm *( fvm::ddt(U1) + fvm::div(phi1, U1) @@ -28,7 +28,7 @@ volScalarField dragCoeff(fluid.dragCoeff()); + fvOptions(alpha1, rho1, U1) ); U1Eqn.relax(); - U1Eqn += fvm::Sp(dragCoeff, U1); + U1Eqn += fvm::Sp(Kd, U1); fvOptions.constrain(U1Eqn); U1.correctBoundaryConditions(); fvOptions.correct(U1); @@ -39,10 +39,10 @@ volScalarField dragCoeff(fluid.dragCoeff()); ( fvm::ddt(alpha2, rho2, U2) + fvm::div(alphaRhoPhi2, U2) - fvm::Sp(contErr2, U2) - + mrfZones(alpha2*rho2 + virtualMassCoeff, U2) + + mrfZones(alpha2*rho2 + Vm, U2) + phase2.turbulence().divDevRhoReff(U2) == - - virtualMassCoeff + - Vm *( fvm::ddt(U2) + fvm::div(phi2, U2) @@ -52,7 +52,7 @@ volScalarField dragCoeff(fluid.dragCoeff()); + fvOptions(alpha2, rho2, U2) ); U2Eqn.relax(); - U2Eqn += fvm::Sp(dragCoeff, U2); + U2Eqn += fvm::Sp(Kd, U2); fvOptions.constrain(U2Eqn); U2.correctBoundaryConditions(); fvOptions.correct(U2); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pU/createDDtU.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/createDDtU.H new file mode 100644 index 0000000000..f71a265817 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/createDDtU.H @@ -0,0 +1,17 @@ + Info<< "Calculating field DDtU1 and DDtU2\n" << endl; + + volVectorField DDtU1 + ( + "DDtU1", + fvc::ddt(U1) + + fvc::div(phi1, U1) + - fvc::div(phi1)*U1 + ); + + volVectorField DDtU2 + ( + "DDtU2", + fvc::ddt(U2) + + fvc::div(phi2, U2) + - fvc::div(phi2)*U2 + ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H similarity index 91% rename from applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H rename to applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H index dba4ad8a36..fe3ba4350a 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H @@ -39,35 +39,26 @@ while (pimple.correct()) // Turbulent diffusion, particle-pressure lift and wall-lubrication fluxes { - volScalarField turbulentDiffusivity(fluid.turbulentDiffusivity()); - volVectorField liftForce(fluid.liftForce()); - volVectorField wallLubricationForce(fluid.wallLubricationForce()); + volScalarField D(fluid.D()); + volVectorField F(fluid.F()); surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf()); phiF1 = ( fvc::interpolate ( - rAU1*(turbulentDiffusivity + phase1.turbulence().pPrime()) + rAU1*(D + phase1.turbulence().pPrime()) )*snGradAlpha1 - - + ( - fvc::interpolate(rAU1*(wallLubricationForce + liftForce)) - & mesh.Sf() - ) + + (fvc::interpolate(rAU1*F) & mesh.Sf()) ); phiF2 = ( - fvc::interpolate ( - rAU2*(turbulentDiffusivity + phase2.turbulence().pPrime()) + rAU2*(D + phase2.turbulence().pPrime()) )*snGradAlpha1 - - - ( - fvc::interpolate(rAU2*(wallLubricationForce + liftForce)) - & mesh.Sf() - ) + - (fvc::interpolate(rAU2*F) & mesh.Sf()) ); } @@ -139,8 +130,8 @@ while (pimple.correct()) ); // Face-drag coefficients - surfaceScalarField D1f(fvc::interpolate(rAU1*dragCoeff)); - surfaceScalarField D2f(fvc::interpolate(rAU2*dragCoeff)); + surfaceScalarField D1f(fvc::interpolate(rAU1*Kd)); + surfaceScalarField D2f(fvc::interpolate(rAU2*Kd)); // Construct the mean predicted flux // including explicit drag contributions based on absolute fluxes @@ -306,8 +297,8 @@ while (pimple.correct()) HbyA2 + fvc::reconstruct(alpharAU2f*mSfGradp - phiF2) ); - volScalarField D1(rAU1*dragCoeff); - volScalarField D2(rAU2*dragCoeff); + volScalarField D1(rAU1*Kd); + volScalarField D2(rAU2*Kd); U = alpha1*(U1s + D1*U2) + alpha2*(U2s + D2*U1); volVectorField Ur(((1 - D2)*U1s - (1 - D1)*U2s)/(1 - D1*D2)); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/DDtU.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/DDtU.H new file mode 100644 index 0000000000..63b61cd61a --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/DDtU.H @@ -0,0 +1,9 @@ + ddtPhi1 = + ( + (phi1 - phi1.oldTime())/runTime.deltaT() + ); + + ddtPhi2 = + ( + (phi2 - phi2.oldTime())/runTime.deltaT() + ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/UEqns.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/UEqns.H new file mode 100644 index 0000000000..36559ade99 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/UEqns.H @@ -0,0 +1,50 @@ +mrfZones.correctBoundaryVelocity(U1); +mrfZones.correctBoundaryVelocity(U2); +mrfZones.correctBoundaryVelocity(U); + +fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime); +fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime); + +{ + volScalarField Vm(fluid.Vm()); + + fvVectorMatrix UgradU1 + ( + fvm::div(phi1, U1) - fvm::Sp(fvc::div(phi1), U1) + + mrfZones(U1) + ); + + fvVectorMatrix UgradU2 + ( + fvm::div(phi2, U2) - fvm::Sp(fvc::div(phi2), U2) + + mrfZones(U2) + ); + + { + U1Eqn = + ( + fvm::div(alphaRhoPhi1, U1) - fvm::Sp(fvc::div(alphaRhoPhi1), U1) + + mrfZones(alpha1*rho1, U1) + + phase1.turbulence().divDevRhoReff(U1) + + Vm*(UgradU1 - (UgradU2 & U2)) + ); + U1Eqn.relax(); + fvOptions.constrain(U1Eqn); + U1.correctBoundaryConditions(); + fvOptions.correct(U1); + } + + { + U2Eqn = + ( + fvm::div(alphaRhoPhi2, U2) - fvm::Sp(fvc::div(alphaRhoPhi2), U2) + + mrfZones(alpha2*rho2, U2) + + phase2.turbulence().divDevRhoReff(U2) + + Vm*(UgradU2 - (UgradU1 & U1)) + ); + U2Eqn.relax(); + fvOptions.constrain(U2Eqn); + U2.correctBoundaryConditions(); + fvOptions.correct(U2); + } +} diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/createDDtU.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/createDDtU.H new file mode 100644 index 0000000000..7445d4af23 --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/createDDtU.H @@ -0,0 +1,9 @@ + surfaceScalarField ddtPhi1 + ( + (phi1 - phi1.oldTime())/runTime.deltaT() + ); + + surfaceScalarField ddtPhi2 + ( + (phi2 - phi2.oldTime())/runTime.deltaT() + ); diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H new file mode 100644 index 0000000000..0e5a47a8ff --- /dev/null +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H @@ -0,0 +1,347 @@ +surfaceScalarField alpha1f("alpha1f", fvc::interpolate(alpha1)); +surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f); + +surfaceScalarField alphaRho1f0 +( + "alphaRho1f0", + fvc::interpolate + ( + max(alpha1.oldTime(), fluid.residualAlpha(phase1)) + *rho1.oldTime() + ) +); + +surfaceScalarField alphaRho2f0 +( + "alphaRho2f0", + fvc::interpolate + ( + max(alpha2.oldTime(), fluid.residualAlpha(phase2)) + *rho2.oldTime() + ) +); + +// Drag coefficient +surfaceScalarField Kdf("Kdf", fluid.Kdf()); + +// Virtual-mass coefficient +surfaceScalarField Vmf("Vmf", fluid.Vmf()); + +// Lift and wall-lubrication forces +surfaceScalarField Ff("Ff", fluid.Ff()); + +tmp snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf()); +tmp D(fluid.D()); + +// Turbulent-dispersion force for phase 1 +surfaceScalarField Ftdf1 +( + IOobject::groupName("Ftdf", phase1.name()), + fvc::interpolate(D() + phase1.turbulence().pPrime())*snGradAlpha1() +); + +// Turbulent-dispersion force for phase 2 +surfaceScalarField Ftdf2 +( + IOobject::groupName("Ftdf", phase2.name()), + fvc::interpolate(D + phase2.turbulence().pPrime())*snGradAlpha1 +); + + +while (pimple.correct()) +{ + // Update continuity errors due to temperature changes + #include "correctContErrs.H" + + surfaceScalarField rho1f(fvc::interpolate(rho1)); + surfaceScalarField rho2f(fvc::interpolate(rho2)); + + // Correct flux BCs to be consistent with the velocity BCs + phi1.boundaryField() == + mrfZones.relative(mesh.Sf().boundaryField() & U1.boundaryField()); + phi2.boundaryField() == + mrfZones.relative(mesh.Sf().boundaryField() & U2.boundaryField()); + + rAU1f = + ( + IOobject::groupName("rAUf", phase1.name()), + 1.0 + /( + (alphaRho1f0 + Vmf)/runTime.deltaT() + + fvc::interpolate(U1Eqn.A()) + + Kdf + ) + ); + + rAU2f = + ( + IOobject::groupName("rAUf", phase2.name()), + 1.0 + /( + (alphaRho2f0 + Vmf)/runTime.deltaT() + + fvc::interpolate(U2Eqn.A()) + + Kdf + ) + ); + + surfaceScalarField rAlphaAU1f + ( + IOobject::groupName("rAlphaAUf", phase1.name()), + max(alpha1f, fluid.residualAlpha(phase1))*rAU1f + ); + + surfaceScalarField rAlphaAU2f + ( + IOobject::groupName("rAlphaAUf", phase2.name()), + max(alpha2f, fluid.residualAlpha(phase2))*rAU2f + ); + + + volScalarField rho("rho", fluid.rho()); + surfaceScalarField ghSnGradRho + ( + "ghSnGradRho", + ghf*fvc::snGrad(rho)*mesh.magSf() + ); + + // Add the phase-1 buoyancy force + surfaceScalarField phiF1 + ( + IOobject::groupName("phiF", phase1.name()), + rAlphaAU1f + *( + ghSnGradRho + - alpha2f*(rho1f - rho2f)*(g & mesh.Sf()) + ) + ); + + // Add the phase-2 buoyancy force + surfaceScalarField phiF2 + ( + IOobject::groupName("phiF", phase2.name()), + rAlphaAU2f + *( + ghSnGradRho + - alpha1f*(rho2f - rho1f)*(g & mesh.Sf()) + ) + ); + + + // Phase-1 predicted flux + surfaceScalarField phiHbyA1 + ( + IOobject::groupName("phiHbyA", phase1.name()), + phi1 + ); + + phiHbyA1 = + rAU1f + *( + (alphaRho1f0 + Vmf) + *mrfZones.absolute(phi1.oldTime())/runTime.deltaT() + + (fvc::interpolate(U1Eqn.H()) & mesh.Sf()) + + Vmf*ddtPhi2 + + Kdf*mrfZones.absolute(phi2) + - Ff + - Ftdf1 + ); + + // Phase-2 predicted flux + surfaceScalarField phiHbyA2 + ( + IOobject::groupName("phiHbyA", phase2.name()), + phi2 + ); + + phiHbyA2 = + rAU2f + *( + (alphaRho2f0 + Vmf) + *mrfZones.absolute(phi2.oldTime())/runTime.deltaT() + + (fvc::interpolate(U2Eqn.H()) & mesh.Sf()) + + Vmf*ddtPhi1 + + Kdf*mrfZones.absolute(phi1) + + Ff + + Ftdf2 + ); + + + surfaceScalarField phiHbyA + ( + "phiHbyA", + alpha1f*(phiHbyA1 - phiF1) + alpha2f*(phiHbyA2 - phiF2) + ); + mrfZones.makeRelative(phiHbyA); + + phiHbyA1 -= phiF1; + phiHbyA2 -= phiF2; + + surfaceScalarField rAUf + ( + "rAUf", + mag(alpha1f*rAlphaAU1f + alpha2f*rAlphaAU2f) + ); + + // Update the fixedFluxPressure BCs to ensure flux consistency + setSnGrad + ( + p_rgh.boundaryField(), + ( + phiHbyA.boundaryField() + - ( + alpha1f.boundaryField()*phi1.boundaryField() + + alpha2f.boundaryField()*phi2.boundaryField() + ) + )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) + ); + + tmp pEqnComp1; + tmp 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_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); + + phi = phiHbyA + pEqnIncomp.flux(); + + surfaceScalarField phi1s + ( + phiHbyA1 + + rAlphaAU1f*mSfGradp + - rAU1f*Kdf*mrfZones.absolute(phi2) + ); + + surfaceScalarField phi2s + ( + phiHbyA2 + + rAlphaAU2f*mSfGradp + - rAU2f*Kdf*mrfZones.absolute(phi1) + ); + + surfaceScalarField phir + ( + ((phi2s + rAU2f*Kdf*phi1s) - (phi1s + rAU1f*Kdf*phi2s)) + /(1.0 - rAU1f*rAU2f*sqr(Kdf)) + ); + + phi1 = phi - alpha2f*phir; + phi2 = phi + alpha1f*phir; + + U1 = fvc::reconstruct(mrfZones.absolute(phi1)); + U1.correctBoundaryConditions(); + fvOptions.correct(U1); + + U2 = fvc::reconstruct(mrfZones.absolute(phi2)); + U2.correctBoundaryConditions(); + fvOptions.correct(U2); + + U = fluid.U(); + + fluid.dgdt() = + ( + alpha1*(pEqnComp2 & p_rgh) + - alpha2*(pEqnComp1 & p_rgh) + ); + } + } + + // Update and limit the static pressure + p = max(p_rgh + rho*gh, pMin); + + // Limit p_rgh + p_rgh = p - rho*gh; + + // Update densities from change in p_rgh + rho1 += psi1*(p_rgh - p_rgh_0); + rho2 += psi2*(p_rgh - p_rgh_0); + + // Correct p_rgh for consistency with p and the updated densities + rho = fluid.rho(); + p_rgh = p - rho*gh; + p_rgh.correctBoundaryConditions(); +} + +// Update the phase kinetic energies +K1 = 0.5*magSqr(U1); +K2 = 0.5*magSqr(U2); + +// Update the pressure time-derivative if required +if (thermo1.dpdt() || thermo2.dpdt()) +{ + dpdt = fvc::ddt(p); +} + +#include "pUf/DDtU.H" diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C index f518de6177..c9679c3439 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C @@ -50,6 +50,11 @@ int main(int argc, char *argv[]) pimpleControl pimple(mesh); + Switch faceMomentum + ( + pimple.dict().lookupOrDefault("faceMomentum", false) + ); + #include "createFields.H" #include "createMRFZones.H" #include "createFvOptions.H" @@ -58,6 +63,9 @@ int main(int argc, char *argv[]) #include "CourantNos.H" #include "setInitialDeltaT.H" + #include "pUf/createDDtU.H" + #include "pU/createDDtU.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -78,10 +86,21 @@ int main(int argc, char *argv[]) fluid.correct(); #include "contErrs.H" - #include "UEqns.H" #include "EEqns.H" - #include "pEqn.H" - #include "DDtU.H" + + if (faceMomentum) + { + Info<< "Constructing face momentum equations" << endl; + #include "pUf/UEqns.H" + #include "pUf/pEqn.H" + } + else + { + Info<< "Constructing momentum equations" << endl; + #include "pU/UEqns.H" + #include "pU/pEqn.H" + #include "pU/DDtU.H" + } if (pimple.turbCorr()) { diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C index 93d9bb4c78..090ffd1bcc 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C @@ -25,14 +25,15 @@ License #include "BlendedInterfacialModel.H" #include "fixedValueFvsPatchFields.H" +#include "surfaceInterpolate.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // template -template +template void Foam::BlendedInterfacialModel::correctFixedFluxBCs ( - GeometricField& field + GeometricField& field ) const { forAll(pair_.phase1().phi().boundaryField(), patchI) @@ -45,7 +46,8 @@ void Foam::BlendedInterfacialModel::correctFixedFluxBCs ) ) { - field.boundaryField()[patchI] = pTraits::zero; + field.boundaryField()[patchI] + = pTraits::zero; } } } @@ -139,9 +141,12 @@ Foam::BlendedInterfacialModel::K() const ( IOobject ( - modelType::typeName + "Coeff", + modelType::typeName + ":K", pair_.phase1().mesh().time().timeName(), - pair_.phase1().mesh() + pair_.phase1().mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false ), pair_.phase1().mesh(), dimensionedScalar("zero", modelType::dimK, 0) @@ -176,6 +181,74 @@ Foam::BlendedInterfacialModel::K() const } +template +Foam::tmp +Foam::BlendedInterfacialModel::Kf() const +{ + tmp f1, f2; + + if (model_.valid() || model1In2_.valid()) + { + f1 = fvc::interpolate + ( + blending_.f1(pair1In2_.dispersed(), pair2In1_.dispersed()) + ); + } + + if (model_.valid() || model2In1_.valid()) + { + f2 = fvc::interpolate + ( + blending_.f2(pair1In2_.dispersed(), pair2In1_.dispersed()) + ); + } + + tmp x + ( + new surfaceScalarField + ( + IOobject + ( + modelType::typeName + ":Kf", + pair_.phase1().mesh().time().timeName(), + pair_.phase1().mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + pair_.phase1().mesh(), + dimensionedScalar("zero", modelType::dimK, 0) + ) + ); + + if (model_.valid()) + { + x() += model_->Kf()*(f1() - f2()); + } + + if (model1In2_.valid()) + { + x() += model1In2_->Kf()*(1 - f1); + } + + if (model2In1_.valid()) + { + x() += model2In1_->Kf()*f2; + } + + if + ( + correctFixedFluxBCs_ + && (model_.valid() || model1In2_.valid() || model2In1_.valid()) + ) + { + correctFixedFluxBCs(x()); + } + + return x; +} + + template template Foam::tmp > @@ -199,9 +272,12 @@ Foam::BlendedInterfacialModel::F() const ( IOobject ( - modelType::typeName + "Coeff", + modelType::typeName + ":F", pair_.phase1().mesh().time().timeName(), - pair_.phase1().mesh() + pair_.phase1().mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false ), pair_.phase1().mesh(), dimensioned("zero", modelType::dimF, pTraits::zero) @@ -236,6 +312,74 @@ Foam::BlendedInterfacialModel::F() const } +template +Foam::tmp +Foam::BlendedInterfacialModel::Ff() const +{ + tmp f1, f2; + + if (model_.valid() || model1In2_.valid()) + { + f1 = fvc::interpolate + ( + blending_.f1(pair1In2_.dispersed(), pair2In1_.dispersed()) + ); + } + + if (model_.valid() || model2In1_.valid()) + { + f2 = fvc::interpolate + ( + blending_.f2(pair1In2_.dispersed(), pair2In1_.dispersed()) + ); + } + + tmp x + ( + new surfaceScalarField + ( + IOobject + ( + modelType::typeName + ":Ff", + pair_.phase1().mesh().time().timeName(), + pair_.phase1().mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + pair_.phase1().mesh(), + dimensionedScalar("zero", modelType::dimF*dimArea, 0) + ) + ); + + if (model_.valid()) + { + x() += model_->Ff()*(f1() - f2()); + } + + if (model1In2_.valid()) + { + x() += model1In2_->Ff()*(1 - f1); + } + + if (model2In1_.valid()) + { + x() -= model2In1_->Ff()*f2; // note : subtraction + } + + if + ( + correctFixedFluxBCs_ + && (model_.valid() || model1In2_.valid() || model2In1_.valid()) + ) + { + correctFixedFluxBCs(x()); + } + + return x; +} + + template Foam::tmp Foam::BlendedInterfacialModel::D() const @@ -258,9 +402,12 @@ Foam::BlendedInterfacialModel::D() const ( IOobject ( - modelType::typeName + "Coeff", + modelType::typeName + ":D", pair_.phase1().mesh().time().timeName(), - pair_.phase1().mesh() + pair_.phase1().mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false ), pair_.phase1().mesh(), dimensionedScalar("zero", modelType::dimD, 0) @@ -295,6 +442,19 @@ Foam::BlendedInterfacialModel::D() const } +template +bool Foam::BlendedInterfacialModel::hasModel +( + const class phaseModel& phase +) const +{ + return + &phase == &(pair_.phase1()) + ? model1In2_.valid() + : model2In1_.valid(); +} + + template const modelType& Foam::BlendedInterfacialModel::phaseModel ( diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H index f56ab025b6..e6493fa1cf 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H @@ -88,11 +88,8 @@ class BlendedInterfacialModel void operator=(const BlendedInterfacialModel&); //- Correct coeff/value on fixed flux boundary conditions - template - void correctFixedFluxBCs - ( - GeometricField& field - ) const; + template + void correctFixedFluxBCs(GeometricField& field) const; public: @@ -117,18 +114,27 @@ public: // Member Functions + //- Return true if a model is specified for the supplied phase + bool hasModel(const phaseModel& phase) const; + + //- Return the model for the supplied phase + const modelType& phaseModel(const phaseModel& phase) const; + //- Return the blended force coefficient tmp K() const; + //- Return the face blended force coefficient + tmp Kf() const; + //- Return the blended force template tmp > F() const; + //- Return the face blended force + tmp Ff() const; + //- Return the blended diffusivity tmp 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 32ca0ace85..391b38a76a 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C @@ -26,7 +26,6 @@ License #include "twoPhaseSystem.H" #include "PhaseCompressibleTurbulenceModel.H" #include "BlendedInterfacialModel.H" -#include "dragModel.H" #include "virtualMassModel.H" #include "heatTransferModel.H" #include "liftModel.H" @@ -48,6 +47,15 @@ License #include "blendingMethod.H" #include "HashPtrTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +Foam::dimensionedScalar Foam::twoPhaseSystem::zeroResidualAlpha_ +( + "zeroResidualAlpha", dimless, 0 +); + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::twoPhaseSystem::twoPhaseSystem @@ -299,46 +307,49 @@ Foam::tmp Foam::twoPhaseSystem::calcPhi() const } -Foam::tmp Foam::twoPhaseSystem::dragCoeff() const +Foam::tmp Foam::twoPhaseSystem::Kd() const { return drag_->K(); } -Foam::tmp Foam::twoPhaseSystem::virtualMassCoeff() const +Foam::tmp Foam::twoPhaseSystem::Kdf() const +{ + return drag_->Kf(); +} + + +Foam::tmp Foam::twoPhaseSystem::Vm() const { return virtualMass_->K(); } -Foam::tmp Foam::twoPhaseSystem::heatTransferCoeff() const +Foam::tmp Foam::twoPhaseSystem::Vmf() const +{ + return virtualMass_->Kf(); +} + + +Foam::tmp Foam::twoPhaseSystem::Kh() const { return heatTransfer_->K(); } -Foam::tmp Foam::twoPhaseSystem::liftForce() const +Foam::tmp Foam::twoPhaseSystem::F() const { - return lift_->F(); + return lift_->F() + wallLubrication_->F(); } -Foam::tmp -Foam::twoPhaseSystem::wallLubricationForce() const +Foam::tmp Foam::twoPhaseSystem::Ff() const { - return wallLubrication_->F(); + return lift_->Ff() + wallLubrication_->Ff(); } -Foam::tmp -Foam::twoPhaseSystem::turbulentDispersionForce() const -{ - return turbulentDispersion_->F(); -} - - -Foam::tmp -Foam::twoPhaseSystem::turbulentDiffusivity() const +Foam::tmp Foam::twoPhaseSystem::D() const { return turbulentDispersion_->D(); } @@ -354,6 +365,13 @@ void Foam::twoPhaseSystem::solve() const surfaceScalarField& phi1 = phase1_.phi(); const surfaceScalarField& phi2 = phase2_.phi(); + Switch faceMomentum + ( + //pimple.dict().lookupOrDefault("faceMomentum", false) + mesh_.solutionDict().subDict("PIMPLE") + .lookupOrDefault("faceMomentum", false) + ); + const dictionary& alphaControls = mesh_.solverDict ( alpha1.name() @@ -381,18 +399,40 @@ void Foam::twoPhaseSystem::solve() if (implicitPhasePressure) { - const volScalarField& rAU1 = mesh_.lookupObject - ( - IOobject::groupName("rAU", phase1_.name()) - ); - const volScalarField& rAU2 = mesh_.lookupObject - ( - IOobject::groupName("rAU", phase2_.name()) - ); + if (faceMomentum) + { + const surfaceScalarField& rAU1f = + mesh_.lookupObject + ( + IOobject::groupName("rAUf", phase1_.name()) + ); + const surfaceScalarField& rAU2f = + mesh_.lookupObject + ( + IOobject::groupName("rAUf", phase2_.name()) + ); - pPrimeByA = - fvc::interpolate(rAU1*phase1_.turbulence().pPrime()) - + fvc::interpolate(rAU2*phase2_.turbulence().pPrime()); + volScalarField D(this->D()); + + pPrimeByA = + rAU1f*fvc::interpolate(D + phase1_.turbulence().pPrime()) + + rAU2f*fvc::interpolate(D + phase2_.turbulence().pPrime()); + } + else + { + const volScalarField& rAU1 = mesh_.lookupObject + ( + IOobject::groupName("rAU", phase1_.name()) + ); + const volScalarField& rAU2 = mesh_.lookupObject + ( + IOobject::groupName("rAU", phase2_.name()) + ); + + pPrimeByA = + fvc::interpolate(rAU1*phase1_.turbulence().pPrime()) + + fvc::interpolate(rAU2*phase2_.turbulence().pPrime()); + } surfaceScalarField phiP ( @@ -596,8 +636,21 @@ bool Foam::twoPhaseSystem::read() } -const Foam::dragModel& -Foam::twoPhaseSystem::drag(const phaseModel& phase) const +const Foam::dimensionedScalar& +Foam::twoPhaseSystem::residualAlpha(const phaseModel& phase) const +{ + if (drag_->hasModel(phase)) + { + return drag_->phaseModel(phase).residualAlpha(); + } + else + { + return zeroResidualAlpha_; + } +} + + +const Foam::dragModel& Foam::twoPhaseSystem::drag(const phaseModel& phase) const { return drag_->phaseModel(phase); } diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H index c0d3facc7c..974af111c1 100644 --- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H +++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H @@ -40,13 +40,13 @@ SourceFiles #include "orderedPhasePair.H" #include "volFields.H" #include "surfaceFields.H" +#include "dragModel.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { -class dragModel; class virtualMassModel; class heatTransferModel; class liftModel; @@ -113,6 +113,9 @@ class twoPhaseSystem autoPtr > turbulentDispersion_; + //- + static dimensionedScalar zeroResidualAlpha_; + // Private member functions @@ -141,26 +144,29 @@ public: tmp U() const; //- Return the drag coefficient - tmp dragCoeff() const; + tmp Kd() const; + + //- Return the face drag coefficient + tmp Kdf() const; //- Return the virtual mass coefficient - tmp virtualMassCoeff() const; + tmp Vm() const; + + //- Return the face virtual mass coefficient + tmp Vmf() const; //- Return the heat transfer coefficient - tmp heatTransferCoeff() const; + tmp Kh() const; - //- Return the lift force - tmp liftForce() const; + //- Return the combined force (lift + wall-lubrication) + tmp F() const; - //- Return the wall lubrication force - tmp wallLubricationForce() const; + //- Return the combined face-force (lift + wall-lubrication) + tmp Ff() const; //- Return the turbulent diffusivity // Multiplies the phase-fraction gradient - tmp turbulentDiffusivity() const; - - //- Return the turbulent dispersion force - tmp turbulentDispersionForce() const; + tmp D() const; //- Solve for the two-phase-fractions void solve(); @@ -176,10 +182,17 @@ public: // Access - //- Return the drag model for the supplied phase + //- Return the residual phase-fraction for given phase + // Used to stabilize the phase momentum as the phase-fraction -> 0 + const dimensionedScalar& residualAlpha + ( + const phaseModel& phase + ) const; + + //- Return the drag model for the given phase const dragModel& drag(const phaseModel& phase) const; - //- Return the virtual mass model for the supplied phase + //- Return the virtual mass model for the given phase const virtualMassModel& virtualMass(const phaseModel& phase) const; //- Return the surface tension coefficient