thixotropicViscosity.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "thixotropicViscosity.H"
27 #include "kinematicSingleLayer.H"
29 
30 #include "fvmDdt.H"
31 #include "fvmDiv.H"
32 #include "fvcDiv.H"
33 #include "fvmSup.H"
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 namespace regionModels
40 {
41 namespace surfaceFilmModels
42 {
43 
44 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
45 
46 defineTypeNameAndDebug(thixotropicViscosity, 0);
47 
49 (
50  filmViscosityModel,
51  thixotropicViscosity,
52  dictionary
53 );
54 
55 
56 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
57 
58 thixotropicViscosity::thixotropicViscosity
59 (
60  surfaceFilmModel& owner,
61  const dictionary& dict,
62  volScalarField& mu
63 )
64 :
65  filmViscosityModel(typeName, owner, dict, mu),
66  a_("a", dimless/dimTime, coeffDict_),
67  b_("b", dimless, coeffDict_),
68  d_("d", dimless, coeffDict_),
69  c_("c", pow(dimTime, d_.value() - scalar(1)), coeffDict_),
70  mu0_("mu0", dimPressure*dimTime, coeffDict_),
71  muInf_("muInf", mu0_.dimensions(), coeffDict_),
72  K_(1.0 - Foam::sqrt(muInf_/mu0_)),
73  lambda_
74  (
75  IOobject
76  (
77  typeName + ":lambda",
78  owner.regionMesh().time().timeName(),
79  owner.regionMesh(),
82  ),
83  owner.regionMesh()
84  )
85 {
86  lambda_.min(1.0);
87  lambda_.max(0.0);
88 
89  // Initialise viscosity to inf value because it cannot be evaluated yet
90  mu_ = muInf_;
91  mu_.correctBoundaryConditions();
92 }
93 
94 
95 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
96 
98 {}
99 
100 
101 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
102 
104 (
105  const volScalarField& p,
106  const volScalarField& T
107 )
108 {
109  const kinematicSingleLayer& film = filmType<kinematicSingleLayer>();
110 
111  // references to film fields
112  const volVectorField& U = film.U();
113  const volVectorField& Uw = film.Uw();
114  const volScalarField& delta = film.delta();
115  const volScalarField& deltaRho = film.deltaRho();
116  const surfaceScalarField& phi = film.phi();
117  const volScalarField& alpha = film.alpha();
118  const Time& runTime = this->owner().regionMesh().time();
119 
120  // Shear rate
121  volScalarField gDot("gDot", alpha*mag(U - Uw)/(delta + film.deltaSmall()));
122 
123  if (debug && runTime.writeTime())
124  {
125  gDot.write();
126  }
127 
128  dimensionedScalar deltaRho0("deltaRho0", deltaRho.dimensions(), ROOTVSMALL);
129  surfaceScalarField phiU(phi/fvc::interpolate(deltaRho + deltaRho0));
130 
131  dimensionedScalar c0("c0", dimless/dimTime, ROOTVSMALL);
132  volScalarField coeff("coeff", -c_*pow(gDot, d_) + c0);
133 
134  // Limit the filmMass and deltaMass to calculate the effect of the added
135  // droplets with lambda = 0 to the film
136  const volScalarField filmMass
137  (
138  "thixotropicViscosity:filmMass",
139  film.netMass() + dimensionedScalar("SMALL", dimMass, ROOTVSMALL)
140  );
141  const volScalarField deltaMass
142  (
143  "thixotropicViscosity:deltaMass",
144  max(dimensionedScalar("zero", dimMass, 0), film.deltaMass())
145  );
146 
147  fvScalarMatrix lambdaEqn
148  (
150  + fvm::div(phiU, lambda_)
151  - fvm::Sp(fvc::div(phiU), lambda_)
152  ==
153  a_*pow((1.0 - lambda_), b_)
154  + fvm::SuSp(coeff, lambda_)
155  - fvm::Sp(deltaMass/(runTime.deltaT()*filmMass), lambda_)
156  );
157 
158  lambdaEqn.relax();
159  lambdaEqn.solve();
160 
161  lambda_.min(1.0);
162  lambda_.max(0.0);
163 
164  mu_ = muInf_/(sqr(1.0 - K_*lambda_) + ROOTVSMALL);
166 }
167 
168 
169 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
170 
171 } // End namespace surfaceFilmModels
172 } // End namespace regionModels
173 } // End namespace Foam
174 
175 // ************************************************************************* //
scalar delta
surfaceScalarField & phi
U
Definition: pEqn.H:83
tmp< volScalarField > deltaMass() const
Return the change in film mass due to sources/sinks.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
tmp< fvMatrix< Type > > Sp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
dimensionedScalar sqrt(const dimensionedScalar &ds)
virtual const volVectorField & U() const
Return the film velocity [m/s].
virtual const volVectorField & Uw() const
Return the film wall velocity [m/s].
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:715
bool writeTime() const
Return true if this is a write time.
Definition: TimeStateI.H:65
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
Macros for easy insertion into run-time selection tables.
const fvMesh & regionMesh() const
Return the region mesh database.
Definition: regionModelI.H:61
addToRunTimeSelectionTable(surfaceFilmModel, kinematicSingleLayer, mesh)
tmp< fvMatrix< Type > > SuSp(const DimensionedField< scalar, volMesh > &, const GeometricField< Type, fvPatchField, volMesh > &)
const volScalarField & delta() const
Return const access to the film thickness / [m].
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
Calulate the matrix for the first temporal derivative.
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
Definition: fvPatchField.H:71
void min(const dimensioned< Type > &)
virtual void correct(const volScalarField &p, const volScalarField &T)
Correct.
const dimensionSet dimPressure
dimensionedScalar muInf_
Limiting viscosity when lambda = 0.
const dimensionSet & dimensions() const
Return dimensions.
Calculate the divergence of the given field.
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
virtual const volScalarField & deltaRho() const
Return the film thickness*density (helper field) [kg/m3].
void relax(const scalar alpha)
Relax matrix (for steady-state solution).
Definition: fvMatrix.C:511
volScalarField & mu_
Reference to the viscosity field.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Calculate the matrix for the divergence of the given field and flux.
tmp< volScalarField > netMass() const
Return the net film mass available over the next integration.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
dimensionedScalar deltaT() const
Return time step.
Definition: TimeStateI.H:53
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const volScalarField & alpha() const
Return the film coverage, 1 = covered, 0 = uncovered / [].
void max(const dimensioned< Type > &)
const surfaceFilmModel & owner() const
Return const access to the owner surface film model.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
virtual bool write() const
Write using setting from DB.
void correctBoundaryConditions()
Correct boundary field.
dimensionedScalar c_
Model `c&#39; coefficient (read after d since dims depend on d value)
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual const surfaceScalarField & phi() const
Return the film flux [kg.m/s].
const dimensionedScalar & deltaSmall() const
Return small delta.
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
volScalarField & p
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
defineTypeNameAndDebug(kinematicSingleLayer, 0)
Calculate the matrix for implicit and explicit sources.
Namespace for OpenFOAM.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243