[anuargimenez]

Hi all,

I am developing my own shear and temperature dependent viscosity model in OpenFoam for a compressible solver (rhoSimpleFoam). For that, I want to modify the WLF transport model (located in src/thermophysicalModels/specie/transport) to add the Cross Power Law equation.

My code is shown below:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2018 OpenFOAM Foundation
-------------------------------------------------------------------------------
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 "specie.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

template <class Thermo>
inline Foam::polymerKCrossWLFTransport<Thermo>::polymerKCrossWLFTransport(
    const word &name,
    const polymerKCrossWLFTransport &wlft)
    : Thermo(name, wlft),
      mu0_(wlft.mu0_),
      Tr_(wlft.Tr_),
      C1_(wlft.C1_),
      C2_(wlft.C2_),
      rPr_(wlft.rPr_),
      Tg_(wlft.Tg_),
      TS_(wlft.TS_),
      MS_(wlft.MS_),
      tauStar_(wlft.tauStar_),
      n_(wlft.n_)
{
}

template <class Thermo>
inline Foam::autoPtr<Foam::polymerKCrossWLFTransport<Thermo>>
Foam::polymerKCrossWLFTransport<Thermo>::clone() const
{
    return autoPtr<polymerKCrossWLFTransport<Thermo>>(
        new polymerKCrossWLFTransport<Thermo>(*this));
}

template <class Thermo>
inline Foam::autoPtr<Foam::polymerKCrossWLFTransport<Thermo>>
Foam::polymerKCrossWLFTransport<Thermo>::New(
    const dictionary &dict)
{
    return autoPtr<polymerKCrossWLFTransport<Thermo>>(
        new polymerKCrossWLFTransport<Thermo>(dict));
}

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

template <class Thermo>
inline Foam::scalar Foam::polymerKCrossWLFTransport<Thermo>::mu(
    const scalar p,
    const scalar T) const
{

HERE I NEED TO CALL STRAINRATE FOR THE CALCULATION OF mu


    
//CROSS WLF EQUATION:
       return (mu0_*exp(-C1_*(T - Tr_)/(C2_ + T - Tr_))) / (1 + pow((mu0_*exp(-C1_*(T - Tr_)/(C2_ + T - Tr_)))*strainRate/tauStar_,1-n_));    

    
   // return mu0_*exp(-C1_*(T - Tr_)/(C2_ + T - Tr_));

    /* WLF with Melting Fraction (MF):
    scalar MF;
    MF=pow((tanh((T - Tg_) * TS_) + 1) / 2, MS_);
    
    return mu0_*exp(((1-MF)*(-C1_*(Tg_ - Tr_)/(C2_ + Tg_ - Tr_))) + (MF*(-C1_*(T - Tr_)/(C2_ + T - Tr_)))); */
}

template <class Thermo>
inline Foam::scalar Foam::polymerKCrossWLFTransport<Thermo>::kappa(
    const scalar p, const scalar T) const
{
// original version:
// return this->Cp(p, T)*mu(p, T)*rPr_;

// new version:
#include "polymerKWLF_Data_ThermalConduct.H"
}

template <class Thermo>
inline Foam::scalar Foam::polymerKCrossWLFTransport<Thermo>::alphah(
    const scalar p,
    const scalar T) const
{
    // ORIGINAL VERSION Pr constant:
    // return mu(p, T)*rPr_;

    // NEW VERSION FOR POLYMER Pr not constant:
    // Pr = mu(p,T)*Cp(p,T)/kappa(p,T)
    // alpha = kappa(p,T)/cp(p,T)*rho(p,T) = kappa(p,T)/Cp(p,T)

    return kappa(p, T) / this->Cp(p, T);

    // end of alpha version implemented for polymer
}

// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //

template <class Thermo>
inline void Foam::polymerKCrossWLFTransport<Thermo>::operator+=(
    const polymerKCrossWLFTransport<Thermo> &wlft)
{
    scalar Y1 = this->Y();

    Thermo::operator+=(wlft);

    if (mag(this->Y()) > SMALL)
    {
        Y1 /= this->Y();
        scalar Y2 = wlft.Y() / this->Y();

        mu0_ = Y1 * mu0_ + Y2 * wlft.mu0_;
        Tr_ = Y1 * Tr_ + Y2 * wlft.Tr_;
        C1_ = Y1 * C1_ + Y2 * wlft.C1_;
        C2_ = Y1 * C2_ + Y2 * wlft.C2_;
        rPr_ = 1.0 / (Y1 / rPr_ + Y2 / wlft.rPr_);
        Tg_ = Y1 * Tg_ + Y2 * wlft.Tg_;
        TS_ = Y1 * TS_ + Y2 * wlft.TS_;
        MS_ = Y1 * MS_ + Y2 * wlft.MS_;
        tauStar_ = Y1 * tauStar_ + Y2 * wlft.tauStar_;
        n_ = Y1 * n_ + Y2 * wlft.n_;        
    }
}

template <class Thermo>
inline void Foam::polymerKCrossWLFTransport<Thermo>::operator*=(
    const scalar s)
{
    Thermo::operator*=(s);
}

// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //

template <class Thermo>
inline Foam::polymerKCrossWLFTransport<Thermo> Foam::operator+(
    const polymerKCrossWLFTransport<Thermo> &wlft1,
    const polymerKCrossWLFTransport<Thermo> &wlft2)
{
    Thermo t(
        static_cast<const Thermo &>(wlft1) + static_cast<const Thermo &>(wlft2));

    if (mag(t.Y()) < SMALL)
    {
        return polymerKCrossWLFTransport<Thermo>(
            t,
            0,
            wlft1.mu0_,
            wlft1.Tr_,
            wlft1.C1_,
            wlft1.C2_,
            wlft1.rPr_,
            wlft1.Tg_,
            wlft1.TS_,
            wlft1.MS_,
            wlft1.tauStar_,
            wlft1.n_);
    }
    else
    {
        scalar Y1 = wlft1.Y() / t.Y();
        scalar Y2 = wlft2.Y() / t.Y();

        return polymerKCrossWLFTransport<Thermo>(
            t,
            Y1 * wlft1.mu0_ + Y2 * wlft2.mu0_,
            Y1 * wlft1.Tr_ + Y2 * wlft2.Tr_,
            Y1 * wlft1.C1_ + Y2 * wlft2.C1_,
            Y1 * wlft1.C2_ + Y2 * wlft2.C2_,
            1.0 / (Y1 / wlft1.rPr_ + Y2 / wlft2.rPr_),
            Y1 * wlft1.Tg_ + Y2 * wlft2.Tg_,
            Y1 * wlft1.TS_ + Y2 * wlft2.TS_,
            Y1 * wlft1.MS_ + Y2 * wlft2.MS_,
            Y1 * wlft1.tauStar_ + Y2 * wlft2.tauStar_,
            Y1 * wlft1.n_ + Y2 * wlft2.n_);
    }
}

template <class Thermo>
inline Foam::polymerKCrossWLFTransport<Thermo> Foam::operator*(
    const scalar s,
    const polymerKCrossWLFTransport<Thermo> &wlft)
{
    return polymerKCrossWLFTransport<Thermo>(
        s * static_cast<const Thermo &>(wlft),
        wlft.mu0_,
        wlft.Tr_,
        wlft.C1_,
        wlft.C2_,
        1.0 / wlft.rPr_,
        wlft.Tg_,
        wlft.TS_,
        wlft.MS_,
        wlft.tauStar_,
        wlft.n_ );
}

// ************************************************************************* //

How could I call the strainRate value from this library (or call the velocity and compute the strainrate)?



Thanks
Anuar R. Giménez El Amrani