# Elastic response of a cylinder loaded by a Hertzian contact pressure and maintained in equilibrium by its inertia

Pierric Mora, Fabien Treyssède, Laurent Laguerre

► To cite this version:

Pierric Mora, Fabien Treyssède, Laurent Laguerre. Elastic response of a cylinder loaded by a Hertzian contact pressure and maintained in equilibrium by its inertia. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2024, 480 (2287), 10.1098/rspa.2023.0891. hal-04532296

HAL Id: hal-04532296

https://hal.science/hal-04532296

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Elastic response of a cylinder loaded by a Hertzian contact pressure and maintained in equilibrium by its inertia

Pierric Mora, Fabien Treyssède, Laurent Laguerre

► To cite this version:

Pierric Mora, Fabien Treyssède, Laurent Laguerre. Elastic response of a cylinder loaded by a Hertzian contact pressure and maintained in equilibrium by its inertia. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2024, 480 (2287), 10.1098/rspa.2023.0891. hal-04532296

HAL Id: hal-04532296

https://hal.science/hal-04532296v1

Submitted on 4 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

**Cite this article:** Mora P, Treyssède F, Laguerre L. 2024 Elastic response of a cylinder loaded by a Hertzian contact pressure and maintained in equilibrium by its inertia. *Proc. R. Soc. A 480:* 20230891.

https://doi.org/10.1098/rspa.2023.0891

Received: 4 December 2023

Accepted: 29 February 2024

**Subject Areas:**

mechanics

**Keywords:**

Hertzian contact, compliance, elastic cylinder, contact resonance, quasi-static approximation

**Author for correspondence:**

Pierric Mora

e-mail: pierric.mora@univ-eiffel.fr

# Elastic response of a cylinder loaded by a Hertzian contact pressure and maintained in equilibrium by its inertia

Pierric Mora, Fabien Treyssède and Laurent Laguerre

GeoEND, GERS-GeoEND, Univ Gustave Eiffel, IFSTTAR, F-44344

Bouguenais, France

<figure><img src="image_5.png" /></figure>

PM, 0000-0002-1013-6863

We derive the compliance of an elastic cylinder submitted to a line Hertzian contact. The cylinder is maintained in static equilibrium by bulk forces that are proportional to rigid body motions. Displacements are measured by setting integral gauges that amount to prescribing zero net linear and angular momentum, if the problem were to depend upon time. Various cases are covered, representing either infinitesimal or finite contact displacements, including partial slip. The developments are illustrated by revisiting a classical example in what could be called *The heavy cylinder on a vibrating soil*. The four contact resonances and forced response of the system are given in closed form in the quasi-static approximation, and compared against a reference numerical solution. The formulae can also be used as building blocks to assemble the compliance matrix of a system comprising several cylinders.

## 1. Introduction

Obtaining the compliance of a Hertzian line contact is a notoriously difficult task. For point contacts, i.e. between two non-conformal surfaces such as two spheres or two non-parallel cylinders, a comprehensive collection of formulae have been derived from the 1950s following Mindlin's milestone article [1], and are now gathered in textbooks [2]. The key property that allows for a general treatment of the point contact

© 2024 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.
<figure><img src="image_2.png" /></figure>

Research

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Check for
updates

**Cite this article:** Mora P, Treyssède F, Laguerre L. 2024 Elastic response of a cylinder loaded by a Hertzian contact pressure and maintained in equilibrium by its inertia. *Proc. R. Soc. A 480:* 20230891.

https://doi.org/10.1098/rspa.2023.0891

Received: 4 December 2023

Accepted: 29 February 2024

**Subject Areas:**

mechanics

**Keywords:**

Hertzian contact, compliance, elastic cylinder, contact resonance, quasi-static approximation

**Author for correspondence:**

Pierric Mora

e-mail: pierric.mora@univ-eiffel.fr

# Elastic response of a cylinder loaded by a Hertzian contact pressure and maintained in equilibrium by its inertia

Pierric Mora, Fabien Treyssède and Laurent Laguerre

GeoEND, GERS-GeoEND, Univ Gustave Eiffel, IFSTTAR, F-44344

Bouguenais, France

<figure><img src="image_5.png" /></figure>

PM, 0000-0002-1013-6863

We derive the compliance of an elastic cylinder submitted to a line Hertzian contact. The cylinder is maintained in static equilibrium by bulk forces that are proportional to rigid body motions. Displacements are measured by setting integral gauges that amount to prescribing zero net linear and angular momentum, if the problem were to depend upon time. Various cases are covered, representing either infinitesimal or finite contact displacements, including partial slip. The developments are illustrated by revisiting a classical example in what could be called *The heavy cylinder on a vibrating soil*. The four contact resonances and forced response of the system are given in closed form in the quasi-static approximation, and compared against a reference numerical solution. The formulae can also be used as building blocks to assemble the compliance matrix of a system comprising several cylinders.

## 1. Introduction

Obtaining the compliance of a Hertzian line contact is a notoriously difficult task. For point contacts, i.e. between two non-conformal surfaces such as two spheres or two non-parallel cylinders, a comprehensive collection of formulae have been derived from the 1950s following Mindlin's milestone article [1], and are now gathered in textbooks [2]. The key property that allows for a general treatment of the point contact

© 2024 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.
compliance problem is the fact that the strains and the elastic displacements decay to 0 as the distance to contact increases—respectively as $s^{-2}$ and $s^{-1}$ if we call s that distance. The problem is therefore purely local and comes down to a half-space problem: it is enough to know the local curvatures to derive the contact area caused by a normal force, and the contact displacement can be measured relatively to any ‘far point’.

Regarding line contacts, only part of the problem can be treated as local: the strains decay to 0 as $s^{-1}$, enabling a general solution for the stresses and the contact width as a function of the local curvatures. However, the ‘far point’ concept does not apply. This is because the elastic displacement does not tend to 0 far from the contact, but diverges as logs. As a consequence, there is no general solution, there are only special cases: the solution depends on the entire geometry of the cross section, on the number and location of multiple contacts, on the presence of other loads such as bulk loads and on the non-trivial choice of a reference to measure displacements—the ‘datum’ in [2], in this article, we speak of ‘gauges’. There is a large literature on line contacts that will not be reviewed here; we rather refer to [2] (in which line contacts are called non-Hertzian because of the ‘datum’ issue), to [3] for a historical perspective and to [4] for the closely related problem of line contacts between three-dimensional bodies.

In this article, we address the problem of a cylinder loaded by a single Hertzian contact, balanced with body forces that are proportional to rigid body motions. This problem arises as an *inner* problem when considering the overall motion and deformation of a cylinder put into a time-varying contact with an external body. We treat the four cases that result into a non-zero net force or moment, namely monopolar loads in the three directions and a dipolar normal load. Closed form solutions are given for a variety of traction distributions representing either small incremental, or finite contact displacements, including partial slip within the Coulomb model. An important aspect in our derivation is how the reference for displacements is set: instead of choosing arbitrarily a point—there is no symmetry, and therefore no reason to choose the axis of the cylinder—we introduce integral gauges. Indeed, these gauges naturally appear when addressing the dynamic behaviour of the cylinder: they translate into prescribing zero linear and angular momentum, and enable separating the total motion into what we call below a *skeleton* and an *inner* field.

The problem addressed can represent a heavy cylinder resting on a vibrating soil, as we call it in an example. It also appears in statics when modelling the contacts between peripheral and core wires of a wire strand (see [5,6])—where the uniform normal body force is caused by a tension applied to the cable and by the curvature of the wires, or with slightly misaligned cylinders pressed together [7]. As a piece of context, our motivation stems from modelling guided wave propagation in wire strands and wire ropes, where contacts are known to play an important role on wave phenomena [8,9], and where current models [10] would certainly benefit from closed form formulae to account for contacts. Let us also cite [11] that develops a semianalytical strategy to model Hertzian stick-slip contacts between cylindrical particles in a discrete element method approach, and [12] that simplifies the modelling of line contacts in wire strands by using the well-known compliance of the symmetrically loaded cylinder.

We start in §2 by defining the problems that aim to be addressed, then sketch in §3 a general methodology to solve them and motivate the quasi-static *inner* problems that are the core of this article. Section 4 is devoted to these *inner* problems. A first step in their resolution is to obtain the related Green's functions. Although the general solutions are known since Michell's work [13] and are tabulated in textbooks [14], and despite the literature on similar problems being very large ('the heavy disc', 'the rotating disc', see [15,16], or [17] for a review with a historical perspective), we did not manage to find a record of the solutions to our four cases including the displacement field—not mentioning the gauges. Therefore, we believe that it should be valuable to report in this article a complete, self-sustained set of formulae for both the stress and the displacement fields. We finally give examples in §5 by revisiting two classical cases: A first example revisits the well-known formula for the 'compression of a cylinder pressed by two others', while a second one revisits 'the heavy disc' by adding a dynamic behaviour.
compliance problem is the fact that the strains and the elastic displacements decay to 0 as the distance to contact increases—respectively as $s^{-2}$ and $s^{-1}$ if we call s that distance. The problem is therefore purely local and comes down to a half-space problem: it is enough to know the local curvatures to derive the contact area caused by a normal force, and the contact displacement can be measured relatively to any ‘far point’.

Regarding line contacts, only part of the problem can be treated as local: the strains decay to 0 as $s^{-1}$, enabling a general solution for the stresses and the contact width as a function of the local curvatures. However, the ‘far point’ concept does not apply. This is because the elastic displacement does not tend to 0 far from the contact, but diverges as logs. As a consequence, there is no general solution, there are only special cases: the solution depends on the entire geometry of the cross section, on the number and location of multiple contacts, on the presence of other loads such as bulk loads and on the non-trivial choice of a reference to measure displacements—the ‘datum’ in [2], in this article, we speak of ‘gauges’. There is a large literature on line contacts that will not be reviewed here; we rather refer to [2] (in which line contacts are called non-Hertzian because of the ‘datum’ issue), to [3] for a historical perspective and to [4] for the closely related problem of line contacts between three-dimensional bodies.

In this article, we address the problem of a cylinder loaded by a single Hertzian contact, balanced with body forces that are proportional to rigid body motions. This problem arises as an *inner* problem when considering the overall motion and deformation of a cylinder put into a time-varying contact with an external body. We treat the four cases that result into a non-zero net force or moment, namely monopolar loads in the three directions and a dipolar normal load. Closed form solutions are given for a variety of traction distributions representing either small incremental, or finite contact displacements, including partial slip within the Coulomb model. An important aspect in our derivation is how the reference for displacements is set: instead of choosing arbitrarily a point—there is no symmetry, and therefore no reason to choose the axis of the cylinder—we introduce integral gauges. Indeed, these gauges naturally appear when addressing the dynamic behaviour of the cylinder: they translate into prescribing zero linear and angular momentum, and enable separating the total motion into what we call below a *skeleton* and an *inner* field.

The problem addressed can represent a heavy cylinder resting on a vibrating soil, as we call it in an example. It also appears in statics when modelling the contacts between peripheral and core wires of a wire strand (see [5,6])—where the uniform normal body force is caused by a tension applied to the cable and by the curvature of the wires, or with slightly misaligned cylinders pressed together [7]. As a piece of context, our motivation stems from modelling guided wave propagation in wire strands and wire ropes, where contacts are known to play an important role on wave phenomena [8,9], and where current models [10] would certainly benefit from closed form formulae to account for contacts. Let us also cite [11] that develops a semianalytical strategy to model Hertzian stick-slip contacts between cylindrical particles in a discrete element method approach, and [12] that simplifies the modelling of line contacts in wire strands by using the well-known compliance of the symmetrically loaded cylinder.

We start in §2 by defining the problems that aim to be addressed, then sketch in §3 a general methodology to solve them and motivate the quasi-static *inner* problems that are the core of this article. Section 4 is devoted to these *inner* problems. A first step in their resolution is to obtain the related Green's functions. Although the general solutions are known since Michell's work [13] and are tabulated in textbooks [14], and despite the literature on similar problems being very large ('the heavy disc', 'the rotating disc', see [15,16], or [17] for a review with a historical perspective), we did not manage to find a record of the solutions to our four cases including the displacement field—not mentioning the gauges. Therefore, we believe that it should be valuable to report in this article a complete, self-sustained set of formulae for both the stress and the displacement fields. We finally give examples in §5 by revisiting two classical cases: A first example revisits the well-known formula for the 'compression of a cylinder pressed by two others', while a second one revisits 'the heavy disc' by adding a dynamic behaviour.
<figure><img src="image_6.png" /><figcaption>Figure 1. (a and b) Scheme of the problem and of the coordinate systems.</figcaption></figure>

## 2. Problem statement

### (a) Summary

#### (i) System

A long elastic cylinder of radius $R$, shear modulus $\mu$, Poisson's ratio $\nu$ and density $\rho$ is put in contact with an external body, which is such that the system remains invariant along the axis direction (z-axis). The characteristics of the contact—efforts, area—are allowed to change over time. The system is represented in figure 1 along with three sets of coordinates: cylinder-centred Cartesian $(x, y, z)$ and polar $(r, \theta, z)$, and contact-centred polar $(s, \beta, z)$.

#### (ii) Boundary condition for the displacement

Under this action, the cylinder deforms and a small contact strip of half-width $a \ll R$ arises, centred around the line whose coordinates in the $(xy)$ plane are $C = (0, -R)^T$. Anticipating that the elastic field will be decomposed into two parts, we refer to it as the total field: the displacement and stress within the cylinder are notated $\mathbf{u}_{\text{tot}}$ and $\sigma_{\text{tot}}$. The displacement satisfies the following boundary condition:

$$ \mathbf{u}_{\text{tot}} = \mathbf{d} \quad \text{along the contact: } r = R, \quad \frac{-a}{R} < \theta < \frac{a}{R}, \qquad (2.1) $$

with $\mathbf{d}$ a given displacement distribution, considered as known. Whether $\mathbf{d}$ is indeed explicitly known (prescribed) or not (e.g. the external body is part of a larger system) is problem-dependent and irrelevant at this stage.

Equation (2.1) translates into a traction distribution $\mathbf{q}$ applied on the contact strip. The symmetric part of $\mathbf{d}$ will give rise to a monopolar-like traction while the antisymmetric part will produce a dipolar-like traction.

#### (iii) Model restriction and degrees of freedom

We restrict to certain contact models by assuming that $\mathbf{q}$ is a weighted sum of a finite number of well known distributions. Because they are caused by smooth surfaces in contact, and because we work in the usual limit $a \ll R$, we will refer to them as Hertzian, although this qualification entails different meanings in the literature: it is often reserved to frictionless contacts, and sometimes excludes line contacts. The traction distributions considered in this article are depicted in figure 2—their expressions and what they physically represent is recalled in §2b. In a nutshell, both small incremental (full stick for normal, shear and tilt) and finite displacements (normal, or shear with partial slip) are covered, although the focus is set on the full stick case in the examples.
<figure><img src="image_6.png" /><figcaption>Figure 1. (a and b) Scheme of the problem and of the coordinate systems.</figcaption></figure>

## 2. Problem statement

### (a) Summary

#### (i) System

A long elastic cylinder of radius $R$, shear modulus $\mu$, Poisson's ratio $\nu$ and density $\rho$ is put in contact with an external body, which is such that the system remains invariant along the axis direction (z-axis). The characteristics of the contact—efforts, area—are allowed to change over time. The system is represented in figure 1 along with three sets of coordinates: cylinder-centred Cartesian $(x, y, z)$ and polar $(r, \theta, z)$, and contact-centred polar $(s, \beta, z)$.

#### (ii) Boundary condition for the displacement

Under this action, the cylinder deforms and a small contact strip of half-width $a \ll R$ arises, centred around the line whose coordinates in the $(xy)$ plane are $C = (0, -R)^T$. Anticipating that the elastic field will be decomposed into two parts, we refer to it as the *total field*: the displacement and stress within the cylinder are notated $\mathbf{u}_{\text{tot}}$ and $\sigma_{\text{tot}}$. The displacement satisfies the following boundary condition:

$$ \mathbf{u}_{\text{tot}} = \mathbf{d} \quad \text{along the contact: } r=R, \quad \frac{-a}{R} < \theta < \frac{a}{R}, \qquad (2.1) $$

with $\mathbf{d}$ a given displacement distribution, considered as known. Whether $\mathbf{d}$ is indeed explicitly known (prescribed) or not (e.g. the external body is part of a larger system) is problem-dependent and irrelevant at this stage.

Equation (2.1) translates into a traction distribution $\mathbf{q}$ applied on the contact strip. The symmetric part of $\mathbf{d}$ will give rise to a monopolar-like traction while the antisymmetric part will produce a dipolar-like traction.

#### (iii) Model restriction and degrees of freedom

We restrict to certain contact models by assuming that $\mathbf{q}$ is a weighted sum of a finite number of well known distributions. Because they are caused by smooth surfaces in contact, and because we work in the usual limit $a \ll R$, we will refer to them as Hertzian, although this qualification entails different meanings in the literature: it is often reserved to frictionless contacts, and sometimes excludes line contacts. The traction distributions considered in this article are depicted in figure 2—their expressions and what they physically represent is recalled in §2b. In a nutshell, both small incremental (full stick for normal, shear and tilt) and finite displacements (normal, or shear with partial slip) are covered, although the focus is set on the full stick case in the examples.
<figure><img src="image_7.png" /><figcaption>**Figure 2.** (a) Force distribution along a line Hertzian contact, representing: a small incremental displacement (b1) or tilt (b2), a finite normal displacement (b3), a finite transverse displacement with partial slip (b4).</figcaption></figure>

We write the traction as the sum of four components:

$$ \mathbf{q} = Q_x q_{Q_x}(x)\mathbf{e}_{x} + Q_y q_{Q_y}(x)\mathbf{e}_{y} + Q_z q_{Q_z}(x)\mathbf{e}_{z} + \frac{M_z}{R} q_{M_z}(x)\mathbf{e}_{y}, \quad (2.2) $$

where $q_{Q_i}(x)$ are monopolar-like distributions while $q_{M_z}(x)$ is a dipolar-like one, and where $\mathbf{e}_i$ is either unit vector of the standard basis ($i = x, y, z$). Note that in equation (2.2) we have chosen to omit $\mathbf{e}_x$ and $\mathbf{e}_z$ dipolar terms. This is because the net force and moment of such distributions is zero: although they do produce a local field, they bring no contribution to the overall cylinder motion and as such are not under interest here. Finally, let us introduce the average traction $\mathbf{Q}$ and moment $\mathbf{M}$. Assuming $q_{Q_i}$ and $q_{M_z}$ are properly normalized, they express as

$$ \mathbf{Q} = \oint_{r=R} \mathbf{q} \cdot \mathbf{R} \, d\theta = Q_x \mathbf{e}_x + Q_y \mathbf{e}_y + Q_z \mathbf{e}_z \quad (2.3a) $$

and

$$ \mathbf{M} = \oint_{r=R} \mathbf{r} \times \mathbf{q} \cdot \mathbf{R} \, d\theta = (RQ_x + M_z)\mathbf{e}_z. \quad (2.3b) $$

### (iv) Problem

The problem is to find the displacement field $\mathbf{u}_{tot}$ resulting from the load, possibly as a function of time for a dynamic problem, and possibly from other loads if the actual system is comprised of several bodies in contact.

### (b) Distributions of load

We here specify the distributions of load $q(x)$ for which explicit formulae for the contact compliance are given in §4. These distributions are well known (e.g. [2]). $q(x)$ is defined for $-a < x < a$, where the contact half-width $a \ll R$ is treated as an arbitrary parameter—its value eventually depends on the radii and elastic constants of both bodies and on the normal pressure, following Hertz's Law. In equation (2.2), we have introduced the notation $q_{Q_i}$ to refer to a monopolar-like load on the $i = x, y, z$ component, and $q_{M_z}$ for a dipolar-like load around the z axis. The formulae below are normalized to $\int_{-a}^{a} q_{Q_i} \, dx = 1$ and $\int_{-a}^{a} xq_{M_z} \, dx = R$ to be compatible with equations (2.2) and (2.3).

— Small incremental uniform displacement:

$$ q_{Q_i}(x) = \frac{1}{\pi a} \frac{1}{\sqrt{1 - x^2/a^2}}. \quad (2.4) $$

Equation (2.4) represents a prescribed displacement that is uniform along the contact, i.e. assuming full-stick (e.g. [2] §7.2). It is well suited to describe a small vibration around a
<figure><img src="image_7.png" /><figcaption>**Figure 2.** (a) Force distribution along a line Hertzian contact, representing: a small incremental displacement (b1) or tilt (b2), a finite normal displacement (b3), a finite transverse displacement with partial slip (b4).</figcaption></figure>

We write the traction as the sum of four components:

$$ \mathbf{q} = Q_x q_{Q_x}(x)\mathbf{e}_{x} + Q_y q_{Q_y}(x)\mathbf{e}_{y} + Q_z q_{Q_z}(x)\mathbf{e}_{z} + \frac{M_z}{R} q_{M_z}(x)\mathbf{e}_{y}, \quad (2.2) $$

where $q_{Q_i}(x)$ are monopolar-like distributions while $q_{M_z}(x)$ is a dipolar-like one, and where $\mathbf{e}_i$ is either unit vector of the standard basis ($i = x, y, z$). Note that in equation (2.2) we have chosen to omit $\mathbf{e}_x$ and $\mathbf{e}_z$ dipolar terms. This is because the net force and moment of such distributions is zero: although they do produce a local field, they bring no contribution to the overall cylinder motion and as such are not under interest here. Finally, let us introduce the average traction $\mathbf{Q}$ and moment $\mathbf{M}$. Assuming $q_{Q_i}$ and $q_{M_z}$ are properly normalized, they express as

$$ \mathbf{Q} = \oint_{r=R} \mathbf{q} \cdot \mathbf{R} \, d\theta = Q_x \mathbf{e}_x + Q_y \mathbf{e}_y + Q_z \mathbf{e}_z \quad (2.3a) $$

and

$$ \mathbf{M} = \oint_{r=R} \mathbf{r} \times \mathbf{q} \cdot \mathbf{R} \, d\theta = (RQ_x + M_z)\mathbf{e}_z. \quad (2.3b) $$

### (iv) Problem

The problem is to find the displacement field $\mathbf{u}_{tot}$ resulting from the load, possibly as a function of time for a dynamic problem, and possibly from other loads if the actual system is comprised of several bodies in contact.

### (b) Distributions of load

We here specify the distributions of load $q(x)$ for which explicit formulae for the contact compliance are given in §4. These distributions are well known (e.g. [2]). $q(x)$ is defined for $-a < x < a$, where the contact half-width $a \ll R$ is treated as an arbitrary parameter—its value eventually depends on the radii and elastic constants of both bodies and on the normal pressure, following Hertz's Law. In equation (2.2), we have introduced the notation $q_{Q_i}$ to refer to a monopolar-like load on the $i = x, y, z$ component, and $q_{M_z}$ for a dipolar-like load around the z axis. The formulae below are normalized to $\int_{-a}^{a} q_{Q_i} dx = 1$ and $\int_{-a}^{a} xq_{M_z} dx = R$ to be compatible with equations (2.2) and (2.3).

— Small incremental uniform displacement:

$$ q_{Q_i}(x) = \frac{1}{\pi a} \frac{1}{\sqrt{1 - x^2/a^2}}. \quad (2.4) $$

Equation (2.4) represents a prescribed displacement that is uniform along the contact, i.e. assuming full-stick (e.g. [2] §7.2). It is well suited to describe a small vibration around a
given static state, assuming that finite displacement effects are of second order: slip effects for the transverse components ($Q_xq_{Q_x}$ and $Q_zq_{Q_z}$) and the nonlinear dependence of the contact half-width upon the normal component $Q_yq_{Q_y}$.

— Small incremental uniform tilt

$$q_{M_z}(x) = \frac{2R}{\pi a^2} \frac{x/a}{\sqrt{1 - x^2/a^2}} \quad (2.5)$$

Equation (2.5) represents a small, shearless prescribed tilt $u_y(x) = \alpha x$, with $\alpha = R(\kappa + 1)/(2\pi \mu a^2)$, assuming full-stick. This distribution is usually called the inclined flat punch (e.g. [2], ch. 2.7).

— Finite normal displacement

$$q_{Q_y}(x) = \frac{2}{\pi a} \sqrt{\frac{1-x^2}{a^2}} \quad (2.6)$$

Equation (2.6) represents a non-uniform deformation along the contact caused by a normal pressure of finite magnitude $Q_y$, resulting for instance into a flat shaped contact zone of half-width $a(Q_y)$ for a contact between two identical cylinders (e.g. [2] §4.2).

— Finite transverse displacement, partial slip

$$q_{Q_{x,z}}(x) = \begin{cases} \frac{2}{\pi(a^2 - c^2)}(\sqrt{a^2 - x^2} - \sqrt{c^2 - x^2}) & \text{in the stick zone } |x| < c, \\ \frac{2}{\pi(a^2 - c^2)}\sqrt{a^2 - x^2} & \text{in the slip zone } c < |x| < a. \end{cases} \quad (2.7)$$

Equation (2.7) represents the distribution caused by a tangential traction of finite magnitude, assuming that the adhesion follows Coulomb's law (e.g. [2] §7.2 and [18,19]). Although this adhesive model is simple and widely used, there is a vast literature on friction modelling, for which [20] can be a good starting point. Here, the contact zone separates into a stick zone of half-width $c$, and a slip zone near the edges. By calling $f$ the coefficient of friction, $P = Q_y$ the finite normal pressure and $Q_{shear} = \sqrt{Q_x^2 + Q_z^2} < fP$ the total shear magnitude, $c$ is given by the relation $c = a\sqrt{1 - Q_{shear}/fP}$. Equation (2.7) further assumes that the normal traction is given by equation (2.6), implying no tilt from equation (2.5), or $M_z = 0$. A way to reduce this hypothesis is shown in [21].

## 3. General principle of resolution

The idea is to describe the motion of the cylinder as the superposition of an average, rigid body like motion—the skeleton field, adopting a naming used in [22]—and an elastic deformation—the inner field. The equations of motion are written in terms of a state vector containing only the motion of the skeleton (the outer problem), while independent equations are obtained for the inner field and solved once for all (the inner problem). Then, these inner solutions are converted into equivalent, possibly nonlinear spring constants and reported into the outer problem and the boundary condition (2.1). The outer problem is thereby reduced to one with a finite number of degrees of freedom and can eventually be solved.

### (a) Skeleton and inner fields

The total field is expanded as follows:

$$\mathbf{u}_{\text{tot}}(\mathbf{r}) = \mathbf{U} + \boldsymbol{\Theta} \times \mathbf{r} + \mathbf{u}(\mathbf{r}) \quad (3.1a)$$

and

$$\boldsymbol{\sigma}_{\text{tot}} = \boldsymbol{\Sigma} + \boldsymbol{\sigma}, \quad (3.1b)$$
given static state, assuming that finite displacement effects are of second order: slip effects for the transverse components ($Q_xq_{Q_x}$ and $Q_zq_{Q_z}$) and the nonlinear dependence of the contact half-width upon the normal component $Q_yq_{Q_y}$.

— Small incremental uniform tilt

$$q_{M_z}(x) = \frac{2R}{\pi a^2} \frac{x/a}{\sqrt{1 - x^2/a^2}} \quad (2.5)$$

Equation (2.5) represents a small, shearless prescribed tilt $u_y(x) = \alpha x$, with $\alpha = R(\kappa + 1)/(2\pi \mu a^2)$, assuming full-stick. This distribution is usually called the inclined flat punch (e.g. [2], ch. 2.7).

— Finite normal displacement

$$q_{Q_y}(x) = \frac{2}{\pi a} \sqrt{\frac{1-x^2}{a^2}} \quad (2.6)$$

Equation (2.6) represents a non-uniform deformation along the contact caused by a normal pressure of finite magnitude $Q_y$, resulting for instance into a flat shaped contact zone of half-width $a(Q_y)$ for a contact between two identical cylinders (e.g. [2] §4.2).

— Finite transverse displacement, partial slip

$$q_{Q_{x,z}}(x) = \begin{cases} \frac{2}{\pi(a^2 - c^2)}(\sqrt{a^2 - x^2} - \sqrt{c^2 - x^2}) & \text{in the stick zone } |x| < c, \\ \frac{2}{\pi(a^2 - c^2)}\sqrt{a^2 - x^2} & \text{in the slip zone } c < |x| < a. \end{cases} \quad (2.7)$$

Equation (2.7) represents the distribution caused by a tangential traction of finite magnitude, assuming that the adhesion follows Coulomb's law (e.g. [2] §7.2 and [18,19]). Although this adhesive model is simple and widely used, there is a vast literature on friction modelling, for which [20] can be a good starting point. Here, the contact zone separates into a stick zone of half-width $c$, and a slip zone near the edges. By calling $f$ the coefficient of friction, $P = Q_y$ the finite normal pressure and $Q_{shear} = \sqrt{Q_x^2 + Q_z^2} < fP$ the total shear magnitude, $c$ is given by the relation $c = a\sqrt{1 - Q_{shear}/fP}$. Equation (2.7) further assumes that the normal traction is given by equation (2.6), implying no tilt from equation (2.5), or $M_z = 0$. A way to reduce this hypothesis is shown in [21].

### 3. General principle of resolution

The idea is to describe the motion of the cylinder as the superposition of an average, rigid body like motion—the skeleton field, adopting a naming used in [22]—and an elastic deformation—the inner field. The equations of motion are written in terms of a state vector containing only the motion of the skeleton (the outer problem), while independent equations are obtained for the inner field and solved once for all (the inner problem). Then, these inner solutions are converted into equivalent, possibly nonlinear spring constants and reported into the outer problem and the boundary condition (2.1). The outer problem is thereby reduced to one with a finite number of degrees of freedom and can eventually be solved.

#### (a) Skeleton and inner fields

The total field is expanded as follows:

$$\mathbf{u}_{\text{tot}}(\mathbf{r}) = \mathbf{U} + \boldsymbol{\Theta} \times \mathbf{r} + \mathbf{u}(\mathbf{r}) \quad (3.1a)$$

and

$$\boldsymbol{\sigma}_{\text{tot}} = \boldsymbol{\Sigma} + \boldsymbol{\sigma}, \quad (3.1b)$$
— **r** = (x, y, z) is the position vector from the axis of the cylinder.

— **U** + **Θ** × **r** is the motion of the skeleton and will be called skeleton field: **U** is a uniform translation, and **Θ** × **r** = **Θ** **e**<sub>z</sub> × **r** is a uniform rotation.

— **u** is an elastic displacement relative to the skeleton and will be called inner field.

— **Σ** = **0** is the (null) skeleton stress field.

— **σ** is the inner stress field.

(b) Gauges for the inner field

Expansion (3.1a) is defined up to integration constants that represent rigid-body motions. Here, the inner displacement field is made unique by requiring zero average uniform translations and rotation about the principal axis

$$ \iint_S \rho \mathbf{u} \, dS = \mathbf{0} \qquad (3.2a) $$

and

$$ \iint_S \rho \mathbf{r} \times \mathbf{u} \, dS = \mathbf{0}, \qquad (3.2b) $$

in which $\iint_S dS$ refers to integrating over the cross section of the cylinder.

Because they come down to requiring zero linear and angular momentum for **u** in the limit of small displacements, these gauges are well suited to both static and dynamic problems.

(c) State vectors

In equation (2.2), we have written the traction as a weighted sum of four distributions after restricting to cases listed in §2b. Likewise, **d** and the boundary condition (2.1) can then also be characterized with four values. For this purpose, it is useful to define $d_\alpha|_C = \partial_x d_y|_C$ the prescribed inclination at contact, and

$$ \alpha|_C = \partial_x u_y|_C \qquad (3.3) $$

the inner field local inclination at contact.

Let us introduce the following state vectors to write the equations in a compact form:

— $\hat{\mathbf{d}}|_C = (d_x, d_y, d_z, R d_\alpha)^T|_C$: prescribed displacement at contact.

— $\hat{\mathbf{Q}} = (Q_x, Q_y, Q_z, R^{-1} M_z)^T$: net efforts at contact.

— $\hat{\mathbf{u}} = (u_x, u_y, u_z, R\alpha)^T$: inner field → $\hat{\mathbf{u}}|_C$: inner field at contact.

— $\hat{\mathbf{U}} = (U_x, U_y, U_z, R\Theta)^T$: skeleton field state variables.

By using the relations $\mathbf{u}_{tot}|_C = \mathbf{U} + R\Theta\mathbf{e}_x + \mathbf{u}|_C$ and $\partial_x u_{tot}|_C = \Theta + \alpha$, equation (2.1) can be rewritten

$$ \hat{\mathbf{R}}\hat{\mathbf{U}} + \hat{\mathbf{u}}|_C = \hat{\mathbf{d}}|_C, \qquad (3.4) $$

with $\hat{\mathbf{R}}$ the relations matrix

$$ \hat{\mathbf{R}} = \begin{bmatrix} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}. \qquad (3.5) $$
— **r** = (x, y, z) is the position vector from the axis of the cylinder.

— **U** + **Θ** × **r** is the motion of the skeleton and will be called skeleton field: **U** is a uniform translation, and **Θ** × **r** = **Θ** **e**<sub>z</sub> × **r** is a uniform rotation.

— **u** is an elastic displacement relative to the skeleton and will be called inner field.

— **Σ** = **0** is the (null) skeleton stress field.

— **σ** is the inner stress field.

(b) Gauges for the inner field

Expansion (3.1a) is defined up to integration constants that represent rigid-body motions. Here, the inner displacement field is made unique by requiring zero average uniform translations and rotation about the principal axis

$$ \iint_S \rho \mathbf{u} \, dS = \mathbf{0} \qquad (3.2a) $$

and

$$ \iint_S \rho \mathbf{r} \times \mathbf{u} \, dS = \mathbf{0}, \qquad (3.2b) $$

in which $\iint_S dS$ refers to integrating over the cross section of the cylinder.

Because they come down to requiring zero linear and angular momentum for **u** in the limit of small displacements, these gauges are well suited to both static and dynamic problems.

(c) State vectors

In equation (2.2), we have written the traction as a weighted sum of four distributions after restricting to cases listed in §2b. Likewise, **d** and the boundary condition (2.1) can then also be characterized with four values. For this purpose, it is useful to define $d_\alpha|_C = \partial_x d_y|_C$ the prescribed inclination at contact, and

$$ \alpha|_C = \partial_x u_y|_C \qquad (3.3) $$

the inner field local inclination at contact.

Let us introduce the following state vectors to write the equations in a compact form:

— $\hat{\mathbf{d}}|_C = (d_x, d_y, d_z, R d_\alpha)^T|_C$: prescribed displacement at contact.

— $\hat{\mathbf{Q}} = (Q_x, Q_y, Q_z, R^{-1} M_z)^T$: net efforts at contact.

— $\hat{\mathbf{u}} = (u_x, u_y, u_z, R\alpha)^T$: inner field → $\hat{\mathbf{u}}|_C$: inner field at contact.

— $\hat{\mathbf{U}} = (U_x, U_y, U_z, R\Theta)^T$: skeleton field state variables.

By using the relations $\mathbf{u}_{tot}|_C = \mathbf{U} + R\Theta\mathbf{e}_x + \mathbf{u}|_C$ and $\partial_x u_{tot}|_C = \Theta + \alpha$, equation (2.1) can be rewritten

$$ \hat{\mathbf{R}}\hat{\mathbf{U}} + \hat{\mathbf{u}}|_C = \hat{\mathbf{d}}|_C, \qquad (3.4) $$

with $\hat{\mathbf{R}}$ the relations matrix

$$ \hat{\mathbf{R}} = \begin{bmatrix} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}. \qquad (3.5) $$
(d) Dynamic inner and outer problems

(i) Balance of forces

The balance of forces for the *total* field reads

$$ \rho \partial_t^2 \mathbf{u}_{\text{tot}} - \nabla \cdot \boldsymbol{\sigma}_{\text{tot}} = \mathbf{0} \quad \text{for } r < R \qquad (3.6a) $$

and

$$ \boldsymbol{\sigma}_{\text{tot}} \cdot \mathbf{e}_r = \mathbf{q} \quad \text{for } r = R. \qquad (3.6b) $$

The next step is to separate equations (3.6) into equations involving only the *skeleton* on the one hand, and only the *inner* field on the other hand.

(ii) The outer problem

We apply the gauge integrals (3.2) to equations (3.6), use the divergence theorem to transform $\iint_S \nabla \cdot \boldsymbol{\sigma}_{\text{tot}} dS = \mathbf{Q}$, and $\iint_S \mathbf{r} \times \nabla \cdot \boldsymbol{\sigma}_{\text{tot}} dS = \mathbf{M}$ by noticing that $\mathbf{r} \times \nabla \cdot \boldsymbol{\sigma}_{\text{tot}} = \partial_i(\mathbf{r} \times \boldsymbol{\sigma}_{i,\text{tot}})$ because $\mathbf{e}_i \times \boldsymbol{\sigma}_{i,\text{tot}} = \mathbf{0}$. We then obtain the law of conservation of momentum

$$ \partial_t^2 \hat{\mathcal{M}} \hat{\mathbf{U}} = \hat{\mathcal{R}}^T \hat{\mathbf{Q}}, \qquad (3.7) $$

with

$$ \hat{\mathcal{M}} = m \operatorname{diag} \left( \left[ 1, 1, 1, \frac{J}{mR^2} \right] \right), \qquad (3.8) $$

the mass matrix, $m = \iint_S \rho dS = \rho\pi R^2$ the mass and $J = \iint_S \rho r^2 dS = (1/2) mR^2$ the moment of inertia with respect to the principal axis.

(iii) The inner problem

Inserting equations (3.7) into (3.6) leads to a balance of forces that does not explicitly involve the *skeleton*

$$ \rho \partial_t^2 \mathbf{u} - \nabla \cdot \boldsymbol{\sigma} = \mathbf{b} \quad \text{for } r < R, \qquad (3.9a) $$

$$ \boldsymbol{\sigma} \cdot \mathbf{e}_r = \mathbf{q} \quad \text{for } r = R \qquad (3.9b) $$

and

$$ \mathbf{b} = -\frac{\rho}{m}\mathbf{Q} - \frac{\rho}{J}\mathbf{M} \times \mathbf{r}. \qquad (3.9c) $$

In equations (3.9), the bulk load **b** balances the surface traction **q** on average, which ensures that the static limit can be taken, as done later on. **b** can be interpreted as the superposition of rigid-body inertia forces

$$ \mathbf{b} = \mathbf{b}^{\text{tr}} + \mathbf{b}^{\text{rot}}, \qquad (3.10) $$

where $\mathbf{b}^{\text{tr}} = -\rho m^{-1} \mathbf{Q}$ can be viewed as the inertia of a uniform translation and $\mathbf{b}^{\text{rot}} = -\rho J^{-1} \mathbf{M} \times \mathbf{r}$ can be viewed as the inertia of a uniform rotation.

(e) Quasi-static compliance of the system

Taking the static limit in equations (3.9) by neglecting $\rho \partial_t^2 \mathbf{u} \approx \mathbf{0}$ leads to a well-posed *inner* problem that can be solved analytically for the traction distributions listed earlier: this solution is given in §4. Formally, let us write

$$ \hat{\mathbf{u}}|_{\mathbf{C}} = \hat{\mathcal{C}}\hat{\mathbf{Q}}, \qquad (3.11) $$

with $\hat{\mathcal{C}}(a)$ the contact compliance matrix, which depends on the contact half-width $a$. For finite displacements (i.e. $q_{Q_y}$ is given by equation (2.6)), we recall that $a$ depends on the normal pressure $Q_y$, which renders the relation (3.11) nonlinear. For small incremental displacements ($q_{Q_y}$ is given by equation (2.4)), relation (3.11) is of course linear, and is defined around the static state given by $a$.
(d) Dynamic inner and outer problems

(i) Balance of forces

The balance of forces for the *total* field reads

$$ \rho \partial_t^2 \mathbf{u}_{\text{tot}} - \nabla \cdot \boldsymbol{\sigma}_{\text{tot}} = \mathbf{0} \quad \text{for } r < R \qquad (3.6a) $$

and

$$ \boldsymbol{\sigma}_{\text{tot}} \cdot \mathbf{e}_r = \mathbf{q} \quad \text{for } r = R. \qquad (3.6b) $$

The next step is to separate equations (3.6) into equations involving only the *skeleton* on the one hand, and only the *inner* field on the other hand.

(ii) The outer problem

We apply the gauge integrals (3.2) to equations (3.6), use the divergence theorem to transform $\iint_S \nabla \cdot \boldsymbol{\sigma}_{\text{tot}} dS = \mathbf{Q}$, and $\iint_S \mathbf{r} \times \nabla \cdot \boldsymbol{\sigma}_{\text{tot}} dS = \mathbf{M}$ by noticing that $\mathbf{r} \times \nabla \cdot \boldsymbol{\sigma}_{\text{tot}} = \partial_i(\mathbf{r} \times \boldsymbol{\sigma}_{i,\text{tot}})$ because $\mathbf{e}_i \times \boldsymbol{\sigma}_{i,\text{tot}} = \mathbf{0}$. We then obtain the law of conservation of momentum

$$ \partial_t^2 \hat{\mathcal{M}} \hat{\mathbf{U}} = \hat{\mathcal{R}}^T \hat{\mathbf{Q}}, \qquad (3.7) $$

with

$$ \hat{\mathcal{M}} = m \operatorname{diag} \left( \left[ 1, 1, 1, \frac{J}{mR^2} \right] \right), \qquad (3.8) $$

the mass matrix, $m = \iint_S \rho dS = \rho\pi R^2$ the mass and $J = \iint_S \rho r^2 dS = (1/2) mR^2$ the moment of inertia with respect to the principal axis.

(iii) The inner problem

Inserting equations (3.7) into (3.6) leads to a balance of forces that does not explicitly involve the *skeleton*

$$ \rho \partial_t^2 \mathbf{u} - \nabla \cdot \boldsymbol{\sigma} = \mathbf{b} \quad \text{for } r < R, \qquad (3.9a) $$

$$ \boldsymbol{\sigma} \cdot \mathbf{e}_r = \mathbf{q} \quad \text{for } r = R \qquad (3.9b) $$

and

$$ \mathbf{b} = -\frac{\rho}{m}\mathbf{Q} - \frac{\rho}{J}\mathbf{M} \times \mathbf{r}. \qquad (3.9c) $$

In equations (3.9), the bulk load **b** balances the surface traction **q** on average, which ensures that the static limit can be taken, as done later on. **b** can be interpreted as the superposition of rigid-body inertia forces

$$ \mathbf{b} = \mathbf{b}^{\text{tr}} + \mathbf{b}^{\text{rot}}, \qquad (3.10) $$

where $\mathbf{b}^{\text{tr}} = -\rho m^{-1} \mathbf{Q}$ can be viewed as the inertia of a uniform translation and $\mathbf{b}^{\text{rot}} = -\rho J^{-1} \mathbf{M} \times \mathbf{r}$ can be viewed as the inertia of a uniform rotation.

(e) Quasi-static compliance of the system

Taking the static limit in equations (3.9) by neglecting $\rho \partial_t^2 \mathbf{u} \approx \mathbf{0}$ leads to a well-posed *inner* problem that can be solved analytically for the traction distributions listed earlier: this solution is given in §4. Formally, let us write

$$ \hat{\mathbf{u}}|_{\mathbf{C}} = \hat{\mathcal{C}}\hat{\mathbf{Q}}, \qquad (3.11) $$

with $\hat{\mathcal{C}}(a)$ the contact compliance matrix, which depends on the contact half-width $a$. For finite displacements (i.e. $q_{Q_y}$ is given by equation (2.6)), we recall that $a$ depends on the normal pressure $Q_y$, which renders the relation (3.11) nonlinear. For small incremental displacements ($q_{Q_y}$ is given by equation (2.4)), relation (3.11) is of course linear, and is defined around the static state given by $a$.
Anticipating the solutions given in the next section, $\hat{\mathcal{C}}$ can be approximated by a diagonal matrix because the coordinates $(x, y)$ align with the tangential and normal directions

$$ \hat{\mathcal{C}} = \mathrm{diag}([C_{xx}, C_{yy}, C_{zz}, C_{\alpha\alpha}]). \qquad (3.12) $$

Indeed, the off-diagonal terms in equation (3.12) $C_{\alpha\alpha} = R u_x |_{C} / M_z$ and $C_{\alpha x} = R \alpha |_{C} / Q_x$ are much smaller than the corresponding diagonal terms and can be neglected: $C_{\alpha\alpha} \ll C_{xx}$ and $C_{\alpha x} \ll C_{\alpha\alpha}$.

### (f) Extension to several contact points: far field compliance

If the overall system is comprised of more than two cylinders, such as one of these undergoes several contacts, then the global compliance matrix will need other terms to be assembled. Indeed, by calling $C_2$ another contact point to the cylinder under consideration, the inner field at that point is also composed of a far field contribution from the efforts at point C—renamed $C_1$ here. Let us then formally extend equation (3.11) as

$$ \hat{\mathbf{u}}|_{C_1} = \hat{\mathcal{C}}(C_1, a_1)\hat{\mathbf{Q}}_1 + \hat{\mathcal{G}}(C_2, R, \theta_{2,1})\hat{\mathbf{Q}}_2 \qquad (3.13a) $$

and

$$ \hat{\mathbf{u}}|_{C_2} = \hat{\mathcal{G}}(C_1, R, \theta_{1,2})\hat{\mathbf{Q}}_1 + \hat{\mathcal{C}}(C_2, a_2)\hat{\mathbf{Q}}_2, \qquad (3.13b) $$

where we have defined $\hat{\mathcal{G}}(C_i, r, \theta)$ the far-field compliance matrix from a source at point $C_i$ ($i=1,2$), $\hat{\mathbf{u}}|_{C_i}(\hat{\mathbf{Q}}_i)$ the inner displacement (efforts) at point $C_i$, $\theta_{1,2}$ the angle from $C_1$ to $C_2$ and $a_i$ the $i$th contact half-width. Note that we have added the $C_i$ argument to $\hat{\mathcal{C}}$ to emphasize that it might in general not be diagonal—to preserve this property one should use polar coordinates for $\hat{\mathbf{u}}_i$ and $\hat{\mathbf{Q}}_i$, which also leads to a more symmetric form for $\hat{\mathcal{G}}$.

Unlike $\hat{\mathcal{C}}, \hat{\mathcal{G}}$ does not depend on $a_i$, however, it contains off-diagonal terms and can be approximated by

$$ \hat{\mathcal{G}} = \begin{bmatrix} \mathcal{G}_{xx} & \mathcal{G}_{xy} & 0 & 0 \\ \mathcal{G}_{yx} & \mathcal{G}_{yy} & 0 & 0 \\ 0 & 0 & \mathcal{G}_{zz} & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}. \qquad (3.14) $$

Here, we have neglected any far field term involving $\alpha$ because of its qualitatively faster decrease with the distance to the source (as $1/||r - C||$) than the other terms (as $\log ||r - C||$).

### (g) Resolution

Equations (3.4), (3.7) and (3.11) (or (3.13)) represent a well-posed system, which can be solved. Let us emphasize again that an essential brick is equations (3.11) and (3.13), which is the object of §4.

## 4. The static inner problems and their solutions

In this section, we treat the inner problem defined by equations (3.9) and (3.2) as a set of four independent problems. We seek to obtain the displacement field $\mathbf{u}$. We solve these four problems analytically in the static limit, i.e. by taking $\rho \partial_t^2 \mathbf{u} \approx 0$.

We first define four Green's problems by considering line loads of unit amplitude ($Q_{x,y,z} = 1$ or $0$, $R^{-1}M_z = 1$ or $0$) and call $\hat{\mathcal{G}}$ their solutions (Green's matrix). We then consider normalized distributions of finite extent (see §2b) and call $\hat{\mathcal{C}}$ the matrix of solutions at the centre of the contact, i.e. for $\mathbf{r} = C$.

**Plane strain:** We assume that the cylinder is long enough for the plane strain hypothesis to be valid. Solutions are expressed using Kolosov's<fnref n="1" /> constant $\kappa = 3 - 4\nu$.

<footnote n="1">¹The formulae are also valid for plane stress by taking $\kappa = (3 - \nu)/(1 + \nu)$.</footnote>
Anticipating the solutions given in the next section, $\hat{\mathcal{C}}$ can be approximated by a diagonal matrix because the coordinates $(x, y)$ align with the tangential and normal directions

$$ \hat{\mathcal{C}} = \mathrm{diag}([C_{xx}, C_{yy}, C_{zz}, C_{\alpha\alpha}]). \qquad (3.12) $$

Indeed, the off-diagonal terms in equation (3.12) $C_{\alpha\alpha} = R u_x |_{C} / M_z$ and $C_{\alpha x} = R \alpha |_{C} / Q_x$ are much smaller than the corresponding diagonal terms and can be neglected: $C_{\alpha\alpha} \ll C_{xx}$ and $C_{\alpha x} \ll C_{\alpha\alpha}$.

### (f) Extension to several contact points: far field compliance

If the overall system is comprised of more than two cylinders, such as one of these undergoes several contacts, then the global compliance matrix will need other terms to be assembled. Indeed, by calling $C_2$ another contact point to the cylinder under consideration, the inner field at that point is also composed of a far field contribution from the efforts at point C—renamed $C_1$ here. Let us then formally extend equation (3.11) as

$$ \hat{\mathbf{u}}|_{C_1} = \hat{\mathcal{C}}(C_1, a_1)\hat{\mathbf{Q}}_1 + \hat{\mathcal{G}}(C_2, R, \theta_{2,1})\hat{\mathbf{Q}}_2 \qquad (3.13a) $$

and

$$ \hat{\mathbf{u}}|_{C_2} = \hat{\mathcal{G}}(C_1, R, \theta_{1,2})\hat{\mathbf{Q}}_1 + \hat{\mathcal{C}}(C_2, a_2)\hat{\mathbf{Q}}_2, \qquad (3.13b) $$

where we have defined $\hat{\mathcal{G}}(C_i, r, \theta)$ the far-field compliance matrix from a source at point $C_i$ ($i=1,2$), $\hat{\mathbf{u}}|_{C_i}(\hat{\mathbf{Q}}_i)$ the inner displacement (efforts) at point $C_i$, $\theta_{1,2}$ the angle from $C_1$ to $C_2$ and $a_i$ the $i$th contact half-width. Note that we have added the $C_i$ argument to $\hat{\mathcal{C}}$ to emphasize that it might in general not be diagonal—to preserve this property one should use polar coordinates for $\hat{\mathbf{u}}_i$ and $\hat{\mathbf{Q}}_i$, which also leads to a more symmetric form for $\hat{\mathcal{G}}$.

Unlike $\hat{\mathcal{C}}, \hat{\mathcal{G}}$ does not depend on $a_i$, however, it contains off-diagonal terms and can be approximated by

$$ \hat{\mathcal{G}} = \begin{bmatrix} \mathcal{G}_{xx} & \mathcal{G}_{xy} & 0 & 0 \\ \mathcal{G}_{yx} & \mathcal{G}_{yy} & 0 & 0 \\ 0 & 0 & \mathcal{G}_{zz} & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}. \qquad (3.14) $$

Here, we have neglected any far field term involving $\alpha$ because of its qualitatively faster decrease with the distance to the source (as $1/||r - C||$) than the other terms (as $\log ||r - C||$).

### (g) Resolution

Equations (3.4), (3.7) and (3.11) (or (3.13)) represent a well-posed system, which can be solved. Let us emphasize again that an essential brick is equations (3.11) and (3.13), which is the object of §4.

## 4. The static inner problems and their solutions

In this section, we treat the inner problem defined by equations (3.9) and (3.2) as a set of four independent problems. We seek to obtain the displacement field $\mathbf{u}$. We solve these four problems analytically in the static limit, i.e. by taking $\rho \partial_t^2 \mathbf{u} \approx 0$.

We first define four Green's problems by considering line loads of unit amplitude ($Q_{x,y,z} = 1$ or $0$, $R^{-1}M_z = 1$ or $0$) and call $\hat{\mathcal{G}}$ their solutions (Green's matrix). We then consider normalized distributions of finite extent (see §2b) and call $\hat{\mathcal{C}}$ the matrix of solutions at the centre of the contact, i.e. for $\mathbf{r} = C$.

**Plane strain:** We assume that the cylinder is long enough for the plane strain hypothesis to be valid. Solutions are expressed using Kolosov's<fnref n="1" /> constant $\kappa = 3 - 4\nu$.

<footnote n="1">¹The formulae are also valid for plane stress by taking $\kappa = (3 - \nu)/(1 + \nu)$.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Load distribution in Green's problems. (a) Normal load, (b) in-plane tangential load, (c) out-of-plane tangential load and (d) dipolar normal load.</figcaption></figure>

(a) The four Green's problems and their solutions

As a first step of resolution, let us address the following four cases of unit line loads, represented in figure 3.

**Definitions.**

(a) normal load: $\mathbf{q} = \delta(\mathbf{r} - \mathbf{C})\mathbf{e}_y$ and $\mathbf{b} = -(\pi R^2)^{-1}\mathbf{e}_y$,

(b) in-plane tangential load: $\mathbf{q} = \delta(\mathbf{r} - \mathbf{C})\mathbf{e}_x$ and $\mathbf{b} = -(\pi R^2)^{-1}(\mathbf{e}_x + 2rR^{-1}\mathbf{e}_\theta)$,

(c) out-of-plane tangential load: $\mathbf{q} = \delta(\mathbf{r} - \mathbf{C})\mathbf{e}_z$ and $\mathbf{b} = -(\pi R^2)^{-1}\mathbf{e}_z$,

(d) dipolar normal load: $\mathbf{q} = R\partial_x\delta(\mathbf{r} - \mathbf{C})\mathbf{e}_y$ and $\mathbf{b} = -2r/(\pi R^3)\mathbf{e}_\theta$,

where $\delta$ is the unit impulse symbol.

(i) Solutions on the edge

We recall in appendix A the well-known methodology to obtain the solutions to these four problems, and provide the expression of Green's functions at any point inside the cylinder. Here, we report only their expression at $r=R$, where simplifications can be done using $\beta = (\theta - \pi)/2$ and $s = 2|\sin\theta/2|$ for $0 \le \theta \le 2\pi$. Indeed, these are the formulae that are required to assemble the compliance matrix of systems of several cylinders in contact as illustrated by the example in §5a.

(a) Normal load

$$4\pi \mu \mathcal{G}_{xy}|_{r=R} = \frac{1}{2}(\kappa - 1)(\theta - \pi) + \frac{1}{2}(\kappa + 1 + 2\cos\theta)\sin\theta \quad (4.1a)$$

and

$$4\pi \mu \mathcal{G}_{yy}|_{r=R} = \left( \frac{1}{2} - (\kappa + 1) \log \left| 2 \sin \frac{\theta}{2} \right| \right) - \frac{1}{2}(\kappa + 1 + 2\cos\theta)\cos\theta. \quad (4.1b)$$

(b) In-plane tangential load

$$4\pi \mu \mathcal{G}_{xx}|_{r=R} = \left( \frac{1}{2} - (\kappa + 1) \log \left| 2 \sin \frac{\theta}{2} \right| \right) - \left( \kappa - \frac{1}{3} - \cos\theta \right) \cos\theta \quad (4.2a)$$

and

$$4\pi \mu \mathcal{G}_{yx}|_{r=R} = -\frac{1}{2}(\kappa - 1)(\theta - \pi) - \left(\kappa - \frac{1}{3} - \cos\theta\right) \sin\theta. \quad (4.2b)$$
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Load distribution in Green's problems. (a) Normal load, (b) in-plane tangential load, (c) out-of-plane tangential load and (d) dipolar normal load.</figcaption></figure>

(a) The four Green's problems and their solutions

As a first step of resolution, let us address the following four cases of unit line loads, represented in figure 3.

**Definitions.**

(a) normal load: $\mathbf{q} = \delta(\mathbf{r} - \mathbf{C})\mathbf{e}_y$ and $\mathbf{b} = -(\pi R^2)^{-1}\mathbf{e}_y$,

(b) in-plane tangential load: $\mathbf{q} = \delta(\mathbf{r} - \mathbf{C})\mathbf{e}_x$ and $\mathbf{b} = -(\pi R^2)^{-1}(\mathbf{e}_x + 2rR^{-1}\mathbf{e}_\theta)$,

(c) out-of-plane tangential load: $\mathbf{q} = \delta(\mathbf{r} - \mathbf{C})\mathbf{e}_z$ and $\mathbf{b} = -(\pi R^2)^{-1}\mathbf{e}_z$,

(d) dipolar normal load: $\mathbf{q} = R\partial_x\delta(\mathbf{r} - \mathbf{C})\mathbf{e}_y$ and $\mathbf{b} = -2r/(\pi R^3)\mathbf{e}_\theta$,

where $\delta$ is the unit impulse symbol.

(i) Solutions on the edge

We recall in appendix A the well-known methodology to obtain the solutions to these four problems, and provide the expression of Green's functions at any point inside the cylinder. Here, we report only their expression at $r=R$, where simplifications can be done using $\beta = (\theta - \pi)/2$ and $s = 2|\sin\theta/2|$ for $0 \le \theta \le 2\pi$. Indeed, these are the formulae that are required to assemble the compliance matrix of systems of several cylinders in contact as illustrated by the example in §5a.

(a) Normal load

$$4\pi \mu \mathcal{G}_{xy}|_{r=R} = \frac{1}{2}(\kappa - 1)(\theta - \pi) + \frac{1}{2}(\kappa + 1 + 2\cos\theta)\sin\theta \quad (4.1a)$$

and

$$4\pi \mu \mathcal{G}_{yy}|_{r=R} = \left( \frac{1}{2} - (\kappa + 1) \log \left| 2 \sin \frac{\theta}{2} \right| \right) - \frac{1}{2}(\kappa + 1 + 2\cos\theta)\cos\theta. \quad (4.1b)$$

(b) In-plane tangential load

$$4\pi \mu \mathcal{G}_{xx}|_{r=R} = \left( \frac{1}{2} - (\kappa + 1) \log \left| 2 \sin \frac{\theta}{2} \right| \right) - \left( \kappa - \frac{1}{3} - \cos\theta \right) \cos\theta \quad (4.2a)$$

and

$$4\pi \mu \mathcal{G}_{yx}|_{r=R} = -\frac{1}{2}(\kappa - 1)(\theta - \pi) - \left(\kappa - \frac{1}{3} - \cos\theta\right) \sin\theta. \quad (4.2b)$$
(c) Out-of-plane tangential load

$$4\pi \mu \mathcal{G}_{zz}|_{r=R} = -4 \log \left| 2 \sin \frac{\theta}{2} \right| + \frac{1}{2}. \quad (4.3)$$

(d) Dipolar normal load

$$4\pi \mu \mathcal{G}_{x\alpha}|_{r=R} = -\frac{1}{2}(\kappa - 3 - 2\cos\theta) - \left(\kappa + \frac{2}{3}\right)\cos\theta \quad (4.4a)$$

and

$$4\pi \mu \mathcal{G}_{y\alpha}|_{r=R} = \frac{1}{2} \cot \frac{\theta}{2} (\kappa + 3 - 2 \cos \theta) - \left(\kappa + \frac{2}{3}\right) \sin \theta. \quad (4.4b)$$

Let us emphasize that equations (4.1)–(4.4) are valid only for $0 \le \theta \le 2\pi$ because of the relation $\beta(\theta)$. For another range of definition of $\theta$ such as $-\pi \le \theta \le \pi$, they have to be adapted by replacing $\theta \to \theta \mod 2\pi$, where mod is the modulus operator.

(b) The Hertzian contact problems and their solutions at the contact point

Let us now consider that the traction distribution $q$ is composed of normalized distributions of finite extent. As in the previous section, we set $Q_{x,y,z} = 1$ or $0$, $R^{-1}M_z = 1$ or $0$.

(i) Method of resolution

Green's functions given in appendix A are convolved with the normalized distributions defined in §2b, yielding the contact compliance matrix $\hat{C}$ when the convolutions are evaluated at point C. The integrals involved are solved analytically in the limit $2a \ll R$ – for the sake of conciseness they are given in appendix Cc.

(ii) Solutions for full-stick, small incremental displacements

We start by giving the compliance corresponding to distributions given by equations (2.4) and (2.5). Only non-zero terms are reported.

(a) Small, incremental normal load

$$4\pi \mu C_{yy} = (\kappa + 1) \log \frac{2R}{a\sqrt{e}} - \frac{1}{2}. \quad (4.5)$$

(b) Small, incremental in-plane tangential load

$$4\pi \mu C_{xx} = (\kappa + 1) \log \frac{2R}{ae} + 3 - \frac{1}{6} \quad (4.6a)$$

and

$$4\pi \mu C_{\alpha x} = \frac{R}{2a}\pi(\kappa - 1). \quad (4.6b)$$

(c) Small, incremental out-of-plane tangential load

$$4\pi \mu C_{zz} = 4 \log \frac{2R}{a} + \frac{1}{2}. \quad (4.7)$$

(d) Small, incremental dipolar normal load

$$4\pi \mu C_{\alpha\alpha} = \frac{2R^2}{a^2}(\kappa + 1) \quad (4.8a)$$

and

$$4\pi \mu C_{x\alpha} = -\frac{3}{2}\kappa + \frac{11}{6}, \quad (4.8b)$$

with $e = \exp(1)$. An application of the normal contact compliance formula (4.5) can be found in [6,23], where contact effects are accounted for in the stiffness matrix of a wire strand.
(c) Out-of-plane tangential load

$$4\pi \mu \mathcal{G}_{zz}|_{r=R} = -4 \log \left| 2 \sin \frac{\theta}{2} \right| + \frac{1}{2}. \qquad (4.3)$$

(d) Dipolar normal load

$$4\pi \mu \mathcal{G}_{x\alpha}|_{r=R} = -\frac{1}{2}(\kappa - 3 - 2\cos\theta) - \left(\kappa + \frac{2}{3}\right)\cos\theta \qquad (4.4a)$$

and

$$4\pi \mu \mathcal{G}_{y\alpha}|_{r=R} = \frac{1}{2} \cot \frac{\theta}{2} (\kappa + 3 - 2 \cos \theta) - \left(\kappa + \frac{2}{3}\right) \sin \theta. \qquad (4.4b)$$

Let us emphasize that equations (4.1)–(4.4) are valid only for $0 \le \theta \le 2\pi$ because of the relation $\beta(\theta)$. For another range of definition of $\theta$ such as $-\pi \le \theta \le \pi$, they have to be adapted by replacing $\theta \to \theta$ mod $2\pi$, where mod is the modulus operator.

(b) The Hertzian contact problems and their solutions at the contact point

Let us now consider that the traction distribution $q$ is composed of normalized distributions of finite extent. As in the previous section, we set $Q_{x,y,z} = 1$ or $0$, $R^{-1}M_z = 1$ or $0$.

(i) Method of resolution

Green's functions given in appendix A are convolved with the normalized distributions defined in §2b, yielding the contact compliance matrix $\hat{C}$ when the convolutions are evaluated at point C. The integrals involved are solved analytically in the limit $2a \ll R$ – for the sake of conciseness they are given in appendix Cc.

(ii) Solutions for full-stick, small incremental displacements

We start by giving the compliance corresponding to distributions given by equations (2.4) and (2.5). Only non-zero terms are reported.

(a) Small, incremental normal load

$$4\pi \mu C_{yy} = (\kappa + 1) \log \frac{2R}{a\sqrt{e}} - \frac{1}{2}. \qquad (4.5)$$

(b) Small, incremental in-plane tangential load

$$4\pi \mu C_{xx} = (\kappa + 1) \log \frac{2R}{ae} + 3 - \frac{1}{6} \qquad (4.6a)$$

and

$$4\pi \mu C_{\alpha x} = \frac{R}{2a}\pi(\kappa - 1). \qquad (4.6b)$$

(c) Small, incremental out-of-plane tangential load

$$4\pi \mu C_{zz} = 4 \log \frac{2R}{a} + \frac{1}{2}. \qquad (4.7)$$

(d) Small, incremental dipolar normal load

$$4\pi \mu C_{\alpha\alpha} = \frac{2R^2}{a^2}(\kappa + 1) \qquad (4.8a)$$

and

$$4\pi \mu C_{x\alpha} = -\frac{3}{2}\kappa + \frac{11}{6}, \qquad (4.8b)$$

with $e = \exp(1)$. An application of the normal contact compliance formula (4.5) can be found in [6,23], where contact effects are accounted for in the stiffness matrix of a wire strand.
(iii) Solutions for finite displacements

The formulae for C<sub>xx</sub>, C<sub>yy</sub> and C<sub>zz</sub> corresponding to distributions given by equations (2.6) and (2.7) read

(a) Finite normal load

$$4\pi \mu C_{yy} = (\kappa + 1) \log \frac{2R}{a} - \frac{1}{2}. \qquad (4.9)$$

(b) Finite in-plane tangential load with partial slip, by calling c the half-width of the stick strip

$$4\pi \mu C_{xx} = (\kappa + 1) \log \frac{2R}{a\sqrt{e}} + \frac{(\kappa + 1)c^2}{a^2 - c^2} \log \frac{c}{a} + 3 - \frac{1}{6} \qquad (4.10a)$$

and

$$4\pi \mu C_{\alpha x} = \text{same as equation (4.6)}. \qquad (4.10b)$$

(c) Finite out-of-plane tangential load, with partial slip

$$4\pi \mu C_{zz} = 4 \log \frac{2R\sqrt{e}}{a} + \frac{4c^2}{a^2 - c^2} \log \frac{c}{a} + \frac{1}{2}. \qquad (4.11)$$

(c) Solution: inner field at and far from contact

From these latter two parts, using the state vectors defined in §3c, the inner field can be expressed as

$$\hat{\mathbf{u}}(\mathbf{r}) = \begin{cases} \hat{\mathcal{G}}(\mathbf{r})\hat{\mathcal{Q}} & \text{far from the contact: } \|\mathbf{r} - \mathbf{C}\| \gg a, \\ \hat{\mathcal{C}}(a)\hat{\mathcal{Q}} & \text{at the centre of the contact: } \mathbf{r} = \mathbf{C}, \end{cases} \qquad (4.12)$$

where the far field terms $||\mathbf{r} - \mathbf{C}|| \gg a$ are approximated by merely multiplying Green's functions by the average loads $Q_{x,y,z}$ and $M_z/R$, and where the intermediate region $s \sim a$ is small and then not of interest.

5. Examples

Let us illustrate how the formulae derived in the previous parts can be applied to obtain analytical solutions to problems involving cylinders in contact. In example 5a, we show how the compliance of a simple system of three cylinders can be properly assembled by accounting for non-local terms (i.e. the far field in Green's solutions), by revisiting an example given in [2]. Although there already exists a general recipe to obtain the stress field within a cylinder under multiple contacts ([24]—§37), example 5a indicates how this could be extended to get the displacement field. In example 5b, we derive the free and forced dynamic response of a single cylinder in contact with a stiff soil.

(a) Example 1: Compression of a cylinder in contact with two others

It may be instructive to show how a well-known formula can be recovered using equations (4.1) and (4.9). Consider a cylinder in contact with two others that are located at opposite ends of a diameter. The two external bodies compress the cylinder at a magnitude $Q_y = P$ per unit axial length. This case is treated in [2] §5.6. We call $C_1 = (0, -R)^T$ and $C_2 = (0, R)^T$ the contact points before any pressure is applied, and $a_1(P)$ and $a_2(P)$ the contact half-widths. This situation can be viewed as a superposition of two cases of type (a) treated in §4, with load distributions following equation (2.6): the body forces cancel and only the Hertzian loads remain. We call $\mathbf{u}^1$ and $\mathbf{u}^2$ the solutions to these two cases, given by equations (4.1) and (4.9) under an appropriate change of
(iii) Solutions for finite displacements

The formulae for $C_{xx}$, $C_{yy}$ and $C_{zz}$ corresponding to distributions given by equations (2.6) and (2.7) read

(a) Finite normal load

$$4\pi \mu C_{yy} = (\kappa + 1) \log \frac{2R}{a} - \frac{1}{2}. \qquad (4.9)$$

(b) Finite in-plane tangential load with partial slip, by calling c the half-width of the stick strip

$$4\pi \mu C_{xx} = (\kappa + 1) \log \frac{2R}{a\sqrt{e}} + \frac{(\kappa + 1)c^2}{a^2 - c^2} \log \frac{c}{a} + 3 - \frac{1}{6} \qquad (4.10a)$$

and

$$4\pi \mu C_{\alpha x} = \text{same as equation (4.6)}. \qquad (4.10b)$$

(c) Finite out-of-plane tangential load, with partial slip

$$4\pi \mu C_{zz} = 4 \log \frac{2R\sqrt{e}}{a} + \frac{4c^2}{a^2 - c^2} \log \frac{c}{a} + \frac{1}{2}. \qquad (4.11)$$

(c) Solution: inner field at and far from contact

From these latter two parts, using the state vectors defined in §3c, the inner field can be expressed as

$$\hat{\mathbf{u}}(\mathbf{r}) = \begin{cases} \hat{\mathcal{G}}(\mathbf{r})\hat{\mathcal{Q}} & \text{far from the contact: } \|\mathbf{r} - \mathbf{C}\| \gg a, \\ \hat{\mathcal{C}}(a)\hat{\mathcal{Q}} & \text{at the centre of the contact: } \mathbf{r} = \mathbf{C}, \end{cases} \qquad (4.12)$$

where the far field terms $\|\mathbf{r} - \mathbf{C}\| \gg a$ are approximated by merely multiplying Green's functions by the average loads $Q_{x,y,z}$ and $M_z/R$, and where the intermediate region $s \sim a$ is small and then not of interest.

5. Examples

Let us illustrate how the formulae derived in the previous parts can be applied to obtain analytical solutions to problems involving cylinders in contact. In example 5a, we show how the compliance of a simple system of three cylinders can be properly assembled by accounting for non-local terms (i.e. the far field in Green's solutions), by revisiting an example given in [2]. Although there already exists a general recipe to obtain the stress field within a cylinder under multiple contacts ([24]—§37), example 5a indicates how this could be extended to get the displacement field. In example 5b, we derive the free and forced dynamic response of a single cylinder in contact with a stiff soil.

(a) Example 1: Compression of a cylinder in contact with two others

It may be instructive to show how a well-known formula can be recovered using equations (4.1) and (4.9). Consider a cylinder in contact with two others that are located at opposite ends of a diameter. The two external bodies compress the cylinder at a magnitude $Q_y = P$ per unit axial length. This case is treated in [2] §5.6. We call $C_1 = (0, -R)^T$ and $C_2 = (0, R)^T$ the contact points before any pressure is applied, and $a_1(P)$ and $a_2(P)$ the contact half-widths. This situation can be viewed as a superposition of two cases of type (a) treated in §4, with load distributions following equation (2.6): the body forces cancel and only the Hertzian loads remain. We call $\mathbf{u}^1$ and $\mathbf{u}^2$ the solutions to these two cases, given by equations (4.1) and (4.9) under an appropriate change of
<figure><img src="image_9.png" /><figcaption>Figure 4. A cylinder in Hertzian contact with a vibrating soil, represented at low frequencies as a mass-spring system with 4 d.f.</figcaption></figure>

coordinate system. The displacement $\mathbf{u} = \mathbf{u}^1 + \mathbf{u}^2$ at point $C_1$ is obtained by subtracting equation (4.1) evaluated at $\theta = \pi$ to equation (4.9):

$$ \frac{4\pi\mu}{P} u_y(C_1) = (\kappa + 1) \log \frac{2R}{a_1} - \frac{1}{2} - \left( -(\kappa + 1) \log 2 - 0 + \frac{1}{2}(\kappa - 1) + \frac{1}{2} \right) \\ = (\kappa + 1) \left( \log \frac{4R}{a_1} - \frac{1}{2} \right). \qquad (5.1) $$

A similar expression can be obtained at point $C_2$. The compression of the diameter $\delta = u_y(C_1) - u_y(C_2)$ is then

$$ \delta = P \frac{(\kappa + 1)}{4\pi\mu} \left( \log \frac{4R}{a_1} + \log \frac{4R}{a_2} - 1 \right). \qquad (5.2) $$

With $\kappa = 3 - 4\nu$ and $\mu = E/2(1 + \nu)$, one recovers the formula given in [2] §5.6:

$$ \delta = 2P \frac{(1 - \nu^2)}{\pi E} \left( \log \frac{4R}{a_1} + \log \frac{4R}{a_2} - 1 \right). \qquad (5.3) $$

## (b) Example 2: The heavy cylinder on a vibrating soil

Here, we go back to the minimalist example of a single cylinder and add a dynamic behaviour. The resolution method is to replace contacts with equivalent springs under a low-frequency approximation, as represented in figure 4. We revisit the classical example of the ‘heavy disc’ or ‘heavy cylinder’ [15] by introducing a small dynamic motion in the stiff soil supporting the cylinder. The original problem in Michell’s paper was to find the stress field caused by the weight and a line-wise contact with the soil. Here, the problem is to find the displacement field as a function of the forcing frequency, under the assumption that the contact follows Hertz’s Law.

### (i) Initial state and perturbed state

A cylinder lays over a stiff soil. Under the action of a normal body force of magnitude $P_0$ (e.g. the weight) and the reaction of the soil $q_0$ (Hertzian contact distribution following equation (2.6)), the cylinder deforms into an initial (static) state ($u_0, \sigma_0$) characterized by a contact strip of width $2a$, where $a = \sqrt{4P_0R(1 - \nu^2)/\pi E}$ according to Hertz’s Law.

Then, let us consider that the soil slightly vibrates harmonically owing to a small prescribed displacement $\mathbf{d} e^{i\omega t}$, with $\omega$ the angular frequency and $t$ the time. The related contact load $\mathbf{q} e^{i\omega t}$ follows equations (2.4) and (2.5). This excitation produces a small incremental dynamic displacement $\mathbf{u}_{dyn} e^{i\omega t}$. The total displacement and stress become

$$ \mathbf{u}_{tot}(t) = \mathbf{u}_0 + \mathbf{u}_{dyn} e^{i\omega t} \qquad (5.4a) $$

and

$$ \sigma_{tot}(t) = \sigma_0 + \sigma_{dyn} e^{i\omega t}. \qquad (5.4b) $$
<figure><img src="image_9.png" /><figcaption>Figure 4. A cylinder in Hertzian contact with a vibrating soil, represented at low frequencies as a mass-spring system with 4 d.f.</figcaption></figure>

coordinate system. The displacement $\mathbf{u} = \mathbf{u}^1 + \mathbf{u}^2$ at point $C_1$ is obtained by subtracting equation (4.1) evaluated at $\theta = \pi$ to equation (4.9):

$$ \begin{aligned} \frac{4\pi\mu}{P} u_y(C_1) &= (\kappa + 1) \log \frac{2R}{a_1} - \frac{1}{2} - \left( -(\kappa + 1) \log 2 - 0 + \frac{1}{2}(\kappa - 1) + \frac{1}{2} \right) \\ &= (\kappa + 1) \left( \log \frac{4R}{a_1} - \frac{1}{2} \right). \end{aligned} \quad (5.1) $$

A similar expression can be obtained at point $C_2$. The compression of the diameter $\delta = u_y(C_1) - u_y(C_2)$ is then

$$ \delta = P \frac{(\kappa + 1)}{4\pi\mu} \left( \log \frac{4R}{a_1} + \log \frac{4R}{a_2} - 1 \right). \quad (5.2) $$

With $\kappa = 3 - 4\nu$ and $\mu = E/2(1 + \nu)$, one recovers the formula given in [2] §5.6:

$$ \delta = 2P \frac{(1 - \nu^2)}{\pi E} \left( \log \frac{4R}{a_1} + \log \frac{4R}{a_2} - 1 \right). \quad (5.3) $$

## (b) Example 2: The heavy cylinder on a vibrating soil

Here, we go back to the minimalist example of a single cylinder and add a dynamic behaviour. The resolution method is to replace contacts with equivalent springs under a low-frequency approximation, as represented in figure 4. We revisit the classical example of the 'heavy disc' or 'heavy cylinder' [15] by introducing a small dynamic motion in the stiff soil supporting the cylinder. The original problem in Michell's paper was to find the stress field caused by the weight and a line-wise contact with the soil. Here, the problem is to find the displacement field as a function of the forcing frequency, under the assumption that the contact follows Hertz's Law.

### (i) Initial state and perturbed state

A cylinder lays over a stiff soil. Under the action of a normal body force of magnitude $P_0$ (e.g. the weight) and the reaction of the soil $q_0$ (Hertzian contact distribution following equation (2.6)), the cylinder deforms into an initial (static) state ($u_0, \sigma_0$) characterized by a contact strip of width $2a$, where $a = \sqrt{4P_0R(1 - \nu^2)/\pi E}$ according to Hertz's Law.

Then, let us consider that the soil slightly vibrates harmonically owing to a small prescribed displacement $\mathbf{d} e^{i\omega t}$, with $\omega$ the angular frequency and $t$ the time. The related contact load $\mathbf{q} e^{i\omega t}$ follows equations (2.4) and (2.5). This excitation produces a small incremental dynamic displacement $\mathbf{u}_{dyn} e^{i\omega t}$. The total displacement and stress become

$$ \mathbf{u}_{tot}(t) = \mathbf{u}_0 + \mathbf{u}_{dyn} e^{i\omega t} \quad (5.4a) $$

and

$$ \sigma_{tot}(t) = \sigma_0 + \sigma_{dyn} e^{i\omega t}. \quad (5.4b) $$
For simplicity, we use a coordinate system centred on the centre of mass of the cylinder in the initial state $\mathbf{u}_0$, such as in the previous sections (figure 1). From here on, only the dynamic part $\mathbf{u}_{dyn}$ is of interest and we apply the notations defined earlier to it: U is the skeleton part of $\mathbf{u}_{dyn}$, u is the inner field part, etc.

## (ii) Solutions: forced motion and contact resonances

Equations (3.4), (3.7) and (3.11) can be combined to obtain an equation on the skeleton field only

$$(-\omega^2 \hat{\mathcal{M}} + \hat{\mathcal{K}}) \hat{U} = \hat{\mathcal{K}} \hat{\mathcal{R}}^{-1} \hat{\mathbf{d}}|_{\mathcal{C}}, \qquad (5.5)$$

with

$$\hat{\mathcal{K}} = \hat{\mathcal{R}}^T \hat{\mathcal{C}}^{-1} \hat{\mathcal{R}}, \qquad (5.6)$$

the (real, symmetric, positive definite) stiffness matrix. The solution of equation (5.5) (forced motion) can be straightforwardly obtained by inverting the left-hand side matrix. It is however more insightful to express the solution as a combination of the eigenmodes of the system: we call $(\omega_n^2, \hat{U}_n)$ the eigenvalues and eigenvectors of the Hermitian generalized eigenvalue problem $\hat{\mathcal{K}}\hat{U}_n = \omega_n^2 \hat{\mathcal{M}}\hat{U}_n$ and require the following normalization: $\hat{U}_n^T \hat{\mathcal{M}} \hat{U}_n = m$. Then,

$$\hat{U} = \sum_{n=1}^{4} d_n \frac{1}{1 - \omega^2 / \omega_n^2} \hat{U}_n, \qquad (5.7)$$

in which $d_n = \hat{\mathbf{d}}_n^T \hat{\mathbf{d}}|_{\mathcal{C}}$ is the participation factor of mode $n$, and with $\hat{\mathbf{d}}_n = m^{-1}(\hat{\mathbf{R}}^T)^{-1} \hat{\mathcal{M}} \hat{U}_n$. The contact forces can be obtained using equation (3.11)

$$\hat{\mathbf{Q}} = \sum_{n=1}^{4} d_n \frac{\omega^2 / \omega_n^2}{1 - \omega^2 / \omega_n^2} \hat{\mathbf{Q}}_n, \qquad (5.8)$$

where we have defined $\hat{\mathbf{Q}}_n = -m\omega_n^2 \hat{\mathbf{d}}_n$ the eigen contact forces. Knowing the contact efforts $\hat{\mathbf{Q}}$ at a given frequency $\omega$, the inner field $\mathbf{u}(\mathbf{r}, \omega)$ can be obtained at any point using equation (4.12). Finally, the total dynamic field $\mathbf{u}_{dyn}(\mathbf{r}, \omega)$ can be obtained using equation (3.1a).

The eigenvalue problem is simple enough for analytical solutions to be derived; indeed only the $U_x$ and $\Theta$ components are coupled. The four eigenmodes are

$$\omega_1 = \sqrt{\frac{1}{C_{\alpha\alpha}m_{\alpha}}}, \quad \hat{U}_1 = \sqrt{\frac{2}{3}(-1,0,0,1)^T}, \quad \hat{\mathbf{d}}_1 = \sqrt{\frac{2}{3}(-1,0,0,\frac{3}{2})^T}, \quad \hat{\mathbf{Q}}_1 = -\frac{2}{3}C_{\alpha\alpha}^{-1}\hat{\mathbf{d}}_1, \qquad (5.9a)$$

$$\omega_2 = \sqrt{\frac{1}{C_{zz}m}}, \quad \hat{U}_2 = (0,0,1,0)^T, \quad \hat{\mathbf{d}}_2 = \hat{U}_2, \quad \hat{\mathbf{Q}}_2 = -C_{zz}^{-1} \hat{\mathbf{d}}_2, \qquad (5.9b)$$

$$\omega_3 = \sqrt{\frac{1}{C_{yy}m}}, \quad \hat{U}_3 = (0,1,0,0)^T, \quad \hat{\mathbf{d}}_3 = \hat{U}_3, \quad \hat{\mathbf{Q}}_3 = -C_{yy}^{-1} \hat{\mathbf{d}}_3 \qquad (5.9c)$$

$$\text{and} \quad \omega_4 = \sqrt{\frac{1}{C_{xx}m_x}}, \quad \hat{U}_4 = \frac{1}{\sqrt{3}}(1,0,0,2)^T, \quad \hat{\mathbf{d}}_4 = \frac{1}{\sqrt{3}}(1,0,0,0)^T, \quad \hat{\mathbf{Q}}_4 = -3C_{xx}^{-1}\hat{\mathbf{d}}_4, \qquad (5.9d)$$

$$\text{with } m_x = (m^{-1} + (J/R^2)^{-1})^{-1} = 1/3 \text{ m and } m_\alpha = 3/2 \text{ m.}$$

**Remark.** Equations (5.9a) and (5.9d) were obtained by noticing that $C_{xx} \ll C_{\alpha\alpha}$, expanding in powers of $C_{xx}/C_{\alpha\alpha}$ and keeping only the leading term.

## (iii) Numerical evaluation, discussion

Let us now apply equations (5.7) and (5.9) with parameters representing a steel cylinder: $E = 217$ GPa, $\nu = 0.28$, $\rho = 7.8$ g cm$^{-3}$, $R = 3.6$ mm. The contact half-width is set to $a/R = 1\%$: this order of magnitude can be reached in a wire strand between a peripheral wire and the core when
For simplicity, we use a coordinate system centred on the centre of mass of the cylinder in the initial state $\mathbf{u}_0$, such as in the previous sections (figure 1). From here on, only the dynamic part $\mathbf{u}_{dyn}$ is of interest and we apply the notations defined earlier to it: U is the skeleton part of $\mathbf{u}_{dyn}$, u is the inner field part, etc.

## (ii) Solutions: forced motion and contact resonances

Equations (3.4), (3.7) and (3.11) can be combined to obtain an equation on the skeleton field only

$$(-\omega^2 \hat{\mathcal{M}} + \hat{\mathcal{K}}) \hat{U} = \hat{\mathcal{K}} \hat{\mathcal{R}}^{-1} \hat{\mathbf{d}}|_{\mathcal{C}}, \quad (5.5)$$

with

$$\hat{\mathcal{K}} = \hat{\mathcal{R}}^T \hat{\mathcal{C}}^{-1} \hat{\mathcal{R}}, \quad (5.6)$$

the (real, symmetric, positive definite) stiffness matrix. The solution of equation (5.5) (forced motion) can be straightforwardly obtained by inverting the left-hand side matrix. It is however more insightful to express the solution as a combination of the eigenmodes of the system: we call $(\omega_n^2, \hat{U}_n)$ the eigenvalues and eigenvectors of the Hermitian generalized eigenvalue problem $\hat{\mathcal{K}}\hat{U}_n = \omega_n^2 \hat{\mathcal{M}}\hat{U}_n$ and require the following normalization: $\hat{U}_n^T \hat{\mathcal{M}} \hat{U}_n = m$. Then,

$$\hat{U} = \sum_{n=1}^{4} d_n \frac{1}{1 - \omega^2 / \omega_n^2} \hat{U}_n, \quad (5.7)$$

in which $d_n = \hat{\mathbf{d}}_n^T \hat{\mathbf{d}}|_{\mathcal{C}}$ is the participation factor of mode $n$, and with $\hat{\mathbf{d}}_n = m^{-1}(\hat{\mathcal{R}}^T)^{-1} \hat{\mathcal{M}} \hat{U}_n$. The contact forces can be obtained using equation (3.11)

$$\hat{\mathbf{Q}} = \sum_{n=1}^{4} d_n \frac{\omega^2 / \omega_n^2}{1 - \omega^2 / \omega_n^2} \hat{\mathbf{Q}}_n, \quad (5.8)$$

where we have defined $\hat{\mathbf{Q}}_n = -m\omega_n^2 \hat{\mathbf{d}}_n$ the eigen contact forces. Knowing the contact efforts $\hat{\mathbf{Q}}$ at a given frequency $\omega$, the inner field $\mathbf{u}(\mathbf{r}, \omega)$ can be obtained at any point using equation (4.12). Finally, the total dynamic field $\mathbf{u}_{dyn}(\mathbf{r}, \omega)$ can be obtained using equation (3.1a).

The eigenvalue problem is simple enough for analytical solutions to be derived; indeed only the $U_x$ and $\Theta$ components are coupled. The four eigenmodes are

$$\omega_1 = \sqrt{\frac{1}{C_{\alpha\alpha} m_{\alpha}}}, \quad \hat{U}_1 = \sqrt{\frac{2}{3}(-1, 0, 0, 1)^T}, \quad \hat{\mathbf{d}}_1 = \sqrt{\frac{2}{3}(-1, 0, 0, \frac{3}{2})^T}, \quad \hat{\mathbf{Q}}_1 = -\frac{2}{3} C_{\alpha\alpha}^{-1} \hat{\mathbf{d}}_1, \quad (5.9a)$$

$$\omega_2 = \sqrt{\frac{1}{C_{zz} m}}, \quad \hat{U}_2 = (0, 0, 1, 0)^T, \quad \hat{\mathbf{d}}_2 = \hat{U}_2, \quad \hat{\mathbf{Q}}_2 = -C_{zz}^{-1} \hat{\mathbf{d}}_2, \quad (5.9b)$$

$$\omega_3 = \sqrt{\frac{1}{C_{yy} m}}, \quad \hat{U}_3 = (0, 1, 0, 0)^T, \quad \hat{\mathbf{d}}_3 = \hat{U}_3, \quad \hat{\mathbf{Q}}_3 = -C_{yy}^{-1} \hat{\mathbf{d}}_3 \quad (5.9c)$$

$$\text{and} \quad \omega_4 = \sqrt{\frac{1}{C_{xx} m_x}}, \quad \hat{U}_4 = \frac{1}{\sqrt{3}}(1, 0, 0, 2)^T, \quad \hat{\mathbf{d}}_4 = \frac{1}{\sqrt{3}}(1, 0, 0, 0)^T, \quad \hat{\mathbf{Q}}_4 = -3C_{xx}^{-1} \hat{\mathbf{d}}_4, \quad (5.9d)$$

$$\text{with } m_x = (m^{-1} + (J/R^2)^{-1})^{-1} = 1/3 \text{ m and } m_\alpha = 3/2 \text{ m.}$$

**Remark.** Equations (5.9a) and (5.9d) were obtained by noticing that $C_{xx} \ll C_{\alpha\alpha}$, expanding in powers of $C_{xx}/C_{\alpha\alpha}$ and keeping only the leading term.

## (iii) Numerical evaluation, discussion

Let us now apply equations (5.7) and (5.9) with parameters representing a steel cylinder: $E = 217 \text{ GPa}$, $\nu = 0.28$, $\rho = 7.8 \text{ g cm}^{-3}$, $R = 3.6 \text{ mm}$. The contact half-width is set to $a/R = 1\%$: this order of magnitude can be reached in a wire strand between a peripheral wire and the core when
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Forced response of a steel heavy cylinder in Hertzian contact with a vibrating soil, for different vibration polarizations, obtained using the analytical solutions in the quasi-static approximation (solid lines) or using finite elements (FE) (dashed lines, see appendix B for details). The first four resonances are the contact resonances of the system, then elastic resonances start to contribute. Contact half-width: $a/R = 1\%$.</figcaption></figure>

applying a tension close to the maximal service load—as a comparison, the action of the weight only yields a much weaker value ($a/R \sim 6.8 \times 10^{-3} \%$). The frequencies will be normalized with a characteristic angular frequency of the cylinder $\omega_{cyl} = c_s/R$, with $c_s = \sqrt{E/(2(1+\nu))\rho}$ the shear wave velocity.

Figure 5 represents the forced response (equation (5.7)) at the top of the cylinder ($u_{dyn}|_{r=R, \theta=\pi}$) as a function of the normalized driving frequency $\omega/\omega_{cyl}$, for different polarizations of the soil motion. A reference solution obtained with FE—see appendix B for details—is superimposed to show the transition from which the low-frequency hypothesis stops being valid: here we observe excellent agreement up to roughly $\omega/\omega_{cyl} \approx 1.5$, after which elastic resonances start to dominate the solution. It is noticeable that the lowest contact resonance frequency is orders of magnitude below the other three. This is because the response to the dipolar normal load is qualitatively different: the fundamental solution (displacement) behaves as $1/s$ near the load instead of logs, leading to a contact compliance in $(a/R)^{-2}$ instead of $-\log a/R + \text{const.}$, hence resulting into a much higher compliance for $a/R \ll 1$.

Figure 6 illustrates further this dependence on the contact width. The eigen angular frequencies $\omega_n$ (equations (5.9)) are represented as a function of $a/R$, emphasizing three families of modes: one pendulum-like contact mode with $\omega_1 \sim a/R$, then three other contact modes with $\omega_{2,3,4} \sim 1/\sqrt{-\log a/R + \text{const.}}$, responding to piston or tangential forces, and finally elastic resonances in infinite number from (depending on $a/R$) half an order of magnitude above the highest contact mode, whose frequencies are almost independent on the contact width in the range compatible with the Hertzian approximation. Figure 6 also illustrates the decomposition of the displacement field into a *skeleton* field that allows for an average description of the motion, and an *inner* field that captures the deformation due to the contact.

## 6. Conclusion

This manuscript gives a derivation of the compliance of an elastic cylinder to surface forces representing a Hertzian contact, that are balanced by bulk forces representing the inertia of rigid body motions. It contains an extensive set of closed form formulae, starting from Green's functions of the problems for the displacement and stress fields, to the value of the response at the contact point. An original aspect that relates to a long standing discussion in the literature is how displacements are made unique: here integral gauges have been introduced whose justification
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Forced response of a steel heavy cylinder in Hertzian contact with a vibrating soil, for different vibration polarizations, obtained using the analytical solutions in the quasi-static approximation (solid lines) or using finite elements (FE) (dashed lines, see appendix B for details). The first four resonances are the contact resonances of the system, then elastic resonances start to contribute. Contact half-width: $a/R = 1\%$.</figcaption></figure>

applying a tension close to the maximal service load—as a comparison, the action of the weight only yields a much weaker value ($a/R \sim 6.8 \times 10^{-3} \%$). The frequencies will be normalized with a characteristic angular frequency of the cylinder $\omega_{cyl} = c_s/R$, with $c_s = \sqrt{E/(2(1+\nu))\rho}$ the shear wave velocity.

Figure 5 represents the forced response (equation (5.7)) at the top of the cylinder ($u_{dyn} |_{r=R, \theta=\pi}$) as a function of the normalized driving frequency $\omega/\omega_{cyl}$, for different polarizations of the soil motion. A reference solution obtained with FE—see appendix B for details—is superimposed to show the transition from which the low-frequency hypothesis stops being valid: here we observe excellent agreement up to roughly $\omega/\omega_{cyl} \approx 1.5$, after which elastic resonances start to dominate the solution. It is noticeable that the lowest contact resonance frequency is orders of magnitude below the other three. This is because the response to the dipolar normal load is qualitatively different: the fundamental solution (displacement) behaves as $1/s$ near the load instead of logs, leading to a contact compliance in $(a/R)^{-2}$ instead of $-\log a/R + \text{const.}$, hence resulting into a much higher compliance for $a/R \ll 1$.

Figure 6 illustrates further this dependence on the contact width. The eigen angular frequencies $\omega_n$ (equations (5.9)) are represented as a function of $a/R$, emphasizing three families of modes: one pendulum-like contact mode with $\omega_1 \sim a/R$, then three other contact modes with $\omega_{2,3,4} \sim 1/\sqrt{-\log a/R + \text{const.}}$, responding to piston or tangential forces, and finally elastic resonances in infinite number from (depending on $a/R$) half an order of magnitude above the highest contact mode, whose frequencies are almost independent on the contact width in the range compatible with the Hertzian approximation. Figure 6 also illustrates the decomposition of the displacement field into a *skeleton* field that allows for an average description of the motion, and an *inner* field that captures the deformation due to the contact.

## 6. Conclusion

This manuscript gives a derivation of the compliance of an elastic cylinder to surface forces representing a Hertzian contact, that are balanced by bulk forces representing the inertia of rigid body motions. It contains an extensive set of closed form formulae, starting from Green's functions of the problems for the displacement and stress fields, to the value of the response at the contact point. An original aspect that relates to a long standing discussion in the literature is how displacements are made unique: here integral gauges have been introduced whose justification
<figure><img src="image_11.png" /><figcaption>**Figure 6.** First eigenmodes of a steel heavy cylinder in Hertzian contact with the soil, obtained using the analytical solutions in the quasi-static approximation or using FE, see appendix B for details. (a) Eigenfrequency as a function of the contact half-width *a*. The dashed lines stand for FE solutions. (b) Mode shapes; analytical solutions are represented by filtering out or including the inner field. Pink-green colourmap depicts out-of-plane modes.</figcaption></figure>

arises when addressing dynamic problems. From these formulae, the forced and free response of a cylinder in contact with a vibrating soil have been derived.

The natural way to use these formulae is to address systems of several—or many—cylinders in contact, either in static equilibrium or in linear or nonlinear dynamic motion around a given static state, and construct the compliance matrix by linear superposition. For statics there already exists a general solution to obtain the stress within a cylinder submitted to multiple contacts, here the benefit would be to also obtain the displacement. For dynamics, one could think of modelling vibrations in a lattice of cylinders, obtaining the normal modes of a stack of cylinders, or, more specifically, modelling the behaviour of wire strands at low frequencies. For this latter perspective, it seems a promising way to couple the present results to beam theory. Another further step could be to model frictional hysteretic damping by employing the contact compliance with partial slip, a case which has been only barely sketched here.

A limitation is the geometry of the cross section, which has to be purely circular. In the context of wire ropes, there exists a wide variety of architectures that also resort to exotic cross sections (e.g. Z-shaped). Devising a numerical strategy to adapt the current formulae to these shapes would also be beneficial.

**Data accessibility.** This article has no additional data.

**Declaration of AI use.** We have not used AI-assisted technologies in creating this article.

**Authors' contributions.** P.M.: conceptualization, formal analysis, investigation, methodology, validation, writing—original draft, writing—review and editing; F.T.: formal analysis, funding acquisition, writing—original draft, writing—review and editing; L.L.: formal analysis, writing—original draft, writing—review and editing.
<figure><img src="image_11.png" /><figcaption>**Figure 6.** First eigenmodes of a steel heavy cylinder in Hertzian contact with the soil, obtained using the analytical solutions in the quasi-static approximation or using FE, see appendix B for details. (a) Eigenfrequency as a function of the contact half-width *a*. The dashed lines stand for FE solutions. (b) Mode shapes; analytical solutions are represented by filtering out or including the inner field. Pink-green colourmap depicts out-of-plane modes.</figcaption></figure>

arises when addressing dynamic problems. From these formulae, the forced and free response of a cylinder in contact with a vibrating soil have been derived.

The natural way to use these formulae is to address systems of several—or many—cylinders in contact, either in static equilibrium or in linear or nonlinear dynamic motion around a given static state, and construct the compliance matrix by linear superposition. For statics there already exists a general solution to obtain the stress within a cylinder submitted to multiple contacts, here the benefit would be to also obtain the displacement. For dynamics, one could think of modelling vibrations in a lattice of cylinders, obtaining the normal modes of a stack of cylinders, or, more specifically, modelling the behaviour of wire strands at low frequencies. For this latter perspective, it seems a promising way to couple the present results to beam theory. Another further step could be to model frictional hysteretic damping by employing the contact compliance with partial slip, a case which has been only barely sketched here.

A limitation is the geometry of the cross section, which has to be purely circular. In the context of wire ropes, there exists a wide variety of architectures that also resort to exotic cross sections (e.g. Z-shaped). Devising a numerical strategy to adapt the current formulae to these shapes would also be beneficial.

**Data accessibility.** This article has no additional data.

**Declaration of AI use.** We have not used AI-assisted technologies in creating this article.

**Authors' contributions.** P.M.: conceptualization, formal analysis, investigation, methodology, validation, writing—original draft, writing—review and editing; F.T.: formal analysis, funding acquisition, writing—original draft, writing—review and editing; L.L.: formal analysis, writing—original draft, writing—review and editing.
All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

**Conflict of interest declaration.** We declare we have no competing interests.

**Funding.** P.M. acknowledges funding of the PULSAR-2022 program, Région Pays de la Loire, grant USonNoLiMatGC.

**Acknowledgements.** We are thankful to Marc Bonnet (CNRS, POems team) for stimulating discussions and for indicating the value of Integrals (C 4), (C 5), to Patrice Cartraud (École Centrale de Nantes, GeM lab.) for bibliographic references, and to Konstantinos Kondylidis (Univ. Eiffel, GeoEND lab.) for proof-reading. We also acknowledge the anonymous reviewers for their help in improving the clarity of the manuscript.

# Appendix A. Full solutions of Green's inner static problems

We here recall how the four Green's problems defined in §4 can be solved. We give an expression of the displacement field that is valid at all points inside the cylinder, from which the simplified formulae (4.1)–(4.4) were derived. We also report the partial formulae from which the total stress field can be obtained.

The solutions are represented in figure 7. It can be observed that while the displacement turns out to be zero on the cylinder axis for case (a), and while in case (d) the far-point concept would also reasonably apply for that point, cases (b) and (c) have no trivial fixed point. A validation of the formulae with FE is shown in figure 8, see details in appendix B.

## (a) Method of resolution

Green's functions can be constructed by superposing:

- the response of a half-space to a load applied on its boundary (the Flamant solution),

- the contribution of balancing body forces (gravitational-like loading, and rotational acceleration),

- solutions of the bi-harmonic equation chosen to satisfy the boundary condition (the Michell solutions),

- rigid-body displacements chosen to satisfy Gauges (3.2).

Corresponding formulae can either be readily found in [14], or be deduced with simple algebra from solutions in stress reported in [14]. The integrals required for the gauges are given in appendix C. Case (c) is simpler, it is amenable to a more straightforward resolution using separation of variables. Case (d) is qualitatively different because the dipolar load leads to a decrease of the fundamental solution as $1/||\mathbf{r} - \mathbf{C}||$ rather than $\log ||\mathbf{r} - \mathbf{C}||$: for our purpose, it would be enough to only account for the local curvature and prescribe zero displacement of a 'far point' (half-space approximation), as done in three dimensions for deriving the compliance of a point-wise contact. However, for the sake of uniformity in the derivations, we treat case (d) in the same manner as the others and give the solution to the entire problem.

## (b) Polar coordinates (s, β)

Parts of the formulae can be most conveniently expressed by using polar coordinates centred on point C, whose definition is recalled here (figure 1): $s = ||\mathbf{r} - \mathbf{C}||$ and $\tan \beta = -r \sin \theta / (R - r \cos \theta)$, or equivalently $x = -s \sin \beta$ and $y + R = s \cos \beta$, with $-\pi/2 \le \beta \le \pi/2$.

## (c) Normal load

We here treat the case depicted in figure 3a, which can represent a cylinder that deforms under its own weight. This case is a classic since [15] (the 'heavy disc' or 'heavy cylinder') and is often given as an exercise in textbooks, where the problem is to find the stress field (e.g. [14] §12). However, in spite of that, it seems hardly possible to find an explicit report of the displacement field caused
All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

**Conflict of interest declaration.** We declare we have no competing interests.

**Funding.** P.M. acknowledges funding of the PULSAR-2022 program, Région Pays de la Loire, grant USonNoLiMatGC.

**Acknowledgements.** We are thankful to Marc Bonnet (CNRS, POems team) for stimulating discussions and for indicating the value of Integrals (C 4), (C 5), to Patrice Cartraud (École Centrale de Nantes, GeM lab.) for bibliographic references, and to Konstantinos Kondylidis (Univ. Eiffel, GeoEND lab.) for proof-reading. We also acknowledge the anonymous reviewers for their help in improving the clarity of the manuscript.

# Appendix A. Full solutions of Green's inner static problems

We here recall how the four Green's problems defined in §4 can be solved. We give an expression of the displacement field that is valid at all points inside the cylinder, from which the simplified formulae (4.1)–(4.4) were derived. We also report the partial formulae from which the total stress field can be obtained.

The solutions are represented in figure 7. It can be observed that while the displacement turns out to be zero on the cylinder axis for case (a), and while in case (d) the far-point concept would also reasonably apply for that point, cases (b) and (c) have no trivial fixed point. A validation of the formulae with FE is shown in figure 8, see details in appendix B.

## (a) Method of resolution

Green's functions can be constructed by superposing:

- the response of a half-space to a load applied on its boundary (the Flamant solution),

- the contribution of balancing body forces (gravitational-like loading, and rotational acceleration),

- solutions of the bi-harmonic equation chosen to satisfy the boundary condition (the Michell solutions),

- rigid-body displacements chosen to satisfy Gauges (3.2).

Corresponding formulae can either be readily found in [14], or be deduced with simple algebra from solutions in stress reported in [14]. The integrals required for the gauges are given in appendix C. Case (c) is simpler, it is amenable to a more straightforward resolution using separation of variables. Case (d) is qualitatively different because the dipolar load leads to a decrease of the fundamental solution as $1/||r - C||$ rather than $\log ||r - C||$: for our purpose, it would be enough to only account for the local curvature and prescribe zero displacement of a 'far point' (half-space approximation), as done in three dimensions for deriving the compliance of a point-wise contact. However, for the sake of uniformity in the derivations, we treat case (d) in the same manner as the others and give the solution to the entire problem.

## (b) Polar coordinates (s, β)

Parts of the formulae can be most conveniently expressed by using polar coordinates centred on point C, whose definition is recalled here (figure 1): $s = ||r - C||$ and $\tan \beta = -r \sin \theta / (R - r \cos \theta)$, or equivalently $x = -s \sin \beta$ and $y + R = s \cos \beta$, with $-\pi/2 \le \beta \le \pi/2$.

## (c) Normal load

We here treat the case depicted in figure 3a, which can represent a cylinder that deforms under its own weight. This case is a classic since [15] (the 'heavy disc' or 'heavy cylinder') and is often given as an exercise in textbooks, where the problem is to find the stress field (e.g. [14] §12). However, in spite of that, it seems hardly possible to find an explicit report of the displacement field caused
<figure><img src="image_12.png" /><figcaption>**Figure 7.** (a–d) Analytical Green’s functions: displacement field for a Poisson’s ratio $\nu = 0.28$. Log scale colours are proportional to the magnitude of the displacement.</figcaption></figure>

by this system of forces. The solution is partially given in [7]—where only the displacement for the Hertzian problem at the contact point is derived. We did find reports of the solution in [23] and [5], however, with typographic errors—the reference quoted in [23] that supposedly contains the proof is a real challenge to find.

The solution can be constructed by adding the contributions of:

— the Flamant solution

— Airy stress function: $\phi^F = -\pi^{-1}s\beta \sin \beta$,

— Stress tensor: $2\pi(\sigma_{ss}^F, \sigma_{\beta\beta}^F, \sigma_{s\beta}^F) = -(4 \cos \beta / s, 0, 0)$,

— the irrotational force field $\mathbf{b}^{tr} = -\nabla V$

— Potential: $V = y/\pi R^2$,

— Stress tensor: $2\pi(\sigma_{rr}^{tr}, \sigma_{\theta\theta}^{tr}, \sigma_{r\theta}^{tr}) = (2\pi V, 2\pi V, 0)$,

— the following two corrective fields:

1. Airy stress function: $\phi^{C1} = r^2/(4\pi R)$,

1. Stress tensor: $2\pi(\sigma_{rr}^{C1}, \sigma_{\theta\theta}^{C1}, \sigma_{r\theta}^{C1}) = (R^{-1}, R^{-1}, 0)$,

2. Airy stress function: $\phi^{C2} = r^3 \cos \theta / (4\pi R^2)$,

2. Stress tensor: $2\pi(\sigma_{rr}^{C2}, \sigma_{\theta\theta}^{C2}, \sigma_{r\theta}^{C2}) = rR^{-2}(\cos \theta, 3 \cos \theta, \sin \theta)$,

giving the following result for the displacement field:

$$4\pi \mu \mathcal{G}_{xy} = -\sin 2\beta + (\kappa - 1)\beta + \frac{r}{2R}(\kappa - 1)\sin \theta + \frac{r^2}{2R^2}\sin 2\theta \quad (A 1a)$$

and

$$4\pi \mu \mathcal{G}_{yy} = -(\kappa + 1) \log \frac{s}{R} - 2 \sin^2 \beta - \frac{r}{2R}(\kappa - 1) \cos \theta + \frac{r^2}{2R^2}(1 + 2 \sin^2 \theta). \quad (A 1b)$$

This solution is represented in figures 7a and 8a. As discussed in the introduction, it can be observed that the displacement is indeed not zero far from the source, which calls for a careful treatment of how references are measured (through gauges in this manuscript). This is in contrast to the case of a point contact, or to the dipolar case (d).

**Remark.** The displacement at point O is zero: this point remains the centre of mass of the cylinder. This is not the case for tangential loads.

### (d) In-plane tangential load

We here treat the case represented in figure 3b. Unlike case (a), this case seems not to appear in the literature, even though the closely related case of a cylinder accelerated by two diametrically opposed sources can be found in [14], §12.
<figure><img src="image_12.png" /><figcaption>**Figure 7.** (a–d) Analytical Green’s functions: displacement field for a Poisson’s ratio $v = 0.28$. Log scale colours are proportional to the magnitude of the displacement.</figcaption></figure>

by this system of forces. The solution is partially given in [7]—where only the displacement for the Hertzian problem at the contact point is derived. We did find reports of the solution in [23] and [5], however, with typographic errors—the reference quoted in [23] that supposedly contains the proof is a real challenge to find.

The solution can be constructed by adding the contributions of:

— the Flamant solution

— Airy stress function: $\phi^F = -\pi^{-1}s\beta \sin \beta$,

— Stress tensor: $2\pi(\sigma_{ss}^F, \sigma_{\beta\beta}^F, \sigma_{s\beta}^F) = -(4 \cos \beta / s, 0, 0)$,

— the irrotational force field $\mathbf{b}^{tr} = -\nabla V$

— Potential: $V = y/\pi R^2$,

— Stress tensor: $2\pi(\sigma_{rr}^{tr}, \sigma_{\theta\theta}^{tr}, \sigma_{r\theta}^{tr}) = (2\pi V, 2\pi V, 0)$,

— the following two corrective fields:

1. Airy stress function: $\phi^{C1} = r^2/(4\pi R)$,

1. Stress tensor: $2\pi(\sigma_{rr}^{C1}, \sigma_{\theta\theta}^{C1}, \sigma_{r\theta}^{C1}) = (R^{-1}, R^{-1}, 0)$,

2. Airy stress function: $\phi^{C2} = r^3 \cos \theta / (4\pi R^2)$,

2. Stress tensor: $2\pi(\sigma_{rr}^{C2}, \sigma_{\theta\theta}^{C2}, \sigma_{r\theta}^{C2}) = rR^{-2}(\cos \theta, 3 \cos \theta, \sin \theta)$,

giving the following result for the displacement field:

$$4\pi \mu \mathcal{G}_{xy} = -\sin 2\beta + (\kappa - 1)\beta + \frac{r}{2R}(\kappa - 1)\sin \theta + \frac{r^2}{2R^2}\sin 2\theta \quad (A 1a)$$

and

$$4\pi \mu \mathcal{G}_{yy} = -(\kappa + 1) \log \frac{s}{R} - 2 \sin^2 \beta - \frac{r}{2R}(\kappa - 1) \cos \theta + \frac{r^2}{2R^2}(1 + 2 \sin^2 \theta). \quad (A 1b)$$

This solution is represented in figures 7a and 8a. As discussed in the introduction, it can be observed that the displacement is indeed not zero far from the source, which calls for a careful treatment of how references are measured (through gauges in this manuscript). This is in contrast to the case of a point contact, or to the dipolar case (d).

**Remark.** The displacement at point O is zero: this point remains the centre of mass of the cylinder. This is not the case for tangential loads.

### (d) In-plane tangential load

We here treat the case represented in figure 3b. Unlike case (a), this case seems not to appear in the literature, even though the closely related case of a cylinder accelerated by two diametrically opposed sources can be found in [14], §12.
The solution can be constructed by adding the contributions of

— the Flamant solution

— Airy stress function: $\phi^F = -\pi^{-1}s\beta \cos \beta$,

— Stress tensor: $2\pi(\sigma_{ss}^F, \sigma_{\beta\beta}^F, \sigma_{s\beta}^F) = (4 \sin \beta / s, 0, 0)$,

— the irrotational force field $\mathbf{b}^{tr} = -\nabla V$

— Potential: $V = x/\pi R^2$,

— Stress tensor: $2\pi(\sigma_{rr}^{tr}, \sigma_{\theta\theta}^{tr}, \sigma_{r\theta}^{tr}) = (2\pi V, 2\pi V, 0)$,

— the solenoidal force field $\mathbf{b}^{rot} = \nabla \times \mathbf{A}$

— Vector potential: $\mathbf{A} = r^2/(\pi R^3)\mathbf{e}_z$,

— Stress tensor: $2\pi(\sigma_{rr}^{rot}, \sigma_{\theta\theta}^{rot}, \sigma_{r\theta}^{rot}) = (0, 0, \pi A)$,

— a corrective field

— Airy stress function: $\phi^C = -r^3 \sin \theta / (4\pi R^2)$,

— Stress tensor: $2\pi(\sigma_{rr}^C, \sigma_{\theta\theta}^C, \sigma_{r\theta}^C) = rR^{-2}(-\sin \theta, -3 \sin \theta, \cos \theta)$,

giving the following result for the displacement field:

$$4\pi \mu \mathcal{G}_{xx} = -(\kappa + 1) \log \frac{s}{R} + 2 \sin^2 \beta + \frac{r^2}{2R^2} (1 + 2 \cos^2 \theta) + \left( \frac{r^3}{R^3} - \left( \kappa + \frac{5}{3} \right) \frac{r}{R} \right) \cos \theta - 1 \quad (A 2a)$$

and

$$4\pi \mu \mathcal{G}_{yx} = -\sin 2\beta - (\kappa - 1)\beta + \frac{r^2}{2R^2} \sin 2\theta + \left( \frac{r^3}{R^3} - \left( \kappa + \frac{5}{3} \right) \frac{r}{R} \right) \sin \theta. \quad (A 2b)$$

This solution is represented in figures 7b and 8b.

**Remark.** In contrast to the case of a normal load, the displacement at point O is here non-zero.

### (e) Out-of-plane tangential load

We here treat the case represented in figure 3c. This case is scalar, and is simple enough to be solved in a tractable manner by separation of variables. The solution is

$$4\pi \mu \mathcal{G}_{zz} = -4 \log \frac{s}{R} + \frac{r^2}{R^2} - \frac{1}{2}. \quad (A 3)$$

This solution is represented in figures 7c and 8c. The shear stresses $\sigma_{rz} = \frac{1}{2}\mu\partial_r\mathcal{G}_{zz}$ and $\sigma_{\theta z} = \frac{1}{2}\mu\partial_\theta\mathcal{G}_{zz}/r$ can be derived using $s\partial_r s = r - R \cos \theta$ and $s\partial_\theta s = rR \sin \theta$—obtained from $s^2 = R^2 - 2rR \cos \theta + r^2$. One obtains: $2\pi\sigma_{rz} = -(r - R \cos \theta)/s^2 + r/(2R^2)$ and $2\pi\sigma_{\theta z} = -R \sin \theta/s^2$.

**Remark.** Here again, in contrast to the case of a normal load, the displacement at point O is non-zero.

### (f) Dipolar normal load

We here treat the case represented in figure 3d. This case is qualitatively different to the others as the fundamental solution has a stronger divergence at the source, and decreases to zero in the far field. In principle, approximating the contact compliance at the leading order could then be done by using the half-space approximation and relying on the 'far point' concept, as classically done for a point-wise (three-dimensional) contact compliance. Nevertheless, for the sake of a uniform presentation, we give here the complete solution of case (d).

The solution can be constructed by adding the contributions of:

— the Flamant solution
The solution can be constructed by adding the contributions of

— the Flamant solution

— Airy stress function: $\phi^F = -\pi^{-1}s\beta \cos \beta$,

— Stress tensor: $2\pi(\sigma_{ss}^F, \sigma_{\beta\beta}^F, \sigma_{s\beta}^F) = (4 \sin \beta / s, 0, 0)$,

— the irrotational force field $\mathbf{b}^{tr} = -\nabla V$

— Potential: $V = x/\pi R^2$,

— Stress tensor: $2\pi(\sigma_{rr}^{tr}, \sigma_{\theta\theta}^{tr}, \sigma_{r\theta}^{tr}) = (2\pi V, 2\pi V, 0)$,

— the solenoidal force field $\mathbf{b}^{rot} = \nabla \times \mathbf{A}$

— Vector potential: $\mathbf{A} = r^2/(\pi R^3)\mathbf{e}_z$,

— Stress tensor: $2\pi(\sigma_{rr}^{rot}, \sigma_{\theta\theta}^{rot}, \sigma_{r\theta}^{rot}) = (0, 0, \pi A)$,

— a corrective field

— Airy stress function: $\phi^C = -r^3 \sin \theta / (4\pi R^2)$,

— Stress tensor: $2\pi(\sigma_{rr}^C, \sigma_{\theta\theta}^C, \sigma_{r\theta}^C) = rR^{-2}(-\sin \theta, -3 \sin \theta, \cos \theta)$,

giving the following result for the displacement field:

$$4\pi \mu \mathcal{G}_{xx} = -(\kappa + 1) \log \frac{s}{R} + 2 \sin^2 \beta + \frac{r^2}{2R^2} (1 + 2 \cos^2 \theta) + \left( \frac{r^3}{R^3} - \left( \kappa + \frac{5}{3} \right) \frac{r}{R} \right) \cos \theta - 1 \quad (A 2a)$$

and

$$4\pi \mu \mathcal{G}_{yx} = -\sin 2\beta - (\kappa - 1)\beta + \frac{r^2}{2R^2} \sin 2\theta + \left( \frac{r^3}{R^3} - \left( \kappa + \frac{5}{3} \right) \frac{r}{R} \right) \sin \theta. \quad (A 2b)$$

This solution is represented in figures 7b and 8b.

**Remark.** In contrast to the case of a normal load, the displacement at point O is here non-zero.

### (e) Out-of-plane tangential load

We here treat the case represented in figure 3c. This case is scalar, and is simple enough to be solved in a tractable manner by separation of variables. The solution is

$$4\pi \mu \mathcal{G}_{zz} = -4 \log \frac{s}{R} + \frac{r^2}{R^2} - \frac{1}{2}. \quad (A 3)$$

This solution is represented in figures 7c and 8c. The shear stresses $\sigma_{rz} = \frac{1}{2}\mu\partial_r\mathcal{G}_{zz}$ and $\sigma_{\theta z} = \frac{1}{2}\mu\partial_\theta\mathcal{G}_{zz}/r$ can be derived using $s\partial_r s = r - R\cos\theta$ and $s\partial_\theta s = rR\sin\theta$—obtained from $s^2 = R^2 - 2rR\cos\theta + r^2$. One obtains: $2\pi\sigma_{rz} = -(r - R\cos\theta)/s^2 + r/(2R^2)$ and $2\pi\sigma_{\theta z} = -R\sin\theta/s^2$.

**Remark.** Here again, in contrast to the case of a normal load, the displacement at point O is non-zero.

### (f) Dipolar normal load

We here treat the case represented in figure 3d. This case is qualitatively different to the others as the fundamental solution has a stronger divergence at the source, and decreases to zero in the far field. In principle, approximating the contact compliance at the leading order could then be done by using the half-space approximation and relying on the 'far point' concept, as classically done for a point-wise (three-dimensional) contact compliance. Nevertheless, for the sake of a uniform presentation, we give here the complete solution of case (d).

The solution can be constructed by adding the contributions of:

— the Flamant solution
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Solutions to the Hertzian problems at the edge, for a contact half-width $a/R = 0.1\%$, a Poisson's ratio $\nu = 0.28$, a shear modulus $\mu = 1$ and unit average traction ((a) $Q_y = 1$, (b) $Q_x = 1$, (c) $Q_z = 1$, (d) $M_z = 1$ and $R = 1$). The traction distribution is taken from equation (2.6) for (a) (finite disp.), from equation (2.4) for (b) and (c), and from equation (2.5) for (d) (small incremental disp. and tilt). Solid lines stand for analytical formulae (see equation (4.12)); dashed lines represent a validation obtained using finite elements, see appendix B for details. Note that (d) is strongly zoomed to see the field at non-zero angles.</figcaption></figure>

— Airy stress function: $-R\partial_x((1/\pi)x\beta) = -(R/\pi)\beta - (R/2\pi)\sin 2\beta$—
using $\partial_x\beta = -\partial_x \arctan x/(y+R) = -\cos\beta/s$,

— Stress tensor: $2\pi(\sigma_{ss}^F, \sigma_{\beta\beta}^F, \sigma_{s\beta}^F) = 4Rs^{-2}(\sin 2\beta, 0, -\cos^2\beta)$,

— the solenoidal force field $\mathbf{b}^{\text{rot}} = \nabla \times \mathbf{A}$

— Vector potential: $\mathbf{A} = r^2/(\pi R^3)\mathbf{e}_z$,

— Stress tensor: $2\pi(\sigma_{rr}^{\text{rot}}, \sigma_{\theta\theta}^{\text{rot}}, \sigma_{r\theta}^{\text{rot}}) = (0, 0, \pi A)$,

— (no corrective field is needed),

giving the following result for the displacement field:

$$4\pi \mu \mathcal{G}_{x\alpha} = \frac{R \cos \beta}{s} (\kappa + 1 - 4 \cos^2 \beta) - (\kappa - 2) + \left( \frac{r^3}{R^3} - \left( \kappa + \frac{5}{3} \right) \frac{r}{R} \right) \cos \theta \quad (A 4a)$$

and

$$4\pi \mu \mathcal{G}_{y\alpha} = -\frac{R \sin \beta}{s} (\kappa + 1 + 4 \cos^2 \beta) + \left( \frac{r^3}{R^3} - \left( \kappa + \frac{5}{3} \right) \frac{r}{R} \right) \sin \theta. \quad (A 4b)$$

This solution is represented in figures 7d and 8d.

## Appendix B. Comparison with finite elements

Finite-elements (FE) reference computations presented in figures 5, 6 and 8 were performed using the FreeFem++ software [25]. In all cases, a mesh refinement loop was used to properly converge to the solution near the contact region. P2 Lagrange elements were used. In neither case, the geometry of the cylinder was deformed prior to applying a load (figure 8) or applying a Dirichlet condition (figures 5 and 6): even if the cases may represent a boundary deformed by an initial pressure, this geometric effect was neglected, as do the analytical approximations.

The static problems represented in figure 8a,b and c are ill behaved due to the difficulty to reach the balance of forces at a required numerical accuracy: a small penalization term was added to remedy this issue. Cases (b) and (c) involve a force distribution that is weakly singular at the tips of the contact strip: the distribution was therefore slightly truncated (and re-normalized). Case (d) involves a much stronger singularity: no static FE computation was intended, however we believe that the dynamic FE computations shown in figures 5 and 6 (forced response and eigenmodes) provide solid validation.
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Solutions to the Hertzian problems at the edge, for a contact half-width $a/R = 0.1\%$, a Poisson's ratio $\nu = 0.28$, a shear modulus $\mu = 1$ and unit average traction ((a) $Q_y = 1$, (b) $Q_x = 1$, (c) $Q_z = 1$, (d) $M_z = 1$ and $R = 1$). The traction distribution is taken from equation (2.6) for (a) (finite disp.), from equation (2.4) for (b) and (c), and from equation (2.5) for (d) (small incremental disp. and tilt). Solid lines stand for analytical formulae (see equation (4.12)); dashed lines represent a validation obtained using finite elements, see appendix B for details. Note that (d) is strongly zoomed to see the field at non-zero angles.</figcaption></figure>

— Airy stress function: $-R\partial_x((1/\pi)x\beta) = -(R/\pi)\beta - (R/2\pi)\sin 2\beta$—
using $\partial_x\beta = -\partial_x \arctan x/(y+R) = -\cos\beta/s$,

— Stress tensor: $2\pi(\sigma_{ss}^F, \sigma_{\beta\beta}^F, \sigma_{s\beta}^F) = 4Rs^{-2}(\sin 2\beta, 0, -\cos^2\beta)$,

— the solenoidal force field $\mathbf{b}^{\text{rot}} = \nabla \times \mathbf{A}$

— Vector potential: $\mathbf{A} = r^2/(\pi R^3)\mathbf{e}_z$,

— Stress tensor: $2\pi(\sigma_{rr}^{\text{rot}}, \sigma_{\theta\theta}^{\text{rot}}, \sigma_{r\theta}^{\text{rot}}) = (0, 0, \pi A)$,

— (no corrective field is needed),

giving the following result for the displacement field:

$$4\pi \mu \mathcal{G}_{x\alpha} = \frac{R \cos \beta}{s} (\kappa + 1 - 4 \cos^2 \beta) - (\kappa - 2) + \left(\frac{r^3}{R^3} - \left(\kappa + \frac{5}{3}\right) \frac{r}{R}\right) \cos \theta \quad (A 4a)$$

and

$$4\pi \mu \mathcal{G}_{y\alpha} = -\frac{R \sin \beta}{s} (\kappa + 1 + 4 \cos^2 \beta) + \left(\frac{r^3}{R^3} - \left(\kappa + \frac{5}{3}\right) \frac{r}{R}\right) \sin \theta. \quad (A 4b)$$

This solution is represented in figures 7d and 8d.

## Appendix B. Comparison with finite elements

Finite-elements (FE) reference computations presented in figures 5, 6 and 8 were performed using the FreeFem++ software [25]. In all cases, a mesh refinement loop was used to properly converge to the solution near the contact region. P2 Lagrange elements were used. In neither case, the geometry of the cylinder was deformed prior to applying a load (figure 8) or applying a Dirichlet condition (figures 5 and 6): even if the cases may represent a boundary deformed by an initial pressure, this geometric effect was neglected, as do the analytical approximations.

The static problems represented in figure 8a,b and c are ill behaved due to the difficulty to reach the balance of forces at a required numerical accuracy: a small penalization term was added to remedy this issue. Cases (b) and (c) involve a force distribution that is weakly singular at the tips of the contact strip: the distribution was therefore slightly truncated (and re-normalized). Case (d) involves a much stronger singularity: no static FE computation was intended, however we believe that the dynamic FE computations shown in figures 5 and 6 (forced response and eigenmodes) provide solid validation.
# Appendix C. Useful integrals

All formulae in this appendix are expressed with dimensionless coordinates, obtained by normalizing with the radius of the cylinder. Equivalently, it comes down to writing $R = 1$ in the notations used above.

Let us start with a few useful relations obtained by introducing complex numbers. The two sets of polar coordinates relate as

$$s e^{i\beta} = 1 - r e^{i\theta}, \qquad (C 1)$$

with $i = \sqrt{-1}$. Taking the logarithm of (C 1), expanding the right-hand side into a Laurent series, and identifying real and imaginary parts yields

$$\log s = - \sum_{n \ge 1} \frac{r^n \cos n\theta}{n} \qquad (C 2)$$

and

$$\beta = - \sum_{n \ge 1} \frac{r^n \sin n\theta}{n}. \qquad (C 3)$$

## (a) Gauges on translations

The integrals listed below relate to gauges (3.2a).

### (i) Non-zero integrals

$$- \iint_S dS = \pi.$$

$$- \iint_S r^n dS = \int_{\theta=0}^{2\pi} \int_{r=0}^{1} r^{n+1} dr d\theta = 2\pi/(n+2).$$

$$- \iint_S \sin^2 \theta dS = \iint_S \cos^2 \theta dS = \pi/2.$$

$$- \iint_S \sin^2 \beta \, dS = \pi/4: \text{ using } \iint_S dS = \int_{\beta=-\pi/2}^{\pi/2} \int_{s=0}^{2 \cos \beta} s \, ds \, d\beta.$$

$$- \iint_S \cos \beta/s \, dS = \pi: \text{ same method.}$$

$$- \iint_S \cos \beta^3/s \, dS = 3\pi/4: \text{ same method.}$$

### (ii) Zero integrals

$$- \iint_S \log s \, dS = 0: \text{ all terms of the series expansion (C 2) give 0.}$$

$$- \text{ Integrals of } r^n \cos m\theta, r^n \sin m\theta, \text{ with } m \ge 1 \text{ an integer.}$$

$$- \text{ Integrals of } \sin \beta, \sin 2\beta \text{ and } \beta: \text{ odd functions of } x.$$

## (b) Gauge on rotation

The integrals listed below relate to gauge (3.2b): $\iint_S r u_\theta \, dS = \iint_S (r u_x \cos \theta + r u_y \sin \theta) \, dS$.

### (i) Non-zero integrals

$$- \iint_S r \cos \theta \log s \, dS = -\pi/4: \text{ using the series expansion in equation (C 2), all terms give 0 except for } n=1.$$

$$- \iint_S \beta r \sin \theta \, dS = -\pi/4: \text{ same method, using equation (C 3).}$$

$$- \iint_S r \cos \theta \sin^2 \beta \, dS = \pi/12: \text{ using } r \cos \theta = 1 - s \cos \beta \text{ and } \iint_S dS = \int_{\beta=-\pi/2}^{\pi/2} \int_{s=0}^{2 \cos \beta} s \, ds \, d\beta.$$

$$- \iint_S r \cos \theta \cos \beta/s \, dS = \pi/4: \text{ same method.}$$

$$- \iint_S r \cos \theta \cos \beta^3/s \, dS = \pi/8: \text{ same method.}$$

$$- \iint_S r \sin \theta \sin 2\beta \, dS = -\pi/3: \text{ same method, using } r \sin \theta = -s \sin \beta.$$

$$- \iint_S r \sin \theta \sin \beta/s \, dS = -\pi/4: \text{ same method.}$$

$$- \iint_S r \sin \theta \cos^2 \beta \sin \beta/s \, dS = -\pi/8: \text{ same method.}$$

$$- \iint_S r^2 \, dS = \pi/2.$$

$$- \iint_S r^4 \, dS = \pi/3.$$
All formulae in this appendix are expressed with dimensionless coordinates, obtained by normalizing with the radius of the cylinder. Equivalently, it comes down to writing $R = 1$ in the notations used above.

Let us start with a few useful relations obtained by introducing complex numbers. The two sets of polar coordinates relate as

$$s e^{i\beta} = 1 - r e^{i\theta}, \qquad (C 1)$$

with $i = \sqrt{-1}$. Taking the logarithm of (C 1), expanding the right-hand side into a Laurent series, and identifying real and imaginary parts yields

$$\log s = - \sum_{n \ge 1} \frac{r^n \cos n\theta}{n} \qquad (C 2)$$

and

$$\beta = - \sum_{n \ge 1} \frac{r^n \sin n\theta}{n}. \qquad (C 3)$$

(a) Gauges on translations

The integrals listed below relate to gauges (3.2a).

(i) Non-zero integrals

$$- \iint_S dS = \pi.$$

$$- \iint_S r^n dS = \int_{\theta=0}^{2\pi} \int_{r=0}^{1} r^{n+1} dr d\theta = 2\pi/(n+2).$$

$$- \iint_S \sin^2 \theta dS = \iint_S \cos^2 \theta dS = \pi/2.$$

$$- \iint_S \sin^2 \beta \, dS = \pi/4: \text{ using } \iint_S dS = \int_{\beta=-\pi/2}^{\pi/2} \int_{s=0}^{2\cos\beta} s \, ds \, d\beta.$$

$$- \iint_S \cos \beta/s \, dS = \pi: \text{ same method.}$$

$$- \iint_S \cos \beta^3/s \, dS = 3\pi/4: \text{ same method.}$$

(ii) Zero integrals

$$- \iint_S \log s \, dS = 0: \text{ all terms of the series expansion (C 2) give 0.}$$

$$- \text{ Integrals of } r^n \cos m\theta, r^n \sin m\theta, \text{ with } m \ge 1 \text{ an integer.}$$

$$- \text{ Integrals of } \sin \beta, \sin 2\beta \text{ and } \beta: \text{ odd functions of } x.$$

(b) Gauge on rotation

The integrals listed below relate to gauge (3.2b): $\iint_S r u_\theta \, dS = \iint_S (r u_x \cos \theta + r u_y \sin \theta) \, dS$.

(i) Non-zero integrals

$$- \iint_S r \cos \theta \log s \, dS = -\pi/4: \text{ using the series expansion in equation (C 2), all terms give 0 except for } n = 1.$$

$$- \iint_S \beta r \sin \theta \, dS = -\pi/4: \text{ same method, using equation (C 3).}$$

$$- \iint_S r \cos \theta \sin^2 \beta \, dS = \pi/12: \text{ using } r \cos \theta = 1 - s \cos \beta \text{ and } \iint_S dS = \int_{\beta=-\pi/2}^{\pi/2} \int_{s=0}^{2\cos \beta} s \, ds \, d\beta.$$

$$- \iint_S r \cos \theta \cos \beta/s \, dS = \pi/4: \text{ same method.}$$

$$- \iint_S r \cos \theta \cos \beta^3/s \, dS = \pi/8: \text{ same method.}$$

$$- \iint_S r \sin \theta \sin 2\beta \, dS = -\pi/3: \text{ same method, using } r \sin \theta = -s \sin \beta.$$

$$- \iint_S r \sin \theta \sin \beta/s \, dS = -\pi/4: \text{ same method.}$$

$$- \iint_S r \sin \theta \cos^2 \beta \sin \beta/s \, dS = -\pi/8: \text{ same method.}$$

$$- \iint_S r^2 \, dS = \pi/2.$$

$$- \iint_S r^4 \, dS = \pi/3.$$
(ii) Zero integrals

— r sin θ log s, r cos θ sin 2 β, r cos θ β, and r sin θ sin² β: odd functions of x.

— r cos θ sin 2θ, r sin θ sin² θ and r sin θ cos² θ: odd functions of x.

— r sin θ sin 2θ, r cos θ sin² θ and r cos θ cos² θ: odd functions of y.

— cos θ sin θ.

(c) Convolutions with Hertzian distributions at the mid-point of contact

Most integrals are either zero (odd functions convolved by even distributions), or straightforward (sin² β = 1 is a constant on y = 0, far field terms r, θ are constant in the limit a → 0). The only tedious convolutions are those involving log s:

— For distributions (2.6) and (2.7): from [26] §4.241, equation (9):

$$
2 \int_{0}^{1} \sqrt{1-x^2} \log x \, dx = -\frac{\pi}{4} - \frac{\pi}{2} \log 2 = -\frac{\pi}{2} \log 2\sqrt{e}. \quad (C 4)
$$

— For distribution (2.4): from [26] §4.241, equation (7):

$$
\int_0^1 \frac{\log x}{\sqrt{1-x^2}} \, dx = -\frac{\pi}{2} \log 2. \qquad (C 5)
$$

References

1. Mindlin RD. 1949 Compliance of elastic bodies in contact. *J. Appl. Mech.* **16**, 259–268. (doi:10.1115/1.4009973)

2. Johnson KL. 1987 *Contact mechanics*. Cambridge, UK: Cambridge University Press.

3. Johnson KL. 1982 One hundred years of hertz contact. *Proc. Inst. Mech. Eng.* **196**, 363–378. (doi:10.1243/PIME_PROC_1982_196_039_02)

4. Kalker JJ. 1972 On elastic line contact. *J. Appl. Mech.* **39**, 1125–1132. (doi:10.1115/1.3422841)

5. Argatov I, Kachanov M, Mishuris G. 2017 On the concept of ‘far points’ in Hertz contact problems. *Int. J. Eng. Sci.* **113**, 20–36. (doi:10.1016/j.ijengsci.2016.11.009)

6. Foti F, Martinelli L. 2019 Modeling the axial-torsional response of metallic strands accounting for the deformability of the internal contact surfaces: derivation of the symmetric stiffness matrix. *Int. J. Solids Struct.* **171**, 30–46. (doi:10.1016/j.ijsolstr.2019.05.008)

7. Castillo J, Barber JR. 1997 Lateral contact of slender prismatic bodies. *Proc. R. Soc. Lond. A* **453**, 2397–2412. (doi:10.1098/rspa.1997.0128)

8. Kwun H, Bartels KA, Hanley JJ. 1998 Effects of tensile loading on the properties of elastic-wave propagation in a strand. *J. Acoust. Soc. Am.* **103**, 3370–3375. (doi:10.1121/1.423051)

9. Laguerre L, Treyssede F. 2011 Non destructive evaluation of seven-wire strands using ultrasonic guided waves. *Eur. J. Environ. Civil Eng.* **15**, 487–500. (doi:10.1080/19648189.2011.9693342)

10. Treyssède F, Frikha A, Cartraud P. 2013 Mechanical modeling of helical structures accounting for translational invariance. Part 2: guided wave propagation under axial loads. *Int. J. Solids Struct.* **50**, 1383–1393. (doi:10.1016/j.ijsolstr.2013.01.006)

11. Lai Z, Chen Q, Huang L. 2021 A semianalytical Hertzian frictional contact model in 2D. *Appl. Math. Modell.* **92**, 546–564. (doi:10.1016/j.apm.2020.11.016)

12. Han Y, Yong H, Zhou Y. 2023 The global mechanical response and local contact in multilevel helical structures under axial tension. *Int. J. Mech. Sci.* **239**, 107886. (doi:10.1016/j.ijmecsci.2022.107886)

13. Michell JH. 1899 On the direct determination of stress in an elastic solid, with application to the theory of plates. *Proc. Lond. Math. Soc.* **s1-31**, 100–124. (doi:10.1112/plms/s1-31.1.100)

14. Barber JR. 2010 *Elasticity*, 3rd edn. Dordrecht, the Netherlands: Springer.

15. Michell JH. 1900 Some elementary distributions of stress in three dimensions. *Proc. Lond. Math. Soc.* **s1-32**, 23–61. (doi:10.1112/plms/s1-32.1.23)

16. Stevenson AC. 1945 Complex potentials in two-dimensional elasticity. *Proc. R. Soc. Lond. A* **184**, 129–179. (doi:10.1098/rspa.1945.0015)
(ii) Zero integrals

— r sin θ log s, r cos θ sin 2 β, r cos θ β, and r sin θ sin² β: odd functions of x.

— r cos θ sin 2θ, r sin θ sin² θ and r sin θ cos² θ: odd functions of x.

— r sin θ sin 2θ, r cos θ sin² θ and r cos θ cos² θ: odd functions of y.

— cos θ sin θ.

(c) Convolutions with Hertzian distributions at the mid-point of contact

Most integrals are either zero (odd functions convolved by even distributions), or straightforward (sin² β = 1 is a constant on y = 0, far field terms r, θ are constant in the limit a → 0). The only tedious convolutions are those involving log s:

— For distributions (2.6) and (2.7): from [26] §4.241, equation (9):

$$
2 \int_{0}^{1} \sqrt{1-x^2} \log x \, dx = -\frac{\pi}{4} - \frac{\pi}{2} \log 2 = -\frac{\pi}{2} \log 2\sqrt{e}. \quad (C 4)
$$

— For distribution (2.4): from [26] §4.241, equation (7):

$$
\int_0^1 \frac{\log x}{\sqrt{1-x^2}} \, dx = -\frac{\pi}{2} \log 2. \qquad (C 5)
$$

References

1. Mindlin RD. 1949 Compliance of elastic bodies in contact. *J. Appl. Mech.* **16**, 259–268. (doi:10.1115/1.4009973)

2. Johnson KL. 1987 *Contact mechanics*. Cambridge, UK: Cambridge University Press.

3. Johnson KL. 1982 One hundred years of hertz contact. *Proc. Inst. Mech. Eng.* **196**, 363–378. (doi:10.1243/PIME_PROC_1982_196_039_02)

4. Kalker JJ. 1972 On elastic line contact. *J. Appl. Mech.* **39**, 1125–1132. (doi:10.1115/1.3422841)

5. Argatov I, Kachanov M, Mishuris G. 2017 On the concept of ‘far points’ in Hertz contact problems. *Int. J. Eng. Sci.* **113**, 20–36. (doi:10.1016/j.ijengsci.2016.11.009)

6. Foti F, Martinelli L. 2019 Modeling the axial-torsional response of metallic strands accounting for the deformability of the internal contact surfaces: derivation of the symmetric stiffness matrix. *Int. J. Solids Struct.* **171**, 30–46. (doi:10.1016/j.ijsolstr.2019.05.008)

7. Castillo J, Barber JR. 1997 Lateral contact of slender prismatic bodies. *Proc. R. Soc. Lond. A* **453**, 2397–2412. (doi:10.1098/rspa.1997.0128)

8. Kwun H, Bartels KA, Hanley JJ. 1998 Effects of tensile loading on the properties of elastic-wave propagation in a strand. *J. Acoust. Soc. Am.* **103**, 3370–3375. (doi:10.1121/1.423051)

9. Laguerre L, Treyssede F. 2011 Non destructive evaluation of seven-wire strands using ultrasonic guided waves. *Eur. J. Environ. Civil Eng.* **15**, 487–500. (doi:10.1080/19648189.2011.9693342)

10. Treyssède F, Frikha A, Cartraud P. 2013 Mechanical modeling of helical structures accounting for translational invariance. Part 2: guided wave propagation under axial loads. *Int. J. Solids Struct.* **50**, 1383–1393. (doi:10.1016/j.ijsolstr.2013.01.006)

11. Lai Z, Chen Q, Huang L. 2021 A semianalytical Hertzian frictional contact model in 2D. *Appl. Math. Modell.* **92**, 546–564. (doi:10.1016/j.apm.2020.11.016)

12. Han Y, Yong H, Zhou Y. 2023 The global mechanical response and local contact in multilevel helical structures under axial tension. *Int. J. Mech. Sci.* **239**, 107886. (doi:10.1016/j.ijmecsci.2022.107886)

13. Michell JH. 1899 On the direct determination of stress in an elastic solid, with application to the theory of plates. *Proc. Lond. Math. Soc.* **s1-31**, 100–124. (doi:10.1112/plms/s1-31.1.100)

14. Barber JR. 2010 *Elasticity*, 3rd edn. Dordrecht, the Netherlands: Springer.

15. Michell JH. 1900 Some elementary distributions of stress in three dimensions. *Proc. Lond. Math. Soc.* **s1-32**, 23–61. (doi:10.1112/plms/s1-32.1.23)

16. Stevenson AC. 1945 Complex potentials in two-dimensional elasticity. *Proc. R. Soc. Lond. A* **184**, 129–179. (doi:10.1098/rspa.1945.0015)
17. Meleshko V. 2003 Selected topics in the history of the two-dimensional biharmonic problem. *Appl. Mech. Rev.* **56**, 33–85. (doi:10.1115/1.1521166)

18. Ciavarella M. 1998 The generalized Cattaneo partial slip plane contact problem. I–Theory. *Int. J. Solids Struct.* **35**, 2349–2362. (doi:10.1016/S0020-7683(97)00154-6)

19. Ciavarella M. 1998 The generalized Cattaneo partial slip plane contact problem. II-examples. *Int. J. Solids Struct.* **35**, 2363–2378. (doi:10.1016/S0020-7683(97)00155-8)

20. Vakis A et al. 2018 Modeling and simulation in tribology across scales: an overview. *Tribol. Int.* **125**, 169–199. (doi:10.1016/j.triboint.2018.02.005)

21. Yan JF, Huang GY. 2019 A double-Hertz model for adhesive contact between cylinders under inclined forces. *Proc. R. Soc. A* **475**, 20180589. (doi:10.1098/rspa.2018.0589)

22. Lurie AI. 2002 *Analytical mechanics*. Foundations of Engineering Mechanics. Berlin, Heidelberg: Springer.

23. Argatov I. 2011 Response of a wire rope strand to axial and torsional loads: asymptotic modeling of the effect of interwire contact deformations. *Int. J. Solids Struct.* **48**, 1413–1423. (doi:10.1016/j.ijsolstr.2011.01.021)

24. Timoshenko S, Goodier J. 1951 *Theory of elasticity*. New York, NY: McGraw-Hill.

25. Hecht F. 2012 New development in FreeFem++. *J. Numer. Math.* **20**, 251–265. (doi:10.1515/jnum-2012-0013)

26. Gradshtein IS, Ryzhik IM. 2007 *Table of integrals, series, and products*, 7th edn. New York, NY: Academic Press.
17. Meleshko V. 2003 Selected topics in the history of the two-dimensional biharmonic problem. *Appl. Mech. Rev.* **56**, 33–85. (doi:10.1115/1.1521166)

18. Ciavarella M. 1998 The generalized Cattaneo partial slip plane contact problem. I–Theory. *Int. J. Solids Struct.* **35**, 2349–2362. (doi:10.1016/S0020-7683(97)00154-6)

19. Ciavarella M. 1998 The generalized Cattaneo partial slip plane contact problem. II-examples. *Int. J. Solids Struct.* **35**, 2363–2378. (doi:10.1016/S0020-7683(97)00155-8)

20. Vakis A et al. 2018 Modeling and simulation in tribology across scales: an overview. *Tribol. Int.* **125**, 169–199. (doi:10.1016/j.triboint.2018.02.005)

21. Yan JF, Huang GY. 2019 A double-Hertz model for adhesive contact between cylinders under inclined forces. *Proc. R. Soc. A* **475**, 20180589. (doi:10.1098/rspa.2018.0589)

22. Lurie AI. 2002 *Analytical mechanics*. Foundations of Engineering Mechanics. Berlin, Heidelberg: Springer.

23. Argatov I. 2011 Response of a wire rope strand to axial and torsional loads: asymptotic modeling of the effect of interwire contact deformations. *Int. J. Solids Struct.* **48**, 1413–1423. (doi:10.1016/j.ijsolstr.2011.01.021)

24. Timoshenko S, Goodier J. 1951 *Theory of elasticity*. New York, NY: McGraw-Hill.

25. Hecht F. 2012 New development in FreeFem++. *J. Numer. Math.* **20**, 251–265. (doi:10.1515/jnum-2012-0013)

26. Gradshtein IS, Ryzhik IM. 2007 *Table of integrals, series, and products*, 7th edn. New York, NY: Academic Press.