[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Getfem-users] Problem with Hertzian contact

From: David Danan
Subject: Re: [Getfem-users] Problem with Hertzian contact
Date: Mon, 9 May 2022 15:26:09 +0200

Dear Kostas,

First, as always, thank you for your dedication! It is very much appreciated!
Let me answer to each point you raised

1)2) Sorry, it seems it seems i was not very clear once again.

When i said negligible i meant that such a modification cannot explain the large gap between the analytical solution and numerical solution in the script i sent in my first email (the full disk), which is exactly the problem i am having (point 2) you raised).
I suspected i made a mistake in the analytical formula for the pressure/max radius but i was unable to find it.

For completeness, please find enclosed the following results in a .ppt:
-The half disk case with the numerical solution, the numerical solution modified with the modification regarding the multiplier value and the analytical solution (with the mesh associated)
-The full disk case with the numerical solution, the numerical solution modified with the modification regarding the multiplier value and the analytical solution (with the mesh associated)
Regarding 3), the mesh is refined enough near the boundary but the load was indeed very low, you are right.

Half disk case:
I am afraid these pictures raised even more questions rather than answer some. The strange behaviour of the modified solution (converted multiplier in my previous email) on specific nodes seems to be related to the interpolation of det(F) equal to 0 on these nodes. But I don't understand how it could be the case.

I have tried to investigate a little:
-Even Norm(F) seems to behave like that on these nodes
-I expected the new multiplier to "explode" on these nodes (because division by 0) but it is 0 instead. The interpolation of det(F) and (1/det(F)) is 0 in both cases on these nodes
-I have tried to plot the value of det(F) on the integration points of the contact boundary; there it seems reasonable (value between 0.9 and 1). What is happening here?
-Increasing the order of integration does not solve the issue, neither does increasing the finite element order (but i can test more thoroughly if you think the explanation lies there)

Full disk case:
I observe the same phenomena for det(F) and still need to find what is wrong with the analytical formula.
At least, the integral of the numerical contact multiplier over the contact boundary is consistent with the load which, I can reasonably assume, is not the case for the analytical formula.

4) I see, it might be interesting to try to apply this to the half disk case to check whether i can improve the performances even further.
For now, when using a neo hookean law instead of a linear elasticity law , it fails to converge. I am fiddling with the solver parameter.

Kind regards and thanks for your help!

Le sam. 7 mai 2022 à 14:05, Konstantinos Poulios <> a écrit :
Dear David,

There are a few points to clarify
1) Regarding the difference of 1st PK traction versus Cauchy traction, your result does not let conclude that is negligible. In fact it looks similar at the top, to the deviation that you showed before.
We will need you to show us a much finer resolution and a valid comparison with the analytical Hertzian pressure curve in order to conclude about this point

2) The last comparison you showed us cannot be valid because the total load of the analytical result is much lower than the total load in the numerical result. So you need to fix that first.

3) The result you showed is very coarse, either your load is too little all the mesh size too large.

4) Regarding the quote from our paper. It makes sense, it says that in general for a very large contact width you should not expect the Hertzian solution to be valid, as Andriy basically pointed out. Moreover, for large contact width, the real pressure distribution will depend on the constitutive law. Only if you use the same hyperelastic law that we used in our example 1, you will get a pressure distribution that is quite close to Hertzian, but how close, I cannot tell you. We do not have any very fine comparison in the paper. And it is actually not so interesting either since the result is anyway material-dependent.

Best regards

On Fri, May 6, 2022 at 7:59 PM David Danan <> wrote:
Dear all,

Following your suggestions, please find enclosed the deformed disk result and the comparison between the analytical solution (max displacement is small compared to the disk radius=200), the original multiplier and the multiplier converted.

I have refined the mesh a bit near the interface to have more points in contact.
I had also checked beforehand the value of " (Det(F)*Norm(Inv(F)'*Normal)) " at the contact boundary Gauss points and it was close to 1 everywhere (which was expected in the small deformation theory). Thus, the multiplier and the multiplier converted are almost the same.

The "converted multiplier" is interpolated on the nodes with the following commands, where lambda is the contact multiplier
    gwflMult = "lambda/(Det(Id(2)+Grad_u)*Norm(Inv(Id(2)+Grad_u)'*Normal))"
    newContactMult=model.interpolation(gwflMult, np.transpose(nodesOnContactBoundary), mesh, boundaryTags["CONTACT_BOUNDARY"])

Also, even if the results were any different, quoting your 2015 article

" It should be noted that the studied problem is not strictly a
Hertzian contact due to the high load resulting to a contact width of considerable size compared
to the disc radius. However, numerical experiments showed that the considered material law (35)
preserves the validity of the Hertzian solution in the finite strain regime quite accurately"

Very strange...


Le ven. 6 mai 2022 à 15:55, Konstantinos Poulios <> a écrit :
Dear Andriy,

Yes I agree this is another possible explanation. Basically when the ratio between contact width and contact radius becomes large, the Hertzian solution is not valid anymore. I had somehow assumed that David's plot was for a reasonably small contact width, but you are of course right one should also check that.


On Fri, May 6, 2022 at 12:53 PM Andriy Andreykiv <> wrote:
Kostas, is it correct to differentiate between 1st Piola Kirchhoff stress and Cauchy stress when we are within small deformation theory?
I suspect that when the strains in the David's simulation become large, Hertz analytical solution doesn't apply anymore because of small deformation assumption in Hertz contact

On Fri, 6 May 2022 at 12:20, Konstantinos Poulios via <> wrote:
ah ok, I think I know what the issue is here. The Lagrange multiplier for imposing the contact condition corresponds to a 1st Piola Kirchhoff stress, while the analytical solution gives you the Cauchy stress. You can use Nanson's formula to convert between the two. In your model the Lagrange multiplier is like q in the upper left case of the following table:

to convert this q to the equivalent one for the upper right case you need to postprocess your Lagrange multiplier and export the following quantity instead:


I hope it works.

It is actually a good observation. I don't recall accounting for this difference in the presentation of the results in our paper from 2015, but maybe our meshes were too coarse to notice the difference, compared to your simulations.


On Fri, May 6, 2022 at 11:52 AM David Danan <> wrote:
Dear Kostas,

first, thank a lot for the feedback, the convergence is much cleaner now!
Thanks also for sharing this script here, it contains several features i had never use so far, it's quite interesting/intriguing.

However, i am afraid i did not properly explained what my issue was.
Somehow, i was able to make the solver converge even without enabling some options for the solve and the result was qualitatively correct, as far as i could tell.
 But it turned out there is a large difference between the pressure predicted by the Hertz theory and the result i get and i was wondering whether i am doing something wrong. I tried to list the potential suspects in the first email

Regarding your remarks:
-Yes,using a punctual load bother me a lot, which is why i wonder if i am supposed to do that to reproduce Hertz analytical case.
-Regarding the unit, i wanted to stay in the small deformation assumption as i was not sure the analytical solution is supposed to hold for large deformation. Therefore, here, the Saint Venant Kirchhoff "becomes" a mere Hooke's law.
Do you have any idea how to reproduce the analytical results?

Thanks to your help regarding the parameters fiddling, i was at least able to reproduce a different case where there is no ponctual force (surfacic load instead) with a half-disk and a different analytical formula and...It worked (see the results enclosed for the pressure comparison on the contact boundary)!
Before, it failed to converge. However, this formula is not really well documented at all, i would rather have a comparison with the classical formula.

Thanks in advance for your help,

PS: sorry for the duplicate, i forgot to put the mailing list in copy previously

Le ven. 6 mai 2022 à 09:23, Konstantinos Poulios <> a écrit :
Dear David,

With the attached script you should also be able to reproduce the results for the first example of our paper
(the script is only slightly modified compared to the original one, so that it can run on python3 and the current version of GetFEM).

Best regards

On Fri, May 6, 2022 at 8:45 AM Konstantinos Poulios <> wrote:

Your script looks very nice and clean. It just needs a bit of fiddling with the augmentation parameter and the solver. The attached version produces this result

Just some general comments:
- Point constraints/loads make almost never sense (unless you know beforehand the reaction force at a constrained point will be zero)
- Maybe there is some issue with your units, I needed to increase the load to a very large value to get some visible deformations
- The Saint-Venant Kirchhoff material law, that we often use for demonstrations, is a "toy" material law. If you do anything with moderate/large strains, use some proper hyperelastic law
- The augmentation parameter should be in the same order of magnitude like the modulus of elasticity, from 1-2 orders of magnitude below to 1-2 orders of magnitude above.

Best regards

On Thu, May 5, 2022 at 8:04 PM David Danan <> wrote:
Dear Getfem-users,

My issue might be utterly simple, but i cannot find what i am doing wrong...

I am trying to compare the solution provided by Getfem to the analytical solution of an Hertzian contact in 2D (namely, a deformable ball against a perfectly rigid plane) but, so far, i do not get correct results.

Please, find enclosed my script. I use gmsh python api to generate the mesh but, if you do not want to use it, you can comment the line to generate the mesh in the script and juste use the mesh also enclosed.

Here are some of my concerns:
-Are the formula used in the script for the analytical contact pressure/contact surface area etc... correct in 2D?
-Is it correct to use a punctual nodal force in such a case? Does it make sense in this context?
-Is it correct to compare the analytical pressure to the lagrange multipliers?

I have also tried to use a testcase provided in Code_Aster documentation and i got much better results, although not perfect but within the tolerance given in the documentation below.
It might be because these formula are for 3D only but i am not sure yet.

Thanks in advance for your help,
kind regards,

Attachment: HertzContact.pptx
Description: MS-Powerpoint 2007 presentation

reply via email to

[Prev in Thread] Current Thread [Next in Thread]