High-Order Finite Element Methods for Computational Contact Mechanics

Project team leader Dipl.-Ing. Stefan Kollmannsberger MSc
Phd students Dipl.-Ing. David Franke
Principal investigators PD Dr.- Ing. Alexander Düster
Prof. Dr. rer.nat. Ernst Rank
Sponsorship Deutsche Forschungsgemeinschaft

Project description

Solving mechanical contact problems with the help of finite element methods is well established and state-of-the-art. However, due to its nonlinear nature, it is still computationally quite costly. Up to now commercial FE programs for contact problems are based on low order finite elements, using mostly linear or quadratic shape functions. In order to approximate the geometry and to increase the solution's accuracy, in the h-version the meshes are refined. Using a different set of shape functions introduced by Szabo and Babuska we make use of high-order (p-FEM) elements with hierarchic shape functions in order to reduce the computational effort. These high-order elements enable exponential convergence rates in energy norm, provided that an increase of polynomial degree is combined with a proper mesh design.

High-order elements:
Hierachic shape functions based on integrated Legendre polynomials are used as a finite element basis. In contrast to Lagrange polynomials, Legendre polynomials lead to a better conditioned stiffness matrix when it comes to higher polynomial Ansatz spaces. Applying the blending function method allows to describe the given (curved) geometry very accurately, even on coarse meshes. Due to the robustness of high-order elements against distortion, large deformations (and therefore high aspect ratios) do not pose a significant threat. Furthermore, it has been shown that high-order elements are robust with respect to locking effects.

Contact Discretization:
In order to treat mechanical contact problems, the virtual work equation has to be solved, where the sum of the internal and external energy has to be minimized over the two contacting bodies. The influence of frictionless contact is applied for by enhancing the virtual work equation with an extra term representing the contact related part. The chosen method to satisfy this enhanced equation by means of the Hertz-Signorini-Moreau conditions for frictionless contact is the penalty method.

Numerical Example:
The well known Herzian contact problem for 2D is used as a numerical example, where an infinitely long cylinder is pressed onto a rigid plate. This example was chosen because an analytic solution for the contact pressure in the contact zone is available.Only one quater of the total Model was used in the numerical solution. The bottom cylinder is supposed to be flat and rigid. The quater cylinder is loaded by a constant displacement load at the top, which is equivalent to a load of F=70000 N. The radius r=10 mm, poisson ratio is 0,3, and E=210000 N/mm²

Hertz_modell Hertz_num_modell

Hertzian contact problem. Numerical model for Hertzian contact problem.


Three different high order methods, namely the p-version, the hp-version, and the rp-version of the FEM, were used to for the simulation. The p-version has problems to cope with the fact the one element has a change in its boundary conditions (from contact to no contact). Better approximations are derived by using the hp-, or the rp-method. The meshes and pressure distributions in the contact zone are plotted in the following figures.

4elem_mesh 4elem_stress

Four element mesh used for the simulation. Contact pressure in the contact zone for non fitting mesh.


4elem_mesh_hp 4elem_stress_hp

Initial four element mesh used for the simulation Contact pressure in the contact zone using the hp-method and a polynomial degree of 8.

after four hp- refinement steps.


4elem_mesh_it 4elem_stress_rp

Four element mesh used for the simulation Contact pressure in the contact zone using the rp-method.

after rp-method is applied.

From this it can be observed that the contact pressure computed with the p-FEM approach exhibits oscillations. This is due to the fact that the exact solution contains a loss of regularity at the end of the contact zone which is located in the interior of the edge of an element. Since the corresponding point with reduced regularity in the exact solution can not be represented by a smooth polynomial, the approximation of the element tends to oscillate.

Results and Outlook:
By means of a numerical example it has been demonstrated that high-order finite elements can yield very accurate results also for problems of computational contact mechanics. To achieve high accuracy, adaptive methods should be applied to account for the point with reduced regularity at the end of the contact zone.



  • Franke, D.; Duester, A.; Rank, E.:
    The p-version of the FEM for computational contact mechanics. In: Proceedings of Applied Mechanics, Bremen, Germany 2008.
  • Szabo, B.; Düster, A.; Rank, E.:
    The p-version of the Finite Element Method. In E. Stein, R. de Borst, T.J.R. Hughes, editors: Encyclopedia of Computational Mechanics. Volume 1, Chapter 5, pp. 119-139. John Wiley & Sons, 2004.

Contact: Dipl.-Ing. David Franke