From 0bf7758dd4471bdeb641ed18c73c870e5dbf86e7 Mon Sep 17 00:00:00 2001 From: Johan Roenby Date: Fri, 10 Dec 2021 16:16:39 +0000 Subject: [PATCH 1/4] ENH: interIsoFoam/isoAdvection: new porosity extension MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Niels Gjøl Jacobsen <> Co-authored-by: Konstantinos Missios <> Co-authored-by: Henning Scheufler Co-authored-by: Johan Roenby --- .../solvers/multiphase/interIsoFoam/UEqn.H | 35 ++++++++ .../multiphase/interIsoFoam/UEqnAddPorosity.H | 48 ++++++++++ .../multiphase/interIsoFoam/alphaEqn.H | 14 ++- .../multiphase/interIsoFoam/createFields.H | 1 + .../multiphase/interIsoFoam/createPorosity.H | 34 ++++++++ .../multiphase/interIsoFoam/interIsoFoam.C | 6 +- .../interIsoFoam/porousAlphaCourantNo.H | 69 +++++++++++++++ .../multiphase/interIsoFoam/porousCourantNo.H | 67 ++++++++++++++ .../isoAdvection/isoAdvection.C | 87 +++++++++++++++---- .../isoAdvection/isoAdvection.H | 33 +++++-- .../isoAdvection/isoAdvectionTemplates.C | 54 ++++++++---- .../geometricVoF/cellCuts/cutFace/cutFace.C | 6 +- 12 files changed, 408 insertions(+), 46 deletions(-) create mode 100644 applications/solvers/multiphase/interIsoFoam/UEqn.H create mode 100644 applications/solvers/multiphase/interIsoFoam/UEqnAddPorosity.H create mode 100644 applications/solvers/multiphase/interIsoFoam/createPorosity.H create mode 100644 applications/solvers/multiphase/interIsoFoam/porousAlphaCourantNo.H create mode 100644 applications/solvers/multiphase/interIsoFoam/porousCourantNo.H diff --git a/applications/solvers/multiphase/interIsoFoam/UEqn.H b/applications/solvers/multiphase/interIsoFoam/UEqn.H new file mode 100644 index 0000000000..1f696f0650 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/UEqn.H @@ -0,0 +1,35 @@ + MRF.correctBoundaryVelocity(U); + + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + fvm::div(rhoPhi, U) + + MRF.DDt(rho, U) + + turbulence->divDevRhoReff(rho, U) + == + fvOptions(rho, U) + ); + + #include "UEqnAddPorosity.H" + + UEqn.relax(); + + fvOptions.constrain(UEqn); + + if (pimple.momentumPredictor()) + { + solve + ( + UEqn + == + fvc::reconstruct + ( + ( + mixture.surfaceTensionForce() + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + ) * mesh.magSf() + ) + ); + + fvOptions.correct(U); + } diff --git a/applications/solvers/multiphase/interIsoFoam/UEqnAddPorosity.H b/applications/solvers/multiphase/interIsoFoam/UEqnAddPorosity.H new file mode 100644 index 0000000000..5e0242966e --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/UEqnAddPorosity.H @@ -0,0 +1,48 @@ +// Including porosity effects in UEqn following: +// Jensen, B., Jacobsen, N. G., & Christensen, E. D. (2014). +// Investigations on the porous media equations and resistance +// coefficients for coastal structures. Coastal Engineering, 84, 56-72. + +if (porosityEnabled) +{ + const volScalarField& porosity = tporosity.cref(); + + const word porosityModel("JensenEtAl2014"); + const dictionary& dict = + porosityProperties.subDict(porosityModel + "Coeffs"); + const dimensionedScalar alpha(dimless/dimArea, dict.get("alpha")); + const dimensionedScalar beta(dimless/dimLength, dict.get("beta")); + const dimensionedScalar d50(dimless, dict.get("d50")); + const dimensionedScalar KC(dimless, dict.get("KC")); + + // Generating Darcy-Forchheimer coefficient: F = rho*U*(a + b*|U|) + // Shoud it be mu or muEff in the equation below? + { + // Darcy term + volScalarField DarcyForchheimerCoeff + ( + alpha*sqr(1 - porosity)*mixture.mu()/sqr(porosity)/sqr(d50) + ); + + // Adding Forchheimer term + DarcyForchheimerCoeff += rho*mag(U) + *beta*(1 + pos(KC)*7.5/KC)*(1 - porosity)/sqr(porosity)/d50; + + // Adding Darcy-Forchheimer term as implicit source term + UEqn += fvm::Sp(DarcyForchheimerCoeff, U); + } + + { + // Generating added mass force coefficient + const dimensionedScalar gamma_p(dimless, dict.get("gamma_p")); + const volScalarField Cm(gamma_p*(1 - porosity)); + + UEqn += Cm*fvm::ddt(rho, U); + UEqn += Cm*MRF.DDt(rho, U); + } + + // Dividing both matrix entries and source term by porosity to compensate + // for the fact that the FVM cell volume averages use division by cell + // volume V whereas only the cell pore volume, porosity*V, is accessible. + UEqn *= scalar(1)/porosity; +} diff --git a/applications/solvers/multiphase/interIsoFoam/alphaEqn.H b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H index b6f8aecb2f..3f3280a652 100644 --- a/applications/solvers/multiphase/interIsoFoam/alphaEqn.H +++ b/applications/solvers/multiphase/interIsoFoam/alphaEqn.H @@ -22,8 +22,20 @@ mixture.correct(); } +scalar domainFraction = 0; +if (porosityEnabled) +{ + const volScalarField& porosity = tporosity.cref(); + rhoPhi *= scalar(1)/fvc::interpolate(porosity); + domainFraction = alpha1.weightedAverage(mesh.Vsc()*porosity).value(); +} +else +{ + domainFraction = alpha1.weightedAverage(mesh.Vsc()).value(); +} + Info<< "Phase-1 volume fraction = " - << alpha1.weightedAverage(mesh.Vsc()).value() + << domainFraction << " Min(" << alpha1.name() << ") = " << min(alpha1).value() << " Max(" << alpha1.name() << ") = " << max(alpha1).value() << endl; diff --git a/applications/solvers/multiphase/interIsoFoam/createFields.H b/applications/solvers/multiphase/interIsoFoam/createFields.H index f3cfac66a0..1fdb0243d0 100644 --- a/applications/solvers/multiphase/interIsoFoam/createFields.H +++ b/applications/solvers/multiphase/interIsoFoam/createFields.H @@ -130,5 +130,6 @@ mesh.setFluxRequired(alpha1.name()); #include "createMRF.H" #include "createFvOptions.H" +#include "createPorosity.H" isoAdvection advector(alpha1, phi, U); diff --git a/applications/solvers/multiphase/interIsoFoam/createPorosity.H b/applications/solvers/multiphase/interIsoFoam/createPorosity.H new file mode 100644 index 0000000000..ca7c52d98a --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/createPorosity.H @@ -0,0 +1,34 @@ +// Reading porosity properties from constant directory +IOdictionary porosityProperties +( + IOobject + ( + "porosityProperties", + runTime.constant(), + runTime, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ) +); + +const bool porosityEnabled +( + porosityProperties.getOrDefault("porosityEnabled", false) +); + +tmp tporosity; +if (porosityEnabled) +{ + tporosity = tmp::New + ( + IOobject + ( + "porosity", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); +} diff --git a/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C index 24f20fd6a4..e9c2f4ca4c 100644 --- a/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C +++ b/applications/solvers/multiphase/interIsoFoam/interIsoFoam.C @@ -92,7 +92,7 @@ int main(int argc, char *argv[]) #include "initCorrectPhi.H" #include "createUfIfPresent.H" - #include "CourantNo.H" + #include "porousCourantNo.H" #include "setInitialDeltaT.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -101,8 +101,8 @@ int main(int argc, char *argv[]) while (runTime.run()) { #include "readDyMControls.H" - #include "CourantNo.H" - #include "alphaCourantNo.H" + #include "porousCourantNo.H" + #include "porousAlphaCourantNo.H" #include "setDeltaT.H" ++runTime; diff --git a/applications/solvers/multiphase/interIsoFoam/porousAlphaCourantNo.H b/applications/solvers/multiphase/interIsoFoam/porousAlphaCourantNo.H new file mode 100644 index 0000000000..c4c1fcbe3c --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/porousAlphaCourantNo.H @@ -0,0 +1,69 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 Johan Roenby +------------------------------------------------------------------------------- +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 + porousAlphaCourantNo + +Description + Calculates and outputs the mean and maximum Courant Numbers. + +\*---------------------------------------------------------------------------*/ + +scalar maxAlphaCo = runTime.controlDict().get("maxAlphaCo"); +scalar alphaCoNum = 0; +scalar meanAlphaCoNum = 0; + +if (mesh.nInternalFaces()) +{ + scalarField sumPhi + ( + mixture.nearInterface()().primitiveField() + *fvc::surfaceSum(mag(phi))().primitiveField() + ); + + if (porosityEnabled) + { + const volScalarField& porosity = tporosity.cref(); + + alphaCoNum = + 0.5*gMax(sumPhi/mesh.V().field()/porosity)*runTime.deltaTValue(); + + meanAlphaCoNum = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()*porosity)) + *runTime.deltaTValue(); + } + else + { + alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); + + meanAlphaCoNum = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue(); + } +} + +Info<< "Interface Courant Number mean: " << meanAlphaCoNum + << " max: " << alphaCoNum << endl; + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interIsoFoam/porousCourantNo.H b/applications/solvers/multiphase/interIsoFoam/porousCourantNo.H new file mode 100644 index 0000000000..aac3f02051 --- /dev/null +++ b/applications/solvers/multiphase/interIsoFoam/porousCourantNo.H @@ -0,0 +1,67 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 Johan Roenby +------------------------------------------------------------------------------- +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 + porousCourantNo + +Description + Calculates and outputs the mean and maximum Courant Numbers. + +\*---------------------------------------------------------------------------*/ + +scalar CoNum = 0; +scalar meanCoNum = 0; + +{ + scalarField sumPhi + ( + fvc::surfaceSum(mag(phi))().primitiveField() + ); + + if (porosityEnabled) + { + const volScalarField& porosity = tporosity.cref(); + + CoNum = + 0.5*gMax(sumPhi/mesh.V().field()/porosity) + *runTime.deltaTValue(); + + meanCoNum = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()*porosity)) + *runTime.deltaTValue(); + } + else + { + 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/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C index 13d80e8250..39d57acdfb 100644 --- a/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C +++ b/src/transportModels/geometricVoF/advectionSchemes/isoAdvection/isoAdvection.C @@ -7,8 +7,8 @@ ------------------------------------------------------------------------------- Copyright (C) 2016-2017 DHI Modified code Copyright (C) 2016-2017 OpenCFD Ltd. - Modified code Copyright (C) 2019 Johan Roenby Modified code Copyright (C) 2019-2020 DLR + Modified code Copyright (C) 2018, 2021 Johan Roenby ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -108,6 +108,10 @@ Foam::isoAdvection::isoAdvection bsn0_(bsFaces_.size()), bsUn0_(bsFaces_.size()), + // Porosity + porosityEnabled_(dict_.getOrDefault("porosityEnabled", false)), + porosityPtr_(nullptr), + // Parallel run data procPatchLabels_(mesh_.boundary().size()), surfaceCellFacesOnProcPatches_(0) @@ -132,6 +136,48 @@ Foam::isoAdvection::isoAdvection // Get boundary mesh and resize the list for parallel comms setProcessorPatches(); } + + // Reading porosity properties from constant directory + IOdictionary porosityProperties + ( + IOobject + ( + "porosityProperties", + mesh_.time().constant(), + mesh_, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ) + ); + + porosityEnabled_ = + porosityProperties.getOrDefault("porosityEnabled", false); + + if (porosityEnabled_) + { + if (mesh_.foundObject("porosity")) + { + porosityPtr_ = mesh_.getObjectPtr("porosity"); + + if + ( + gMin(porosityPtr_->primitiveField()) <= 0 + || gMax(porosityPtr_->primitiveField()) > 1 + SMALL + ) + { + FatalErrorInFunction + << "Porosity field has values <= 0 or > 1" + << exit(FatalError); + } + } + else + { + FatalErrorInFunction + << "Porosity enabled in constant/porosityProperties " + << "but no porosity field is found in object registry." + << exit(FatalError); + } + } } @@ -222,6 +268,26 @@ void Foam::isoAdvection::timeIntegratedFlux() DynamicList> isoFacePts; const DynamicField