[oort]

Hello

In papers I see Herchel-Bulkley model given by:

tau = tau0 + k * (shearRate)^n

if

tau = viscEff * (shearRate)

then

viscEff = (tau0 + k * (shearRate)^n) / (shearRate)

So why does OpenFOAM implement Herschel-Bulkley model with this equation (in /src/transportModels/incompressible/viscosityModels/HerschelBulkley/HerschelBulkley.C
):

{
dimensionedScalar tone("tone", dimTime, 1.0);
dimensionedScalar rtone("rtone", dimless/dimTime, 1.0);
tmp<volScalarField> sr(strainRate());

return (min(nu0_,(tau0_ + k_* rtone *( pow(tone * sr(), n_)
+ pow(tone*tau0_/nu0_,n_))) / (max(sr(), dimensionedScalar
("VSMALL", dimless/dimTime, VSMALL)))));
}

from where comes:

pow(tone*tau0_/nu0_,n_) ?

what are the values "tone", "rtone" and "VSMALL"?

Thanks

[MartinB]

Hi Carlos,

I played around with Herschel-Bulkley some time ago. To implement the (more common?) version presented at Wikipedia (http://en.wikipedia.org/wiki/Herschel–Bulkley_fluid) I used a new viscosity model named "HerschelBulkleySimple".

Here is my code, some comments are added to the .C file:

First the "HerschelBulkleySimple.H":

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
     \\/     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 2 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, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Class
    Foam::viscosityModels::HerschelBulkleySimple

Description
     Herschel-Bulkley non-Newtonian viscosity model.
     The wikipedia.org model is realized.

Authors
    Martin Becker
    becker_martin@t-online.de

    Ahmed El Sawi
    ahmed.el.sawi@gmail.com

SourceFiles
    HerschelBulkleySimple.C

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

#ifndef HerschelBulkleySimple_H
#define HerschelBulkleySimple_H

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

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

namespace Foam
{
namespace viscosityModels
{

/*---------------------------------------------------------------------------*\
                           Class HerschelBulkleySimple Declaration
\*---------------------------------------------------------------------------*/

class HerschelBulkleySimple
:
    public viscosityModel
{
    // Private data

        dictionary HerschelBulkleySimpleCoeffs_;

        dimensionedScalar k_;
        dimensionedScalar n_;
        dimensionedScalar tau0_;
        dimensionedScalar nu0_;

        volScalarField nu_;

    // Private Member Functions

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


public:

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


    // Constructors

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


    // Destructor

        ~HerschelBulkleySimple()
        {}


    // Member Functions

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

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

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


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

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

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

#endif

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

The "HerschelBulkleySimple.C":

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
     \\/     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 2 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, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

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

#include "HerschelBulkleySimple.H"
#include "addToRunTimeSelectionTable.H"
#include "surfaceFields.H"

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

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

    addToRunTimeSelectionTable
    (
        viscosityModel,
        HerschelBulkleySimple,
        dictionary
    );
}
}


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

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::HerschelBulkleySimple::calcNu() const
{
    // used to make the dimension [-] to [s]:
    dimensionedScalar tone("tone", dimTime, 1.0);
    // used to make the dimension [-] to [1/s]:
    dimensionedScalar rtone("rtone", dimless/dimTime, 1.0);
    // Pointer to the strainRate (shear rate gamma):
    tmp<volScalarField> sr(strainRate());

    // For debugging purpose the constants that where read from the dicts and
    // interesting values can be written out:
    //Info << "max(strainRate): " << max(sr()).value()
    //     << ", min(strainRate): " << min(sr()).value()
    //     << endl;
    //
    //Info << "nu0_: " << nu0_ << endl;
    //Info << "tau_: " << tau0_ << endl;
    //Info << "k_: " << k_ << endl;
    //Info << "n_: " << n_ << endl;
    //
    //Info << "tone: " << tone
    //     << ",          rtone: " << rtone << endl;
    //
    //Info << "nu(max(sr)) = " << (min(nu0_, (tau0_ + k_* rtone *
    // ( pow(tone * max(sr()), n_))) / (max(max(sr()), dimensionedScalar
    // ("VSMALL", dimless/dimTime, VSMALL))))) << endl
    //     << "nu(min(sr)) = " << (min(nu0_, (tau0_ + k_* rtone *
    // ( pow(tone * min(sr()), n_))) / (max(min(sr()), dimensionedScalar
    // ("VSMALL", dimless/dimTime, VSMALL))))) << endl
    //     << endl;

    return
    (
      min
      (
        nu0_,
        (
          tau0_             // dimension is [m^2 / s^2]
          + k_* rtone       // term is set to [m^2 / s] * [1 / s] = [m^2 / s^2]
          * ( pow(tone * sr(), n_))  // strain rate dimension is set to
        )                   // [s]*[1/s]=[-]
      / (max(sr(), dimensionedScalar ("VSMALL", dimless/dimTime, VSMALL)))
           // avoid division by zero if strain rate is too small or zero
      ) // return value has dimension [m^2 / s^2] / [1/s] = [m^2 / s]
    );
}


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

Foam::viscosityModels::HerschelBulkleySimple::HerschelBulkleySimple
(
    const word& name,
    const dictionary& viscosityProperties,
    const volVectorField& U,
    const surfaceScalarField& phi
)
:
    viscosityModel(name, viscosityProperties, U, phi),
    HerschelBulkleySimpleCoeffs_(viscosityProperties.subDict(typeName + "Coeffs")),
    k_(HerschelBulkleySimpleCoeffs_.lookup("k")),
    n_(HerschelBulkleySimpleCoeffs_.lookup("n")),
    tau0_(HerschelBulkleySimpleCoeffs_.lookup("tau0")),
    nu0_(HerschelBulkleySimpleCoeffs_.lookup("nu0")),
    nu_
    (
        IOobject
        (
            name,
            U_.time().timeName(),
            U_.db(),
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        calcNu()
    )
{}


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

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

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

    HerschelBulkleySimpleCoeffs_.lookup("k") >> k_;
    HerschelBulkleySimpleCoeffs_.lookup("n") >> n_;
    HerschelBulkleySimpleCoeffs_.lookup("tau0") >> tau0_;
    HerschelBulkleySimpleCoeffs_.lookup("nu0") >> nu0_;

    return true;
}

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

what are the values "tone", "rtone" and "VSMALL"?

The "tone" and "rtone" remove the dimensions, so it can be seen as a "to dimensionless number one", I suppose.

The "VSMALL" gives a very small value, which is used to avoid divisions by zero.

If not already done you should patch the definition of strainRate, as mentioned here:
http://www.cfd-online.com/Forums/ope...1-vs-15-a.html

To add the HerschelBulkleySimple model, you must edit the file
"OpenFOAM-1.6.x/src/transportModels/incompressible/Make/files".

Call "wclean" in "OpenFOAM-1.6.x/src/transportModels/incompressible" and "./Allwmake" in "OpenFOAM-1.6.x/src/transportModels/" to build the library.

Hope this info is useful,

Martin

p.s.:
An example for a "constant/transportProperties" file:

/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.6                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

transportModel  HerschelBulkleySimple;//Newtonian;
nu             nu [ 0 2 -1 0 0 0 0 ] 0.01;

HerschelBulkleySimpleCoeffs
{
    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 18.0; // mu0 / rho; 
    tau0           tau0 [ 0 2 -2 0 0 0 0 ] 3.0; // tau0 / rho;
    k               k [ 0 2 -1 0 0 0 0 ] 0.3; // k / rho;
    n               n [ 0 0 0 0 0 0 0 ] 0.7;
}

rho               rho [ 1 -3 0 0 0 0 0 ] 1350.0;

[oort]

Originally Posted by MartinB

Hi Martin.
Thank you very much for your reply.

Since I have the experimental Herschel-Bulkley parameters in SI units I think that I can assume:

tau0 = tau0_experimental / density

k = k_experimental / density

(this is also what is comented in transportProperties file)

I think that is an "abuse" say that k has m^2/s dimensions because it should be m^2*s^n/s^2, but since the power was made dimensionless numerically is the same.

So if I print VSMALL in any solver this value is printed in the screen. VSMALL is only a defined constant with a very small value right?

Why there is the parameter "rho" in transportProperties? I think is unecessary if original SI tau0 and k parameters are divided by density.

About the patch of strainRate... the original post is from version 1.5. Version 1.7 still have

return mag(symm(fvc::grad(U_)));

If this was wrong wouldn't be corrected after 2 years?

Thanks,
carlos

[MartinB]

Hi Carlos,

tau0 = tau0_experimental / density

k = k_experimental / density

(this is also what is comented in transportProperties file)

I do assume the same.

So if I print VSMALL in any solver this value is printed in the screen. VSMALL is only a defined constant with a very small value right?

Right!

Why there is the parameter "rho" in transportProperties? I think is unecessary if original SI tau0 and k parameters are divided by density.

It's unnecessary for most standard solvers of OpenFOAM, right. I defined it there because I usually write out results for pressure in [Bar], as it is usual practice in my job. Furthermore it's for documentation purposes so that the material can be identified easily at a later point of time.

About the patch of strainRate... the original post is from version 1.5. Version 1.7 still have

return mag(symm(fvc::grad(U_)));

If this was wrong wouldn't be corrected after 2 years?

Well, it's not wrong since there exist some material models using the strain rate in this definition. But to get the same results as from other FEM software with the same material parameters, it was necessary for us to change the strain rate definition.

Martin

[oort]

Hi.

I think I have some problem...

I created the folder

~/OpenFOAM/carlos-1.7.x/run/HerschelBulkleySimple

and I copied the files "HerschelBulkleySimple.H" and "HerschelBulkleySimple.C" to this folder.
Then I copied the files "files" and "options" to

~/OpenFOAM/carlos-1.7.x/run/HerschelBulkleySimple/Make.

file "files":

HerschelBulkleySimple.C

LIB = $(FOAM_USER_LIBBIN)/HerschelBulkleySimple

file "options":

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/transportModels/incompressible/lnInclude
    
LIB_LIBS = \
    -lfiniteVolume

Then I compiled with "wmake libso" and i think that the compilation is OK.

Now I have also the folders "lnInclude" and "Make/linux64GccDPOpt".

I also have the file "HerschelBulkleySimple.so" in the folder

/home/carlos/OpenFOAM/carlos-1.7.x/lib/linux64GccDPOpt

My problem is that when I run blockMesh and checkMesh appears a warning in beggining:

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  * * //
Create time

--> FOAM Warning : 
    From function dlLibraryTable::open(const fileName&  functionLibName)
    in file db/dlLibraryTable/dlLibraryTable.C at line 78
    could not load  /home/carlos/OpenFOAM/carlos-1.7.x/lib/linux64GccDPOpt/HerschelBulkleySimple.so:  undefined symbol: _ZN4Foam14viscosityModel8typeNameE

Creating block mesh from
     "/home/carlos/OpenFOAM/carlos-1.7.x/run/canalRectangular_WG35X002_1.296_25.211/constant/polyMesh/blockMeshDict"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  * * //
Create time

--> FOAM Warning : 
    From function dlLibraryTable::open(const fileName&  functionLibName)
    in file db/dlLibraryTable/dlLibraryTable.C at line 78
    could not load  /home/carlos/OpenFOAM/carlos-1.7.x/lib/linux64GccDPOpt/HerschelBulkleySimple.so:  undefined symbol: _ZN4Foam14viscosityModel8typeNameE

Creating block mesh from
     "/home/carlos/OpenFOAM/carlos-1.7.x/run/canalRectangular_WG35X002_1.296_25.211/constant/polyMesh/blockMeshDict"

However when I run the nonNewtonianIcoFoam solver it gives no errors or warnings.

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh, no clear-out for time = 0

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model HerschelBulkleySimple

Starting time loop

Time = 0.001

Courant Number mean: 0 max: 0.360239
DILUPBiCG:  Solving for Ux, Initial residual = 1, Final residual = 8.14969e-06, No Iterations 8
DILUPBiCG:  Solving for Uy, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for Uz, Initial residual = 0, Final residual = 0, No Iterations 0
DICPCG:  Solving for p, Initial residual = 1, Final residual = 9.89034e-07, No Iterations 822
time step continuity errors : sum local = 8.91534e-10, global = 5.93651e-12, cumulative = 5.93651e-12
DICPCG:  Solving for p, Initial residual = 0.00243736, Final residual = 9.80171e-07, No Iterations 668
time step continuity errors : sum local = 3.54787e-07, global = 1.42402e-09, cumulative = 1.42996e-09
ExecutionTime = 193.75 s  ClockTime = 200 s

Time = 0.002

This is my setup:

system/controlDict:

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

application     nonNewtonianIcoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         3;

deltaT          0.001;

writeControl    runTime;;

writeInterval   0.1;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression uncompressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

libs ("HerschelBulkleySimple.so");


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

constant/transportProperties:

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

// nu(water+glycerol(35%)+xanthan(0.02%)) @ t = 30 ºC

transportModel  HerschelBulkleySimple;

HerschelBulkleySimpleCoeffs
{
    nu0             nu0 [ 0 2 -1 0 0 0 0 ] 0.000076940; // mu0 / rho; 
    tau0            tau0 [ 0 2 -2 0 0 0 0 ] 0.000033814; // tau0 / rho;
    k               k [ 0 2 -1 0 0 0 0 ] 0.000008426; // k / rho;
    n               n [ 0 0 0 0 0 0 0 ] 0.855757552;
}

rho               rho [ 1 -3 0 0 0 0 0 ] 1087.13909;

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

Are the warnings in blockMesh and checkMesh very important and conduct to errors in the solver?

Thanks

[oort]

Well... I don't know if the error messages produces any error in the solving, however the results are similar to the ones expected: flattened profile in the center and a slightly raise in the sides.

[TemC]

Hi foamers,

I'am working with simpleFoam, and I want to implement a new Herschel-Bulkley model based on what Sir "MartinB" proposed at the beginning of this post. I have already done all the modifications that he suggested, an now I'am just stucked at the final step, meaning:

- Call "wclean" in "OpenFOAM-1.6.x/src/transportModels/incompressible"
- and "./Allwmake" in "OpenFOAM-1.6.x/src/transportModels/" to build the library.

When I try to build the library, I got the error depicted in the attachment. Does somebody got the same message?

Thanks in advance, and have a nice week.

Regards.