diff --git a/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/dispersedDragModel/dispersedDragModel.H b/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/dispersedDragModel/dispersedDragModel.H index b38d4179d3..78ac88075e 100644 --- a/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/dispersedDragModel/dispersedDragModel.H +++ b/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/dispersedDragModel/dispersedDragModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -44,11 +44,14 @@ SourceFiles namespace Foam { + +class swarmCorrection; + namespace dragModels { /*---------------------------------------------------------------------------*\ - Class dispersedDragModel Declaration + Class dispersedDragModel Declaration \*---------------------------------------------------------------------------*/ class dispersedDragModel diff --git a/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/dragModel/dragModel.H b/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/dragModel/dragModel.H index 737812d064..267b4cda85 100644 --- a/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/dragModel/dragModel.H +++ b/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/dragModel/dragModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -46,8 +46,6 @@ SourceFiles namespace Foam { -class swarmCorrection; - /*---------------------------------------------------------------------------*\ Class dragModel Declaration \*---------------------------------------------------------------------------*/ @@ -83,8 +81,8 @@ public: //- Coefficient dimensions static const dimensionSet dimK; - //- Does this model require correcting on fixed flux boundaries? - static const bool correctFixedFluxBCs = true; + //- This model MUST NOT be set to 0 on fixed flux boundaries + static const bool correctFixedFluxBCs = false; // Constructors diff --git a/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/segregated/segregated.C b/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/segregated/segregated.C index 70674b2710..ca6ca45c86 100644 --- a/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/segregated/segregated.C +++ b/applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/segregated/segregated.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2014-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -72,29 +72,16 @@ Foam::tmp Foam::dragModels::segregated::K() const const volScalarField& alpha1(interface_.phase1()); const volScalarField& alpha2(interface_.phase2()); - const volScalarField& rho1(interface_.phase1().rho()); - const volScalarField& rho2(interface_.phase2().rho()); + const volScalarField::Internal& rho1(interface_.phase1().rho()); + const volScalarField::Internal& rho2(interface_.phase2().rho()); tmp tnu1(interface_.phase1().thermo().nu()); tmp tnu2(interface_.phase2().thermo().nu()); - const volScalarField& nu1(tnu1()); - const volScalarField& nu2(tnu2()); + const volScalarField::Internal& nu1(tnu1()); + const volScalarField::Internal& nu2(tnu2()); - volScalarField L - ( - IOobject - ( - "L", - mesh.time().name(), - mesh - ), - mesh, - dimensionedScalar(dimLength, 0), - zeroGradientFvPatchField::typeName - ); - L.primitiveFieldRef() = cbrt(mesh.V()); - L.correctBoundaryConditions(); + const volScalarField::Internal L(cbrt(mesh.V())); const dimensionedScalar residualAlpha ( @@ -104,44 +91,67 @@ Foam::tmp Foam::dragModels::segregated::K() const )/2 ); - const volScalarField I1(alpha1/max(alpha1 + alpha2, residualAlpha)); - const volScalarField I2(alpha2/max(alpha1 + alpha2, residualAlpha)); - const volScalarField magGradI + const volScalarField I1 + ( + alpha1/max(alpha1 + alpha2, residualAlpha) + ); + const volScalarField I2 + ( + alpha2/max(alpha1 + alpha2, residualAlpha) + ); + const volScalarField::Internal magGradI ( max ( - (rho2*mag(fvc::grad(I1)) + rho1*mag(fvc::grad(I2)))/(rho1 + rho2), + ( + rho2*mag(fvc::grad(I1)()()) + + rho1*mag(fvc::grad(I2)()()) + )/(rho1 + rho2), residualAlpha/2/L ) ); - const volScalarField muI(rho1*nu1*rho2*nu2/(rho1*nu1 + rho2*nu2)); + const volScalarField::Internal muI(rho1*nu1*rho2*nu2/(rho1*nu1 + rho2*nu2)); - const volScalarField limitedAlpha1 + const volScalarField::Internal limitedAlpha1 ( max(alpha1, interface_.phase1().residualAlpha()) ); - const volScalarField limitedAlpha2 + const volScalarField::Internal limitedAlpha2 ( max(alpha2, interface_.phase2().residualAlpha()) ); - const volScalarField muAlphaI + const volScalarField::Internal muAlphaI ( limitedAlpha1*rho1*nu1*limitedAlpha2*rho2*nu2 /(limitedAlpha1*rho1*nu1 + limitedAlpha2*rho2*nu2) ); - const volScalarField ReI + const volScalarField::Internal ReI ( - interface_.rho()*interface_.magUr() + (interface_.rho()()()*interface_.magUr()()()) /(magGradI*limitedAlpha1*limitedAlpha2*muI) ); - const volScalarField lambda(m_*ReI + n_*muAlphaI/muI); + const volScalarField::Internal lambda(m_*ReI + n_*muAlphaI/muI); - return lambda*sqr(magGradI)*muI; + tmp tK + ( + volScalarField::New + ( + "K", + mesh, + dimensionedScalar(dimK, 0), + zeroGradientFvPatchField::typeName + ) + ); + + tK.ref().ref() = lambda*sqr(magGradI)*muI; + tK.ref().correctBoundaryConditions(); + + return tK; } diff --git a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/dispersedVirtualMassModel/dispersedVirtualMassModel.C b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/dispersedVirtualMassModel/dispersedVirtualMassModel.C index b27e669fe8..867107ddb3 100644 --- a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/dispersedVirtualMassModel/dispersedVirtualMassModel.C +++ b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/dispersedVirtualMassModel/dispersedVirtualMassModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -58,14 +58,12 @@ Foam::virtualMassModels::dispersedVirtualMassModel::Ki() const Foam::tmp Foam::virtualMassModels::dispersedVirtualMassModel::K() const { - return interface_.dispersed()*Ki(); -} - - -Foam::tmp -Foam::virtualMassModels::dispersedVirtualMassModel::Kf() const -{ - return fvc::interpolate(interface_.dispersed())*fvc::interpolate(Ki()); + return + max + ( + interface_.dispersed(), + interface_.dispersed().residualAlpha() + )*Ki(); } diff --git a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/dispersedVirtualMassModel/dispersedVirtualMassModel.H b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/dispersedVirtualMassModel/dispersedVirtualMassModel.H index df18b8da0b..2ef3a45ad0 100644 --- a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/dispersedVirtualMassModel/dispersedVirtualMassModel.H +++ b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/dispersedVirtualMassModel/dispersedVirtualMassModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -96,10 +96,6 @@ public: // 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; }; diff --git a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/noVirtualMass/noVirtualMass.C b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/noVirtualMass/noVirtualMass.C index 9afb7f0fcb..272360da2e 100644 --- a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/noVirtualMass/noVirtualMass.C +++ b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/noVirtualMass/noVirtualMass.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2014-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -72,16 +72,4 @@ Foam::virtualMassModels::noVirtualMass::K() const } -Foam::tmp -Foam::virtualMassModels::noVirtualMass::Kf() const -{ - return surfaceScalarField::New - ( - "Kf", - interface_.mesh(), - dimensionedScalar(dimK, 0) - ); -} - - // ************************************************************************* // diff --git a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/noVirtualMass/noVirtualMass.H b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/noVirtualMass/noVirtualMass.H index f4d0f32b14..c2cb36686d 100644 --- a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/noVirtualMass/noVirtualMass.H +++ b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/noVirtualMass/noVirtualMass.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2014-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -85,10 +85,6 @@ public: //- Return the virtual mass coefficient K // used in the momentum equations virtual tmp K() const; - - //- Return the virtual mass coefficient Kf - // used in the face-momentum equations - virtual tmp Kf() const; }; diff --git a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.C b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.C index 83570ddd05..bfda371122 100644 --- a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.C +++ b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2014-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -81,10 +81,4 @@ Foam::tmp Foam::blendedVirtualMassModel::K() const } -Foam::tmp Foam::blendedVirtualMassModel::Kf() const -{ - return evaluate(&virtualMassModel::Kf, "Kf", virtualMassModel::dimK, false); -} - - // ************************************************************************* // diff --git a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.H b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.H index 627368f6dd..2ce96ca2ae 100644 --- a/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.H +++ b/applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2014-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -81,8 +81,8 @@ public: //- Coefficient dimensions static const dimensionSet dimK; - //- Does this model require correcting on fixed flux boundaries? - static const bool correctFixedFluxBCs = true; + //- This model should not be set to 0 on fixed flux boundaries + static const bool correctFixedFluxBCs = false; // Constructors @@ -119,10 +119,6 @@ public: // ddt(alpha2*rho2*U2) + ... = ... K*(DU1_Dt - DU2_Dt) virtual tmp K() const = 0; - //- Return the virtual mass coefficient Kf - // used in the face-momentum equations - virtual tmp Kf() const = 0; - // Dummy write for regIOobject bool writeData(Ostream& os) const; }; @@ -159,9 +155,6 @@ public: //- Return the lift coefficient K tmp K() const; - - //- Return the lift coefficient Kf - tmp Kf() const; }; diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C index c14e1f5e89..d2c02f836e 100644 --- a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C +++ b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2022-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -58,30 +58,39 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() rAUs.clear(); rAUs.setSize(phases.size()); + PtrList rAUfs(phases.size()); + + PtrList Kds(fluid.Kds()); + forAll(fluid.movingPhases(), movingPhasei) { phaseModel& phase = fluid.movingPhases()[movingPhasei]; const volScalarField& alpha = phase; + const volScalarField AU + ( + UEqns[phase.index()].A() + Kds[phase.index()] + + byDt + ( + max(phase.residualAlpha() - alpha, scalar(0)) + *phase.rho() + ) + ); + rAUs.set ( phase.index(), new volScalarField ( IOobject::groupName("rAU", phase.name()), - 1.0 - /( - UEqns[phase.index()].A() - + byDt - ( - max(phase.residualAlpha() - alpha, scalar(0)) - *phase.rho() - ) - ) + 1/AU ) ); + + rAUfs.set(phase.index(), 1/fvc::interpolate(AU)); } fluid.fillFields("rAU", dimTime/dimDensity, rAUs); + fluid.fillFields("rAUf", dimTime/dimDensity, rAUfs); // Phase diagonal coefficients PtrList alpharAUfs(phases.size()); @@ -94,13 +103,14 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() ( phasei, ( - fvc::interpolate(max(alpha, phase.residualAlpha())*rAUs[phasei]) + fvc::interpolate(max(alpha, phase.residualAlpha())) + *rAUfs[phasei] ).ptr() ); } // Explicit force fluxes - PtrList phiFs(fluid.phiFs(rAUs)); + PtrList Fs(fluid.Fs()); // Mass transfer rates PtrList dmdts(fluid.dmdts()); @@ -134,8 +144,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() ( phasei, ( - alpharAUfs[phasei] - *( + ( ghSnGradRho - (fvc::interpolate(phase.rho() - rho)) *(buoyancy.g & mesh.Sf()) @@ -143,11 +152,6 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() ) ).ptr() ); - - if (phiFs.set(phasei)) - { - phigFs[phasei] += phiFs[phasei]; - } } } @@ -156,21 +160,22 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() PtrList phiHbyAs(phases.size()); { // Correction force fluxes - PtrList ddtCorrByAs(fluid.ddtCorrByAs(rAUs)); + PtrList ddtCorrs(fluid.ddtCorrs()); forAll(fluid.movingPhases(), movingPhasei) { phaseModel& phase = fluid.movingPhases()[movingPhasei]; const volScalarField& alpha = phase; + const label phasei = phase.index(); HbyAs.set ( - phase.index(), + phasei, constrainHbyA ( - rAUs[phase.index()] + rAUs[phasei] *( - UEqns[phase.index()].H() + UEqns[phasei].H() + byDt ( max(phase.residualAlpha() - alpha, scalar(0)) @@ -185,13 +190,17 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() phiHbyAs.set ( - phase.index(), + phasei, new surfaceScalarField ( IOobject::groupName("phiHbyA", phase.name()), - fvc::flux(HbyAs[phase.index()]) - - phigFs[phase.index()] - - ddtCorrByAs[phase.index()] + rAUfs[phasei] + *( + fvc::flux(HbyAs[phasei]/rAUs[phasei]) + + ddtCorrs[phasei] + ) + - alpharAUfs[phasei]*phigFs[phasei] + - rAUfs[phasei]*Fs[phasei] ) ); } @@ -200,19 +209,19 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs); // Add explicit drag forces and fluxes - PtrList KdUByAs(fluid.KdUByAs(rAUs)); - PtrList phiKdPhis(fluid.phiKdPhis(rAUs)); + PtrList KdUs(fluid.KdUs()); + PtrList KdPhis(fluid.KdPhis()); forAll(phases, phasei) { - if (KdUByAs.set(phasei)) + if (KdUs.set(phasei)) { - HbyAs[phasei] -= KdUByAs[phasei]; + HbyAs[phasei] -= rAUs[phasei]*KdUs[phasei]; } - if (phiKdPhis.set(phasei)) + if (KdPhis.set(phasei)) { - phiHbyAs[phasei] -= phiKdPhis[phasei]; + phiHbyAs[phasei] -= rAUfs[phasei]*KdPhis[phasei]; } } @@ -255,6 +264,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() forAll(phases, phasei) { rAUf += alphafs[phasei]*alpharAUfs[phasei]; + // rAUf += alphafs[phasei]*alphafs[phasei]*rAUfs[phasei]; } rAUf = mag(rAUf); @@ -346,22 +356,54 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector() mSfGradp = pEqnIncomp.flux()/rAUf; - forAll(fluid.movingPhases(), movingPhasei) + if (!dragCorrection) { - phaseModel& phase = fluid.movingPhases()[movingPhasei]; + forAll(fluid.movingPhases(), movingPhasei) + { + phaseModel& phase = fluid.movingPhases()[movingPhasei]; + const label phasei = phase.index(); - phase.URef() = - HbyAs[phase.index()] - + fvc::reconstruct - ( - alpharAUfs[phase.index()]*mSfGradp - - phigFs[phase.index()] - ); + phase.URef() = + HbyAs[phasei] + + fvc::reconstruct + ( + alpharAUfs[phasei]*(mSfGradp - phigFs[phasei]) + - rAUfs[phasei]*Fs[phasei] + ); + } + } + else + { + PtrList dragCorrs(phases.size()); + PtrList dragCorrfs(phases.size()); + fluid.dragCorrs(dragCorrs, dragCorrfs); + + forAll(fluid.movingPhases(), movingPhasei) + { + phaseModel& phase = fluid.movingPhases()[movingPhasei]; + const label phasei = phase.index(); + + phase.URef() = + HbyAs[phasei] + + fvc::reconstruct + ( + alpharAUfs[phasei]*(mSfGradp - phigFs[phasei]) + + rAUfs[phasei]*(dragCorrfs[phasei] - Fs[phasei]) + ) + - rAUs[phasei]*dragCorrs[phasei]; + } } if (partialElimination) { - fluid.partialElimination(rAUs, KdUByAs, alphafs, phiKdPhis); + fluid.partialElimination + ( + rAUs, + KdUs, + alphafs, + rAUfs, + KdPhis + ); } else { diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C index c7059c851a..8067a4eca3 100644 --- a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C +++ b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2022-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -71,7 +71,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() rAUfs.clear(); rAUfs.setSize(phases.size()); { - PtrList AFfs(fluid.AFfs()); + PtrList KdVmfs(fluid.KdVmfs()); forAll(fluid.movingPhases(), movingPhasei) { @@ -87,7 +87,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() /( byDt(alphaRho0fs[phase.index()]) + fvc::interpolate(UEqns[phase.index()].A()) - + AFfs[phase.index()] + + KdVmfs[phase.index()] ) ) ); @@ -111,7 +111,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() } // Explicit force fluxes - PtrList phiFfs(fluid.phiFfs(rAUfs)); + PtrList Ffs(fluid.Ffs()); // Mass transfer rates PtrList dmdts(fluid.dmdts()); @@ -155,9 +155,9 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() ).ptr() ); - if (phiFfs.set(phasei)) + if (Ffs.set(phasei)) { - phigFs[phasei] += phiFfs[phasei]; + phigFs[phasei] += rAUfs[phasei]*Ffs[phasei]; } } } @@ -193,13 +193,13 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs); // Add explicit drag forces and fluxes - PtrList phiKdPhifs(fluid.phiKdPhifs(rAUfs)); + PtrList KdPhifs(fluid.KdPhifs()); forAll(phases, phasei) { - if (phiKdPhifs.set(phasei)) + if (KdPhifs.set(phasei)) { - phiHbyAs[phasei] -= phiKdPhifs[phasei]; + phiHbyAs[phasei] -= rAUfs[phasei]*KdPhifs[phasei]; } } @@ -330,7 +330,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector() if (partialElimination) { - fluid.partialEliminationf(rAUfs, alphafs, phiKdPhifs); + fluid.partialEliminationf(rAUfs, alphafs, KdPhifs); } else { diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/momentumPredictor.C b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/momentumPredictor.C index 97b1439212..a59e1b47f3 100644 --- a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/momentumPredictor.C +++ b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/momentumPredictor.C @@ -24,6 +24,7 @@ License \*---------------------------------------------------------------------------*/ #include "multiphaseEuler.H" +#include "fvmSup.H" // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // @@ -37,6 +38,8 @@ void Foam::solvers::multiphaseEuler::cellMomentumPredictor() phaseSystem::momentumTransferTable& momentumTransfer(momentumTransferPtr()); + const PtrList Kds(fluid.Kds()); + forAll(fluid.movingPhases(), movingPhasei) { phaseModel& phase = fluid.movingPhases()[movingPhasei]; @@ -54,6 +57,7 @@ void Foam::solvers::multiphaseEuler::cellMomentumPredictor() == *momentumTransfer[phase.name()] + fvModels().source(alpha, rho, U) + // - fvm::Sp(Kds[phase.index()], U) ) ); diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C index 9cfeb86a0e..c653cd3388 100644 --- a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C +++ b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C @@ -52,6 +52,9 @@ void Foam::solvers::multiphaseEuler::readControls() faceMomentum = pimple.dict().lookupOrDefault("faceMomentum", false); + dragCorrection = + pimple.dict().lookupOrDefault("dragCorrection", false); + partialElimination = pimple.dict().lookupOrDefault("partialElimination", false); @@ -97,6 +100,11 @@ Foam::solvers::multiphaseEuler::multiphaseEuler(fvMesh& mesh) pimple.dict().lookupOrDefault("faceMomentum", false) ), + dragCorrection + ( + pimple.dict().lookupOrDefault("dragCorrection", false) + ), + partialElimination ( pimple.dict().lookupOrDefault("partialElimination", false) diff --git a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H index 1ab0c4593b..a213c6f8b3 100644 --- a/applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H +++ b/applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H @@ -81,6 +81,10 @@ protected: // Defaults to false, i.e. uses the cell momentum equation Switch faceMomentum; + //- Cell/face drag correction for cell momentum corrector + // Defaults to false + Switch dragCorrection; + //- Partial elimination drag contribution optimisation // Defaults to false Switch partialElimination; diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C index 8bc482af64..bbcca7429a 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C @@ -31,18 +31,18 @@ License #include "wallLubricationModel.H" #include "turbulentDispersionModel.H" -#include "HashPtrTable.H" - #include "fvmDdt.H" #include "fvmDiv.H" #include "fvmSup.H" -#include "fvcAverage.H" + #include "fvcDdt.H" #include "fvcDiv.H" #include "fvcFlux.H" #include "fvcSnGrad.H" #include "fvcMeshPhi.H" -#include "fvMatrix.H" +#include "fvcReconstruct.H" + +#include "pimpleNoLoopControl.H" // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // @@ -197,30 +197,6 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() *Kds_[dragModelIter.key()] = dragModelIter()->K(); } - // Add the implicit part of the drag force - forAllConstIter(KdTable, Kds_, KdIter) - { - const volScalarField& K(*KdIter()); - const phaseInterface interface(*this, KdIter.key()); - - forAllConstIter(phaseInterface, interface, iter) - { - if (!iter().stationary()) - { - fvVectorMatrix& eqn = *eqns[iter().name()]; - - const volScalarField phaseK - ( - iter.otherPhase() - /max(iter.otherPhase(), iter.otherPhase().residualAlpha()) - *K - ); - - eqn -= fvm::Sp(phaseK, eqn.psi()); - } - } - } - // Initialise Vms table if required if (!Vms_.size()) { @@ -283,14 +259,20 @@ Foam::MomentumTransferPhaseSystem::momentumTransfer() const tmp taphi(fvc::absolute(phi, U)); const surfaceScalarField& aphi(taphi()); + const volScalarField VmPhase + ( + (otherPhase/max(otherPhase, otherPhase.residualAlpha())) + *Vm + ); + eqn -= - Vm + VmPhase *( fvm::ddt(U) + fvm::div(aphi, U) - fvm::Sp(fvc::div(aphi), U) - otherPhase.DUDt() ) - + this->MRF_.DDt(Vm, U - otherPhase.U()); + + this->MRF_.DDt(VmPhase, U - otherPhase.U()); } } } @@ -363,7 +345,7 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() *Kdfs_[dragModelIter.key()] = dragModelIter()->Kf(); } - // Initialise Vms and Vmfs tables if required + // Initialise Vms table if required if (!Vms_.size()) { forAllConstIter @@ -391,22 +373,6 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() dimensionedScalar(virtualMassModel::dimK, 0) ) ); - - Vmfs_.insert - ( - interface, - new surfaceScalarField - ( - IOobject - ( - IOobject::groupName("Vmf", interface.name()), - this->mesh().time().name(), - this->mesh() - ), - this->mesh(), - dimensionedScalar(virtualMassModel::dimK, 0) - ) - ); } } @@ -419,7 +385,6 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() ) { *Vms_[virtualMassModelIter.key()] = virtualMassModelIter()->K(); - *Vmfs_[virtualMassModelIter.key()] = virtualMassModelIter()->Kf(); } // Create U & grad(U) fields @@ -460,8 +425,14 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() if (!phase.stationary()) { + const volScalarField VmPhase + ( + (otherPhase/max(otherPhase, otherPhase.residualAlpha())) + *Vm + ); + *eqns[phase.name()] -= - Vm + VmPhase *( UgradUs[phase.index()] - (UgradUs[otherPhase.index()] & otherPhase.U()) @@ -476,9 +447,9 @@ Foam::MomentumTransferPhaseSystem::momentumTransferf() template Foam::PtrList -Foam::MomentumTransferPhaseSystem::AFfs() const +Foam::MomentumTransferPhaseSystem::KdVmfs() const { - PtrList AFfs(this->phaseModels_.size()); + PtrList KdVmfs(this->phaseModels_.size()); // Add the implicit part of the drag force forAllConstIter(KdfTable, Kdfs_, KdfIter) @@ -488,48 +459,57 @@ Foam::MomentumTransferPhaseSystem::AFfs() const forAllConstIter(phaseInterface, interface, iter) { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + const surfaceScalarField alphaf ( - fvc::interpolate(iter.otherPhase()) + fvc::interpolate(otherPhase) ); addField ( - iter(), - "AFf", - alphaf/max(alphaf, iter.otherPhase().residualAlpha()) + phase, + "KdVmf", + (alphaf/max(alphaf, otherPhase.residualAlpha())) *Kf, - AFfs + KdVmfs ); } } // Add the implicit part of the virtual mass force - forAllConstIter(VmfTable, Vmfs_, VmfIter) + forAllConstIter(VmTable, Vms_, VmIter) { - const surfaceScalarField& Vmf(*VmfIter()); - const phaseInterface interface(*this, VmfIter.key()); + const volScalarField& Vm(*VmIter()); + const phaseInterface interface(*this, VmIter.key()); forAllConstIter(phaseInterface, interface, iter) { - addField(iter(), "AFf", byDt(Vmf), AFfs); + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + const volScalarField VmPhase + ( + (otherPhase/max(otherPhase, otherPhase.residualAlpha())) + *Vm + ); + + addField(phase, "KdVmf", byDt(fvc::interpolate(VmPhase)), KdVmfs); } } - this->fillFields("AFf", dimDensity/dimTime, AFfs); + this->fillFields("KdVmf", dimDensity/dimTime, KdVmfs); - return AFfs; + return KdVmfs; } template Foam::PtrList -Foam::MomentumTransferPhaseSystem::phiFs -( - const PtrList& rAUs -) +Foam::MomentumTransferPhaseSystem::Fs() { - PtrList phiFs(this->phaseModels_.size()); + PtrList Fs(this->phaseModels_.size()); // Add the lift force forAllConstIter @@ -546,16 +526,16 @@ Foam::MomentumTransferPhaseSystem::phiFs addField ( interface.phase1(), - "phiF", - fvc::flux(rAUs[interface.phase1().index()]*F), - phiFs + "F", + fvc::flux(F), + Fs ); addField ( interface.phase2(), - "phiF", - -fvc::flux(rAUs[interface.phase2().index()]*F), - phiFs + "F", + -fvc::flux(F), + Fs ); } @@ -575,16 +555,16 @@ Foam::MomentumTransferPhaseSystem::phiFs addField ( interface.phase1(), - "phiF", - fvc::flux(rAUs[interface.phase1().index()]*F), - phiFs + "F", + fvc::flux(F), + Fs ); addField ( interface.phase2(), - "phiF", - -fvc::flux(rAUs[interface.phase2().index()]*F), - phiFs + "F", + -fvc::flux(F), + Fs ); } @@ -598,10 +578,10 @@ Foam::MomentumTransferPhaseSystem::phiFs addField ( phase, - "phiF", - fvc::interpolate(rAUs[phasei]*pPrime(), pPrime().name()) + "F", + fvc::interpolate(pPrime(), pPrime().name()) *fvc::snGrad(phase)*this->mesh_.magSf(), - phiFs + Fs ); } @@ -617,15 +597,7 @@ Foam::MomentumTransferPhaseSystem::phiFs turbulentDispersionModelIter()->interface(); const volScalarField D(turbulentDispersionModelIter()->D()); - - const surfaceScalarField DByA1f - ( - fvc::interpolate(rAUs[interface.phase1().index()]*D) - ); - const surfaceScalarField DByA2f - ( - fvc::interpolate(rAUs[interface.phase2().index()]*D) - ); + const surfaceScalarField Df(fvc::interpolate(D)); const volScalarField alpha12(interface.phase1() + interface.phase2()); const surfaceScalarField snGradAlpha1By12 @@ -645,41 +617,47 @@ Foam::MomentumTransferPhaseSystem::phiFs )*this->mesh_.magSf() ); - addField(interface.phase1(), "phiF", DByA1f*snGradAlpha1By12, phiFs); - addField(interface.phase2(), "phiF", DByA2f*snGradAlpha2By12, phiFs); + addField(interface.phase1(), "F", Df*snGradAlpha1By12, Fs); + addField(interface.phase2(), "F", Df*snGradAlpha2By12, Fs); } if (this->fillFields_) { - this->fillFields("phiF", dimForce/dimDensity/dimVelocity, phiFs); + this->fillFields("F", dimArea*dimDensity*dimAcceleration, Fs); } - return phiFs; + return Fs; } template Foam::PtrList -Foam::MomentumTransferPhaseSystem::phiFfs -( - const PtrList& rAUfs -) +Foam::MomentumTransferPhaseSystem::Ffs() { - PtrList phiFfs(this->phaseModels_.size()); + PtrList Ffs(this->phaseModels_.size()); // Add the explicit part of the virtual mass force - forAllConstIter(VmfTable, Vmfs_, VmfIter) + forAllConstIter(VmTable, Vms_, VmIter) { - const surfaceScalarField& Vmf(*VmfIter()); - const phaseInterface interface(*this, VmfIter.key()); + const volScalarField& Vm(*VmIter()); + const phaseInterface interface(*this, VmIter.key()); forAllConstIter(phaseInterface, interface, iter) { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + const volScalarField VmPhase + ( + (otherPhase/max(otherPhase, otherPhase.residualAlpha())) + *Vm + ); + addField ( - iter(), - "phiFf", - -rAUfs[iter().index()]*Vmf + phase, + "Ff", + -fvc::interpolate(VmPhase) *( byDt ( @@ -689,9 +667,9 @@ Foam::MomentumTransferPhaseSystem::phiFfs iter().U() ) ) - + iter.otherPhase().DUDtf() + + otherPhase.DUDtf() ), - phiFfs + Ffs ); } } @@ -711,16 +689,16 @@ Foam::MomentumTransferPhaseSystem::phiFfs addField ( interface.phase1(), - "phiFf", - rAUfs[interface.phase1().index()]*Ff, - phiFfs + "Ff", + Ff, + Ffs ); addField ( interface.phase2(), - "phiFf", - -rAUfs[interface.phase2().index()]*Ff, - phiFfs + "Ff", + -Ff, + Ffs ); } @@ -740,16 +718,16 @@ Foam::MomentumTransferPhaseSystem::phiFfs addField ( interface.phase1(), - "phiFf", - rAUfs[interface.phase1().index()]*Ff, - phiFfs + "Ff", + Ff, + Ffs ); addField ( interface.phase2(), - "phiFf", - -rAUfs[interface.phase2().index()]*Ff, - phiFfs + "Ff", + -Ff, + Ffs ); } @@ -761,10 +739,10 @@ Foam::MomentumTransferPhaseSystem::phiFfs addField ( phase, - "phiFf", - rAUfs[phasei]*fvc::interpolate(phase.pPrime()) + "Ff", + fvc::interpolate(phase.pPrime()) *fvc::snGrad(phase)*this->mesh_.magSf(), - phiFfs + Ffs ); } @@ -784,9 +762,6 @@ Foam::MomentumTransferPhaseSystem::phiFfs fvc::interpolate(turbulentDispersionModelIter()->D()) ); - const surfaceScalarField DByA1f(rAUfs[interface.phase1().index()]*Df); - const surfaceScalarField DByA2f(rAUfs[interface.phase2().index()]*Df); - const volScalarField alpha12(interface.phase1() + interface.phase2()); const surfaceScalarField snGradAlpha1By12 ( @@ -805,27 +780,24 @@ Foam::MomentumTransferPhaseSystem::phiFfs )*this->mesh_.magSf() ); - addField(interface.phase1(), "phiF", DByA1f*snGradAlpha1By12, phiFfs); - addField(interface.phase2(), "phiF", DByA2f*snGradAlpha2By12, phiFfs); + addField(interface.phase1(), "F", Df*snGradAlpha1By12, Ffs); + addField(interface.phase2(), "F", Df*snGradAlpha2By12, Ffs); } if (this->fillFields_) { - this->fillFields("phiFf", dimForce/dimDensity/dimVelocity, phiFfs); + this->fillFields("Ff", dimArea*dimDensity*dimAcceleration, Ffs); } - return phiFfs; + return Ffs; } template Foam::PtrList -Foam::MomentumTransferPhaseSystem::phiKdPhis -( - const PtrList& rAUs -) const +Foam::MomentumTransferPhaseSystem::KdPhis() const { - PtrList phiKdPhis(this->phaseModels_.size()); + PtrList KdPhis(this->phaseModels_.size()); // Add the explicit part of the drag force forAllConstIter(KdTable, Kds_, KdIter) @@ -835,22 +807,24 @@ Foam::MomentumTransferPhaseSystem::phiKdPhis forAllConstIter(phaseInterface, interface, iter) { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + addField ( - iter(), - "phiKdPhi", + phase, + "KdPhi", fvc::interpolate ( - -iter.otherPhase() - /max(iter.otherPhase(), iter.otherPhase().residualAlpha()) - *rAUs[iter().index()]*K + -(otherPhase/max(otherPhase, otherPhase.residualAlpha())) + *K ) *fvc::absolute ( - this->MRF().absolute(iter.otherPhase().phi()), - iter.otherPhase().U() + this->MRF().absolute(otherPhase.phi()), + otherPhase.U() ), - phiKdPhis + KdPhis ); } } @@ -859,24 +833,21 @@ Foam::MomentumTransferPhaseSystem::phiKdPhis { this->fillFields ( - "phiKdPhi", - dimForce/dimDensity/dimVelocity, - phiKdPhis + "KdPhi", + dimArea*dimDensity*dimAcceleration, + KdPhis ); } - return phiKdPhis; + return KdPhis; } template Foam::PtrList -Foam::MomentumTransferPhaseSystem::phiKdPhifs -( - const PtrList& rAUfs -) const +Foam::MomentumTransferPhaseSystem::KdPhifs() const { - PtrList phiKdPhifs(this->phaseModels_.size()); + PtrList KdPhifs(this->phaseModels_.size()); // Add the explicit part of the drag force forAllConstIter(KdfTable, Kdfs_, KdfIter) @@ -886,24 +857,26 @@ Foam::MomentumTransferPhaseSystem::phiKdPhifs forAllConstIter(phaseInterface, interface, iter) { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + const surfaceScalarField alphaf ( - fvc::interpolate(iter.otherPhase()) + fvc::interpolate(otherPhase) ); addField ( - iter(), - "phiKdPhif", - -rAUfs[iter().index()] - *alphaf/max(alphaf, iter.otherPhase().residualAlpha()) + phase, + "KdPhif", + -(alphaf/max(alphaf, otherPhase.residualAlpha())) *Kf *fvc::absolute ( - this->MRF().absolute(iter.otherPhase().phi()), - iter.otherPhase().U() + this->MRF().absolute(otherPhase.phi()), + otherPhase.U() ), - phiKdPhifs + KdPhifs ); } } @@ -912,24 +885,57 @@ Foam::MomentumTransferPhaseSystem::phiKdPhifs { this->fillFields ( - "phiKdPhif", - dimForce/dimDensity/dimVelocity, - phiKdPhifs + "KdPhif", + dimArea*dimDensity*dimAcceleration, + KdPhifs ); } - return phiKdPhifs; + return KdPhifs; +} + + +template +Foam::PtrList +Foam::MomentumTransferPhaseSystem::Kds() const +{ + PtrList Kds(this->phaseModels_.size()); + + forAllConstIter(KdTable, Kds_, KdIter) + { + const volScalarField& K(*KdIter()); + const phaseInterface interface(*this, KdIter.key()); + + forAllConstIter(phaseInterface, interface, iter) + { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + addField + ( + phase, + "Kd", + (otherPhase/max(otherPhase, otherPhase.residualAlpha())) + *K, + Kds + ); + } + } + + if (this->fillFields_) + { + this->fillFields("Kd", dimDensity/dimTime, Kds); + } + + return Kds; } template Foam::PtrList -Foam::MomentumTransferPhaseSystem::KdUByAs -( - const PtrList& rAUs -) const +Foam::MomentumTransferPhaseSystem::KdUs() const { - PtrList KdUByAs(this->phaseModels_.size()); + PtrList KdUs(this->phaseModels_.size()); // Add the explicit part of the drag force forAllConstIter(KdTable, Kds_, KdIter) @@ -939,25 +945,26 @@ Foam::MomentumTransferPhaseSystem::KdUByAs forAllConstIter(phaseInterface, interface, iter) { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + addField ( - iter(), - "KdUByA", - -rAUs[iter().index()] - *iter.otherPhase() - /max(iter.otherPhase(), iter.otherPhase().residualAlpha()) - *K*iter.otherPhase().U(), - KdUByAs + phase, + "KdU", + -(otherPhase/max(otherPhase, otherPhase.residualAlpha())) + *K*otherPhase.U(), + KdUs ); } } if (this->fillFields_) { - this->fillFields("KdUByA", dimVelocity, KdUByAs); + this->fillFields("KdU", dimDensity*dimAcceleration, KdUs); } - return KdUByAs; + return KdUs; } @@ -1082,76 +1089,70 @@ Foam::MomentumTransferPhaseSystem::alphaDByAf template Foam::PtrList -Foam::MomentumTransferPhaseSystem::ddtCorrByAs +Foam::MomentumTransferPhaseSystem::ddtCorrs ( - const PtrList& rAUs, const bool includeVirtualMass ) const { - PtrList ddtCorrByAs(this->phaseModels_.size()); - - // Construct phi differences - PtrList phiCorrs(this->phaseModels_.size()); - - forAll(this->movingPhases(), movingPhasei) - { - const phaseModel& phase = this->movingPhases()[movingPhasei]; - const label phasei = phase.index(); - - phiCorrs.set - ( - phasei, - ( - phase.Uf().valid() - ? (this->mesh_.Sf() & phase.Uf()().oldTime())() - : phase.phi()().oldTime() - ) - - fvc::flux(phase.U()().oldTime()) - ); - } + PtrList ddtCorrs(this->phaseModels_.size()); // Add correction forAll(this->movingPhases(), movingPhasei) { const phaseModel& phase = this->movingPhases()[movingPhasei]; - const label phasei = phase.index(); - const volScalarField& alpha = phase; - - // Apply ddtPhiCorr filter in pure(ish) phases - surfaceScalarField alphafBar - ( - fvc::interpolate(fvc::average(fvc::interpolate(alpha))) - ); - - tmp phiCorrCoeff = pos0(alphafBar - 0.99); - - surfaceScalarField::Boundary& phiCorrCoeffBf = - phiCorrCoeff.ref().boundaryFieldRef(); - - forAll(this->mesh_.boundary(), patchi) - { - // Set ddtPhiCorr to 0 on non-coupled boundaries - if (!this->mesh_.boundary()[patchi].coupled()) - { - phiCorrCoeffBf[patchi] = 0; - } - } addField ( phase, - "ddtCorrByA", - -phiCorrCoeff*phiCorrs[phasei]*fvc::interpolate + "ddtCorr", + fvc::ddtCorr ( - byDt(alpha.oldTime()*phase.rho()().oldTime()*rAUs[phasei]) + phase, + phase.rho()(), + phase.U()(), + phase.phi()(), + phase.Uf() ), - ddtCorrByAs + ddtCorrs ); } - // Add virtual mass correction - if (includeVirtualMass) + + const pimpleNoLoopControl& pimple = this->pimple(); + const Switch VmDdtCorr + ( + pimple.dict().lookupOrDefault("VmDdtCorrection", false) + ); + + // Optionally ddd virtual mass correction + if (VmDdtCorr) { + PtrList VmDdtCoeffs(this->phaseModels_.size()); + PtrList VmDdtCorrs(this->phaseModels_.size()); + + forAll(this->movingPhases(), movingPhasei) + { + const phaseModel& phase = this->movingPhases()[movingPhasei]; + const label phasei = phase.index(); + + VmDdtCoeffs.set + ( + phasei, + fvm::ddt(phase.U()())().A() + ); + + VmDdtCorrs.set + ( + phasei, + fvc::ddtCorr + ( + phase.U()(), + phase.phi()(), + phase.Uf() + ) + ); + } + forAllConstIter(VmTable, Vms_, VmIter) { const volScalarField& Vm(*VmIter()); @@ -1160,30 +1161,146 @@ Foam::MomentumTransferPhaseSystem::ddtCorrByAs forAllConstIter(phaseInterface, interface, iter) { const phaseModel& phase = iter(); + const label phasei = phase.index(); const phaseModel& otherPhase = iter.otherPhase(); + const label otherPhasei = otherPhase.index(); + + const volScalarField VmPhase + ( + (otherPhase/max(otherPhase, otherPhase.residualAlpha())) + *Vm + ); addField ( - iter(), - "ddtCorrByA", - -fvc::interpolate(Vm*byDt(rAUs[phase.index()])) + phase, + "ddtCorr", + fvc::interpolate(VmPhase) *( - phiCorrs[phase.index()] - + fvc::absolute - ( - this->MRF().absolute(otherPhase.phi()), - otherPhase.U() - ) - - fvc::flux(otherPhase.U()) - - phiCorrs[otherPhase.index()] + VmDdtCorrs[phasei] + + ( + fvc::interpolate(VmDdtCoeffs[otherPhasei]) + *( + otherPhase.Uf().valid() + ? (this->mesh_.Sf() & otherPhase.Uf()())() + : otherPhase.phi()() + ) + - fvc::flux(VmDdtCoeffs[otherPhasei]*otherPhase.U()) + ) + - VmDdtCorrs[otherPhase.index()] ), - ddtCorrByAs + ddtCorrs ); } } } - return ddtCorrByAs; + return ddtCorrs; +} + + +template +void Foam::MomentumTransferPhaseSystem::dragCorrs +( + PtrList& dragCorrs, + PtrList& dragCorrfs +) const +{ + const phaseSystem::phaseModelList& phases = this->phaseModels_; + + PtrList Uphis(phases.size()); + + forAll(phases, i) + { + if (!phases[i].stationary()) + { + Uphis.set + ( + i, + fvc::reconstruct(phases[i].phi()) + ); + } + } + + forAllConstIter(KdTable, Kds_, KdIter) + { + const volScalarField& K(*KdIter()); + const phaseInterface interface(*this, KdIter.key()); + + const phaseModel& phase1 = interface.phase1(); + const phaseModel& phase2 = interface.phase2(); + + const label phase1i = phase1.index(); + const label phase2i = phase2.index(); + + if (!phase1.stationary()) + { + const volScalarField K1 + ( + (phase2/max(phase2, phase2.residualAlpha()))*K + ); + + addField + ( + phase1, + "dragCorr", + K1 + *( + phase2.stationary() + ? -Uphis[phase1i] + : (Uphis[phase2i] - Uphis[phase1i]) + ), + dragCorrs + ); + + addField + ( + phase1, + "dragCorrf", + fvc::interpolate(K1) + *( + phase2.stationary() + ? -phase1.phi() + : (phase2.phi() - phase1.phi()) + ), + dragCorrfs + ); + } + + if (!phase2.stationary()) + { + const volScalarField K2 + ( + (phase1/max(phase1, phase1.residualAlpha()))*K + ); + + addField + ( + phase2, + "dragCorr", + K2 + *( + phase1.stationary() + ? -Uphis[phase2i] + : (Uphis[phase1i] - Uphis[phase2i]) + ), + dragCorrs + ); + + addField + ( + phase2, + "dragCorrf", + fvc::interpolate(K2) + *( + phase1.stationary() + ? -phase2.phi() + : (phase1.phi() - phase2.phi()) + ), + dragCorrfs + ); + } + } } @@ -1191,9 +1308,10 @@ template void Foam::MomentumTransferPhaseSystem::partialElimination ( const PtrList& rAUs, - const PtrList& KdUByAs, + const PtrList& KdUs, const PtrList& alphafs, - const PtrList& phiKdPhis + const PtrList& rAUfs, + const PtrList& KdPhis ) { Info<< "Inverting drag systems: "; @@ -1222,15 +1340,15 @@ void Foam::MomentumTransferPhaseSystem::partialElimination { if (!phases[i].stationary()) { - phases[i].URef() += KdUByAs[i]; - phases[i].phiRef() += phiKdPhis[i]; + phases[i].URef() += rAUs[i]*KdUs[i]; + phases[i].phiRef() += rAUfs[i]*KdPhis[i]; } } { // Create drag coefficient matrices PtrList> KdByAs(phases.size()); - PtrList> phiKds(phases.size()); + PtrList> KdByAfs(phases.size()); forAll(phases, i) { @@ -1240,7 +1358,7 @@ void Foam::MomentumTransferPhaseSystem::partialElimination new PtrList(phases.size()) ); - phiKds.set + KdByAfs.set ( i, new PtrList(phases.size()) @@ -1255,50 +1373,58 @@ void Foam::MomentumTransferPhaseSystem::partialElimination const label phase1i = interface.phase1().index(); const label phase2i = interface.phase2().index(); + const volScalarField K1 + ( + interface.phase2() + /max(interface.phase2(), interface.phase2().residualAlpha()) + *K + ); + + const volScalarField K2 + ( + interface.phase1() + /max(interface.phase1(), interface.phase1().residualAlpha()) + *K + ); + addField ( interface.phase2(), "KdByA", - -rAUs[phase1i] - *interface.phase2() - /max(interface.phase2(), interface.phase2().residualAlpha()) - *K, + -rAUs[phase1i]*K1, KdByAs[phase1i] ); addField ( interface.phase1(), "KdByA", - -rAUs[phase2i] - *interface.phase1() - /max(interface.phase1(), interface.phase1().residualAlpha()) - *K, + -rAUs[phase2i]*K2, KdByAs[phase2i] ); addField ( interface.phase2(), - "phiKd", - fvc::interpolate(KdByAs[phase1i][phase2i]), - phiKds[phase1i] + "KdByAf", + -rAUfs[phase1i]*fvc::interpolate(K1), + KdByAfs[phase1i] ); addField ( interface.phase1(), - "phiKd", - fvc::interpolate(KdByAs[phase2i][phase1i]), - phiKds[phase2i] + "KdByAf", + -rAUfs[phase2i]*fvc::interpolate(K2), + KdByAfs[phase2i] ); } forAll(phases, i) { this->fillFields("KdByAs", dimless, KdByAs[i]); - this->fillFields("phiKds", dimless, phiKds[i]); + this->fillFields("KdByAfs", dimless, KdByAfs[i]); KdByAs[i][i] = 1; - phiKds[i][i] = 1; + KdByAfs[i][i] = 1; } // Decompose @@ -1307,23 +1433,23 @@ void Foam::MomentumTransferPhaseSystem::partialElimination for (label j = i + 1; j < phases.size(); j++) { KdByAs[i][j] /= KdByAs[i][i]; - phiKds[i][j] /= phiKds[i][i]; + KdByAfs[i][j] /= KdByAfs[i][i]; for (label k = i + 1; k < phases.size(); ++ k) { KdByAs[j][k] -= KdByAs[j][i]*KdByAs[i][k]; - phiKds[j][k] -= phiKds[j][i]*phiKds[i][k]; + KdByAfs[j][k] -= KdByAfs[j][i]*KdByAfs[i][k]; } } } { volScalarField detKdByAs(KdByAs[0][0]); - surfaceScalarField detPhiKdfs(phiKds[0][0]); + surfaceScalarField detPhiKdfs(KdByAfs[0][0]); for (label i = 1; i < phases.size(); i++) { detKdByAs *= KdByAs[i][i]; - detPhiKdfs *= phiKds[i][i]; + detPhiKdfs *= KdByAfs[i][i]; } Info<< "Min cell/face det = " << gMin(detKdByAs.primitiveField()) @@ -1338,7 +1464,7 @@ void Foam::MomentumTransferPhaseSystem::partialElimination for (label j = 0; j < i; j ++) { phases[i].URef() -= KdByAs[i][j]*phases[j].U(); - phases[i].phiRef() -= phiKds[i][j]*phases[j].phi(); + phases[i].phiRef() -= KdByAfs[i][j]*phases[j].phi(); } } } @@ -1349,10 +1475,10 @@ void Foam::MomentumTransferPhaseSystem::partialElimination for (label j = phases.size() - 1; j > i; j--) { phases[i].URef() -= KdByAs[i][j]*phases[j].U(); - phases[i].phiRef() -= phiKds[i][j]*phases[j].phi(); + phases[i].phiRef() -= KdByAfs[i][j]*phases[j].phi(); } phases[i].URef() /= KdByAs[i][i]; - phases[i].phiRef() /= phiKds[i][i]; + phases[i].phiRef() /= KdByAfs[i][i]; } } } @@ -1367,7 +1493,7 @@ void Foam::MomentumTransferPhaseSystem::partialEliminationf ( const PtrList& rAUfs, const PtrList& alphafs, - const PtrList& phiKdPhifs + const PtrList& KdPhifs ) { Info<< "Inverting drag system: "; @@ -1380,7 +1506,7 @@ void Foam::MomentumTransferPhaseSystem::partialEliminationf { if (!phases[i].stationary()) { - phases[i].phiRef() += phiKdPhifs[i]; + phases[i].phiRef() += rAUfs[i]*KdPhifs[i]; } } diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H index afd965ad99..d3f4d0595e 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2015-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2015-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -63,8 +63,6 @@ class MomentumTransferPhaseSystem : public BasePhaseSystem { -private: - // Private typedefs typedef HashPtrTable @@ -88,13 +86,6 @@ private: phaseInterfaceKey::hash > VmTable; - typedef HashPtrTable - < - surfaceScalarField, - phaseInterfaceKey, - phaseInterfaceKey::hash - > VmfTable; - typedef HashTable < autoPtr, @@ -142,9 +133,6 @@ private: //- Virtual mass coefficients VmTable Vms_; - //- Face virtual mass coefficients - VmfTable Vmfs_; - // Sub Models @@ -206,72 +194,33 @@ public: //- Return implicit force coefficients on the faces, for the face-based // algorithm. - virtual PtrList AFfs() const; + virtual PtrList KdVmfs() const; //- Return the explicit force fluxes for the cell-based algorithm, that // do not depend on phase mass/volume fluxes, and can therefore be // evaluated outside the corrector loop. This includes things like // lift, turbulent dispersion, and wall lubrication. - virtual PtrList phiFs - ( - const PtrList& rAUs - ); + virtual PtrList Fs(); - //- As phiFs, but for the face-based algorithm - virtual PtrList phiFfs - ( - const PtrList& rAUfs - ); + //- As Fs, but for the face-based algorithm + virtual PtrList Ffs(); //- Return the explicit drag force fluxes for the cell-based algorithm. // These depend on phase mass/volume fluxes, and must therefore be // evaluated inside the corrector loop. - virtual PtrList phiKdPhis - ( - const PtrList& rAUs - ) const; + virtual PtrList KdPhis() const; - //- As phiKdPhis, but for the face-based algorithm - virtual PtrList phiKdPhifs - ( - const PtrList& rAUfs - ) const; + //- As KdPhis, but for the face-based algorithm + virtual PtrList KdPhifs() const; + + //- Return the implicit part of the drag force + virtual PtrList Kds() const; //- Return the explicit part of the drag force for the cell-based - // algorithm. This is the cell-equivalent of phiKdPhis. These depend on + // algorithm. This is the cell-equivalent of KdPhis. These depend on // phase velocities, and must therefore be evaluated inside the // corrector loop. - virtual PtrList KdUByAs - ( - const PtrList& rAUs - ) const; - - //- Solve the drag system for the velocities and fluxes - virtual void partialElimination - ( - const PtrList& rAUs, - const PtrList& KdUByAs, - const PtrList& alphafs, - const PtrList& phiKdPhis - ); - - //- As partialElimination, but for the face-based algorithm. Only solves - // for the fluxes. - virtual void partialEliminationf - ( - const PtrList& rAUfs, - const PtrList& alphafs, - const PtrList& phiKdPhifs - ); - - //- Return the flux corrections for the cell-based algorithm. These - // depend on phase mass/volume fluxes, and must therefore be evaluated - // inside the corrector loop. - virtual PtrList ddtCorrByAs - ( - const PtrList& rAUs, - const bool includeVirtualMass = false - ) const; + virtual PtrList KdUs() const; //- Returns true if the phase pressure is treated implicitly // in the phase fraction equation @@ -289,6 +238,40 @@ public: const PtrList& rAUfs ) const; + //- Return the flux corrections for the cell-based algorithm. These + // depend on phase mass/volume fluxes, and must therefore be evaluated + // inside the corrector loop. + virtual PtrList ddtCorrs + ( + const bool includeVirtualMass = false + ) const; + + //- Set the cell and faces drag correction fields + virtual void dragCorrs + ( + PtrList& dragCorrs, + PtrList& dragCorrf + ) const; + + //- Solve the drag system for the velocities and fluxes + virtual void partialElimination + ( + const PtrList& rAUs, + const PtrList& KdUs, + const PtrList& alphafs, + const PtrList& rAUfs, + const PtrList& KdPhis + ); + + //- As partialElimination, but for the face-based algorithm. Only solves + // for the fluxes. + virtual void partialEliminationf + ( + const PtrList& rAUfs, + const PtrList& alphafs, + const PtrList& KdPhifs + ); + //- Read base phaseProperties dictionary virtual bool read(); }; diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/diameterModels/velocityGroup/velocityGroup.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/diameterModels/velocityGroup/velocityGroup.C index d45e511a39..52770a136b 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/diameterModels/velocityGroup/velocityGroup.C +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/diameterModels/velocityGroup/velocityGroup.C @@ -212,7 +212,7 @@ void Foam::diameterModels::velocityGroup::correct() const populationBalanceModel& popBal = phase().mesh().lookupObject(popBalName_); - if (!popBal.solveOnFinalIterOnly() || popBal.pimple().finalIter()) + if (!popBal.solveOnFinalIterOnly() || popBal.fluid().pimple().finalIter()) { forAll(sizeGroups_, i) { diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C index dabbdc669c..d866bdedde 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C @@ -141,9 +141,7 @@ Foam::MovingPhaseModel::MovingPhaseModel ( IOobject::groupName("alphaPhi", this->name()), fluid.mesh().time().name(), - fluid.mesh(), - IOobject::READ_IF_PRESENT, - IOobject::NO_WRITE + fluid.mesh() ), fluid.mesh(), dimensionedScalar(dimensionSet(0, 3, -1, 0, 0), 0) @@ -154,9 +152,7 @@ Foam::MovingPhaseModel::MovingPhaseModel ( IOobject::groupName("alphaRhoPhi", this->name()), fluid.mesh().time().name(), - fluid.mesh(), - IOobject::READ_IF_PRESENT, - IOobject::AUTO_WRITE + fluid.mesh() ), fluid.mesh(), dimensionedScalar(dimensionSet(1, 0, -1, 0, 0), 0) diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C index a3f09fbd66..ff5e65c3a0 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C @@ -225,6 +225,8 @@ Foam::phaseSystem::phaseSystem mesh_(mesh), + pimple_(mesh_.lookupObject("solutionControl")), + MRF_(mesh_), referencePhaseName_(lookupOrDefault("referencePhase", word::null)), diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H index 09a00702cb..821c86e577 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H @@ -45,6 +45,8 @@ SourceFiles #include "PtrListDictionary.H" #include "hashedWordList.H" +#include "pimpleNoLoopControl.H" + #include "IOMRFZoneList.H" #include "fvModels.H" #include "fvConstraints.H" @@ -126,6 +128,9 @@ protected: //- Reference to the mesh const fvMesh& mesh_; + //- Reference to pimpleNoLoopControl + const pimpleNoLoopControl& pimple_; + //- Optional MRF zones IOMRFZoneList MRF_; @@ -273,6 +278,9 @@ public: //- Return the mesh inline const fvMesh& mesh() const; + //- Return pimpleNoLoopControl + inline const pimpleNoLoopControl& pimple() const; + //- Return the phase models inline const phaseModelList& phases() const; @@ -515,37 +523,25 @@ public: //- Return the implicit force coefficients for the face-based // algorithm - virtual PtrList AFfs() const = 0; + virtual PtrList KdVmfs() const = 0; //- Return the force fluxes for the cell-based algorithm - virtual PtrList phiFs - ( - const PtrList& rAUs - ) = 0; + virtual PtrList Fs() = 0; //- Return the force fluxes for the face-based algorithm - virtual PtrList phiFfs - ( - const PtrList& rAUfs - ) = 0; + virtual PtrList Ffs() = 0; //- Return the force fluxes for the cell-based algorithm - virtual PtrList phiKdPhis - ( - const PtrList& rAUs - ) const = 0; + virtual PtrList KdPhis() const = 0; //- Return the force fluxes for the face-based algorithm - virtual PtrList phiKdPhifs - ( - const PtrList& rAUfs - ) const = 0; + virtual PtrList KdPhifs() const = 0; + + //- Return the implicit part of the drag force + virtual PtrList Kds() const = 0; //- Return the explicit part of the drag force - virtual PtrList KdUByAs - ( - const PtrList& rAUs - ) const = 0; + virtual PtrList KdUs() const = 0; //- Returns true if the phase pressure is treated implicitly // in the phase fraction equation @@ -563,13 +559,27 @@ public: const PtrList& rAUfs ) const = 0; + //- Return the flux corrections for the cell-based algorithm + virtual PtrList ddtCorrs + ( + const bool includeVirtualMass = false + ) const = 0; + + //- Set the cell and faces drag correction fields + virtual void dragCorrs + ( + PtrList& dragCorrs, + PtrList& dragCorrf + ) const = 0; + //- Solve the drag system for the new velocities and fluxes virtual void partialElimination ( const PtrList& rAUs, - const PtrList& KdUByAs, + const PtrList& KdUs, const PtrList& alphafs, - const PtrList& phiKdPhis + const PtrList& rAUfs, + const PtrList& KdPhis ) = 0; //- Solve the drag system for the new fluxes @@ -577,7 +587,7 @@ public: ( const PtrList& rAUfs, const PtrList& alphafs, - const PtrList& phiKdPhifs + const PtrList& KdPhifs ) = 0; //- Re-normalise the flux of the phases @@ -588,13 +598,6 @@ public: const surfaceScalarField& phim ); - //- Return the flux corrections for the cell-based algorithm - virtual PtrList ddtCorrByAs - ( - const PtrList& rAUs, - const bool includeVirtualMass = false - ) const = 0; - //- Return the heat transfer matrices virtual autoPtr heatTransfer() const = 0; diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemI.H b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemI.H index 1f41050559..bb8966b07e 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemI.H +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2014-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2014-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -31,6 +31,12 @@ inline const Foam::fvMesh& Foam::phaseSystem::mesh() const } +inline const Foam::pimpleNoLoopControl& Foam::phaseSystem::pimple() const +{ + return pimple_; +} + + inline const Foam::phaseSystem::phaseModelList& Foam::phaseSystem::phases() const { diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C index 09c324d462..1bbd8749c6 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C @@ -797,7 +797,6 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel ( fluid_.subDict("populationBalanceCoeffs").subDict(name_) ), - pimple_(mesh_.lookupObject("solutionControl")), continuousPhase_ ( mesh_.lookupObject @@ -1130,7 +1129,7 @@ Foam::diameterModels::populationBalanceModel::continuousTurbulence() const void Foam::diameterModels::populationBalanceModel::solve() { - if (!solveOnFinalIterOnly() || pimple_.finalIter()) + if (!solveOnFinalIterOnly() || fluid_.pimple().finalIter()) { const label nCorr = this->nCorr(); const scalar tolerance = diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H b/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H index f551cc5c37..87e3b05b92 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2017-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -206,7 +206,6 @@ SourceFiles #include "sizeGroup.H" #include "phaseSystem.H" -#include "pimpleNoLoopControl.H" #include "phaseCompressibleMomentumTransportModel.H" #include "HashPtrTable.H" #include "Pair.H" @@ -232,8 +231,6 @@ class populationBalanceModel : public regIOobject { -private: - // Private Data //- Reference to the phaseSystem @@ -251,9 +248,6 @@ private: //- Dictionary dictionary dict_; - //- Reference to pimpleNoLoopControl - const pimpleNoLoopControl& pimple_; - //- Continuous phase const phaseModel& continuousPhase_; @@ -438,9 +432,6 @@ public: //- Return populationBalanceCoeffs dictionary inline const dictionary& dict() const; - //- Return pimpleNoLoopControl - const pimpleNoLoopControl& pimple() const; - //- Return the number of corrections inline label nCorr() const; diff --git a/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H b/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H index 4a0afb691a..81efbe0a8c 100644 --- a/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H +++ b/applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2017-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2023 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -64,13 +64,6 @@ Foam::diameterModels::populationBalanceModel::dict() const } -inline const Foam::pimpleNoLoopControl& -Foam::diameterModels::populationBalanceModel::pimple() const -{ - return pimple_; -} - - inline Foam::label Foam::diameterModels::populationBalanceModel::nCorr() const { return mesh_.solution().solverDict(name_).lookup