Thick-Walled Cylinder Stress Analysis — Analytical Pipeline (Croccolo 2009)

The classical Lamé equations handle a single thick-walled cylinder under internal and external pressure. Real shaft–hub assemblies are rarely that simple: multiple coaxial rings of different materials, interference fits at each interface, and a rotating load case that adds centrifugal body forces to the stress field. Croccolo & Vincenzi (2009) derived a closed-form analytical solution that covers exactly this configuration — nn n isotropic coaxial cylinders with arbitrary interference and rotation — and validated it against finite element analysis to within 2% on all stress components.

This tool implements that solution as a seven-node calculation pipeline. The pipeline takes geometry, material properties, interferences, boundary pressures, and angular velocity as input, solves the tridiagonal compatibility system at the interfaces, and returns contact pressures and full stress distributions \(\sigma_r(r)\) , \(\sigma_\theta(r)\) , \(\sigma_{VM}(r)\) across all elements. A Python notebook ready to run on Google Colab is provided alongside the derivation, so the result of every intermediate step can be verified by hand or extended to a new case.

Quick Example

Three coaxial cylinders — steel outer ring, aluminium middle ring, steel inner shaft — assembled with interference fits \(Z_{tot,1} = 0.025\,\text{mm}\) and \(Z_{tot,2} = 0.010\,\text{mm}\) , subject to external pressure \(p_0 = 150\,\text{MPa}\) , bore pressure \(p_3 = 30\,\text{MPa}\) , and angular velocity \(\omega = 420\,\text{rad/s}\) .

Contact pressures: \(p_1 = 161.4\,\text{MPa}\) , \(p_2 = 175.0\,\text{MPa}\)

All equivalent stresses remain below typical yield strengths (\(\sigma_{y,\text{steel}} \approx 350\,\text{MPa}\) for structural steel, \(\sigma_{y,\text{Al}} \approx 280\,\text{MPa}\) for 6061-T6), except element 3 inner surface (\(\sigma_{VM} = 537.7\,\text{MPa}\) ) which exceeds the yield limit — redesign needed.

Pipeline summary:

Pipeline Overview

The pipeline solves a tridiagonal compatibility system (one equation per interface) reduced to bidiagonal form by the forward sweep (Node 4) and resolved back-to-front (Node 5). Lamé equations and centrifugal body-force corrections (Node 6) are then applied element by element.

Node 1 — Geometric Ratios

Purpose. Express all radial positions as dimensionless ratios relative to the outer radius of each element.

For a solid shaft (\(r_n = 0\) ): \(Q_n = 0\) .

Node 2 — Compatibility Coefficients K, L, N

Purpose. Assemble the coefficients of the displacement-compatibility equation at each of the \(n-1\) interfaces.

Special case: if \(Q_{i+1} = 0\) (solid shaft bore), \(N_{i+1} = 0\) and the second bracket of \(K_i\) reduces to \((1 - \nu_{i+1}) / E_{i+1}\) .

Node 3 — Dynamic Compatibility Contribution \(V_i\)

Purpose. Account for the centrifugal body-force term that shifts the effective interference at each interface when the assembly rotates.

Unit note: with \(\rho\) in kg/m³, \(\omega\) in rad/s, \(r\) in mm, and \(E\) in MPa, multiply the numerator by \(10^{-12}\) to obtain \(V_i\) in the same units as \(Z_{tot}/r\) (dimensionless). \(V_i = 0\) for all interfaces when \(\omega = 0\) .

Node 4 — Forward Sweep: H, I, J

Purpose. Reduce the \((n-1) \times (n-1)\) tridiagonal compatibility system to upper-bidiagonal form by propagating boundary conditions from the outer surface inward.

Initialisation (\(i = 1\) ):

Recursion (\(i = 2, \ldots, n-1\) ):

Node 5 — Interface Pressures

Purpose. Solve the reduced bidiagonal system by back-substitution from the innermost interface outward.

Node 6 — Stress Distributions

Purpose. Compute \(\sigma_r\) , \(\sigma_\theta\) , and \(\sigma_{VM}\) at any radius within each element by superimposing static Lamé and centrifugal contributions.

Static Lamé (eqs. 2):

Centrifugal (eq. 25), with units as in Node 3:

Von Mises (plane stress, \(\sigma_z = 0\) ):

Numerical Case — Croccolo & Vincenzi (2009), Tables 1–3

Input

\(Z_{tot,1} = 0.025\,\text{mm}\) , \(Z_{tot,2} = 0.010\,\text{mm}\) , \(p_0 = 150\,\text{MPa}\) , \(p_3 = 30\,\text{MPa}\) , \(\omega = 420\,\text{rad/s}\) .

Node 1 — Geometric ratios

Node 2 — Coefficients at interface 1 (\(i=1\) )

At interface 2 (\(i=2\) ):

Node 3 — Dynamic contributions

For interface 1, with \(\rho_2 = 2770\) , \(r_1 = 24.988\,\text{mm}\) , \(\rho_1 = 7850\) , \(r_0 = 35.000\,\text{mm}\) :

Node 4 — Forward sweep

Initialisation:

Recursion at \(i=2\) :

Node 5 — Interface pressures

Node 6 — Stresses at element boundaries

All values in MPa. For brevity only boundary radii (\(r = r_{int}\) and \(r = r_{ext}\) ) are shown.

FEM validation (Table 3 of the paper): maximum deviation on all stress components is below 2%.

Stress Distribution

The curves are continuous across interfaces. The sharp change in slope at each interface reflects the material discontinuity. The hoop stress in element 3 is dominated by the centrifugal term at the inner bore, reaching −552 MPa — the governing failure mode for this configuration.

How to Use the Notebook

The Python notebook implements the full 7-node pipeline and runs directly on Google Colab with no installation required (numpy and matplotlib are preinstalled).

Open in Colab

Modify only the cell marked ✏️ Input. The parameters are:

Known limitations. The model assumes isotropic linear-elastic materials (\(\sigma_{VM} < \sigma_y\) must be verified after the fact), axial symmetry, plane stress (\(\sigma_z = 0\) , accurate for short couplings with \(L \sim D\) ), identical axial length for all elements, and a single common \(\omega\) . No thermal loading, friction, or axial force is accounted for.

References

[1] Croccolo, D. & Vincenzi, N. (2009). A generalized theory for shaft–hub couplings. Proc. IMechE Part C: J. Mechanical Engineering Science, 223, 2231–2239.

[2] Croccolo, D. & De Agostinis, M. (2013). Analytical solution of stress and strain distributions in press fitted orthotropic cylinders. Int. J. Mechanical Sciences, 71, 21–29.