ENH: Surface film modelling - updated thixotropic viscosity model

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

@@ -91,12 +91,12 @@ thixotropicViscosity::thixotropicViscosity
 )
 :
     filmViscosityModel(typeName, owner, dict, mu),
-    a_(coeffDict_.lookup("a")),
-    b_(coeffDict_.lookup("b")),
-    c_(coeffDict_.lookup("c")),
-    d_(coeffDict_.lookup("d")),
-    mu0_(coeffDict_.lookup("mu0")),
-    muInf_(coeffDict_.lookup("muInf")),
+    a_("a", dimless/dimTime, coeffDict_.lookup("a")),
+    b_("b", dimless, coeffDict_.lookup("b")),
+    d_("d", dimless, coeffDict_.lookup("d")),
+    c_("c", pow(dimTime, d_.value() - scalar(1)), coeffDict_.lookup("c")),
+    mu0_("mu0", dimPressure*dimTime, coeffDict_.lookup("mu0")),
+    muInf_("muInf", mu0_.dimensions(), coeffDict_.lookup("muInf")),
     K_(1.0 - Foam::sqrt(muInf_/mu0_)),
     lambda_
     (
@@ -144,27 +144,33 @@ void thixotropicViscosity::correct
     const volScalarField& delta = film.delta();
     const volScalarField& deltaRho = film.deltaRho();
     const surfaceScalarField& phi = film.phi();
+    const volScalarField& alpha = film.alpha();
 
-    // gamma-dot (shear rate) raised to the power d
-    volScalarField gDotPowD
-    (
-        "gDotPowD",
-        pow(mag(U - Uw)/(delta + film.deltaSmall()), d_)
-    );
+    // gamma-dot (shear rate)
+    volScalarField gDot("gDot", alpha*mag(U - Uw)/(delta + film.deltaSmall()));
+
+    if (debug && this->owner().regionMesh().time().outputTime())
+    {
+        gDot.write();
+    }
 
-    dimensionedScalar c0("SMALL", dimMass/sqr(dimLength)/dimTime, SMALL);
-    volScalarField coeff(-deltaRho*c_*gDotPowD + c0);
+    dimensionedScalar deltaRho0("deltaRho0", deltaRho.dimensions(), ROOTVSMALL);
+    surfaceScalarField phiU(phi/fvc::interpolate(deltaRho + deltaRho0));
+
+    dimensionedScalar c0("c0", dimless/dimTime, ROOTVSMALL);
+    volScalarField coeff("coeff", -c_*pow(gDot, d_) + c0);
 
     fvScalarMatrix lambdaEqn
     (
-        fvm::ddt(deltaRho, lambda_)
-      + fvm::div(phi, lambda_)
-      - fvm::Sp(fvc::div(phi), lambda_)
+        fvm::ddt(lambda_)
+      + fvm::div(phiU, lambda_)
+      - fvm::Sp(fvc::div(phiU), lambda_)
       ==
-        deltaRho*a_*pow((1.0 - lambda_), b_)
+        a_*pow((1.0 - lambda_), b_)
       + fvm::SuSp(coeff, lambda_)
     );
 
+
     lambdaEqn.relax();
 
     lambdaEqn.solve();

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,38 +105,18 @@ 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_;
 
+        //- Model `c' coefficient (read after d since dims depend on d value)
+        dimensionedScalar c_;
+
         //- Limiting viscosity when lambda = 1
         dimensionedScalar mu0_;