From 5bdc07e9339078f0198979b657b8538fe4141164 Mon Sep 17 00:00:00 2001 From: "Schlegel, Dr. Fabian (FWDC) - 107795" Date: Tue, 8 Dec 2020 15:31:26 +0100 Subject: [PATCH] multiphaseEulerFoam: Add cipsaMultiphaseEulerFoam solver Solver, which makes use of the compact momentum interpolation method, according to: Cubero, A. and Fueyo, N. (2007). A compact momentum interpolation procedure for unsteady flows and relaxation. Numerical Heat Transfer, Part B: Fundamentals, 52(6):507-529. * A partial elimination sum formulation is added to approximately resolve phase coupling due to the large drag coefficient for an arbitrary number of phases * Interfacial momentum and mass exchanges are handled consistently * Necessary infrastructure in phaseSystems is added * New drag models are added: StrubeljTiselj, IshiiZuberSparse * A new blending method is added: definitePhases * ResolvedPhaseSystem added: calculation of surfaceTension and interfaceCompression * InterfaceCompression becomes a blendedInterfacialModel * InterfaceIndicatorMethods added for detection of resolved interfaces * Added SGSdeltas: copy of LESdeltas, which can be accessed by interfaceIndicatorMethods * Support of PID time step control for cipsaMultiphaseEulerFoam solver * Mixture density of phasePair: formulation for arbitrary number of phases --- .../multiphase/multiphaseEulerFoam/Allwclean | 1 + .../multiphase/multiphaseEulerFoam/Allwmake | 1 + .../cipsaMultiphaseEulerFoam/CourantNo.H | 60 +++ .../cipsaMultiphaseEulerFoam/EEqns.H | 61 +++ .../cipsaMultiphaseEulerFoam/Make/files | 3 + .../cipsaMultiphaseEulerFoam/Make/options | 25 ++ .../cipsaMultiphaseEulerFoam.C | 135 +++++++ .../createFieldRefs.H | 49 +++ .../cipsaMultiphaseEulerFoam/createFields.H | 69 ++++ .../initZeroPtrList.H | 82 +++++ .../cipsaMultiphaseEulerFoam/pU/HUEqn.H | 84 +++++ .../cipsaMultiphaseEulerFoam/pU/UEqns.H | 234 ++++++++++++ .../cipsaMultiphaseEulerFoam/pU/ddtUEqn.H | 47 +++ .../cipsaMultiphaseEulerFoam/pU/dmdtUEqn.H | 74 ++++ .../cipsaMultiphaseEulerFoam/pU/dragDiag.H | 71 ++++ .../cipsaMultiphaseEulerFoam/pU/dragSource.H | 75 ++++ .../cipsaMultiphaseEulerFoam/pU/dragSourcef.H | 84 +++++ .../pU/forcesSourcef.H | 36 ++ .../pU/fvOptionsEqn.H | 18 + .../cipsaMultiphaseEulerFoam/pU/pEqn.H | 347 ++++++++++++++++++ .../pU/stabilizeHUEqn.H | 43 +++ .../cipsaMultiphaseEulerFoam/pU/vmDiag.H | 16 + .../cipsaMultiphaseEulerFoam/pU/vmSource.H | 22 ++ .../cipsaMultiphaseEulerFoam/pU/vmSourcef.H | 23 ++ .../interfacialCompositionModels/Make/files | 5 + .../SGSdeltas/SGSdelta/SGSdelta.C | 93 +++++ .../SGSdeltas/SGSdelta/SGSdelta.H | 143 ++++++++ .../cubeRootVolDelta/cubeRootVolDelta.C | 126 +++++++ .../cubeRootVolDelta/cubeRootVolDelta.H | 119 ++++++ .../SGSdeltas/maxDeltaxyz/maxDeltaxyz.C | 131 +++++++ .../SGSdeltas/maxDeltaxyz/maxDeltaxyz.H | 124 +++++++ .../oneDimensionalDelta/oneDimensionalDelta.C | 125 +++++++ .../oneDimensionalDelta/oneDimensionalDelta.H | 112 ++++++ .../interfacialModels/Make/files | 2 + .../IshiiZuberSparse/IshiiZuberSparse.C | 83 +++++ .../IshiiZuberSparse/IshiiZuberSparse.H | 105 ++++++ .../StrubeljTiselj/StrubeljTiselj.C | 94 +++++ .../StrubeljTiselj/StrubeljTiselj.H | 121 ++++++ .../multiphaseSystems/multiphaseSystems.C | 26 +- .../definitePhases/definitePhases.C | 140 +++++++ .../definitePhases/definitePhases.H | 142 +++++++ .../phaseSystems/Make/files | 16 + .../MomentumTransferPhaseSystem.C | 118 ++++++ .../MomentumTransferPhaseSystem.H | 80 ++++ .../PhaseTransferPhaseSystem.H | 5 + .../ResolvedPhaseSystem/ResolvedPhaseSystem.C | 220 +++++++++++ .../ResolvedPhaseSystem/ResolvedPhaseSystem.H | 141 +++++++ .../constantInterfaceCompression.C | 91 +++++ .../constantInterfaceCompression.H | 100 +++++ .../interfaceCompressionModel.C | 102 +++++ .../interfaceCompressionModel.H | 138 +++++++ .../interfaceCompressionModelNew.C | 60 +++ .../noInterfaceCompression.C | 116 ++++++ .../noInterfaceCompression.H | 96 +++++ .../gradAlpha/gradAlpha.C | 107 ++++++ .../gradAlpha/gradAlpha.H | 142 +++++++ .../interfaceIndicatorMethod.C | 118 ++++++ .../interfaceIndicatorMethod.H | 137 +++++++ .../interfaceIndicatorMethodNew.C | 112 ++++++ .../useBlendingMethod/useBlendingMethod.C | 106 ++++++ .../useBlendingMethod/useBlendingMethod.H | 134 +++++++ .../continuumSurfaceForce.C | 241 ++++++++++++ .../continuumSurfaceForce.H | 119 ++++++ .../noSurfaceTensionForce.C | 102 +++++ .../noSurfaceTensionForce.H | 96 +++++ .../surfaceTensionForceModel.C | 82 +++++ .../surfaceTensionForceModel.H | 129 +++++++ .../surfaceTensionForceModelNew.C | 78 ++++ .../MovingPhaseModel/MovingPhaseModel.C | 8 + .../MovingPhaseModel/MovingPhaseModel.H | 3 + .../StationaryPhaseModel.C | 15 + .../StationaryPhaseModel.H | 3 + .../phaseModel/phaseModel/phaseModel.C | 9 +- .../phaseModel/phaseModel/phaseModel.H | 12 + .../phasePair/phasePair/phasePair.C | 19 +- .../phaseSystems/phaseSystem/phaseSystem.C | 303 +++++---------- .../phaseSystems/phaseSystem/phaseSystem.H | 180 ++++++--- .../phaseSystems/phaseSystem/phaseSystemI.H | 7 + .../phaseSystem/phaseSystemSolve.C | 90 +---- .../general/include/createPIDTimeControls.H | 60 +++ .../general/include/readPIDTimeControls.H | 91 +++++ .../general/include/setInitialPIDDeltaT.H | 53 +++ .../cfdTools/general/include/setPIDDeltaT.H | 82 +++++ 83 files changed, 6935 insertions(+), 337 deletions(-) create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/EEqns.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/files create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/options create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFields.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dmdtUEqn.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/fvOptionsEqn.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/maxDeltaxyz/maxDeltaxyz.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/maxDeltaxyz/maxDeltaxyz.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/IshiiZuberSparse/IshiiZuberSparse.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/dragModels/StrubeljTiselj/StrubeljTiselj.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/BlendedInterfacialModel/blendingMethods/definitePhases/definitePhases.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/ResolvedPhaseSystem.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/ResolvedPhaseSystem.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/constantInterfaceCompression/constantInterfaceCompression.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModel.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/interfaceCompressionModel/interfaceCompressionModelNew.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceCompressionModels/noInterfaceCompression/noInterfaceCompression.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/gradAlpha/gradAlpha.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethod.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/interfaceIndicatorMethod/interfaceIndicatorMethodNew.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/interfaceIndicatorMethods/useBlendingMethod/useBlendingMethod.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/continuumSurfaceForce/continuumSurfaceForce.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/noSurfaceTensionForce/noSurfaceTensionForce.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.C create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModel.H create mode 100644 applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/ResolvedPhaseSystem/surfaceTensionForceModels/surfaceTensionForceModel/surfaceTensionForceModelNew.C create mode 100644 src/finiteVolume/cfdTools/general/include/createPIDTimeControls.H create mode 100644 src/finiteVolume/cfdTools/general/include/readPIDTimeControls.H create mode 100644 src/finiteVolume/cfdTools/general/include/setInitialPIDDeltaT.H create mode 100644 src/finiteVolume/cfdTools/general/include/setPIDDeltaT.H diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean b/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean index b8fe60678..858f85f76 100755 --- a/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean +++ b/applications/solvers/multiphase/multiphaseEulerFoam/Allwclean @@ -7,6 +7,7 @@ wclean libso interfacialCompositionModels wclean libso multiphaseCompressibleMomentumTransportModels wclean libso multiphaseThermophysicalTransportModels multiphaseEulerFoam/Allwclean +wclean libso cipsaMultiphaseEulerFoam wclean libso functionObjects #------------------------------------------------------------------------------ diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake b/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake index db3b01550..9c4618de9 100755 --- a/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake +++ b/applications/solvers/multiphase/multiphaseEulerFoam/Allwmake @@ -10,6 +10,7 @@ wmake $targetType interfacialCompositionModels wmake $targetType multiphaseCompressibleMomentumTransportModels wmake $targetType multiphaseThermophysicalTransportModels multiphaseEulerFoam/Allwmake $targetType $* +wmake $targetType cipsaMultiphaseEulerFoam wmake $targetType functionObjects #------------------------------------------------------------------------------ diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H new file mode 100644 index 000000000..d1d28ab5c --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/CourantNo.H @@ -0,0 +1,60 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Global + CourantNo + +Description + Calculates and outputs the mean and maximum Courant Numbers. + +\*---------------------------------------------------------------------------*/ + +scalar CoNum = 0.0; +scalar meanCoNum = 0.0; + +if (mesh.nInternalFaces()) +{ + scalarField sumPhi + ( + fvc::surfaceSum(mag(phi))().primitiveField() + ); + + forAll(phases, phasei) + { + sumPhi = max + ( + sumPhi, + fvc::surfaceSum(mag(phases[phasei].phi()))().primitiveField() + ); + } + + CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); + + meanCoNum = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue(); +} + +Info<< "Courant Number mean: " << meanCoNum + << " max: " << CoNum << endl; + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/EEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/EEqns.H new file mode 100644 index 000000000..a6e4c5256 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/EEqns.H @@ -0,0 +1,61 @@ +for (int Ecorr=0; Ecorr + heatTransferPtr(fluid.heatTransfer()); + + phaseSystem::heatTransferTable& heatTransfer = heatTransferPtr(); + + forAll(fluid.anisothermalPhases(), anisothermalPhasei) + { + phaseModel& phase = fluid.anisothermalPhases()[anisothermalPhasei]; + + const volScalarField& alpha = phase; + tmp tRhoAlpha = phase.rho(); + const volScalarField& rhoAlpha = tRhoAlpha(); + tmp tUAlpha = phase.U(); + const volVectorField& UAlpha = tUAlpha(); + + fvScalarMatrix EEqn + ( + phase.heEqn() + == + *heatTransfer[phase.name()] + + alpha*rhoAlpha*(UAlpha & g) + + fvOptions(alpha, rhoAlpha, phase.thermoRef().he()) + ); + + EEqn += + fvm::ddt + ( + max(phase.residualAlpha() - alpha, Zero)*rhoAlpha, + phase.thermoRef().he() + ) + - fvc::ddt + ( + max(phase.residualAlpha() - alpha, Zero)*rhoAlpha, + phase.thermoRef().he() + ); + + EEqn.relax(); + fvOptions.constrain(EEqn); + EEqn.solve(); + fvOptions.correct(phase.thermoRef().he()); + } + + fluid.correctThermo(); + fluid.correctContinuityError(); +} + + +forAll(phases, phasei) +{ + phaseModel& phase = phases[phasei]; + + Info<< phase.name() << " min/max T " + << min(phase.thermo().T()).value() + << " - " + << max(phase.thermo().T()).value() + << endl; +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/files b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/files new file mode 100644 index 000000000..3299c2137 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/files @@ -0,0 +1,3 @@ +cipsaMultiphaseEulerFoam.C + +EXE = $(FOAM_APPBIN)/cipsaMultiphaseEulerFoam diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/options b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/options new file mode 100644 index 000000000..915392055 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/Make/options @@ -0,0 +1,25 @@ +EXE_INC = \ + -I../multiphaseEulerFoam/multiphaseSystems \ + -I../phaseSystems/lnInclude \ + -I../interfacialModels/lnInclude \ + -I../interfacialCompositionModels/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude + +EXE_LIBS = \ + -lphaseSystem \ + -lmultiphaseSystems \ + -leulerianInterfacialModels \ + -leulerianInterfacialCompositionModels \ + -lmultiphaseMomentumTransportModels \ + -lmultiphaseThermophysicalTransportModels \ + -lthermophysicalTransportModels \ + -lfiniteVolume \ + -lfvOptions \ + -lmeshTools \ + -lsampling diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C new file mode 100644 index 000000000..2a4bb8d72 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/cipsaMultiphaseEulerFoam.C @@ -0,0 +1,135 @@ +/*---------------------------------------------------------------------------*\ + == == ====== ==== ==== | + \\ || | Unsupported Contributions for OpenFOAM + ====== // || || ===// | + || || // || // || \\ | Copyright (C) 2018-2020 OpenFOAM Foundation + == == ====== ==== == == | and Helmholtz-Zentrum Dresden-Rossendorf +------------------------------------------------------------------------------- +License + This file is a derivative work of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Application + cipsaMultiphaseEulerFoam + +Description + Solver for a system of any number of incompressible fluid phases with a + common pressure, but otherwise separate properties. The solver uses the + consistent momentum interpolation method of Cubero et al. (2014) and a + partial-elimination algorithm. + + Reference: + \verbatim + Cubero, A., Sánchez-Insa, A., & Fueyo, N. (2014). + A consistent momentum interpolation method for steady and unsteady + multiphase flows. + Computers & Chemical Engineering, 62, 96-107. + \endverbatim + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "phaseSystem.H" +#include "phaseCompressibleMomentumTransportModel.H" +#include "pimpleControl.H" +#include "initZeroPtrList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + + #include "setRootCaseLists.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "createPIDTimeControls.H" + #include "createFields.H" + #include "createFieldRefs.H" + #include "CourantNo.H" + #include "setInitialPIDDeltaT.H" + + Switch faceMomentum + ( + pimple.dict().lookupOrDefault("faceMomentum", false) + ); + + Switch partialElimination + ( + pimple.dict().lookupOrDefault("partialElimination", true) + ); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (pimple.run(runTime)) + { + #include "readPIDTimeControls.H" + + int nEnergyCorrectors + ( + pimple.dict().lookupOrDefault("nEnergyCorrectors", 1) + ); + + #include "CourantNo.H" + #include "setPIDDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // --- Pressure-velocity corrector loop (PIMPLE) + while (pimple.loop()) + { + fluid.solve(rAUs, rHAUfs); + fluid.correct(); + // Update continuityError in MovingPhaseModel for heEqn in + // AnisothermalPhaseModel + fluid.correctContinuityError(); + + #include "pU/UEqns.H" + #include "EEqns.H" + + // --- Pressure corrector loop (PISO) + while (pimple.correct()) + { + #include "pU/pEqn.H" + } + + // Update material derivatives + fluid.correctKinematics(); + + if (pimple.turbCorr()) + { + fluid.correctTurbulence(); + } + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H new file mode 100644 index 000000000..1803f8d75 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFieldRefs.H @@ -0,0 +1,49 @@ +surfaceScalarField& phi = fluid.phi(); + +const IOMRFZoneList& MRF = fluid.MRF(); +const fv::options& fvOptions = fluid.fvOptions(); + +List relaxU(phases.size()); +PtrList alphafs(phases.size()); + +// Equations +PtrList HUEqns(phases.size()); +PtrList ddtUEqns(phases.size()); +PtrList dmdtUEqns(phases.size()); +PtrList fvOptionsEqns(phases.size()); +PtrList ddtPhiCoeffs(phases.size()); +PtrList dmdtPhiCoeffs(phases.size()); + +// Cell-based coefficients +PtrList rHUs(phases.size()); +PtrList rAs(phases.size()); +PtrList rAUs; // Dummy argument for solveAlphas() + +PtrList dragDiag(phases.size()); +PtrList dragSource(phases.size()); + +PtrList vmDiag(phases.size()); +PtrList vmSource(phases.size()); + +// Face-based coefficients +PtrList rHUfs(phases.size()); +PtrList rAfs(phases.size()); +PtrList rAUfs(phases.size()); +PtrList rHAUfs(phases.size()); // Argument for solveAlphas() +PtrList alpharAUfs(phases.size()); + +PtrList dragDiagf(phases.size()); +PtrList dragSourcef(phases.size()); +PtrList dragPCoefff(phases.size()); + +PtrList vmDiagf(phases.size()); +PtrList vmSourcef(phases.size()); + +PtrList forcesSourcef(phases.size()); + +// Cell velocities +PtrList HbyAs(phases.size()); + +// Face fluxes +PtrList phigs(phases.size()); +PtrList phiHbyAs(phases.size()); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFields.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFields.H new file mode 100644 index 000000000..d24413a6f --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/createFields.H @@ -0,0 +1,69 @@ +#include "readGravitationalAcceleration.H" +#include "readhRef.H" + +Info<< "Creating phaseSystem\n" << endl; + +autoPtr fluidPtr +( + phaseSystem::New(mesh) +); +phaseSystem& fluid = fluidPtr(); +phaseSystem::phaseModelList& phases = fluid.phases(); +phaseSystem::phaseModelPartialList& movingPhases = fluid.movingPhases(); + +// Create mixture density and velocity +volScalarField rho("rho", fluid.rho()); +volVectorField U("U", fluid.U()); + +dimensionedScalar pMin +( + "pMin", + dimPressure, + fluid +); + +#include "gh.H" + +volScalarField& p = phases[0].thermoRef().p(); + +Info<< "Reading field p_rgh\n" << endl; +volScalarField p_rgh +( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +label pRefCell = 0; +scalar pRefValue = 0.0; +if (fluid.incompressible()) +{ + p = max(p_rgh + fluid.rho()*gh, pMin); + + if (p_rgh.needReference()) + { + setRefCell + ( + p, + p_rgh, + pimple.dict(), + pRefCell, + pRefValue + ); + + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); + p_rgh = p - fluid.rho()*gh; + } +} +mesh.setFluxRequired(p_rgh.name()); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H new file mode 100644 index 000000000..c7816d470 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/initZeroPtrList.H @@ -0,0 +1,82 @@ +/*---------------------------------------------------------------------------*\ + == == ====== ==== ==== | + \\ || | Unsupported Contributions for OpenFOAM + ====== // || || ===// | + || || // || // || \\ | Copyright (C) 2018-2020 OpenFOAM Foundation + == == ====== ==== == == | and Helmholtz-Zentrum Dresden-Rossendorf +------------------------------------------------------------------------------- +License + This file is a derivative work of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +InNamespace + Foam + +Description + Initialize the fields of a PtrList with zero + +\*---------------------------------------------------------------------------*/ + +#ifndef initZeroPtrList_H +#define initZeroPtrList_H + +#include "volFieldsFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template class PatchField, class GeoMesh> +void initZeroPtrList +( + PtrList>& a, + const fvMesh& mesh, + const dimensionSet& ds, + const word& name +) +{ + using T = Foam::GeometricField; + + forAll(a, i) + { + if (a.set(i)) + { + a[i] == Zero; + } + else + { + a.set + ( + i, + T::New(name, mesh, dimensioned(ds, Zero)) + ); + } + } + + return; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H new file mode 100644 index 000000000..d7ad17416 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/HUEqn.H @@ -0,0 +1,84 @@ +forAll(movingPhases, movingPhasei) +{ + phaseModel& phase = movingPhases[movingPhasei]; + volVectorField& UAlpha = phase.URef(); + const volScalarField& alpha = phase; + const volScalarField& rhoAlpha = phase.thermoRef().rho(); + + HUEqns.set + ( + phase.index(), + ( + fvm::div(phase.alphaRhoPhi(), UAlpha) + + phase.divDevTau(UAlpha) + + fvm::SuSp(- fvc::div(phase.alphaRhoPhi()), UAlpha) + + MRF.DDt(alpha*rhoAlpha, UAlpha) + ) + ); +} + +// Create U & grad(U) fields +PtrList UgradUs(phases.size()); +forAll(movingPhases, movingPhasei) +{ + const phaseModel& phase = movingPhases[movingPhasei]; + const volVectorField& UAlpha = phase.U(); + const surfaceScalarField& phiAlpha = phase.phi(); + + UgradUs.set + ( + phase.index(), + new fvVectorMatrix + ( + fvm::div(phiAlpha, UAlpha) + + fvm::SuSp(- fvc::div(phiAlpha), UAlpha) + ) + ); +} + +// Add convective part of the virtual mass force +forAllConstIter(phaseSystem::VmTable, fluid.Vms(), VmIter) +{ + const volScalarField& Vm(*VmIter()); + const phasePair& pair(fluid.phasePairs()[VmIter.key()]); + + forAllConstIter(phasePair, pair, iter) + { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + HUEqns[phase.index()] += + Vm + *( + UgradUs[phase.index()] + - (UgradUs[otherPhase.index()] & otherPhase.U()) + ); + } +} + +// Relaxation factor +forAll(movingPhases, movingPhasei) +{ + const phaseModel& phase = movingPhases[movingPhasei]; + + word name = HUEqns[phase.index()].psi().select(pimple.finalPimpleIter()); + + relaxU[phase.index()] = 1; + + if (HUEqns[phase.index()].psi().mesh().relaxField(name)) + { + relaxU[phase.index()] = + HUEqns[phase.index()].psi().mesh().fieldRelaxationFactor(name); + + if + ( + pimple.finalPimpleIter() + && HUEqns[phase.index()].psi().mesh().fieldRelaxationFactor(name) != 1 + ) + { + Warning << "Field relaxation factor of " << relaxU[phase.index()] + << " for solution variable '" << name + << "' is not recommended." << endl << endl; + } + } +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H new file mode 100644 index 000000000..90a7d8cf5 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/UEqns.H @@ -0,0 +1,234 @@ +// Update mixture properties +rho = fluid.rho(); +U = fluid.U(); + +// Store velocity and flux of previous iteration +forAll(movingPhases, movingPhasei) +{ + movingPhases[movingPhasei].U()().storePrevIter(); + movingPhases[movingPhasei].phi()().storePrevIter(); +} + +// Update drag and virtual mass coefficient table and interfacial forces tables +fluid.momentumTransfer(); + +// Time derivative equation +#include "pU/ddtUEqn.H" + +// Momentum transfer by mass transfer equation +#include "pU/dmdtUEqn.H" + +// Convection, diffusion equation +#include "pU/HUEqn.H" +#include "pU/stabilizeHUEqn.H" + +// fvOptions for momentum equation +#include "pU/fvOptionsEqn.H" + +// Face volume fractions +forAll(phases, phasei) +{ + const phaseModel& phase = phases[phasei]; + const volScalarField& alpha = phase; + + alphafs.set(phasei, fvc::interpolate(alpha)); +} + +// rHUs, rHUfs, rAs, rAfs +forAll(movingPhases, movingPhasei) +{ + const phaseModel& phase = movingPhases[movingPhasei]; + + rHUs.set(phase.index(), 1.0/HUEqns[phase.index()].A()); + rHUfs.set(phase.index(), fvc::interpolate(rHUs[phase.index()])); + + // Diagonal coefficients for partial elimination + rAs.set + ( + phase.index(), + 1.0 + /( + HUEqns[phase.index()].A() + + ddtUEqns[phase.index()].A() + + fvOptionsEqns[phase.index()].A() + ) + ); + + // Face diagonal coefficients for partial elimination + rAfs.set + ( + phase.index(), + 1.0 + /( + fvc::interpolate(HUEqns[phase.index()].A()) + + fvc::interpolate(ddtUEqns[phase.index()].A()) + + fvc::interpolate(fvOptionsEqns[phase.index()].A()) + ) + ); +} +fluid.fillFields("rHUs", dimless, rHUs); +fluid.fillFields("rHUfs", dimless, rHUfs); +fluid.fillFields("rAs", dimTime/dimDensity, rAs); +fluid.fillFields("rAfs", dimTime/dimDensity, rAfs); + +// Implicit and pressure part of the drag +#include "pU/dragDiag.H" + +// Implicit part of the time derivative of the virtual mass +#include "pU/vmDiag.H" + +// Momentum predictor +if (pimple.momentumPredictor()) +{ + // Correct p_rgh for consistency with p + p_rgh = p - rho*gh; + + // Explicit part of the drag terms + #include "pU/dragSource.H" + + // Explicit part of the time derivative of the virtual mass + #include "pU/vmSource.H" + + // Add lift, wall lubrication and turbulent dispersion forces + #include "pU/forcesSourcef.H" + + // Mixture density gradient times gravitational acceleration + const surfaceScalarField ghSnGradRho + ( + "ghSnGradRho", + ghf*fvc::snGrad(rho)*mesh.magSf() + ); + + // rAUfs, alpharAUfs, phigs WITHOUT field relaxation of UAlpha + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + const volScalarField& rhoAlpha = phase.rho(); + volVectorField& UAlpha = phase.URef(); + + rAUfs.set + ( + phase.index(), + 1.0 + /( + 1.0 + + rHUfs[phase.index()] + *( + fvc::interpolate(ddtUEqns[phase.index()].A()) + + fvc::interpolate(dmdtUEqns[phase.index()].A()) + + fvc::interpolate(fvOptionsEqns[phase.index()].A()) + + dragDiagf[phase.index()] + + vmDiagf[phase.index()] + ) + ) + ); + + // Phase diagonal coefficients + alpharAUfs.set + ( + phase.index(), + rHUfs[phase.index()]*rAUfs[phase.index()] + *(alphafs[phase.index()] + dragPCoefff[phase.index()]) + ); + + // Combined buoyancy and force fluxes + phigs.set + ( + phase.index(), + alpharAUfs[phase.index()] + *( + - ghSnGradRho + + fluid.surfaceTension(phase)*mesh.magSf() + ) + ); + + // Momentum predictor equation + fvVectorMatrix UEqns + ( + ddtUEqns[phase.index()] + + dmdtUEqns[phase.index()] + + HUEqns[phase.index()] + + fvm::Sp(dragDiag[phase.index()], UAlpha) + + fvm::Sp(vmDiag[phase.index()], UAlpha) + + fvOptionsEqns[phase.index()] + == + dragSource[phase.index()] + + alpha*(rhoAlpha - rho)*g + + ( + ddtUEqns[phase.index()].A() + + dmdtUEqns[phase.index()].A() + + HUEqns[phase.index()].A() + + dragDiag[phase.index()] + + fvOptionsEqns[phase.index()].A() + ) + *fvc::reconstruct + ( + phigs[phase.index()] + - alpharAUfs[phase.index()]*fvc::snGrad(p_rgh)*mesh.magSf() + + rHUfs[phase.index()]*rAUfs[phase.index()] + *forcesSourcef[phase.index()] + ) + + vmSource[phase.index()] + ); + + UEqns.relax(); + fvOptions.constrain(UEqns); + UEqns.solve(); + fvOptions.correct(UAlpha); + } +} + +// Construct pressure "diffusivity" +surfaceScalarField rAUf +( + IOobject + ( + "rAUf", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar(dimensionSet(-1, 3, 1, 0, 0), Zero) +); + +// rAUfs, rHAUfs, alpharAUfs, rAUf WITH field relaxation of UAlpha +initZeroPtrList(rAUfs, mesh, dimless, "rAUfs"); +initZeroPtrList(rHAUfs, mesh, dimTime/dimDensity, "rHAUfs"); +initZeroPtrList(alpharAUfs, mesh, dimTime/dimDensity, "alpharAUfs"); +forAll(movingPhases, movingPhasei) +{ + const phaseModel& phase = movingPhases[movingPhasei]; + + rAUfs.set + ( + phase.index(), + 1.0 + /( + 1.0/max(relaxU[phase.index()], small) + + rHUfs[phase.index()] + *( + fvc::interpolate(ddtUEqns[phase.index()].A()) + + fvc::interpolate(dmdtUEqns[phase.index()].A()) + + fvc::interpolate(fvOptionsEqns[phase.index()].A()) + + dragDiagf[phase.index()] + + vmDiagf[phase.index()] + ) + ) + ); + + rHAUfs.set(phase.index(), rHUfs[phase.index()]*rAUfs[phase.index()]); + + // Phase diagonal coefficients + alpharAUfs.set + ( + phase.index(), + rHUfs[phase.index()]*rAUfs[phase.index()] + *(alphafs[phase.index()] + dragPCoefff[phase.index()]) + ); + + // Pressure "diffusivity" + rAUf += alphafs[phase.index()]*alpharAUfs[phase.index()]; +} + +rAUf = mag(rAUf); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H new file mode 100644 index 000000000..fb81d8e0d --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/ddtUEqn.H @@ -0,0 +1,47 @@ +forAll(movingPhases, movingPhasei) +{ + const phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + const volScalarField& rhoAlpha = phase.rho(); + const volVectorField& UAlpha = phase.U(); + + ddtUEqns.set + ( + phase.index(), + fvm::ddt(alpha, rhoAlpha, UAlpha) + - fvm::Sp(fvc::ddt(alpha, rhoAlpha), UAlpha) + ); + + word ddtScheme + ( + mesh.ddtScheme + ( + "ddt(" + + alpha.name() + + ',' + rhoAlpha.name() + ',' + + UAlpha.name() + ')' + ).last().wordToken() + ); + + if (ddtScheme == "Euler") + { + ddtPhiCoeffs.set + ( + phase.index(), + ( + byDt + ( + fvc::interpolate(alpha.oldTime()) + *fvc::interpolate(rhoAlpha.oldTime()) + *MRF.absolute(phase.phi()().oldTime()) + ) + ) + ); + } + else + { + FatalErrorInFunction + << "Only Euler is an allowed time scheme " + << "for this solver." << exit(FatalError); + } +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dmdtUEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dmdtUEqn.H new file mode 100644 index 000000000..b5b1d6bac --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dmdtUEqn.H @@ -0,0 +1,74 @@ +forAll(movingPhases, movingPhasei) +{ + const phaseModel& phase = movingPhases[movingPhasei]; + const volVectorField& UAlpha = phase.U(); + const PtrList dmdts = fluid.dmdts(); + + const dimensionedScalar ZERO(dimensionSet(1, -3, -1, 0, 0), Zero); + + dmdtUEqns.set + ( + phase.index(), + ( + fvm::Sp(ZERO, UAlpha) + + fvm::Su(Zero, UAlpha) + ) + ); + + volScalarField source + ( + volScalarField::New + ( + IOobject::groupName("source", phase.name()), + mesh, + dimensionedScalar(dimDensity/dimTime, 0) + ) + ); + + if (fvOptions.appliesToField(rho.name())) + { + source += fvOptions(phase, rho)ρ + } + + if (dmdts.set(phase.index())) + { + source += dmdts[phase.index()]; + } + + dmdtUEqns[phase.index()] += fvm::SuSp(source, UAlpha); +} + +initZeroPtrList +( + dmdtPhiCoeffs, + mesh, + dimensionSet(1, 0, -2, 0, 0), + "dmdtPhiCoeff" +); + +// Add mass sources +// Inspired by MomentumTransferPhaseSystem::addDmdtUfs() but with fvm::SuSp() +// instead of fvm::Sp() +forAllConstIter(phaseSystem::dmdtfTable, fluid.dmdtfsPhaseTransfer(), dmdtfIter) +{ + const phasePairKey& key = dmdtfIter.key(); + const phasePair& pair(fluid.phasePairs()[key]); + + const volScalarField dmdtf(Pair::compare(pair, key)**dmdtfIter()); + + phaseModel& phase1 = phases[pair.phase1().name()]; + phaseModel& phase2 = phases[pair.phase2().name()]; + + const volScalarField dmdtf21(posPart(dmdtf)); + const volScalarField dmdtf12(negPart(dmdtf)); + + dmdtUEqns[pair.phase1().index()] -= + dmdtf21*phase2.U() + fvm::SuSp(dmdtf12, phase1.URef()); + dmdtPhiCoeffs[pair.phase1().index()] += + fvc::interpolate(dmdtf21)*MRF.absolute(phase2.phi()); + + dmdtUEqns[pair.phase2().index()] += + dmdtf12*phase1.U() - fvm::SuSp(-dmdtf21, phase2.URef()); + dmdtPhiCoeffs[pair.phase2().index()] -= + fvc::interpolate(dmdtf12)*MRF.absolute(phase1.phi()); +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H new file mode 100644 index 000000000..ce183e3bf --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragDiag.H @@ -0,0 +1,71 @@ +initZeroPtrList(dragDiag, mesh, dimensionSet(1, -3, -1, 0, 0), "dragDiag"); +initZeroPtrList(dragDiagf, mesh, dimensionSet(1, -3, -1, 0, 0), "dragDiagf"); +initZeroPtrList(dragPCoefff, mesh, dimensionSet(0, 0, 0, 0, 0), "dragPCoefff"); + +forAllConstIter(phaseSystem::KdTable, fluid.Kds(), KdIter) +{ + const volScalarField& K(*KdIter()); + const phasePair& pair(fluid.phasePairs()[KdIter.key()]); + + forAllConstIter(phasePair, pair, iter) + { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + if (phase.stationary()) continue; + + if (partialElimination) + { + dragDiag[phase.index()] += + K + *( + 1 + + rAs[otherPhase.index()] + *( + HUEqns[phase.index()].A() + + ddtUEqns[phase.index()].A() + + fvOptionsEqns[phase.index()].A() + ) + ); + } + else + { + dragDiag[phase.index()] += K; + } + } +} + +forAllConstIter(phaseSystem::KdfTable, fluid.Kdfs(), KdfIter) +{ + const surfaceScalarField& Kf(*KdfIter()); + const phasePair& pair(fluid.phasePairs()[KdfIter.key()]); + + forAllConstIter(phasePair, pair, iter) + { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + if (phase.stationary()) continue; + + if (partialElimination) + { + dragDiagf[phase.index()] += + Kf + *( + 1 + + rAfs[otherPhase.index()] + *( + fvc::interpolate(HUEqns[phase.index()].A()) + + fvc::interpolate(ddtUEqns[phase.index()].A()) + + fvc::interpolate(fvOptionsEqns[phase.index()].A()) + ) + ); + + dragPCoefff[phase.index()] += Kf*rAfs[otherPhase.index()]; + } + else + { + dragDiagf[phase.index()] += Kf; + } + } +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H new file mode 100644 index 000000000..7eace2f3d --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSource.H @@ -0,0 +1,75 @@ +initZeroPtrList(dragSource, mesh, dimensionSet(1, -2, -2, 0, 0), "dragSource"); + +{ + // Summed RHS of cell-based momentum balance equation + volVectorField HSum + ( + IOobject + ( + "HSum", + runTime.timeName(), + mesh + ), + mesh, + dimensionedVector(dimensionSet(1, -2, -2, 0, 0), Zero) + ); + + if (partialElimination) + { + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + HSum += + HUEqns[phase.index()].H() + + ddtUEqns[phase.index()].H() + + fvOptionsEqns[phase.index()].H(); + } + } + + forAllConstIter(phaseSystem::KdTable, fluid.Kds(), KdIter) + { + const volScalarField& K(*KdIter()); + const phasePair& pair(fluid.phasePairs()[KdIter.key()]); + + forAllConstIter(phasePair, pair, iter) + { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + if (phase.stationary() || otherPhase.stationary()) continue; + + if (partialElimination) + { + dragSource[phase.index()] += K*rAs[otherPhase.index()]*HSum; + + forAll(movingPhases, movingPhasei) + { + const phaseModel& remainingPhase = + movingPhases[movingPhasei]; + + if + ( + (remainingPhase.index() == phase.index()) + || (remainingPhase.index() == otherPhase.index()) + ) continue; + + dragSource[phase.index()] -= + ( + K*rAs[otherPhase.index()] + *( + HUEqns[remainingPhase.index()].A() + + ddtUEqns[remainingPhase.index()].A() + + fvOptionsEqns[remainingPhase.index()].A() + ) + *phases[remainingPhase.index()].U() + ); + } + } + else + { + dragSource[phase.index()] += K*phases[otherPhase.index()].U(); + } + } + } +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H new file mode 100644 index 000000000..5f24aa1e6 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/dragSourcef.H @@ -0,0 +1,84 @@ +initZeroPtrList(dragSourcef, mesh, dimensionSet(0, 3, -1, 0, 0), "dragSourcef"); + +{ + // Summed RHS of face-based momentum balance equation + surfaceScalarField HSumf + ( + IOobject + ( + "HSumf", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar(dimensionSet(1, 0, -2, 0, 0), Zero) + ); + + if (partialElimination) + { + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + HSumf += + fvc::flux(HUEqns[phase.index()].H()) + + ddtPhiCoeffs[phase.index()] + + fvc::flux(fvOptionsEqns[phase.index()].H()); + } + } + + forAllConstIter(phaseSystem::KdfTable, fluid.Kdfs(), KdfIter) + { + const surfaceScalarField& Kf(*KdfIter()); + const phasePair& pair(fluid.phasePairs()[KdfIter.key()]); + + forAllConstIter(phasePair, pair, iter) + { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + if (phase.stationary() || otherPhase.stationary()) continue; + + if (partialElimination) + { + dragSourcef[phase.index()] + += Kf*rAfs[otherPhase.index()]*rHUfs[phase.index()]*HSumf; + + forAll(movingPhases, movingPhasei) + { + const phaseModel& remainingPhase = + movingPhases[movingPhasei]; + + if + ( + (remainingPhase.index() == phase.index()) + || (remainingPhase.index() == otherPhase.index()) + ) continue; + + dragSourcef[phase.index()] -= + ( + Kf*rAfs[otherPhase.index()]*rHUfs[phase.index()] + *( + fvc::interpolate(HUEqns[remainingPhase.index()].A()) + + fvc::interpolate + ( + ddtUEqns[remainingPhase.index()].A() + ) + + fvc::interpolate + ( + fvOptionsEqns[remainingPhase.index()].A() + ) + ) + *MRF.absolute(phases[remainingPhase.index()].phi()) + ); + } + } + else + { + dragSourcef[phase.index()] + += Kf*rHUfs[phase.index()] + *MRF.absolute(phases[otherPhase.index()].phi()); + } + } + } +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H new file mode 100644 index 000000000..2c674a319 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/forcesSourcef.H @@ -0,0 +1,36 @@ +initZeroPtrList +( + forcesSourcef, mesh, dimensionSet(1, 0, -2, 0, 0), "forcesSourcef" +); + +// Add lift force +// Inspired by MomentumTransferPhaseSystem::phiFfs() +forAllConstIter(phaseSystem::FlfTable, fluid.Flfs(), FlfIter) +{ + const surfaceScalarField& Ff(*FlfIter()); + const phasePair& pair(fluid.phasePairs()[FlfIter.key()]); + + forcesSourcef[pair.phase1().index()] -= Ff; + forcesSourcef[pair.phase2().index()] += Ff; +} + +// Add wall lubrication force +// Inspired by MomentumTransferPhaseSystem::phiFfs() +forAllConstIter(phaseSystem::FwfTable, fluid.Fwfs(), FwfIter) +{ + const surfaceScalarField& Ff(*FwfIter()); + const phasePair& pair(fluid.phasePairs()[FwfIter.key()]); + + forcesSourcef[pair.phase1().index()] -= Ff; + forcesSourcef[pair.phase2().index()] += Ff; +} + +// Add turbulent dispersion force +forAllConstIter(phaseSystem::FtfTable, fluid.Ftfs(), FtfIter) +{ + const surfaceScalarField& Ff(*FtfIter()); + const phasePair& pair(fluid.phasePairs()[FtfIter.key()]); + + forcesSourcef[pair.phase1().index()] -= Ff; + forcesSourcef[pair.phase2().index()] += Ff; +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/fvOptionsEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/fvOptionsEqn.H new file mode 100644 index 000000000..bab37091e --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/fvOptionsEqn.H @@ -0,0 +1,18 @@ +forAll(movingPhases, movingPhasei) +{ + phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + volVectorField& UAlpha = phase.URef(); + volScalarField& rhoAlpha = phase.thermoRef().rho(); + + fvOptionsEqns.set + ( + phase.index(), + fvm::Sp(dimensionedScalar(dimDensity/dimTime, Zero), UAlpha) + ); + + if(fvOptions.appliesToField(UAlpha.name())) + { + fvOptionsEqns[phase.index()] -= fvOptions(alpha, rhoAlpha, UAlpha); + } +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H new file mode 100644 index 000000000..7f2dc1e87 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/pEqn.H @@ -0,0 +1,347 @@ +rho = fluid.rho(); + +// Correct p_rgh for consistency with p and the updated densities due to +// updated temperature +p_rgh = p - rho*gh; + +// Correct fixed-flux BCs to be consistent with the velocity BCs +forAll(fluid.movingPhases(), movingPhasei) +{ + phaseModel& phase = fluid.movingPhases()[movingPhasei]; + MRF.correctBoundaryFlux(phase.U(), phase.phiRef()); +} + +{ + // Mixture density gradient times gravitational acceleration + const surfaceScalarField ghSnGradRho + ( + "ghSnGradRho", + ghf*fvc::snGrad(rho)*mesh.magSf() + ); + + // phigs WITH field relaxation of UAlpha + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + // Combined buoyancy and force fluxes + phigs.set + ( + phase.index(), + alpharAUfs[phase.index()] + *( + - ghSnGradRho + + fluid.surfaceTension(phase)*mesh.magSf() + ) + ); + } +} + +// Update explicit drag terms for PISO loop +#include "pU/dragSource.H" +#include "pU/dragSourcef.H" + +// Explicit part of the time derivative of the virtual mass for PISO loop +#include "pU/vmSource.H" +#include "pU/vmSourcef.H" + +// Update lift, wall lubrication and turbulent dispersion forces for PISO loop +#include "pU/forcesSourcef.H" + +// Construct Total predicted flux +surfaceScalarField phiHbyA +( + IOobject + ( + "phiHbyA", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar(dimArea*dimVelocity, Zero) +); + +forAll(movingPhases, movingPhasei) +{ + const phaseModel& phase = fluid.movingPhases()[movingPhasei]; + const volScalarField& rhoAlpha = phase.rho(); + const volVectorField& UAlpha = phase.U(); + + // Pseudo convected velocities + HbyAs.set + ( + phase.index(), + constrainHbyA + ( + rHUs[phase.index()]*HUEqns[phase.index()].H(), + UAlpha, + p_rgh + ) + ); + + // Predicted fluxes for each phase + phiHbyAs.set + ( + phase.index(), + constrainPhiHbyA + ( + rAUfs[phase.index()] + *( + fvc::flux(HbyAs[phase.index()]) + + rHUfs[phase.index()] + *( + ddtPhiCoeffs[phase.index()] + + dmdtPhiCoeffs[phase.index()] + + alphafs[phase.index()]*(g & mesh.Sf()) + *( + fvc::interpolate(rhoAlpha) - fvc::interpolate(rho) + ) + + forcesSourcef[phase.index()] + + fvc::flux(fvOptionsEqns[phase.index()].H()) + ) + + (1 - relaxU[phase.index()])/max(relaxU[phase.index()], small) + *MRF.absolute(phases[phase.index()].phi()().prevIter()) + + dragSourcef[phase.index()] + + vmSourcef[phase.index()] + ), + UAlpha, + p_rgh + ) + + phigs[phase.index()] + ); + + phiHbyA += alphafs[phase.index()]*phiHbyAs[phase.index()]; +} + +MRF.makeRelative(phiHbyA); + +// Update the fixedFluxPressure BCs to ensure flux consistency +{ + surfaceScalarField::Boundary phib(phi.boundaryField()); + phib = 0; + forAll(phases, phasei) + { + const phaseModel& phase = phases[phasei]; + phib += alphafs[phasei].boundaryField()*phase.phi()().boundaryField(); + } + + setSnGrad + ( + p_rgh.boundaryFieldRef(), + ( + phiHbyA.boundaryField() - phib + )/(mesh.magSf().boundaryField()*rAUf.boundaryField()) + ); +} + +// Compressible pressure equations +PtrList pEqnComps(phases.size()); +PtrList dmdts(fluid.dmdts()); +forAll(phases, phasei) +{ + phaseModel& phase = phases[phasei]; + const volScalarField& alpha = phase; + volScalarField& rhoAlpha = phase.thermoRef().rho(); + + pEqnComps.set(phasei, new fvScalarMatrix(p_rgh, dimVolume/dimTime)); + fvScalarMatrix& pEqnComp = pEqnComps[phasei]; + + // Density variation + if (!phase.isochoric() || !phase.pure()) + { + pEqnComp += + ( + fvc::ddt(alpha, rhoAlpha) + fvc::div(phase.alphaRhoPhi()) + - fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rhoAlpha) + )/rhoAlpha; + } + + // Compressibility + if (!phase.incompressible()) + { + if (pimple.transonic()) + { + const surfaceScalarField phid + ( + IOobject::groupName("phid", phase.name()), + fvc::interpolate(phase.thermo().psi())*phase.phi() + ); + + pEqnComp += + correction + ( + (alpha/rhoAlpha)* + ( + phase.thermo().psi()*fvm::ddt(p_rgh) + + fvm::div(phid, p_rgh) + - fvm::Sp(fvc::div(phid), p_rgh) + ) + ); + + deleteDemandDrivenData + ( + pEqnComps[phasei].faceFluxCorrectionPtr() + ); + + pEqnComps[phasei].relax(); + } + else + { + pEqnComp += + (alpha*phase.thermo().psi()/rhoAlpha) + *correction(fvm::ddt(p_rgh)); + } + } + + // Option sources + if (fvOptions.appliesToField(rhoAlpha.name())) + { + pEqnComp -= (fvOptions(alpha, rhoAlpha) & rhoAlpha)/rhoAlpha; + } + + // Mass transfer + if (dmdts.set(phasei)) + { + pEqnComp -= dmdts[phasei]/rhoAlpha; + } +} + +// Cache p prior to solve for density update +volScalarField p_rgh_0(p_rgh); + +// Iterate over the pressure equation to correct for non-orthogonality +while (pimple.correctNonOrthogonal()) +{ + // Construct the transport part of the pressure equation + fvScalarMatrix pEqnIncomp + ( + fvc::div(phiHbyA) + - fvm::laplacian(rAUf, p_rgh) + ); + + // Solve + { + fvScalarMatrix pEqn(pEqnIncomp); + + forAll(phases, phasei) + { + pEqn += pEqnComps[phasei]; + } + + if (fluid.incompressible()) + { + pEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); + } + + pEqn.solve(); + } + + // Correct fluxes and velocities on last non-orthogonal iteration + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqnIncomp.flux(); + + surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf); + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases[movingPhasei]; + + phase.phiRef() = + phiHbyAs[phase.index()] + + alpharAUfs[phase.index()]*mSfGradp; + + // Set the phase dilatation rate + phase.divU(-pEqnComps[phase.index()] & p_rgh); + } + + // Optionally relax pressure for velocity correction + p_rgh.relax(); + + mSfGradp = pEqnIncomp.flux()/rAUf; + + // Momentum corrector + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + const volScalarField& rhoAlpha = phase.rho(); + volVectorField& UAlpha = phase.URef(); + + if (faceMomentum) + { + UAlpha = fvc::reconstruct(phase.phi()); + } + else + { + UAlpha = + 1.0 + /( + ddtUEqns[phase.index()].A() + + dmdtUEqns[phase.index()].A() + + HUEqns[phase.index()].A() + /max(relaxU[phase.index()], small) + + dragDiag[phase.index()] + + vmDiag[phase.index()] + + fvOptionsEqns[phase.index()].A() + ) + *( + HUEqns[phase.index()].A()*HbyAs[phase.index()] + + ddtUEqns[phase.index()].H() + + dmdtUEqns[phase.index()].H() + + dragSource[phase.index()] + + alpha*(rhoAlpha - rho)*g + + (1 - relaxU[phase.index()]) + /max(relaxU[phase.index()], small) + *HUEqns[phase.index()].A()*UAlpha.prevIter() + + vmSource[phase.index()] + + fvOptionsEqns[phase.index()].H() + ) + + fvc::reconstruct + ( + phigs[phase.index()] + + alpharAUfs[phase.index()]*mSfGradp + + rHUfs[phase.index()]*rAUfs[phase.index()] + *forcesSourcef[phase.index()] + ); + } + + UAlpha.correctBoundaryConditions(); + fvOptions.correct(UAlpha); + + MRF.makeRelative(phase.phiRef()); + } + } +} + +// Update the mixture properties and static pressure +U = fluid.U(); +p = p_rgh + rho*gh; + +// Account for static pressure reference +if (p_rgh.needReference() && fluid.incompressible()) +{ + p += dimensionedScalar + ( + "p", + p.dimensions(), + pRefValue - getRefCellValue(p, pRefCell) + ); +} + +// Limit p_rgh +p_rgh = p - rho*gh; + +// Update densities from change in p_rgh +forAll(phases, phasei) +{ + phaseModel& phase = phases[phasei]; + phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0); +} + +// Correct p_rgh for consistency with p and the updated densities +rho = fluid.rho(); +p_rgh = p - rho*gh; +p_rgh.correctBoundaryConditions(); diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H new file mode 100644 index 000000000..22fb7747b --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/stabilizeHUEqn.H @@ -0,0 +1,43 @@ +forAll(movingPhases, movingPhasei) +{ + phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + const volScalarField& rhoAlpha = phase.rho(); + volVectorField& UAlpha = phase.URef(); + + // Stabilize due to zero flux + const dimensionedScalar residualDivU(dimensionSet(1, -3, -1, 0, 0), vSmall); + + HUEqns[phase.index()] += + fvm::Sp + ( + max + ( + residualDivU - fvc::div(phase.alphaRhoPhi()), + dimensionedScalar(dimensionSet(1, -3, -1, 0, 0), Zero) + ), + UAlpha + ) + - fvc::Sp + ( + max + ( + residualDivU - fvc::div(phase.alphaRhoPhi()), + dimensionedScalar(dimensionSet(1, -3, -1, 0, 0), Zero) + ), + UAlpha + ); + + // Stabilize due to zero phase fraction + HUEqns[phase.index()] += + fvm::ddt + ( + max(phase.residualAlpha() - alpha, scalar(0))*rhoAlpha, + UAlpha + ) + - fvc::ddt + ( + max(phase.residualAlpha() - alpha, scalar(0))*rhoAlpha, + UAlpha + ); +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H new file mode 100644 index 000000000..1a6f80659 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmDiag.H @@ -0,0 +1,16 @@ +initZeroPtrList(vmDiag, mesh, dimensionSet(1, -3, -1, 0, 0), "vmDiag"); +initZeroPtrList(vmDiagf, mesh, dimensionSet(1, -3, -1, 0, 0), "vmDiagf"); + +forAllConstIter(phaseSystem::VmTable, fluid.Vms(), VmIter) +{ + const volScalarField& Vm(*VmIter()); + const phasePair& pair(fluid.phasePairs()[VmIter.key()]); + + forAllConstIter(phasePair, pair, iter) + { + const phaseModel& phase = iter(); + + vmDiag[phase.index()] += byDt(Vm); + vmDiagf[phase.index()] += byDt(fvc::interpolate(Vm)); + } +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H new file mode 100644 index 000000000..ddaabab31 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSource.H @@ -0,0 +1,22 @@ +initZeroPtrList(vmSource, mesh, dimensionSet(1, -2, -2, 0, 0), "vmSource"); + +forAllConstIter(phaseSystem::VmTable, fluid.Vms(), VmIter) +{ + const volScalarField& Vm(*VmIter()); + const phasePair& pair(fluid.phasePairs()[VmIter.key()]); + + forAllConstIter(phasePair, pair, iter) + { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + vmSource[phase.index()] += + Vm + *byDt + ( + phase.U()().oldTime() + + otherPhase.U() + - otherPhase.U()().oldTime() + ); + } +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H new file mode 100644 index 000000000..eed96da3b --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/cipsaMultiphaseEulerFoam/pU/vmSourcef.H @@ -0,0 +1,23 @@ +initZeroPtrList(vmSourcef, mesh, dimensionSet(0, 3, -1, 0, 0), "vmSourcef"); + +forAllConstIter(phaseSystem::VmTable, fluid.Vms(), VmIter) +{ + const volScalarField& Vm(*VmIter()); + const phasePair& pair(fluid.phasePairs()[VmIter.key()]); + + forAllConstIter(phasePair, pair, iter) + { + const phaseModel& phase = iter(); + const phaseModel& otherPhase = iter.otherPhase(); + + vmSourcef[phase.index()] += + rHUfs[phase.index()] + *fvc::interpolate(Vm) + *byDt + ( + MRF.absolute(phase.phi()().oldTime()) + + MRF.absolute(otherPhase.phi()) + - MRF.absolute(otherPhase.phi()().oldTime()) + ); + } +} diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/Make/files b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/Make/files index 24000605a..8850485af 100644 --- a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/Make/files +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/Make/files @@ -22,4 +22,9 @@ saturationModels/polynomial/polynomial.C saturationModels/function1/function1.C saturationModels/constantSaturationConditions/constantSaturationConditions.C +SGSdeltas/SGSdelta/SGSdelta.C +SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.C +SGSdeltas/maxDeltaxyz/maxDeltaxyz.C +SGSdeltas/oneDimensionalDelta/oneDimensionalDelta.C + LIB = $(FOAM_LIBBIN)/libeulerianInterfacialCompositionModels diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.C new file mode 100644 index 000000000..fedaf01a2 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.C @@ -0,0 +1,93 @@ +/*---------------------------------------------------------------------------*\ + == == ====== ==== ==== | + \\ || | Unsupported Contributions for OpenFOAM + ====== // || || ===// | + || || // || // || \\ | (C) 2019 Helmholtz-Zentrum Dresden-Rossendorf + == == ====== ==== == == | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "SGSdelta.H" +#include "calculatedFvPatchFields.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(SGSdelta, 0); + defineRunTimeSelectionTable(SGSdelta, dictionary); +} + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::SGSdelta::SGSdelta +( + const word& name, + const fvMesh& mesh +) +: + mesh_(mesh), + delta_ + ( + IOobject + ( + name, + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar(name, dimLength, small), + calculatedFvPatchScalarField::typeName + ) +{} + + +// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * // + +Foam::autoPtr Foam::SGSdelta::New +( + const word& name, + const fvMesh& mesh, + const dictionary& dict +) +{ + const word deltaType(dict.lookup("delta")); + + Info<< "Selecting SGS delta type " << deltaType << endl; + + dictionaryConstructorTable::iterator cstrIter = + dictionaryConstructorTablePtr_->find(deltaType); + + if (cstrIter == dictionaryConstructorTablePtr_->end()) + { + FatalErrorInFunction + << "Unknown SGSdelta type " + << deltaType << nl << nl + << "Valid SGSdelta types are :" << endl + << dictionaryConstructorTablePtr_->sortedToc() + << exit(FatalError); + } + + return autoPtr(cstrIter()(name, mesh, dict)); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.H b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.H new file mode 100644 index 000000000..d8faefa03 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/SGSdelta/SGSdelta.H @@ -0,0 +1,143 @@ +/*---------------------------------------------------------------------------*\ + == == ====== ==== ==== | + \\ || | Unsupported Contributions for OpenFOAM + ====== // || || ===// | + || || // || // || \\ | (C) 2019 Helmholtz-Zentrum Dresden-Rossendorf + == == ====== ==== == == | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::SGSdelta + +Description + Abstract base class for SGS deltas + +See also + Foam::LESdelta + +SourceFiles + SGSdelta.C + +\*---------------------------------------------------------------------------*/ + +#ifndef SGSdelta_H +#define SGSdelta_H + +#include "volFields.H" +#include "runTimeSelectionTables.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class SGSdelta Declaration +\*---------------------------------------------------------------------------*/ + +class SGSdelta +{ +protected: + + // Protected data + + //- Reference to the mesh + const fvMesh& mesh_; + + volScalarField delta_; + + +public: + + //- Runtime type information + TypeName("SGSdelta"); + + + // Declare run-time constructor selection table + + declareRunTimeSelectionTable + ( + autoPtr, + SGSdelta, + dictionary, + ( + const word& name, + const fvMesh& mesh, + const dictionary& dict + ), + (name, mesh, dict) + ); + + + // Constructors + + //- Construct from name and mesh + SGSdelta + ( + const word& name, + const fvMesh& mesh + ); + + //- Disallow default bitwise copy construction + SGSdelta(const SGSdelta&) = delete; + + // Selectors + + //- Return a reference to the selected SGS delta + static autoPtr New + ( + const word& name, + const fvMesh& mesh, + const dictionary& dict + ); + + + //- Destructor + virtual ~SGSdelta() + {} + + + // Member Functions + + //- Read the SGSdelta dictionary + virtual void read(const dictionary&) = 0; + + // Correct values + virtual void correct() = 0; + + + // Member Operators + + void operator=(const SGSdelta&) = delete; + + virtual operator const volScalarField&() const + { + return delta_; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.C b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.C new file mode 100644 index 000000000..c5ed71b21 --- /dev/null +++ b/applications/solvers/multiphase/multiphaseEulerFoam/interfacialCompositionModels/SGSdeltas/cubeRootVolDelta/cubeRootVolDelta.C @@ -0,0 +1,126 @@ +/*---------------------------------------------------------------------------*\ + == == ====== ==== ==== | + \\ || | Unsupported Contributions for OpenFOAM + ====== // || || ===// | + || || // || // || \\ | (C) 2019 Helmholtz-Zentrum Dresden-Rossendorf + == == ====== ==== == == | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "cubeRootVolDelta.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace SGSModels +{ + defineTypeNameAndDebug(cubeRootVolDelta, 0); + addToRunTimeSelectionTable(SGSdelta, cubeRootVolDelta, dictionary); +} +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::SGSModels::cubeRootVolDelta::calcDelta() +{ + label nD = mesh_.nGeometricD(); + + if (nD == 3) + { + delta_.primitiveFieldRef() = deltaCoeff_*pow(mesh_.V(), 1.0/3.0); + } + else if (nD == 2) + { + WarningInFunction + << "Case is 2D, SGS is not strictly applicable\n" + << endl; + + const Vector