From a8cb8a61dabddb54009131de901ece9e9b39e33f Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Thu, 30 Mar 2023 12:27:48 +0100 Subject: [PATCH] solvers::multiphaseEuler: New cell momentum/pressure algorithm The cell-base momentum/pressure algorithm in the multiphaseEuler solver module has been substantially updated to improve consistency, conservation and reduce drag generated staggering patterns at sharp interfaces and the boundaries with stationary phases. For most if not all cases this new algorithm can be used to provide well resolved and reliable solutions where the faceMomentum algorithm would have been chosen previously in order to obtain sufficiently smooth solutions but at the expense of a noticeable loss in accuracy and resolution. The first significant change in the momentum/pressure algorithm is in the interpolation practice used to construct the flux predictor equation from the cell momentum equation: rather than interpolating the H/A ratio to the faces i.e. (H/A)_f the terms in the momentum equation are interpolated separately so that H_f/A_f is used. The same approach is used for the drag i.e. (D_f/A_f) and virtual mass contributions. The advantage of this change is that the phase forces are now consistent in both the momentum and flux equations, i.e. sum to zero for each pair of phases. The second significant change is in the handling of ddtCorr which converts the old-time time-derivative contributions in H from velocity to flux which is now consistent due to the change to H/A interpolation and also generalised to use the fvc::ddtCorr function which has been updated for multiphase. Additionally ddtCorr may optionally be applied to the time-derivative in the virtual mass term in a consistent manner so that the contributions to the flux equation sum to zero for each pair of phases. The third significant change is the addition of an optional drag correction term to the momentum corrector to reduce the staggering patters generated in the velocity field due to sudden changes in drag force between phase, e.g. at sharp interfaces between phases or at the boundaries with stationary phases. This is particularly beneficial for fluidised bed simulations. However this correction is not and cannot be phase consistent, i.e. the correction does not sum to zero for pairs of phases it is applied to so a small drag error is introduced, but tests so far have shown that the error is small and outweighed by the benefit in the reduction in numerical artefacts in the solution. The final significant change is in the handling of residualAlpha for drag and virtual mass to provide stable and physical phase velocities in the limit of the phase-fraction -> 0. The new approach is phase asymmetric such that the residual drag is applied only to the phase with a phase-fraction less than residualAlpha and not to the carrier phase. This change ensures that the flow of a pure phase is unaffected by the residualAlpha and residual drag of the other phases that are stabilised in pure phase region. There are now four options in the PIMPLE section of the fvSolutions dictionary relating to the multiphase momentum/pressure algorithm: PIMPLE { faceMomentum no; VmDdtCorrection yes; dragCorrection yes; partialElimination no; } faceMomentum: Switches between the cell and face momentum equation algorithms. Provides much smoother and reliable solutions for even the most challenging multiphase cases at the expense of a noticeable loss in accuracy and resolution. Defaults to 'no'. VmDdtCorrection: Includes the ddtCorr correction term to the time-derivative part of the virtual-mass term in the flux equation which ensures consistency between the phase virtual mass force on the faces but generates solutions which are slightly less smooth and more likely to contain numerical artefacts. Defaults to 'no'. Testing so far has shown that the loss in smoothness is small and there is some noticeable improvement is some cases so in the future the default may be changed to 'yes'. dragCorrection: Includes the momentum corrector drag correction term to reduce the staggering patters generated in the velocity field due to sudden changes in drag force at the expense of a small error in drag consistency. Defaults to 'no' partialElimination: Switches the partial-elimination momentum corrector which inverts the drag matrix for both the momentum equations and/or flux equations to provide a drag implicit correction to the phase velocity and flux fields. The algorithm is the same as previously but updated for the new consistent drag interpolation. All the tutorials/modules/multiphaseEuler tutorial cases have been updated and tested with the above developments and the four options set appropriately for each. --- .../dispersedDragModel/dispersedDragModel.H | 7 +- .../dragModels/dragModel/dragModel.H | 8 +- .../dragModels/segregated/segregated.C | 72 +- .../dispersedVirtualMassModel.C | 16 +- .../dispersedVirtualMassModel.H | 6 +- .../noVirtualMass/noVirtualMass.C | 14 +- .../noVirtualMass/noVirtualMass.H | 6 +- .../virtualMassModel/virtualMassModel.C | 8 +- .../virtualMassModel/virtualMassModel.H | 13 +- .../multiphaseEuler/cellPressureCorrector.C | 128 ++-- .../multiphaseEuler/facePressureCorrector.C | 20 +- .../multiphaseEuler/momentumPredictor.C | 4 + .../multiphaseEuler/multiphaseEuler.C | 8 + .../multiphaseEuler/multiphaseEuler.H | 4 + .../MomentumTransferPhaseSystem.C | 672 +++++++++++------- .../MomentumTransferPhaseSystem.H | 111 ++- .../velocityGroup/velocityGroup.C | 2 +- .../MovingPhaseModel/MovingPhaseModel.C | 8 +- .../phaseSystems/phaseSystem/phaseSystem.C | 2 + .../phaseSystems/phaseSystem/phaseSystem.H | 65 +- .../phaseSystems/phaseSystem/phaseSystemI.H | 8 +- .../populationBalanceModel.C | 3 +- .../populationBalanceModel.H | 11 +- .../populationBalanceModelI.H | 9 +- .../Grossetete/system/fvSolution | 5 +- .../LBend/constant/phaseProperties | 8 +- .../multiphaseEuler/LBend/system/fvSolution | 6 +- .../multiphaseEuler/bubbleColumn/0/p_rgh | 2 +- .../bubbleColumn/system/fvSolution | 5 + .../bubbleColumnEvaporating/system/fvSolution | 5 + .../system/fvSolution | 5 + .../system/fvSolution | 5 + .../bubbleColumnIATE/system/fvSolution | 4 + .../bubbleColumnLES/system/fvSolution | 5 + .../bubbleColumnLaminar/system/fvSolution | 4 + .../bubblePipe/system/fvSolution | 6 +- .../damBreak4phase/system/controlDict | 6 +- .../damBreak4phase/system/fvSolution | 5 + .../fluidisedBed/system/fvSolution | 5 +- .../constant/phaseProperties | 10 +- .../fluidisedBedLaminar/system/fvSolution | 9 +- .../hydrofoil/system/fvSolution | 3 + .../injection/constant/phaseProperties | 11 +- .../injection/system/fvSolution | 4 + .../mixerVessel2D/constant/phaseProperties | 4 +- .../constant/physicalProperties.mercury | 2 +- .../constant/physicalProperties.oil | 2 +- .../constant/physicalProperties.water | 2 +- .../mixerVessel2D/system/fvSolution | 5 + .../mixerVessel2DMRF/system/fvSolution | 5 + .../pipeBend/system/fvSolution | 23 +- .../multiphaseEuler/steamInjection/0/T.steam | 12 +- .../multiphaseEuler/steamInjection/0/T.water | 12 +- .../multiphaseEuler/steamInjection/0/U.steam | 6 +- .../multiphaseEuler/steamInjection/0/U.water | 2 + .../steamInjection/0/alpha.steam | 8 +- .../steamInjection/0/alpha.water | 12 +- .../steamInjection/0/nut.water | 2 + .../multiphaseEuler/steamInjection/0/p | 2 + .../multiphaseEuler/steamInjection/0/p_rgh | 9 +- .../steamInjection/0/water.steam | 41 -- .../steamInjection/0/water.water | 41 -- .../steamInjection/constant/phaseProperties | 12 +- .../steamInjection/system/fvSolution | 49 +- .../titaniaSynthesis/system/fvSolution | 5 +- .../titaniaSynthesisSurface/system/fvSolution | 5 +- .../wallBoilingIATE/system/fvSolution | 6 +- .../wallBoilingPolydisperse/system/fvSolution | 6 +- .../system/fvSolution | 6 +- 69 files changed, 863 insertions(+), 734 deletions(-) delete mode 100644 tutorials/modules/multiphaseEuler/steamInjection/0/water.steam delete mode 100644 tutorials/modules/multiphaseEuler/steamInjection/0/water.water 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