songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ENH: regionFaModels: new filmSeparation models for finite-area framework

Browse Source
This commit is contained in:
Kutalmis Bercin 2024-12-06 18:07:40 +00:00 committed by Andrew Heather
parent 83b9c60a88
commit 51d1050339
14 changed files with 1843 additions and 407 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
src/regionFaModels/Make/files
View File

@ -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
kinematic = liquidFilm/subModels/kinematic
liquidFilm/subModels/kinematic/injectionModel/injectionModelList/injectionModelList.C
liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModel.C
liquidFilm/subModels/kinematic/injectionModel/injectionModel/injectionModelNew.C
$(kinematic)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
$(kinematic)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModelNew.C
$(kinematic)/filmTurbulenceModel/laminar/laminar.C
liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C
liquidFilm/subModels/kinematic/injectionModel/BrunDrippingInjection/BrunDrippingInjection.C
$(kinematic)/injectionModel/injectionModelList/injectionModelList.C
$(kinematic)/injectionModel/injectionModel/injectionModel.C
$(kinematic)/injectionModel/injectionModel/injectionModelNew.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)/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

201
src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.H
View File

@ -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
// ************************************************************************* //

127
src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparation.C Normal file
View File

@ -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();
}
}
}
// ************************************************************************* //

153
src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparation.H Normal file
View File

@ -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
// ************************************************************************* //

593
src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.C Normal file
View File

@ -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
// ************************************************************************* //

288
src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/FriedrichModel/FriedrichModel.H Normal file
View File

@ -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
// ************************************************************************* //

356
src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/curvatureSeparation/curvatureSeparation.C → src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.C
View File

@ -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"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace regionModels
{
namespace areaSurfaceFilmModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(curvatureSeparation, 0);
addToRunTimeSelectionTable
(
injectionModel,
curvatureSeparation,
dictionary
);
namespace Foam
{
namespace filmSeparationModels
{
defineTypeNameAndDebug(OwenRyleyModel, 0);
addToRunTimeSelectionTable(filmSeparationModel, OwenRyleyModel, dictionary);
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
tmp<areaScalarField> curvatureSeparation::calcInvR1
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
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,176 +96,59 @@ 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();
}
}
return clamp(cosAngle, scalarMinMax(-1, 1));
}
// * * * * * * * * * * * * * * * * 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
)
{
const faMesh& mesh = film().regionMesh();
const areaScalarField& delta = 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));
// calculate force balance
scalarField Fnet(mesh.nFaces(), Zero);
scalarField separated(mesh.nFaces(), Zero);
forAll(invR1, i)
{
if ((invR1[i] > minInvR1_) && (delta[i]*invR1[i] > deltaByR1Min_))
{
const scalar R1 = 1.0/(invR1[i] + ROOTVSMALL);
const scalar R2 = R1 + delta[i];
// inertial force
const scalar Fi = -delta[i]*rho[i]*magSqrU[i]*72.0/60.0*invR1[i];
// body force
const scalar Fb =
- 0.5*rho[i]*magG_*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i];
// surface force
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
(
mesh.newIOobject("Fnet"),
mesh,
dimensionedScalar(dimForce, Zero)
);
volFnet.primitiveFieldRef() = Fnet;
volFnet.write();
}
{
areaScalarField areaSeparated
(
mesh.newIOobject("separated"),
mesh,
dimensionedScalar(dimMass, Zero)
);
areaSeparated.primitiveFieldRef() = separated;
areaSeparated.write();
}
{
areaScalarField areaMassToInject
(
mesh.newIOobject("massToInject"),
mesh,
dimensionedScalar(dimMass, Zero)
);
areaMassToInject.primitiveFieldRef() = massToInject;
areaMassToInject.write();
areaCosAngle.primitiveFieldRef() = cosAngle;
areaCosAngle.correctBoundaryConditions();
areaCosAngle.write();
}
{
@ -290,14 +163,123 @@ void curvatureSeparation::correct
}
}
injectionModel::correct();
return clamp(cosAngle, scalarMinMax(-1, 1));
}
tmp<scalarField> OwenRyleyModel::netForce() const
{
const faMesh& mesh = film().regionMesh();
const areaScalarField& h = film().h();
const areaVectorField& U = film().Uf();
const edgeScalarField& phi = film().phi2s();
const areaScalarField& rho = film().rho();
const areaScalarField& sigma = film().sigma();
const scalarField magSqrU(magSqr(U));
const scalarField invR1(calcInvR1(U));
const scalarField cosAngle(calcCosAngle(phi, invR1));
const scalar magG = mag(film().g().value());
// Initialize the net-force magnitude
auto tFnet = tmp<scalarField>::New(mesh.nFaces(), Zero);
auto& Fnet = tFnet.ref();
forAll(Fnet, i)
{
if ((invR1[i] > minInvR1()) && (h[i]*invR1[i] > minHbyR1()))
{
const scalar R1 = 1.0/(invR1[i] + ROOTVSMALL);
const scalar R2 = R1 + h[i];
// Inertial force (OR:Eq. 11; '2' is an exponent in the Eq.)
const scalar Fi = -72.0/60.0*h[i]*rho[i]*magSqrU[i]*invR1[i];
// Body force (OR:Eq. 11)
const scalar Fb =
-0.5*rho[i]*magG*invR1[i]*(sqr(R1) - sqr(R2))*cosAngle[i];
// Surface force (OR:Eq. 11)
const scalar Fs = sigma[i]/R2;
Fnet[i] = Fi + Fb + Fs;
}
}
if (debug && mesh.time().writeTime())
{
{
areaScalarField areaFnet
(
mesh.newIOobject("Fnet"),
mesh,
dimensionedScalar(dimForce, Zero)
);
areaFnet.primitiveFieldRef() = Fnet;
areaFnet.write();
}
}
return tFnet;
}
// * * * * * * * * * * * * * * * * 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)
{
separated[i] = 1;
}
}
return tseparated;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace areaSurfaceFilmModels
} // End namespace regionModels
} // End namespace filmSeparationModels
} // End namespace Foam
// ************************************************************************* //

221
src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/OwenRyleyModel/OwenRyleyModel.H Normal file
View File

@ -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
// ************************************************************************* //

51
src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.C Normal file
View File

@ -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)
{}
// ************************************************************************* //

149
src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModel.H Normal file
View File

@ -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
// ************************************************************************* //

60
src/regionFaModels/liquidFilm/subModels/kinematic/injectionModel/filmSeparation/filmSeparationModels/filmSeparationModel/filmSeparationModelNew.C Normal file
View File

@ -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));
}
// ************************************************************************* //

8
tutorials/lagrangian/kinematicParcelFoam/drippingChair/0/U
View File

@ -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

5
tutorials/lagrangian/kinematicParcelFoam/pitzDailyWithSprinklers/0/U
View File

@ -65,13 +65,14 @@ boundaryField
injectionModels
(
curvatureSeparation
filmSeparation
);
forces ();
curvatureSeparationCoeffs
filmSeparationCoeffs
{
model OwenRyley;
definedPatchRadii 0;
minInvR1 0;
deltaByR1Min 0; // h*invRi = 0.004*10

5
tutorials/lagrangian/reactingParcelFoam/liquidFilmStepWithSprinklers/0/U
View File

@ -68,13 +68,14 @@ boundaryField
injectionModels
(
curvatureSeparation
filmSeparation
);
forces ();
curvatureSeparationCoeffs
filmSeparationCoeffs
{
model OwenRyley;
definedPatchRadii 0;
}