# 2D analytical solution of a multilayered hollow cylinder subjected to non-uniform loads

Vincent Renaud

► To cite this version:

Vincent Renaud. 2D analytical solution of a multilayered hollow cylinder subjected to non-uniform loads. European Journal of Mechanics - A/Solids, 2026, 116, pp.105926. 10.1016/j.euromechsol.2025.105926. ineris-05339306

HAL Id: ineris-05339306

https://ineris.hal.science/ineris-05339306v1

Submitted on 3 Nov 2025

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
Full Length Article

2D analytical solution of a multilayered hollow cylinder subjected to
non-uniform loads

Vincent Renaud

Institut national de L'environnement Industriel et des Risques (Ineris) - ARTEM, 54042, Nancy, France

ARTICLE INFO

Keywords:

Hollow cylinder

Multilayer

Analytical solution

Non-uniform loading

A B S T R A C T

In this paper, a new analytical solution has been developed using the recursive method to assess stresses, strains and displacements in the cross-section of an elastic multilayered hollow cylinder subjected to non-uniform internal and external loads. Using Michell's analysis of the stress function and Hooke's law, we derived the analytical solution and compared it to the numerical results obtained using FDM software. The analytical and numerical results are in perfect agreement and lead to the conclusion that the non-uniformity of the loading significantly influences the amplitudes of the orthradial stresses. These orthradial stresses may not always be equal to the maximum principal stress, depending on the position within the layer. This result suggests that it is not possible to find an analytical solution in the elastoplastic framework where the fracture criterion is expressed in terms of principal stresses, σ₁ = σrr or σθθ, depending on the position considered. The present study makes it possible to deduce the non-uniform loading at the interface of 2 layers by minimizing the deviations with stress and/or strain measurements in one of the layers of the hollow cylinder. This study should also help to understand the influence of a low stiffness material interspersed between rigid materials.

1. Introduction

The stability of large cylindrical structures is of interest in the mining and geotechnical industries. The analysis of a lining loaded with a solid is essential not only to reduce manufacturing costs, but above all to guarantee the stability of the structures. In order to better quantify the phenomena, a better understanding of the mechanical behavior of the lining is essential.

The analytical solution developed in this article is part of the general problem of a stressed hollow structure (the load is often purely me- chanical). The geometry of this structure (lining, support) may have central or planar symmetries, while the loading may or may not be assumed to be uniform (symmetrical). In addition, boundary conditions can be taken at infinity or on the edge of the structure and the host medium may or may not be isotropic (Hefny and Lo, 1999). Following Lamé's solution (Lamé, 1852) addressing a hollow cylinder under uni- form load, and Kirsch's work (Kirsch, 1898) on stresses around a circular hole in an infinite isotropic plate subjected to uniaxial stress, authors have proposed solutions for more complex cases. These extended cases consider factors such as thermal coupling (Yeo et al., 2017), plasticity under uniform loading (Sharan, 2005; Park and Kim, 2006; Masoudian and Hashemi, 2016), viscoelasticity (Van Sambeek, 1986), contact

interaction (Qi et al., 2023) and radial variation of mechanical proper-
ties, whether in a continuous or multilayered structure (Hongjun et al.,
2006; Tokovyy and Ma, 2018). Since the 1950s, authors have proposed
solutions for noncircular holes using the theory of complex variables
(Savin, 1951; Savin and Guz, 1964). More recently, Li et al. (2022)
published an elastoplastic analytical solution with biaxial loading.

In some innovative cases, such as the compressive monoblock
voussoir (Zghondi et al., 2018; Stasiak, 2019), measuring the load on the
inner material (concrete) in a reinforced concrete voussoir proves to be
challenging for at least 2 reasons:

1) because stress sensors are large, they cannot be placed at the interface between concrete and the compressible material. Moreover, they are less reliable than vibrating-wire sensors that measure strains; 2) in some particular cases linked to material properties, the experimental behavior of the compressible material (with a load applied to its extrados) at the cylinder scale is sufficiently heterogeneous to induce a non-uniform loading of the inner concrete voussoirs.

Modelling the actual behavior of the compressible material in a
hollow cylinder subjected to an external load with a continuous model
cannot give satisfactory results, even under uniform loading. On the
other hand, leveraging readings from various sensors strategically
placed in the voussoirs allows the analytical solution for the cross-

<footnote>E-mail address: vincent.renaud@ineris.fr.</footnote>
section of an isotropic elastic hollow cylinder under non-uniform loading to approximate the load at the interface between concrete and the compressible material (material discontinuity). This is an inexpensive and practical way to study and improve our understanding of the behavior of segmental prototypes.

There is an advantage in accelerating numerical calculations by initializing the stresses and displacements using the exact solution to the relevant equations. This acceleration is particularly crucial because the significant difference in Young's moduli between two adjacent layers poses a challenge for some computational software. Moreover, these analytical solutions are a useful support for comparing theoretical results with the numerical results provided by different methods (FDM, FEM, etc.) and software.

In the literature on hollow cylinders, 3 theoretical frameworks are used: thin ring theory, Winkler-Bach theory for curved beams, and elasticity theory. As the subject of this article concerns rigid and thick cylindrical structures, we have adopted the theory of elasticity, the other 2 theories being more suitable for thin-walled cylinders. The elastic analytical solution has an obvious interest for the reinforced concrete voussoirs sized to avoid breakage.

For all cases involving a hollow cylindrical structure whose behaviour is elastic (or in a state comparable to elasticity), analytical solutions for cases involving non-uniform internal and/or external pressure have not yet been developed. This is probably because common applications often involve uniform loading, otherwise they assume elliptical loading as an approximation of non-uniform loading. These simple solutions are found in all rock or soil mechanics books (Jaeger and Cook, 1968; Timoshenko and Goodier, 1951; Barber, 2010; Sadd, 2014).

Beyond the mechanical response of the voussoir, it can also be useful to consider other materials enveloping the voussoirs, and this requires extending analytical solutions from the monolayered cylinder to the multilayered cylinder. While some authors have provided analytical solutions for multilayered hollow cylinders, these developments are constrained by limitations. For instance, Qiu and Zhou (2016) focused on boundary conditions under uniform loading, specifically addressing assembly interference in parts. In a related context, Atangana Nkene et al. (2018) aimed to optimize the production costs of rotating cylindrical components and Khechai et al. (2023) developed an analytical solution for functionally graded composite laminated plates with circular cutouts under various in-plane loading conditions. However, these existing analyses do not comprehensively address the broader range of complexities that may be encountered in diverse loading conditions and structural configurations of multilayered hollow cylinders.

The main objective of this work is therefore to find analytical solutions for $\sigma$, $\epsilon$ and $u$ for the case of a multilayered elastic cylindrical hollow cylinder under non-uniform loading (internal and external, radial and tangential). The exact theoretical relations are then numerically compared with finite difference calculations.

## 2. Analytical approaches

### 2.1. Single-layer case

#### 2.1.1. Solutions of the biharmonic equation

According to elasticity theory, if the volume forces are harmonic or zero, the stress function (or Airy function) satisfies the biharmonic equation (which reflects both the static admissibility of the stress field and the kinematics admissibility of the displacement field): $\Delta^2\Phi = 0$.

It is a homogeneous partial differential equation of order 4, which is written more naturally in polar coordinates for cylindrical structures (Fig. 1):

<figure><img src="image_2.png" /><figcaption>Fig. 1. Polar coordinates for the hollow cylinder.</figcaption></figure>

$$ \begin{aligned} \Delta^2 \Phi(r, \theta) &= \nabla^4 \Phi(r, \theta) = \frac{\partial^4 \Phi}{\partial r^4} + \frac{2}{r^2} \frac{\partial^4 \Phi}{\partial r^2 \partial \theta^2} + \frac{1}{r^4} \frac{\partial^4 \Phi}{\partial \theta^4} + \frac{2}{r} \frac{\partial^3 \Phi}{\partial r^3} - \frac{2}{r^3} \frac{\partial^3 \Phi}{\partial r \partial \theta^2} \\ & \quad - \frac{1}{r^2} \frac{\partial^2 \Phi}{\partial r^2} + \frac{4}{r^4} \frac{\partial^2 \Phi}{\partial \theta^2} + \frac{1}{r^3} \frac{\partial \Phi}{\partial r} \\ &= 0 \end{aligned} $$

The complete solution of this equation is not known. The most common approach is by separation of variables, which postulates that $\Phi(r,\theta) = f(r)g(\theta)$. However, one cannot take it for granted that $\Phi(r,\theta)$ can always be decoupled in this way.

Michell (1899) proposed a solution written in the form of a Fourier series in $\theta$, and this approach was later taken up and improved by Obert and Duvall (1967):

$$ \begin{aligned} \Phi &= A_0 r^2 + B_0 r^2 \ln r + C_0 \ln r + D_0 + (E_0 r^2 + F_0 r^2 \ln r + G_0 \ln r \\ & \quad + H_0)\theta + \left(A_1 r + \frac{B_1}{r} + C_1 r^3 + D_1 r \ln r\right) \sin \theta \\ & \quad + \left(E_1 r + \frac{F_1}{r} + G_1 r^3 + H_1 r \ln r\right) \cos \theta + (J_1 r + K_1 r \ln r)\theta \sin \theta \\ & \quad + (L_1 r + M_1 r \ln r)\theta \cos \theta + \end{aligned} \tag{1} $$

$$ \begin{aligned} & \sum_{n=2}^{\infty} \left( A_n r^n + \frac{B_n}{r^n} + C_n r^{2-n} + D_n r^{n+2} \right) \sin n\theta \\ & + \sum_{n=2}^{\infty} \left( E_n r^n + \frac{F_n}{r^n} + G_n r^{2-n} + H_n r^{n+2} \right) \cos n\theta \end{aligned} $$

Because Michell sought a solution in the form of a Fourier series (i.e.: $g(\theta) = \sin n\theta$ or $\cos n\theta$), other permissible forms of $g(\theta)$ were not examined. To address these omissions, Michell's general solution can be expanded using the following functions which are also solutions of the biharmonic equation:

$$ \Phi_1(r, \theta) = e^{in\theta} \left[ (P_1 + Q_1 r^2) \sin(\ln r^n) + (R_1 + S_1 r^2) \cos(\ln r^n) \right] $$

$$ \Phi_2(r, \theta) = r^n e^{\theta} [P_2 \sin(\ln r - n\theta) + Q_2 \sin(\ln r + (2-n)\theta) \\ + R_2 \cos(\ln r - n\theta) + S_2 \cos(\ln r + (2-n)\theta)] $$

$$ \Phi_3(r, \theta) = re^{n\theta} [P_3 \sin(\ln r^n + \theta) + Q_3 \sin(\ln r^n - \theta) \\ + R_3 \cos(\ln r^n + \theta) + S_3 \cos(\ln r^n - \theta)] $$
Moreover, it should be noted that the terms in sin $n\theta$ and cos $n\theta$ of (1) do not require that $n$ be an integer. Indeed, the Airy function which is the solution for the problem of a crack in an infinite medium is of the form:

$$ \Phi_4(r, \theta) = r^3(A \cos \theta / 2 + B \sin \theta / 2 + C \cos 3\theta / 2 + D \sin 3\theta / 2) \quad (\text{with } \Delta^2\Phi_4 = 0). $$

However, these last 4 functions have the disadvantage of being multivalued at $\theta = 2k\pi$ both for $\Phi$ and the stresses themselves (which derive from $\Phi$). Thus, for the analysis of a complete annular structure (i. e. without boundary conditions in $\theta$), the Airy function to consider is solution (1) with $n$ being an integer. We also note that for this type of structure, the Airy function is always defined since the inner radius of the hollow cylinder is a positive real: the singularity $r=0$ does not have to be considered.

### 2.1.2. Stress tensor in polar coordinates

For problems without axial symmetry (as a result of applying a non-uniform load), the stress components are defined in polar coordinates as

$$ \begin{aligned} \sigma_{\theta\theta} = & 2A_0 + 3B_0 - \frac{C_0}{r^2} + \left(\frac{2}{r^3}F_1 + \frac{H_1}{r} + 6G_1r\right) \cos \theta \\ & + \left(\frac{2}{r^3}B_1 + \frac{D_1}{r} + 6C_1r\right) \sin \theta + \\ & \ln r(2B_0 + 2F_0\theta) + \theta \left(2E_0 + 3F_0 - \frac{G_0}{r^2} + \frac{M_1 \cos \theta + K_1 \sin \theta}{r}\right) \\ & + \sum_{n=2}^{\infty} \sin n\theta \left((n-1)nr^{n-2}A_n + n(n+1)\frac{B_n}{r^{n+2}} + (n-2)(n-1)\frac{C_n}{r^n}\right) \\ & + (n+1)(n+2)r^n D_n \right) + \\ & \sum_{n=2}^{\infty} \cos n\theta \left((n-1)nr^{n-2}E_n + n(n+1)\frac{F_n}{r^{n+2}} + (n-2)(n-1)\frac{G_n}{r^n}\right) \\ & + \frac{F_n}{r^{n+2}} + (n-2)(n-1)\frac{G_n}{r^n} + (n+2)(n+1)r^n H_n \end{aligned} \quad (3) $$

$$ \begin{aligned} \sigma_{r\theta} = & \frac{H_0 - G_0}{r^2} - E_0 - F_0 - \left(\frac{D_1 + M_1}{r} - \frac{2}{r^3}B_1 + 2C_1r\right) \cos \theta + \left(\frac{H_1 - K_1}{r} - \frac{2}{r^3}F_1 + 2G_1r\right) \sin \theta + \ln r \left(\frac{G_0}{r^2} - F_0\right) \\ & + \theta \left(\frac{M_1 \sin \theta - K_1 \cos \theta}{r}\right) + \sum_{n=2}^{\infty} n \cos n\theta \left(-(n-1)r^{n-2}A_n + (n+1)\frac{B_n}{r^{n+2}} + (n-1)\frac{C_n}{r^n} - (n+1)r^n D_n\right) \\ & - \sum_{n=2}^{\infty} n \sin n\theta \left(-(n-1)r^{n-2}E_n + (n+1)\frac{F_n}{r^{n+2}} + (n-1)\frac{G_n}{r^n} - (n+1)r^n H_n\right) \end{aligned} \quad (4) $$

follows:

$$ \sigma_{rr} = \frac{1}{r} \frac{\partial \Phi}{\partial r} + \frac{1}{r^2} \frac{\partial^2 \Phi}{\partial \theta^2} \quad \sigma_{\theta\theta} = \frac{\partial^2 \Phi}{\partial r^2} \quad \sigma_{r\theta} = \frac{1}{r^2} \frac{\partial \Phi}{\partial \theta} - \frac{1}{r} \frac{\partial^2 \Phi}{\partial r \partial \theta} $$

It is therefore possible to calculate the expressions for these 3 components (Table 8.1 of Barber, 2010) by noting that these stresses decompose into an axisymmetric part, and a "deviatoric" (non-uniform) part that corresponds to the terms that depend on $\theta$:

$$ \begin{aligned} \sigma_{rr} = & \frac{C_0}{r^2} + 2A_0 + B_0(1 + 2 \ln r) + \left(\frac{H_1 + 2J_1 + 2 \ln r K_1}{r} - \frac{2}{r^3} F_1 + 2G_1 r\right) \\ & \cos \theta + \left(\frac{D_1 - 2L_1 - 2 \ln r M_1}{r} - \frac{2}{r^3} B_1 + 2C_1 r\right) \sin \theta + \theta \\ & \left(\frac{G_0}{r^2} + 2E_0 + F_0(1 + 2 \ln r) + \frac{K_1 \sin \theta + M_1 \cos \theta}{r}\right) \\ & - \sum_{n=2}^{\infty} \sin n\theta \left((n-1)nr^{n-2}A_n + n(n+1)\frac{B_n}{r^{n+2}} + (n+2)(n-1)\frac{C_n}{r^n}\right) \\ & + (n-2)(n+1)r^n D_n \right) - \sum_{n=2}^{\infty} \cos n\theta \left((n-1)nr^{n-2}E_n + n(n+1)\frac{F_n}{r^{n+2}} + (n-2)(n+1)r^n H_n\right) \end{aligned} \quad (2) $$

Under plane strain conditions without axial symmetry (as a result of non-uniform loading), the orthoradial displacement can be written as:

$$ \begin{aligned} u_{\theta} = & \frac{1+\nu}{E} \int \left( r((1-\nu)\sigma_{\theta\theta} - \nu\sigma_{rr}) - \int ((1-\nu)\sigma_{rr} - \nu\sigma_{\theta\theta})dr \right) d\theta \\ = & \frac{1+\nu}{E} \left( 2r(1-\nu)(2\theta B_0 + \theta^2 F_0) + (\ln r - 1)(1-2\nu)(K_1\theta \cos\theta - \theta \sin\theta) \right) \\ & + Z_1(r)\cos\theta + Z_2(r)\sin\theta + Z_3(r, n)\cos n\theta + Z_4(r, n)\sin n\theta \end{aligned} \quad (5) $$

The terms are those that prevent equations (2)–(5) from being single-valued.

At this stage, it is not useful to pursue the 4 functions $Z_1$ to $Z_4$ in more detail.

The studied hollow cylinder is a problem with 4 boundary conditions: 2 stress boundary conditions (those with a radial component: $\sigma_{r\theta}$ and $\sigma_{rr}$) applied to each of the 2 radial limits (internal and external, Fig. 2). To avoid overly complex derivations, only radial loading is considered in the following ($f_3(\theta) = f_4(\theta) = 0$). Derivations considering tangential loads (internal and/or external) can be consulted in the appendix.

The 2 non-zero conditions are the internal and external radial stresses (with loading functions: $f_1(\theta)$ and $f_2(\theta)$). The expressions for the stresses, (2) through (4), are dependent on the boundary conditions of the problem. After defining these boundary conditions, it is then possible to simplify the Airy function (1) and the stresses (2) to (4) derived from it.

### 2.1.3. Simplifying Airy's function

Many simplifications can be made to the Airy function ($\Phi$) by examining the expressions for the stresses: $\sigma_{rr}$, $\sigma_{\theta\theta}$, and $\sigma_{r\theta}$.

First, we note that the constants $A_1$, $E_1$, and $D_0$ do not intervene in equations (2)–(4). They can therefore be eliminated from $\Phi$.

After that, and for a complete hollow cylinder, equations (2)–(5) must be single-valued at $\theta = 2k\pi$. The constants associated with the terms in $\theta$ and $\theta \cos \theta$ in expressions (2) to (5) should therefore be eliminated to avoid $f(\theta) \neq f(\theta+2k\pi)$. This pertains to the constants: $E_0 =$
$$F_0 = G_0 = M_1 = K_1 = B_0 = 0.$$

In expression (4), only one constant remains. Moreover, $\sigma_{r\theta}$ must be zero at $r = r_1$ or $r_2$, and this implies that $H_0 = 0$.

After these first simplifications, the expressions for $\sigma_{rr}$ and $\sigma_{r\theta}$ become:

$$\sigma_{rr} = \dots + 2 \frac{(r^2 - r_1^2)(r^2 - r_2^2)}{r^3} (G_1 \cos \theta + C_1 \sin \theta) + \frac{2}{r} (J_1 \cos \theta - L_1 \sin \theta) + \dots$$

$$\sigma_{rr} = \frac{C_0}{r^2} + 2A_0 + \left( \frac{H_1 + 2J_1}{r} - \frac{2}{r^3} F_1 + 2G_1 r \right) \cos \theta + \left( \frac{D_1 - 2L_1}{r} - \frac{2}{r^3} B_1 + 2C_1 r \right) \sin \theta \\ - \sum_{n=2}^{\infty} \sin n\theta \left( (n-1)nr^{n-2}A_n + (n+1)\frac{B_n}{r^{2+n}} + (n+2)(n-1)\frac{C_n}{r^n} + (n-2)(n+1)r^nD_n \right) \\ - \sum_{n=2}^{\infty} \cos n\theta \left( (n-1)nr^{n-2}E_n + (n+1)\frac{F_n}{r^{2+n}} + (n+2)(n-1)\frac{G_n}{r^n} + (n-2)(n+1)r^nH_n \right)$$

$$\sigma_{r\theta} = - \left( \frac{D_1}{r} - \frac{2}{r^3} B_1 + 2C_1 r \right) \cos \theta + \left( \frac{H_1}{r} - \frac{2}{r^3} F_1 + 2G_1 r \right) \sin \theta \\ + \sum_{n=2}^{\infty} n \cos n\theta \left( -(n-1)r^{n-2}A_n + (n+1)\frac{B_n}{r^{2+n}} + (n-1)\frac{C_n}{r^n} - (n+1)r^nD_n \right) \\ - \sum_{n=2}^{\infty} n \sin n\theta \left( -(n-1)r^{n-2}E_n + (n+1)\frac{F_n}{r^{2+n}} + (n-1)\frac{G_n}{r^n} - (n+1)r^nH_n \right) \quad (6)$$

Finally, for the terms in $\cos\theta$ and $\sin\theta$, the 2 boundary conditions $\sigma_{r\theta} = 0$ for $r = r_1$ and $r_2$ result in $\frac{H_1}{r_1} - \frac{2}{r_1^3}F_1 + 2G_1r_1 = 0$, $\frac{H_1}{r_2} - \frac{2}{r_2^3}F_1 + 2G_1r_2 = 0$, $\frac{D_1}{r_1} - \frac{2}{r_1^3}B_1 + 2C_1r_1 = 0$ and $\frac{D_1}{r_2} - \frac{2}{r_2^3}B_1 + 2C_1r_2 = 0$. This is only possible if

$$\sigma_{r\theta} = 2 \frac{(r^2 - r_1^2)(r^2 - r_2^2)}{r^3} (G_1 \sin \theta - C_1 \cos \theta) + \dots \quad (7)$$

The terms multiplying $G_1$ and $C_1$ ("multiplicative factors") are therefore zero at the boundary conditions (at $r = r_1$ and $r = r_2$): they can then take any value without affecting the loads at $r = r_1$ and $r = r_2$. This means that the same load (internal and/or external), can correspond to an infinity of stress fields: a situation to be avoided by setting to zero 6 constants: $B_1 = C_1 = D_1 = F_1 = G_1 = H_1 = 0$. This leaves the constants $J_1$ and $L_1$, related to the coefficients of $\cos\theta$ and $\sin\theta$ (loading functions):

$$2J_1 / r_1 = \beta_1, 2J_1 / r_2 = \beta_2, -2L_1 / r_1 = \gamma_1 \text{ and } -2L_1 / r_2 = \gamma_2 \rightarrow \beta_1 = \beta_2 r_2 / r_1 \text{ and } \gamma_1 = \gamma_2 r_2 / r_1$$

$H_1 = -2G_1(r_1^2 + r_2^2)$, $F_1 = -G_1r_1^2r_2^2$, $D_1 = -2C_1(r_1^2 + r_2^2)$ and $B_1 = -C_1r_1^2r_2^2$. The coefficients of $\sin\theta$ and $\cos\theta$ for the stresses related to the direction $r$, are then written as follows:

For a general solution where there are no objective reasons to establish a dependency between the internal and external loading functions, this amounts to setting to zero these 2 constants: $L_1 = J_1 = 0$. There is, therefore, no term in $\cos\theta$ and $\sin\theta$ for $\sigma_{rr}$ and thus also none for the load functions $f_1(\theta)$ and $f_2(\theta)$. This is related to the fact that the

<figure><img src="image_3.png" /><figcaption>Fig. 2. Boundary conditions of the annular structure; example of radial loading, loading is proportional to length of the arrows whose colour varies from red (maximum value) to blue (minimum value). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>
solutions in cosθ and sinθ in the general solution produced by Michell (1899), correspond to degenerate cases. Furthermore, as highlighted by Barber (2010), it is essential for the boundary conditions to fulfill moment equilibrium requirements to ensure the static equilibrium of the hollow cylinder:

$$ 
\begin{aligned}
\int_{0}^{2\pi} f_1(\theta)r_1 \cos \theta d\theta &= \int_{0}^{2\pi} f_2(\theta)r_2 \cos \theta d\theta \quad \text{and} \\
&= \int_{0}^{2\pi} f_2(\theta)r_2 \sin \theta d\theta
\end{aligned}
\qquad (8) $$

These conditions are only true if $f_i(\theta)$ is of the form $p_i + b_i \cos n\theta + c_i \sin n\theta$ with $n \ge 2$ ($n = 1$ is here also a degenerate case for moment equilibrium). Without this load balance, the hollow cylinder will be subjected to a rigid displacement defined by $u_r = A \cos\theta + B \sin\theta$ and $u_\theta = -A \sin\theta + B \cos\theta + C r$ ($A, B$ and $C$ being arbitrary constants). This situation should be avoided because it can lead to numerical instabilities.

Finally, the Airy function for a complete hollow cylinder under solely internal and external (non-uniform) radial loads is:

$$ 
\begin{aligned}
\Phi(r, \theta) ={}& A_0 r^2 + C_0 \ln r + \sum_{n=2}^{\infty} \left( A_n r^n + \frac{B_n}{r^n} + C_n r^{2-n} + D_n r^{2+n} \right) \sin n\theta \\
& + \sum_{n=2}^{\infty} \left( E_n r^n + \frac{F_n}{r^n} + G_n r^{2-n} + H_n r^{2+n} \right) \cos n\theta
\end{aligned}
\qquad (9) $$

The radial stress derived from (9) is:

The boundary conditions for $r_1$ and $r_2$ imply loading functions compatible with $\sigma_{rr}(r_1)$ and $\sigma_{rr}(r_2)$, i.e. $f_1(\theta) = p_1 + \sum_{n=2}^{m} \beta_{1,n} \cos n\theta + \sum_{n=2}^{m} \gamma_{1,n} \sin n\theta$ and $f_2(\theta) = p_2 + \sum_{n=2}^{m} \beta_{2,n} \cos n\theta + \sum_{n=2}^{m} \gamma_{2,n} \sin n\theta$ where $p_1$ and $p_2$ correspond to the internal and external mean loads (axisymmetric case) while the coefficients $\beta_{i,n}$ and $\gamma_{i,n}$ represent the "deviatoric" (non-uniform) part of internal and external loads.

### 2.1.4. Determination of constants

For sake of simplicity, some parameters whose formulations differ only by a sign will be expressed more succinctly using the symbols $\pm$ and $\mp$. The upper sign indicates the top element in column vector notation. The constants of the problem can be determined by solving a system of $2+2 \times 4(m-1) = 8m-6$ equations ($n$ runs from 2 to $m$) in as many unknowns. The resolution can be simplified by examining the stress expressions:

$$ 
\begin{aligned}
\sigma_{rr} ={}& \frac{C_0}{r^2} + 2A_0 - \sum_{n=2}^{\infty} \sin n\theta \left( (n-1)nr^{n-2}A_n + n(n+1)\frac{B_n}{r^{2+n}} + (n+2)(n-1)\frac{C_n}{r^n} + (n-2)(n+1)r^nD_n \right) \\
& - \sum_{n=2}^{\infty} \cos n\theta \left( (n-1)nr^{n-2}E_n + n(n+1)\frac{F_n}{r^{2+n}} + (n+2)(n-1)\frac{G_n}{r^n} + (n-2)(n+1)r^nH_n \right)
\end{aligned}
\qquad (10) $$

$$ 
\begin{aligned}
\left( \begin{matrix} \sigma_{rr} \\ \sigma_{\theta\theta} \end{matrix} \right) ={}& 2A_0 \pm \frac{C_0}{r^2} \mp \sum_{n=2}^{\infty} \sin n\theta \left( (n-1)nr^{n-2}A_n + n(n+1)\frac{B_n}{r^{2+n}} + (n \pm 2)(n-1)\frac{C_n}{r^n} + (n \mp 2)(n+1)r^nD_n \right) \\
& \mp \sum_{n=2}^{\infty} \cos n\theta \left( (n-1)nr^{n-2}E_n + n(n+1)\frac{F_n}{r^{2+n}} + (n \pm 2)(n-1)\frac{G_n}{r^n} + (n \mp 2)(n+1)r^nH_n \right) \\
\sigma_{r\theta} ={}& \sum_{n=2}^{\infty} n \cos n\theta \left( -(n-1)r^{n-2}A_n + (n+1)\frac{B_n}{r^{2+n}} + (n-1)\frac{C_n}{r^n} - (n+1)r^nD_n \right) \\
& - \sum_{n=2}^{\infty} n \sin n\theta \left( -(n-1)r^{n-2}E_n + (n+1)\frac{F_n}{r^{2+n}} + (n-1)\frac{G_n}{r^n} - (n+1)r^nH_n \right)
\end{aligned}
\qquad (11) $$

Note that most terms will consist of a constant to be determined, a
multiplicative factor, and either cosnθ or sinnθ for n ≥ 2. To exploit the obvious symmetry of the terms multiplying cosnθ and sinnθ, a vertical bar has been introduced in Eqs (13) and (14) to separate these terms, which makes the formulation more compact. Also note that the system of equations can be decoupled.

The (axisymmetric) average load case is quickly solved (within the theoretical framework of continuum mechanics):

$$
\begin{cases}
C_0/r_1^2 + 2A_0 = p_1 & \text{thus } A_0 = \frac{1}{2} \frac{p_2 r_2^2 - p_1 r_1^2}{r_2^2 - r_1^2} \text{ and } C_0 = (p_2 - p_1) \frac{r_1^2 r_2^2}{r_1^2 - r_2^2} \\
C_0/r_2^2 + 2A_0 = p_2 &
\end{cases}
\tag{12}
$$

For the deviatoric part of the loading:

$$
\begin{pmatrix} E_n | A_n \\ F_n | B_n \\ G_n | C_n \\ H_n | D_n \end{pmatrix} = - \begin{pmatrix} n(n-1)r_2^{n-2} & \frac{n(n+1)}{r_2^{n+2}} & \frac{(n+2)(n-1)}{r_2^n} & (n-2)(n+1)r_2^n \\ n(n-1)r_1^{n-2} & \frac{n(n+1)}{r_1^{n+2}} & \frac{(n+2)(n-1)}{r_1^n} & (n-2)(n+1)r_1^n \\ -n(n-1)r_1^{n-2} & \frac{n(n+1)}{r_1^{n+2}} & \frac{n(n-1)}{r_1^n} & -n(n+1)r_1^n \\ -n(n-1)r_2^{n-2} & \frac{n(n+1)}{r_2^{n+2}} & \frac{n(n-1)}{r_2^n} & -n(n+1)r_2^n \end{pmatrix}^{-1} \times \begin{pmatrix} \beta_{2,n} | \gamma_{2,n} \\ \beta_{1,n} | \gamma_{1,n} \\ 0 \\ 0 \end{pmatrix} \tag{13}
$$

The solution of this subsystem of equations is:

2.1.5. Analytical solutions of stresses, strains and displacements in polar coordinates

Introducing the variable $\chi_i$ (also valid for layer $i$ of the multilayer case: §2.2) and the functions $h_i$ (of order $2n+2$) below enables a more concise representation of equations of stresses, strains and displacements:

$$
\begin{gather*}
\chi_i = \frac{1}{2r^{n+2}} \frac{1}{(r_{i+1}^{2n} - r_i^{2n})^2 - n^2 r_i^{2n-2} r_{i+1}^{2n-2} (r_{i+1}^2 - r_i^2)^2} \\
h_1(x, y) = (x^{2n} - r^{2n})(r^2 - y^2)
\end{gather*}
$$

$$
h_2(x, y) = n(x^{2n} + r^{2n})(r^2 - y^2)
$$

$$
h_3(x) = 2r^2(x^{2n} - r^{2n})
$$

$$
h_4(x, y) = n x^2 (y^{2n}(1 - n) + r^{2n}(n + 1))
$$

$$
h_5(x) = r^2(r^{2n}(1-n)(n+4\nu-2) + x^{2n}(n+1)(n-4\nu+2))
$$

$$
h_6(x) = r^2(r^{2n}(n-1)(n-4\nu+4) + x^{2n}(n+1)(n+4\nu-4))
$$

$$
h_7(x, y) = n x^2 (y^{2n}(n-1) + r^{2n}(n+1))
$$

Substituting equations (12) and (14) in equation (11) yields, after simplifications, the following expressions for the stresses:

$$
\begin{pmatrix} E_n | A_n \\ F_n | B_n \\ G_n | C_n \\ H_n | D_n \end{pmatrix} = - \frac{r_2^n}{2((r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2)} \begin{pmatrix} \frac{r_2^2(r_2^{2n} - r_1^{2n}) + n r_1^{2n} (r_2^2 - r_1^2)}{n-1} \\ r_1^{2n} r_2^2 \frac{r_2^{2n} - r_1^{2n} + n r_2^{2n-2} (r_2^2 - r_1^2)}{n+1} \\ r_1^{2n-2} r_2^2 \frac{r_2^{2n} - r_1^{2n} + n r_2^{2n} (r_2^2 - r_1^2)}{1-n} \\ \frac{r_2^{2n} - r_1^{2n} + n r_1^{2n-2} (r_2^2 - r_1^2)}{n+1} \end{pmatrix} (\beta_{2,n} | \gamma_{2,n}) \\
+ \frac{r_1^n}{2((r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2)} \begin{pmatrix} \frac{r_1^2(r_2^{2n} - r_1^{2n}) + n r_2^{2n} (r_2^2 - r_1^2)}{n-1} \\ r_2^{2n} r_1^2 \frac{r_2^{2n} - r_1^{2n} + n r_1^{2n-2} (r_2^2 - r_1^2)}{n+1} \\ r_1^{2n-2} r_2^2 \frac{r_2^{2n} - r_1^{2n} + n r_1^{2n} (r_2^2 - r_1^2)}{1-n} \\ \frac{r_2^{2n} - r_1^{2n} + n r_2^{2n-2} (r_2^2 - r_1^2)}{n+1} \end{pmatrix} (\beta_{1,n} | \gamma_{1,n}) \tag{14}
$$

$$
\left( \frac{\sigma_{\theta\theta}}{\sigma_{rr}} \right) = \frac{p_2 r_2^2 (r^2 \pm r_1^2) - p_1 r_1^2 (r^2 \pm r_2^2)}{r^2 (r_2^2 - r_1^2)} + \sum_{n=2}^{m \to \infty} \chi_1 \left[ \begin{array}{l} r_2^n (\beta_{2,n} \cos n\theta + \gamma_{2,n} \sin n\theta) \left( \frac{(r_2^{2n} - r_1^{2n})(-h_3(r_1) \pm h_2(r_1, r_2))}{n(r_2^2 - r_1^2)r_1^{2n-2}(-h_3(r_2) \pm h_2(r_2, r_1))} \right) \\ -r_1^n (\beta_{1,n} \cos n\theta + \gamma_{1,n} \sin n\theta) \left( \frac{(r_2^{2n} - r_1^{2n})(-h_3(r_2) \pm h_2(r_2, r_1))}{n(r_2^2 - r_1^2)r_2^{2n-2}(-h_3(r_1) \pm h_2(r_1, r_2))} \right) \end{array} \right] \quad (15)
$$
$$ \sigma_{r\theta} = \sum_{n=2}^{m \to \infty} n \chi_1 \left[ \begin{array}{l} r_2^n ((r_2^{2n} - r_1^{2n})h_1(r_1, r_2) + nr_1^{2n-2}(r_2^2 - r_1^2)h_1(r_2, r_1)) (\gamma_{2,n} \cos n\theta - \beta_{2,n} \sin n\theta) \\ -r_1^n ((r_2^{2n} - r_1^{2n})h_1(r_2, r_1) + nr_2^{2n-2}(r_2^2 - r_1^2)h_1(r_1, r_2)) (\gamma_{1,n} \cos n\theta - \beta_{1,n} \sin n\theta) \end{array} \right] \quad (16) $$

It can be verified (in equation (15)) that $\sigma_{rr} = p_1 + \sum_{n=2}^{m \to \infty} (\beta_{1,n} \cos n\theta + \gamma_{1,n} \sin n\theta) = f_1(\theta)$ if $r = r_1$ and that $\sigma_{rr} = p_2 + \sum_{n=2}^{m \to \infty} (\beta_{2,n} \cos n\theta + \gamma_{2,n} \sin n\theta) = f_2(\theta)$ if $r = r_2$. Similarly, it is easy to see (using (16)) that $\sigma_{r\theta} = 0$ if $r = r_1$ or $r_2$.

The expressions for the strains follow directly from Hooke's relation for plane strains: $\varepsilon_{rr} = (1 + \nu)((1 - \nu)\sigma_{rr} - \nu\sigma_{\theta\theta})/E$, $\varepsilon_{\theta\theta} = (1 + \nu)((1 - \nu)\sigma_{\theta\theta} - \nu\sigma_{rr})/E$ and $\varepsilon_{r\theta} = (1 + \nu)\sigma_{r\theta}/E$.

This yields the following relationships:

The most complex expression concerns the orthoradial displacement since it requires integrating the orthoradial strain (17), from which the

$$ \left( \frac{\varepsilon_{\theta\theta}}{\varepsilon_{rr}} \right) = \frac{1 + \nu}{E} \left\{ \sum_{n=2}^{m \to \infty} \chi_1 \left[ \begin{array}{l} \frac{(p_2 r_2^2 - p_1 r_1^2) r^2 (1 - 2\nu) \pm (p_2 - p_1) r_2^2 r_1^2}{r^2 (r_2^2 - r_1^2)} + \\ r_2^n (\beta_{2,n} \cos n\theta + \gamma_{2,n} \sin n\theta) \left( \begin{array}{l} (r_2^{2n} - r_1^{2n})((2\nu - 1)h_3(r_1) \pm h_2(r_1, r_2)) - \\ n(r_2^2 - r_1^2)r_1^{2n-2}((1 - 2\nu)h_3(r_2) \mp h_2(r_2, r_1)) \end{array} \right) \\ -r_1^n (\beta_{1,n} \cos n\theta + \gamma_{1,n} \sin n\theta) \left( \begin{array}{l} (r_2^{2n} - r_1^{2n})((2\nu - 1)h_3(r_2) \pm h_2(r_2, r_1)) - \\ n(r_2^2 - r_1^2)r_2^{2n-2}((1 - 2\nu)h_3(r_1) \mp h_2(r_1, r_2)) \end{array} \right) \end{array} \right\} \quad (17) $$

$$ \varepsilon_{r\theta} = \frac{1}{2G} \sum_{n=2}^{m \to \infty} n \chi_1 \left[ \begin{array}{l} r_2^n (\gamma_{2,n} \cos n\theta - \beta_{2,n} \sin n\theta) \left( \begin{array}{l} (r_2^{2n} - r_1^{2n})(r_2^2 - r_1^2) + \\ n r_1^{2n-2}(r_2^2 - r_1^2)(r_2^2 - r_1^2) \end{array} \right) \\ -r_1^n (\gamma_{1,n} \cos n\theta - \beta_{1,n} \sin n\theta) \left( \begin{array}{l} (r_2^{2n} - r_1^{2n})(r_2^2 - r_1^2) + \\ n r_2^{2n-2}(r_2^2 - r_1^2)(r_2^2 - r_1^2) \end{array} \right) \end{array} \right] \quad (18) $$

radial displacement (17) is subtracted: $u_\theta = \int (r\varepsilon_{\theta\theta} - u_r)d\theta$, therefore:

where the shear modulus G is equal to $E/2/(1+\nu)$.

The radial displacement is obtained by integrating the radial strain (17): $u_r = \int \varepsilon_{rr} dr$, and yields:

$$ u_r = r \frac{1 + \nu}{E} \left[ \sum_{n=2}^{m \to \infty} \chi_1 \left[ \begin{array}{l} \frac{(p_2 - p_1) \frac{r_2^2 r_1^2}{r^2(r_2^2 - r_1^2)} + \frac{p_2 r_2^2 - p_1 r_1^2}{r^2 - r_1^2}}{1 - 2\nu} + \\ r_2^n \frac{\beta_{2,n} \cos n\theta + \gamma_{2,n} \sin n\theta}{n^2 - 1} \left\{ (r_2^{2n} - r_1^{2n})(h_5(r_1) + h_4(r_2, r_1)) + n(r_2^2 - r_1^2)r_1^{2n-2}(h_5(r_2) + h_4(r_1, r_2)) \right\} \\ -r_1^n \frac{\beta_{1,n} \cos n\theta + \gamma_{1,n} \sin n\theta}{n^2 - 1} \left\{ (r_2^{2n} - r_1^{2n})(h_5(r_2) + h_4(r_1, r_2)) + n(r_2^2 - r_1^2)r_2^{2n-2}(h_5(r_1) + h_4(r_2, r_1)) \right\} \end{array} \right] \quad (19) $$
$$ u_{\theta} = r \frac{1 + \nu}{E} \sum_{n=2}^{m-\infty} \chi_1 \left[ \frac{r_2^{n} \beta_{2,n} \sin n\theta - \gamma_{2,n} \cos n\theta}{n^2 - 1} \left\{ (r_2^{2n} - r_1^{2n})(h_6(r_1) - h_7(r_2, r_1)) + n(r_2^2 - r_1^2)r_1^{2n-2}(h_6(r_2) - h_7(r_1, r_2)) \right\} \right. \\ \left. - r_1^n \frac{\beta_{1,n} \sin n\theta - \gamma_{1,n} \cos n\theta}{n^2 - 1} \left\{ (r_2^{2n} - r_1^{2n})(h_6(r_2) - h_7(r_1, r_2)) + n(r_2^2 - r_1^2)r_2^{2n-2}(h_6(r_1) - h_7(r_2, r_1)) \right\} \right] \quad (20) $$

<figure><img src="image_4.png" /><figcaption>Fig. 3. Distribution of $\sigma_{rr}$ in a single-layer hollow cylinder under a non-uniform tangential load. Arrows on the lower surface and upper surface of the ring visually demonstrate the non-uniform nature of the load, colors and lengths of these arrows correspond to the magnitude of $\sigma_{rr}$. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Comparison of analytical (plain lines) and numerical (diamonds) stress profiles (scale for $\sigma_{\theta\theta}$ on the left, scale for $\sigma_{rr}$ and $\sigma_{r\theta}$ on the right) in cross-section A ($\theta = 2.65$ rad).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Geometric representation of the continuity of the orthoradial strain at interface i between layers i and i+1.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Geometric representation of a trilayered hollow cylinder.</figcaption></figure>

Equations (15)–(20) are the elastic analytical solution of a complete hollow cylinder under internal and external loads of the form:

$$f_1(\theta) = p_1 + \sum_{n=2}^{m-\infty} (\beta_{1,n} \cos n\theta + \gamma_{1,n} \sin n\theta) \text{ and } f_2(\theta) = p_2 + \sum_{n=2}^{m-\infty} (\beta_{2,n} \cos n\theta + \gamma_{2,n} \sin n\theta)$$

with m the interpolation order of the load.

### 2.1.6. Analytical solution compared with a numerical model for a tangentially loaded hollow cylinder

Equations (15)–(20) and equations (15b) to (20b) from the appendix were verified by adjusting the biharmonic equation. In doing so, it was possible to express it in terms of the stresses rather than the stress function, enabling the concurrent verification of the three expressions for $\sigma_{rr}$, $\sigma_{\theta\theta}$ and $\sigma_{r\theta}$. It can also be expressed as a function of strains using Hooke's law with plane strain conditions.

Finally, knowing that $\Delta\Phi = \sigma_{rr} + \sigma_{\theta\theta} = \text{tr}(\sigma)$ i.e. $\Delta^2\Phi = \Delta(\sigma_{rr} + \sigma_{\theta\theta}) = 0$ and noting that $\Delta(\sigma_{rr} + \sigma_{\theta\theta}) = \frac{E}{(1+\nu)(1-2\nu)}\Delta(\epsilon_{rr} + \epsilon_{\theta\theta}) = 0$, we can also transform the Laplacian of the trace of the strain tensor into a relation involving only radial and tangential displacements:

$$\begin{aligned} \Delta^2\Phi &= \frac{\partial^2\sigma_{\theta\theta}}{\partial r^2} + \frac{2}{r}\frac{\partial\sigma_{\theta\theta}}{\partial r} + \frac{1}{r^2}\left(\frac{\partial^2\sigma_{rr}}{\partial\theta^2} + 2\frac{\partial^2\sigma_{\theta\theta}}{\partial\theta^2} + \sigma_{rr} - \sigma_{\theta\theta} + 3\frac{\partial\sigma_{r\theta}}{\partial\theta}\right) = 0 \\ \Delta^2\Phi &= (1-\nu)\frac{\partial^2\epsilon_{\theta\theta}}{\partial r^2} + \nu\frac{\partial^2\epsilon_{rr}}{\partial r^2} + \frac{2}{r}\left((1-\nu)\frac{\partial\epsilon_{\theta\theta}}{\partial r} + \nu\frac{\partial\epsilon_{rr}}{\partial r}\right) + \dots \\ &\quad \dots\frac{1}{r^2}\left((1+\nu)\frac{\partial^2\epsilon_{rr}}{\partial\theta^2} + (2-\nu)\frac{\partial^2\epsilon_{\theta\theta}}{\partial\theta^2} + (1-2\nu)(\epsilon_{rr} - \epsilon_{\theta\theta}) + 3(1-2\nu)\frac{\partial\epsilon_{r\theta}}{\partial\theta}\right) = 0 \end{aligned} \quad (21)$$

$$\begin{aligned} \Delta\text{tr}(\epsilon) = & \frac{\partial^3 u_r}{\partial r^3} + \frac{1}{r} \left( \frac{\partial^3 u_\theta}{\partial r^2 \partial \theta} + 2 \frac{\partial^2 u_r}{\partial r^2} \right) + \frac{1}{r^2} \left( \frac{\partial^3 u_r}{\partial \theta^2 \partial r} - \frac{\partial^2 u_\theta}{\partial r \partial \theta} - \frac{\partial u_r}{\partial r} \right) \\ & + \frac{1}{r^3} \left( \frac{\partial^3 u_\theta}{\partial \theta^3} + \frac{\partial^2 u_r}{\partial \theta^2} + \frac{\partial u_\theta}{\partial \theta} + u_r \right) = 0 \end{aligned}$$

After conducting these checks, the derived relationships for the single-layer hollow cylinder are compared with an FDM numerical model (Flac3D: Itasca, 2023) under both internal and external tangential loads (with m = 4). The relevant equations, specifically (15b) to (20b), are provided in the appendix. These solutions are valid within an elastic continuous medium and can be extended to the multilayer case, following the approach outlined in §2.2. Analytical solutions for radial loads are further compared with a numerical model in §2.2 featuring an illustrative example of a multilayered hollow cylinder.

Fig. 3 illustrates the distribution of $\sigma_{rr}$ in a single-layer ring under the influence of internal and external tangential loads $f_3(\theta)$ and $f_4(\theta)$, respectively. For all 28320 zones in the model, the analytical/numerical
relative deviation remains below 0.5 %. In Fig. 4, a remarkably good fit is observed between the numerical results and the analytical solution for the stress profiles (radial, orthoradial, and tangential) within section A (refer to Fig. 3).

factors correspond to the layer number and the interpolation order of the load: $\beta_{layer,order}$, $\gamma_{layer,order}$, $\kappa_{layer,order}$, $\lambda_{layer,order}$.

In the inner layer $i$, and at $r = r_{i+1}$ (extrados), the orthoradial strain is:

$$ \epsilon_{\theta\theta}^{i} = \frac{1 + \nu_{i}}{E_{i}} \left\{ \sum_{n=2}^{m \to \infty} \left[ \frac{(p_{i+1}r_{i+1}^{2} - p_{i}r_{i}^{2})(1 - 2\nu_{i}) + (p_{i+1} - p_{i})r_{i}^{2}}{r_{i+1}^{2} - r_{i}^{2}} + \frac{(r_{i+1}^{2n} - r_{i}^{2n})^{2}(1 - 2\nu_{i}) + n^{2}r_{i}^{2n-2}r_{i+1}^{2n-2}(r_{i+1}^{2} - r_{i}^{2})^{2}}{(r_{i+1}^{2n} - r_{i}^{2n})^{2} - n^{2}r_{i}^{2n-2}r_{i+1}^{2n-2}(r_{i+1}^{2} - r_{i}^{2})^{2}}(\beta_{i+1,n} \cos n\theta + \gamma_{i+1,n} \sin n\theta) \right. \\ \left. - \frac{2n(r_{i+1}^{2} - r_{i}^{2})(r_{i+1}^{2n} - r_{i}^{2n})r_{i}^{n}r_{i+1}^{n-2}(1 - \nu_{i})}{(r_{i+1}^{2n} - r_{i}^{2n})^{2} - n^{2}r_{i}^{2n-2}r_{i+1}^{2n-2}(r_{i+1}^{2} - r_{i}^{2})^{2}}(\beta_{i,n} \cos n\theta + \gamma_{i,n} \sin n\theta) \right] \right\} \quad (22) $$

## 2.2. Multilayer case

Here we use the solution of a single-layer hollow cylinder, which is transposed to a layer $i$ (by adapting the indices of the variables con-

In the upper layer $i+1$, and at $r = r_{i+1}$ (intrados) the orthoradial strain is:

$$ \epsilon_{\theta\theta}^{i+1} = \frac{1 + \nu_{i+1}}{E_{i+1}} \left\{ \sum_{n=2}^{m \to \infty} \left[ \frac{(p_{i+2}r_{i+2}^{2} - p_{i+1}r_{i+1}^{2})(1 - 2\nu_{i+1}) + (p_{i+2} - p_{i+1})r_{i+2}^{2} + \frac{r_{i+2}^{2} - r_{i+1}^{2}}{r_{i+2}^{2} - r_{i+1}^{2}}}{\frac{r_{i+2}^{2n} - r_{i+1}^{2n}}{(r_{i+2}^{2n} - r_{i+1}^{2n})^{2} - n^{2}r_{i+1}^{2n-2}r_{i+2}^{2n-2}(r_{i+2}^{2} - r_{i+1}^{2})^{2}}(\beta_{i+1,n} \cos n\theta + \gamma_{i+1,n} \sin n\theta) \right. \right. \\ \left. \left. - \frac{2n(r_{i+2}^{2} - r_{i+1}^{2})(r_{i+2}^{2n} - r_{i+1}^{2n})r_{i+1}^{n-2}r_{i+2}^{n-2}(1 - \nu_{i+1})}{(r_{i+2}^{2n} - r_{i+1}^{2n})^{2} - n^{2}r_{i+1}^{2n-2}r_{i+2}^{2n-2}(r_{i+2}^{2} - r_{i+1}^{2})^{2}}(\beta_{i+2,n} \cos n\theta + \gamma_{i+2,n} \sin n\theta) \right] \right\} \quad (23) $$

cerned). It is then sufficient to check the continuity of the radial and orthoradial displacement fields at the interface $i$ separating layer $i$ from layer $i+1$. Note that the continuity of $u_{\theta}$ is automatically checked with a uniform load since in this case $u_{\theta} = 0$.

The continuity condition for orthoradial displacement is based on the

The continuity in orthoradial strain at the interface no. $i$ between the layers $i$ and $i+1$ (Fig. 5) yields the transmission coefficients $\eta$ and $\lambda$ for the load on the extrados of layer $i+1$.

For the axisymmetric (mean) part of the load:

$$ p_{i+2} = \eta_{i+1,0} p_{i+1} + \lambda_{i+1,0} p_{i} \\ \eta_{i+1,0} = \frac{r_{i+2}^{2} + r_{i+1}^{2}(1 - 2\nu_{i+1})}{2r_{i+2}^{2}(1 - \nu_{i+1})} - \frac{G_{i}}{G_{i}} \frac{r_{i+1}^{2} - r_{i+2}^{2}}{r_{i+1}^{2} - r_{i+2}^{2}} \frac{r_{i}^{2} + r_{i+1}^{2}(1 - 2\nu_{i})}{2r_{i+2}^{2}(1 - \nu_{i+1})} \quad \lambda_{i+1,0} = \frac{1 - \nu_{i}^{2}}{1 - \nu_{i+1}^{2}} \frac{E_{i+1}}{E_{i}} \frac{r_{i+2}^{2} - r_{i+1}^{2}}{r_{i}^{2} - r_{i+1}^{2}} \frac{r_{i}^{2}}{r_{i+2}^{2}} \quad (24) $$

orthoradial strain for a layer $i$ between the inner and outer radii $r_i$ and $r_{i+1}$ (case of the hollow cylinder under internal load $f_i(\theta, p_b, \beta_{i,n}, \gamma_{i,n})$ and external load $f_{i+1}(\theta, p_{i+1}, \beta_{i+1,n}, \gamma_{i+1,n})$). Since $u_i$ is naturally continuous at the interface (by integrating solely over $r$), if follows that for $u_{\theta}$ to be continuous in $r = r_{i+1}$, $\epsilon_{\theta\theta}^{i} = \epsilon_{\theta\theta}^{i+1}$ must be verified.

For the deviatoric part of the load ($n \ge 2$):

* Layer $i+1$ ($r_{i+1} < r < r_{i+2}$): $u_{\theta}^{i+1} = r \int_{0}^{2\pi} \epsilon_{\theta\theta}^{i+1} d\theta - \int_{0}^{2\pi} u_{r}^{i+1} d\theta$

* Layer $i$ ($r_i < r < r_{i+1}$): $u_{\theta} = r \int_{0}^{2\pi} \epsilon_{\theta\theta}^{i} d\theta - \int_{0}^{2\pi} u_{r}^{i} d\theta$

The indices of the load function coefficients and the transmission
$$ 
\begin{gathered}
\begin{pmatrix} \beta_{i+2,n} \\ \gamma_{i+2,n} \end{pmatrix} = \eta_{i+1,n} \begin{pmatrix} \beta_{i+1,n} \\ \gamma_{i+1,n} \end{pmatrix} + \lambda_{i+1,n} \begin{pmatrix} \beta_{i,n} \\ \gamma_{i,n} \end{pmatrix} \\
\eta_{i+1,n} = \frac{(r_{i+2}^{2n} - r_{i+1}^{2n})^2 (1 - 2\nu_{i+1}) + n^2 r_{i+1}^{2n-2} r_{i+2}^{2n-2} (r_{i+2}^2 - r_{i+1}^2)^2}{2n(r_{i+2}^2 - r_{i+1}^2)(r_{i+2}^{2n} - r_{i+1}^{2n}) r_{i+2}^{n-2}(1 - \nu_{i+1})} \\
\frac{G_{i+1}}{G_i} \frac{(r_{i+2}^{2n} - r_{i+1}^{2n})^2 - n^2 r_{i+1}^{2n-2} r_{i+2}^{2n-2} (r_{i+2}^2 - r_{i+1}^2)^2}{2n(r_{i+2}^2 - r_{i+1}^2)(r_{i+2}^{2n} - r_{i+1}^{2n}) r_{i+2}^{n-2}(1 - \nu_{i+1})} \frac{(r_{i+1}^{2n} - r_i^{2n})^2 (1 - 2\nu_i) + n^2 r_i^{2n-2} r_{i+1}^{2n-2} (r_{i+1}^2 - r_i^2)^2}{(r_{i+1}^{2n} - r_i^{2n})^2 - n^2 r_i^{2n-2} r_{i+1}^{2n-2} (r_{i+1}^2 - r_i^2)^2} \\
\lambda_{i+1,n} = \frac{1 - \nu_i^2}{1 - \nu_{i+1}^2} \frac{E_{i+1}}{E_i} \frac{(r_{i+2}^{2n} - r_{i+1}^{2n})^2 - n^2 r_{i+1}^{2n-2} r_{i+2}^{2n-2} (r_{i+2}^2 - r_{i+1}^2)^2}{(r_{i+1}^{2n} - r_i^{2n})^2 - n^2 r_i^{2n-2} r_{i+1}^{2n-2} (r_{i+1}^2 - r_i^2)^2} \frac{(r_{i+1}^2 - r_i^2)(r_{i+1}^{2n} - r_i^{2n}) r_i^n}{(r_{i+2}^2 - r_{i+1}^2)(r_{i+2}^{2n} - r_{i+1}^{2n}) r_{i+2}^n}
\end{gathered}
\quad (25)
$$

<figure><img src="image_8.png" /><figcaption>Fig. 7. Representation of internal (data) and external (calculated from the load given in interface 1) applied to a trilayered hollow cylinder at the position specified by cross-section B, loading is proportional to length of the arrows whose colour varies from red (maximum value) to blue (minimum value). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

To find the load functions of the upper layers, we use these recurrence relations expressing the internal loading of layer i+2 (external loading of layer i+1) as a function of the internal loads of the 2 lower layers.

By recurrence, it is possible to express the internal and external loads of the multilayered hollow cylinder in terms of the known internal and external loads for a single layer. It thus becomes sufficient to calculate the transmission coefficients $\eta$ and $\lambda$ for the layers where the loads are unknown and to introduce them into the next iteration of the recurrence to calculate the distributions of $\sigma_{rr}$, $\sigma_{\theta\theta}$, $\sigma_{r\theta}$, $\varepsilon_{rr}$, $\varepsilon_{\theta\theta}$, $\varepsilon_{r\theta}$, $u_r$ and $u_\theta$ (up to one integration constant), the continuity of $\sigma_{rr}$ being automatically ensured at the interfaces by its definition (see equation (22)). The continuity of tangential quantities is trivial since $\sigma_{r\theta}$ and $\varepsilon_{r\theta}$ are zero at $r = r_i$ and $r = r_{i+1}$. Note that $u_r$ and $u_\theta$ are defined as integral expressions requiring the

calculation of integration constants (except for layer 1 where these constants are zero).

The relations (17), (19) and (20) can be transformed to comply with the plane stress assumption by replacing $E$ by $E(1-\nu^2)$ and $\nu$ by $\nu/(1+\nu)$.

These relations can also be made more universal by introducing the Kolosov constant which can take the following 2 values:

$$ \kappa_i = 3 - 4\nu_i \quad : \text{plane strains} $$

$$ \kappa_i = \frac{3 - \nu_i}{1 + \nu_i} \quad : \text{plane stresses} $$

The functions $h_5$ and $h_6$ defined in §2.1.5 are rewritten here to consider the Kolosov constant for a layer of index i ($\kappa_i$):

$$ h_5^\kappa(x) = r^2 (r^{2n}(1-n)(n+1-\kappa_i) + x^{2n}(n+1)(n-1+\kappa_i)) $$

$$ h_6^\kappa(x) = r^2 (r^{2n}(n-1)(n+1+\kappa_i) + x^{2n}(n+1)(n-1-\kappa_i)) $$
$$ \begin{pmatrix} \sigma_{\theta\theta}^i \\ \sigma_{rr}^i \end{pmatrix} = \frac{p_{i+1} r_{i+1}^2 (r^2 \pm r_i^2) - p_i r_i^2 (r^2 \pm r_{i+1}^2)}{r^2 (r_{i+1}^2 - r_i^2)} + \sum_{n=2}^{m \to \infty} \left[ \begin{array}{l} r_{i+1}^n \chi_i (\beta_{i+1,n} \cos n\theta + \gamma_{i+1,n} \sin n\theta) \left( \begin{array}{l} (r_{i+1}^{2n} - r_i^{2n})(-h_3(r_i) \pm h_2(r_i, r_{i+1})) + \\ n(r_{i+1}^2 - r_i^2)r_i^{2n-2}(-h_3(r_{i+1}) \pm h_2(r_{i+1}, r_i)) \end{array} \right) \\ -r_i^n \chi_i (\beta_{i,n} \cos n\theta + \gamma_{i,n} \sin n\theta) \left( \begin{array}{l} (r_{i+1}^{2n} - r_i^{2n})(-h_3(r_{i+1}) \pm h_2(r_{i+1}, r_i)) + \\ n(r_{i+1}^2 - r_i^2)r_{i+1}^{2n-2}(-h_3(r_i) \pm h_2(r_i, r_{i+1})) \end{array} \right) \end{array} \right] \quad (26) $$

$$ \sigma_{r\theta}^i = \sum_{n=2}^{m \to \infty} \left[ \begin{array}{l} n r_{i+1}^n \chi_i ((r_{i+1}^{2n} - r_i^{2n})h_1(r_i, r_{i+1}) + n r_i^{2n-2} (r_{i+1}^2 - r_i^2)h_1(r_{i+1}, r_i)) (\gamma_{i+1,n} \cos n\theta - \beta_{i+1,n} \sin n\theta) \\ -n r_i^n \chi_i ((r_{i+1}^{2n} - r_i^{2n})h_1(r_{i+1}, r_i) + n r_{i+1}^{2n-2} (r_{i+1}^2 - r_i^2)h_1(r_i, r_{i+1})) (\gamma_{i,n} \cos n\theta - \beta_{i,n} \sin n\theta) \end{array} \right] \quad (27) $$

$$ \begin{pmatrix} \varepsilon_{\theta\theta}^i \\ \varepsilon_{rr}^i \end{pmatrix} = \frac{1}{4G_i} \left\{ \sum_{n=2}^{m \to \infty} \left[ \begin{array}{l} \frac{(\kappa_i - 1)(p_{i+1}r_{i+1}^2 - p_i r_i^2)r^2 \pm 2(p_{i+1} - p_i)r_{i+1}^2 r_i^2}{r^2(r_{i+1}^2 - r_i^2)} + \\ r_{i+1}^n \chi_i (\beta_{i+1,n} \cos n\theta + \gamma_{i+1,n} \sin n\theta) \left( \begin{array}{l} (r_{i+1}^{2n} - r_i^{2n})((1 - \kappa_i)h_3(r_i) \pm 2h_2(r_i, r_{i+1})) - \\ n(r_{i+1}^2 - r_i^2)r_i^{2n-2}((\kappa_i - 1)h_3(r_{i+1}) \mp 2h_2(r_{i+1}, r_i)) \end{array} \right) \\ -r_i^n \chi_i (\beta_{i,n} \cos n\theta + \gamma_{i,n} \sin n\theta) \left( \begin{array}{l} (r_{i+1}^{2n} - r_i^{2n})((1 - \kappa_i)h_3(r_{i+1}) \pm 2h_2(r_{i+1}, r_i)) - \\ n(r_{i+1}^2 - r_i^2)r_{i+1}^{2n-2}((\kappa_i - 1)h_3(r_i) \mp 2h_2(r_i, r_{i+1})) \end{array} \right) \end{array} \right\} \quad (28) $$

$$ \bar{u}_r^i = \frac{r}{2G_i} \left[ \sum_{n=2}^{m \to \infty} \left[ \begin{array}{l} \frac{p_{i+1}r_{i+1}^2 - p_i r_i^2}{2(r_{i+1}^2 - r_i^2)} (\kappa_i - 1) + \frac{r_{i+1}^2 r_i^2}{r^2(r_{i+1}^2 - r_i^2)} (p_{i+1} - p_i) + \\ r_{i+1}^n \chi_i \frac{\beta_{i+1,n} \cos n\theta + \gamma_{i+1,n} \sin n\theta}{n^2 - 1} \left\{ (r_{i+1}^{2n} - r_i^{2n})(h_5^k(r_i) + h_4(r_{i+1}, r_i)) + n(r_{i+1}^2 - r_i^2)r_i^{2n-2}(h_5^k(r_{i+1}) + h_4(r_i, r_{i+1})) \right\} \\ -r_i^n \chi_i \frac{\beta_{i,n} \cos n\theta + \gamma_{i,n} \sin n\theta}{n^2 - 1} \left\{ (r_{i+1}^{2n} - r_i^{2n})(h_5^k(r_{i+1}) + h_4(r_i, r_{i+1})) + n(r_{i+1}^2 - r_i^2)r_{i+1}^{2n-2}(h_5^k(r_i) + h_4(r_{i+1}, r_i)) \right\} \end{array} \right] \quad (29) $$

$$ \bar{u}_{\theta}^i = \frac{r}{2G_i} \left[ \sum_{n=2}^{m \to \infty} \left[ \begin{array}{l} r_{i+1}^n \chi_i \frac{\beta_{i+1,n} \sin n\theta - \gamma_{i+1,n} \cos n\theta}{n^2 - 1} \left\{ (r_{i+1}^{2n} - r_i^{2n})(h_6^k(r_i) - h_7(r_{i+1}, r_i)) + n(r_{i+1}^2 - r_i^2)r_i^{2n-2}(h_6^k(r_{i+1}) - h_7(r_i, r_{i+1})) \right\} \\ -r_i^n \chi_i \frac{\beta_{i,n} \sin n\theta - \gamma_{i,n} \cos n\theta}{n^2 - 1} \left\{ (r_{i+1}^{2n} - r_i^{2n})(h_6^k(r_{i+1}) - h_7(r_i, r_{i+1})) + n(r_{i+1}^2 - r_i^2)r_{i+1}^{2n-2}(h_6^k(r_i) - h_7(r_{i+1}, r_i)) \right\} \end{array} \right] \quad (30) $$
For $i \ge 2$, accounting for integration constants becomes essential to obtain accurate values for radial and orthoradial displacements. This requirement leads to the formulation of the following recurrence relations:

$$ 
\begin{aligned}
u_r^i &= u_r^{i-1}(r=r_l) - \bar{u}_r^i(r=r_l) + \bar{u}_r^i \\
u_\theta^i &= u_\theta^{i-1}(r=r_l) - \bar{u}_\theta^i(r=r_l) + \bar{u}_\theta^i
\end{aligned} 
\quad (31) $$

Equations (23) and (24) can also be made more universal using the Kolosov constant. The expressions for the transmission coefficients $\lambda$ of the axisymmetric part can be made more compact by introducing the constant $\alpha^D$ defined by Dundurs (1969): $\alpha_l^D = \frac{G_{i+1}(1+k_i)-G_i(1+k_{i+1})}{G_{i+1}(1+k_i)+G_i(1+k_{i+1})}$

$$ 
\begin{aligned}
\eta_{i+1,0} &= \frac{2}{\kappa_{i+1} + 1} + \frac{G_{i+1} r_{i+2}^2 - r_{i+1}^2}{G_i r_{i+1}^2 - r_i^2} \frac{2r_i^2 + r_{i+1}^2(k_i - 1)}{r_{i+2}^2(k_{i+1} + 1)} + \frac{r_{i+1}^2(k_{i+1} - 1)}{r_{i+2}^2(k_{i+1} + 1)} \\
\lambda_{i+1,0} &= \frac{\alpha_i^D + 1}{\alpha_i^D - 1} \frac{r_{i+2}^2 - r_{i+1}^2}{r_{i+1}^2 - r_i^2} \frac{r_i^2}{r_{i+2}^2}
\end{aligned} 
\quad (32) $$

$$ 
\begin{aligned}
p_4 &= (\eta_{3,0}\eta_{2,0} + \lambda_{3,1})p_2 + \eta_{3,0}\lambda_{2,0}p_1 \\
\beta_{4,n} &= (\eta_{3,n}\eta_{2,n} + \lambda_{3,n})\beta_{2,n} + \eta_{3,n}\lambda_{2,n}\beta_{1,n} \\
\gamma_{4,n} &= (\eta_{3,n}\eta_{2,n} + \lambda_{3,n})\gamma_{2,n} + \eta_{3,n}\lambda_{2,n}\gamma_{1,n}
\end{aligned} 
\quad (34) $$

Displacements and stresses in the 3 layers of the hollow cylinder were calculated in an FDM model without initialization. The search for static equilibrium evidently demanded many iterations, yet no iteration would be necessary if these quantities were initialized using the analytical solutions. Fig. 7 illustrates the perfect fit of the numerical results with the analytical solution for the stress profiles (radial, orthoradial and tangential) of section B (Fig. 8). The deviations (analytical/FDM calculation, Figs. 9–11) in all displacements and stresses in the 35,840 meshes of the 3 layers of the hollow cylinder are less than $10^{-14}$ units (MPa or m), thus validating the analytical derivations.

In this example, where the intermediate layer is 400 times less rigid than the inner layer, the initialization of displacements and stresses using the analytical solution values helps avoid an additional computational cost of approximately 500,000 Flac3D iterations (FDM software). The manual accompanying this software acknowledges a decline in efficiency when solving problems involving large differences in rigidity (Itasca, 2023).

Figs. 9–11 below show the distribution of stresses, strains and dis-

$$ 
\begin{aligned}
\eta_{i+1,n} &= \frac{(r_{i+2}^{2n} - r_{i+1}^{2n})^2 (k_{i+1} - 1) + 2n^2 r_{i+1}^{2n-2} r_{i+2}^{2n-2} (r_{i+2}^2 - r_{i+1}^2)^2}{2n(r_{i+2}^2 - r_{i+1}^2)(r_{i+2}^{2n} - r_{i+1}^{2n})r_{i+1}^n r_{i+2}^{n-2}(k_{i+1} + 1)} - \dots \\
&\quad \frac{G_{i+1} (r_{i+2}^{2n} - r_{i+1}^{2n})^2 - n^2 r_{i+1}^{2n-2} r_{i+2}^{2n-2} (r_{i+2}^2 - r_{i+1}^2)^2}{2n(r_{i+2}^2 - r_{i+1}^2)(r_{i+2}^{2n} - r_{i+1}^{2n})r_{i+1}^n r_{i+2}^{n-2}(k_{i+1} + 1)} \frac{(r_{i+1}^{2n} - r_i^{2n})^2 (k_i - 1) + 2n^2 r_{i+1}^{2n-2} r_{i+2}^{2n-2} (r_{i+1}^2 - r_i^2)^2}{(r_{i+1}^{2n} - r_i^{2n})^2 - n^2 r_{i+1}^{2n-2} r_{i+2}^{2n-2} (r_{i+1}^2 - r_i^2)^2} \\
\lambda_{i+1,n} &= -\frac{\alpha_i^D + 1}{\alpha_i^D - 1} \frac{(r_{i+2}^{2n} - r_{i+1}^{2n})^2 - n^2 r_{i+1}^{2n-2} r_{i+2}^{2n-2} (r_{i+2}^2 - r_{i+1}^2)^2}{(r_{i+2}^{2n} - r_{i+1}^{2n})^2 - n^2 r_{i+1}^{2n-2} r_{i+2}^{2n-2} (r_{i+1}^2 - r_i^2)^2} \frac{(r_{i+1}^2 - r_i^2)(r_{i+1}^{2n} - r_i^{2n})r_i^n}{(r_{i+2}^2 - r_{i+1}^2)(r_{i+2}^{2n} - r_{i+1}^{2n})r_{i+1}^n}
\end{aligned} 
\quad (33) $$

placements in a trilayered hollow cylinder. Their analysis leads to several remarks.

1. The continuity of $\sigma_{rr}$, $\varepsilon_{\theta\theta}$, $u_r$, and $u_\theta$ is verified, as well as the discontinuous distribution of $\sigma_{\theta\theta}$ and $\varepsilon_{rr}$;

2. By expressing the orthoradial stress (plane strains) as a function of the radially continuous quantities $\sigma_{rr}$ and $\varepsilon_{\theta\theta}$, we obtain $\sigma_{\theta\theta} = E/(1 - \nu^2)\varepsilon_{\theta\theta} + \nu/(1 - \nu)\sigma_{rr}$. Considering the minimal variation in Poisson's ratio across different materials (layers), the primary factor leading to a discontinuity is Young's modulus. Specifically, at $r = r_2$ on the concrete side, $\sigma_{\theta\theta}$ is high because $E_1/(1 - \nu_1^2)\varepsilon_{\theta\theta}$ is significant, while on the side of the compressible material, $\sigma_{\theta\theta}$ is limited to $E_2/(1 - \nu_2^2)\varepsilon_{\theta\theta}$ because the term $\nu_2/(1 - \nu_2)\sigma_{rr} \approx \sigma_{rr}/(3 \text{ to } 4)$ is negligible;

1. in the absence of internal pressure ($p_1 = 0$) and assuming that the deviatoric part of the load can be neglected (amplitude lower than the average pressure), the distribution of $\sigma_{rr}$ in the compressible material (no. 2) is mainly governed by $\eta_{2,0} = \frac{r_2^2(1-2\nu_2)+r_3^2+\frac{E_2}{r_2^2}r_1^{2n-1}+\frac{E_3}{r_2^2}r_1^{2n-2}}{2r_3^2(1-\nu_2)}(r_1^2+r_2^2(1-2\nu_1))$, which can be less than 1 when $\frac{E_1}{E_2} > \frac{3\nu^2-4\nu\omega+\omega^2\nu^2}{\omega(2\omega-\omega)}$, with $\nu = \frac{r_3-r_2}{r_2-r_1}$, $\omega = r_3/r_2 - 1$, assuming that $\nu_1 = \nu_2 = 0.25$. With the geometry adopted in Fig. 4, $\eta_{2,0} < 1$ if $E_1 > 8.3 E_2$.

### 3. Comparison with a trilayered example and discussion

Relationships (25) through (29) were initially validated by verifying $\frac{\partial u_r}{\partial r} - \varepsilon_{rr} = 0$ and $\frac{1}{r} \left( \frac{\partial u_\theta}{\partial \theta} + u_r \right) - \varepsilon_{\theta\theta} = 0$, followed by utilisation of the adjusted biharmonic equation (21).

After these checks, these analytical solutions were compared with a numerical model (in Flac3D Itasca, 2023) of a trilayered hollow cylinder (Fig. 6) under an internal and external radial loading of layer 1 (internal layer).

For a trilayered hollow cylinder, adapting the indices of the relations (24) and (25) (or (32) and (33)), yields the values for $\eta_{3,0}$, $\lambda_{3,0}$, $\eta_{3,n}$ and $\lambda_{3,n}$ and consequently the extrados loading function of the third layer applied in the model (Fig. 7):

$$ p_4 = \eta_{3,0}p_3 + \lambda_{3,0}p_2, \beta_{4,n} = \eta_{3,n}\beta_{3,n} + \lambda_{3,n}\beta_{2,n} \text{ and } \gamma_{4,n} = \eta_{3,n}\gamma_{3,n} + \lambda_{3,n}\gamma_{2,n} $$

Since the internal loading in layer 3 is identical to the external loading in layer 2, the relationships derived depend exclusively on the loads applied to the first layer:
<figure><img src="image_9.png" /><figcaption>Fig. 8. Analytical-numerical comparison of stress profiles ($\sigma_{rr}$, $\sigma_{\theta\theta}$ and $\sigma_{r\theta}$) in cross-section B ($\theta$ = 2,2 rad).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Distribution of radial and orthoradial stresses in a trilayered hollow cylinder.</figcaption></figure>

This means that $\sigma_{rr}(r_3)<\sigma_{rr}(r_2)$ for a material 8.3 times less rigid than the underlying material.

In scenarios involving non-uniform loads, the derivation of analytical solutions for elastoplastic materials becomes impractical, if not impossible. Under uniform loading, it is often workable to express the principal stresses in polar coordinates (one often ascertains that $\sigma_1 = \sigma_{\theta\theta}$, $\sigma_3 = \sigma_{rr}$ and $\sigma_2 = \sigma_{zz}$). Failure criterion is typically expressed in principal stresses, often using $\sigma_1$ and $\sigma_3$, and sometimes using $\sigma_1$, $\sigma_2$ and $\sigma_3$. Under uniform loading, it is straightforward to rewrite the equilibrium equation (simplified due to axis symmetry) by incorporating the failure criterion (Li et al., 2022; Chen et al., 1999; Fritz, 1984).

Unfortunately, in the case of non-uniform loading, the simplifications allowed by axis symmetry are no longer valid. Importantly, the assumption that the minimum and maximum principal stresses are respectively radial and orthoradial stresses can no longer be maintained. The non-uniformity of the load implies that within the same layer, the orthoradial stress (as well as the radial stress) can be maximum in one location but minimum in another sector. Moreover, although these derivations are strictly valid in elastic conditions, it is recommended that the software compute an elastic equilibrium before introducing elastoplastic properties for the concerned materials. This precaution helps prevent potential errors arising from an overly approximate stress state during the equilibrium phase. Despite being applicable only in
<figure><img src="image_11.png" /><figcaption>Fig. 10. Distribution of radial and orthoradial strains in a trilayered hollow cylinder.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>a)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Distribution of radial and orthoradial displacements in a trilayered hollow cylinder.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>b)</figcaption></figure>

elastic scenarios, these analytical solutions provide a reliable means to expedite calculations—even in elastic cases where such solutions are no longer valid, as illustrated in Fig. 12 where a flaw is simulated by an empty space in the structure.

Furthermore, the relationships developed in this article can be extended to the thermoelastic case. This adaptation involves incorporating temperature and heat flux limits at the intrados and extrados of the layers and adjusting the expressions that connect strains to stresses when thermal coupling is considered.

It should be noted that non-uniform radial loading must be of the form $p + \sum_{n=2}^{m} (\beta \cos n\theta + \gamma \sin n\theta)$, this series expansion allows any form of loading to be represented by applying the Fourier method (Sadd, 2014). There may, however, be a limitation in the case of loading

calculated by minimization of $n_c$ measurements in the layers of the hollow cylinder (ideally, with equal angular and radial distributions). Indeed, for reasons of numerical stability, it is not advisable to have more unknowns than measurements, which means that if the number of sensors is small, the value of $m$ will also be low. This approximation occurs because the small number of sensors constrain the accuracy of the calculated loading, which would have been more accurate with a greater number of sensors.

## 4. Conclusion

A simple recursive analytical solution for predicting stress, strain and displacement fields in a multilayered elastic hollow cylinder under non-
<figure><img src="image_15.png" /><figcaption>Fig. 12. Distribution of orthoradial stress variations ($\Delta\sigma_{\theta\theta} = \sigma_{\theta\theta}^{\text{with defect}} - \sigma_{\theta\theta}^{\text{without defect}}$) in a trilayered hollow cylinder.</figcaption></figure>

uniform internal and external loads has been proposed. In particular, the proposed analytical solution can be used as an efficient, fast, and inex- pensive means to evaluate, for example, the loading at the interface of layers in instrumented annular structures (assessing stress, strain, and/ or displacement at the interface). The value of the parameter m (the interpolation order of the load in the analytical solution) will then be determined to best match the number of sensors ($n_c$) arranged in the hollow cylinder. It is numerically most efficient to set $2m-1 \approx n_c$ to pre- vent excessive constraints on the loading function. The results of this analytical solution are in perfect agreement with FDM calculations and can be easily extended to the thermoelastic case. Being able to calculate a non-uniform load at the interface of a layer of a cylindrical hollow

cylinder can help in assessing the contact between the roadway lining and the rock mass. The amplitude of the deviatoric part of the load can be seen as an over-contact or a defect of contact according to its sign. Finally, these solutions are not only appropriate for Civil Engineering but also for other disciplines, using composite tubular structures for instance.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A

If the boundary conditions (internal and/or external) of the hollow cylinder involve tangential stresses (non-zero functions $f_3(\theta)$ and $f_4(\theta)$), then equations (13)–(20) are modified to accommodate the coefficients $b_{1,n}, b_{2,n}, c_{1,n}$, and $c_{2,n}$ of the tangential load. These boundary conditions in $r_1$ and $r_2$ imply loading functions compatible with $\sigma_{r\theta}(r_1)$ and $\sigma_{r\theta}(r_2)$, i.e. $f_3(\theta) = \sum_{n=2}^{m} b_{1,n} \cos n\theta + \sum_{n=2}^{m} c_{1,n} \sin n\theta$ and $f_4(\theta) = \sum_{n=2}^{m} b_{2,n} \cos n\theta + \sum_{n=2}^{m} c_{2,n} \sin n\theta$, with the $b_{i,n}$ and $c_{i,n}$ coefficients representing internal and external non-uniform loads.

Note that the Airy function (9) is not modified by changing boundary conditions, especially for non-zero tangential stresses applied at $r = r_1$ and/or $r = r_2$. Indeed, the remaining constant $H_0$ in expression (4): $\sigma_{r\theta} = H_0/r^2 + \dots$, implies the dependence $\sigma_{r\theta}(r_1, \theta)/\sigma_{r\theta}(r_2, \theta) = r_2^2/r_1^2$, which is not an acceptable choice for a general analytic solution. Thus, tangential loading functions have no constant terms, consequently $H_0 = 0$, and these expressions become: $f_3(\theta) = \sum_{n=2}^{m} b_{1,n} \cos n\theta + \sum_{n=2}^{m} c_{1,n} \sin n\theta$ and $f_4(\theta) = \sum_{n=2}^{m} b_{2,n} \cos n\theta + \sum_{n=2}^{m} c_{2,n} \sin n\theta$. Moreover, even with non-zero tangential stresses applied to the hollow cylinder boundaries, the 8 constants occurring in the $\cos\theta$ and $\sin\theta$ terms of $\sigma_{rr}$ and $\sigma_{r\theta}$ must be zero to avoid dependencies between the radial and tangential loads: $B_1=C_1=D_1=F_1=G_1=H_1=J_1=L_1=0$.

The system to be solved becomes:

$$ \begin{pmatrix} E_n | A_n \\ F_n | B_n \\ G_n | C_n \\ H_n | D_n \end{pmatrix} = - \begin{pmatrix} n(n-1)r_2^{n-2} & \frac{n(n+1)}{r_2^{n+2}} & \frac{(n+2)(n-1)}{r_2^n} & (n-2)(n+1)r_2^n \\ n(n-1)r_1^{n-2} & \frac{n(n+1)}{r_1^{n+2}} & \frac{(n+2)(n-1)}{r_1^n} & (n-2)(n+1)r_1^n \\ -n(n-1)r_1^{n-2} & \frac{n(n+1)}{r_1^{n+2}} & \frac{n(n-1)}{r_1^n} & -n(n+1)r_1^n \\ -n(n-1)r_2^{n-2} & \frac{n(n+1)}{r_2^{n+2}} & \frac{n(n-1)}{r_2^n} & -n(n+1)r_2^n \end{pmatrix}^{-1} \times \begin{pmatrix} \beta_{2,n} | \gamma_{2,n} \\ \beta_{1,n} | \gamma_{1,n} \\ c_{1,n} | b_{1,n} \\ c_{2,n} | b_{2,n} \end{pmatrix} \quad (13b) $$

The solution of this system of equations is:
$$
\begin{pmatrix} E_n | A_n \\ F_n | B_n \\ G_n | C_n \\ H_n | D_n \end{pmatrix} = - \frac{r_2^n}{2((r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2)} \begin{pmatrix} \frac{r_2^2 (r_2^{2n} - r_1^{2n}) + n r_1^{2n} (r_2^2 - r_1^2)}{n-1} \\ r_1^{2n} r_2^2 \frac{r_2^{2n} - r_1^{2n} + n r_2^{2n-2} (r_2^2 - r_1^2)}{n+1} \\ r_1^{2n-2} r_1^2 \frac{r_2^{2n} - r_1^{2n} + n r_2^{2n} (r_2^2 - r_1^2)}{1-n} \\ - \frac{r_2^{2n} - r_1^{2n} + n r_1^{2n-2} (r_2^2 - r_1^2)}{n+1} \end{pmatrix} (\beta_{2,n} | \gamma_{2,n}) \\
+ \frac{r_1^n}{2((r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2)} \begin{pmatrix} \frac{r_1^2 (r_2^{2n} - r_1^{2n}) + n r_2^{2n} (r_2^2 - r_1^2)}{n-1} \\ r_1^{2n} r_2^2 \frac{r_2^{2n} - r_1^{2n} + n r_1^{2n-2} (r_2^2 - r_1^2)}{n+1} \\ r_2^{2n-2} r_2^2 \frac{r_2^{2n} - r_1^{2n} + n r_1^{2n} (r_2^2 - r_1^2)}{1-n} \\ - \frac{r_2^{2n} - r_1^{2n} + n r_2^{2n-2} (r_2^2 - r_1^2)}{n+1} \end{pmatrix} (\beta_{1,n} | \gamma_{1,n}) \\
+ \frac{r_2^n}{2((r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2)} \begin{pmatrix} \frac{(2-n) r_2^2 (r_2^{2n} - r_1^{2n}) + n^2 r_1^2 (r_2^2 - r_1^2)}{n(n-1)} \\ r_1^{2n} r_2^2 \frac{(n+2)(r_2^{2n} - r_1^{2n}) - n^2 r_2^{2n-2} (r_2^2 - r_1^2)}{n(n+1)} \\ r_1^{2n-2} r_1^2 \frac{r_2^{2n} - r_1^{2n}}{1-n} \\ - \frac{r_2^{2n} - r_1^{2n} - (n+2) r_1^{2n-2} (r_2^2 - r_1^2)}{n+1} \end{pmatrix} (b_{2,n} | c_{2,n}) \\
+ \frac{r_1^n}{2((r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2)} \begin{pmatrix} \frac{(n-2) r_1^2 (r_2^{2n} - r_1^{2n}) - n^2 r_2^{2n} (r_2^2 - r_1^2)}{n(n-1)} \\ r_2^{2n} r_1^2 \frac{-(n+2)(r_2^{2n} - r_1^{2n}) + n^2 r_1^{2n-2} (r_2^2 - r_1^2)}{n(n+1)} \\ r_2^{2n-2} r_2^2 \frac{r_2^{2n} - r_1^{2n}}{n-1} \\ - \frac{r_2^{2n} - r_1^{2n} - (n+2) n r_2^{2n-2} (r_2^2 - r_1^2)}{n+1} \end{pmatrix} (b_{1,n} | c_{1,n})
\qquad (14b)
$$

Equations (15b) to (20b) below correspond to analytical solutions with only tangential loading. They are added to equations (15)–(20) in the case of a load combining both radial and tangential components.

$$
\sigma_{rr} = \dots + \sum_{n=2}^{m-\infty} \left[ \frac{r_2^n}{2r^{n+2}} \frac{c_{2,n} \cos n\theta - b_{2,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2} \begin{pmatrix} (r_2^{2n} - r_1^{2n}) (r_2^2 - r_1^2) ((n-2)r^{2n} - (n+2)r_1^{2n}) + \\ (r_2^2 - r_1^2) r_1^{2n-2} (r_2^{2n} - r_1^{2n}) (n^2 r_1^2 - (n^2-4)r_1^2) \\ (r_2^{2n} - r_1^{2n}) (r_2^2 - r_1^2) ((n-2)r^{2n} - (n+2)r_2^{2n}) + \\ (r_2^2 - r_1^2) r_2^{2n-2} (r_2^{2n} - r_1^{2n}) (n^2 r_2^2 - (n^2-4)r_2^2) \end{pmatrix} \right] \quad (15b)
$$
$$
\sigma_{ij} = \cdots +
$$

$$
\sum_{n=2}^{m \to \infty} \left[
\begin{aligned}
& \frac{r_2^n}{2r^{n+2}} \frac{c_{2,n} \cos n\theta - b_{2,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n} r_2^{2n-2} (r_2^2 - r_1^2)^2} \left( (r_2^{2n} - r_1^{2n}) (2(r_2^{2n} + r_1^{2n}) (r_2^2 + r^2) - n(r_2^{2n} - r_1^{2n}) (r_2^2 - r^2)) + (r_2^2 - r_1^2) r_1^{2n-2} \left( 4r^2((1-n)r_2^{2n} - (n+1)r_2^{2n}) + n^2(r_2^2 - r_1^2)(r_2^{2n} - r^{2n}) \right) \right) \\
& - \frac{r_1^n}{2r^{n+2}} \frac{c_{1,n} \cos n\theta - b_{1,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n} r_2^{2n-2} (r_2^2 - r_1^2)^2} \left( (r_2^{2n} - r_1^{2n}) (2(r_2^{2n} + r_2^{2n}) (r_1^2 + r^2) - n(r_2^{2n} - r^{2n}) (r_2^2 - r_1^2)) + (r_2^2 - r_1^2) r_2^{2n-2} \left( 4r^2((1-n)r_1^{2n} - (n+1)r_1^{2n}) + n^2(r_2^2 - r^2)(r_2^{2n} - r_1^{2n}) \right) \right)
\end{aligned}
\right]
$$

$$
\sigma_{r\theta} = \dots + \sum_{n=2}^{m \to \infty} \left[ \frac{r_2^n}{2r^{n+2}} \frac{b_{2,n} \cos n\theta + c_{2,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n} r_2^{2n-2} (r_2^2 - r_1^2)^2} \left( (r_2^{2n} - r_1^{2n}) (2(r_2^{2n} - r_1^{2n})r_2^2 - n(r_2^{2n} + r_1^{2n})(r_2^2 - r^2)) + nr_1^{2n-2}(r_2^2 - r_1^2)(2r^2(r_2^{2n} - r^{2n}) - n(r_2^2 - r_1^2)(r_2^{2n} + r^{2n})) \right) \right. \\
\left. + \frac{r_1^n}{2r^{n+2}} \frac{b_{1,n} \cos n\theta + c_{1,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n} r_2^{2n-2} (r_2^2 - r_1^2)^2} \left( (r_2^{2n} - r_1^{2n}) (2(r_2^{2n} - r^{2n})r_1^2 - n(r_2^{2n} + r_2^{2n})(r_2^2 - r_1^2)) + nr_2^{2n-2}(r_2^2 - r_1^2)(2r^2(r_2^{2n} - r_1^{2n}) - n(r_2^2 - r_1^2)(r_1^{2n} + r^{2n})) \right) \right] \quad (16b)
$$

It can be verified that (equation (16b)) that $\sigma_{r\theta} = \sum_{n=2}^{m\to\infty} (b_{1,n} \cos n\theta + c_{1,n} \sin n\theta) = f_3(\theta)$ if $r=r_1$ and that
$\sigma_{rr} = \sum_{n=2}^{m\to\infty} (b_{2,n} \cos n\theta + c_{2,n} \sin n\theta) = f_4(\theta)$ if $r=r_2$.

$$
\epsilon_{rr} = \cdots +
$$

$$
\frac{1+\nu}{E} \sum_{n=2}^{m \to \infty} \left[ \begin{array}{l} \frac{r_2^n}{2r^{n+2}} \frac{c_{2,n} \cos n\theta - b_{2,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n} r_2^{2n-2} (r_1^2 - r_2^2)^2} \left( (r_2^{2n} - r_1^{2n}) \right. \\ \left. - (r_2^2 - r_1^2) r_1^{2n-2} \left( r^2 \left( \begin{array}{l} (n+2)(n-2+4\nu)r_1^{2n} \\ -(n-2)(n+2-4\nu)r_2^{2n} \end{array} \right) + n^2 r_1^2 (r_2^{2n} - r^{2n}) \end{array} \right) \\[1em] \frac{r_1^n}{2r^{n+2}} \frac{c_{1,n} \cos n\theta - b_{1,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n} r_2^{2n-2} (r_1^2 - r_2^2)^2} \left( (r_2^{2n} - r_1^{2n}) \right. \\ \left. - (r_2^2 - r_1^2) r_2^{2n-2} \left( r^2 \left( \begin{array}{l} (n+2)(n-2+4\nu)r_2^{2n} \\ -(n-2)(n+2-4\nu)r_1^{2n} \end{array} \right) + n^2 r_2^2 (r_1^{2n} - r^{2n}) \end{array} \right) \end{array} \right] \quad (17b)
$$
$$ \epsilon_{\theta\theta} = \dots + $$

$$ \frac{1 + \nu}{E} \sum_{n=2}^{m \to \infty} \left\{ \begin{array}{l} \frac{r_2^n}{2r^{n+2}} \frac{c_{2,n} \cos n\theta - b_{2,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_1^2 - r_2^2)^2} \left( \begin{array}{l} (r_2^{2n} - r_1^{2n}) \left( r_2^2((n+2)r_1^{2n} - (n-2)r_2^{2n}) + r_1^2((n+2-4\nu)r_2^{2n} - (n-2+4\nu)r_1^{2n}) \right) \\ + (r_2^2 - r_1^2)r_1^{2n-2} \begin{array}{l} r^2 \left( \begin{array}{l} (n-2)(n-2+4\nu)r_2^{2n} - \\ (n+2)(n+2-4\nu)r_1^{2n} \end{array} \right) \\ - n^2 r_1^2 (r_2^{2n} - r_1^{2n}) \end{array} \right) \\ \\ - \frac{r_1^n}{2r^{n+2}} \frac{c_{1,n} \cos n\theta - b_{1,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_1^2 - r_2^2)^2} \left( \begin{array}{l} (r_2^{2n} - r_1^{2n}) \left( r_1^2((n+2)r_2^{2n} - (n-2)r_1^{2n}) + r_2^2((n+2-4\nu)r_1^{2n} - (n-2+4\nu)r_2^{2n}) \right) \\ + (r_2^2 - r_1^2)r_2^{2n-2} \begin{array}{l} r^2 \left( \begin{array}{l} (n-2)(n-2+4\nu)r_1^{2n} - \\ (n+2)(n+2-4\nu)r_2^{2n} \end{array} \right) \\ - n^2 r_2^2 (r_1^{2n} - r_2^{2n}) \end{array} \right) \end{array} \right\} $$

$$ \epsilon_{r\theta} = \dots + $$

$$ \frac{1}{4G} \sum_{n=2}^{m \to \infty} \left[ \frac{r_2^n}{r^{n+2}} \frac{b_{2,n} \cos n\theta + c_{2,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_1^2 - r_2^2)^2} \left( \begin{array}{l} (r_2^{2n} - r_1^{2n}) (2(r_2^{2n} - r_1^{2n})r_2^2 - n(r_2^{2n} + r_1^{2n})(r_2^2 - r_1^2)) + \\ nr_1^{2n-2} (r_2^2 - r_1^2) (2r_2^2(r_2^{2n} - r_1^{2n}) - n(r_2^2 - r_1^2)(r_2^{2n} + r_1^{2n})) \end{array} \right) + \\ + \frac{r_1^n}{r^{n+2}} \frac{b_{1,n} \cos n\theta + c_{1,n} \sin n\theta}{(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_1^2 - r_2^2)^2} \left( \begin{array}{l} (r_2^{2n} - r_1^{2n}) (2(r_2^{2n} - r_2^{2n})r_1^2 - n(r_2^{2n} + r_2^{2n})(r_1^2 - r_2^2)) + \\ nr_2^{2n-2} (r_1^2 - r_2^2) (2r_2^2(r_1^{2n} - r_2^{2n}) - n(r_2^2 - r_2^2)(r_1^{2n} + r_2^{2n})) \end{array} \right) \right] \quad (18b) $$

$$ u_r = \dots + $$

$$ \frac{1 + \nu}{2E} \sum_{n=2}^{m \to \infty} \left[ \frac{r_2^n}{r^{n+1}} \frac{c_{2,n} \cos n\theta - b_{2,n} \sin n\theta}{(n^2 - 1) ((r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2)} \left\{ \begin{array}{l} (r_2^{2n-2} (n^2 r_1^2 - (n-2)(n+1)r_2^2) - (n+2)r_1^{2n})(n-1)r_2^2 r_1^{2n} + \\ (r_1^{2n-2} ((n+2)r_2^2 - (n+1)r_1^2) - r_2^{2n})(n-1)(n-2+4\nu)r_1^{2n+2} + \\ + (r_1^{2n+2} + r_2^{2n} ((n-2)r_2^2 - (n-1)r_1^2))(n+1)(n+2-4\nu)r_1^{2n-2} \\ + ((n-2)r_2^{2n+2} + r_1^{2n} (n^2 r_1^2 - (n+2)(n-1)r_2^2))(n+1)r_2^{2n} \end{array} \right\} \\ \\ + \frac{r_1^n}{r^{n+1}} \frac{c_{1,n} \cos n\theta - b_{1,n} \sin n\theta}{(n^2 - 1) ((r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2)} \left\{ \begin{array}{l} (r_1^{2n-2} (n^2 r_2^2 - (n-2)(n+1)r_1^2) - (n+2)r_2^{2n})(n-1)r_1^2 r_2^{2n} + \\ (r_2^{2n-2} ((n+2)r_1^2 - (n+1)r_2^2) - r_1^{2n})(n-1)(n-2+4\nu)r_2^{2n+2} + \\ + (r_2^{2n+2} + r_1^{2n} ((n-2)r_1^2 - (n-1)r_2^2))(n+1)(n+2-4\nu)r_2^{2n-2} \\ + ((n-2)r_1^{2n+2} + r_2^{2n} (n^2 r_2^2 - (n+2)(n-1)r_1^2))(n+1)r_1^{2n} \end{array} \right\} \end{array} \right] $$

(19b)
$$u_{\theta} = \cdots +$$

$$ \frac{1 + \nu}{2E} \sum_{n=2}^{m \to \infty} \left[ \frac{r_2^n}{r^{n+1}} \frac{c_{2,n} \cos n\theta - b_{2,n} \sin n\theta}{(n^2 - 1)(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2} + \frac{r_1^n}{r^{n+1}} \frac{c_{1,n} \cos n\theta - b_{1,n} \sin n\theta}{(n^2 - 1)(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2} \left\{ \begin{array}{l} r_2^{2n-2} (n^2 r_1^2 - (n-2)(n+1)r_2^2) - (n+2)r_1^{2n} (n-1)r_1^{2n} r_2^2 + \\ r_1^{2n} ((n+2)(n-1)r_2^2 - n^2 r_1^2) - (n-2)r_2^{2n+2} (n+1)r_2^{2n} + \\ r_2^{2n} + r_1^{2n-2} ((n+1)r_1^2 - (n+2)r_2^2) (n-1)(n+4-4\nu)r^{2+2n} + \\ r_1^{2n+2} + r_2^{2n} ((n-2)r_2^2 - (n-1)r_1^2) (n+1)(n-4+4\nu)r^2 r_1^{2n-2} \end{array} \right\} + \frac{r_1^n}{r^{n+1}} \frac{c_{1,n} \cos n\theta - b_{1,n} \sin n\theta}{(n^2 - 1)(r_2^{2n} - r_1^{2n})^2 - n^2 r_1^{2n-2} r_2^{2n-2} (r_2^2 - r_1^2)^2} \left\{ \begin{array}{l} r_2^{2n-2} (n^2 r_2^2 - (n-2)(n+1)r_1^2) - (n+2)r_2^{2n} (n-1)r_2^{2n} r_1^2 + \\ r_2^{2n} ((n+2)(n-1)r_1^2 - n^2 r_2^2) - (n-2)r_1^{2n+2} (n+1)r_1^{2n} + \\ r_1^{2n} + r_2^{2n-2} ((n+1)r_2^2 - (n+2)r_1^2) (n-1)(n+4-4\nu)r^{2+2n} + \\ r_2^{2n+2} + r_1^{2n} ((n-2)r_1^2 - (n-1)r_2^2) (n+1)(n-4+4\nu)r^2 r_2^{2n-2} \end{array} \right\} \right] \quad (20b) $$

## Appendix B. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.euromechsol.2025.105926.

## Data availability

No data was used for the research described in the article.

## References

Atangana Nkene, E.R., Mambou Ngueyep, L.L., Ndop, J., Djiokeng, E.S., Ndjaka, J.M.B., 2018. Displacements, strains, and stresses investigations in an inhomogeneous rotating hollow cylinder made of functionally graded materials under an axisymmetric radial loading. World J. Mech. 8 (3), 59–72. https://doi.org/10.4236/wjm.2018.83005.

Barber, J.R., 2010. Elasticity, third ed. Springer, Dordrecht Heidelberg London New York. https://doi.org/10.1007/978-90-481-3809-8.

Chen, X., Tan, C.P., Haberfield, M., 1999. Solutions for the deformations and stability of elastoplastic hollow cylinders subjected to boundary pressures. Int. J. Numer. Anal. Methods GeoMech. 23, 779–800. https://doi.org/10.1002/(SICI)1096-9853 (199907)23:8%3C779::AID-NAG10%3E3.0.CO, 2-0.

Dundurs, J., 1969. Discussion on 'Edge bonded dissimilar orthogonal elastic wedges under normal and shear loading'. J. Appl. Mech. 36 (3), 650–652. https://doi.org/10.1115/1.3564739.

Fritz, P., 1984. An analytical solution for axisymmetric tunnel problems in elastoviscoplastic media. Int. J. Numer. Anal. Methods GeoMech. 8, 325–342.

Hefny, A.M., Lo, K.Y., 1999. Analytical solutions for stresses and displacements around tunnels driven in cross-anisotropic rocks. Int. J. Numer. Anal. Methods GeoMech. 23 (2), 161–177. https://doi.org/10.1002/(SICI)1096-9853(199902)23:2<161::AID-NAG963>3.0.CO;2-B.

Hongjun, X., Zhifei, S., Taotao, Z., 2006. Elastic analyses of heterogeneous hollow cylinders. Mech. Res. Commun. 33 (5), 681–691. https://doi.org/10.1016/j.mechrescom.2006.01.005.

Itasca, 2023. FLAC3D (Fast Lagrangian Analysis of Continua in Three Dimensions) Version 9.0, User's Guide and Command Reference. Itasca Consulting Group Inc., Minneapolis, USA.

Jaeger, J.C., Cook, N.G.W., 1968. Fundamentals of Rock Mechanics. Methuen & Co, London. Published by.

Khechai, A., Belarbi, M.O., Bouaziz, A., 2023. A general analytical solution of stresses around circular holes in functionally graded plates under various in-plane loading conditions. Acta Mech. 234, 671–691. https://doi.org/10.1007/s00707-022-03413-1.

Kirsch, E.G., 1898. Die Theorie der Elastizität und die Bedürfnisse der Festigkeitslehre. Z. Des. Vereines Dtsch. Ingenieure 42, 797–807 (in german).

Lamé, G., 1852. Leçons Sur la Théorie Mathématique De L'Élasticité Des Corps Solides. impr.-libr. Bachelier, Paris, p. 335 (in French).

Li, G., Li, N., Bai, Y., 2022. A new elastic-plastic analytical solution of circular tunnel under non-axisymmetric conditions. Sci. Rep. 12, 4367. https://doi.org/10.1038/s41598-022-08353-3.

Masoudian, M.S., Hashemi, M.A., 2016. Analytical solution of a circular opening in an axisymmetric elastic-brittle-plastic swelling rock. J. Nat. Gas Sci. Eng. 35, 483–496. https://doi.org/10.1016/j.jngse.2016.08.076.

Michell, J.H., 1899. On the direct determination of stress in an elastic solid, with application to the theory of plates. Proc. Lond. Math. Soc. 31, 100–124. https://doi.org/10.1112/plms/s1-31.1.100.

Obert, L., Duvall, W.I., 1967. Rock Mechanics and the Design of Structure in Rock. J. Wiley & Sons Inc., New York London Sydney.

Park, K.H., Kim, Y.J., 2006. Analytical solution for a circular opening in an elastic-brittle-plastic rock. Int. J. Rock Mech. Min. Sci. 43, 616–622. https://doi.org/10.1016/j.ijrmms.2005.11.004.

Qi, G., Liu, C., Feng, K., Ma, L., Schroder, K.U., 2023. Analytical solutions of stress distribution within a hollow cylinder under contact interactions. Int. J. Mech. Sci. 239, 107897. https://doi.org/10.1016/j.jimcsci.2022.107897.

Qiu, J., Zhou, M., 2016. Analytical solution for interference fit for multi-layer thick-walled cylinders and the application in crankshaft bearing design. Appl. Sci. 6 (6), 167. https://doi.org/10.3390/app6060167.

Sadd, M.H., 2014. Elasticity. Theory, Applications, and Numerics, third ed. Elsevier. https://doi.org/10.1016/C2012-0-06981-5.

Savin, G.N., Guz, A.N., 1964. On the stress state near curvilinear holes in shells. Izv. SSSR, Mekhanika i Mashinostroenie (6).

Savin, G.N., 1951. Stress Concentrations near Holes. GTITL, Moscow [in Russian].

Sharan, S.K., 2005. Exact and approximate solutions for displacements around circular openings in elastic-brittle-plastic Hoek-Brown rock. Int. J. Rock Mech. Min. Sci. 42, 542–549. https://doi.org/10.1016/j.ijrmms.2005.03.019.

Stasiak, M., 2019. Compression uniaxiale d'un matériau granulaire fragile très compressible - une étude par éléments discrets 3D. PhD of Université Grenoble Alpes (in french). https://www.theses.fr/2019GREAI041.

Timoshenko, S., Goodier, J.N., 1951. Theory of Elasticity, second ed. McGraw-Hill.

Tokovyy, Y.V., Ma, C.C., 2018. Analysis of 2D non-axisymmetric elasticity and thermoelasticity problems for radially inhomogeneous hollow cylinders. J. Eng. Math. 61, 171–184. https://doi.org/10.1007/s10665-007-9154-6.

Van Sambeek, L.L., 1986. Creep of Rock Salt Under Inhomogeneous Stress Conditions. Ph.D. Thesis. Colorado School of Mines.

Yeo, W.H., Purbolaksono, J., Aliabadi, M.H., Ramesh, S., Liew, H.L., 2017. Exact solution for stresses/displacements in a multilayered hollow cylinder under thermo-mechanical loading. Int. J. Pres. Ves. Pip. 151, 45–53. https://doi.org/10.1016/j.ijpvp.2017.01.003.