CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   OpenFOAM141dev A new implementation of dynamicKistlerGammaContactAngle (

eberberovic November 13, 2008 06:56

Dear all, I would like to s
Dear all,

I would like to share the implementation of a new boundary condition of dynamic contact angle for interFoam aplication. It is based on the investigations of Hoffman, firstly introduced by Kistler(1993) through a Hoffman function. Details can be found in:

author = {Stephan F. Kistler},
title = {Hydrodynymics of Wetting},
booktitle = {Wettability},
pages = {311--429},
publisher = {Marcel Dekker Inc.},
year = {1993},
editor = {John C. Berg},
chapter = {6},

Basically, one of the problems in assigning an expression for the dynamic contact angle (in addition to modelling the contact line velocity) seems to be the appearance of a static contact angle. The static contact angle shows a hysteresis, and is therefore not single valued. In order to overcome this difficulty, a term is introduced in the function for the dynamic contact angle, through which the influence of the static contact angle is absorbed. The function is of the form, Kistler(1993):

thetaD = F_hoff(Ca + F_hoff<sup>-1</sup>(theta0)),

where the symbols are :
thetaD &ndash; dynamic contact angle,
F_hoff &ndash; the function of Hoffman,
Ca &ndash; Capillary number based on the contact line velocity,
F_hoff<sup>-1</sup> &ndash; inverse of the function of Hoffman,
theta0 - advancing or receding contact angle (depending on the advancing or receding motion).

The function must be solved iteratively due to the appearance of the inverse of the function as one of its own arguments. In the implementation, a C++ function object is introduced defining the inverse Hoffman function, which is solved using RiddersRoot (src/ODE/findRoot). ThetaD is then evaluated using Ca number and the value of F_hoff<sup>-1</sup>. The Ca number is evaluated using the existing model for contact line velocity from dynamicContactAngle in OpenFOAM.

The interFoam application is slightly modified. The createFields.H containes initialization of mixture viscosity and surface tension coefficient, which are passed as volScalarFields to the dynamicKistlerContactAngle class. The mixture viscosity is updated during time steps. The application with the new dynamicKistlerContactAngle class and a sample damBreak case are here:

and a couple of snapshots from damBreak case (on a still relatively coarse grid):

This model for the dynamic contact angle (but with another simpler model for contact line velocity) has been successively used in other numerical codes in predictions of drop impact. Examples can be found in:

author = {S. Sikalo and H. D. Wilhelm and I. V. Roisman and S. Jakirlic and C. Tropea},
title = {Dynamic contact angle of spreading droplets: experiments and simulations},
journal = {Physics of Fluids},
year = {2005},
volume = {17},
number = {6},
pages = {62103},

author = {I. V. Roisman and L. Opfer and C. Tropea and M. Raessi and J. Mostaghimi and S. Chandra},
title = {Drop impact onto a dry surface: Role of the dynamic contact angle},
journal = {Colloids and Surfaces A},
year = {2008},
volume = {322},
pages = {183--191},

I am aware that this might not be the most elegant or general way of implementation. Therefore any further suggestions for improvement are welcome.

With kind regards,
Edin Berberovic.

eberberovic November 13, 2008 07:10

Sorry, the files were too larg
Sorry, the files were too large for upload.
Here is the interFoam solver with only the part containing the dynamicKistlerGammaContactAngle class. It should be put into transportModels/interfaceProperties/gammaContactAngle/dynamicKistlerGammaContactAngle interFoamKistler.tgz damBreak_thetaA_115.tgz

gerado November 13, 2008 09:50

Hi Edin , I am also implement
Hi Edin ,
I am also implementing some models of dynamic contact angle in Openfoam . why you didn t use interfaceproperties.C in your model


eberberovic November 13, 2008 10:24

Dear Gerard, I didn't want
Dear Gerard,

I didn't want to create any objects of type interfaceProperties within the class of type ...FvPatchScalarField just in order to get some data from them. In such case, I believe I would have to use also the twoPhaseProperties objects.

I wanted to pass the values for mu and sigma from the application itself (interFoam).
Therefore I used as a basis the implementation of totalPressure BC.


gerado November 13, 2008 10:57

Dear Edin , I thank you for y
Dear Edin ,
I thank you for your fast answer . you know in interfaceproperties.C the contact angle ist corrected (correctgammacontactangle)at each time step . where did you correct it in your inplementation wenn you don t use interfaceproperties.



eberberovic November 13, 2008 12:00

Dear Gerard, OK, now I see
Dear Gerard,

OK, now I see what you are asking.

Yes, I do use the same function from interfaceProperties.C to correct the contact angle. What I did is the following: I derived the dynamicKistlerGammaContactAngleFvPatchScalarField from gammaContactAngleFvPatchScalarField, the same way as the existing dynamicGammaContactAngleFvPatchScalarField was derived. The same function correctContactAngle from interfaceProperties.c is used. Only the way how thetaD is calculated has changed (I had to put in a function object and make use of RiddersRoot) as well as the way of passing the arguments to the constructor of dynamicKistlerGammaContactAngleFvPatchScalarField.

I did not make it global, I copied the transportModels folder from src into my interFoam application folder, and made an additional folder called dynamicKistlerGammaContactAngle within.

This is what I wanted to post here, but it was to large, so if you want I can email it to you.


gerado November 13, 2008 12:25

Dear Edin, I am really happy
Dear Edin,
I am really happy to have discussed it with you . it would be very nice if you can sent it to me . I have a last question to you . in your code you say
if uwall>0 theta0 = thetaR
if uwall<0 theta0 = thetaA
is that not the opposite?



asaha November 14, 2008 00:17

Hello Edin, Thanks for shar
Hello Edin,

Thanks for sharing your implementation on dynamic contact angle in interFoam. I have the following question:

Will it be possible to modify the equation for dyanamic contact angle in your implementation with the following equation

|tan(theta)| = a*Ca^(1/3) -
b*(lambda)^0.04*Ca^0.293; a = 7.48, b = 3.28, lambda = 10E-8, Ca = Capillary number

I would be extremely happy to hear from you in this regard.

Thanks again.

A A Saha.

eberberovic November 14, 2008 04:46

Gerard, No, I believe that

No, I believe that this is correct. This is the most general approach I could make using the exsisting model for contact line velocity in OpenFOAM.
uwall (small u) is evaluated as a dot product of nWall and Uwall. nWall is calculated using nHat by subtracting from it the component normal to the wall. So nWall is parallel to the wall and points into liquid (due to normalized grad(gamma) in nHat). Similarly, Uwall is calculated by taking U from the first cell near the wall and subtracting its normal component. So Uwall is also parallel to the wall, but points in the direction of motion.
Now if the liquid advances, these two vectors point in opposite sides, and if the liquid recedes the point to the same side. Therefore the advancing motion is determined by uwall<0>0.
I also put in some lines to print out all these values at the wall. Uncomment them and you will see that this is correct at the wall.


I think this should not be a problem. After you calculate the Ca number, you can simply use your expression for theta. In this case you will not need any function objects, since, as far as I can see, this is an ordinary linear equation which you can evaluate directly (i.e. by theta = arctan(...)).
So, remove the function object I put in, and simply replace the expression for theta.


eberberovic November 14, 2008 04:49

Sorry for the mistype: Now if
Sorry for the mistype:
Now if the liquid advances, these two vectors point in opposite sides, and if the liquid recedes they point to the same side. Therefore the advancing motion is determined by uwall<0>0.

eberberovic November 14, 2008 04:51

The advancing motion is determ
The advancing motion is determined by uwall<0.
The receding motion is determined by uwall>0.

asaha November 14, 2008 09:49

Hello Edin, I would appreci
Hello Edin,

I would appreciate if you can share you code and send the same to


a a saha

feijooos October 12, 2009 22:51

First, this is a nice addition to the OF code. Second, I am trying to get this code to work with OF 1.5. I am getting some errors regarding the -llduSolver. Did anybody successfully implement this?

Thanks for the help!

farhagim July 20, 2010 14:59

Hello Edin,

thanks for sharing your information. I could not open your files. Would you please share your code again or is it possible to help me to implement the new one?


Originally Posted by eberberovic (Post 188399)
The advancing motion is determined by uwall<0.
The receding motion is determined by uwall>0.

benni November 4, 2014 06:45

I try to implement Kistler's dynamic contact angle model in OpenFOAM 2.3.0. With a few changes in Edin's files, I almost could get it work.
Can anybody send me his implementation and particularly the files of a case?

Thanks a lot!

All times are GMT -4. The time now is 03:27.