ENH: Surface film - added new thixotropic viscosity model

src/regionModels/surfaceFilmModels/Make/files

@@ -24,8 +24,9 @@ $(KINEMATICMODELS)/injectionModel/drippingInjection/drippingInjection.C
 $(KINEMATICMODELS)/injectionModel/removeInjection/removeInjection.C
 $(KINEMATICMODELS)/injectionModel/curvatureSeparation/curvatureSeparation.C
 
-$(KINEMATICMODELS)/filmThermoModel/constantFilmThermo/constantFilmThermo.C
 $(KINEMATICMODELS)/filmThermoModel/filmThermoModel/filmThermoModel.C
+$(KINEMATICMODELS)/filmThermoModel/filmThermoModel/filmThermoModelNew.C
+$(KINEMATICMODELS)/filmThermoModel/constantFilmThermo/constantFilmThermo.C
 $(KINEMATICMODELS)/filmThermoModel/liquidFilmThermo/liquidFilmThermo.C
 
 $(KINEMATICMODELS)/filmTurbulenceModel/filmTurbulenceModel/filmTurbulenceModel.C
@@ -50,6 +51,12 @@ $(THERMOMODELS)/filmRadiationModel/constantRadiation/constantRadiation.C
 $(THERMOMODELS)/filmRadiationModel/primaryRadiation/primaryRadiation.C
 $(THERMOMODELS)/filmRadiationModel/standardRadiation/standardRadiation.C
 
+$(THERMOMODELS)/filmViscosityModel/filmViscosityModel/filmViscosityModel.C
+$(THERMOMODELS)/filmViscosityModel/filmViscosityModel/filmViscosityModelNew.C
+$(THERMOMODELS)/filmViscosityModel/constantViscosity/constantViscosity.C
+$(THERMOMODELS)/filmViscosityModel/liquidViscosity/liquidViscosity.C
+$(THERMOMODELS)/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
+
 
 /* Boundary conditions */
 PATCHFIELDS=derivedFvPatchFields

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

@@ -438,6 +438,8 @@ kinematicSingleLayer::kinematicSingleLayer
 
     cumulativeContErr_(0.0),
 
+    deltaSmall_("deltaSmall", dimLength, SMALL),
+
     rho_
     (
         IOobject

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.H

@@ -97,6 +97,9 @@ protected:
             //- Cumulative continuity error
             scalar cumulativeContErr_;
 
+            //- Small delta
+            const dimensionedScalar deltaSmall_;
+
 
         // Thermo properties
 
@@ -324,6 +327,9 @@ public:
             //- Return the number of non-orthogonal correctors
             inline label nNonOrthCorr() const;
 
+            //- Return small delta
+            inline const dimensionedScalar& deltaSmall() const;
+
 
         // Thermo properties
 

src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayerI.H

@@ -62,6 +62,12 @@ inline label kinematicSingleLayer::nNonOrthCorr() const
 }
 
 
+inline const dimensionedScalar& kinematicSingleLayer::deltaSmall() const
+{
+    return deltaSmall_;
+}
+
+
 inline const volScalarField& kinematicSingleLayer::mu() const
 {
     return mu_;

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C

@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 "thixotropicViscosity.H"
+#include "kinematicSingleLayer.H"
+#include "addToRunTimeSelectionTable.H"
+
+#include "fvmDdt.H"
+#include "fvmDiv.H"
+#include "fvcDiv.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(thixotropicViscosity, 0);
+
+addToRunTimeSelectionTable
+(
+    filmViscosityModel,
+    thixotropicViscosity,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void thixotropicViscosity::updateMu()
+{
+    const kinematicSingleLayer& film = filmType<kinematicSingleLayer>();
+
+    // blend based on mass fraction of added- to existing film mass
+    const dimensionedScalar m0("zero", dimMass, 0.0);
+    const dimensionedScalar mSMALL("SMALL", dimMass, ROOTVSMALL);
+    const volScalarField deltaMass("deltaMass", max(m0, film.deltaMass()));
+    const volScalarField filmMass("filmMass", film.netMass() + mSMALL);
+
+    // weighting field to blend new and existing mass contributions
+    const volScalarField w("w", max(0.0, min(1.0, deltaMass/filmMass)));
+
+    // evaluate thixotropic viscosity
+    volScalarField muThx("muThx", muInf_/(sqr(1.0 - K_*lambda_) + ROOTVSMALL));
+
+    // set new viscosity based on weight field
+    mu_ = w*muInf_ + (1.0 - w)*muThx;
+    mu_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+thixotropicViscosity::thixotropicViscosity
+(
+    const surfaceFilmModel& owner,
+    const dictionary& dict,
+    volScalarField& mu
+)
+:
+    filmViscosityModel(typeName, owner, dict, mu),
+    a_(coeffs().lookup("a")),
+    b_(coeffs().lookup("b")),
+    c_(coeffs().lookup("c")),
+    d_(coeffs().lookup("d")),
+    mu0_(coeffs().lookup("mu0")),
+    muInf_(coeffs().lookup("muInf")),
+    K_(1.0 - Foam::sqrt(muInf_/mu0_)),
+    lambda_
+    (
+        IOobject
+        (
+            typeName + ":lambda",
+            owner.regionMesh().time().timeName(),
+            owner.regionMesh(),
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::AUTO_WRITE
+        ),
+        owner.regionMesh()
+    )
+{
+    lambda_.min(1.0);
+    lambda_.max(0.0);
+
+    // initialise viscosity to inf value
+    // - cannot call updateMu() since this calls film.netMass() which
+    //   cannot be evaluated yet (still in construction)
+    mu_ = muInf_;
+    mu_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+thixotropicViscosity::~thixotropicViscosity()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void thixotropicViscosity::correct
+(
+    const volScalarField& p,
+    const volScalarField& T
+)
+{
+    const kinematicSingleLayer& film = filmType<kinematicSingleLayer>();
+
+    // references to film fields
+    const volVectorField& U = film.U();
+    const volVectorField& Uw = film.Uw();
+    const volScalarField& delta = film.delta();
+    const volScalarField& deltaRho = film.deltaRho();
+    const surfaceScalarField& phi = film.phi();
+
+    // gamma-dot (shear rate) raised to the power d
+    volScalarField gDotPowD
+    (
+        "gDotPowD",
+        pow(mag(U - Uw)/(delta + film.deltaSmall()), d_)
+    );
+
+    dimensionedScalar c0("SMALL", dimMass/sqr(dimLength)/dimTime, SMALL);
+    volScalarField coeff(-deltaRho*c_*gDotPowD + c0);
+
+    fvScalarMatrix lambdaEqn
+    (
+        fvm::ddt(deltaRho, lambda_)
+      + fvm::div(phi, lambda_)
+      - fvm::Sp(fvc::div(phi), lambda_)
+      ==
+        deltaRho*a_*pow((1.0 - lambda_), b_)
+      + fvm::SuSp(coeff, lambda_)
+    );
+
+    lambdaEqn.relax();
+
+    lambdaEqn.solve();
+
+    lambda_.min(1.0);
+    lambda_.max(0.0);
+
+    updateMu();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.H

@@ -0,0 +1,205 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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::thixotropicViscosity
+
+Description
+    Thixotropic viscosity model based on the evolution of the structural
+    parameter \f$ \lambda \f$:
+
+        \f[
+            \lambda = a(1 - \lambda)^b - c \lambda \dot{\gamma}^d
+        \f]
+
+    The viscosity is then calculated using the expression
+
+        \f[
+            \mu = \frac{\mu_{inf}}{{1 - K \lambda}^2}
+        \f]
+
+    Where the parameter K is given by:
+
+        \f[
+            K = 1 - \frac{\mu_{\inf}}{\mu_{0}}^0.5
+        \f]
+
+    Here:
+    \vartable
+        \lambda         | structural parameter
+        a               | model coefficient
+        b               | model coefficient
+        c               | model coefficient
+        d               | model coefficient
+        \dot{\gamma}    | stress rate [1/s]
+        \mu_{0}         | limiting viscosity when \f$ \lambda = 1 \f$
+        \mu_{\inf}      | limiting viscosity when \f$ \lambda = 0 \f$
+    \endvartable
+
+
+    Reference:
+    \verbatim
+        Barnes H A, 1997.  Thixotropy - a review.  J. Non-Newtonian Fluid
+        Mech 70, pp 1-33
+    \endverbatim
+
+SourceFiles
+    thixotropicViscosity.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef thixotropicViscosity_H
+#define thixotropicViscosity_H
+
+#include "filmViscosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class thixotropicViscosity Declaration
+\*---------------------------------------------------------------------------*/
+
+class thixotropicViscosity
+:
+    public filmViscosityModel
+{
+private:
+
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        thixotropicViscosity(const thixotropicViscosity&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const thixotropicViscosity&);
+
+
+protected:
+
+    // Protected data
+
+/*
+        //- Model `a' coefficient
+        scalar a_;
+
+        //- Model `b' coefficient
+        scalar b_;
+
+        //- Model `c' coefficient
+        scalar c_;
+
+        //- Model `d' coefficient
+        scalar d_;
+
+        //- Limiting viscosity when lambda = 1
+        scalar mu0_;
+
+        //- Limiting viscosity when lambda = 0
+        scalar muInf_;
+*/
+
+        //- Model `a' coefficient
+        dimensionedScalar a_;
+
+        //- Model `b' coefficient
+        dimensionedScalar b_;
+
+        //- Model `c' coefficient
+        dimensionedScalar c_;
+
+        //- Model `d' coefficient
+        dimensionedScalar d_;
+
+        //- Limiting viscosity when lambda = 1
+        dimensionedScalar mu0_;
+
+        //- Limiting viscosity when lambda = 0
+        dimensionedScalar muInf_;
+
+        //- Model coeffiicient
+        dimensionedScalar K_;
+
+        //- Structural parameter
+        //  0 = freestream value (most liquid)
+        //  1 = fully built (most solid)
+        volScalarField lambda_;
+
+
+    // Protected Member Functions
+
+        //- Update the viscosity
+        void updateMu();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("thixotropic");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        thixotropicViscosity
+        (
+            const surfaceFilmModel& owner,
+            const dictionary& dict,
+            volScalarField& mu
+        );
+
+
+    //- Destructor
+    virtual ~thixotropicViscosity();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct
+            virtual void correct
+            (
+                const volScalarField& p,
+                const volScalarField& T
+            );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //