diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 5247c979..425ae77a 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -33,7 +33,7 @@ jobs: os: ['ubuntu-latest'] defaults: run: - shell: bash -l {0} + shell: bash -leo pipefail {0} working-directory: ${{github.workspace}} steps: - uses: actions/checkout@v4 @@ -67,7 +67,7 @@ jobs: os: ['ubuntu-latest', 'macos-latest'] defaults: run: - shell: bash -l {0} + shell: bash -leo pipefail {0} working-directory: ${{github.workspace}} steps: @@ -106,7 +106,7 @@ jobs: os: ['ubuntu-latest', 'macos-latest'] defaults: run: - shell: bash -l {0} + shell: bash -leo pipefail {0} working-directory: ${{github.workspace}} steps: @@ -119,7 +119,7 @@ jobs: cache: true cache-write: ${{ github.event_name == 'push' && github.ref == 'refs/heads/main' }} - - name: Formatting and sorting import + - name: Test run: | pixi run test @@ -133,7 +133,7 @@ jobs: os: ['ubuntu-latest'] defaults: run: - shell: bash -l {0} + shell: bash -leo pipefail {0} working-directory: ${{github.workspace}} steps: - uses: actions/checkout@v4 @@ -162,8 +162,8 @@ jobs: - name: Test run: nox -s tests -- no-reports - Regression-Test-OF: - name: Regression-Test-OF (${{ matrix.python-version }}, ${{ matrix.os }}) + Regression-Test-OF9: + name: Regression-Test-OF9 (${{ matrix.python-version }}, ${{ matrix.os }}) runs-on: ${{ matrix.os }} strategy: fail-fast: false @@ -172,7 +172,7 @@ jobs: os: ['ubuntu-22.04'] defaults: run: - shell: bash -l {0} + shell: bash -leo pipefail {0} working-directory: ${{github.workspace}} steps: - uses: actions/checkout@v4 @@ -202,64 +202,109 @@ jobs: - name: Compile solver run: | - cd OFsolvers/birdmultiphaseEulerFoam + cd OFsolvers/OF9/birdmultiphaseEulerFoam export WM_COMPILE_OPTION=Debug ./Allwmake - cd ../../ - name: Run deckwer17 PBE run: | - cd experimental_cases/deckwer17 + cd experimental_cases/OF9/deckwer17 bash run.sh - cd ../../ - name: Run deckwer17 constantD run: | - cd experimental_cases/deckwer17 + cd experimental_cases/OF9/deckwer17 cp constant/phaseProperties_constantd constant/phaseProperties bash run.sh - cd ../../ - name: Run deckwer19 PBE run: | - cd experimental_cases/deckwer19 + cd experimental_cases/OF9/deckwer19 bash run.sh - cd ../../ - name: Run side sparger tutorial run: | - cd tutorial_cases/side_sparger + cd tutorial_cases/OF9/side_sparger bash run.sh - cd ../../ - name: Run bubble column tutorial run: | - cd tutorial_cases/bubble_column_20L + cd tutorial_cases/OF9/bubble_column_20L bash run.sh - cd ../../ - name: Run stirred-tank tutorial run: | - cd tutorial_cases/stirred_tank + cd tutorial_cases/OF9/stirred_tank bash run.sh - cd ../../ - name: Run reactive loop reactor tutorial run: | - cd tutorial_cases/loop_reactor_reacting + cd tutorial_cases/OF9/loop_reactor_reacting bash run.sh - cd ../../ - name: Run mixing loop reactor tutorial run: | - cd tutorial_cases/loop_reactor_mixing + cd tutorial_cases/OF9/loop_reactor_mixing bash run.sh - cd ../../ - name: Run airlift reactor tutorial run: | - cd tutorial_cases/airlift_40m + cd tutorial_cases/OF9/airlift_40m bash run.sh - cd ../../ - name: Run flat panel reactor tutorial run: | - cd tutorial_cases/FlatPanel_250L_ASU + cd tutorial_cases/OF9/FlatPanel_250L_ASU bash run.sh - cd ../../ + + Regression-Test-OF13: + name: Regression-Test-OF13 (${{ matrix.python-version }}, ${{ matrix.os }}) + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + python-version: ['3.14'] + os: ['ubuntu-22.04'] + defaults: + run: + shell: bash -leo pipefail {0} + working-directory: ${{github.workspace}} + steps: + - uses: actions/checkout@v4 + - uses: actions/setup-python@v5 + with: + python-version: ${{matrix.python-version}} + - uses: gerlero/setup-openfoam@v1 + with: + openfoam-version: 13 + + - name: Set up micromamba + uses: mamba-org/setup-micromamba@v1 + with: + environment-name: test-env + create-args: >- + python=${{ matrix.python-version }} + channels: conda-forge + channel-priority: strict + cache-downloads: true + cache-env: true + + - name: Install dependencies + run: | + micromamba install --yes -n test-env -c conda-forge paraview + pip install --upgrade pip + pip install . + + - name: Compile solver + run: | + cd OFsolvers/OF13/birdmultiphaseEuler + export WM_COMPILE_OPTION=Debug + ./Allwmake + - name: Run bubble column + run: | + cd tutorial_cases/OF13/bubble_column + ./Allrun --test + - name: STR MRF + run: | + cd tutorial_cases/OF13/STR_MRF + ./Allrun --test + - name: STR rotating mesh + run: | + cd tutorial_cases/OF13/STR_rotatingMesh + ./Allrun --test - Regression-Test-OF-pixi: - name: Regression-Test-OF-pixi (${{ matrix.os }}) + Regression-Test-OF9-pixi: + name: Regression-Test-OF9-pixi (${{ matrix.os }}) runs-on: ${{ matrix.os }} strategy: fail-fast: false @@ -267,7 +312,7 @@ jobs: os: ['ubuntu-22.04'] defaults: run: - shell: bash -l {0} + shell: bash -leo pipefail {0} working-directory: ${{github.workspace}} steps: - uses: actions/checkout@v4 @@ -282,13 +327,11 @@ jobs: cache-write: ${{ github.event_name == 'push' && github.ref == 'refs/heads/main' }} - name: Compile solver run: | - cd OFsolvers/birdmultiphaseEulerFoam + cd OFsolvers/OF9/birdmultiphaseEulerFoam export WM_COMPILE_OPTION=Debug ./Allwmake - cd ../../ - name: Run bubble column tutorial run: | - cd tutorial_cases/bubble_column_20L + cd tutorial_cases/OF9/bubble_column_20L bash run_pixi.sh - cd ../../ diff --git a/OFsolvers/OF13/birdmultiphaseEuler/Allwclean b/OFsolvers/OF13/birdmultiphaseEuler/Allwclean new file mode 100755 index 00000000..f9e3cd1f --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/Allwclean @@ -0,0 +1,14 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory + +wclean libso phaseSystem +wclean libso populationBalance +wclean libso thermophysicalTransportModels +wclean libso reactions +wclean libso momentumTransportModels +wclean + +wclean libso fvModels +wclean libso functionObjects + +#------------------------------------------------------------------------------ diff --git a/OFsolvers/OF13/birdmultiphaseEuler/Allwmake b/OFsolvers/OF13/birdmultiphaseEuler/Allwmake new file mode 100755 index 00000000..fc432ad1 --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/Allwmake @@ -0,0 +1,17 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory + +# Parse arguments for library compilation +. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments + +wmake -j 20 $targetType phaseSystem +wmake -j 20 $targetType populationBalance +wmake -j 20 $targetType thermophysicalTransportModels +wmake -j 20 $targetType reactions +wmake -j 20 $targetType momentumTransportModels +wmake $targetType + +wmake -j 20 $targetType fvModels +wmake -j 20 $targetType functionObjects + +#------------------------------------------------------------------------------ diff --git a/OFsolvers/OF13/birdmultiphaseEuler/Make/files b/OFsolvers/OF13/birdmultiphaseEuler/Make/files new file mode 100644 index 00000000..5f5c35c6 --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/Make/files @@ -0,0 +1,11 @@ +setRDeltaT.C +moveMesh.C +momentumPredictor.C +thermophysicalPredictor.C +compressibilityEqns.C +facePressureCorrector.C +cellPressureCorrector.C +pressureCorrector.C +birdmultiphaseEuler.C + +LIB = $(FOAM_USER_LIBBIN)/libbirdmultiphaseEuler diff --git a/OFsolvers/OF13/birdmultiphaseEuler/Make/options b/OFsolvers/OF13/birdmultiphaseEuler/Make/options new file mode 100644 index 00000000..2d25e36f --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/Make/options @@ -0,0 +1,24 @@ +EXE_INC = \ + -I$(FOAM_MODULES)/fluidSolver/lnInclude \ + -IphaseSystem/lnInclude \ + -IpopulationBalance/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/physicalProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \ + -I$(LIB_SRC)/twoPhaseModels/compressibleTwoPhases/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude + +LIB_LIBS = \ + -L$(FOAM_USER_LIBBIN) \ + -lfluidSolver \ + -lbirdpopulationBalance \ + -lbirdmultiphaseEulerThermophysicalTransportModels \ + -lbirdmultiphaseEulerMomentumTransportModels \ + -lbirdmultiphaseEulerReactions \ + -lfvModels \ + -lfvConstraints diff --git a/OFsolvers/OF13/birdmultiphaseEuler/birdmultiphaseEuler.C b/OFsolvers/OF13/birdmultiphaseEuler/birdmultiphaseEuler.C new file mode 100644 index 00000000..6bbc419f --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/birdmultiphaseEuler.C @@ -0,0 +1,319 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2022-2025 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 . + +\*---------------------------------------------------------------------------*/ + +#include "birdmultiphaseEuler.H" +#include "localEulerDdtScheme.H" +#include "surfaceFields.H" +#include "fvcDiv.H" +#include "fvcSurfaceIntegrate.H" +#include "fvcMeshPhi.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solvers +{ + defineTypeNameAndDebug(birdmultiphaseEuler, 0); + addToRunTimeSelectionTable(solver, birdmultiphaseEuler, fvMesh); +} +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +bool Foam::solvers::birdmultiphaseEuler::read() +{ + fluidSolver::read(); + + predictMomentum = + pimple.dict().lookupOrDefault("momentumPredictor", false); + + faceMomentum = + pimple.dict().lookupOrDefault("faceMomentum", false); + + dragCorrection = + pimple.dict().lookupOrDefault("dragCorrection", false); + + nEnergyCorrectors = + pimple.dict().lookupOrDefault("nEnergyCorrectors", 1); + + nTurbulenceCorrectors = + pimple.dict().lookupOrDefault("nTurbulenceCorrectors", 1); + + alphaControls.read(mesh.solution().solverDict("alpha")); + + return true; +} + + +void Foam::solvers::birdmultiphaseEuler::correctCoNum() +{ + scalarField sumPhi + ( + fvc::surfaceSum(mag(phi))().primitiveField() + ); + + forAll(movingPhases, movingPhasei) + { + sumPhi = max + ( + sumPhi, + fvc::surfaceSum(mag(movingPhases[movingPhasei].phi()))() + .primitiveField() + ); + } + + CoNum_ = 0.5*gMax(sumPhi/mesh.V().primitiveField())*runTime.deltaTValue(); + + const scalar meanCoNum = + 0.5 + *(gSum(sumPhi)/gSum(mesh.V().primitiveField())) + *runTime.deltaTValue(); + + Info<< "Courant Number mean: " << meanCoNum + << " max: " << CoNum << endl; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::solvers::birdmultiphaseEuler::birdmultiphaseEuler(fvMesh& mesh) +: + fluidSolver(mesh), + + predictMomentum + ( + pimple.dict().lookupOrDefault("momentumPredictor", false) + ), + + faceMomentum + ( + pimple.dict().lookupOrDefault("faceMomentum", false) + ), + + dragCorrection + ( + pimple.dict().lookupOrDefault("dragCorrection", false) + ), + + interphaseHeatCorrectionScheme + ( + pimple.dict().lookupOrDefault("interphaseHeatCorrectionScheme", false) + ), + + nEnergyCorrectors + ( + pimple.dict().lookupOrDefault("nEnergyCorrectors", 1) + ), + + nTurbulenceCorrectors + ( + pimple.dict().lookupOrDefault("nTurbulenceCorrectors", 1) + ), + + trDeltaT + ( + LTS + ? new volScalarField + ( + IOobject + ( + fv::localEulerDdt::rDeltaTName, + runTime.name(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar(dimless/dimTime, 1), + extrapolatedCalculatedFvPatchScalarField::typeName + ) + : nullptr + ), + + trDeltaTf + ( + LTS && faceMomentum + ? new surfaceScalarField + ( + IOobject + ( + fv::localEulerDdt::rDeltaTfName, + runTime.name(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar(dimless/dimTime, 1) + ) + : nullptr + ), + + buoyancy(mesh), + + fluid_(mesh), + + phases_(fluid_.phases()), + + movingPhases_(fluid_.movingPhases()), + + phi_(fluid_.phi()), + + momentumTransferSystem_(fluid_), + + heatTransferSystem_(fluid_), + + populationBalanceSystem_(fluid_), + + p_(movingPhases_[0].fluidThermo().p()), + + p_rgh_(buoyancy.p_rgh), + + pressureReference + ( + p_, + p_rgh_, + pimple.dict(), + fluid_.incompressible() + ), + + MRF(fluid_.MRF()), + + fluid(fluid_), + phases(phases_), + movingPhases(movingPhases_), + momentumTransfer(momentumTransferSystem_), + heatTransfer(heatTransferSystem_), + p(p_), + p_rgh(p_rgh_), + phi(phi_), + rhoCpvs(fluid.thermalPhases().size()) +{ + // Read the controls + read(); + + mesh.schemes().setFluxRequired(p_rgh.name()); + + if (transient()) + { + correctCoNum(); + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::solvers::birdmultiphaseEuler::~birdmultiphaseEuler() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::solvers::birdmultiphaseEuler::preSolve() +{ + if (transient()) + { + correctCoNum(); + } + else if (LTS) + { + setRDeltaT(); + } + + // Store divU from the previous mesh so that it can be + // mapped and used in correctPhi to ensure the corrected phi + // has the same divergence + if (correctPhi || mesh.topoChanging()) + { + // Construct and register divU for mapping + divU = new volScalarField + ( + "divU0", + fvc::div(fvc::absolute(phi, movingPhases[0].U())) + ); + } + + fvModels().preUpdateMesh(); + + // Update the mesh for topology change, mesh to mesh mapping + mesh_.update(); +} + + +void Foam::solvers::birdmultiphaseEuler::prePredictor() +{ + if (pimple.thermophysics() || pimple.flow()) + { + alphaControls.correct(CoNum); + + fluid_.solve(alphaControls, rAs, momentumTransferSystem_); + populationBalanceSystem_.solve(); + + fluid_.correct(); + populationBalanceSystem_.correct(); + + fluid_.correctContinuityError(populationBalanceSystem_.dmdts()); + } +} + + +void Foam::solvers::birdmultiphaseEuler::momentumTransportPredictor() +{ + fluid_.predictMomentumTransport(); +} + + +void Foam::solvers::birdmultiphaseEuler::thermophysicalTransportPredictor() +{ + // Moved inside the nEnergyCorrectors loop in thermophysicalPredictor() + // fluid_.predictThermophysicalTransport(); +} + + +void Foam::solvers::birdmultiphaseEuler::momentumTransportCorrector() +{ + for (int i = 0; i < nTurbulenceCorrectors; i++) + { + fluid_.correctMomentumTransport(); + } +} + + +void Foam::solvers::birdmultiphaseEuler::thermophysicalTransportCorrector() +{ + fluid_.correctThermophysicalTransport(); +} + + +void Foam::solvers::birdmultiphaseEuler::postSolve() +{ + divU.clear(); +} + + +// ************************************************************************* // diff --git a/OFsolvers/OF13/birdmultiphaseEuler/birdmultiphaseEuler.H b/OFsolvers/OF13/birdmultiphaseEuler/birdmultiphaseEuler.H new file mode 100644 index 00000000..1c06c3c2 --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/birdmultiphaseEuler.H @@ -0,0 +1,325 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2022-2025 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 . + +Class + Foam::solvers::birdmultiphaseEuler + +Description + Solver module for a system of any number of compressible fluid phases with a + common pressure, but otherwise separate properties. The type of phase model + is run time selectable and can optionally represent multiple species and + in-phase reactions. The phase system is also run time selectable and can + optionally represent different types of momentum, heat and mass transfer. + + Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and + pseudo-transient and steady simulations. + + Optional fvModels and fvConstraints are provided to enhance the simulation + in many ways including adding various sources, Lagrangian + particles, surface film etc. and constraining or limiting the solution. + +SourceFiles + birdmultiphaseEuler.C + +See also + Foam::solvers::compressibleVoF + Foam::solvers::fluidSolver + Foam::solvers::incompressibleFluid + +\*---------------------------------------------------------------------------*/ + +#ifndef birdmultiphaseEuler_H +#define birdmultiphaseEuler_H + +#include "fluidSolver.H" +#include "momentumTransferSystem.H" +#include "heatTransferSystem.H" +#include "populationBalanceSystem.H" +#include "phaseCompressibleMomentumTransportModel.H" +#include "buoyancy.H" +#include "pressureReference.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace solvers +{ + +/*---------------------------------------------------------------------------*\ + Class birdmultiphaseEuler Declaration +\*---------------------------------------------------------------------------*/ + +class birdmultiphaseEuler +: + public fluidSolver +{ + +protected: + + // Controls + + //- Momentum equation predictor switch + // Defaults to false + Switch predictMomentum; + + //- Cell/face momentum equation switch + // 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; + + Switch interphaseHeatCorrectionScheme; + + //- Number of energy correctors + // Used to improve stability of phase-change simulations + // Defaults to 1 + int nEnergyCorrectors; + + // Number of turbulence correctors + int nTurbulenceCorrectors; + + phaseSystem::alphaControl alphaControls; + + + //- Optional LTS reciprocal time-step field + tmp trDeltaT; + + //- Optional LTS reciprocal face time-step field + tmp trDeltaTf; + + //- Buoyancy force + solvers::buoyancy buoyancy; + + + // Phase properties + + phaseSystem fluid_; + + phaseSystem::phaseModelList& phases_; + + phaseSystem::phaseModelPartialList& movingPhases_; + + surfaceScalarField& phi_; + + momentumTransferSystem momentumTransferSystem_; + + heatTransferSystem heatTransferSystem_; + + populationBalanceSystem populationBalanceSystem_; + + + // Thermophysical properties + + volScalarField& p_; + + //- Reference to the buoyant pressure for buoyant cases + // otherwise to the pressure + volScalarField& p_rgh_; + + + // Pressure reference + + //- Pressure reference + Foam::pressureReference pressureReference; + + + // Optional models + + const IOMRFZoneList& MRF; + + + // Cached temporary fields + + //- Temporary phase momentum matrices + // shared between the momentum predictor and pressure corrector + PtrList UEqns; + + //- Temporary storage for the reciprocal momentum equation diagonal + // Used by the phase-fraction predictor and pressure corrector + PtrList rAs; + + //- Stored divU from the previous mesh so that it can be + // mapped and used in correctPhi to ensure the corrected phi + // has the same divergence + autoPtr divU; + + //- Read controls + virtual bool read(); + + +private: + + // Private Member Functions + + //- Set rDeltaT for LTS + virtual void setRDeltaT(); + + //- Correct the cached Courant numbers + void correctCoNum(); + + //- Composition predictor + // called in prePredictor() after phase-fraction prediction + void compositionPredictor(); + + //- Energy and temperature predictor + void energyPredictor(); + + //- Cell-based momentum equation predictor + void cellMomentumPredictor(); + + //- Face-based momentum equation predictor + void faceMomentumPredictor(); + + //- Cell-based pressure equation predictor + void cellPressureCorrector(); + + //- Face-based pressure equation predictor + void facePressureCorrector(); + + //- Return the list of pressure equation compressibility contributions + PtrList compressibilityEqns + ( + const PtrList& dmdts + ) const; + + bool solveComposition() + { + return pimple.dict().lookupOrDefault("solveComposition",true); + } + + bool solveEnergy() + { + return pimple.dict().lookupOrDefault("solveEnergy",true); + } + + +public: + + // Public Data + + //- Reference to the multiphase fluid + const phaseSystem& fluid; + + //- Reference to the phases + const phaseSystem::phaseModelList& phases; + + //- Reference to the moving phases + const phaseSystem::phaseModelPartialList& movingPhases; + + //- Reference to the momentum transfer system + const momentumTransferSystem& momentumTransfer; + + //- Reference to the heat transfer system + const heatTransferSystem& heatTransfer; + + //- Reference to the pressure field + const volScalarField& p; + + //- Reference to the buoyant pressure for buoyant cases + // otherwise to the pressure + const volScalarField& p_rgh; + + //- Reference to the mass-flux field + const surfaceScalarField& phi; + + PtrList rhoCpvs; + + + //- Runtime type information + TypeName("birdmultiphaseEuler"); + + + // Constructors + + //- Construct from region mesh + birdmultiphaseEuler(fvMesh& mesh); + + //- Disallow default bitwise copy construction + birdmultiphaseEuler(const birdmultiphaseEuler&) = delete; + + + //- Destructor + virtual ~birdmultiphaseEuler(); + + + // Member Functions + + //- Called at the start of the time-step, before the PIMPLE loop + virtual void preSolve(); + + //- Called at the start of the PIMPLE loop to move the mesh + virtual void moveMesh(); + + //- Corrections that follow mesh motion + virtual void motionCorrector(); + + //- Called at the start of the PIMPLE loop + virtual void prePredictor(); + + //- Predict the momentum transport + virtual void momentumTransportPredictor(); + + //- Predict thermophysical transport + virtual void thermophysicalTransportPredictor(); + + //- Construct and optionally solve the momentum equation + virtual void momentumPredictor(); + + //- Construct and solve the energy equation, + // convert to temperature + // and update thermophysical and transport properties + virtual void thermophysicalPredictor(); + + //- Construct and solve the pressure equation in the PISO loop + virtual void pressureCorrector(); + + //- Correct the momentum transport + virtual void momentumTransportCorrector(); + + //- Correct the thermophysical transport + virtual void thermophysicalTransportCorrector(); + + //- Called after the PIMPLE loop at the end of the time-step + virtual void postSolve(); + + + // Member Operators + + //- Disallow default bitwise assignment + void operator=(const birdmultiphaseEuler&) = delete; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace solvers +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/OFsolvers/OF13/birdmultiphaseEuler/cellPressureCorrector.C b/OFsolvers/OF13/birdmultiphaseEuler/cellPressureCorrector.C new file mode 100644 index 00000000..5faf715a --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/cellPressureCorrector.C @@ -0,0 +1,504 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2022-2025 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 . + +\*---------------------------------------------------------------------------*/ + +#include "birdmultiphaseEuler.H" +#include "constrainHbyA.H" +#include "constrainPressure.H" +#include "findRefCell.H" +#include "fvcDdt.H" +#include "fvcDiv.H" +#include "fvcSup.H" +#include "fvcSnGrad.H" +#include "fvmDdt.H" +#include "fvmDiv.H" +#include "fvmLaplacian.H" +#include "fvmSup.H" +#include "fvcFlux.H" +#include "fvcMeshPhi.H" +#include "fvcReconstruct.H" + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +void Foam::solvers::birdmultiphaseEuler::cellPressureCorrector() +{ + volScalarField& p(p_); + volScalarField& p_rgh = p_rgh_; + + volScalarField rho("rho", fluid.rho()); + + // Correct p_rgh for consistency with the current density + p_rgh = p - rho*buoyancy.gh - buoyancy.pRef; + + // Face volume fractions + PtrList alphafs(phases.size()); + forAll(phases, phasei) + { + const phaseModel& phase = phases[phasei]; + const volScalarField& alpha = phase; + + alphafs.set(phasei, fvc::interpolate(max(alpha, scalar(0))).ptr()); + alphafs[phasei].rename("pEqn" + alphafs[phasei].name()); + } + + // Diagonal coefficients + rAs.clear(); + if (fluid.implicitPhasePressure()) + { + rAs.setSize(phases.size()); + } + + PtrList HVms(movingPhases.size()); + PtrList> invADVs; + PtrList> invADVfs; + { + PtrList As(movingPhases.size()); + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + + As.set + ( + movingPhasei, + UEqns[phase.index()].A() + + byDt + ( + max(phase.residualAlpha() - alpha, scalar(0)) + *phase.rho() + ) + ); + + if (fluid.implicitPhasePressure()) + { + rAs.set + ( + phase.index(), + new volScalarField + ( + IOobject::groupName("rA", phase.name()), + 1/As[movingPhasei] + ) + ); + } + } + + momentumTransferSystem_.invADVs(As, HVms, invADVs, invADVfs); + } + + // Explicit force fluxes + PtrList alphaByADfs; + PtrList FgByADfs; + { + PtrList Ffs(momentumTransferSystem_.Fs()); + + const surfaceScalarField ghSnGradRho + ( + "ghSnGradRho", + buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf() + ); + + UPtrList movingAlphafs(movingPhases.size()); + PtrList Fgfs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + movingAlphafs.set(movingPhasei, &alphafs[phase.index()]); + + Fgfs.set + ( + movingPhasei, + Ffs[phase.index()] + + alphafs[phase.index()] + *( + ghSnGradRho + - fluid.surfaceTension(phase)*mesh.magSf() + ) + - fvc::interpolate(max(phase, phase.residualAlpha())) + *fvc::interpolate(phase.rho() - rho)*(buoyancy.g & mesh.Sf()) + ); + } + + alphaByADfs = invADVfs & movingAlphafs; + FgByADfs = invADVfs & Fgfs; + } + + // Mass transfer rates + PtrList dmdts(populationBalanceSystem_.dmdts()); + + // --- Optional momentum predictor + if (predictMomentum) + { + PtrList HbyADs; + { + PtrList Hs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + + Hs.set + ( + movingPhasei, + UEqns[phase.index()].H() + + byDt + ( + max(phase.residualAlpha() - alpha, scalar(0)) + *phase.rho() + ) + *phase.U()().oldTime() + ); + + if (HVms.set(movingPhasei)) + { + Hs[movingPhasei] += HVms[movingPhasei]; + } + } + + HbyADs = invADVs & Hs; + } + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + constrainHbyA(HbyADs[movingPhasei], phase.U(), p_rgh); + } + + const surfaceScalarField mSfGradp(-mesh.magSf()*fvc::snGrad(p_rgh)); + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + + phase.URef() = + HbyADs[movingPhasei] + + fvc::reconstruct + ( + alphaByADfs[movingPhasei]*mSfGradp + - FgByADfs[movingPhasei] + ); + + phase.URef().correctBoundaryConditions(); + fvConstraints().constrain(phase.URef()); + } + } + + // --- Pressure corrector loop + while (pimple.correct()) + { + // Correct fixed-flux BCs to be consistent with the velocity BCs + fluid_.correctBoundaryFlux(); + + PtrList HbyADs; + PtrList phiHbyADs; + { + // Predicted velocities and fluxes for each phase + PtrList Hs(movingPhases.size()); + PtrList phiHs(movingPhases.size()); + + // Correction force fluxes + PtrList ddtCorrs + ( + momentumTransferSystem_.ddtCorrs() + ); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + + Hs.set + ( + movingPhasei, + UEqns[phase.index()].H() + + byDt + ( + max(phase.residualAlpha() - alpha, scalar(0)) + *phase.rho() + ) + *phase.U()().oldTime() + ); + + if (HVms.set(movingPhasei)) + { + Hs[movingPhasei] += HVms[movingPhasei]; + } + + phiHs.set + ( + movingPhasei, + fvc::flux(Hs[movingPhasei]) + ddtCorrs[phase.index()] + ); + } + + HbyADs = invADVs & Hs; + phiHbyADs = invADVfs & phiHs; + } + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + constrainHbyA(HbyADs[movingPhasei], phase.U(), p_rgh); + constrainPhiHbyA(phiHbyADs[movingPhasei], phase.U(), p_rgh); + + phiHbyADs[movingPhasei] -= FgByADfs[movingPhasei]; + } + + // Total predicted flux + surfaceScalarField phiHbyA + ( + IOobject + ( + "phiHbyA", + runTime.name(), + mesh + ), + mesh, + dimensionedScalar(dimVolumetricFlux, 0) + ); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + phiHbyA += alphafs[phase.index()]*phiHbyADs[movingPhasei]; + } + + MRF.makeRelative(phiHbyA); + fvc::makeRelative(phiHbyA, movingPhases[0].U()); + + // Pressure "diffusivity" + surfaceScalarField rAf + ( + IOobject + ( + "rAf", + runTime.name(), + mesh + ), + mesh, + dimensionedScalar(dimTime/dimDensity, 0) + ); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + rAf += alphafs[phase.index()]*alphaByADfs[movingPhasei]; + } + + // Update the fixedFluxPressure BCs to ensure flux consistency + { + surfaceScalarField::Boundary phib + ( + surfaceScalarField::Internal::null(), + phi.boundaryField() + ); + phib = 0; + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + + phib += + alphafs[phase.index()].boundaryField() + *phase.phi()().boundaryField(); + } + + setSnGrad + ( + p_rgh.boundaryFieldRef(), + ( + phiHbyA.boundaryField() - phib + )/(mesh.magSf().boundaryField()*rAf.boundaryField()) + ); + } + + // Compressible pressure equations + PtrList pEqnComps(compressibilityEqns(dmdts)); + + // 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(rAf, p_rgh) + ); + + // Solve + { + fvScalarMatrix pEqn(pEqnIncomp); + + forAll(phases, phasei) + { + pEqn += pEqnComps[phasei]; + } + + if (fluid.incompressible()) + { + pEqn.setReference + ( + pressureReference.refCell(), + pressureReference.refValue() + ); + } + + fvConstraints().constrain(pEqn); + + pEqn.solve(); + } + + // Correct fluxes and velocities on last non-orthogonal iteration + if (pimple.finalNonOrthogonalIter()) + { + phi_ = phiHbyA + pEqnIncomp.flux(); + + surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAf); + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + + phase.phiRef() = + phiHbyADs[movingPhasei] + + alphaByADfs[movingPhasei]*mSfGradp; + + // Set the phase dilatation rate + phase.divU(-pEqnComps[phase.index()] & p_rgh); + + MRF.makeRelative(phase.phiRef()); + fvc::makeRelative(phase.phiRef(), phase.U()); + } + + // Optionally relax pressure for velocity correction + p_rgh.relax(); + + mSfGradp = pEqnIncomp.flux()/rAf; + + if (dragCorrection) + { + PtrList dragCorrs(movingPhases.size()); + PtrList dragCorrfs(movingPhases.size()); + momentumTransferSystem_.dragCorrs(dragCorrs, dragCorrfs); + + PtrList dragCorrByADs + ( + invADVs & dragCorrs + ); + + PtrList dragCorrByADfs + ( + invADVfs & dragCorrfs + ); + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + + phase.URef() = + HbyADs[movingPhasei] + + fvc::reconstruct + ( + alphaByADfs[movingPhasei]*mSfGradp + - FgByADfs[movingPhasei] + + dragCorrByADfs[movingPhasei] + ) + - dragCorrByADs[movingPhasei]; + } + } + else + { + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + + phase.URef() = + HbyADs[movingPhasei] + + fvc::reconstruct + ( + alphaByADfs[movingPhasei]*mSfGradp + - FgByADfs[movingPhasei] + ); + } + } + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + + phase.URef().correctBoundaryConditions(); + phase.correctUf(); + fvConstraints().constrain(phase.URef()); + } + } + } + + // Update and limit the static pressure + p_ = p_rgh + rho*buoyancy.gh + buoyancy.pRef; + fvConstraints().constrain(p_); + + // Account for static pressure reference + if (p_rgh.needReference() && fluid.incompressible()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pressureReference.refValue() + - getRefCellValue(p, pressureReference.refCell()) + ); + } + + // Limit p_rgh + p_rgh = p - rho*buoyancy.gh - buoyancy.pRef; + + // Update densities from change in p_rgh + forAll(phases, phasei) + { + phaseModel& phase = phases_[phasei]; + if (!phase.incompressible()) + { + phase.rho() += phase.fluidThermo().psi()*(p_rgh - p_rgh_0); + } + } + + // Correct p_rgh for consistency with p and the updated densities + rho = fluid.rho(); + p_rgh = p - rho*buoyancy.gh - buoyancy.pRef; + } + + UEqns.clear(); +} + + +// ************************************************************************* // diff --git a/OFsolvers/OF13/birdmultiphaseEuler/compressibilityEqns.C b/OFsolvers/OF13/birdmultiphaseEuler/compressibilityEqns.C new file mode 100644 index 00000000..652efa5a --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/compressibilityEqns.C @@ -0,0 +1,120 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2022-2025 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 . + +\*---------------------------------------------------------------------------*/ + +#include "birdmultiphaseEuler.H" +#include "fvcDdt.H" +#include "fvcDiv.H" +#include "fvcSup.H" +#include "fvmDdt.H" +#include "fvmDiv.H" +#include "fvmSup.H" + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +Foam::PtrList +Foam::solvers::birdmultiphaseEuler::compressibilityEqns +( + const PtrList& dmdts +) const +{ + volScalarField& p_rgh = p_rgh_; + + PtrList pEqnComps(phases.size()); + + forAll(phases_, phasei) + { + phaseModel& phase = phases_[phasei]; + const volScalarField& alpha = phase; + volScalarField& rho = phase.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, rho) + fvc::div(phase.alphaRhoPhi()) + - fvc::Sp(fvc::ddt(alpha) + fvc::div(phase.alphaPhi()), rho) + )/rho; + } + + // Mesh dilatation correction + if (mesh.moving()) + { + pEqnComp += fvc::div(mesh.phi())*alpha; + } + + // Compressibility + if (!phase.incompressible()) + { + if (pimple.transonic()) + { + const surfaceScalarField phid + ( + IOobject::groupName("phid", phase.name()), + fvc::interpolate(phase.fluidThermo().psi())*phase.phi() + ); + + pEqnComp += + correction + ( + (alpha/rho)* + ( + phase.fluidThermo().psi()*fvm::ddt(p_rgh) + + fvm::div(phid, p_rgh) + - fvm::Sp(fvc::div(phid), p_rgh) + ) + ); + + pEqnComps[phasei].relax(); + } + else + { + pEqnComp += + (alpha*phase.fluidThermo().psi()/rho) + *correction(fvm::ddt(p_rgh)); + } + } + + // Option sources + if (fvModels().addsSupToField(rho.name())) + { + pEqnComp -= fvModels().sourceProxy(alpha, rho, p_rgh)/rho; + } + + // Mass transfer + if (dmdts.set(phasei)) + { + pEqnComp -= dmdts[phasei]/rho; + } + } + + return pEqnComps; +} + + +// ************************************************************************* // diff --git a/OFsolvers/OF13/birdmultiphaseEuler/facePressureCorrector.C b/OFsolvers/OF13/birdmultiphaseEuler/facePressureCorrector.C new file mode 100644 index 00000000..d7022611 --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/facePressureCorrector.C @@ -0,0 +1,401 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2022-2025 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 . + +\*---------------------------------------------------------------------------*/ + +#include "birdmultiphaseEuler.H" +#include "constrainHbyA.H" +#include "constrainPressure.H" +#include "findRefCell.H" +#include "fvcDdt.H" +#include "fvcDiv.H" +#include "fvcSup.H" +#include "fvcSnGrad.H" +#include "fvmDdt.H" +#include "fvmDiv.H" +#include "fvmLaplacian.H" +#include "fvmSup.H" +#include "fvcFlux.H" +#include "fvcMeshPhi.H" +#include "fvcReconstruct.H" + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +void Foam::solvers::birdmultiphaseEuler::facePressureCorrector() +{ + volScalarField& p(p_); + volScalarField& p_rgh = p_rgh_; + + // Face volume fractions + PtrList alphafs(phases.size()); + forAll(phases, phasei) + { + const phaseModel& phase = phases[phasei]; + const volScalarField& alpha = phase; + + alphafs.set(phasei, fvc::interpolate(alpha).ptr()); + alphafs[phasei].rename("pEqn" + alphafs[phasei].name()); + } + + // Diagonal coefficients + rAs.clear(); + if (fluid.implicitPhasePressure()) + { + rAs.setSize(phases.size()); + } + + PtrList HVmfs(movingPhases.size()); + PtrList> invADVfs; + { + PtrList Afs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + + const volScalarField A + ( + byDt + ( + max(alpha.oldTime(), phase.residualAlpha()) + *phase.rho().oldTime() + ) + + UEqns[phase.index()].A() + ); + + if (fluid.implicitPhasePressure()) + { + rAs.set + ( + phase.index(), + new volScalarField + ( + IOobject::groupName("rA", phase.name()), + 1/A + ) + ); + } + + Afs.set + ( + movingPhasei, + new surfaceScalarField + ( + IOobject::groupName("rAf", phase.name()), + fvc::interpolate(A) + ) + ); + } + + invADVfs = momentumTransferSystem_.invADVfs(Afs, HVmfs); + } + + volScalarField rho("rho", fluid.rho()); + + // Phase diagonal coefficients + PtrList alphaByADfs; + PtrList FgByADfs; + { + // Explicit force fluxes + PtrList Ffs(momentumTransferSystem_.Ffs()); + + const surfaceScalarField ghSnGradRho + ( + "ghSnGradRho", + buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf() + ); + + UPtrList movingAlphafs(movingPhases.size()); + PtrList Fgfs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + movingAlphafs.set(movingPhasei, &alphafs[phase.index()]); + + Fgfs.set + ( + movingPhasei, + Ffs[phase.index()] + + alphafs[phase.index()] + *( + ghSnGradRho + - fluid.surfaceTension(phase)*mesh.magSf() + ) + - max(alphafs[phase.index()], phase.residualAlpha()) + *fvc::interpolate(phase.rho() - rho)*(buoyancy.g & mesh.Sf()) + ); + } + + alphaByADfs = invADVfs & movingAlphafs; + FgByADfs = invADVfs & Fgfs; + } + + // Mass transfer rates + PtrList dmdts(populationBalanceSystem_.dmdts()); + + // --- Pressure corrector loop + while (pimple.correct()) + { + // Correct fixed-flux BCs to be consistent with the velocity BCs + fluid_.correctBoundaryFlux(); + + // Predicted fluxes for each phase + PtrList phiHbyADs; + { + PtrList phiHs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + + phiHs.set + ( + movingPhasei, + ( + fvc::interpolate + ( + max(alpha.oldTime(), phase.residualAlpha()) + *phase.rho().oldTime() + ) + *byDt + ( + phase.Uf().valid() + ? (mesh.Sf() & phase.Uf()().oldTime()) + : MRF.absolute(phase.phi()().oldTime()) + ) + + fvc::flux(UEqns[phase.index()].H()) + ) + ); + + if (HVmfs.set(movingPhasei)) + { + phiHs[movingPhasei] += HVmfs[movingPhasei]; + } + } + + phiHbyADs = invADVfs & phiHs; + } + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + constrainPhiHbyA(phiHbyADs[movingPhasei], phase.U(), p_rgh); + + phiHbyADs[movingPhasei] -= FgByADfs[movingPhasei]; + } + + // Total predicted flux + surfaceScalarField phiHbyA + ( + IOobject + ( + "phiHbyA", + runTime.name(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar(dimVolumetricFlux, 0) + ); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + phiHbyA += alphafs[phase.index()]*phiHbyADs[movingPhasei]; + } + + MRF.makeRelative(phiHbyA); + fvc::makeRelative(phiHbyA, movingPhases[0].U()); + + // Construct pressure "diffusivity" + surfaceScalarField rAf + ( + IOobject + ( + "rAf", + runTime.name(), + mesh + ), + mesh, + dimensionedScalar(dimensionSet(-1, 3, 1, 0, 0), 0) + ); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + rAf += alphafs[phase.index()]*alphaByADfs[movingPhasei]; + } + + // Update the fixedFluxPressure BCs to ensure flux consistency + { + surfaceScalarField::Boundary phib + ( + surfaceScalarField::Internal::null(), + phi.boundaryField() + ); + phib = 0; + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + + phib += + alphafs[phase.index()].boundaryField() + *phase.phi()().boundaryField(); + } + + setSnGrad + ( + p_rgh.boundaryFieldRef(), + ( + phiHbyA.boundaryField() - phib + )/(mesh.magSf().boundaryField()*rAf.boundaryField()) + ); + } + + // Compressible pressure equations + PtrList pEqnComps(compressibilityEqns(dmdts)); + + // 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(rAf, p_rgh) + ); + + // Solve + { + fvScalarMatrix pEqn(pEqnIncomp); + + forAll(phases, phasei) + { + pEqn += pEqnComps[phasei]; + } + + if (fluid.incompressible()) + { + pEqn.setReference + ( + pressureReference.refCell(), + pressureReference.refValue() + ); + } + + fvConstraints().constrain(pEqn); + + pEqn.solve(); + } + + // Correct fluxes and velocities on last non-orthogonal iteration + if (pimple.finalNonOrthogonalIter()) + { + phi_ = phiHbyA + pEqnIncomp.flux(); + + surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAf); + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + const label phasei = phase.index(); + + phase.phiRef() = + phiHbyADs[movingPhasei] + + alphaByADfs[movingPhasei]*mSfGradp; + + // Set the phase dilatation rate + phase.divU(-pEqnComps[phasei] & p_rgh); + } + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + + MRF.makeRelative(phase.phiRef()); + fvc::makeRelative(phase.phiRef(), phase.U()); + + phase.URef() = fvc::reconstruct + ( + fvc::absolute(MRF.absolute(phase.phi()), phase.U()) + ); + + phase.URef().correctBoundaryConditions(); + phase.correctUf(); + fvConstraints().constrain(phase.URef()); + } + } + } + + // Update and limit the static pressure + p = p_rgh + rho*buoyancy.gh + buoyancy.pRef; + fvConstraints().constrain(p); + + // Account for static pressure reference + if (p_rgh.needReference() && fluid.incompressible()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pressureReference.refValue() + - getRefCellValue(p, pressureReference.refCell()) + ); + } + + // Limit p_rgh + p_rgh = p - rho*buoyancy.gh - buoyancy.pRef; + + // Update densities from change in p_rgh + forAll(phases, phasei) + { + phaseModel& phase = phases_[phasei]; + if (!phase.incompressible()) + { + phase.rho() += phase.fluidThermo().psi()*(p_rgh - p_rgh_0); + } + } + + // Correct p_rgh for consistency with p and the updated densities + rho = fluid.rho(); + p_rgh = p - rho*buoyancy.gh - buoyancy.pRef; + p_rgh.correctBoundaryConditions(); + } + + UEqns.clear(); +} + + +// ************************************************************************* // diff --git a/OFsolvers/OF13/birdmultiphaseEuler/facePressureTurbulenceCorrector.C b/OFsolvers/OF13/birdmultiphaseEuler/facePressureTurbulenceCorrector.C new file mode 100644 index 00000000..db71a6c7 --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/facePressureTurbulenceCorrector.C @@ -0,0 +1,401 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2022-2025 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 . + +\*---------------------------------------------------------------------------*/ + +#include "birdmultiphaseEuler.H" +#include "constrainHbyA.H" +#include "constrainPressure.H" +#include "findRefCell.H" +#include "fvcDdt.H" +#include "fvcDiv.H" +#include "fvcSup.H" +#include "fvcSnGrad.H" +#include "fvmDdt.H" +#include "fvmDiv.H" +#include "fvmLaplacian.H" +#include "fvmSup.H" +#include "fvcFlux.H" +#include "fvcMeshPhi.H" +#include "fvcReconstruct.H" + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +void Foam::solvers::birdmultiphaseEuler::facePressureTurbulenceCorrector() +{ + volScalarField& p(p_); + volScalarField& p_rgh = p_rgh_; + + // Face volume fractions + PtrList alphafs(phases.size()); + forAll(phases, phasei) + { + const phaseModel& phase = phases[phasei]; + const volScalarField& alpha = phase; + + alphafs.set(phasei, fvc::interpolate(alpha).ptr()); + alphafs[phasei].rename("pEqn" + alphafs[phasei].name()); + } + + // Diagonal coefficients + rAs.clear(); + if (fluid.implicitPhasePressure()) + { + rAs.setSize(phases.size()); + } + + PtrList HVmfs(movingPhases.size()); + PtrList> invADVfs; + { + PtrList Afs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + + const volScalarField A + ( + byDt + ( + max(alpha.oldTime(), phase.residualAlpha()) + *phase.rho().oldTime() + ) + + UEqns[phase.index()].A() + ); + + if (fluid.implicitPhasePressure()) + { + rAs.set + ( + phase.index(), + new volScalarField + ( + IOobject::groupName("rA", phase.name()), + 1/A + ) + ); + } + + Afs.set + ( + movingPhasei, + new surfaceScalarField + ( + IOobject::groupName("rAf", phase.name()), + fvc::interpolate(A) + ) + ); + } + + invADVfs = momentumTransferSystem_.invADVfs(Afs, HVmfs); + } + + volScalarField rho("rho", fluid.rho()); + + // Phase diagonal coefficients + PtrList alphaByADfs; + PtrList FgByADfs; + { + // Explicit force fluxes + PtrList Ffs(momentumTransferSystem_.Ffs()); + + const surfaceScalarField ghSnGradRho + ( + "ghSnGradRho", + buoyancy.ghf*fvc::snGrad(rho)*mesh.magSf() + ); + + UPtrList movingAlphafs(movingPhases.size()); + PtrList Fgfs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + movingAlphafs.set(movingPhasei, &alphafs[phase.index()]); + + Fgfs.set + ( + movingPhasei, + Ffs[phase.index()] + + alphafs[phase.index()] + *( + ghSnGradRho + - fluid.surfaceTension(phase)*mesh.magSf() + ) + - max(alphafs[phase.index()], phase.residualAlpha()) + *fvc::interpolate(phase.rho() - rho)*(buoyancy.g & mesh.Sf()) + ); + } + + alphaByADfs = invADVfs & movingAlphafs; + FgByADfs = invADVfs & Fgfs; + } + + // Mass transfer rates + PtrList dmdts(populationBalanceSystem_.dmdts()); + + // --- Pressure corrector loop + while (pimple.correct()) + { + // Correct fixed-flux BCs to be consistent with the velocity BCs + fluid_.correctBoundaryFlux(); + + // Predicted fluxes for each phase + PtrList phiHbyADs; + { + PtrList phiHs(movingPhases.size()); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + const volScalarField& alpha = phase; + + phiHs.set + ( + movingPhasei, + ( + fvc::interpolate + ( + max(alpha.oldTime(), phase.residualAlpha()) + *phase.rho().oldTime() + ) + *byDt + ( + phase.Uf().valid() + ? (mesh.Sf() & phase.Uf()().oldTime()) + : MRF.absolute(phase.phi()().oldTime()) + ) + + fvc::flux(UEqns[phase.index()].H()) + ) + ); + + if (HVmfs.set(movingPhasei)) + { + phiHs[movingPhasei] += HVmfs[movingPhasei]; + } + } + + phiHbyADs = invADVfs & phiHs; + } + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + constrainPhiHbyA(phiHbyADs[movingPhasei], phase.U(), p_rgh); + + phiHbyADs[movingPhasei] -= FgByADfs[movingPhasei]; + } + + // Total predicted flux + surfaceScalarField phiHbyA + ( + IOobject + ( + "phiHbyA", + runTime.name(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar(dimVolumetricFlux, 0) + ); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + phiHbyA += alphafs[phase.index()]*phiHbyADs[movingPhasei]; + } + + MRF.makeRelative(phiHbyA); + fvc::makeRelative(phiHbyA, movingPhases[0].U()); + + // Construct pressure "diffusivity" + surfaceScalarField rAf + ( + IOobject + ( + "rAf", + runTime.name(), + mesh + ), + mesh, + dimensionedScalar(dimensionSet(-1, 3, 1, 0, 0), 0) + ); + + forAll(movingPhases, movingPhasei) + { + const phaseModel& phase = movingPhases[movingPhasei]; + + rAf += alphafs[phase.index()]*alphaByADfs[movingPhasei]; + } + + // Update the fixedFluxPressure BCs to ensure flux consistency + { + surfaceScalarField::Boundary phib + ( + surfaceScalarField::Internal::null(), + phi.boundaryField() + ); + phib = 0; + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + + phib += + alphafs[phase.index()].boundaryField() + *phase.phi()().boundaryField(); + } + + setSnGrad + ( + p_rgh.boundaryFieldRef(), + ( + phiHbyA.boundaryField() - phib + )/(mesh.magSf().boundaryField()*rAf.boundaryField()) + ); + } + + // Compressible pressure equations + PtrList pEqnComps(compressibilityEqns(dmdts)); + + // 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(rAf, p_rgh) + ); + + // Solve + { + fvScalarMatrix pEqn(pEqnIncomp); + + forAll(phases, phasei) + { + pEqn += pEqnComps[phasei]; + } + + if (fluid.incompressible()) + { + pEqn.setReference + ( + pressureReference.refCell(), + pressureReference.refValue() + ); + } + + fvConstraints().constrain(pEqn); + + pEqn.solve(); + } + + // Correct fluxes and velocities on last non-orthogonal iteration + if (pimple.finalNonOrthogonalIter()) + { + phi_ = phiHbyA + pEqnIncomp.flux(); + + surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAf); + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + const label phasei = phase.index(); + + phase.phiRef() = + phiHbyADs[movingPhasei] + + alphaByADfs[movingPhasei]*mSfGradp; + + // Set the phase dilatation rate + phase.divU(-pEqnComps[phasei] & p_rgh); + } + + forAll(movingPhases, movingPhasei) + { + phaseModel& phase = movingPhases_[movingPhasei]; + + MRF.makeRelative(phase.phiRef()); + fvc::makeRelative(phase.phiRef(), phase.U()); + + phase.URef() = fvc::reconstruct + ( + fvc::absolute(MRF.absolute(phase.phi()), phase.U()) + ); + + phase.URef().correctBoundaryConditions(); + phase.correctUf(); + fvConstraints().constrain(phase.URef()); + } + } + } + + // Update and limit the static pressure + p = p_rgh + rho*buoyancy.gh + buoyancy.pRef; + fvConstraints().constrain(p); + + // Account for static pressure reference + if (p_rgh.needReference() && fluid.incompressible()) + { + p += dimensionedScalar + ( + "p", + p.dimensions(), + pressureReference.refValue() + - getRefCellValue(p, pressureReference.refCell()) + ); + } + + // Limit p_rgh + p_rgh = p - rho*buoyancy.gh - buoyancy.pRef; + + // Update densities from change in p_rgh + forAll(phases, phasei) + { + phaseModel& phase = phases_[phasei]; + if (!phase.incompressible()) + { + phase.rho() += phase.fluidThermo().psi()*(p_rgh - p_rgh_0); + } + } + + // Correct p_rgh for consistency with p and the updated densities + rho = fluid.rho(); + p_rgh = p - rho*buoyancy.gh - buoyancy.pRef; + p_rgh.correctBoundaryConditions(); + } + + UEqns.clear(); +} + + +// ************************************************************************* // diff --git a/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/Make/files b/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/Make/files new file mode 100644 index 00000000..4e0cb801 --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/Make/files @@ -0,0 +1,16 @@ +populationBalanceMoments/populationBalanceMoments.C +populationBalanceSizeDistribution/populationBalanceSizeDistribution.C +populationBalanceSetSizeDistribution/populationBalanceSetSizeDistribution.C + +phaseForces/phaseForces.C +phaseMap/phaseMap.C + +wallBoilingProperty/wallBoilingProperty.C + +adjustTimeStepToNucleation/adjustTimeStepToNucleation.C + +disengagement/disengagement.C + +kLA/kLA.C + +LIB = $(FOAM_USER_LIBBIN)/libbirdmultiphaseEulerFunctionObjects diff --git a/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/Make/options b/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/Make/options new file mode 100644 index 00000000..2fc81d3d --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/Make/options @@ -0,0 +1,26 @@ +EXE_INC = \ + -I../lnInclude \ + -I../phaseSystem/lnInclude \ + -I../populationBalance/lnInclude \ + -I../fvModels/lnInclude \ + -I$(LIB_SRC)/fileFormats/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/functionObjects/field/lnInclude \ + -I$(LIB_SRC)/physicalProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/twoPhaseModels/twoPhaseMixture/lnInclude \ + -I$(LIB_SRC)/twoPhaseModels/compressibleTwoPhases/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \ + -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \ + -I$(LIB_SRC)/fvModels/general/lnInclude + +LIB_LIBS = \ + -L$(FOAM_USER_LIBBIN) \ + -lbirdphaseSystem \ + -lbirdpopulationBalance \ + -lbirdmultiphaseEulerFvModels \ + -lfileFormats \ + -lfieldFunctionObjects diff --git a/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/adjustTimeStepToNucleation/adjustTimeStepToNucleation.C b/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/adjustTimeStepToNucleation/adjustTimeStepToNucleation.C new file mode 100644 index 00000000..98fd564b --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/adjustTimeStepToNucleation/adjustTimeStepToNucleation.C @@ -0,0 +1,152 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2025 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 . + +\*---------------------------------------------------------------------------*/ + +#include "adjustTimeStepToNucleation.H" +#include "fvModels.H" +#include "nucleation.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace functionObjects +{ + defineTypeNameAndDebug(adjustTimeStepToNucleation, 0); + + addToRunTimeSelectionTable + ( + functionObject, + adjustTimeStepToNucleation, + dictionary + ); +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::functionObjects::adjustTimeStepToNucleation::adjustTimeStepToNucleation +( + const word& name, + const Time& runTime, + const dictionary& dict +) +: + fvMeshFunctionObject(name, runTime, dict), + modelName_(word::null), + maxCo_(NaN) +{ + read(dict); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::functionObjects::adjustTimeStepToNucleation::~adjustTimeStepToNucleation() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::functionObjects::adjustTimeStepToNucleation::read +( + const dictionary& dict +) +{ + modelName_ = dict.lookupOrDefault("model", word::null); + maxCo_ = dict.lookupOrDefault("maxCo", 1); + + return true; +} + + +bool Foam::functionObjects::adjustTimeStepToNucleation::execute() +{ + return true; +} + + +bool Foam::functionObjects::adjustTimeStepToNucleation::write() +{ + return true; +} + + +Foam::scalar +Foam::functionObjects::adjustTimeStepToNucleation::maxDeltaT() const +{ + if (!time_.controlDict().lookupOrDefault("adjustTimeStep", false)) + { + return vGreat; + } + + const Foam::fvModels& fvModels = Foam::fvModels::New(mesh_); + + if (modelName_ == word::null) + { + bool found = false; + + tmp tTau = + volScalarField::Internal::New + ( + typedName("tau"), + mesh(), + dimensionedScalar(dimTime, vGreat) + ); + + forAll(fvModels, fvModeli) + { + if (isA(fvModels[fvModeli])) + { + found = true; + + const fv::nucleation& nucleationModel = + refCast(fvModels[fvModeli]); + + tTau = min(tTau, nucleationModel.tau()); + } + } + + if (!found) + { + WarningInFunction + << "No nucleation models found" + << exit(FatalError); + } + + return gMin(tTau()); + } + else + { + const fv::nucleation& nucleationModel = + refCast(fvModels[modelName_]); + + return gMin(nucleationModel.tau()()); + } +} + + +// ************************************************************************* // diff --git a/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/adjustTimeStepToNucleation/adjustTimeStepToNucleation.H b/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/adjustTimeStepToNucleation/adjustTimeStepToNucleation.H new file mode 100644 index 00000000..40de6228 --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/adjustTimeStepToNucleation/adjustTimeStepToNucleation.H @@ -0,0 +1,148 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2025 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 . + +Class + Foam::functionObjects::adjustTimeStepToNucleation + +Description + Controls the time-step to the minimum nucleation time scale. A specific + nucleation model can be specified, otherwise it will scan for and use all + the nucleation fvModels in operation. + + Note that this function only does anything if time step adjustment is + enabled in the controlDict. + + Example of function object specification: + \verbatim + adjustTimeStepToNucleation + { + type adjustTimeStepToNucleation; + libs ("libbirdmultiphaseEulerFunctionObjects.so"); + maxCo 0.1; + } + \endverbatim + +Usage + \table + Property | Description | Required | Default value + type | type name: adjustTimeStepToNucleation | yes | + model | name of the nucleation fvModel | no | + maxCo | maximum nucleation "Courant" number | no | 1 + \endtable + +SourceFiles + adjustTimeStepToNucleation.C + +\*---------------------------------------------------------------------------*/ + +#ifndef adjustTimeStepToNucleation_functionObject_H +#define adjustTimeStepToNucleation_functionObject_H + +#include "fvMeshFunctionObject.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace functionObjects +{ + +/*---------------------------------------------------------------------------*\ + Class adjustTimeStepToNucleation Declaration +\*---------------------------------------------------------------------------*/ + +class adjustTimeStepToNucleation +: + public fvMeshFunctionObject +{ + // Private Data + + //- Name of the nucleation fvModel + word modelName_; + + //- Maximum nucleation "Courant" number + scalar maxCo_; + + +public: + + //- Runtime type information + TypeName("adjustTimeStepToNucleation"); + + + // Constructors + + //- Construct from components + adjustTimeStepToNucleation + ( + const word& name, + const Time& runTime, + const dictionary& dict + ); + + //- Disallow default bitwise copy construction + adjustTimeStepToNucleation(const adjustTimeStepToNucleation&) = delete; + + + // Destructor + virtual ~adjustTimeStepToNucleation(); + + + // Member Functions + + //- Read and reset the timeStep Function1 + virtual bool read(const dictionary&); + + //- Return the list of fields required + virtual wordList fields() const + { + return wordList::null(); + } + + //- Do nothing + virtual bool execute(); + + //- Do nothing + virtual bool write(); + + //- Return the minimum chemistry.deltaTChem() + virtual scalar maxDeltaT() const; + + + // Member Operators + + //- Disallow default bitwise assignment + void operator=(const adjustTimeStepToNucleation&) = delete; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace functionObjects +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/disengagement/disengagement.C b/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/disengagement/disengagement.C new file mode 100644 index 00000000..31df68e3 --- /dev/null +++ b/OFsolvers/OF13/birdmultiphaseEuler/functionObjects/disengagement/disengagement.C @@ -0,0 +1,233 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "disengagement.H" +#include "addToRunTimeSelectionTable.H" +#include "fvcVolumeIntegrate.H" +#include "fixedValueFvPatchField.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace functionObjects +{ + defineTypeNameAndDebug(disengagement, 0); + addToRunTimeSelectionTable(functionObject, disengagement, dictionary); +} +} + +//- Add an element to the back of a list and shift all elements left. +// Removes first element. This allows to only keep the information needed. +template +void add_to_list_like_static_queue(Foam::List& list, const Type& value) +{ + for (Foam::label i = 0; i < list.size() - 1; i++) + { + //- Movel all elements + list[i] = list[i+1]; + } + + //- emplace last element + list[list.size()-1] = value; + +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::functionObjects::disengagement::disengagement +( + const word& name, + const Time& runTime, + const dictionary& dict +) +: + fvMeshFunctionObject(name, runTime, dict), + logFiles(obr_,name), + phases_(mesh_.lookupObject("phaseProperties").phases()), + phaseName_(dict.lookup("phase")), + inletPhaseName_(dict.lookup("inletPhase")), + inletPatch_(dict.lookup("inlet")), + tolerance_(dict.lookup("tolerance")), + nsamples_(dict.lookup