diff --git a/src/regionFaModels/Make/files b/src/regionFaModels/Make/files index 4e761c227dd8a241d70a83a771b021413ca113f2..cd0717dbd83b2fc4ffd6e125939af84db33f3251 100644 --- a/src/regionFaModels/Make/files +++ b/src/regionFaModels/Make/files @@ -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 diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C index 32ff65e6f25775d4e42959688354d654f9a753c3..5dc0e4b8d10353bab81ff6e4ca3326c4db7b8360 100644 --- a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.C @@ -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(); - tmp<areaScalarField> talpha = film().alpha(); - const areaScalarField& sigma = film().sigma(); + 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& 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; } diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H index f99aab9ff2ac30e1b3133a1cd09dbab38cbdef56..d763d37e98d339b430a90520a069e65d1f06f900 100644 --- a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/contactAngleForce/contactAngleForce.H @@ -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 diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/perturbedTemperatureDependent/perturbedTemperatureDependentContactAngleForce.C b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C similarity index 61% rename from src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/perturbedTemperatureDependent/perturbedTemperatureDependentContactAngleForce.C rename to src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C index 85cf19f2c5d5ac40272fbfd8d9db29e3dd38d3ce..4119e50a984a6efe8804acc470280971ba4bec0a 100644 --- a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/perturbedTemperatureDependent/perturbedTemperatureDependentContactAngleForce.C +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.C @@ -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 diff --git a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/perturbedTemperatureDependent/perturbedTemperatureDependentContactAngleForce.H b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.H similarity index 57% rename from src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/perturbedTemperatureDependent/perturbedTemperatureDependentContactAngleForce.H rename to src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.H index ae165b90292e02837d751c9f64a5f49729250a61..8f2ec84b8813ac112b133ce4ec661acd3aa50751 100644 --- a/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/perturbedTemperatureDependent/perturbedTemperatureDependentContactAngleForce.H +++ b/src/regionFaModels/liquidFilm/subModels/kinematic/force/contactAngleForces/dynamicContactAngleForce/dynamicContactAngleForce.H @@ -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; }; diff --git a/tutorials/incompressible/pimpleFoam/laminar/filmPanel0/0.orig/U b/tutorials/incompressible/pimpleFoam/laminar/filmPanel0/0.orig/U index 820885bc3606ee4eb6d0dece8bf5c0f39ea836bb..c1d3599d0c17fdc7f3639384e3fb9b9d39b4d272 100644 --- a/tutorials/incompressible/pimpleFoam/laminar/filmPanel0/0.orig/U +++ b/tutorials/incompressible/pimpleFoam/laminar/filmPanel0/0.orig/U @@ -68,12 +68,12 @@ boundaryField injectionModels (); - forces (perturbedTemperatureDependentContactAngle); + forces (dynamicContactAngle); - perturbedTemperatureDependentContactAngleCoeffs + dynamicContactAngleCoeffs { Ccf 0.4; - theta constant 0; + Ttheta constant 0; distribution { type normal; diff --git a/tutorials/incompressible/pimpleFoam/laminar/inclinedPlaneFilm/0/U b/tutorials/incompressible/pimpleFoam/laminar/inclinedPlaneFilm/0/U index d93910cbe5d35e61d6b0523f99639363d9102be2..2c2f5c41226c1ed6217fbff04b4fb64274cb9e5a 100644 --- a/tutorials/incompressible/pimpleFoam/laminar/inclinedPlaneFilm/0/U +++ b/tutorials/incompressible/pimpleFoam/laminar/inclinedPlaneFilm/0/U @@ -59,12 +59,12 @@ boundaryField injectionModels (); - forces (perturbedTemperatureDependentContactAngle); + forces (dynamicContactAngle); - perturbedTemperatureDependentContactAngleCoeffs + dynamicContactAngleCoeffs { Ccf 0.085; - theta constant 45; + Ttheta constant 45; distribution { type normal; diff --git a/tutorials/lagrangian/reactingParcelFoam/splashPanelFilm/0.orig/U b/tutorials/lagrangian/reactingParcelFoam/splashPanelFilm/0.orig/U index 032780506ce18aae7e013e3dd10a23c633624de9..fbf0aa9b9d67770e1b24a7d69872ef308ecc374c 100644 --- a/tutorials/lagrangian/reactingParcelFoam/splashPanelFilm/0.orig/U +++ b/tutorials/lagrangian/reactingParcelFoam/splashPanelFilm/0.orig/U @@ -54,12 +54,12 @@ boundaryField injectionModels (); - forces (perturbedTemperatureDependentContactAngle); + forces (dynamicContactAngle); - perturbedTemperatureDependentContactAngleCoeffs + dynamicContactAngleCoeffs { Ccf 0.4; - theta constant 0; + Ttheta constant 0; distribution { type normal;