CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM

Herschel-Bulkley model definition in HerschelBulkley.C

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By TemC

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 12, 2010, 20:03
Default Herschel-Bulkley model definition in HerschelBulkley.C
  #1
New Member
 
alex
Join Date: Jun 2009
Posts: 17
Rep Power: 16
oort is on a distinguished road
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
oort is offline   Reply With Quote

Old   November 13, 2010, 06:51
Default
  #2
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 255
Rep Power: 21
MartinB will become famous soon enough
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":
Code:
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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":
Code:
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  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;
}

// ************************************************************************* //
Quote:
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:
Code:
/*--------------------------------*- 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;

Last edited by MartinB; November 13, 2010 at 07:01. Reason: Added a constant/transportProperties example
MartinB is offline   Reply With Quote

Old   November 13, 2010, 17:00
Default
  #3
New Member
 
alex
Join Date: Jun 2009
Posts: 17
Rep Power: 16
oort is on a distinguished road
Quote:
Originally Posted by MartinB View Post
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
oort is offline   Reply With Quote

Old   November 13, 2010, 19:08
Default
  #4
Senior Member
 
Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 255
Rep Power: 21
MartinB will become famous soon enough
Hi Carlos,

Quote:
tau0 = tau0_experimental / density

k = k_experimental / density

(this is also what is comented in transportProperties file)
I do assume the same.

Quote:
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!

Quote:
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.

Quote:
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
MartinB is offline   Reply With Quote

Old   November 18, 2010, 11:45
Default
  #5
New Member
 
alex
Join Date: Jun 2009
Posts: 17
Rep Power: 16
oort is on a distinguished road
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":

Code:
HerschelBulkleySimple.C

LIB = $(FOAM_USER_LIBBIN)/HerschelBulkleySimple
file "options":

Code:
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:

Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  * * //
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"
Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  * * //
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.

Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
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:

Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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:

Code:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// 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 is offline   Reply With Quote

Old   December 27, 2010, 11:14
Default
  #6
New Member
 
alex
Join Date: Jun 2009
Posts: 17
Rep Power: 16
oort is on a distinguished road
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.



oort is offline   Reply With Quote

Old   February 27, 2017, 11:48
Default
  #7
Member
 
Emery
Join Date: Feb 2017
Location: France.
Posts: 33
Rep Power: 9
TemC is on a distinguished road
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.
Attached Images
File Type: png library.PNG (20.5 KB, 38 views)
ksrinivasasagar likes this.

Last edited by TemC; March 3, 2017 at 04:28.
TemC is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Herschel Bulkley model Problem emad FLUENT 6 September 12, 2018 14:35
Superlinear speedup in OpenFOAM 13 msrinath80 OpenFOAM Running, Solving & CFD 18 March 3, 2015 06:36
Low Reynolds k-epsilon model YJZ ANSYS 1 August 20, 2010 14:57
species transport model or mixture model? achaokaoyan Main CFD Forum 0 July 10, 2010 11:52
multi fluid mixture model issue rystokes CFX 3 August 9, 2009 20:13


All times are GMT -4. The time now is 09:32.