Project file on GitLab

## Motivation of this test case

The aim of this test is:

• to show that it is possible to apply Dirichlet BCs at a boundary that changes over the course of the simulation
• to give an advanced use-case of the Python BC. Here essentially an iterative procedure for a contact problem is implemented within the Python BC.

## Problem description

Two elastic spheres of same radius $R$ are brought into contact. The sphere centers are displaced towards each other by $w_0$, with increasing values in every load step. Due to symmetry reasons a flat circular contact area of radius $a$ forms.

The contact between the two spheres is modelled as a Dirichlet BC on a varying boundary. The exact boundary and Dirichlet values for the $y$ displacements are determined in a Python script. Compared to the sketch above the prescribed $y$ displacements correspond to $w_0/2$, because due to symmetry only half of the system (a section of the lower sphere) is simulated.

## Analytical solution

The derivation of the formulae below can be found, e.g., in this book (in German).

The radius of the contact area is given by $$\begin{equation} a = \sqrt{\frac{w_0 R}{2}} \end{equation}$$

The average pressure $\bar p$ over a the secant with distance $\xi$ to the center of the contact area (cf. vertical dashed line in the sketch above) is assumed to be $$\begin{equation} \bar p(\xi) = \kappa \sqrt{a^2 - \xi^2} \label{eq:bar-p} \end{equation}$$ with the prefactor $\kappa$ given by $$\begin{equation} \kappa = \frac{G}{R \cdot (1-\nu)} ,. \end{equation}$$

The total force $F$ exerted on the contact area is given by $$\begin{equation} F = \frac{8 a^3}{3\kappa} ,. \end{equation}$$

## Results

Average pressure $\bar{p}$:
Total force $F$:
The simulation results for contact radii and total force reproduce the analytical square root and cubic laws, respectively. For the average pressure $\bar p$ the analytical form of Eq. $(\ref{eq:bar-p})$ is reproduced, too. But for $\bar p$ there are rather strong deviations between the numerical and analytical results, which might be due to deviations in the contact radii $a$, due to insufficient mesh resolution or due to the chosen linear order of FEM ansatz functions.