How to deal with penetration when modeling indentation test?

Please login with a confirmed email address before reporting spam

Hi,

I am working on the modeling of a elastic homogeous material indented by a rigid indenter, spherical or conical.

The indenter is model by rigid material, and the contact is defined through contact pair. The indenter and the elastic material is put together through 'form assembly'. The 'segregated augmented Lagrangian' method is employed, with preset penalty factor control tunned for speed. The prescribed displacement increases cubically (the model file of an isotropic elastic trail run is attached)

However, severe penetration occurs slightly away from the indenter tip (see the attached figure).

I have beening reading related threads in this forum, and tried

1)select the same boundaries as both source and destination in a contact pair as posted by @Henri in this thread. I switch the source and destination of a contact pair (p1) and define a new one (p2). Both p1 and p2 are selected as pairs in the definition of 'contact'. This causes error 'Rigid domain adjacent to contact destination boundary'.

2)reduce the increase step path. The step number is increased from 100 to 1000.

3)refine the mesh. I have tried to refine either side of the contact pair.

4)increase the scale factor of contact pressure as described by an example 'models.sme.cylinder_roller_contact' and mentioned by @John in this post. I increased this factor from the default 1e8 to 1e14. Only the error of lumped step decreases, but penetration still occurred.

I would really appreciate it if someone could shed light on this.

Best,

HC Liu



4 Replies Last Post Dec 1, 2023, 4:44 a.m. EST
Henrik Sönnerlind COMSOL Employee

Please login with a confirmed email address before reporting spam

Posted: 5 months ago Nov 27, 2023, 10:22 a.m. EST

First some comments to what you have tried:

  1. Is used when there is a risk of self-contact, that is a certain boundary comes into contact with itself. This is not the case here. The error is caused by the fact that it is required that at least one of the two contacting objects is elastic (and this should be the 'destination' side).
  2. Can be used if there is divergence in the iterations. Else, it should not be needed.
  3. Can be a good idea. In particular, you want to resolve the curvature of the indentor.
  4. Doing that is catastrophic. It essentially means that you remove all convergence checks on the contact pressure. The scale factor should be of the same order as the expected contact pressure. The default 1E8 (=100 MPa) is characteristic for steel-to-steel contact.

The augmented Lagrangion method that you are using is 'exact' in the sense that there will be no overlap for a converged solution. What you need to do is to make sure that you get more iterations. You can either decrease the solver tolerance or the scalings of the degrees of freedom.

Note that the 'penalty factor' in the Contact settings strongly affects the convergence rate of the contact iterations. It also affects the estimate of the error. A low penalty factor will give meany iterations, and also a risk that the convergence criteria look fulfilled while there is still an overlap.

Another hint: Using a rigid domain / rigid material for the indentor is overkill. You can just use a Prescribed Displacement on that domain.

-------------------
Henrik Sönnerlind
COMSOL
First some comments to what you have tried: 1. Is used when there is a risk of self-contact, that is a certain boundary comes into contact with itself. This is not the case here. The error is caused by the fact that it is required that at least one of the two contacting objects is elastic (and this should be the 'destination' side). 2. Can be used if there is divergence in the iterations. Else, it should not be needed. 3. Can be a good idea. In particular, you want to resolve the curvature of the indentor. 4. Doing that is catastrophic. It essentially means that you remove all convergence checks on the contact pressure. The scale factor should be of the same order as the expected contact pressure. The default 1E8 (=100 MPa) is characteristic for steel-to-steel contact. The augmented Lagrangion method that you are using is 'exact' in the sense that there will be no overlap *for a converged solution*. What you need to do is to make sure that you get more iterations. You can either decrease the solver tolerance or the scalings of the degrees of freedom. Note that the 'penalty factor' in the Contact settings strongly affects the convergence rate of the contact iterations. It also affects the estimate of the error. A low penalty factor will give meany iterations, and also a risk that the convergence criteria look fulfilled while there is still an overlap. Another hint: Using a rigid domain / rigid material for the indentor is overkill. You can just use a Prescribed Displacement on that domain.

Please login with a confirmed email address before reporting spam

Posted: 5 months ago Nov 27, 2023, 11:13 p.m. EST
Updated: 5 months ago Nov 28, 2023, 7:15 a.m. EST

First some comments to what you have tried:

Thank you so much @Henrik for your detailed explanation.

Your advice definitely make sense to me. So, I tried to follow it, but failed to get expected results. I am not sure if I did it right. I was wondering if you could give furher more guidance.

About the Comments,

The default 1E8 (=100 MPa) is characteristic for steel-to-steel contact.

Is this true only for the SI unit system or for all unit systems? I mean, does the meaning of this value keep consistant with the unit system?

I am used to reset the 'Unit System' to 'None' like in other FEM softwares. But I keep it in my mind that the unit system is 'mm for length, N for force, t for mass, MPa for pressure, ...'. Under this circumstance, does the default 1E8 represent 100MPa or 1E8 MPa?

(refine the mesh) Can be a good idea. In particular, you want to resolve the curvature of the indentor.

About the mesh refinement, I encountered another problem. Sometimes, I have to refine the mesh. Otherwise, the derived 'indentation modulus v.s. depth' curve isn't smooth.

But when the mesh is refined, the distribution of contact pressure becomes bumpy near the contact edge (see the figure below or the attached graph). The contact pressure increases and decreased suddently by a vary large extent. Could you give me a hint about this?

attached figure

Here are some trail run following your advice.

What you need to do is to make sure that you get more iterations. You can either decrease the solver tolerance or the scalings of the degrees of freedom.

I changed the 'Relative tolerance' from default 0.001 to 0.0001 under the subnode 'Study 5-->Solver Configurations-->Solution5 -->Stationary Solver 1 -->General'. And the 'Maximum number of refinements' from 15 to 30 under '...-->Stationary Solver1 --> Direct-->Error-->Iterative refinement'. The 'Maximum number of iterations' in the setting for 'segregated 1' is also change from 15 to 30.

The calculation time apprximately doubled. I can see from the Convergence Plot that the Error is below 1E-4 now, but the Graphics view still gives a penetration plot and no appreciable change can be observed.

A low penalty factor will give meany iterations, and also a risk that the convergence criteria look fulfilled while there is still an overlap.

I use Segrated Augmented Lagrangian Method. In previous simulation, I used the preset penalty factor control tunned for speed, since the indenter initially touches the structure through no more than 1 point and there is presumably no contact pressure.

I tried to manually tune the penalty factor multipiler. The default 1 makes no difference. So, I changed it to 10 (which is the maximun contact pressure gained from previous simulation). This caused divergence.

How to set the penalty factor control so that better performance than the preset parameters can be achieved? Could you shed more light on this?

Another hint: Using a rigid domain / rigid material for the indentor is overkill. You can just use a Prescribed Displacement on that domain.

I use Prescribed Displacement this time in substitution of the previous Rigid Material to model the rigid indenter. The displacement is prescribed as (u0r,u0z)=(0,-hi*pw1(par)). The solver configurations is reset to default before running the computation.

The deformation plot doesn't change. The penetration still occurs.

I have got loads of questions. Thank you so~ much for your time and patience.

Best, HC Liu

>First some comments to what you have tried: Thank you so much @Henrik for your detailed explanation. Your advice definitely make sense to me. So, I tried to follow it, but failed to get expected results. I am not sure if I did it right. **I was wondering if you could give furher more guidance**. **About the Comments**, > The default 1E8 (=100 MPa) is characteristic for steel-to-steel contact. **Is this true only for the SI unit system or for all unit systems**? I mean, does the meaning of this value keep consistant with the unit system? I am used to reset the 'Unit System' to 'None' like in other FEM softwares. But I keep it in my mind that the unit system is 'mm for length, N for force, t for mass, MPa for pressure, ...'. Under this circumstance, does the default 1E8 represent 100MPa or 1E8 MPa? > (refine the mesh) Can be a good idea. In particular, you want to resolve the curvature of the indentor. About the mesh refinement, I encountered another problem. Sometimes, I have to refine the mesh. Otherwise, the derived 'indentation modulus v.s. depth' curve isn't smooth. But when the mesh is refined, **the distribution of contact pressure becomes bumpy near the contact edge** (see the figure below or the attached graph). The contact pressure **increases and decreased suddently by a vary large extent**. Could you give me a hint about this? ![attached figure](https://img.chkaja.com/b132e40659e1dc27.png) **Here are some trail run following your advice**. > What you need to do is to make sure that you get more iterations. You can either decrease the solver tolerance or the scalings of the degrees of freedom. I changed the '**Relative tolerance**' from default 0.001 to 0.0001 under the subnode 'Study 5-->Solver Configurations-->Solution5 -->Stationary Solver 1 -->General'. And the '**Maximum number of refinements**' from 15 to 30 under '...-->Stationary Solver1 --> Direct-->Error-->Iterative refinement'. The '**Maximum number of iterations**' in the setting for 'segregated 1' is also change from 15 to 30. The calculation time apprximately doubled. I can see from the Convergence Plot that the Error is below 1E-4 now, but **the Graphics view still gives a penetration plot** and no appreciable change can be observed. > A low penalty factor will give meany iterations, and also a risk that the convergence criteria look fulfilled while there is still an overlap. I use **Segrated Augmented Lagrangian** Method. In previous simulation, I used the **preset** penalty factor control tunned for **speed**, since the indenter initially touches the structure through no more than 1 point and there is presumably no contact pressure. I tried to **manually tune** the **penalty factor multipiler**. The default 1 makes no difference. So, I changed it to 10 (which is the maximun contact pressure gained from previous simulation). This caused divergence. How to set the **penalty factor control** so that better performance than the **preset** parameters can be achieved? Could you shed more light on this? > Another hint: Using a rigid domain / rigid material for the indentor is overkill. You can just use a Prescribed Displacement on that domain. I use **Prescribed Displacement** this time in substitution of the previous **Rigid Material** to model the rigid indenter. The displacement is prescribed as (u0r,u0z)=(0,-hi\*pw1(par)). The solver configurations is reset to default before running the computation. The deformation plot doesn't change. The penetration still occurs. I have got loads of questions. **Thank you so~ much for your time and patience**. Best, HC Liu


Henrik Sönnerlind COMSOL Employee

Please login with a confirmed email address before reporting spam

Posted: 5 months ago Nov 29, 2023, 9:00 a.m. EST

Generally speaking, it is difficult to tell why you see an overlap. I would still suspect that you are looking at an insufficiently converged solution.

Variable scaling: If you switch off the unit handling, you are on your own. I would not do that unless necessary. In your case, the contact pressure scaling should then be in terms of MPa, so 10-100 is probably a good number. This could still be a source of your problem. Also, the displacement scale should then be a relevant number of mm. BTW: There is a built-in unit system with your favorite units: The MPa system.

The bumpy result: This is may be a plotting artifact, but it also possible that there is an unconverged solution. Start by switching off 'Refinement' and 'Smoothing' in the plot so that you can see the raw data. Zoom in on the region with the bump. The contact pressure is represented by a linear shape function over each element.

Penalty factor choice: A too large value will cause divergence. It is not easy to determine an optimal value, but it seems that the defaults are OK for you.

Iterative refinements: If iterative refinements are needed, that may indicate some problem in the physics. Iterative refinements are only needed when the stiffness matrix is ill-conditioned. Double-check boundary conditions and material data.

-------------------
Henrik Sönnerlind
COMSOL
Generally speaking, it is difficult to tell why you see an overlap. I would still suspect that you are looking at an insufficiently converged solution. Variable scaling: If you switch off the unit handling, you are on your own. I would not do that unless necessary. In your case, the contact pressure scaling should then be in terms of MPa, so 10-100 is probably a good number. This could still be a source of your problem. Also, the displacement scale should then be a relevant number of mm. BTW: There is a built-in unit system with your favorite units: The MPa system. The bumpy result: This is may be a plotting artifact, but it also possible that there is an unconverged solution. Start by switching off 'Refinement' and 'Smoothing' in the plot so that you can see the raw data. Zoom in on the region with the bump. The contact pressure is represented by a linear shape function over each element. Penalty factor choice: A too large value will cause divergence. It is not easy to determine an optimal value, but it seems that the defaults are OK for you. Iterative refinements: If iterative refinements are needed, that may indicate some problem in the physics. Iterative refinements are only needed when the stiffness matrix is ill-conditioned. Double-check boundary conditions and material data.

Please login with a confirmed email address before reporting spam

Posted: 5 months ago Dec 1, 2023, 4:44 a.m. EST
Updated: 5 months ago Dec 1, 2023, 4:45 a.m. EST

Generally speaking, it is difficult to tell why you see an overlap. I would still suspect that you are looking at an insufficiently converged solution.

Thank you @Henrik Sönnerlind for your reply.

I checked the overall results, e.g., the derivated indentation modulus. It matches the prediction of related theories, though there are penetration and bumpy contact pressure.

About the penetration, is it possible that the peneration is also a plotting artifact or that it doesn't matter so much? If so, I can live with that.

To be honest, the distribution of contact pressure concerns me more, because I might need to gain insights about the contact process through it.

Variable scaling: If you switch off the unit handling, you are on your own. I would not do that unless necessary. In your case, the contact pressure scaling should then be in terms of MPa, so 10-100 is probably a good number. This could still be a source of your problem. Also, the displacement scale should then be a relevant number of mm. BTW: There is a built-in unit system with your favorite units: The MPa system.

I'll leave this to the default.

The bumpy result: This is may be a plotting artifact, but it also possible that there is an unconverged solution. Start by switching off 'Refinement' and 'Smoothing' in the plot so that you can see the raw data. Zoom in on the region with the bump. The contact pressure is represented by a linear shape function over each element.

Not sure if I did it right. In the settings of contact pressure (comp1.solid.Tn) Line Graph, I chose No refinement for the Resolution and None for the Smoothing. But the plot dosen't change.

The zoomed in graph is as follows (or see the attached picture) Zoomed in of contact pressure

Note: the x-axis is r re-normalized by the element size, which can be regarded as the element id on the contact surface.

Penalty factor choice: A too large value will cause divergence. It is not easy to determine an optimal value, but it seems that the defaults are OK for you.

I understand that the penalty factor is hard to determine a priori without enough experience, so I'll leave this to the default.

Iterative refinements: If iterative refinements are needed, that may indicate some problem in the physics. Iterative refinements are only needed when the stiffness matrix is ill-conditioned. Double-check boundary conditions and material data.

I'll leave this to the default, too.

Thank you again for your detailed comments and explanation.

Best,

HCL

>Generally speaking, it is difficult to tell why you see an overlap. I would still suspect that you are looking at an insufficiently converged solution. Thank you @Henrik Sönnerlind for your reply. I checked the overall results, e.g., the derivated indentation modulus. It matches the prediction of related theories, though there are penetration and bumpy contact pressure. About the penetration, is it possible that the peneration is also a plotting artifact or that it doesn't matter so much? If so, I can live with that. To be honest, the distribution of contact pressure concerns me more, because I might need to gain insights about the contact process through it. > Variable scaling: If you switch off the unit handling, you are on your own. I would not do that unless necessary. In your case, the contact pressure scaling should then be in terms of MPa, so 10-100 is probably a good number. This could still be a source of your problem. Also, the displacement scale should then be a relevant number of mm. BTW: There is a built-in unit system with your favorite units: The MPa system. I'll leave this to the default. >The bumpy result: This is may be a plotting artifact, but it also possible that there is an unconverged solution. Start by switching off 'Refinement' and 'Smoothing' in the plot so that you can see the raw data. Zoom in on the region with the bump. The contact pressure is represented by a linear shape function over each element. Not sure if I did it right. In the settings of contact pressure (comp1.solid.Tn) Line Graph, I chose **No refinement** for the **Resolution** and **None** for the **Smoothing**. But the plot dosen't change. The zoomed in graph is as follows (or see the attached picture) ![Zoomed in of contact pressure](https://img.chkaja.com/3c8e039e0bcac06f.png) Note: the x-axis is r re-normalized by the element size, which can be regarded as the element id on the contact surface. >Penalty factor choice: A too large value will cause divergence. It is not easy to determine an optimal value, but it seems that the defaults are OK for you. I understand that the penalty factor is hard to determine a priori without enough experience, so I'll leave this to the default. >Iterative refinements: If iterative refinements are needed, that may indicate some problem in the physics. Iterative refinements are only needed when the stiffness matrix is ill-conditioned. Double-check boundary conditions and material data. I'll leave this to the default, too. Thank you again for your detailed comments and explanation. Best, HCL

Reply

Please read the discussion forum rules before posting.

Please log in to post a reply.

Note that while COMSOL employees may participate in the discussion forum, COMSOL® software users who are on-subscription should submit their questions via the Support Center for a more comprehensive response from the Technical Support team.