momentumTransportModels::lambdaThixotropic: New momentum transport model for thixotropic liquids

Description
    Thixotropic viscosity momentum transport model based on the evolution of
    the structural parameter \f$ \lambda \f$:

        \f[
            \frac{D\lambda}{Dt} = a(1 - \lambda)^b - c \lambda \dot{\gamma}^d
        \f]

    The viscosity is then calculated using the expression

        \f[
            \nu = \frac{\nu_{\infty}}{{1 - K \lambda}^2}
        \f]

    Where the parameter K is given by:

        \f[
            K = 1 - \sqrt{\frac{\nu_{\infty}}{\nu_{0}}}
        \f]

    Here:
    \vartable
        \lambda         | structural parameter
        a               | model coefficient
        b               | model coefficient
        c               | model coefficient
        d               | model coefficient
        \dot{\gamma}    | stress rate [1/s]
        \nu_{0}         | limiting viscosity when \f$ \lambda = 1 \f$
        \nu_{\infty}    | 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

applications/solvers/multiphase/compressibleInterFoam/VoFphaseCompressibleMomentumTransportModels/VoFphaseCompressibleMomentumTransportModels.C

@@ -67,6 +67,12 @@ makeBaseMomentumTransportModel
 #include "Stokes.H"
 makeLaminarModel(Stokes);
 
+#include "generalizedNewtonian.H"
+makeLaminarModel(generalizedNewtonian);
+
+#include "lambdaThixotropic.H"
+makeLaminarModel(lambdaThixotropic);
+
 #include "Maxwell.H"
 makeLaminarModel(Maxwell);
 

src/MomentumTransportModels/compressible/fluidThermoMomentumTransportModels/fluidThermoMomentumTransportModels.C

@@ -47,6 +47,9 @@ makeLaminarModel(Stokes);
 #include "generalizedNewtonian.H"
 makeLaminarModel(generalizedNewtonian);
 
+#include "lambdaThixotropic.H"
+makeLaminarModel(lambdaThixotropic);
+
 #include "Maxwell.H"
 makeLaminarModel(Maxwell);
 

src/MomentumTransportModels/incompressible/kinematicMomentumTransportModels/kinematicMomentumTransportModels.C

@@ -47,6 +47,9 @@ makeLaminarModel(Stokes);
 #include "generalizedNewtonian.H"
 makeLaminarModel(generalizedNewtonian);
 
+#include "lambdaThixotropic.H"
+makeLaminarModel(lambdaThixotropic);
+
 #include "Maxwell.H"
 makeLaminarModel(Maxwell);
 

src/MomentumTransportModels/momentumTransportModels/laminar/lambdaThixotropic/lambdaThixotropic.C

@@ -0,0 +1,213 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2020 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 "lambdaThixotropic.H"
+#include "fvOptions.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicMomentumTransportModel>
+lambdaThixotropic<BasicMomentumTransportModel>::lambdaThixotropic
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport
+)
+:
+    linearViscousStress<laminarModel<BasicMomentumTransportModel>>
+    (
+        typeName,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport
+    ),
+
+    a_("a", dimless/dimTime, this->coeffDict_),
+    b_("b", dimless, this->coeffDict_),
+    d_("d", dimless, this->coeffDict_),
+    c_("c", pow(dimTime, d_.value() - scalar(1)), this->coeffDict_),
+    nu0_("nu0", dimViscosity, this->coeffDict_),
+    nuInf_("nuInf", dimViscosity, this->coeffDict_),
+    K_(1 - sqrt(nuInf_/nu0_)),
+
+    lambda_
+    (
+        IOobject
+        (
+            IOobject::groupName
+            (
+                IOobject::modelName("lambda", typeName),
+                alphaRhoPhi.group()
+            ),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+
+    nu_
+    (
+        IOobject
+        (
+            IOobject::groupName
+            (
+                IOobject::modelName("nu", typeName),
+                alphaRhoPhi.group()
+            ),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        calcNu()
+    )
+{}
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class BasicMomentumTransportModel>
+tmp<volScalarField>
+lambdaThixotropic<BasicMomentumTransportModel>::calcNu() const
+{
+    return nuInf_/(sqr(1 - K_*lambda_) + rootVSmall);
+}
+
+
+template<class BasicMomentumTransportModel>
+tmp<volScalarField::Internal>
+lambdaThixotropic<BasicMomentumTransportModel>::strainRate() const
+{
+    return sqrt(2.0)*mag(symm(fvc::grad(this->U())()()));
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicMomentumTransportModel>
+bool lambdaThixotropic<BasicMomentumTransportModel>::read()
+{
+    if (laminarModel<BasicMomentumTransportModel>::read())
+    {
+        a_.read(this->coeffDict());
+        b_.read(this->coeffDict());
+        d_.read(this->coeffDict());
+
+        c_ = dimensionedScalar
+        (
+            "c",
+            pow(dimTime, d_.value() - scalar(1)),
+            this->coeffDict_
+        );
+
+        nu0_.read(this->coeffDict());
+        nuInf_.read(this->coeffDict());
+
+        K_ = (1 - sqrt(nuInf_/nu0_));
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicMomentumTransportModel>
+tmp<volScalarField>
+lambdaThixotropic<BasicMomentumTransportModel>::nuEff() const
+{
+    return volScalarField::New
+    (
+        IOobject::groupName("nuEff", this->alphaRhoPhi_.group()),
+        nu_
+    );
+}
+
+
+template<class BasicMomentumTransportModel>
+tmp<scalarField>
+lambdaThixotropic<BasicMomentumTransportModel>::nuEff
+(
+    const label patchi
+) const
+{
+    return nu_.boundaryField()[patchi];
+}
+
+
+template<class BasicMomentumTransportModel>
+void lambdaThixotropic<BasicMomentumTransportModel>::correct()
+{
+    // Local references
+    const surfaceScalarField& phi = this->phi_;
+    const fv::options& fvOptions(fv::options::New(this->mesh_));
+
+    tmp<fvScalarMatrix> lambdaEqn
+    (
+        fvm::ddt(lambda_) + fvm::div(phi, lambda_)
+      - fvm::Sp(fvc::div(phi), lambda_)
+     ==
+        a_*pow(1 - lambda_(), b_)
+      - fvm::Sp(c_*pow(strainRate(), d_), lambda_)
+      + fvOptions(lambda_)
+    );
+
+    lambdaEqn.ref().relax();
+    fvOptions.constrain(lambdaEqn.ref());
+    solve(lambdaEqn);
+    fvOptions.correct(lambda_);
+
+    lambda_.maxMin(scalar(0), scalar(1));
+
+    nu_ = calcNu();
+
+    laminarModel<BasicMomentumTransportModel>::correct();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace laminarModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/MomentumTransportModels/momentumTransportModels/laminar/lambdaThixotropic/lambdaThixotropic.H

@@ -0,0 +1,195 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2020 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::laminarModels::lambdaThixotropic
+
+Description
+    Thixotropic viscosity momentum transport model based on the evolution of
+    the structural parameter \f$ \lambda \f$:
+
+        \f[
+            \frac{D\lambda}{Dt} = a(1 - \lambda)^b - c \lambda \dot{\gamma}^d
+        \f]
+
+    The viscosity is then calculated using the expression
+
+        \f[
+            \nu = \frac{\nu_{\infty}}{{1 - K \lambda}^2}
+        \f]
+
+    Where the parameter K is given by:
+
+        \f[
+            K = 1 - \sqrt{\frac{\nu_{\infty}}{\nu_{0}}}
+        \f]
+
+    Here:
+    \vartable
+        \lambda         | structural parameter
+        a               | model coefficient
+        b               | model coefficient
+        c               | model coefficient
+        d               | model coefficient
+        \dot{\gamma}    | stress rate [1/s]
+        \nu_{0}         | limiting viscosity when \f$ \lambda = 1 \f$
+        \nu_{\infty}    | 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
+    lambdaThixotropic.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef lambdaThixotropic_H
+#define lambdaThixotropic_H
+
+#include "laminarModel.H"
+#include "linearViscousStress.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace laminarModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class lambdaThixotropic Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicMomentumTransportModel>
+class lambdaThixotropic
+:
+    public linearViscousStress<laminarModel<BasicMomentumTransportModel>>
+{
+    // Private data
+
+        //- Model a coefficient
+        dimensionedScalar a_;
+
+        //- Model b coefficient
+        dimensionedScalar b_;
+
+        //- Model d coefficient
+        dimensionedScalar d_;
+
+        //- Model c coefficient (read after d since dims depend on d value)
+        dimensionedScalar c_;
+
+        //- Limiting viscosity when lambda = 1
+        dimensionedScalar nu0_;
+
+        //- Limiting viscosity when lambda = 0
+        dimensionedScalar nuInf_;
+
+        //- Model coefficient
+        dimensionedScalar K_;
+
+        //- Structural parameter
+        //  0 = freestream value (most liquid)
+        //  1 = fully built (most solid)
+        volScalarField lambda_;
+
+        //- The non-Newtonian viscosity field
+        volScalarField nu_;
+
+
+    // Private Member Functions
+
+        //- Calculates and returns the viscosity from the current lambda
+        tmp<volScalarField> calcNu() const;
+
+        //- Returns the current strain rate from the velocity field
+        tmp<volScalarField::Internal> strainRate() const;
+
+
+public:
+
+    typedef typename BasicMomentumTransportModel::alphaField alphaField;
+    typedef typename BasicMomentumTransportModel::rhoField rhoField;
+    typedef typename BasicMomentumTransportModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("lambdaThixotropic");
+
+
+    // Constructors
+
+        //- Construct from components
+        lambdaThixotropic
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport
+        );
+
+
+    //- Destructor
+    virtual ~lambdaThixotropic()
+    {}
+
+
+    // Member Functions
+
+        //- Read momentumTransport dictionary
+        virtual bool read();
+
+        //- Return the effective viscosity
+        //  i.e. the lambdaThixotropic viscosity
+        virtual tmp<volScalarField> nuEff() const;
+
+        //- Return the effective viscosity on patch
+        virtual tmp<scalarField> nuEff(const label patchi) const;
+
+        //- Correct the lambdaThixotropic viscosity
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace laminarModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "lambdaThixotropic.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //