Closed-Form Analysis of Stress and Deformation in Functionally Graded Multi-Layer Hyperelastic Cylinders under Internal Pressure

1. Introduction

Exact, large-deformation predictions of pressure-driven stresses in incompressible cylindrical vessels are essential for safe and efficient design under service loads. Finite-strain, calibration-based formulations outperform linearized surrogates, and analytic benchmarks remain vital because volumetric locking and stabilization methods can bias high-pressure stress predictions [1,2,3]. Recent exact and semi-analytical vessel studies show peak stresses are highly sensitive to constitutive model and boundary conditions, driving the need for parameter-transparent predictions [4,5]. Biomedical applications heighten these requirements, while studies of graded elastomers demonstrate that smooth through-thickness property variation reshapes pressure-driven stress profiles, motivating exact treatments of gradation effects [6,7].

Despite established large-deformation vessel studies, a clear gap remains: an exact hyperelastic formulation that predicts pressure–stress fields in continuously graded, multi-layer hyperelastic cylinders subjected to internal pressure. Analytical (ANL) and semi-analytical multilayer treatments provide closed-form relations for layered vessels, but their piecewise-constant property descriptions can obscure how truly continuous gradation redistributes pressure-driven stresses and masks gradation-specific effects [8]. Semi-analytical finite-length and three-dimensional treatments demonstrate the importance of strictly enforcing boundary conditions and incompressibility to control peak stresses, yet most do not couple those enforcement details to a continuous through-thickness gradation [9,10,11]. Independent analyses of FG hyperelastic cylinders show that smooth radial variation in constitutive parameters reshapes radial and hoop stress profiles, motivating continuous gradation laws rather than stepwise layering when accurate, parameter-transparent predictions are required [12]. Benchmark ring and tube studies further show that abrupt stiffness jumps introduce artificial stress discontinuities that complicate comparisons with graded limits, underscoring the need for an exact formulation that preserves incompressibility while eliminating interface artifacts when passing to FG descriptions [13]. Together, these works set a clear target: an exact, closed-form FG finite-strain theory for internally pressurized cylinders that yields parameter-transparent pressure–stress predictions and recovers multilayer and homogeneous results as piecewise-constant and zero-gradient limits, respectively [8,14].

Assuming axisymmetry with plane strain fixes the kinematics and allows exact enforcement of incompressibility. Recent mixed and stabilized formulations demonstrate that substituting an artificially large bulk modulus remains prone to volumetric locking and pressure bias, whereas mixed displacement–pressure schemes preserve reliable stress predictions at finite strain [15,16]. Benchmarks for large-strain hyperelastic elements further show that predicted stress peaks and through-thickness gradients depend critically on how incompressibility and boundary conditions are enforced as well as on the chosen strain-energy form (e.g., Neo-Hookean versus Mooney–Rivlin), reinforcing the need to specify the constitutive class explicitly and to audit enforcement details [17]. In numerical practice, contact algorithms or artificial prescriptions (for example, virtual thermal strains used to mimic assembly loads) can conflate loading kinematics with solution procedures and thereby obscure parameter transparency at interfaces [18,19]. In this study we instead prescribe boundary tractions and geometric constraints directly, enforce exact incompressibility, and require continuity of displacement and stress across layers. Functional gradation is introduced as a smooth radial variation of the stored-energy parameters; recent analyses of graded tubes under coupled inflation–torsion illustrate how continuous laws reshape hoop and axial responses but do not address the combined effect of continuous gradation and explicit boundary-traction prescriptions [20,21,22,23]. These considerations define our modeling choices: axisymmetric, plane-strain, incompressible hyperelasticity with direct enforcement of incompressibility and a smooth FG law tied to the stored-energy parameters [15,16,17,18,20].

Radial equilibrium together with exact incompressibility reduces the boundary-value problem to a single scalar condition for an integration constant that is fixed by the prescribed radial tractions. Once that constant is found, the radial mapping and closed-form expressions for radial displacement and for all stress components follow algebraically. Multilayer analytical treatments make this scalar structure explicit for piecewise-constant parameter distributions [24]. Analytical solutions for FG cylinders (including studies of inflation and coupled loading) preserve the same scalar organization after replacing piecewise constants with smooth radial functions, although many such works treat different loading modes or simplifying assumptions and do not emphasize a traction-based scalar closure for graded materials [25] Other graded-cylinder analyses likewise exploit the scalar reduction but focus on alternate loadings or numerical strategies rather than on an explicit, closed-form traction closure [26,27]. By contrast, finite-element approximations generally require mixed or stabilized formulations to control volumetric locking and spurious pressure oscillations; their accuracy depends sensitively on how near-incompressibility and boundary data are enforced, and an explicit scalar closure typically emerges only through post-processing rather than directly from the discrete formulation [28,29]. This scalar framework yields an analytically tractable basis for quantifying how continuous gradation modifies stress distributions and for producing design-relevant, parameter-transparent maps [24,25,26,28,29].

We present an exact, closed-form solution for the large-deformation radial mapping, displacements, and full stress field of a thick-walled, incompressible, FG multi-layer hyperelastic cylinder under internal pressure. The formulation enforces incompressibility exactly, admits general smooth through-thickness gradation laws (here illustrated with a Mooney–Rivlin model and an exponential–logarithmic gradation), and reduces the boundary-value problem to a single scalar closure whose solution yields algebraic expressions for radial, hoop and axial stresses. The theory recovers piecewise-constant multilayer and homogeneous limits and is applied to representative bi-layer and tri-layer PVC [30] cases with two gradation variants. Analytical predictions are benchmarked against mixed/hybrid finite-element solutions. Parametric studies quantify how grading strength and layer contrasts govern stress smoothing and peak magnitudes, producing directly usable maps for design of graded elastomer pressure vessels. The closed-form base state also enables straightforward extensions (incremental stability, non-axisymmetric or finite-length numerical studies, and thermo-mechanical coupling).

2. Materials and Methods

2.1. Materials

2.1.1. Material Model and Properties

In continuum mechanics, hyperelastic materials are characterized by an energy density function. In this study, the Mooney-Rivlin model was employed, which consists of two material parameters, denoted as C10 and C01.

To incorporate FG distribution of these material constants, two cases were studied: one with increasing and the other with decreasing material properties as the radius increases. Multiplication factors (MFs) of 1.2 and 0.8 were used to scale the properties to their maximum and minimum values from the base values. The main values, which maximization and minimization were performed on, were chosen by the properties belong to PVC, which are reported in Table 1.

In each case study, regardless of the number of layers and whether the material properties increase or decrease with the radius, two values were defined for each parameter within a layer. For instance, C_{10, max} and C_{10, min} refer to the maximum and minimum values of C₁₀ within the layer. Similarly, C_{01, max} and C_{01, min} refer to the maximum and minimum values of C₀₁ within the same layer. The C₁₀ and C₀₁ properties of PVC were utilized for the first layer.

Additionally, in this study, two cases of bi-layer and multi-layer cylinders were analyzed. In both cases, the material properties of the inner layer were higher than those of the outer layers. This is because the inner layer is in direct contact with the location where the pressure is applied, requiring to be stronger. In the bi-layer cylinder, the properties of the second layer were 50% lower than those of the first layer. The material properties of the bi-layer cylinder with two MFs of 0.8 and 1.2 are provided in Table 2 and Table 3, respectively.

In the multi-layer cylinder, the properties of the second and third layers were 30% and 60% lower than those of the first layer, respectively. The material properties of the multi-layer cylinder, with two MFs of 0.8 and 1.2 are provided in Table 4 and Table 5, respectively.

2.2. Methods

2.2.1. Analytical Solution

Figure 1 illustrates the schematic of a thick, long, multi-layer cylinder under radial pressure. Due to the infinite length of the cylinder, plain strain assumption is made to simplify the analytical approach, as the longitudinal variation of stress and displacement fields can be neglected. A distinct material is assigned to each layer of the cylinder, and the thickness can vary from one layer to another. However, in this study, equal thickness is assumed to all layers. Given the theory of hyperelasticity, which is a branch of nonlinear mechanics, the proposed analytical solution must be capable of capturing large deformation and inflation caused by the applied radial pressure.

The inner and outer radial pressures are denoted as $p_{in}$ and $p_{out}$, respectively. Each layer consists of two radii: the inner and outer radii, denoted as $R_{in}and{R}_{out}$, respectively. The radii $R_1,R_2,R_3,\mathrm{a}\mathrm{n}\mathrm{d}{R}_{n+1}$ correspond to the inner radius of the first, second, third, and the last layer, respectively.

The general form of the deformation gradient in the cylindrical coordinate system [33] is as follows:

$$F=\left[\begin{array}{ccc}{\displaystyle \frac{\partial r}{\partial R}}& {\displaystyle \frac{1}{R}}{\displaystyle \frac{\partial r}{\partial \mathsf{\Theta }}}& {\displaystyle \frac{\partial r}{\partial Z}}\\ r{\displaystyle \frac{\partial \theta }{\partial R}}& {\displaystyle \frac{r}{R}}{\displaystyle \frac{\partial \theta }{\partial \mathsf{\Theta }}}& r{\displaystyle \frac{\partial \theta }{\partial \mathrm{Z}}}\\ {\displaystyle \frac{\partial z}{\partial R}}& {\displaystyle \frac{1}{R}}{\displaystyle \frac{\partial z}{\partial \mathsf{\Theta }}}& {\displaystyle \frac{\partial z}{\partial \mathrm{Z}}}\end{array}\right].$$

(1)

In this problem, all fields depend only on the reference radial coordinate, so the deformation is represented by the monotone radial mapping,

$$r=r\left(R\right), \theta =\mathsf{\Theta },\mathrm{z}=\mathrm{Z}.$$

(2)

where $r$ and $R$ are the current and reference radii. Using this assumption, the deformation gradient tensor is defined as:

$$F=\left[\begin{array}{ccc}{\displaystyle \frac{\partial r}{\partial R}}& 0& 0\\ 0& {\displaystyle \frac{r}{R}}& 0\\ 0& 0& 1\end{array}\right].$$

(3)

The left Cauchy-Green deformation tensor is given as follows:

$$b=F{F}^T=\left[\begin{array}{ccc}{\left({\displaystyle \frac{\partial r}{\partial R}}\right)}^2& 0& 0\\ 0& {\left({\displaystyle \frac{r}{R}}\right)}^2& 0\\ 0& 0& 1\end{array}\right].$$

(4)

The first three principal invariants of the left Cauchy-Green deformation tensors are defined as $I_1,I_2,and{I}_3$ as follows:

$$I_1=tr\left(b\right)={\left({\displaystyle \frac{\partial r}{\partial R}}\right)}^2+{\left({\displaystyle \frac{r}{R}}\right)}^2+1,$$

(5)

$${ I}_2={\displaystyle \frac{1}{2}}\left[{\left(trb\right)}^2-tr\left(b^2\right)\right]=\left[{\left({\displaystyle \frac{\partial r}{\partial R}}\right)}^2{\left({\displaystyle \frac{r}{R}}\right)}^2+{\left({\displaystyle \frac{\partial r}{\partial R}}\right)}^2+{\left({\displaystyle \frac{r}{R}}\right)}^2\right],$$

(6)

$$I_3=\left|b\right|={\left|F\right|}^2=J^2=1.$$

(7)

Exact incompressibility closes the kinematics through the radial identity as follows:

$$J=1\to \left({\displaystyle \frac{\partial r}{\partial R}}\right)\left({\displaystyle \frac{r}{R}}\right)=1\to r{\displaystyle \frac{\partial r}{\partial R}}=R.$$

(8)

The displacement field is defined as:

$$u=\mathrm{r}-\mathrm{R}.$$

(9)

Integration of the incompressibility relation yields an explicit primitive for mapping,

$$r^2\left(R\right)=R^2+c.$$

(10)

where $c$ is the integration constant.

To impose the continuity condition on the radial deformation of the successive layers, the following steps were taken:

$$at R=R_{n+1}\to u_{r, n}=u_{r,n+1}.$$

(11)

$u$ was substituted from Eq. (9) into Eq. (11):

$$\sqrt{R_{n+1}^2+c_n}-R_{n+1}=\sqrt{R_{n+1}^2+c_{n+1}}-R_{n+1}.$$

(12)

which concludes:

$$c_n=c_{n+1} ;n=1, 2, \dots , n-1.$$

(13)

This implies that $c_n$ does not change from one layer to another.

For an incompressible hyperelastic cylinder, the Mooney-Rivlin model [34] is defined as:

$$\psi =C_{10} \left(I_1-3\right)+ C_{01} \left(I_2-3\right)$$

(14)

where $C_{10} \mathrm{a}\mathrm{n}\mathrm{d} C_{01}$ are the material constants.

To obtain an expression for the different components of Cauchy stress, the second Piola stress was first defined as:

$$S=2{\displaystyle \frac{\partial \psi }{\partial C}}=2\left[{\displaystyle \frac{\partial \psi }{\partial I_1}}{\displaystyle \frac{\partial I_1}{\partial C}}+{\displaystyle \frac{\partial \psi }{\partial I_2}}{\displaystyle \frac{\partial I_2}{\partial C}}+{\displaystyle \frac{\partial \psi }{\partial I_3}}{\displaystyle \frac{\partial I_3}{\partial C}}\right]$$

(15)

The derivatives of the principal invariants with respect to the right Cauchy-Green deformation tensor can be calculated as:

$${\displaystyle \frac{\partial I_1}{\partial C}}=\mathrm{I}; {\displaystyle \frac{\partial I_2}{\partial C}}=I_1\mathrm{I}-C; {\displaystyle \frac{\partial I_3}{\partial C}}=I_3{C}^{-1}.$$

(16)

By substituting Eq. (16) into Eq. (15), the expression for the second Piola stress tensor [35] is as follows:

$$S=2\left[{\displaystyle \frac{\partial \psi }{\partial I_1}}I+{\displaystyle \frac{\partial \psi }{\partial I_2}}\left(I_1\mathrm{I}-C\right)+{\displaystyle \frac{\partial \psi }{\partial I_3}}\left(I_3{C}^{-1}\right)\right].$$

(17)

The Cauchy stress formula can be written using its relation with the second Piola stress as follows:

$$\begin{gathered} \sigma ={\displaystyle \frac{1}{J}}FS{F}^T={\displaystyle \frac{2}{J}}F\left[{\displaystyle \frac{\partial \psi }{\partial I_1}}I+{\displaystyle \frac{\partial \psi }{\partial I_2}}\left(I_1\mathrm{I}-C\right)+{\displaystyle \frac{\partial \psi }{\partial I_3}}\left(I_3{C}^{-1}\right)\right]F^T \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{2}{\sqrt{I_3}}}\left[{\displaystyle \frac{\partial \psi }{\partial I_1}}b+{\displaystyle \frac{\partial \psi }{\partial I_2}}\left(I_1\mathrm{B}-{\mathrm{F}\left(F^T\mathrm{F}\right)F}^T\right)+{\displaystyle \frac{\partial \psi }{\partial I_3}}\left(I_{3}F\left(F^{-1}{F}^{-T}\right)F^T\right)\right] \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\displaystyle \frac{2}{\sqrt{I_3}}}\left[{\displaystyle \frac{\partial \psi }{\partial I_1}}b+{\displaystyle \frac{\partial \psi }{\partial I_2}}\left(I_1\mathrm{b}-b^2\right)+{\displaystyle \frac{\partial \psi }{\partial I_3}}\left(I_3\mathrm{I}\right)\right]. \end{gathered}$$

(18)

The index notation for Eq. (18) is:

$${\sigma }_{ij}={\displaystyle \frac{2}{\sqrt{I_3}}}\left[\left({\displaystyle \frac{\partial \psi }{\partial I_1}}+I_1{\displaystyle \frac{\partial \psi }{\partial I_2}}\right){\mathrm{b}}_{ij}-{\displaystyle \frac{\partial \psi }{\partial I_2}}\sum _{m=1}^3{\mathrm{b}}_{im}{\mathrm{b}}_{mj}\right]+2\sqrt{I_3}{\displaystyle \frac{\partial \psi }{\partial I_3}}{\delta }_{ij}.$$

(19)

An incompressible material is only deformed when subjected to isochoric stress, which is denoted as ${\overline{\sigma }}_{ij}.$Incompressibility constraint does not allow the material to change its volume, so volumetric stress, which is denoted as $-p$ plays no role in material deformation.

The total stress, which is the summation of isochoric and volumetric stress, is defined as follows:

$${\sigma }_{ij}={\overline{\sigma }}_{ij}+\left(-p\right){\delta }_{ij}.$$

(20)

where ${\delta }_{ij}$ is the kronecker delta, defined as:

$${\delta }_{ij}=\left\{\begin{array}{c}0\\ 1\end{array}\begin{array}{c} if i\ne j\\ if i=j\end{array}\right.$$

(21)

The principal stress is defined as:

$${\sigma }_i\equiv {\sigma }_{ii}={\overline{\sigma }}_i-p$$

(22)

Eq. (18) is simplified by using the Eq. (22), $I_3=1,\mathrm{a}\mathrm{n}\mathrm{d}{\displaystyle \frac{\partial \psi }{\partial I_3}}=0$ as:

$${\sigma }_i=2\left[\left({\displaystyle \frac{\partial \psi }{\partial I_1}}+I_1{\displaystyle \frac{\partial \psi }{\partial I_2}}\right){\mathrm{b}}_{ii}-{\displaystyle \frac{\partial \psi }{\partial I_2}}\sum _{m=1}^3{\mathrm{b}}_{im}{\mathrm{b}}_{mi}\right]-p$$

(23)

The radial, tangential, and axial stress components, denoted as ${\sigma }_r,{\sigma }_{\theta },and{\sigma }_z$, respectively, are defined using Eq. (23):

$$\begin{gathered} {\sigma }_r=2{\displaystyle \frac{\partial \psi }{\partial I_1}}{\mathrm{b}}_{11}+2{\displaystyle \frac{\partial \psi }{\partial I_2}}\left(I_1{\mathrm{b}}_{11}-{\mathrm{b}}_{11}^2\right)-p \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=2C_{10}\left({\displaystyle \frac{R^2}{r^2}}\right)+2C_{01}\left(\left(\left({\displaystyle \frac{R^2}{r^2}}\right)+\left({\displaystyle \frac{r^2}{R^2}}\right)+1\right)\left({\displaystyle \frac{R^2}{r^2}}\right)-\left({\displaystyle \frac{R^4}{r^4}}\right)\right)-p \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\stackrel{1+\left({\displaystyle \frac{R^2}{r^2}}\right)=-{\displaystyle \frac{r^2}{R^2}}}{\to } {\sigma }_r=2C_{10}{\left({\displaystyle \frac{R}{r}}\right)}^2-2C_{01}{\left({\displaystyle \frac{r}{R}}\right)}^2-p \end{gathered}$$

(24)

$${\sigma }_{\theta }=2{\displaystyle \frac{\partial \psi }{\partial I_1}}{\mathrm{b}}_{22}+2{\displaystyle \frac{\partial \psi }{\partial I_2}}\left(I_1{\mathrm{b}}_{22}-{\mathrm{b}}_{22}^2\right)-p=2C_{10}\left({\displaystyle \frac{r^2}{R^2}}\right)+2C_{01}\left(\left(\left({\displaystyle \frac{R^2}{r^2}}\right)+\left({\displaystyle \frac{r^2}{R^2}}\right)+1\right)\left({\displaystyle \frac{r^2}{R^2}}\right)-\left({\displaystyle \frac{r^4}{R^4}}\right)\right)-p\stackrel{1+\left({\displaystyle \frac{r^2}{R^2}}\right)=-{\displaystyle \frac{R^2}{r^2}}}{\to }{ \sigma }_{\theta }=2C_{10}{\left({\displaystyle \frac{r}{R}}\right)}^2-2C_{01}{\left({\displaystyle \frac{R}{r}}\right)}^2-p$$

(25)

$${\sigma }_z=2{\displaystyle \frac{\partial \psi }{\partial I_1}}{\mathrm{b}}_{33}+2{\displaystyle \frac{\partial \psi }{\partial I_2}}\left(I_1{\mathrm{b}}_{33}-{\mathrm{b}}_{33}^2\right)-p=2{\displaystyle \frac{\partial \psi }{\partial I_1}}-2{\displaystyle \frac{\partial \psi }{\partial I_2}}-p={\sigma }_r-2{\displaystyle \frac{\partial \psi }{\partial I_1}}\left[{\left({\displaystyle \frac{R}{r}}\right)}^2-1\right]+2{\displaystyle \frac{\partial \psi }{\partial I_2}}\left[{\left({\displaystyle \frac{r}{R}}\right)}^2-1\right]-p={\sigma }_r+2\left(C_{01}{\displaystyle \frac{c}{R^2}}+C_{10}{\displaystyle \frac{c}{R^2+c}}\right)-p$$

(26)

In polar coordinates, in the absence of body forces, the radial stress equilibrium is defined as follows:

$${\displaystyle \frac{\mathrm{d}{\mathsf{\sigma }}_r}{\mathrm{d}\mathrm{r}}}+{\displaystyle \frac{{\mathsf{\sigma }}_r-{\mathsf{\sigma }}_{\theta }}{r}}=0$$

(27)

Now, the expression related to ${\mathsf{\sigma }}_r$ and ${\mathsf{\sigma }}_{\theta }$ can be substituted from Eq. (24) and Eq. (25), respectively into Eq. (27):

$${\displaystyle \frac{\mathrm{d}{\mathsf{\sigma }}_r}{\mathrm{d}\mathrm{r}}}=-{\displaystyle \frac{2}{r}}\left(C_{10}+C_{01}\right)\left[{\left({\displaystyle \frac{R}{r}}\right)}^2-{\left({\displaystyle \frac{r}{R}}\right)}^2\right]$$

(28)

Using Eq. (8), Eq. (28) becomes:

$${\displaystyle \frac{\mathrm{d}{\mathsf{\sigma }}_r}{\mathrm{d}\mathrm{R}}}=-{\displaystyle \frac{2R}{r^2}}\left(C_{10}+C_{01}\right)\left[{\left({\displaystyle \frac{R}{r}}\right)}^2-{\left({\displaystyle \frac{r}{R}}\right)}^2\right]$$

(29)

Then, integration is applied to both sides of Eq. (29):

$${\mathsf{\sigma }}_r={\int }_{R_n}^R{\displaystyle \frac{d{\mathsf{\sigma }}_r}{dR}} dR={\int }_{R_n}^R\left[-{\displaystyle \frac{2R}{r^2}}\left(C_{10}+C_{01}\right)\left[{\left({\displaystyle \frac{R}{r}}\right)}^2-{\left({\displaystyle \frac{r}{R}}\right)}^2\right]\right] dR$$

(30)

By using Eq. (10), Eq. (30) is rewritten as follows:

$${\mathsf{\sigma }}_r={\int }_{R_n}^R{\displaystyle \frac{d{\mathsf{\sigma }}_r}{dR}} dR={\int }_{R_n}^R\left[-2\left(C_{10}+C_{01}\right)\left[{\displaystyle \frac{R^3}{{\left(R^2+c_n\right)}^2}}-{\displaystyle \frac{1}{R}}\right]\right] dR$$

(31)

Eq. (31) is suitable for problems where material properties are assumed to be homogenous. To generalize the solution for heterogeneous material distributions, a FG approach is employed for the Mooney-Rivlin material parameters, C10 and C01.

To account for the FG distribution, an exponential-logarithmic distribution form of C10 and C01 is defined in Eq. (32) and Eq. (33), respectively:

$$C_{10}\left(R\right)={C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right).$$

(32)

$$C_{01}\left(R\right)={C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right).$$

(33)

The homogeneous limit is recovered when ${\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}=1\mathrm{a}\mathrm{n}\mathrm{d}{\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}=1.$The expressions for $C_{10}\left(R\right)\mathrm{a}\mathrm{n}\mathrm{d}{C}_{01}\left(R\right)$ are required to be substituted into Eq. (31) to employ the integration process. To facilitate integration, the Taylor expansion is used for the exponential part of Eq. (32) and Eq. (33), expanded up to four sentences. The expansion is used as follows:

$$C_{10}\left(R\right)={C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(1+\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^2+{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^3\right).$$

(34)

$$C_{01}\left(R\right)={C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(1+\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^2+{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^3\right).$$

(35)

Now, $C_{10}\left(R\right)\mathrm{a}\mathrm{n}\mathrm{d}{C}_{01}\left(R\right)$ are substituted into Eq. (31):

$${\mathsf{\sigma }}_r={\int }_{R_n}^R\left[-2\left(\left({C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(1+\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^2+{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^3\right)\right)+\left({C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(1+\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^2+{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^3\right)\right)\right)\left[{\displaystyle \frac{R^3}{{\left(R^2+c_n\right)}^2}}-{\displaystyle \frac{1}{R}}\right]\right] dR-p_n.$$

(36)

where $p_n$ denotes the pressure applied on the radius of $R_n$. The above integral requires a long-hand solution. The final closed-form solution is reported in Appendix A.1.

By considering ${R=R}_{n+1}$, Eq. (36) becomes as follows:

$$p_n-p_{n+1} ={\int }_{R_n}^R\left[-2\left(\left({C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(1+\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^2+{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^3\right)\right)+\left({C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(1+\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^2+{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{R-R_{out}}{R_{out}-R_{in}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}}}\right)\right)}^3\right)\right)\right)\left[{\displaystyle \frac{R^3}{{\left(R^2+c_n\right)}^2}}-{\displaystyle \frac{1}{R}}\right]\right] dR.$$

(37)

where $p_{n+1}$ denotes the pressure applied on the radius of $R_{n+1}$, which is equal to ${\mathsf{\sigma }}_{n+1}$.

To obtain $c_n$, which is the same for all layers within the cylinder, Eq. (37) must be solved successively, as shown in Eq. (38):

$$p_{in}-p_{out}=\sum _{m=1}^{m=n}{p}_m-p_{m+1}=2\sum _{m=1}^n\left({\int }_{R_n}^R\left[-2\left(\left({C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{m}}\left(1+\left({\displaystyle \frac{R-R_{out,\mathrm{m}}}{R_{out,\mathrm{m}}-R_{in,\mathrm{m}}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{m}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{m}}}}\right)\right)+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{R-R_{out,\mathrm{m}}}{R_{out,\mathrm{m}}-R_{in,\mathrm{m}}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{m}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{m}}}}\right)\right)}^2+{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{R-R_{out,\mathrm{m}}}{R_{out,\mathrm{m}}-R_{in,\mathrm{m}}}}\ln \left({\displaystyle \frac{{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{m}}}{{C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{m}}}}\right)\right)}^3\right)\right)+\left({C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(1+\left({\displaystyle \frac{R-R_{out,\mathrm{m}}}{R_{out,\mathrm{m}}-R_{in,\mathrm{m}}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{m}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{m}}}}\right)\right)+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{R-R_{out,\mathrm{m}}}{R_{out,\mathrm{m}}-R_{in,\mathrm{m}}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{m}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{m}}}}\right)\right)}^2+{\displaystyle \frac{1}{6}}{\left({\displaystyle \frac{R-R_{out,\mathrm{m}}}{R_{out,\mathrm{m}}-R_{in,\mathrm{m}}}}\ln \left({\displaystyle \frac{{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{m}}}{{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{m}}}}\right)\right)}^3\right)\right)\right)\left[{\displaystyle \frac{R^3}{{\left(R^2+c_n\right)}^2}}-{\displaystyle \frac{1}{R}}\right]\right] dR\right)$$

(38)

where ${C_{10}}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{m}},{C_{10}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{m}},{C_{01}}_{\mathrm{m}\mathrm{i}\mathrm{n},\mathrm{m}}\mathrm{a}\mathrm{n}\mathrm{d}{C_{01}}_{\mathrm{m}\mathrm{a}\mathrm{x},\mathrm{m}}$ are the minimum and maximum values of $C_{10}$ and $C_{01}$in the m-th layer, respectively. $R_{in,\mathrm{m}}\mathrm{a}\mathrm{n}\mathrm{d}{R}_{out,\mathrm{m}}$, denote the inner and outer radii of the m-th layer. Also, to implement the traction free of outermost surface, $p_{out}$ is considered zero.

By starting from the outer radius and solving Eq. (37) successively, the pressure applied on each boundary radius is obtained.

After solving Eq. (36), the tangential and axial stresses are calculated in Eq. (39) and Eq. (40), respectively:

$${\sigma }_{\theta }={\sigma }_r+2{\left(C_{10}+C_{01}\right)}_k\left({\displaystyle \frac{R^2+c_n}{R^2}}-{\displaystyle \frac{R^2}{R^2+c_n}}\right).$$

(39)

$${\sigma }_z={\sigma }_r+2\left[{\left(C_{10}\right)}_k{\displaystyle \frac{c_n}{R^2}}+{\left(C_{01}\right)}_k{\displaystyle \frac{c_1}{R^2+c_n}}\right].$$

(40)

2.2.2. Finite Element Solution

To verify the analytical solution, finite element analysis of bi-layer and multi-layer cylinders subjected to internal pressure was carried out using COMSOL Multiphysics software. To capture hyperelasticity and large deformation, a Mooney-Rivlin material model was used. This model includes two material parameters, denoted as C₁₀ and C_01.

2.2.2.1. Geometry

One quarter of each cylinder was modeled to reduce the simulation runtime. The geometries of the bi-layer and multi-layer cylinders are illustrated in Figure 2.a and Figure 2.b, respectively. No gap or interference was assumed between the layers. The outer layer is highlighted in blue in each case.

The inner and outer radii of each layer for the bi-layer and multi-layer cases are presented in Table 6 and Table 7, respectively.

2.2.2.2. Boundary Conditions

To prevent any undesirable movement of the cylinders, two displacement-type boundary conditions were imposed. The vertical displacements of radii that lie horizontally in parallel to the x-axis were prescribed as zero, and the horizontal displacements of radii that lie vertically in parallel to the y-axis were prescribed as zero. Radial pressure is directly applied to the inner radius of the cylinder. To ensure symmetric results, a symmetry boundary condition was defined for the radii that are parallel to the x-axis and y-axis

3. Results

3.1. The Bi-Layer Cylinder

3.1.1. Increasing Material Properties with Increasing Radius

The radial displacement and stress distributions of the bi-layer cylinder in the case of increasing material properties with increasing radius are shown in Figure 3a, 3.b, 3.c, and 3.d, respectively.

The error between the analytical and finite element solutions in the bi-layer cylinder, in the case of increasing material properties with increasing radius, is reported in Table 8.

3.1.2. Decreasing Material Properties with Increasing Radius

The radial displacement and stress distributions of the bi-layer cylinder in the case of decreasing material properties with increasing radius are shown in Figure 4a, 4b, 4c, and 4d, respectively.

The error between the analytical and finite element solutions in the bi-layer cylinder, in the case of decreasing material properties with increasing radius, is reported in Table 9.

3.2. The Multi-Layer Cylinder

3.2.1. Increasing Material Properties with Increasing Radius

The radial displacement and stress distributions of the multi-layer cylinder in the case of increasing material properties with increasing radius are shown in Figure 5a, 5b, 5c, and 5d, respectively.

The error between the analytical and finite element solutions in the multi-layer cylinder, in the case of increasing material properties with increasing radius, is reported in Table 10.

3.1.2. Decreasing Material Properties with Increasing Radius

The radial displacement and stress distributions of the multi-layer cylinder in the case of decreasing material properties with increasing radius are shown in Figure 6.a, 6.b, 6.c, and 6.d, respectively.

The error between the analytical and finite element solutions in the multi-layer cylinder, in the case of decreasing material properties with increasing radius, is reported in Table 11.

4. Discussion

In this study, two cases of bi-layered and multi-layered thick-walled, functionally graded, hyperelastic cylinders subjected to internal pressure are investigated. The Mooney-Rivlin material model is used to account for the hyperelasticity and large deformation. Distinct material properties are assigned to each layer of the cylinder. The analytical formulation ensures the incorporation of incompressibility conditions, which is not guaranteed when using numerical software such as COMSOL. The continuity condition of the radial displacement is also defined as a constraint to ensure that its values at the shared edges are compatible with each other.

According to Figure 3Figures toFigure 6, the comparison between the distributions of radial displacement and stress components, including radial, tangential, and axial, derived from both the numerical and analytical solutions, shows excellent agreement, which indicates that the proposed analytical solution can properly predict the fields across the layers of the cylinder.

According to the second subfigures of Figure 3Figures toFigure 6, the inner layer, which is in direct contact with the applied pressure, demonstrates the highest amount of stress compared to the outer layers, which are significantly smaller due to their lower material properties. This behavior is observed in both the bi-layer and multi-layer cylinders, indicating that, regardless of the number of layers, more compliant materials lead to lower stress concentration. Both the tangential and axial stresses are tensile, which is typical for pressurized cylinders. Both stress components show discontinuities at the interfaces between layers due to the differences in material properties.

In this work, three values of 0.2, 0.15, and 0.1 MPa for internal pressure are investigated for each material gradation case. Regardless of the number of layers, excellent agreement between the analytical and FEM solutions, with errors below 1% is recorded. As the internal pressure increases, the magnitude of the stress components follow the same trend and increase, and the errors increase slightly. This observation highlights the sensitivity of the error margin to higher internal pressures, but the discrepancies remain small and within an acceptable range.

The findings can provide valuable insights for designing FG pressurized cylinders, as the stress and displacement distributions can be optimized by selecting appropriate material properties across the thickness of the cylinder, which leads to more efficient and durable designs. The proposed closed-form solution can be used by engineers and designers to predict the behavior of multi-layer, functionally graded, hyperelastic cylinders under pressure, without the need for complex and time-consuming numerical simulations.

Engineers should also consider the material gradation to achieve a more uniform stress distribution, which is essential for the structural integrity of the vessels. Multi-layer cylinders with smoothly varying properties across the radius demonstrate reduced stress concentrations and improved overall performance when subjected to internal pressure.

Future research could explore the addition of viscoelastic behavior, where time-dependent deformations and energy dissipation play an important role under cyclic loading conditions. In order to further provide insight into this prediction, incorporation of damage mechanics into the model could be a promising approach for failure prediction. particularly when fatigue or localized material degradation are concerns for designers. Additionally, the use of shape memory behavior could be an interesting avenue for future work, particularly in systems where self-healing or adaptive responses to external loads are desired.

Appendix A.1

The closed-form analytical solution obtained after integrating Eq. (36), as expressed in Eq. (A1), involves the parameters $M={\displaystyle \frac{C_{{10}_{max}}}{C_{{10}_{min}}}}$, $N={\displaystyle \frac{C_{{01}_{max}}}{C_{{01}_{min}}}}$, $C_1=C_{{10}_{max}}$, $C_2=C_{{10}_{min}}$, $C_3=C_{{01}_{max}}$, $C_4=C_{{01}_{min}}$, with $a=R_{in}$, $b=R_{out}$, and $c$ representing the integration constant.