CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM (http://www.cfd-online.com/Forums/openfoam/)
-   -   Herschel-Bulkley model definition in HerschelBulkley.C (http://www.cfd-online.com/Forums/openfoam/82002-herschel-bulkley-model-definition-herschelbulkley-c.html)

 oort November 12, 2010 20:03

Herschel-Bulkley model definition in HerschelBulkley.C

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 November 13, 2010 06:51

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;```

 oort November 13, 2010 17:00

Quote:
 Originally Posted by MartinB (Post 283313)
Hi Martin.

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

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

Thanks,
carlos

 MartinB November 13, 2010 19:08

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

 oort November 18, 2010 11:45

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 December 27, 2010 11:14

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.

http://www.peniche-online.com/images/hb_01.jpg

http://www.peniche-online.com/images/hb_02.jpg

 All times are GMT -4. The time now is 19:08.