From fc6b44ee3cd8ff7ea89b8edb4f4470247612fffb Mon Sep 17 00:00:00 2001 From: Henry Date: Mon, 27 Apr 2015 21:33:58 +0100 Subject: [PATCH] twoPhaseEulerFoam: Added experimental face-based momentum equation formulation This formulation provides C-grid like pressure-flux staggering on an unstructured mesh which is hugely beneficial for Euler-Euler multiphase equations as it allows for all forces to be treated in a consistent manner on the cell-faces which provides better balance, stability and accuracy. However, to achieve face-force consistency the momentum transport terms must be interpolated to the faces reducing accuracy of this part of the system but this is offset by the increase in accuracy of the force-balance. Currently it is not clear if this face-based momentum equation formulation is preferable for all Euler-Euler simulations so I have included it on a switch to allow evaluation and comparison with the previous cell-based formulation. To try the new algorithm simply switch it on, e.g.: PIMPLE { nOuterCorrectors 3; nCorrectors 1; nNonOrthogonalCorrectors 0; faceMomentum yes; } It is proving particularly good for bubbly flows, eliminating the staggering patterns often seen in the air velocity field with the previous algorithm, removing other spurious numerical artifacts in the velocity fields and improving stability and allowing larger time-steps For particle-gas flows the advantage is noticeable but not nearly as pronounced as in the bubbly flow cases. Please test the new algorithm on your cases and provide feedback. Henry G. Weller CFD Direct --- .../multiphase/twoPhaseEulerFoam/EEqns.H | 14 +- .../twoPhaseEulerFoam/createFields.H | 103 +++--- .../dragModels/dragModel/dragModel.C | 20 +- .../dragModels/dragModel/dragModel.H | 22 +- .../dragModels/segregated/segregated.C | 7 + .../dragModels/segregated/segregated.H | 5 +- .../liftModels/liftModel/liftModel.C | 21 +- .../liftModels/liftModel/liftModel.H | 12 +- .../liftModels/noLift/noLift.C | 77 +++- .../liftModels/noLift/noLift.H | 5 +- .../constantVirtualMassCoefficient.C | 28 +- .../virtualMassModel/virtualMassModel.C | 18 +- .../virtualMassModel/virtualMassModel.H | 17 +- .../wallLubricationModels/Antal/Antal.C | 3 +- .../wallLubricationModels/Antal/Antal.H | 6 +- .../wallLubricationModels/Frank/Frank.C | 3 +- .../wallLubricationModels/Frank/Frank.H | 6 +- .../TomiyamaWallLubrication.C | 3 +- .../TomiyamaWallLubrication.H | 6 +- .../noWallLubrication/noWallLubrication.C | 58 ++- .../noWallLubrication/noWallLubrication.H | 5 +- .../wallLubricationModel.C | 25 +- .../wallLubricationModel.H | 12 +- .../twoPhaseEulerFoam/{ => pU}/DDtU.H | 0 .../twoPhaseEulerFoam/{ => pU}/UEqns.H | 16 +- .../twoPhaseEulerFoam/pU/createDDtU.H | 17 + .../twoPhaseEulerFoam/{ => pU}/pEqn.H | 29 +- .../multiphase/twoPhaseEulerFoam/pUf/DDtU.H | 9 + .../multiphase/twoPhaseEulerFoam/pUf/UEqns.H | 50 +++ .../twoPhaseEulerFoam/pUf/createDDtU.H | 9 + .../multiphase/twoPhaseEulerFoam/pUf/pEqn.H | 347 ++++++++++++++++++ .../twoPhaseEulerFoam/twoPhaseEulerFoam.C | 25 +- .../BlendedInterfacialModel.C | 178 ++++++++- .../BlendedInterfacialModel.H | 22 +- .../twoPhaseSystem/twoPhaseSystem.C | 115 ++++-- .../twoPhaseSystem/twoPhaseSystem.H | 41 ++- 36 files changed, 1098 insertions(+), 236 deletions(-) rename applications/solvers/multiphase/twoPhaseEulerFoam/{ => pU}/DDtU.H (100%) rename applications/solvers/multiphase/twoPhaseEulerFoam/{ => pU}/UEqns.H (78%) create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/pU/createDDtU.H rename applications/solvers/multiphase/twoPhaseEulerFoam/{ => pU}/pEqn.H (91%) create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/pUf/DDtU.H create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/pUf/UEqns.H create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/pUf/createDDtU.H create mode 100644 applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H 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