Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
Open AccessResearch articles

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

Pierric Mora

Pierric Mora

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

[email protected]

Contribution: Conceptualization, Formal analysis, Investigation, Methodology, Validation, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed

Fabien Treyssède

Fabien Treyssède

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

Contribution: Formal analysis, Funding acquisition, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed

Laurent Laguerre

Laurent Laguerre

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

Contribution: Formal analysis, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed


    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 compliance problem is the fact that the strains and the elastic displacements decay to 0 as the distance to contact increases—respectively as s2 and s1 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 s1, 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.

    2. Problem statement

    (a) Summary

    (i) System

    A long elastic cylinder of radius R, shear modulus μ, Poisson’s ratio ν and density ρ 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,θ,z), and contact-centred polar (s,β,z).

    Figure 1.

    Figure 1. (a and b) Scheme of the problem and of the coordinate systems.

    (ii) Boundary condition for the displacement

    Under this action, the cylinder deforms and a small contact strip of half-width aR 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 utot and σtot. The displacement satisfies the following boundary condition:

    utot=dalong the contact:r=R,aR<θ<aR,2.1
    with d a given displacement distribution, considered as known. Whether 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 q applied on the contact strip. The symmetric part of 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 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 aR, 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. We write the traction as the sum of four components:

    where qQi(x) are monopolar-like distributions while qMz(x) is a dipolar-like one, and where ei is either unit vector of the standard basis (i=x,y,z). Note that in equation (2.2) we have chosen to omit ex and ez 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 Q and moment M. Assuming qQi and qMz are properly normalized, they express as
    Figure 2.

    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).

    (iv) Problem

    The problem is to find the displacement field utot 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 aR 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 qQi to refer to a monopolar-like load on the i=x,y,z component, and qMz for a dipolar-like load around the z axis. The formulae below are normalized to aaqQidx=1 and aaxqMzdx=R to be compatible with equations (2.2) and (2.3).

    Small incremental uniform displacement:

    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 (QxqQx and QzqQz) and the nonlinear dependence of the contact half-width upon the normal component QyqQy.

    Small incremental uniform tilt

    Equation (2.5) represents a small, shearless prescribed tilt uy(x)=αx, with α=R(κ+1)/(2πμa2), assuming full-stick. This distribution is usually called the inclined flat punch (e.g. [2], ch. 2.7).

    Finite normal displacement

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

    Finite transverse displacement, partial slip

    qQx,z(x)={2π(a2c2)(a2x2c2x2)in the stick zone|x|<c,2π(a2c2)a2x2in the slip zonec<|x|<a.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=Qy the finite normal pressure and Qshear=Qx2+Qz2<fP the total shear magnitude, c is given by the relation c=a1Qshear/fP. Equation (2.7) further assumes that the normal traction is given by equation (2.6), implying no tilt from equation (2.5), or Mz=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:


    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=Θez×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

    in which SdS 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α|C=xdy|C the prescribed inclination at contact, and

    the inner field local inclination at contact.

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

    d^|C=(dx,dy,dz,Rdα)T|C: prescribed displacement at contact.

    Q^=(Qx,Qy,Qz,R1Mz)T: net efforts at contact.

    u^=(ux,uy,uz,Rα)T: inner field u^|C: inner field at contact.

    U^=(Ux,Uy,Uz,RΘ)T: skeleton field state variables.

    By using the relations utot|C=U+RΘex+u|C and xuy,tot|C=Θ+α, equation (2.1) can be re-written
    with R^ the relations matrix

    (d) Dynamic inner and outer problems

    (i) Balance of forces

    The balance of forces for the total field reads

    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 SσtotdS=Q, and Sr×σtotdS=M by noticing that r×σtot=i(r×σi,tot) because ei×σi,tot=0. We then obtain the law of conservation of momentum

    the mass matrix, m=SρdS=ρπR2 the mass and J=Sρr2dS=(1/2)mR2 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

    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
    where btr=ρm1Q can be viewed as the inertia of a uniform translation and brot=ρJ1M×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 ρt2u0 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

    with C^(a) the contact compliance matrix, which depends on the contact half-width a. For finite displacements (i.e. qQy is given by equation (2.6)), we recall that a depends on the normal pressure Qy, which renders the relation (3.11) nonlinear. For small incremental displacements (qQy 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, C^ can be approximated by a diagonal matrix because the coordinates (x,y) align with the tangential and normal directions

    Indeed, the off-diagonal terms in equation (3.12) Cxα=Rux|C/Mz and Cαx=Rα|C/Qx are much smaller than the corresponding diagonal terms and can be neglected: CxαCxx and CαxCαα.

    (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 C2 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 C1 here. Let us then formally extend equation (3.11) as

    where we have defined G^(Ci,r,θ) the far-field compliance matrix from a source at point Ci (i=1,2), u^|Ci (Q^i) the inner displacement (efforts) at point Ci, θ1,2 the angle from C1 to C2 and ai the ith contact half-width. Note that we have added the Ci argument to C^ to emphasize that it might in general not be diagonal—to preserve this property one should use polar coordinates for u^i and Q^i, which also leads to a more symmetric form for G^.

    Unlike C^, G^ does not depend on ai, however, it contains off-diagonal terms and can be approximated by

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

    (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 u. We solve these four problems analytically in the static limit, i.e. by taking ρt2u0.

    We first define four Green’s problems by considering line loads of unit amplitude (Qx,y,z=1or0, R1Mz=1or0) and call G^ their solutions (Green’s matrix). We then consider normalized distributions of finite extent (see §2b) and call C^ the matrix of solutions at the centre of the contact, i.e. for 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’s1 constant κ=34ν.

    (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.

    Figure 3.

    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.



    normal load: q=δ(rC)ey and b=(πR2)1ey,


    in-plane tangential load: q=δ(rC)ex and b=(πR2)1(ex+2rR1eθ),


    out-of-plane tangential load: q=δ(rC)ez and b=(πR2)1ez,


    dipolar normal load: q=Rxδ(rC)ey and b=2r/(πR3)eθ,

    where δ 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 β=(θπ)/2 and s=2|sinθ/2| for 0θ2π. 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.


    Normal load



    In-plane tangential load



    Out-of-plane tangential load



    Dipolar normal load


    Let us emphasize that equations (4.1)–(4.4) are valid only for 0θ2π because of the relation β(θ). For another range of definition of θ such as πθπ, they have to be adapted by replacing θθmod2π, 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 Qx,y,z=1or0, R1Mz=1or0.

    (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 C^ when the convolutions are evaluated at point C. The integrals involved are solved analytically in the limit 2aR – 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.


    Small, incremental normal load



    Small, incremental in-plane tangential load



    Small, incremental out-of-plane tangential load



    Small, incremental dipolar normal load


    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 Cxx, Cyy and Czz corresponding to distributions given by equations (2.6) and (2.7) read


    Finite normal load



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

    4πμCαx=same as equation (4.6).4.10b


    Finite out-of-plane tangential load, with partial slip


    (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

    u^(r)={G^(r)Q^far from the contact:||rC||a,C^(a)Q^at the centre of the contact:r=C,4.12
    where the far field terms ||rC||a are approximated by merely multiplying Green’s functions by the average loads Qx,y,z and Mz/R, and where the intermediate region sa 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 Qy=P per unit axial length. This case is treated in [2] §5.6. We call C1=(0,R)T and C2=(0,R)T the contact points before any pressure is applied, and a1(P) and a2(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 u1 and u2 the solutions to these two cases, given by equations (4.1) and (4.9) under an appropriate change of coordinate system. The displacement u=u1+u2 at point C1 is obtained by subtracting equation (4.1) evaluated at θ=π to equation (4.9):

    A similar expression can be obtained at point C2. The compression of the diameter δ=uy(C1)uy(C2) is then
    With κ=34ν and μ=E/2(1+ν), one recovers the formula given in [2] §5.6:

    (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.

    Figure 4.

    Figure 4. A cylinder in Hertzian contact with a vibrating soil, represented at low frequencies as a mass-spring system with 4 d.f.

    (i) Initial state and perturbed state

    A cylinder lays over a stiff soil. Under the action of a normal body force of magnitude P0 (e.g. the weight) and the reaction of the soil q0 (Hertzian contact distribution following equation (2.6)), the cylinder deforms into an initial (static) state (u0,σ0) characterized by a contact strip of width 2a, where a=4P0R(1ν2)/πE according to Hertz’s Law.

    Then, let us consider that the soil slightly vibrates harmonically owing to a small prescribed displacement deiωt, with ω the angular frequency and t the time. The related contact load qeiωt follows equations (2.4) and (2.5). This excitation produces a small incremental dynamic displacement udyneiωt. The total displacement and stress become

    For simplicity, we use a coordinate system centred on the centre of mass of the cylinder in the initial state u0, such as in the previous sections (figure 1). From here on, only the dynamic part udyn is of interest and we apply the notations defined earlier to it: U is the skeleton part of udyn, 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

    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 (ωn2,U^n) the eigenvalues and eigenvectors of the Hermitian generalized eigenvalue problem K^U^n=ωn2M^U^n and require the following normalization: U^nTM^U^n=m. Then,
    in which dn=d^nTd^|C is the participation factor of mode n, and with d^n=m1(R^T)1M^U^n. The contact forces can be obtained using equation (3.11)
    where we have defined Q^n=mωn2d^n the eigen contact forces. Knowing the contact efforts Q^ at a given frequency ω, the inner field u(r,ω) can be obtained at any point using equation (4.12). Finally, the total dynamic field udyn(r,ω) can be obtained using equation (3.1a).

    The eigenvalue problem is simple enough for analytical solutions to be derived; indeed only the Ux and Θ components are coupled. The four eigenmodes are

    with mx=(m1+(J/R2)1) 1=1/3m and mα=3/2m.

    Remark. Equations (5.9a) and (5.9d) were obtained by noticing that CxxCαα, expanding in powers of Cxx/Cαα 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=217GPa, ν=0.28, ρ=7.8gcm3, R=3.6mm. 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 applying a tension close to the maximal service load—as a comparison, the action of the weight only yields a much weaker value (a/R6.8×103%). The frequencies will be normalized with a characteristic angular frequency of the cylinder ωcyl=cS/R, with cS=E/2(1+ν)ρ the shear wave velocity.

    Figure 5 represents the forced response (equation (5.7)) at the top of the cylinder (udyn|r=R,θ=π) as a function of the normalized driving frequency ω/ω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 ω/ωcyl1.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 loga/R+const., hence resulting into a much higher compliance for a/R1.

    Figure 5.

    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%.

    Figure 6 illustrates further this dependence on the contact width. The eigen angular frequencies ωn (equations (5.9)) are represented as a function of a/R, emphasizing three families of modes: one pendulum-like contact mode with ω1a/R, then three other contact modes with ω2,3,41/loga/R+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.

    Figure 6.

    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.

    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 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.


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


    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.


    1 The formulae are also valid for plane stress by taking κ=(3ν)/(1+ν).

    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/||rC|| rather than log||rC||: 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.
    Figure 7.

    Figure 7. (ad) Analytical Green’s functions: displacement field for a Poisson’s ratio ν=0.28. Log scale colours are proportional to the magnitude of the displacement.

    (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=||rC|| and tanβ=rsinθ/(Rrcosθ), or equivalently x=ssinβ and y+R=scosβ, with π/2βπ/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 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: ϕF=π1sβsinβ,

    Stress tensor: 2π(σssF,σββF,σsβF)=(4cosβ/s,0,0),

    the irrotational force field btr=V

    Potential: V=y/πR2,

    Stress tensor: 2π(σrrtr,σθθtr,σrθtr)=(2πV,2πV,0),

    the following two corrective fields:


    Airy stress function: ϕC1=r2/(4πR),


    Stress tensor: 2π(σrrC1,σθθC1,σrθC1)=(R1,R1,0),


    Airy stress function: ϕC2=r3cosθ/(4πR2),


    Stress tensor: 2π(σrrC2,σθθC2,σrθC2)=rR2(cosθ,3cosθ,sinθ),

    giving the following result for the displacement field:
    4πμGxy=sin2β+(κ1)β+r2R(κ1)sinθ+r22R2sin2θA 1a
    4πμGyy=(κ+1)logsR2sin2βr2R(κ1)cosθ+r22R2(1+2sin2θ).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: ϕF=π1sβcosβ,

    Stress tensor: 2π(σssF,σββF,σsβF)=(4sinβ/s,0,0),

    the irrotational force field btr=V

    Potential: V=x/πR2,

    Stress tensor: 2π(σrrtr,σθθtr,σrθtr)=(2πV,2πV,0),

    the solenoidal force field brot=×A

    Vector potential: A=r2/(πR3)ez,

    Stress tensor: 2π(σrrrot,σθθrot,σrθrot)=(0,0,πA),

    a corrective field

    Airy stress function: ϕC=r3sinθ/(4πR2),

    Stress tensor: 2π(σrrC,σθθC,σrθC)=rR2(sinθ,3sinθ,cosθ),

    giving the following result for the displacement field:
    4πμGxx=(κ+1)logsR+2sin2β+r22R2(1+2cos2θ)+(r3R3(κ+53)rR)cosθ1A 2a
    4πμGyx=sin2β(κ1)β+r22R2sin2θ+(r3R3(κ+53)rR)sinθ.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πμGzz=4logsR+r2R212.A 3
    This solution is represented in figures 7c and 8c. The shear stresses σrz=12μrGzz and σθz=12μθGzz/r can be derived using srs=rRcosθ and sθs=rRsinθ—obtained from s2=R22rRcosθ+r2. One obtains: 2πσrz=(rRcosθ)/s2+r/(2R2) and 2πσθz=Rsinθ/s2.

    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

    Airy stress function: Rx((1/π)xβ)=(R/π)β(R/2π)sin2β—using xβ=xarctanx/(y+R)=cosβ/s,

    Stress tensor: 2π(σssF,σββF,σsβF)=4Rs2(sin2β,0,cos2β),

    the solenoidal force field brot=×A

    Vector potential: A=r2/(πR3)ez,

    Stress tensor: 2π(σrrrot,σθθrot,σrθrot)=(0,0,πA),

    (no corrective field is needed),

    giving the following result for the displacement field:
    4πμGxα=Rcosβs(κ+14cos2β)(κ2)+(r3R3(κ+53)rR)cosθA 4a
    4πμGyα=Rsinβs(κ+1+4cos2β)+(r3R3(κ+53)rR)sinθ.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.

    Figure 8.

    Figure 8. Solutions to the Hertzian problems at the edge, for a contact half-width a/R=0.1%, a Poisson’s ratio ν=0.28, a shear modulus μ=1 and unit average traction ((a) Qy=1, (b) Qx=1, (c) Qz=1, (d) Mz=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.

    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

    seiβ=1reiθ,C 1
    with i=1. Taking the logarithm of (C 1), expanding the right-hand side into a Laurent series, and identifying real and imaginary parts yields
    logs=n1rncosnθnC 2
    β=n1rnsinnθn.C 3

    (a) Gauges on translations

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

    (i) Non-zero integrals




    Ssin2βdS=π/4: using SdS=β=π/2π/2s=02cosβsdsdβ.

    Scosβ/sdS=π: same method.

    Scosβ3/sdS=3π/4: same method.

    (ii) Zero integrals

    SlogsdS=0: all terms of the series expansion (C 2) give 0.

    Integrals of rncosmθ, rnsinmθ, with m1 an integer.

    Integrals of sinβ, sin2β and β: odd functions of x.

    (b) Gauge on rotation

    The integrals listed below relate to gauge (3.2b): SruθdS=S(ruxcosθ+ruysinθ)dS.

    (i) Non-zero integrals

    SrcosθlogsdS=π/4: using the series expansion in equation (C 2), all terms give 0 except for n=1.

    SβrsinθdS=π/4: same method, using equation (C 3).

    Srcosθsin2βdS=π/12: using rcosθ=1scosβ and SdS=β=π/2π/2s=02cosβsdsdβ.

    Srcosθcosβ/sdS=π/4: same method.

    Srcosθcosβ3/sdS=π/8: same method.

    Srsinθsin2βdS=π/3: same method, using rsinθ=ssinβ.

    Srsinθsinβ/sdS=π/4: same method.

    Srsinθcos2βsinβ/sdS=π/8: same method.



    (ii) Zero integrals

    rsinθlogs, rcosθsin2β, rcosθβ, and rsinθsin2β: odd functions of x.

    rcosθsin2θ, rsinθsin2θ and rsinθcos2θ: odd functions of x.

    rsinθsin2θ, rcosθsin2θ and rcosθcos2θ: odd functions of y.


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

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

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

    2011x2logxdx=π4π2log2=π2log2e.C 4

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

    01logx1x2dx=π2log2.C 5

    Published by the Royal Society under the terms of the Creative Commons Attribution License, which permits unrestricted use, provided the original author and source are credited.