Timo,

the problem with numerical contact is, as others have mentioned, the overclosure. Theory formulations consider a perfect geometry, but when dealing with a discretization you will have a node penetrating the opposite element. This is where the overclosure/contact formulation comes into play. Penalty methods are easier to program, but suffer from numerical instability: the ideal penalty is infinite, but then you don't get a solution (infinite times zero is indetermined), and if you relax the penalty value, you might not capture the mechanics you are looking for.

An approach that has worked for me in the past is to start with a somewhat relaxed penalty value (e.g. equal to the bulk modulus or a small integer fraction of that). Once you get convergence, increase the penalty incrementally and refine your results, starting from the last converged solution. This way you can increase the penalty to values that would otherwise lead to divergence.

Depending on the software you are using, you may consider using a different contact algorithm, e.g. Lagrange multipliers, which are a much more elegant formulation and do not suffer from this instability...

Hope this helps,

Fernando

In my experience, no other method than 'trial and error' checking results against experimental evidence. Elastic properties are not sufficient: also friction and adhesion play a role.

Only by a contact interaction, perhaps you can adjust the penality value of the contact to fit an experiemental result

Clear trial and error but one can try \epsilon = 1/E

Agreed with Frédéric in case of very good surface quality, else stiffness increases with overclosure.

So, You do not think that it would be possible to apply the Hertz contact theory?

I could calculate the contact stiffness using the Hertz contact theory, see equations (5) and (6) in the attached report. This approach would give for the slope k of the pressure-overclosure curve

k = (2/pi) * Eeff / sqrt(A)

where Eeff is the effective elastic modulus and A the projected contact area. The term sqrt(A) would be replaced with the average element size h. For two steel surfaces in contact, this approach would result in k = 130 GPa / h.

Any comments?

Timo,

the problem with numerical contact is, as others have mentioned, the overclosure. Theory formulations consider a perfect geometry, but when dealing with a discretization you will have a node penetrating the opposite element. This is where the overclosure/contact formulation comes into play. Penalty methods are easier to program, but suffer from numerical instability: the ideal penalty is infinite, but then you don't get a solution (infinite times zero is indetermined), and if you relax the penalty value, you might not capture the mechanics you are looking for.

An approach that has worked for me in the past is to start with a somewhat relaxed penalty value (e.g. equal to the bulk modulus or a small integer fraction of that). Once you get convergence, increase the penalty incrementally and refine your results, starting from the last converged solution. This way you can increase the penalty to values that would otherwise lead to divergence.

Depending on the software you are using, you may consider using a different contact algorithm, e.g. Lagrange multipliers, which are a much more elegant formulation and do not suffer from this instability...

Hope this helps,

Fernando

Timo,

you can apply the Hertz theory. However, utilizing the penalty method you only have perfect contact if the penalty parameter is infinite. Lagrange multipliers are an alternative but you will get more degrees of freedoms. Also mixed methods are possible. Which method you have to use depends on the problem you want to solve.

I agree with Elizabeth and Brett, a trial and error method is required. A starting point for the trial and error method for estimating the penalty stiffness is discussed in Chapter 3, Chandrupatla and Beleguntu, Introduction to Finite Elements in Engg.

In the stiffness matrix that has been obtained, say K, the value that has maximum magnitude is taken, and the contact stiffness is assumed to be 10000 times of that.

So if the maximum stiffness in the stiffness matrix is obtained to be Kij, then

C = Kij * (10^4)

Reason for taking such a high value is that at the point of contact, there are two nodes, one belonging to the contacting member and one belonging to the target member. The displacements of the two nodes are related to each other through a multi point constraint equation. In order to effectively transmit the displacements of the nodes over the contacts, the stiffness between them needs to be high.

The authors i referred to above say that this choice has been found to be satisfactory on most of the computers.

However please note that this choice will also depend on the type of the contact between the parts being modelled.

It is unclear from the original question whether it is a theoretical contact mechanics problem (for example, a perfectly smooth elastic sylinder on perfect flat rigid or eleastic plane), or it is about real bodies in contact.

For the first problem, I agree with Elisabetta Zanetti and Fernando Cacho-Nerin.

For the second problem, one needs consider the surface roughness and friction. Surface roughness can be measured to various degrees of accuracy and included in a FE model (I did this in real contact problem). Friction should be considered as there is local slip even if the bodies are stationary. In the second problem, trail-and-error is still needed and the contact stiffness is largely related to the local roughness profile.

As one's experience in dealing with a particualr type of contact problems increases, it becomes easier to select the right penalty value.

It would seem reasonable to relate the pressure - overclosure function to the bulk modulus, but how?

just a comment

The answer about the value of penalty is of course bouncing between very high value causing ill-conditioning of the matrices and a small value leading to the large penetration.

There is a trial to start with a small value of penalty and then to correct it based on the overclosure (penetration) and the value of pressure (seems like exact as you are asking). The method is originated to R.Taylor as the large penetration method. It has been generalized by us in a series of publications - last one is

Article: 3D Frictionless contact problems with large load-steps based on the covariant description for higher order approximation

Ridvan Izi, Alexander Konyukhov, Karl Schweizerhof

[]

Engineering Structures 01/2013; 50:107-114.

can we use this approach for low velocity impact problems? I mean, is penalty method good approach to simulate low velocity impact problem?

The penalty method and augmented Lagrangian methods are used in the most FEM packages (ANSYS, LS-DYNA, ABAQUS) for contact simulations. It is independent of which simulation is involved high- , or low velocity impact.

But there is a clear preference which time integration scheme should be used during the simulation. Thus, for low- velocity is both implicit and explicit can be applied, however, explicit scheme is preferred for the high-velocity impact simulation due to the problem with convergence within the implicit scheme.

Dear Dr. Alexander Konyukhov 

I am modelling behavior of rock joints with contact mechanics view, so I want to have a basic penalty or lagrange contact friction code to develope to my aim.

So , I want to help me,please.

The following answer was found in the ansys-blog. Hope it is useful.

The penalty method introduces a force at the contact detection point(s) that has penetrated across the target surface with the express purpose of eliminating the penetration. This method uses very simple formulas:

Fc = kc*Dp

for contact detection points that penetrate across the target, and

Fc = zero

for open contact detection points. kc is the contact stiffness (also called penalty stiffness); it is a predetermined property of the contact element. Dp is the penetration at the contact element. Hence, the larger the penetration, the greater the calculated force. The challenge here is that the magnitude of the force necessary to prevent penetration is completely unknown beforehand. Obviously, the force needs to be large enough to push the contact surface back to the target surface and, thereby, eliminate unwanted penetration — but not so large that it pushes the contact completely off the target surface, causing error and instability. A positive aspect of the penalty method is that it is elegantly simple. The negative is that you end up with a finite amount of penetration at the end of the load step. Of course, this penetration is necessary for a contact force to be generated. It is important, therefore, that the contact stiffness be large enough so that the resulting penetration is negligibly small; but the contact stiffness cannot be so large as to cause instability and nonconvergence. This same strategy is used both in the opposite direction to prevent separation (with bonded or no-separation behaviors) and in the tangential direction to enforce frictional resistance and no-sliding behavior.

A stiffness relationship between two bodies must be established for contact to occur. Without a contact stiffness, bodies will pass through one another. The relationship is generated through an 'elastic spring' that is put between the two bodies, where the contact force is equal to the product of the contact stiffness (k) and the penetration (δ). The amount of penetration (δ), or incompatibility, between the two bodies is therefore dependent on the stiffness k. Ideally, there should be no penetration, but this implies that k =  , which will lead to numerical instabilities. The value of k that is used depends on the relative stiffness of the bodies in contact.

https://www.sharcnet.ca/Software/Ansys/17.0/en-us/help/ans_lsd/Hlp_L_contact6.html

Any MotionSolve users here ?! This solver is providing only Penalty contact and there is no such an option as Trial and Error

I am simulating a bolt pulling test. For which I have used linear pressure overclosure. I tried using 21000 and 30000 as the contact stiffness values, but it changes the resultant mises stress. How much is the accurate value I should use as contact stiffness. Material is steel.

A stiffness relationship between two bodies must be established for contact to occur. Without a contact stiffness, bodies will pass through one another. The relationship is generated through an 'elastic spring' that is put between the two bodies, where the contact force is equal to the product of the contact stiffness (k) and the penetration (δ). The amount of penetration (δ), or incompatibility, between the two bodies is therefore dependent on the stiffness k. Ideally, there should be no penetration, but this implies that k =  , which will lead to numerical instabilities. The value of k that is used depends on the relative stiffness of the bodies in contact.

Hello I have a question about penalty contact stiffness. I'm simulating a roll forming process for two layer of steel. I used penalty method and I need to add it . is there any hand book to find this values?