CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   A newbie question (

janusz January 7, 2008 16:00

Dear all: Any advice on the
Dear all:

Any advice on the topic would be greatly appreciated especially since I am new to OpenFoam.

The problem I am trying to solve is to model the behavior of a visco-plastic material (e.g. a toothpaste) using a variation of a Bingham model, where typically:

shear_stress = yield_stress + nu*shear_rate

In the past I used to approach such tasks in FLUENT coding an appropriate UDF for apparent viscocity.

In the problem yield stress is spatially dependent (x,y,z), therefore this information needs to be acquired by the model .

A question then:

How can I add a new material model to OpenFoam?
What solver would be best suited for this task? This is a relatively low Re task (>3000).

Best regards & appreciate any response,
Janusz Goldasz

hjasak January 7, 2008 16:17

Welcome. First, have a look a
Welcome. First, have a look at simpleFoam to see how it picks the viscosity model - it is called transportModel. Then, go to:

OpenFOAM-1.4.1-dev/src/transportModels/incompressible/viscosityModels and you will see a family of viscosity models like the one you need. Start from the most similar one, e.g. BirdCarreau, copy it, change all instances of "BirdCarreau" in the file and in the business bit calculate your viscosity. Thus (I'm looking at the existing Bird-Carreau):

//- Calculate and return the laminar viscosity
tmp<volscalarfield> BirdCarreau::calcNu() const
+ (nu0_ - nuInf_)
*pow(scalar(1) + sqr(k_*strainRate()), (n_ - 1.0)/2.0);


This should be clear...

Enjoy OpenFOAM,


hjasak January 7, 2008 16:19

P.S. A more descriptive thread
P.S. A more descriptive thread name in the future please - people need to be able to search for this. Thus:

How to implement own viscosity model?


visco-plastic material, Bingham model


annikagram April 16, 2008 04:44

On the Bingham model: non-newt
On the Bingham model: non-newtonian flow

Has anyone succeeded in implementing the Bingham model shear_stress = yield_stress + nu*shear_rate or maybe Hershel-bulkley at src/transportModels/incompressible/viscosityModels? It is a visco plastic material model including yield stress, which I am not quite sure how to implement into existing models, where no shear stress is present. Thank you!

eugene April 16, 2008 05:08

We are in the process of imple
We are in the process of implementing a Hershel-Bulkley model. If you are willing to wait a few months I can send you a copy when when the job is done.

janusz April 16, 2008 05:12

The subject was discussed here
The subject was discussed here:

Janusz Goldasz

annikagram April 16, 2008 05:34

Dear Eugene, sounds great!
Dear Eugene,

sounds great! Where and for what use are you implementing the Hershel-Bulkley model? I am currently modelling fresh concrete flow, so Hershel-Bulkley would fit rather nicely. My e-mail address is:

annikagram April 16, 2008 05:37

Dear Janusz, thank you. Did
Dear Janusz,

thank you. Did you succeed with the Bingham model? I am rather curious how you introduced the yield stress into the existing Coss or BirdCarreau model. I was thinking they are not stress based, but viscosity based material models?

eugene April 16, 2008 06:08

I'm not sure what it is for. P
I'm not sure what it is for. Please mail me directly at the end of May for an update, otherwise I will forget.

annikagram April 16, 2008 06:44

Dear Eugene, what e-mail ad
Dear Eugene,

what e-mail address do you have?

andersking April 17, 2008 02:48

Hi, Please find attached a

Please find attached a Herschel-Bulkley model, and test case. I wrote this for OF-1.3, but never got around to using it/validating it. I've updated it for OF-1.4, but probably still needs to be validated (and checked for errors).

The implementation is based on the one described in the fluent user manual, ie

nu = (tau0 + k * (pow(strainRate(),n) - pow(tau0/nu0,n)))/strainRate();

the parameters are
\nbsp \nbsp tau0 - yield stress
\nbsp \nbsp nu0 - pre-yield viscosity
\nbsp \nbsp k - consistency index
\nbsp \nbsp n - power index

To imitate a bingham plastic, k becomes the viscosity after yield, n=1 and nu0 should be set to a very high value.

I've attached a diff against the OF-extend version of OF (should work on the standard I think), or alternatively just the required source files.

to use the source files HerschelBulkley.tgz needs to be extracted in the src/transportModels/incompressible/viscosityModels directory.

cd OF-1.4.1/src/transportModels/incompressible/viscosityModels
tar -zxvf HerschelBulkley.tgz

and then the file src/transportModels/incompressible/Make/files needs to be modified to add the model by adding
viscosityModels/HerschelBulkley/HerschelBulkley.C after the line

then you will need to rebuild libincompressibleTransportModels, ie.

cd OF-1.4.1/src/transportModels/incompressible
wmake libso .

I have also attached an example case (pipe flow)

Andrew HerschelBulkley.tgz HerschelBulkley.diff HerschelBulkleyCase.tgz

andersking April 17, 2008 02:51

PS. sorry, \nbsp should be \
PS. sorry, \nbsp should be \ ch { nbsp }

andersking April 17, 2008 02:55

PPS. I just noticed that I hav
PPS. I just noticed that I have nu0 as mu0 in the code. the case should work as is though.

gopala June 25, 2008 07:53

Hello, Has any one succeede

Has any one succeeded in validating the HerschelBulkley model of Andrew?

I implemented the Bingham model in OpenFOAM-1.4 and while validating I see a strong dependence on the time step size. For a simple channel flow case, I get different pressure fields for different values of time step size. Any suggestions?

rvajapey November 30, 2008 19:14

Hello OpenFOAM experts I am
Hello OpenFOAM experts

I am a beginner and have been trying to implement these models into OpenFOAM-1.4.1-dev version. I have a few questions about these models. I can see the implementation of Herschel Bulkley in OpenFOAM-1.5, what happens when the strainRate() becomes zero? Does the model have an artificial cap on the maximum value of nu so that it does not blow up?

Sorry if the question sounds stupid.

Thank you very much for your time....

rvajapey November 30, 2008 19:24

Hello OpenFOAM experts I am
Hello OpenFOAM experts

I am also trying to implement the Giesekus model. The expression for the model is given:

(taup_+lambda_*(fvm::ddt(taup_)+fvm::div(phi_,taup _)
-((gradU_.T()&taup_)+(taup_&gradU_)))+ (alpha_*lambda_*(taup_&taup_)/etap_)) == 2*etap_*gradU_);

All my tensors here are symmetric. I have noticed that the term ((gradU_.T()&taup_)+(taup_&gradU_))+(alpha_*lambda _*(taup_&taup_)/etap_));
is blowing up and I do not understand the reason. I have tried the linear Maxwell model which is the same as above without the ((gradU_.T()&taup_)+(taup_&gradU_))+(alpha_*lambda _*(taup_&taup_)/etap_));
term and that seems to work fine. I cannot use Sp, SuSp because they cannot work on symmetric tensors.

I would be really happy if anyone could suggest what is wrong with the equation.

Thank you very much for your time....


rvajapey December 1, 2008 16:24

Hello OpenFOAM experts I ha
Hello OpenFOAM experts

I have been trying to model the Extended form of Herschel Bulkley equation. The equation has the form as given

nu = min(nu0, (tau0 + k * (pow(strainRate(),n)+nuInf)/strainRate())

The term pow(tau0/nu0,n) (present in the FOAM-1.5 version) is missing here. Can anybody explain to me the significance of this term?

Thank you very much for your time.


prjohnston January 19, 2009 22:15

Hello, I am trying to imple

I am trying to implement a generalised power law model for blood viscosity (Ballyk et al, Biorheology, 31 (5), pp 565-586 (1994). In this model the exponent in the power law is a function of the local strain rate. My implemented calcNu() funciton is:

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

Foam::viscosityModels::GeneralisedPowerLaw::calcNu () const
tmp<volscalarfield> sr(strainRate());
return ((nu0_*exp(-(scalar(1)+(sr())/a_)*exp(-b_/(sr())))+nuInf_)*pow(sr(),(nInf_-n0_*e xp(-(scalar(1)+sr()/c_)*exp(-d_/sr()))-scalar(1))));

Unfortunately, this fails to compile giving the following error messages:
Making dependency list for source file viscosityModels/GeneralisedPowerLaw/GeneralisedPowerLaw.C
SOURCE=viscosityModels/GeneralisedPowerLaw/GeneralisedPowerLaw.C ; g++ -m32 -Dlinux -DDP -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast -O3 -DNoRepository -ftemplate-depth-40 -I.. -I/home/peter/OpenFOAM/OpenFOAM-1.5/src/finiteVolume/lnInclude -IlnInclude -I. -I/home/peter/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude -I/home/peter/OpenFOAM/OpenFOAM-1.5/src/OSspecific/Unix/lnInclude -fPIC -pthread -c $SOURCE -o Make/linuxGccDPOpt/GeneralisedPowerLaw.o
/home/peter/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/GeometricScalarField.C: In function 'Foam::tmp<foam::geometricfield<double,> > Foam::pow(const Foam::GeometricField<double,>&, const Foam::tmp<foam::geometricfield<double,> >&) [with PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]':
viscosityModels/GeneralisedPowerLaw/GeneralisedPowerLaw.C:61: instantiated from here
/home/peter/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/GeometricScalarField.C: 220: error: no matching function for call to 'Foam::dimensioned<double>::dimensioned(const char [2], double, const Foam::dimensionSet&)'
/home/peter/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/dimensionedType.C:112: note: candidates are: Foam::dimensioned<type>::dimensioned(const Foam::word&, const Foam::dimensionSet&, Foam::Istream&) [with Type = double]
/home/peter/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/dimensionedType.C:98: note: Foam::dimensioned<type>::dimensioned(const Foam::word&, Foam::Istream&) [with Type = double]
/home/peter/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/dimensionedType.C:85: note: Foam::dimensioned<type>::dimensioned(Foam::Istream &) [with Type = double]
/home/peter/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/dimensionedType.H:95: note: Foam::dimensioned<type>::dimensioned(const Type&) [with Type = double]
/home/peter/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/dimensionedType.C:73: note: Foam::dimensioned<type>::dimensioned(const Foam::word&, const Foam::dimensioned<type>&) [with Type = double]
/home/peter/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/dimensionedType.C:60: note: Foam::dimensioned<type>::dimensioned(const Foam::word&, const Foam::dimensionSet&, Type) [with Type = double]
/home/peter/OpenFOAM/OpenFOAM-1.5/src/OpenFOAM/lnInclude/dimensionedScalarFwd.H: 42: note: Foam::dimensioned<double>::dimensioned(const Foam::dimensioned<double>&)
make: *** [Make/linuxGccDPOpt/GeneralisedPowerLaw.o] Error 1

I believe the error is due to the fact that in the function pow(a,b), b must be a scalar. In my model b would be a volScalarField (I think).

Is there some way to overcome this problem? Perhaps I have misinterpreted the error messages and there is something else wrong. Any insights would be appreciated.



eldudereeno January 21, 2009 10:11

Hello everybody, I am very
Hello everybody,

I am very new with OpenFOAM and having loosy beginners questions :-)
My questions are about the HerschelBuckley viscosity model.

I always get errors with my dimensions.
Could anybody explain me what i have to insert for "nu0,tau0,k,n"

At the moment i have:
0 2 -1 0 0 0 0 for nu0
0 2 -1 0 0 for tau0
0 2 -1 0 0 0 0 for k
0 0 0 0 0 0 0 for n

which actually does not make sence but its the combination with the least errors.

Thanks a lot

gwierink January 21, 2009 10:38

Hi Peter, I think the dimen
Hi Peter,

I think the dimensions should be:

[0 2 -1 0 0 0 0] for nu0
[0 2 -2 0 0 0 0] for tau0
[0 2 -2 0 0 0 0] for k
[0 0 0 0 0 0 0] for n

What solver are you using?

Rgds, Gijsbert

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