zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'feature-contact-line-movement' into 'develop'

Browse Source 
ENH: dynamicContactAngleForce: new finite-area contact-angle force model

See merge request Development/openfoam!584
This commit is contained in:
Andrew Heather
2022-12-09 13:09:21 +00:00
parent 8f63fd5230 91a62ca05a
commit 36f66606a2
8 changed files with 264 additions and 145 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/regionFaModels/Make/files
View File

@ -30,7 +30,7 @@ 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/perturbedTemperatureDependent/perturbedTemperatureDependentContactAngleForce.C
liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C
liquidFilm/subModels/filmSubModelBase.C

124
src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C
View File

@ -44,13 +44,6 @@ namespace areaSurfaceFilmModels
defineTypeNameAndDebug(contactAngleForce, 0);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void contactAngleForce::initialise()
{
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
contactAngleForce::contactAngleForce
@ -62,27 +55,7 @@ contactAngleForce::contactAngleForce
:
force(typeName, film, dict),
Ccf_(coeffDict_.get<scalar>("Ccf")),
mask_
(
IOobject
(
typeName + ":fContactForceMask",
film.primaryMesh().time().timeName(),
film.primaryMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
film.regionMesh(),
dimensionedScalar("mask", dimless, 1.0)
)
{
initialise();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
contactAngleForce::~contactAngleForce()
hCrit_(coeffDict_.getOrDefault<scalar>("hCrit", GREAT))
{}
@ -90,84 +63,103 @@ contactAngleForce::~contactAngleForce()
tmp<faVectorMatrix> contactAngleForce::correct(areaVectorField& U)
{
tmp<areaVectorField> tForce
auto tforce = tmp<areaVectorField>::New
(
new areaVectorField
IOobject
(
IOobject
(
typeName + ":contactForce",
film().primaryMesh().time().timeName(),
film().primaryMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
film().regionMesh(),
dimensionedVector(dimForce/dimDensity/dimArea, Zero)
)
IOobject::scopedName(typeName, "contactForce"),
film().primaryMesh().time().timeName(),
film().primaryMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
film().regionMesh(),
dimensionedVector(dimForce/dimDensity/dimArea, Zero)
);
vectorField& force = tForce.ref().primitiveFieldRef();
vectorField& force = tforce.ref().primitiveFieldRef();
const labelUList& own = film().regionMesh().owner();
const labelUList& nbr = film().regionMesh().neighbour();
const labelUList& nei = film().regionMesh().neighbour();
const DimensionedField<scalar, areaMesh>& magSf = film().regionMesh().S();
const scalarField& magSff = magSf.field();
const edgeScalarField& invDx = film().regionMesh().deltaCoeffs();
const areaScalarField& sigma = film().sigma();
const areaScalarField& mu = film().mu();
const areaScalarField& rhof = film().rho();
const areaVectorField& Uf = film().Uf();
const areaScalarField& hf = film().h();
tmp<areaScalarField> talpha = film().alpha();
const areaScalarField& sigma = film().sigma();
const areaScalarField& rhof = film().rho();
const areaScalarField& alpha = talpha();
tmp<areaScalarField> ttheta = theta();
const areaScalarField& theta = ttheta();
const areaVectorField gradAlpha(fac::grad(talpha()));
tmp<areaVectorField> tgradAlpha = fac::grad(alpha);
const areaVectorField& gradAlpha = tgradAlpha();
forAll(nbr, edgei)
forAll(nei, edgei)
{
const label faceO = own[edgei];
const label faceN = nbr[edgei];
const label faceN = nei[edgei];
label facei = -1;
if ((talpha()[faceO] > 0.5) && (talpha()[faceN] < 0.5))
if ((alpha[faceO] > 0.5) && (alpha[faceN] < 0.5))
{
facei = faceO;
}
else if ((talpha()[faceO] < 0.5) && (talpha()[faceN] > 0.5))
else if ((alpha[faceO] < 0.5) && (alpha[faceN] > 0.5))
{
facei = faceN;
}
if (facei != -1 && mask_[facei] > 0.5)
if (facei != -1)
{
const scalar invDx = film().regionMesh().deltaCoeffs()[edgei];
const vector n
(
gradAlpha[facei]/(mag(gradAlpha[facei]) + ROOTVSMALL)
);
const vector n(normalised(gradAlpha[facei]));
const scalar cosTheta = cos(degToRad(theta[facei]));
// (MHDX:Eq. 13)
force[facei] +=
Ccf_*n*sigma[facei]*(1 - cosTheta)/invDx/rhof[facei];
Ccf_*n*sigma[facei]*(1 - cosTheta)
/invDx[edgei]/rhof[facei]/magSff[facei];
// (NDPC:Eq. 11)
if (hf[facei] > hCrit_)
{
force[facei] -= mu[facei]*Uf[facei]/hCrit_;
}
}
}
for (const faPatchScalarField& sigmaBf : sigma.boundaryField())
{
const faPatch& p = sigmaBf.patch();
force /= magSf.field();
if (!p.coupled())
{
const labelUList& faces = p.edgeFaces();
forAll(sigmaBf, edgei)
{
const label face0 = faces[edgei];
force[face0] = Zero;
}
}
}
if (film().regionMesh().time().writeTime())
{
tForce().write();
tforce().write();
gradAlpha.write();
}
tmp<faVectorMatrix> tfvm
(
new faVectorMatrix(U, dimForce/dimDensity)
);
auto tfvm = tmp<faVectorMatrix>::New(U, dimForce/dimDensity);
tfvm.ref() += tForce;
tfvm.ref() += tforce;
return tfvm;
}

84
src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -32,13 +32,82 @@ Description
The effect of the contact angle force can be ignored over a specified
distance from patches.
The contact-angle force is implemented as follows:
\f[
\vec{f}_\theta =
\beta \frac{\sigma (1 - cos(\theta))}{\Delta_{cl}} \vec{n}_{cl}
\frac{1}{\rho_f|\vec{S}_f|}
- \vec{f}_{cl}
\f]
with if \f$ h > h_{crit} \f$:
\f[
\vec{f}_{cl} = \nu_f \frac{\vec{U}_f}{h_{crit}}
\f]
where
\vartable
\vec{f}_\theta | Contact-angle force [N/m^2/rho]
\beta | Empirical constant [-]
\sigma | Liquid-air interfacial surface tension [kg/s^2]
\theta | Contact angle [rad]
\Delta_{cl} | Inverse width in direction normal to contact line [1/m]
\vec{n}_{cl} | Direction normal to contact line [-]
\rho_f | Film density [kg/m^3]
\vec{S}_f | Face-area vector [m^2]
\vec{f}_{cl} | Contact-line movement force [N/m^2/rho]
\nu_f | Film kinematic viscosity [m^2/s]
\vec{U}_f | Film velocity [m/s]
h_{crit} | Critical film height [m]
\endvartable
Reference:
\verbatim
Governing equations (tag:MHDX):
Meredith, K. V., Heather, A., De Vries, J., & Xin, Y. (2011).
A numerical model for partially-wetted flow of thin liquid films.
Computational Methods in Multiphase Flow VI, 70, 239.
Contact line movement (tag:NDPC):
Novák, M., Devaradja, R., Papper, J., & Černý, M. (2020).
Efficient CFD methods for assessment of water management.
In 20. Internationales Stuttgarter Symposium (pp. 151-170).
Springer Vieweg, Wiesbaden.
\endverbatim
Usage
Minimal example:
\verbatim
{
// Mandatory entries
Ccf <scalar>;
// Optional entries
hCrit <scalar>;
// Inherited entries
...
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
Ccf | Empirical coefficient | scalar | yes | -
hCrit | Critical film height [m] | scalar | no | GREAT
\endtable
The inherited entries are elaborated in:
- \link force.H \endlink
SourceFiles
contactAngleForce.C
\*---------------------------------------------------------------------------*/
#ifndef contactAngleForce_H
#define contactAngleForce_H
#ifndef areaSurfaceFilmModels_contactAngleForce_H
#define areaSurfaceFilmModels_contactAngleForce_H
#include "force.H"
@ -64,15 +133,12 @@ class contactAngleForce
//- Coefficient applied to the contact angle force
scalar Ccf_;
//- Mask over which force is applied
areaScalarField mask_;
//- Critical film thickness [m]
scalar hCrit_;
// Private Member Functions
//- Initialise
void initialise();
//- No copy construct
contactAngleForce(const contactAngleForce&) = delete;
@ -104,7 +170,7 @@ public:
//- Destructor
virtual ~contactAngleForce();
virtual ~contactAngleForce() = default;
// Member Functions

93
src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/perturbedTemperatureDependent/perturbedTemperatureDependentContactAngleForce.C → src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -25,7 +25,7 @@ License
\*---------------------------------------------------------------------------*/
#include "perturbedTemperatureDependentContactAngleForce.H"
#include "dynamicContactAngleForce.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -39,26 +39,44 @@ namespace areaSurfaceFilmModels
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(perturbedTemperatureDependentContactAngleForce, 0);
defineTypeNameAndDebug(dynamicContactAngleForce, 0);
addToRunTimeSelectionTable
(
force,
perturbedTemperatureDependentContactAngleForce,
dynamicContactAngleForce,
dictionary
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
perturbedTemperatureDependentContactAngleForce::
perturbedTemperatureDependentContactAngleForce
dynamicContactAngleForce::dynamicContactAngleForce
(
liquidFilmBase& film,
const dictionary& dict
)
:
contactAngleForce(typeName, film, dict),
thetaPtr_(Function1<scalar>::New("theta", coeffDict_, &film.primaryMesh())),
U_vs_thetaPtr_
(
Function1<scalar>::NewIfPresent
(
"Utheta",
coeffDict_,
word::null,
&film.primaryMesh()
)
),
T_vs_thetaPtr_
(
Function1<scalar>::NewIfPresent
(
"Ttheta",
coeffDict_,
word::null,
&film.primaryMesh()
)
),
rndGen_(label(0)),
distribution_
(
@ -68,43 +86,53 @@ perturbedTemperatureDependentContactAngleForce
rndGen_
)
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
perturbedTemperatureDependentContactAngleForce::
~perturbedTemperatureDependentContactAngleForce()
{}
{
if (U_vs_thetaPtr_ && T_vs_thetaPtr_)
{
FatalIOErrorInFunction(dict)
<< "Entries Utheta and Ttheta could not be used together"
<< abort(FatalIOError);
}
}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
tmp<areaScalarField> perturbedTemperatureDependentContactAngleForce::
theta() const
tmp<areaScalarField> dynamicContactAngleForce::theta() const
{
tmp<areaScalarField> ttheta
auto ttheta = tmp<areaScalarField>::New
(
new areaScalarField
IOobject
(
IOobject
(
typeName + ":theta",
film().primaryMesh().time().timeName(),
film().primaryMesh()
),
film().regionMesh(),
dimensionedScalar(dimless, Zero)
)
IOobject::scopedName(typeName, "theta"),
film().primaryMesh().time().timeName(),
film().primaryMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
film().regionMesh(),
dimensionedScalar(dimless, Zero)
);
areaScalarField& theta = ttheta.ref();
scalarField& thetai = theta.ref();
const areaScalarField& T = film().Tf();
// Initialize with the function of temperature
thetai = thetaPtr_->value(T());
if (U_vs_thetaPtr_)
{
// Initialize with the function of film speed
const areaVectorField& U = film().Uf();
thetai = U_vs_thetaPtr_->value(mag(U()));
}
if (T_vs_thetaPtr_)
{
// Initialize with the function of film temperature
const areaScalarField& T = film().Tf();
thetai = T_vs_thetaPtr_->value(T());
}
// Add the stochastic perturbation
forAll(thetai, facei)
@ -115,6 +143,7 @@ theta() const
return ttheta;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace surfaceFilmModels

88
src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/perturbedTemperatureDependent/perturbedTemperatureDependentContactAngleForce.H → src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -24,31 +24,66 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::regionModels::areaSurfaceFilmModels::
perturbedTemperatureDependentContactAngleForce
Foam::regionModels::areaSurfaceFilmModels::dynamicContactAngleForce
Description
Temperature dependent contact angle force with a stochastic perturbation.
Film-speed or film-temperature dependent
contact-angle force with a stochastic perturbation.
The contact angle in degrees is specified as a \c Foam::Function1 type,
to enable the use of, e.g. \c constant, \c polynomial, \c table values
and the stochastic perturbation obtained from a
\c Foam::distributionModels::distributionModel.
See also
- Foam::regionModels::areaSurfaceFilmModels::contactAngleForce
- areaSurfaceFilmModels::temperatureDependentContactAngleForce
- Foam::regionModels::areaSurfaceFilmModels::distributionContactAngleForce
- Foam::Function1Types
- Foam::distributionModel
Usage
Minimal example:
\verbatim
forces
(
dynamicContactAngle
);
dynamicContactAngleForceCoeffs
{
// Mandatory entries
distribution <subDict>;
// Conditional entries
// Option-1
Utheta <Function1<scalar>>;
// Option-2
Ttheta <Function1<scalar>>;
// Inherited entries
...
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
dynamicContactAngle | Type name | word | yes | -
Utheta | Contact angle as a function of film speed <!--
--> | \<Function1\<scalar\> | choice | -
Ttheta | Contact angle as a function of film temperature <!--
--> | \<Function1\<scalar\> | choice | -
distribution | Probability distribution model | subDict | yes | -
\endtable
The inherited entries are elaborated in:
- \link contactAngleForce.H \endlink
- \link Function1.H \endlink
- \link distributionModel.H \endlink
SourceFiles
perturbedTemperatureDependentContactAngleForce.C
dynamicContactAngleForce.C
\*---------------------------------------------------------------------------*/
#ifndef perturbedTemperatureDependentContactAngleForce_H
#define perturbedTemperatureDependentContactAngleForce_H
#ifndef areaSurfaceFilmModels_dynamicContactAngleForce_H
#define areaSurfaceFilmModels_dynamicContactAngleForce_H
#include "contactAngleForce.H"
#include "Function1.H"
@ -65,17 +100,20 @@ namespace areaSurfaceFilmModels
{
/*---------------------------------------------------------------------------*\
Class perturbedTemperatureDependentContactAngleForce Declaration
Class dynamicContactAngleForce Declaration
\*---------------------------------------------------------------------------*/
class perturbedTemperatureDependentContactAngleForce
class dynamicContactAngleForce
:
public contactAngleForce
{
// Private Data
//- Contact angle function
autoPtr<Function1<scalar>> thetaPtr_;
//- Contact angle as a function of film speed
autoPtr<Function1<scalar>> U_vs_thetaPtr_;
//- Contact angle as a function of film temperature
autoPtr<Function1<scalar>> T_vs_thetaPtr_;
//- Random number generator
Random rndGen_;
@ -87,16 +125,10 @@ class perturbedTemperatureDependentContactAngleForce
// Private Member Functions
//- No copy construct
perturbedTemperatureDependentContactAngleForce
(
const perturbedTemperatureDependentContactAngleForce&
) = delete;
dynamicContactAngleForce(const dynamicContactAngleForce&) = delete;
//- No copy assignment
void operator=
(
const perturbedTemperatureDependentContactAngleForce&
) = delete;
void operator=(const dynamicContactAngleForce&) = delete;
protected:
@ -108,13 +140,13 @@ protected:
public:
//- Runtime type information
TypeName("perturbedTemperatureDependentContactAngle");
TypeName("dynamicContactAngle");
// Constructors
//- Construct from surface film model
perturbedTemperatureDependentContactAngleForce
dynamicContactAngleForce
(
liquidFilmBase& film,
const dictionary& dict
@ -122,7 +154,7 @@ public:
//- Destructor
virtual ~perturbedTemperatureDependentContactAngleForce();
virtual ~dynamicContactAngleForce() = default;
};

6
tutorials/incompressible/pimpleFoam/laminar/filmPanel0/0.orig/U
View File

@ -68,12 +68,12 @@ boundaryField
injectionModels ();
forces (perturbedTemperatureDependentContactAngle);
forces (dynamicContactAngle);
perturbedTemperatureDependentContactAngleCoeffs
dynamicContactAngleCoeffs
{
Ccf 0.4;
theta constant 0;
Ttheta constant 0;
distribution
{
type normal;

6
tutorials/incompressible/pimpleFoam/laminar/inclinedPlaneFilm/0/U
View File

@ -59,12 +59,12 @@ boundaryField
injectionModels ();
forces (perturbedTemperatureDependentContactAngle);
forces (dynamicContactAngle);
perturbedTemperatureDependentContactAngleCoeffs
dynamicContactAngleCoeffs
{
Ccf 0.085;
theta constant 45;
Ttheta constant 45;
distribution
{
type normal;

6
tutorials/lagrangian/reactingParcelFoam/splashPanelFilm/0.orig/U
View File

@ -54,12 +54,12 @@ boundaryField
injectionModels ();
forces (perturbedTemperatureDependentContactAngle);
forces (dynamicContactAngle);
perturbedTemperatureDependentContactAngleCoeffs
dynamicContactAngleCoeffs
{
Ccf 0.4;
theta constant 0;
Ttheta constant 0;
distribution
{
type normal;