diff --git a/src/regionFaModels/Make/files b/src/regionFaModels/Make/files index 375a0c49d386c1d214c7ba7a06b6e03acb005c41..a8875dd6ed0677f7a9ae26758f6ba0454e0d0855 100644 --- a/src/regionFaModels/Make/files +++ b/src/regionFaModels/Make/files @@ -15,22 +15,29 @@ derivedFvPatchFields/vibrationShell/vibrationShellFvPatchScalarField.C /* Sub-Model */ -liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C -liquidFilm/subModels/kinematic/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C -liquidFilm/subModels/kinematic/filmTurbulenceModel/laminar/laminar.C - -liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C -liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C -liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C - -liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C -liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C - -liquidFilm/subModels/kinematic/force/forceList/forceList.C -liquidFilm/subModels/kinematic/force/force/force.C -liquidFilm/subModels/kinematic/force/force/forceNew.C -liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C -liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C +kinematic = liquidFilm/subModels/kinematic + +$(kinematic)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C +$(kinematic)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C +$(kinematic)/filmTurbulenceModel/laminar/laminar.C + +$(kinematic)/injectionModel/injectionModelList/injectionModelList.C +$(kinematic)/injectionModel/injectionModel/injectionModel.C +$(kinematic)/injectionModel/injectionModel/injectionModelNew.C + +$(kinematic)/injectionModel/filmSeparation/filmSeparation.C +$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.C +$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModelNew.C +$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.C +$(kinematic)/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C + +$(kinematic)/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C + +$(kinematic)/force/forceList/forceList.C +$(kinematic)/force/force/force.C +$(kinematic)/force/force/forceNew.C +$(kinematic)/force/contactAngleForces/contactAngleForce/contactAngleForce.C +$(kinematic)/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C liquidFilm/subModels/filmSubModelBase.C diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H deleted file mode 100644 index 07dbef9e18a2f8129b333a44fdb066a03ac5cefa..0000000000000000000000000000000000000000 --- a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H +++ /dev/null @@ -1,201 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | www.openfoam.com - \\/ M anipulation | -------------------------------------------------------------------------------- - Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2021 OpenCFD Ltd. -------------------------------------------------------------------------------- -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 <http://www.gnu.org/licenses/>. - -Class - Foam::regionModels::areaSurfaceFilmModels::curvatureSeparation - -Description - Curvature film separation model - - Assesses film curvature via the mesh geometry and calculates a force - balance of the form: - - F_sum = F_inertial + F_body + F_surface - - If F_sum < 0, the film separates. Similarly, if F_sum > 0 the film will - remain attached. - - Reference: - \verbatim - Owen, I., & Ryley, D. J. (1985). - The flow of thin liquid films around corners. - International journal of multiphase flow, 11(1), 51-62. - \endverbatim - -Usage - Example of the model specification: - \verbatim - injectionModels - ( - curvatureSeparation - ); - - curvatureSeparationCoeffs - { - // Optional entries - deltaByR1Min <scalar>; - definedPatchRadii <scalar>; - fThreshold <scalar>; - minInvR1 <scalar>; - - // Inherited entries - ... - } - \endverbatim - - where the entries mean: - \table - Property | Description | Type | Reqd | Deflt - deltaByR1Min | Minimum gravity driven film thickness | scalar | no | 0 - definedPatchRadii | Patch radius i | scalar | no | 0 - fThreshold | Threshold force for separation | scalar | no | 1e-8 - minInvR1 | Minimum inv R1 for separation | scalar | no | 5 - \endtable - - The inherited entries are elaborated in: - - \link injectionModel.H \endlink - -SourceFiles - curvatureSeparation.C - -\*---------------------------------------------------------------------------*/ - -#ifndef curvatureSeparation_H -#define curvatureSeparation_H - -#include "injectionModel.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ -namespace regionModels -{ -namespace areaSurfaceFilmModels -{ - -/*---------------------------------------------------------------------------*\ - Class curvatureSeparation Declaration -\*---------------------------------------------------------------------------*/ - -class curvatureSeparation -: - public injectionModel -{ - // Private Member Functions - - //- No copy construct - curvatureSeparation(const curvatureSeparation&) = delete; - - //- No copy assignment - void operator=(const curvatureSeparation&) = delete; - - -protected: - - // Protected Data - - //- Gradient of surface normals - areaTensorField gradNHat_; - - //- Minimum gravity driven film thickness (non-dimensionalised delta/R1) - scalar deltaByR1Min_; - - //- Patch radius - scalar definedPatchRadii_; - - //- Magnitude of gravity vector - scalar magG_; - - //- Direction of gravity vector - vector gHat_; - - //- Threshold force for separation - scalar fThreshold_; - - //- Minimum inv R1 for separation - scalar minInvR1_; - - - // Protected Member Functions - - //- Calculate local (inverse) radius of curvature - tmp<areaScalarField> calcInvR1 - ( - const areaVectorField& U, - const scalarField& calcCosAngle - ) const; - - //- Calculate the cosine of the angle between gravity vector and - //- cell out flow direction - tmp<scalarField> calcCosAngle(const edgeScalarField& phi) const; - - -public: - - //- Runtime type information - TypeName("curvatureSeparation"); - - - // Constructors - - //- Construct from surface film model - curvatureSeparation - ( - liquidFilmBase& film, - const dictionary& dict - ); - - - //- Destructor - virtual ~curvatureSeparation() = default; - - - // Member Functions - - // Evolution - - //- Correct - virtual void correct - ( - scalarField& availableMass, - scalarField& massToInject, - scalarField& diameterToInject - ); -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace areaSurfaceFilmModels -} // End namespace regionModels -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparation.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparation.C new file mode 100644 index 0000000000000000000000000000000000000000..105092d9644e528a8703fe0ce277c262c3a629d4 --- /dev/null +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparation.C @@ -0,0 +1,127 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2011-2017 OpenFOAM Foundation + Copyright (C) 2020-2024 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "filmSeparation.H" +#include "filmSeparationModel.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace regionModels +{ +namespace areaSurfaceFilmModels +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(filmSeparation, 0); +addToRunTimeSelectionTable +( + injectionModel, + filmSeparation, + dictionary +); +} +} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::regionModels::areaSurfaceFilmModels::filmSeparation::filmSeparation +( + liquidFilmBase& film, + const dictionary& dict +) +: + injectionModel(type(), film, dict), + filmSeparationModelPtr_(filmSeparationModel::New(film, coeffDict_)) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::regionModels::areaSurfaceFilmModels::filmSeparation::~filmSeparation() +{} // filmSeparationModel was forward declared + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::regionModels::areaSurfaceFilmModels::filmSeparation::correct +( + scalarField& availableMass, + scalarField& massToInject, + scalarField& diameterToInject +) +{ + const faMesh& mesh = film().regionMesh(); + + // Calculate the mass ratio of film separation + tmp<scalarField> tmassRatio = filmSeparationModelPtr_->separatedMassRatio(); + const auto& massRatio = tmassRatio.cref(); + + // Update various film properties based on the mass ratio + massToInject = massRatio*availableMass; + diameterToInject = massRatio*film().h(); + availableMass -= massRatio*availableMass; + + addToInjectedMass(sum(massToInject)); + + injectionModel::correct(); + + + if (debug && mesh.time().writeTime()) + { + { + areaScalarField areaSeparated + ( + mesh.newIOobject("separated"), + mesh, + dimensionedScalar(dimMass, Zero) + ); + areaSeparated.primitiveFieldRef() = massRatio; + areaSeparated.write(); + } + + { + areaScalarField areaMassToInject + ( + mesh.newIOobject("massToInject"), + mesh, + dimensionedScalar(dimMass, Zero) + ); + areaMassToInject.primitiveFieldRef() = massToInject; + areaMassToInject.write(); + } + } +} + + +// ************************************************************************* // diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparation.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparation.H new file mode 100644 index 0000000000000000000000000000000000000000..178a70558d6b8bebc39ad41a90b8fec35740270a --- /dev/null +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparation.H @@ -0,0 +1,153 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2011-2017 OpenFOAM Foundation + Copyright (C) 2021-2024 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +Class + Foam::regionModels::areaSurfaceFilmModels::filmSeparation + +Description + The \c filmSeparation is a collection of curvature thin-film separation + models designed to predict the onset of film separation and mass separation + in geometries featuring sharp and/or rounded corners. + +Usage + Minimal example by using boundary-condition files: + \verbatim + injectionModels + { + // Mandatory entries + filmSeparation + } + + filmSeparationCoeffs + { + model <word>; + + // Conditional entries + + // Option-1: when model == OwenRyley + + // Option-2: when model == Friedrich + + // Inherited entries + ... + } + \endverbatim + + where the entries mean: + \table + Property | Description | Type | Reqd | Deflt + model | Name of the filmSeparation model | word | yes | - + \endtable + + Options for the \c model entry: + \verbatim + OwenRyley | Model proposed by Owen-Ryley (1985) + Friedrich | Model proposed by Friedrich et al. (2008) + \endverbatim + + The inherited entries are elaborated in: + - \link injectionModel.H \endlink + +SourceFiles + filmSeparation.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Foam_filmSeparation_H +#define Foam_filmSeparation_H + +#include "injectionModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward Declarations +class filmSeparationModel; + +namespace regionModels +{ +namespace areaSurfaceFilmModels +{ + +/*---------------------------------------------------------------------------*\ + Class filmSeparation Declaration +\*---------------------------------------------------------------------------*/ + +class filmSeparation +: + public injectionModel +{ + // Private Data + + //- Film-separation model + autoPtr<filmSeparationModel> filmSeparationModelPtr_; + + +public: + + //- Runtime type information + TypeName("filmSeparation"); + + + // Constructors + + //- Construct from base film model and dictionary + filmSeparation + ( + liquidFilmBase& film, + const dictionary& dict + ); + + + //- Destructor + virtual ~filmSeparation(); + + + // Member Functions + + //- Correct film properties due to the film separation + virtual void correct + ( + scalarField& availableMass, + scalarField& massToInject, + scalarField& diameterToInject + ); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace areaSurfaceFilmModels +} // End namespace regionModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C new file mode 100644 index 0000000000000000000000000000000000000000..b9490ebc672677458c13ea84000510d17865e4dc --- /dev/null +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C @@ -0,0 +1,593 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2024 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "FriedrichModel.H" +#include "processorFaPatch.H" +#include "unitConversion.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace filmSeparationModels +{ + defineTypeNameAndDebug(FriedrichModel, 0); + addToRunTimeSelectionTable(filmSeparationModel, FriedrichModel, dictionary); + + +const Foam::Enum +< + FriedrichModel::separationType +> +FriedrichModel::separationTypeNames +({ + { separationType::FULL, "full" }, + { separationType::PARTIAL , "partial" }, +}); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +bitSet FriedrichModel::calcCornerEdges() const +{ + bitSet cornerEdges(mesh().nEdges(), false); + + const areaVectorField& faceCentres = mesh().areaCentres(); + const areaVectorField& faceNormals = mesh().faceAreaNormals(); + + const labelUList& own = mesh().edgeOwner(); + const labelUList& nbr = mesh().edgeNeighbour(); + + // Check if internal face-normal vectors diverge (no separation) + // or converge (separation may occur) + forAll(nbr, edgei) + { + const label faceO = own[edgei]; + const label faceN = nbr[edgei]; + + cornerEdges[edgei] = isCornerEdgeSharp + ( + faceCentres[faceO], + faceCentres[faceN], + faceNormals[faceO], + faceNormals[faceN] + ); + } + + + // Skip the rest of the routine if the simulation is a serial run + if (!Pstream::parRun()) return cornerEdges; + + // Check if processor face-normal vectors diverge (no separation) + // or converge (separation may occur) + const faBoundaryMesh& patches = mesh().boundary(); + + for (const faPatch& fap : patches) + { + if (isA<processorFaPatch>(fap)) + { + const label patchi = fap.index(); + const auto& edgeFaces = fap.edgeFaces(); + const label internalEdgei = fap.start(); + + const auto& faceCentresp = faceCentres.boundaryField()[patchi]; + const auto& faceNormalsp = faceNormals.boundaryField()[patchi]; + + forAll(faceNormalsp, bndEdgei) + { + const label faceO = edgeFaces[bndEdgei]; + const label meshEdgei = internalEdgei + bndEdgei; + + cornerEdges[meshEdgei] = isCornerEdgeSharp + ( + faceCentres[faceO], + faceCentresp[bndEdgei], + faceNormals[faceO], + faceNormalsp[bndEdgei] + ); + } + } + } + + return cornerEdges; +} + + +bool FriedrichModel::isCornerEdgeSharp +( + const vector& faceCentreO, + const vector& faceCentreN, + const vector& faceNormalO, + const vector& faceNormalN +) const +{ + // Calculate the relative position of centres of faces sharing an edge + const vector relativePosition(faceCentreN - faceCentreO); + + // Calculate the relative normal of faces sharing an edge + const vector relativeNormal(faceNormalN - faceNormalO); + + // Return true if the face normals converge, meaning that the edge is sharp + return ((relativeNormal & relativePosition) < -1e-8); +} + + +scalarList FriedrichModel::calcCornerAngles() const +{ + scalarList cornerAngles(mesh().nEdges(), Zero); + + const areaVectorField& faceNormals = mesh().faceAreaNormals(); + + const labelUList& own = mesh().edgeOwner(); + const labelUList& nbr = mesh().edgeNeighbour(); + + // Process internal edges + forAll(nbr, edgei) + { + if (!cornerEdges_[edgei]) continue; + + const label faceO = own[edgei]; + const label faceN = nbr[edgei]; + + cornerAngles[edgei] = calcCornerAngle + ( + faceNormals[faceO], + faceNormals[faceN] + ); + } + + + // Skip the rest of the routine if the simulation is a serial run + if (!Pstream::parRun()) return cornerAngles; + + // Process processor edges + const faBoundaryMesh& patches = mesh().boundary(); + + for (const faPatch& fap : patches) + { + if (isA<processorFaPatch>(fap)) + { + const label patchi = fap.index(); + const auto& edgeFaces = fap.edgeFaces(); + const label internalEdgei = fap.start(); + + const auto& faceNormalsp = faceNormals.boundaryField()[patchi]; + + forAll(faceNormalsp, bndEdgei) + { + const label faceO = edgeFaces[bndEdgei]; + const label meshEdgei = internalEdgei + bndEdgei; + + if (!cornerEdges_[meshEdgei]) continue; + + cornerAngles[meshEdgei] = calcCornerAngle + ( + faceNormals[faceO], + faceNormalsp[bndEdgei] + ); + } + } + } + + return cornerAngles; +} + + +scalar FriedrichModel::calcCornerAngle +( + const vector& faceNormalO, + const vector& faceNormalN +) const +{ + const scalar magFaceNormal = mag(faceNormalO)*mag(faceNormalN); + + // Avoid any potential exceptions during the cosine calculations + if (magFaceNormal < SMALL) return 0; + + scalar cosAngle = (faceNormalO & faceNormalN)/magFaceNormal; + cosAngle = clamp(cosAngle, -1, 1); + + return std::acos(cosAngle); +} + + +bitSet FriedrichModel::calcSeparationFaces() const +{ + bitSet separationFaces(mesh().faces().size(), false); + + const edgeScalarField& phis = film().phi2s(); + + const labelUList& own = mesh().edgeOwner(); + const labelUList& nbr = mesh().edgeNeighbour(); + + // Process internal faces + forAll(nbr, edgei) + { + if (!cornerEdges_[edgei]) continue; + + const label faceO = own[edgei]; + const label faceN = nbr[edgei]; + + isSeparationFace + ( + separationFaces, + phis[edgei], + faceO, + faceN + ); + } + + + // Skip the rest of the routine if the simulation is a serial run + if (!Pstream::parRun()) return separationFaces; + + // Process processor faces + const faBoundaryMesh& patches = mesh().boundary(); + + for (const faPatch& fap : patches) + { + if (isA<processorFaPatch>(fap)) + { + const label patchi = fap.index(); + const auto& edgeFaces = fap.edgeFaces(); + const label internalEdgei = fap.start(); + + const auto& phisp = phis.boundaryField()[patchi]; + + forAll(phisp, bndEdgei) + { + const label faceO = edgeFaces[bndEdgei]; + const label meshEdgei(internalEdgei + bndEdgei); + + if (!cornerEdges_[meshEdgei]) continue; + + isSeparationFace + ( + separationFaces, + phisp[bndEdgei], + faceO + ); + } + } + } + + return separationFaces; +} + + +void FriedrichModel::isSeparationFace +( + bitSet& separationFaces, + const scalar phiEdge, + const label faceO, + const label faceN +) const +{ + const scalar tol = 1e-8; + + // Assuming there are no sources/sinks at the edge + if (phiEdge > tol) // From owner to neighbour + { + separationFaces[faceO] = true; + } + else if ((phiEdge < -tol) && (faceN != -1)) // From neighbour to owner + { + separationFaces[faceN] = true; + } +} + + +scalarList FriedrichModel::calcSeparationAngles +( + const bitSet& separationFaces +) const +{ + scalarList separationAngles(mesh().faces().size(), Zero); + + const labelUList& own = mesh().edgeOwner(); + const labelUList& nbr = mesh().edgeNeighbour(); + + // Process internal faces + forAll(nbr, edgei) + { + if (!cornerEdges_[edgei]) continue; + + const label faceO = own[edgei]; + const label faceN = nbr[edgei]; + + if (separationFaces[faceO]) + { + separationAngles[faceO] = cornerAngles_[edgei]; + } + + if (separationFaces[faceN]) + { + separationAngles[faceN] = cornerAngles_[edgei]; + } + } + + + // Skip the rest of the routine if the simulation is a serial run + if (!Pstream::parRun()) return separationAngles; + + // Process processor faces + const edgeScalarField& phis = film().phi2s(); + const faBoundaryMesh& patches = mesh().boundary(); + + for (const faPatch& fap : patches) + { + if (isA<processorFaPatch>(fap)) + { + const label patchi = fap.index(); + const auto& edgeFaces = fap.edgeFaces(); + const label internalEdgei = fap.start(); + + const auto& phisp = phis.boundaryField()[patchi]; + + forAll(phisp, bndEdgei) + { + const label faceO = edgeFaces[bndEdgei]; + const label meshEdgei(internalEdgei + bndEdgei); + + if (!cornerEdges_[meshEdgei]) continue; + + if (separationFaces[faceO]) + { + separationAngles[faceO] = cornerAngles_[meshEdgei]; + } + } + } + } + + return separationAngles; +} + + +tmp<scalarField> FriedrichModel::Fratio() const +{ + const areaVectorField Up(film().Up()); + const areaVectorField& Uf = film().Uf(); + const areaScalarField& h = film().h(); + const areaScalarField& rho = film().rho(); + const areaScalarField& mu = film().mu(); + const areaScalarField& sigma = film().sigma(); + + // Identify the faces where separation may occur + const bitSet separationFaces(calcSeparationFaces()); + + // Calculate the corner angles corresponding to the separation faces + const scalarList separationAngles(calcSeparationAngles(separationFaces)); + + // Initialize the force ratio + auto tFratio = tmp<scalarField>::New(mesh().faces().size(), Zero); + auto& Fratio = tFratio.ref(); + + // Process internal faces + forAll(separationFaces, i) + { + // Skip the routine if the face is not a candidate for separation + if (!separationFaces[i]) continue; + + // Calculate the corner-angle trigonometric values + const scalar sinAngle = std::sin(separationAngles[i]); + const scalar cosAngle = std::cos(separationAngles[i]); + + // Reynolds number (FLW:Eq. 16) + const scalar Re = h[i]*mag(Uf[i])*rho[i]/mu[i]; + + // Weber number (FLW:Eq. 17) + const scalar We = + h[i]*rhop_*sqr(mag(Up[i]) - mag(Uf[i]))/(2.0*sigma[i]); + + // Characteristic breakup length (FLW:Eq. 15) + const scalar Lb = + 0.0388*Foam::sqrt(h[i])*Foam::pow(Re, 0.6)*Foam::pow(We, -0.5); + + // Force ratio - denominator (FLW:Eq. 20) + const scalar den = + sigma[i]*(sinAngle + 1.0) + rho[i]*magG_*h[i]*Lb*cosAngle; + + if (mag(den) > 0) + { + // Force ratio (FLW:Eq. 20) + Fratio[i] = rho[i]*sqr(mag(Uf[i]))*h[i]*sinAngle/den; + } + } + + + // Skip the rest of the routine if the simulation is a serial run + if (!Pstream::parRun()) return tFratio; + + // Process processor faces + const faBoundaryMesh& patches = mesh().boundary(); + + for (const faPatch& fap : patches) + { + if (isA<processorFaPatch>(fap)) + { + const label patchi = fap.index(); + const label internalEdgei = fap.start(); + + const auto& hp = h.boundaryField()[patchi]; + const auto& Ufp = Uf.boundaryField()[patchi]; + const auto& Upp = Up.boundaryField()[patchi]; + const auto& rhop = rho.boundaryField()[patchi]; + const auto& sigmap = sigma.boundaryField()[patchi]; + const auto& mup = mu.boundaryField()[patchi]; + + forAll(hp, i) + { + // Skip the routine if the face is not a candidate for separation + if (!separationFaces[i]) continue; + + const label meshEdgei = internalEdgei + i; + + // Calculate the corner-angle trigonometric values + const scalar sinAngle = std::sin(cornerAngles_[meshEdgei]); + const scalar cosAngle = std::cos(cornerAngles_[meshEdgei]); + + // Reynolds number (FLW:Eq. 16) + const scalar Re = hp[i]*mag(Ufp[i])*rhop[i]/mup[i]; + + // Weber number (FLW:Eq. 17) + const scalar We = + hp[i]*rhop_*sqr(mag(Upp[i]) - mag(Ufp[i]))/(2.0*sigmap[i]); + + // Characteristic breakup length (FLW:Eq. 15) + const scalar Lb = + 0.0388*Foam::sqrt(hp[i]) + *Foam::pow(Re, 0.6)*Foam::pow(We, -0.5); + + // Force ratio - denominator (FLW:Eq. 20) + const scalar den = + sigmap[i]*(sinAngle + 1.0) + + rhop[i]*magG_*hp[i]*Lb*cosAngle; + + if (mag(den) > 0) + { + // Force ratio (FLW:Eq. 20) + Fratio[i] = rhop[i]*sqr(mag(Ufp[i]))*hp[i]*sinAngle/den; + } + } + } + } + + return tFratio; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +FriedrichModel::FriedrichModel +( + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film, + const dictionary& dict +) +: + filmSeparationModel(film, dict), + separation_ + ( + separationTypeNames.getOrDefault + ( + "separationType", + dict, + separationType::FULL + ) + ), + rhop_(dict.getScalar("rhop")), + magG_(mag(film.g().value())), + C0_(dict.getOrDefault<scalar>("C0", 0.882)), + C1_(dict.getOrDefault<scalar>("C1", -1.908)), + C2_(dict.getOrDefault<scalar>("C2", 1.264)), + cornerEdges_(calcCornerEdges()), + cornerAngles_(calcCornerAngles()) +{ + if (rhop_ < VSMALL) + { + FatalIOErrorInFunction(dict) + << "Primary-phase density, rhop: " << rhop_ << " must be non-zero." + << abort(FatalIOError); + } + + if (mag(C2_) < VSMALL) + { + FatalIOErrorInFunction(dict) + << "Empirical constant, C2 = " << C2_ << "cannot be zero." + << abort(FatalIOError); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +tmp<scalarField> FriedrichModel::separatedMassRatio() const +{ + tmp<scalarField> tFratio = Fratio(); + const auto& Fratio = tFratio.cref(); + + // Initialize the mass ratio of film separation + auto tseparated = tmp<scalarField>::New(mesh().faces().size(), Zero); + auto& separated = tseparated.ref(); + + + switch (separation_) + { + case separationType::FULL: + { + forAll(Fratio, i) + { + if (Fratio[i] > 1) + { + separated[i] = 1; + } + } + break; + } + case separationType::PARTIAL: + { + forAll(Fratio, i) + { + if (Fratio[i] > 1) + { + // (ZJD:Eq. 16) + separated[i] = C0_ + C1_*Foam::exp(-Fratio[i]/C2_); + } + } + break; + } + default: + break; // This should not happen. + } + + if (debug && mesh().time().writeTime()) + { + { + areaScalarField areaFratio + ( + mesh().newIOobject("Fratio"), + mesh(), + dimensionedScalar(dimForce, Zero) + ); + areaFratio.primitiveFieldRef() = Fratio; + areaFratio.write(); + } + } + + + return tseparated; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace filmSeparationModels +} // End namespace Foam + + +// ************************************************************************* // + diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.H new file mode 100644 index 0000000000000000000000000000000000000000..93aeba20916592e91fc6b0f20b21df7bafb30fd6 --- /dev/null +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.H @@ -0,0 +1,288 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2024 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +Class + Foam::filmSeparationModels::FriedrichModel + +Description + Computes film-separation properties from sharp edges for full separation + (Friedrich et al., 2008) and partial separation (Zhang et al., 2018). + + The governing equations for full separation (Friedrich et al., 2008): + + \f[ + F_{ratio} = + \frac{\rho \, |\u|^2 \, h \, sin(\theta)} + {\sigma (1 + sin(\theta)) + \rho \, \mathbf{g}\, h\,L_b cos(\theta)} + \f] + + with: + + \f[ + L_b = 0.0388 h^{0.5} \mathrm{Re}^{0.6} \mathrm{We}^{-0.5} + \f] + + \f[ + \mathrm{Re} = \frac{h \, |\u| \, \rho}{\mu} + \f] + + \f[ + \mathrm{We} = \frac{h \, \rho_p (\u_p - \u)^2}{2 \sigma} + \f] + + where: + + \vartable + F_{ratio} | Force ratio + h | Film thickness + \rho | Film density + \rho_p | Primary-phase (gas) density + \sigma | Film surface tension + \mu | Film dynamic viscosity + \mathbf{u} | Film velocity + \mathbf{g} | Gravitational accelaration + \theta | Sharp-corner angle + L_b | Characteristic breakup length + \endvartable + + The onset of film separation is trigerred and the film is assumed + fully separated when \f$F_{ratio}>1\f$; otherwise, it remains attached. + + + The governing equations for partial separation (Zhang et al., 2018): + + \f[ + m_{ratio} = C_0 + C_1 \, exp\left(-\frac{F_{ratio}}{C_2}\right) + \f] + + where: + + \vartable + m_{ratio} | Mass fraction of separated film mass + C_0 | Empirical constant (0.882) + C_1 | Empirical constant (-1.908) + C_2 | Empirical constant (1.264) + \endvartable + + With the above model modification, the film separation begins when + \f$F_{ratio}>1\f$; however, only the portion with \f$m_{ratio}\f$ + separates while the rest stays attached. + + + Reference: + \verbatim + Governing equations for the full film-separation model (tag:FLW): + Friedrich, M. A., Lan, H., Wegener, J. L., + Drallmeier, J. A., & Armaly, B. F. (2008). + A separation criterion with experimental + validation for shear-driven films in separated flows. + J. Fluids Eng. May 2008, 130(5): 051301. + DOI:10.1115/1.2907405 + + Governing equations for the partial film-separation model (tag:ZJD): + Zhang, Y., Jia, M., Duan, H., Wang, P., + Wang, J., Liu, H., & Xie, M. (2018). + Numerical and experimental study of spray impingement and liquid + film separation during the spray/wall interaction at expanding corners. + International Journal of Multiphase Flow, 107, 67-81. + DOI:10.1016/j.ijmultiphaseflow.2018.05.016 + \endverbatim + +Usage + Minimal example in boundary-condition files: + \verbatim + filmSeparationCoeffs + { + // Mandatory entries + model Friedrich; + rhop <scalar>; + + // Optional entries + separationType <word>; + } + \endverbatim + + where the entries mean: + \table + Property | Description | Type | Reqd | Deflt + model | Model name: Friedrich | word | yes | - + rhop | Primary-phase density | scalar | yes | - + separationType | Film separation type | word | no | full + \endtable + + Options for the \c separationType entry: + \verbatim + full | Full film separation (Friedrich et al., 2008) + partial | Partial film separation (Zhang et al., 2018) + \endverbatim + +SourceFiles + FriedrichModel.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Foam_filmSeparationModels_FriedrichModel_H +#define Foam_filmSeparationModels_FriedrichModel_H + +#include "filmSeparationModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace filmSeparationModels +{ + +/*---------------------------------------------------------------------------*\ + Class FriedrichModel Declaration +\*---------------------------------------------------------------------------*/ + +class FriedrichModel +: + public filmSeparationModel +{ + // Private Enumerations + + //- Options for the film separation type + enum separationType : char + { + FULL = 0, //!< "Full film separation" + PARTIAL //!< "Partial film separation" + }; + + //- Names for separationType + static const Enum<separationType> separationTypeNames; + + + // Private Data + + //- Film separation type + enum separationType separation_; + + //- Approximate uniform mass density of primary phase + scalar rhop_; + + //- Magnitude of the gravitational acceleration + scalar magG_; + + //- Empirical constant for the partial separation model + scalar C0_; + + //- Empirical constant for the partial separation model + scalar C1_; + + //- Empirical constant for the partial separation model + scalar C2_; + + //- List of flags identifying sharp-corner edges where separation + //- may occur + bitSet cornerEdges_; + + //- Corner angles of sharp-corner edges where separation may occur + scalarList cornerAngles_; + + + // Private Member Functions + + //- Return Boolean list of identified sharp-corner edges + bitSet calcCornerEdges() const; + + //- Return true if the given edge is identified as sharp + bool isCornerEdgeSharp + ( + const vector& faceCentreO, + const vector& faceCentreN, + const vector& faceNormalO, + const vector& faceNormalN + ) const; + + //- Return the list of sharp-corner angles for each edge + scalarList calcCornerAngles() const; + + //- Return the sharp-corner angle for a given edge + scalar calcCornerAngle + ( + const vector& faceNormalO, + const vector& faceNormalN + ) const; + + //- Return Boolean list of identified separation faces + bitSet calcSeparationFaces() const; + + //- Return true if the given face is identified as a separation face + void isSeparationFace + ( + bitSet& separationFaces, + const scalar phiEdge, + const label faceO, + const label faceN = -1 + ) const; + + //- Return the list of sharp-corner angles for each face + scalarList calcSeparationAngles(const bitSet& separationFaces) const; + + //- Return the film-separation force ratio per face + tmp<scalarField> Fratio() const; + + +public: + + //- Runtime type information + TypeName("Friedrich"); + + + // Constructors + + //- Construct from the base film model and dictionary + FriedrichModel + ( + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film, + const dictionary& dict + ); + + + // Destructor + virtual ~FriedrichModel() = default; + + + // Member Functions + + // Evaluation + + //- Calculate the mass ratio of film separation + virtual tmp<scalarField> separatedMassRatio() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace filmSeparationModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.C similarity index 50% rename from src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C rename to src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.C index dd3f931c23a8a0ead04fa5215e434c305cff7dbf..3db260f1b8293f5b748257430186d602e434b25c 100644 --- a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2020-2024 OpenCFD Ltd. + Copyright (C) 2021-2024 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -26,63 +26,50 @@ License \*---------------------------------------------------------------------------*/ -#include "curvatureSeparation.H" +#include "OwenRyleyModel.H" +#include "processorFaPatch.H" #include "addToRunTimeSelectionTable.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { -namespace regionModels -{ -namespace areaSurfaceFilmModels +namespace filmSeparationModels { + defineTypeNameAndDebug(OwenRyleyModel, 0); + addToRunTimeSelectionTable(filmSeparationModel, OwenRyleyModel, dictionary); -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -defineTypeNameAndDebug(curvatureSeparation, 0); -addToRunTimeSelectionTable -( - injectionModel, - curvatureSeparation, - dictionary -); - -// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -tmp<areaScalarField> curvatureSeparation::calcInvR1 +tmp<areaScalarField> OwenRyleyModel::calcInvR1 ( - const areaVectorField& U, - const scalarField& calcCosAngle + const areaVectorField& U ) const { const dimensionedScalar smallU(dimVelocity, ROOTVSMALL); const areaVectorField UHat(U/(mag(U) + smallU)); + auto tinvR1 = areaScalarField::New ( "invR1", IOobjectOption::NO_REGISTER, (UHat & (UHat & -gradNHat_)) ); - scalarField& invR1 = tinvR1.ref().primitiveFieldRef(); - // apply defined patch radii - const scalar rMin = 1e-6; - const scalar definedInvR1 = 1.0/max(rMin, definedPatchRadii_); - + // Apply defined patch radii if (definedPatchRadii_ > 0) { - invR1 = definedInvR1; + invR1 = 1.0/max(1e-6, definedPatchRadii_); } - // filter out large radii - const scalar rMax = 1e6; - forAll(invR1, i) + // Filter out large radii + for (auto& x : invR1) { - if ((mag(invR1[i]) < 1/rMax)) + if (mag(x) < 1e-6) { - invR1[i] = -1.0; + x = -1; } } @@ -90,15 +77,18 @@ tmp<areaScalarField> curvatureSeparation::calcInvR1 } -tmp<scalarField> curvatureSeparation::calcCosAngle +tmp<scalarField> OwenRyleyModel::calcCosAngle ( - const edgeScalarField& phi + const edgeScalarField& phi, + const scalarField& invR1 ) const { const areaVectorField& U = film().Uf(); const dimensionedScalar smallU(dimVelocity, ROOTVSMALL); const areaVectorField UHat(U/(mag(U) + smallU)); + const vector gHat(normalised(film().g().value())); + const faMesh& mesh = film().regionMesh(); const labelUList& own = mesh.edgeOwner(); const labelUList& nbr = mesh.edgeNeighbour(); @@ -106,40 +96,70 @@ tmp<scalarField> curvatureSeparation::calcCosAngle scalarField phiMax(mesh.nFaces(), -GREAT); scalarField cosAngle(UHat.size(), Zero); - const scalarField invR1(calcInvR1(U, cosAngle)); - + // Internal edges forAll(nbr, edgei) { const label cellO = own[edgei]; const label cellN = nbr[edgei]; + const scalar phiEdge = phi[edgei]; - if (phi[edgei] > phiMax[cellO]) + if (phiEdge > phiMax[cellO]) { - phiMax[cellO] = phi[edgei]; - cosAngle[cellO] = -gHat_ & UHat[cellN]; + phiMax[cellO] = phiEdge; + cosAngle[cellO] = -gHat & UHat[cellN]; } - if (-phi[edgei] > phiMax[cellN]) + if (-phiEdge > phiMax[cellN]) { - phiMax[cellN] = -phi[edgei]; - cosAngle[cellN] = -gHat_ & UHat[cellO]; + phiMax[cellN] = -phiEdge; + cosAngle[cellN] = -gHat & UHat[cellO]; + } + } + + // Processor edges + for (const auto& phip : phi.boundaryField()) + { + if (isA<processorFaPatch>(phip.patch())) + { + const auto& edgeFaces = phip.patch().edgeFaces(); + const auto& UHatp = UHat.boundaryField()[phip.patch().index()]; + forAll(phip, edgei) + { + const label cellO = edgeFaces[edgei]; + if (phip[edgei] > phiMax[cellO]) + { + phiMax[cellO] = phip[edgei]; + cosAngle[cellO] = -gHat & UHatp[edgei]; + } + } } } cosAngle *= pos(invR1); - // checks + if (debug && mesh.time().writeTime()) { { - areaScalarField volCosAngle + areaScalarField areaCosAngle ( mesh.newIOobject("cosAngle"), mesh, dimensionedScalar(dimless, Zero) ); - volCosAngle.primitiveFieldRef() = cosAngle; - volCosAngle.correctBoundaryConditions(); - volCosAngle.write(); + areaCosAngle.primitiveFieldRef() = cosAngle; + areaCosAngle.correctBoundaryConditions(); + areaCosAngle.write(); + } + + { + areaScalarField areaInvR1 + ( + mesh.newIOobject("InvR1"), + mesh, + dimensionedScalar(inv(dimLength), Zero) + ); + areaInvR1.primitiveFieldRef() = invR1; + areaInvR1.write(); } } @@ -147,157 +167,119 @@ tmp<scalarField> curvatureSeparation::calcCosAngle } -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -curvatureSeparation::curvatureSeparation -( - liquidFilmBase& film, - const dictionary& dict -) -: - injectionModel(type(), film, dict), - gradNHat_(fac::grad(film.regionMesh().faceAreaNormals())), - deltaByR1Min_(coeffDict_.getOrDefault<scalar>("deltaByR1Min", 0)), - definedPatchRadii_ - ( - coeffDict_.getOrDefault<scalar>("definedPatchRadii", 0) - ), - magG_(mag(film.g().value())), - gHat_(Zero), - fThreshold_ - ( - coeffDict_.getOrDefault<scalar>("fThreshold", 1e-8) - ), - minInvR1_ - ( - coeffDict_.getOrDefault<scalar>("minInvR1", 5) - ) -{ - if (magG_ < ROOTVSMALL) - { - FatalErrorInFunction - << "Acceleration due to gravity must be non-zero" - << exit(FatalError); - } - - gHat_ = film.g().value()/magG_; -} - - -// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // - -void curvatureSeparation::correct -( - scalarField& availableMass, - scalarField& massToInject, - scalarField& diameterToInject -) +tmp<scalarField> OwenRyleyModel::netForce() const { const faMesh& mesh = film().regionMesh(); - const areaScalarField& delta = film().h(); + const areaScalarField& h = film().h(); const areaVectorField& U = film().Uf(); const edgeScalarField& phi = film().phi2s(); const areaScalarField& rho = film().rho(); - const scalarField magSqrU(magSqr(film().Uf())); const areaScalarField& sigma = film().sigma(); - const scalarField cosAngle(calcCosAngle(phi)); - const scalarField invR1(calcInvR1(U, cosAngle)); + const scalarField magSqrU(magSqr(U)); + const scalarField invR1(calcInvR1(U)); + const scalarField cosAngle(calcCosAngle(phi, invR1)); + const scalar magG = mag(film().g().value()); - // calculate force balance - scalarField Fnet(mesh.nFaces(), Zero); - scalarField separated(mesh.nFaces(), Zero); + // Initialize the net-force magnitude + auto tFnet = tmp<scalarField>::New(mesh.nFaces(), Zero); + auto& Fnet = tFnet.ref(); - forAll(invR1, i) + forAll(Fnet, i) { - if ((invR1[i] > minInvR1_) && (delta[i]*invR1[i] > deltaByR1Min_)) + if ((invR1[i] > minInvR1()) && (h[i]*invR1[i] > minHbyR1())) { const scalar R1 = 1.0/(invR1[i] + ROOTVSMALL); - const scalar R2 = R1 + delta[i]; + const scalar R2 = R1 + h[i]; - // inertial force - const scalar Fi = -delta[i]*rho[i]*magSqrU[i]*72.0/60.0*invR1[i]; + // Inertial force (OR:Eq. 11; '2' is an exponent in the Eq.) + const scalar Fi = -4.0/3.0*h[i]*rho[i]*magSqrU[i]*invR1[i]; - // body force + // Body force (OR:Eq. 11) const scalar Fb = - - 0.5*rho[i]*magG_*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i]; + -0.5*rho[i]*magG*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i]; - // surface force + // Surface force (OR:Eq. 11) const scalar Fs = sigma[i]/R2; Fnet[i] = Fi + Fb + Fs; - - if (Fnet[i] + fThreshold_ < 0) - { - separated[i] = 1.0; - } } } - // inject all available mass - massToInject = separated*availableMass; - diameterToInject = separated*delta; - availableMass -= separated*availableMass; - - addToInjectedMass(sum(massToInject)); - if (debug && mesh.time().writeTime()) { { - areaScalarField volFnet + areaScalarField areaFnet ( mesh.newIOobject("Fnet"), mesh, dimensionedScalar(dimForce, Zero) ); - volFnet.primitiveFieldRef() = Fnet; - volFnet.write(); + areaFnet.primitiveFieldRef() = Fnet; + areaFnet.write(); } + } - { - areaScalarField areaSeparated - ( - mesh.newIOobject("separated"), - mesh, - dimensionedScalar(dimMass, Zero) - ); - areaSeparated.primitiveFieldRef() = separated; - areaSeparated.write(); - } + return tFnet; +} - { - areaScalarField areaMassToInject - ( - mesh.newIOobject("massToInject"), - mesh, - dimensionedScalar(dimMass, Zero) - ); - areaMassToInject.primitiveFieldRef() = massToInject; - areaMassToInject.write(); - } +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +OwenRyleyModel::OwenRyleyModel +( + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film, + const dictionary& dict +) +: + filmSeparationModel(film, dict), + fThreshold_(dict.getOrDefault<scalar>("fThreshold", 1e-8)), + definedPatchRadii_(dict.getOrDefault<scalar>("definedPatchRadii", 0)), + minHbyR1_(dict.getOrDefault<scalar>("deltaByR1Min", 0)), + minInvR1_(dict.getOrDefault<scalar>("minInvR1", 5)), + gradNHat_(fac::grad(film.regionMesh().faceAreaNormals())) +{ + if (mag(film.g().value()) < ROOTVSMALL) + { + FatalErrorInFunction + << "Gravitational acceleration must be non-zero." + << exit(FatalError); + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +tmp<scalarField> OwenRyleyModel::separatedMassRatio() const +{ + const faMesh& mesh = film().regionMesh(); + + tmp<scalarField> tFnet = netForce(); + const auto& Fnet = tFnet.cref(); + + // Initialize the mass ratio of film separation + auto tseparated = tmp<scalarField>::New(mesh.nFaces(), Zero); + auto& separated = tseparated.ref(); + + forAll(Fnet, i) + { + if ((Fnet[i] + fThreshold_) < 0) { - areaScalarField areaInvR1 - ( - mesh.newIOobject("InvR1"), - mesh, - dimensionedScalar(inv(dimLength), Zero) - ); - areaInvR1.primitiveFieldRef() = invR1; - areaInvR1.write(); + separated[i] = 1; } } - injectionModel::correct(); + return tseparated; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -} // End namespace areaSurfaceFilmModels -} // End namespace regionModels +} // End namespace filmSeparationModels } // End namespace Foam + // ************************************************************************* // + diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.H new file mode 100644 index 0000000000000000000000000000000000000000..0dcedb5d27bacede97dbf42ee3eb507b4aaae239 --- /dev/null +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.H @@ -0,0 +1,221 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2011-2017 OpenFOAM Foundation + Copyright (C) 2021-2024 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +Class + Foam::filmSeparationModels::OwenRyleyModel + +Description + Computes film-separation properties from round edges for full separation + (Owen & Ryley, 1983). The model assesses the film curvature based on the + mesh geometry and flow field, and then calculates a force balance of the + following form: + + \f[ + F_{net} = F_{inertial} + F_{body} + F_{surface} + \f] + + with: + + \f[ + F_{inertial} = -\frac{4}{3} \frac{h \, \rho \, |\mathbf{u}|^2}{R_1} + \f] + + \f[ + F_{body} = + -\frac{\rho |\mathbf{g}| (R_1^2 - R_2^2) cos(\theta)}{2 \, R_1} + \f] + + \f[ + F_{surface} = \frac{\sigma}{R_2} + \f] + + where: + + \vartable + F_{net} | Magnitude of net force on the film + F_{inertial} | Magnitude of inertial forces on the film + F_{body} | Magnitude of body forces on the film + F_{surface} | Magnitude of surface forces on the film + h | Film thickness + \rho | Film density + \sigma | Film surface tension + \mathbf{u} | Film velocity + \mathbf{g} | Gravitational accelaration + \theta | Approximate angle between gravity vector and outflow direction + R_1 | Radius of curvature + R_2 | Sum of the radius of curvature and film thickness + \endvartable + + The film is considered separated when \f$F_{net}<0\f$; otherwise, it + remains attached. + + Reference: + \verbatim + Governing equations (tag:OR): + Owen, I., & Ryley, D. J. (1983). + The flow of thin liquid films around corners. + International journal of multiphase flow, 11(1), 51-62. + DOI:10.1016/0301-9322(85)90005-9 + \endverbatim + +Usage + Minimal example in boundary-condition files: + \verbatim + filmSeparationCoeffs + { + // Mandatory entries + model OwenRyley; + + // Optional entries + fThreshold <scalar>; + definedPatchRadii <scalar>; + deltaByR1Min <scalar>; + minInvR1 <scalar>; + } + \endverbatim + + where the entries mean: + \table + Property | Description | Type | Reqd | Deflt + model | Model name: OwenRyley | word | yes | - + fThreshold | Threshold force for separation | scalar | no | 1e-8 + definedPatchRadii | Patch radius | scalar | no | 0 + deltaByR1Min | Minimum gravity driven film thickness | scalar | no | 0 + minInvR1 | Minimum inv R1 for separation | scalar | no | 5 + \endtable + +Note + - The assumptions underlying the derivation of the model indicate that + the model is better suited for round corners than for sharp corners. + +SourceFiles + OwenRyleyModel.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Foam_filmSeparationModels_OwenRyleyModel_H +#define Foam_filmSeparationModels_OwenRyleyModel_H + +#include "filmSeparationModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace filmSeparationModels +{ + +/*---------------------------------------------------------------------------*\ + Class OwenRyleyModel Declaration +\*---------------------------------------------------------------------------*/ + +class OwenRyleyModel +: + public filmSeparationModel +{ + // Private Data + + //- Threshold force for separation + scalar fThreshold_; + + //- Patch radius + scalar definedPatchRadii_; + + //- Minimum gravity driven film thickness (non-dimensionalised h/R1) + scalar minHbyR1_; + + //- Minimum inv R1 of film velocity for film separation + scalar minInvR1_; + + //- Gradient of surface normals + areaTensorField gradNHat_; + + + // Private Member Functions + + //- Calculate inverse of (local) radius of curvature of film velocity + tmp<areaScalarField> calcInvR1(const areaVectorField& U) const; + + //- Calculate the cosine of the angle between gravity vector and + //- cell-out flow direction + tmp<scalarField> calcCosAngle + ( + const edgeScalarField& phi, + const scalarField& invR1 + ) const; + + //- Calculate the magnitude of net force on the film + tmp<scalarField> netForce() const; + + +public: + + //- Runtime type information + TypeName("OwenRyley"); + + + // Constructors + + //- Construct from the base film model and dictionary + OwenRyleyModel + ( + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film, + const dictionary& dict + ); + + + // Destructor + virtual ~OwenRyleyModel() = default; + + + // Member Functions + + // Access + + // Return access to minimum dimensionless gravity driven film thickness + scalar minHbyR1() const noexcept { return minHbyR1_; } + + // Return access to minimum inv R1 for film separation + scalar minInvR1() const noexcept { return minInvR1_; } + + + // Evaluation + + //- Calculate the mass ratio of film separation + virtual tmp<scalarField> separatedMassRatio() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace filmSeparationModels +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.C new file mode 100644 index 0000000000000000000000000000000000000000..ac27612dc8b47681947d98e224a524ef8144c0a6 --- /dev/null +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.C @@ -0,0 +1,51 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2024 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "filmSeparationModel.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(filmSeparationModel, 0); + defineRunTimeSelectionTable(filmSeparationModel, dictionary); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::filmSeparationModel::filmSeparationModel +( + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film, + const dictionary& dict +) +: + film_(film) +{} + + +// ************************************************************************* // diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.H b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.H new file mode 100644 index 0000000000000000000000000000000000000000..18d402ff64db3e459296bae565697f8473a995fc --- /dev/null +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.H @@ -0,0 +1,149 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2024 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +Namespace + Foam::filmSeparationModels + +Description + A namespace for various \c filmSeparation model implementations. + +Class + Foam::filmSeparationModel + +Description + A base class for \c filmSeparation models. + +SourceFiles + filmSeparationModel.C + filmSeparationModelNew.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Foam_filmSeparationModel_H +#define Foam_filmSeparationModel_H + +#include "liquidFilmBase.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class filmSeparationModel Declaration +\*---------------------------------------------------------------------------*/ + +class filmSeparationModel +{ + // Private Data + + //- Const reference to the film properties + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film_; + + +public: + + //- Runtime type information + TypeName("filmSeparationModel"); + + + // Declare runtime constructor selection table + + declareRunTimeSelectionTable + ( + autoPtr, + filmSeparationModel, + dictionary, + ( + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film, + const dictionary& dict + ), + (film, dict) + ); + + + // Selectors + + //- Return a reference to the selected filmSeparation model + static autoPtr<filmSeparationModel> New + ( + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film, + const dictionary& dict + ); + + + // Generated Methods + + //- No copy construct + filmSeparationModel(const filmSeparationModel&) = delete; + + //- No copy assignment + void operator=(const filmSeparationModel&) = delete; + + + // Constructors + + //- Construct from the base film model and dictionary + filmSeparationModel + ( + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film, + const dictionary& dict + ); + + + //- Destructor + virtual ~filmSeparationModel() = default; + + + // Member Functions + + // Access + + //- Return const access to the film properties + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film() const + { + return film_; + } + + //- Return const access to the finite-area mesh + const faMesh& mesh() const noexcept { return film_.regionMesh(); } + + + // Evaluation + + //- Calculate the mass ratio of film separation + virtual tmp<scalarField> separatedMassRatio() const = 0; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModelNew.C b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModelNew.C new file mode 100644 index 0000000000000000000000000000000000000000..2dcff7fb2224073705d122286ccbea29a99096be --- /dev/null +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModelNew.C @@ -0,0 +1,60 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2024 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "filmSeparationModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::autoPtr<Foam::filmSeparationModel> Foam::filmSeparationModel::New +( + const regionModels::areaSurfaceFilmModels::liquidFilmBase& film, + const dictionary& dict +) +{ + const word modelType(dict.getWord("model")); + + Info<< " Selecting film-separation model: " + << modelType << nl << endl; + + auto* ctorPtr = dictionaryConstructorTable(modelType); + + if (!ctorPtr) + { + FatalIOErrorInLookup + ( + dict, + "filmSeparationModel", + modelType, + *dictionaryConstructorTablePtr_ + ) << exit(FatalIOError); + } + + return autoPtr<filmSeparationModel>(ctorPtr(film, dict)); +} + + +// ************************************************************************* // diff --git a/tutorials/lagrangian/kinematicParcelFoam/drippingChair/0/U b/tutorials/lagrangian/kinematicParcelFoam/drippingChair/0/U index 645aa805d882aa78e85a216c4e77ee36b26bf938..1f280276f689c7c62f7e22f816426288e4d35cd5 100644 --- a/tutorials/lagrangian/kinematicParcelFoam/drippingChair/0/U +++ b/tutorials/lagrangian/kinematicParcelFoam/drippingChair/0/U @@ -73,14 +73,18 @@ boundaryField injectionModels ( - // curvatureSeparation + // filmSeparation // BrunDrippingInjection ); forces (); - curvatureSeparationCoeffs + filmSeparationCoeffs { + // Mandatory entries + model OwenRyley; + + // Optional entries definedPatchRadii 0; minInvR1 0; deltaByR1Min 0; // h*invRi = 0.004*10 diff --git a/tutorials/lagrangian/kinematicParcelFoam/pitzDailyWithSprinklers/0/U b/tutorials/lagrangian/kinematicParcelFoam/pitzDailyWithSprinklers/0/U index 309578ad576262274f6080252f8c1d3f7be1d90b..6c0a5414e4d9c01662630b130ea730bf60b0c767 100644 --- a/tutorials/lagrangian/kinematicParcelFoam/pitzDailyWithSprinklers/0/U +++ b/tutorials/lagrangian/kinematicParcelFoam/pitzDailyWithSprinklers/0/U @@ -65,13 +65,14 @@ boundaryField injectionModels ( - curvatureSeparation + filmSeparation ); forces (); - curvatureSeparationCoeffs + filmSeparationCoeffs { + model OwenRyley; definedPatchRadii 0; minInvR1 0; deltaByR1Min 0; // h*invRi = 0.004*10 diff --git a/tutorials/lagrangian/reactingParcelFoam/liquidFilmStepWithSprinklers/0/U b/tutorials/lagrangian/reactingParcelFoam/liquidFilmStepWithSprinklers/0/U index cf80a66fa7880a0fe445c5c563f1bd7b8e03cfe8..baa0f348367568558bda863a1f87aa250bf20626 100644 --- a/tutorials/lagrangian/reactingParcelFoam/liquidFilmStepWithSprinklers/0/U +++ b/tutorials/lagrangian/reactingParcelFoam/liquidFilmStepWithSprinklers/0/U @@ -68,13 +68,14 @@ boundaryField injectionModels ( - curvatureSeparation + filmSeparation ); forces (); - curvatureSeparationCoeffs + filmSeparationCoeffs { + model OwenRyley; definedPatchRadii 0; }