New Viscosity Model - Best Practice for Temporary Values

MaSch · August 23, 2016, 10:52

Hello FOAMers,
currently I am implementing a temperature-dependant viscosity model, to be concise, the Carreau-WLF-Model, for polymer melts. It is very similiar to a Cross-WLF-model. (See equations 2.29 and 2.32 in this preview for more details: http://files.hanser.de/hanser/docs/2..._Leseprobe.pdf)

I have referred to the following materials:

http://www.cfd-online.com/Forums/ope...time-loop.html

http://www.tfd.chalmers.se/~hani/kur...nFoam%20v2.pdf

http://www.cfd-online.com/Forums/ope...mperature.html

Given the fact that the equation is - while not complex - rather long, I tried to split it up across several lines. See this snippet from GermanCarreau.C.

Code:

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

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::GermanCarreau::calcNu() const
{
const volScalarField& Tcurrent=U_.mesh().lookupObject<volScalarField>("T");


Foam::dimensionedScalar Tsp1 = Te_+scalar(50);
Foam::dimensionedScalar termA = (scalar(8.86)*(Tmeasure_ - Tsp1)) / (scalar(101.6)+(Tmeasure_ -Tsp1));
Foam::dimensionedScalar termB = (scalar(8.86)*(Tcurrent - Tsp1)) / (scalar(101.6)+(Tcurrent -Tsp1));
Foam::dimensionedScalar alpha_shift = pow(scalar(10), (termA - termB));
Foam::dimensionedScalar eta_shifted = ( (alpha_shift*A_) / (pow (eta_shifted*scalar(1)*strainRate()*B_, C_) ));

    return (eta_shifted/rho_.value()

    );
}

This fails to compile, resulting in the following errors:

Code:

viscosityModels/GermanCarreau/GermanCarreau.C:58:82: error: conversion from ‘Foam::tmp<Foam::DimensionedField<double, Foam::volMesh> >’ to non-scalar type ‘Foam::dimensionedScalar {aka Foam::dimensioned<double>}’ requested
 Foam::dimensionedScalar termB = (scalar(8.86)*(Tcurrent.internalField() - Tsp1)) / (scalar(101.6)+(Tcurrent.internalField()  -Tsp1));
                                                                                  ^
viscosityModels/GermanCarreau/GermanCarreau.C:60:58: error: conversion from ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’ to non-scalar type ‘Foam::dimensionedScalar {aka Foam::dimensioned<double>}’ requested
 Foam::dimensionedScalar eta_shifted = ( (alpha_shift*A_) / (pow (eta_shifted*scalar(1)*strainRate()*B_, C_) ));
                                                          ^
viscosityModels/GermanCarreau/GermanCarreau.C:62:24: error: could not convert ‘Foam::operator/(const Foam::dimensioned<Type>&, const Foam::dimensioned<double>&) [with Type = double](Foam::dimensioned<double>((* &((const Foam::viscosityModels::GermanCarreau*)this)->Foam::viscosityModels::GermanCarreau::rho_.Foam::dimensioned<Type>::value<double>())))’ from ‘Foam::dimensioned<double>’ to ‘Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> >’
     return (eta_shifted/rho_.value()
                        ^
viscosityModels/GermanCarreau/GermanCarreau.C: In constructor ‘Foam::viscosityModels::GermanCarreau::GermanCarreau(const Foam::word&, const Foam::dictionary&, const volVectorField&, const surfaceScalarField&)’:
viscosityModels/GermanCarreau/GermanCarreau.C:81:5: error: class ‘Foam::viscosityModels::GermanCarreau’ does not have any field named ‘GermanCarreauCoeffs_’
     GermanCarreauCoeffs_(viscosityProperties.subDict(typeName + "Coeffs")),
     ^
viscosityModels/GermanCarreau/GermanCarreau.C:82:22: error: ‘GermanCarreauCoeffs_’ was not declared in this scope
     A_("A", dimless, GermanCarreauCoeffs_),
                      ^
viscosityModels/GermanCarreau/GermanCarreau.C:83:21: error: expected ‘)’ before ‘GermanCarreauCoeffs_’
     B_("B", dimless GermanCarreauCoeffs_),

As you can see, Scalars and volScalarFields cannot be used together.

As far as I understand, the transport equation is evaluated for each point in the flow field - so how would I go about describing my temporary variables in terms of OpenFOAM data types? How would you approach such a problem?

Is my approach of compartmentalizing the equation fundamentally sound?

Can you point to a implented transport model that utilizes a very long formula in a similar manner?

Thank you very much!

Attached, you'll find GermanCarreau.C and GermanCarreau.H:

GermanCarreau.C

Code:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 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 "GermanCarreau.H"
#include "addToRunTimeSelectionTable.H"
#include "surfaceFields.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
namespace viscosityModels
{
    defineTypeNameAndDebug(GermanCarreau, 0);

    addToRunTimeSelectionTable
    (
        viscosityModel,
        GermanCarreau,
        dictionary
    );
}
}


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

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::GermanCarreau::calcNu() const
{
const volScalarField& Tcurrent=U_.mesh().lookupObject<volScalarField>("T");


Foam::dimensionedScalar Tsp1 = Te_+scalar(50);
Foam::dimensionedScalar termA = (scalar(8.86)*(Tmeasure_ - Tsp1)) / (scalar(101.6)+(Tmeasure_ -Tsp1));
Foam::dimensionedScalar termB = (scalar(8.86)*(Tcurrent - Tsp1)) / (scalar(101.6)+(Tcurrent -Tsp1));
Foam::dimensionedScalar alpha_shift = pow(scalar(10), (termA - termB));
Foam::dimensionedScalar eta_shifted = ( (alpha_shift*A_) / (pow (eta_shifted*scalar(1)*strainRate()*B_, C_) ));

    return (eta_shifted/rho_.value()



    );
}


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

Foam::viscosityModels::GermanCarreau::GermanCarreau
(
    const word& name,
    const dictionary& viscosityProperties,
    const volVectorField& U,
    const surfaceScalarField& phi
)
:
    viscosityModel(name, viscosityProperties, U, phi),
    GermanCarreauCoeffs_(viscosityProperties.subDict(typeName + "Coeffs")),
    A_("A", dimless, GermanCarreauCoeffs_),
    B_("B", dimless GermanCarreauCoeffs_),
    C_("C", dimless, GermanCarreauCoeffs_),
    Te_(GermanCarreauCoeffs_.lookup("Te"),
    rho_("rho", dimless, GermanCarreauCoeffs_),
    Tmeasure_("Tmeasure", dimless, GermanCarreauCoeffs_),
    pFactor_("pFactor", dimless, GermanCarreauCoeffs_),



    nu_
    (
        IOobject
        (
            name,
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        calcNu()
    )
{}


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

bool Foam::viscosityModels::GermanCarreau::read
(
    const dictionary& viscosityProperties
)
{
    viscosityModel::read(viscosityProperties);

    GermanCarreau_ = viscosityProperties.subDict(typeName + "Coeffs");

    GermanCarreauCoeffs_.lookup("A") >> A_;
    GermanCarreauCoeffs_.lookup("B") >> B_;
    GermanCarreauCoeffs_.lookup("C") >> C_;
    GermanCarreauCoeffs_.lookup("Te") >> Te_;
    GermanCarreauCoeffs_.lookup("rho") >> rho_;
    GermanCarreauCoeffs_.lookup("Tmeasure") >> Tmeasure_;
    GermanCarreauCoeffs_.lookup("pFactor") >> pFactor_;


    return true;
}


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

GermanCarreau.H

Code:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 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::viscosityModels::GermanCarreau

Description
    An incompressible isothermal Carreau Three-factor viscosity model.

SourceFiles
    GermanCarreau.C

\*---------------------------------------------------------------------------*/

#ifndef GermanCarreau_H
#define GermanCarreau_H

#include "viscosityModel.H"
#include "dimensionedScalar.H"
#include "volFields.H"

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

namespace Foam
{
namespace viscosityModels
{

/*---------------------------------------------------------------------------*\
                           Class GermanCarreau Declaration
\*---------------------------------------------------------------------------*/

class GermanCarreau
:
    public viscosityModel
{
    // Private data

        dictionary GermanCarreau_;

        dimensionedScalar A_;
        dimensionedScalar B_;
        dimensionedScalar C_;
        dimensionedScalar rho_;
        dimensionedScalar Te_;
        dimensionedScalar Tmeasure_;
        dimensionedScalar pFactor_;

        dimensionedScalar Tsp2_;
        dimensionedScalar termA_;
        dimensionedScalar termB_;
        dimensionedScalar eta_shifted_;
        dimensionedScalar alpha_shift_;

        volScalarField nu_;

    // Private Member Functions

        //- Calculate and return the laminar viscosity
        tmp<volScalarField> calcNu() const;


public:

    //- Runtime type information
    TypeName("GermanCarreau");


    // Constructors

        //- Construct from components
        GermanCarreau
        (
            const word& name,
            const dictionary& viscosityProperties,
            const volVectorField& U,
            const surfaceScalarField& phi
        );


    //- Destructor
    ~GermanCarreau()
    {}


    // Member Functions

        //- Return the laminar viscosity
        tmp<volScalarField> nu() const
        {
            return nu_;
        }

        //- Return the laminar viscosity for patch
        tmp<scalarField> nu(const label patchi) const
        {
            return nu_.boundaryField()[patchi];
        }

        //- Correct the laminar viscosity
        void correct()
        {
            nu_ = calcNu();
        }

        //- Read transportProperties dictionary
        bool read(const dictionary& viscosityProperties);
};


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

} // End namespace viscosityModels
} // End namespace Foam

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

#endif

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

GerhardHolzinger · August 24, 2016, 03:36

Hello,

in Line 58 you try to compute a dimensionedScalar with a formula containing a volScalarField. This is equivalent to compute a scalar directly from a vector. You can compute the magnitude of the vector or the inner product to gain a scalar, but otherwise you won't get a scalar out of a vector. The same logic applies to scalar fields.

Code:

Foam::dimensionedScalar termB = (scalar(8.86)*(Tcurrent - Tsp1)) / (scalar(101.6)+(Tcurrent -Tsp1));

Thus, in your model the variables termB, alpha_shift and eta_shifted need to be volScalarFields, unless you do something about Tcurrent, e.g. apply spatial averaging to make a scalar from the field.

MaSch · August 24, 2016, 03:54

First off: Thank you a lot for your very useful feedback, Gerhard.

There is - however - a simple question remaining:

As far as I understand, the viscosity model is evaluated for each point the flow channel. Otherwise, shear thinning [that is, the reduction of the kinematic viscosity nu in areas of high shear rate] could not be observed.

How would I modify

Code:

const volScalarField& Tcurrent=U_.mesh().lookupObject<volScalarField>("T");{

so that Tcurrent is not the entire volScalarField, but rather the temperature of the node for which nu is currently being calculated?

Thanks again!

GerhardHolzinger · August 24, 2016, 04:36

If the variables termB, alpha_shift and eta_shifted are volScalarFields, then your formula is automatically evaluated for all cells of the mesh. There is no need to change the expression for Tcurrent.

MaSch · August 25, 2016, 09:18

Hello Foamers, in the end I did manage to make it work. I used a hack described here:

http://www.cfd-online.com/Forums/ope...tml#post595370

so I could avoid having to manually create new volScalarFields.

Also, I left out the calculation of Tstandard. Now, Ts is defined in the transportProperties file.

A short test revealed that the results are accurate to within ca. 10%, which I blame on shear heating and my godawful mesh.

Attached you'll find the compileable files GermanCarreau.C and GermanCarreau.H. However, the material model only runs with a non-isothermal solver such as the one found here: http://www.dhcae-tools.com/docs/visousHeatingSolver.pdf

I hope you'll find my results useful!

Code:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 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 "GermanCarreau.H"
#include "addToRunTimeSelectionTable.H"
#include "surfaceFields.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
namespace viscosityModels
{
    defineTypeNameAndDebug(GermanCarreau, 0);

    addToRunTimeSelectionTable
    (
        viscosityModel,
        GermanCarreau,
        dictionary
    );
}
}


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

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::GermanCarreau::calcNu() const
{
const volScalarField& Tcurrent=U_.mesh().lookupObject<volScalarField>("T");

volScalarField alpha_shift =Tcurrent;            //instead of constructing a new volScalarField, we "copy" the old one            
alpha_shift *= scalar(0.0);                //we then set all values of the volScalarField to zero
alpha_shift.dimensions().reset(dimless);        //also, since T is a Value in Kelvin, we make the new Field dimensionsless

volScalarField termA =Tcurrent;
termA *= scalar(0.0);
termA.dimensions().reset(dimless);

volScalarField termB =Tcurrent;
termB *= scalar(0.0);
termB.dimensions().reset(dimless);

volScalarField Tnew =Tcurrent;
Tnew.dimensions().reset(dimless);

//alpha_shift += scalar(1.2);
//Tsp1.dimensions().reset(dimLength*dimLength/dimTime);  //turns into a kinematic viscosity

termB = (scalar(8.86)*(Tnew - Ts_)) / (scalar(101.6)+(Tnew -Ts_));
termA = (scalar(8.86)*(Tmeasure_ - Ts_)) / (scalar(101.6)+(Tmeasure_ -Ts_));
alpha_shift = pow(scalar(10), (termA - termB));


     return (nuA_*alpha_shift)/(scalar(1) + pow(alpha_shift*B_*strainRate(), C_));
    
}


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

Foam::viscosityModels::GermanCarreau::GermanCarreau
(
    const word& name,
    const dictionary& viscosityProperties,
    const volVectorField& U,
    const surfaceScalarField& phi
)
:
    viscosityModel(name, viscosityProperties, U, phi),
    GermanCarreauCoeffs_(viscosityProperties.subDict(typeName + "Coeffs")),
    nuA_("nuA", dimViscosity, GermanCarreauCoeffs_),
    B_("B", dimTime, GermanCarreauCoeffs_),
    C_("C", dimless, GermanCarreauCoeffs_),
    Ts_("Ts", dimless, GermanCarreauCoeffs_),
    Tmeasure_("Tmeasure", dimless, GermanCarreauCoeffs_),
    nu_
    (
        IOobject
        (
            name,
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        
        ),    
               calcNu()
    )
{}

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

bool Foam::viscosityModels::GermanCarreau::read
(
    const dictionary& viscosityProperties
)
{
    viscosityModel::read(viscosityProperties);

    GermanCarreauCoeffs_ = viscosityProperties.subDict(typeName + "Coeffs");

    GermanCarreauCoeffs_.lookup("nuA") >> nuA_;
    GermanCarreauCoeffs_.lookup("B") >> B_;
    GermanCarreauCoeffs_.lookup("C") >> C_;
    GermanCarreauCoeffs_.lookup("Ts") >> Ts_;
    GermanCarreauCoeffs_.lookup("Tmeasure") >> Tmeasure_;



    return true;
}


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

Code:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 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::viscosityModels::GermanCarreau

Description
    An incompressible isothermal Carreau Three-factor viscosity model with WLF correction for temperatures

SourceFiles
    GermanCarreau.C

\*---------------------------------------------------------------------------*/

#ifndef GermanCarreau_H
#define GermanCarreau_H

#include "viscosityModel.H"
#include "dimensionedScalar.H"
#include "volFields.H"

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

namespace Foam
{
namespace viscosityModels
{

/*---------------------------------------------------------------------------*\
                           Class GermanCarreau Declaration
\*---------------------------------------------------------------------------*/

class GermanCarreau
:
    public viscosityModel
{
    // Private data

        dictionary GermanCarreauCoeffs_;

        dimensionedScalar nuA_;
        dimensionedScalar B_;
        dimensionedScalar C_;
        dimensionedScalar Ts_;
        dimensionedScalar Tmeasure_;

        volScalarField nu_;



    // Private Member Functions

        //- Calculate and return the laminar viscosity
        tmp<volScalarField> calcNu() const;


public:

    //- Runtime type information
    TypeName("GermanCarreau");


    // Constructors

        //- Construct from components
        GermanCarreau
        (
            const word& name,
            const dictionary& viscosityProperties,
            const volVectorField& U,
            const surfaceScalarField& phi
        );


    //- Destructor
    ~GermanCarreau()
    {}


    // Member Functions

        //- Return the laminar viscosity
        tmp<volScalarField> nu() const
        {
            return nu_;
        }

        //- Return the laminar viscosity for patch
        tmp<scalarField> nu(const label patchi) const
        {
            return nu_.boundaryField()[patchi];
        }

        //- Correct the laminar viscosity
        void correct()
        {
            nu_ = calcNu();
        }

        //- Read transportProperties dictionary
        bool read(const dictionary& viscosityProperties);
};


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

} // End namespace viscosityModels
} // End namespace Foam

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

#endif

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

malv83 · March 6, 2018, 02:41

Hello everyone,

I have seen some tutorials on how to modify the (incompressible) viscosity models (Powerlaw, Casson, etc.) to make them temperature dependent.

I am trying to modify a compressible thermo physical model, such as Sutherland model (for viscosity) which depends on Pressure and Temperature, but not on the strainRate.

The strainRate() is directly available for incompressible viscosityModels, but it is not available for compressible thermo physical model, such as Sutherland.

Does any body know how to make the Sutherland model strainRate dependent? in addition to the Pressure and Temperature dependency?.

I am using OpenFoam 5.x in Windows by BlueCFD.

Thanks in advanced

MaSch · March 16, 2018, 03:51

Dear malv,

while I'm not certain that it's going to work in the compressible flow model, you might try to calculate the strain rate directly from the flow field.

See this useful thread:
StrainRate definition 141 vs 15

which eventually arrives at:

Code:

strainrate= sqrt(2) * mag(symm(fvc::grad(U)))

March 6, 2018, 02:41	Compressible thermophysical transport model	#6
malv83 Member Alberto Join Date: Sep 2013 Posts: 37 Rep Power: 12	Hello everyone, I have seen some tutorials on how to modify the (incompressible) viscosity models (Powerlaw, Casson, etc.) to make them temperature dependent. I am trying to modify a compressible thermo physical model, such as Sutherland model (for viscosity) which depends on Pressure and Temperature, but not on the strainRate. The strainRate() is directly available for incompressible viscosityModels, but it is not available for compressible thermo physical model, such as Sutherland. Does any body know how to make the Sutherland model strainRate dependent? in addition to the Pressure and Temperature dependency?. I am using OpenFoam 5.x in Windows by BlueCFD. Thanks in advanced

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Viscosity model - Compressible flow	Aurelien Thinat	OpenFOAM Programming & Development	3	March 6, 2018 07:08
Adding new temperature dependent viscosity model	vabishek	OpenFOAM Programming & Development	3	May 15, 2016 21:05
coefficient of power-Law viscosity Model in OpenFOAM	Daniel_Khazaei	OpenFOAM Running, Solving & CFD	6	April 5, 2016 04:23
How to modify the viscosity model	mpml	OpenFOAM Running, Solving & CFD	4	October 13, 2010 07:44
Regarding Implementation of Hershel Bulkley Model in Viscosity correlation	Pathak	FLOW-3D	0	August 24, 2009 06:34

August 24, 2016, 03:36		#2
GerhardHolzinger Senior Member Gerhard Holzinger Join Date: Feb 2012 Location: Austria Posts: 339 Rep Power: 28	Hello, in Line 58 you try to compute a dimensionedScalar with a formula containing a volScalarField. This is equivalent to compute a scalar directly from a vector. You can compute the magnitude of the vector or the inner product to gain a scalar, but otherwise you won't get a scalar out of a vector. The same logic applies to scalar fields. Code: Foam::dimensionedScalar termB = (scalar(8.86)*(Tcurrent - Tsp1)) / (scalar(101.6)+(Tcurrent -Tsp1)); Thus, in your model the variables termB, alpha_shift and eta_shifted need to be volScalarFields, unless you do something about Tcurrent, e.g. apply spatial averaging to make a scalar from the field.

August 24, 2016, 03:54		#3
MaSch New Member Join Date: Aug 2016 Posts: 4 Rep Power: 9	First off: Thank you a lot for your very useful feedback, Gerhard. There is - however - a simple question remaining: As far as I understand, the viscosity model is evaluated for each point the flow channel. Otherwise, shear thinning [that is, the reduction of the kinematic viscosity nu in areas of high shear rate] could not be observed. How would I modify Code: const volScalarField& Tcurrent=U_.mesh().lookupObject<volScalarField>("T");{ so that Tcurrent is not the entire volScalarField, but rather the temperature of the node for which nu is currently being calculated? Thanks again!

August 24, 2016, 04:36		#4
GerhardHolzinger Senior Member Gerhard Holzinger Join Date: Feb 2012 Location: Austria Posts: 339 Rep Power: 28	If the variables termB, alpha_shift and eta_shifted are volScalarFields, then your formula is automatically evaluated for all cells of the mesh. There is no need to change the expression for Tcurrent.

March 16, 2018, 03:51		#7
MaSch New Member Join Date: Aug 2016 Posts: 4 Rep Power: 9	Dear malv, while I'm not certain that it's going to work in the compressible flow model, you might try to calculate the strain rate directly from the flow field. See this useful thread: StrainRate definition 141 vs 15 which eventually arrives at: Code: strainrate= sqrt(2) * mag(symm(fvc::grad(U)))