Commit 23d58e54 authored by andy's avatar andy
Browse files

ENH: Updated and restructured film force models, new contact angle force

parent 3d3cfa91
......@@ -10,6 +10,13 @@ thermoSingleLayer/thermoSingleLayer.C
submodels/subModelBase.C
KINEMATICMODELS=submodels/kinematic
$(KINEMATICMODELS)/force/force/force.C
$(KINEMATICMODELS)/force/force/forceNew.C
$(KINEMATICMODELS)/force/forceList/forceList.C
$(KINEMATICMODELS)/force/contactAngleForce/contactAngleForce.C
$(KINEMATICMODELS)/force/surfaceShearForce/surfaceShearForce.C
$(KINEMATICMODELS)/force/thermocapillaryForce/thermocapillaryForce.C
$(KINEMATICMODELS)/injectionModel/injectionModel/injectionModel.C
$(KINEMATICMODELS)/injectionModel/injectionModel/injectionModelNew.C
$(KINEMATICMODELS)/injectionModel/injectionModelList/injectionModelList.C
......
......@@ -34,6 +34,10 @@ License
#include "directMappedWallPolyPatch.H"
#include "mapDistribute.H"
#include "cachedRandom.H"
#include "normal.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
......@@ -61,8 +65,6 @@ bool kinematicSingleLayer::read()
solution.lookup("nCorr") >> nCorr_;
solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
coeffs_.lookup("Cf") >> Cf_;
return true;
}
else
......@@ -76,9 +78,10 @@ void kinematicSingleLayer::correctThermoFields()
{
if (thermoModel_ == tmConstant)
{
rho_ == dimensionedScalar(coeffs_.lookup("rho0"));
mu_ == dimensionedScalar(coeffs_.lookup("mu0"));
sigma_ == dimensionedScalar(coeffs_.lookup("sigma0"));
const dictionary& constDict(coeffs_.subDict("constantThermoCoeffs"));
rho_ == dimensionedScalar(constDict.lookup("rho0"));
mu_ == dimensionedScalar(constDict.lookup("mu0"));
sigma_ == dimensionedScalar(constDict.lookup("sigma0"));
}
else
{
......@@ -273,25 +276,6 @@ void kinematicSingleLayer::updateSurfaceVelocities()
}
tmp<fvVectorMatrix> kinematicSingleLayer::tau(volVectorField& U) const
{
// Calculate shear stress
volScalarField Cs("Cs", rho_*Cf_*mag(Us_ - U));
volScalarField Cw
(
"Cw",
mu_/(0.3333*(delta_ + dimensionedScalar("SMALL", dimLength, SMALL)))
);
Cw.min(1.0e+06);
return
(
- fvm::Sp(Cs, U) + Cs*Us_ // surface contribution
- fvm::Sp(Cw, U) + Cw*Uw_ // wall contribution
);
}
tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
(
const volScalarField& pu,
......@@ -312,9 +296,8 @@ tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
+ fvm::div(phi_, U_)
==
- USp_
+ tau(U_)
+ fvc::grad(sigma_)
- fvm::SuSp(rhoSp_, U_)
+ forces_.correct(U_)
);
fvVectorMatrix& UEqn = tUEqn();
......@@ -459,8 +442,6 @@ kinematicSingleLayer::kinematicSingleLayer
cumulativeContErr_(0.0),
Cf_(readScalar(coeffs().lookup("Cf"))),
rho_
(
IOobject
......@@ -773,6 +754,8 @@ kinematicSingleLayer::kinematicSingleLayer
injection_(*this, coeffs_),
forces_(*this, coeffs_),
addedMassTotal_(0.0)
{
if (readFields)
......
......@@ -42,6 +42,7 @@ SourceFiles
#include "fvMatrices.H"
#include "injectionModelList.H"
#include "forceList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -93,12 +94,6 @@ protected:
scalar cumulativeContErr_;
// Model parameters
//- Skin frition coefficient for film/primary region interface
scalar Cf_;
// Thermo properties
// Fields
......@@ -199,6 +194,9 @@ protected:
//- Cloud injection
injectionModelList injection_;
//- List of film forces
forceList forces_;
// Checks
......@@ -238,9 +236,6 @@ protected:
//- Update film surface velocities
virtual void updateSurfaceVelocities();
//- Return the stress term for the momentum equation
virtual tmp<fvVectorMatrix> tau(volVectorField& dU) const;
//- Constrain a film region master/slave boundaries of a field to a
// given value
template<class Type>
......@@ -314,12 +309,6 @@ public:
inline label nNonOrthCorr() const;
// Model parameters
//- Return the skin friction coefficient
inline scalar Cf() const;
// Thermo properties
//- Return const access to the dynamic viscosity / [Pa.s]
......
......@@ -61,12 +61,6 @@ inline label kinematicSingleLayer::nNonOrthCorr() const
}
inline scalar kinematicSingleLayer::Cf() const
{
return Cf_;
}
inline const volScalarField& kinematicSingleLayer::mu() const
{
return mu_;
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "contactAngleForce.H"
#include "addToRunTimeSelectionTable.H"
#include "fvcGrad.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace regionModels
{
namespace surfaceFilmModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(contactAngleForce, 0);
addToRunTimeSelectionTable(force, contactAngleForce, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
contactAngleForce::contactAngleForce
(
const surfaceFilmModel& owner,
const dictionary& dict
)
:
force(typeName, owner, dict),
deltaWet_(readScalar(coeffs_.lookup("deltaWet"))),
Ccf_(readScalar(coeffs_.lookup("Ccf"))),
rndGen_(label(0), -1),
distribution_
(
distributionModels::distributionModel::New
(
coeffs_.subDict("contactAngleDistribution"),
rndGen_
)
)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
contactAngleForce::~contactAngleForce()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
tmp<fvVectorMatrix> contactAngleForce::correct(volVectorField& U)
{
tmp<volVectorField> tForce
(
new volVectorField
(
IOobject
(
"contactForce",
owner_.time().timeName(),
owner_.regionMesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
owner_.regionMesh(),
dimensionedVector("zero", dimForce/dimArea, vector::zero)
)
);
vectorField& force = tForce().internalField();
const labelUList& own = owner_.regionMesh().owner();
const labelUList& nbr = owner_.regionMesh().neighbour();
const scalarField& magSf = owner_.magSf();
const volScalarField& delta = owner_.delta();
const volScalarField& sigma = owner_.sigma();
volScalarField alpha
(
"alpha",
pos(delta - dimensionedScalar("deltaWet", dimLength, deltaWet_))
);
volVectorField gradAlpha(fvc::grad(alpha));
scalarField nHits(force.size(), 0.0);
forAll(nbr, faceI)
{
const label cellO = own[faceI];
const label cellN = nbr[faceI];
label cellI = -1;
if ((delta[cellO] > deltaWet_) && (delta[cellN] < deltaWet_))
{
cellI = cellO;
}
else if ((delta[cellO] < deltaWet_) && (delta[cellN] > deltaWet_))
{
cellI = cellN;
}
if (cellI != -1)
{
// const scalar dx = Foam::sqrt(magSf[cellI]);
const scalar dx = owner_.regionMesh().deltaCoeffs()[faceI];
const vector n =
gradAlpha[cellI]/(mag(gradAlpha[cellI]) + ROOTVSMALL);
scalar theta = cos(degToRad(distribution_->sample()));
force[cellI] += Ccf_*n*sigma[cellI]*(1.0 - theta)/dx;
nHits[cellI]++;
}
}
nHits = max(nHits, 1.0);
force /= (nHits*magSf);
if (owner_.regionMesh().time().outputTime())
{
tForce().write();
}
tmp<fvVectorMatrix>
tfvm(new fvVectorMatrix(U, dimForce/dimArea*dimVolume));
tfvm() += tForce;
return tfvm;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace surfaceFilmModels
} // End namespace regionModels
} // End namespace Foam
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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::contactAngleForce
Description
Film contact angle force
SourceFiles
contactAngleForce.C
\*---------------------------------------------------------------------------*/
#ifndef contactAngleForce_H
#define contactAngleForce_H
#include "force.H"
#include "distributionModel.H"
#include "cachedRandom.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace regionModels
{
namespace surfaceFilmModels
{
/*---------------------------------------------------------------------------*\
Class contactAngleForce Declaration
\*---------------------------------------------------------------------------*/
class contactAngleForce
:
public force
{
private:
// Private Data
//- Threshold film thickness beyon which the film is 'wet'
scalar deltaWet_;
//- Coefficient applied to the contact angle force
scalar Ccf_;
//- Random number generator
cachedRandom rndGen_;
//- Parcel size PDF model
const autoPtr<distributionModels::distributionModel> distribution_;
// Private member functions
//- Disallow default bitwise copy construct
contactAngleForce(const contactAngleForce&);
//- Disallow default bitwise assignment
void operator=(const contactAngleForce&);
public:
//- Runtime type information
TypeName("contactAngle");
// Constructors
//- Construct from surface film model
contactAngleForce
(
const surfaceFilmModel& owner,
const dictionary& dict
);
//- Destructor
virtual ~contactAngleForce();
// Member Functions
// Evolution
//- Correct
virtual tmp<fvVectorMatrix> correct(volVectorField& U);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace surfaceFilmModels
} // End namespace regionModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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 "force.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace regionModels
{
namespace surfaceFilmModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(force, 0);
defineRunTimeSelectionTable(force, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
force::force(const surfaceFilmModel& owner)
:
subModelBase(owner)
{}
force::force
(
const word& type,
const surfaceFilmModel& owner,
const dictionary& dict
)
:
subModelBase(type, owner, dict)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
force::~force()
{}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace surfaceFilmModels
} // End namespace regionModels
} // End namespace Foam
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
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::force
Description
Base class for film (stress-based) force models
SourceFiles
force.C
forceNew.C
\*---------------------------------------------------------------------------*/
#ifndef force_H
#define force_H
#include "subModelBase.H"
#include "runTimeSelectionTables.H"
#include "fvMatrices.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace regionModels
{
namespace surfaceFilmModels
{
/*---------------------------------------------------------------------------*\
Class force Declaration
\*---------------------------------------------------------------------------*/
class force
:
public subModelBase
{
private: