Herschel-Bulkley model definition in HerschelBulkley.C

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

 November 12, 2010, 20:03 Herschel-Bulkley model definition in HerschelBulkley.C #1 New Member   alex Join Date: Jun 2009 Posts: 17 Rep Power: 10 Sponsored Links 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 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

November 13, 2010, 06:51
#2
Senior Member

Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 253
Rep Power: 13
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  |
-------------------------------------------------------------------------------
This file is part of OpenFOAM.

OpenFOAM is free software; you can redistribute it and/or modify it
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();
}

};

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

} // 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  |
-------------------------------------------------------------------------------
This file is part of OpenFOAM.

OpenFOAM is free software; you can redistribute it and/or modify it
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 "surfaceFields.H"

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

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

(
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::AUTO_WRITE
),
calcNu()
)
{}

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

(
const dictionary& 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

November 13, 2010, 17:00
#3
New Member

alex
Join Date: Jun 2009
Posts: 17
Rep Power: 10
Quote:
 Originally Posted by MartinB
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

November 13, 2010, 19:08
#4
Senior Member

Martin
Join Date: Oct 2009
Location: Aachen, Germany
Posts: 253
Rep Power: 13
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

 November 18, 2010, 11:45 #5 New Member   alex Join Date: Jun 2009 Posts: 17 Rep Power: 10 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

 December 27, 2010, 11:14 #6 New Member   alex Join Date: Jun 2009 Posts: 17 Rep Power: 10 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.

February 27, 2017, 11:48
#7
New Member

Térence
Join Date: Feb 2017
Location: Bordeaux, France.
Posts: 27
Rep Power: 2
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
 library.PNG (20.5 KB, 8 views)

Last edited by TemC; March 3, 2017 at 04:28.

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post msrinath80 OpenFOAM Running, Solving & CFD 18 March 3, 2015 06:36 emad FLUENT 5 May 18, 2012 13:25 YJZ ANSYS 1 August 20, 2010 13:57 achaokaoyan Main CFD Forum 0 July 10, 2010 10:52 rystokes CFX 3 August 9, 2009 19:13