Analytical and Finite Element Solution for Functionally Graded Pressure Vessels Subjected to Finite Strain Coupled Axial and Torsional Deformations

TL;DR

This paper provides an analytical and numerical solution for the mechanical behavior of a hyperelastic cylinder under combined axial and torsional deformations.

Contribution

The study introduces a novel analytical solution and validates it with finite element simulations for functionally graded hyperelastic materials.

Findings

Analytical and numerical results show strong agreement for various stretch and twist levels.

Torsion causes complex longitudinal stress distribution with two extrema under large twisting.

A critical axial stretch determines the transition of axial force from tension to compression during twisting.

Abstract

This study presents an analytical solution to examine the mechanical behavior of an incompressible, functionally graded hyperelastic cylinder under combined extension and torsion. The exp-exp strain energy density function characterizes the hyperelastic material, with parameters varying exponentially along the radial direction. To validate the solution, finite element simulations using a custom UHYPER in ABAQUS are performed. The analytical and numerical results show strong agreement across different stretch and twist levels. The stress distribution and maximum stress are significantly influenced by the exponential parameter governing material gradients. Unlike axial stretch, torsion induces a more intricate longitudinal stress distribution, with large twisting producing two extrema that shift toward the cylinder’s center and outer surface. Longitudinal stress primarily governs von Mises stress and strain energy density variations across the radial direction. A critical axial stretch is identified, below which torsion-induced axial force transitions to compression, elongating the cylinder during twisting. Beyond this stretch, the axial force shifts from tensile to compressive with increasing twist, causing initial shortening before further elongation.

Figures

Full text

1. Introduction

Rubber and rubber-like materials have recently become pivotal in various industrial applications [1,2,3]. They include the development of sensors and actuators [4,5,6,7], soft microfluidic systems [8,9,10], conductive liquid microchannels [11], and tissue engineering scaffolds [12]. The construction industry is also an interesting field of application for soft materials [13,14,15,16]. Due to the extensive use of stretchable materials, such as elastomers, hydrogel-based systems, and rubber-like materials, they have been the focus of numerous studies [17,18,19]. Notably, the combined extension–torsion behavior of nonlinear elastic materials was first explored in foundational work by Rivlin [20,21]. The widespread use of hyperelastic materials, coupled with their mechanical similarities to soft biological tissues, has further garnered significant interest within the research community. Furthermore, the extension–torsion problem is considered a critical experimental method for characterizing the mechanical behavior of materials [22], particularly soft biological tissues, such as the papillary muscles of the heart [23,24,25].

In nonlinear mechanics, the behavior of rubber-like materials is typically described using hyperelastic models, which employ strain energy density functions. Numerous constitutive models, ranging from micro-mechanical frameworks to phenomenological approaches, have been proposed to characterize the mechanical properties of nonlinear elastic materials [26,27,28,29,30]. For instance, Pan and Zhong [31] introduced a viscoelastic constitutive model for rubber materials, which integrates neo-Hookean elasticity and convolution-based viscous components to account for phenomena such as the Payne effect.

In general, an effective strain energy density function should possess a relatively simple mathematical structure, an optimal number of material parameters, compatibility with experimental data, and the ability to predict material behavior under complex, multi-axial loading conditions [32]. In recent years, inspired by pioneering research, numerous studies have focused on proposing or refining strain energy functions to align with the experimental observations of both homogeneous and composite elastomeric materials.

Bechir et al. [33] introduced a modified neo-Hookean model formulated based on the Seth–Hill invariants, supported by experimental investigations to establish correlations between model parameters. Khajehsaeid et al. [34] developed two constitutive models for rubber-like materials, characterized by strain energy functions exhibiting logarithmic and exponential behavior. Building upon the concept of representing principal stretches as exponential series, Darijani et al. [35] formulated a strain energy function tailored for hyperelastic materials. Subsequently, Mansouri and Darijani [36] proposed a hyperelastic constitutive model utilizing strain invariants expressed through exponential functions. This model is particularly noteworthy for its simplicity, incorporation of the second strain invariant, and its ability to align with experimental data for biological tissues and hyperelastic materials. Compared to traditional models such as the neo-Hookean and Mooney–Rivlin formulations, the exp–exp strain energy function offers several key advantages under extreme loading conditions. While neo-Hookean and Mooney–Rivlin models are limited by their polynomial dependence on the strain invariants, the exp–exp model incorporates exponential terms involving both the first and second strain invariants. This structure enhances its ability to capture the pronounced strain-stiffening behavior often observed in elastomeric and biological materials under large deformations. Furthermore, the inclusion of the second invariant allows the exp–exp model to more accurately represent multi-axial deformation states. Studies such as those by Mansouri and Darijani [36] have demonstrated the superior accuracy and stability of the exp–exp model across a wide range of deformation levels, where simpler models often fail to reflect the nonlinear mechanical response. This makes the exp–exp model particularly suitable for applications involving severe mechanical loading, such as those explored in the present study. Additionally, the nonlinear load–stretch relationships observed in other constitutive models [37] further highlight the versatility of these approaches.

Functionally graded materials (FGMs) have recently garnered significant attention and have found widespread application in various engineering systems due to their adaptability and multifaceted functionality. The unique capability of FGMs to exhibit a gradual variation in material properties, along with the ability to control this gradient, makes them an effective solution for optimizing stress distribution under combined thermal and mechanical loads. These attributes have positioned FGMs as a promising technology for addressing complex engineering challenges. Moreover, the exceptional mechanical performance of hyperelastic materials and their extensive applications across scales, from nano to macro, have driven substantial research into understanding the mechanical behavior of FGMs composed of hyperelastic materials. For instance, Anani and Rahimi [38] examined the stress distribution within a functionally graded hyperelastic pressurized thick spherical shell using a modified neo-Hookean model. Similarly, Bilgili [39] investigated the stress–strain inhomogeneity in a functionally graded hyperelastic hollow tube, incorporating models that accounted for temperature effects, strain stiffening, and the radial variation in the shear modulus within the frameworks of the generalized neo-Hookean and Gent models. Building upon these efforts, Anani and Rahimi [40] derived a closed-form solution for stress distribution in functionally graded hyperelastic spherical shells, utilizing a reinforced neo-Hookean model calibrated against experimental data. In subsequent work, they analyzed the behavior of rotating cylindrical shells using a power-law strain energy function tailored to functionally graded hyperelastic materials [41]. Furthermore, Anani and Rahimi [42] proposed a novel constitutive model for describing the visco-hyperelastic behavior of transversely isotropic functionally graded rubber materials, incorporating an equilibrium spring and Maxwell model. Experimental validations demonstrated the model’s efficacy in predicting the time-dependent deformation of rubbers. Moallemi et al. [43] conducted stress and stability analyses of functionally graded (FG) pressure vessels composed of hyperelastic materials under large deformations, using an exponential–exponential strain energy function. Additionally, the analytical solution offers clear, closed-form expressions that enable rapid assessment and deeper insight into the mechanical behavior under various loading conditions, making it computationally efficient compared [44,45] to full-scale numerical simulations. Almasi et al. [46] developed an analytical solution to describe the behavior of FG hyperelastic cylinders subjected to thermomechanical loading, focusing on pressurized vessels. These studies collectively highlight the growing interest in and the advanced modeling techniques like machine learning methods [47] for characterizing functionally graded hyperelastic materials, emphasizing their importance in modern engineering applications.

The extension–torsion problem has been extensively investigated in the field of mechanics. In contrast to uniaxial extension or compression alone [48,49,50,51], combined extension–torsion loading activates multiple deformation modes, providing critical insights for robust material characterization. Humphrey et al. [52] analyzed the large deformations resulting from the combined effects of extension and torsion in a cylindrical structure composed of a homogeneous, transversely isotropic, hyperelastic material. Their study utilized experimental data from papillary muscles to calibrate the constitutive model. Ogden and Chadwick [53] also addressed this problem within the framework of the well-established Ogden model. Focusing on the strain-stiffening behavior of rubber-like materials, Kanner and Horgan [54] explored the extension–torsion of a cylindrical body and conducted comparative analyses of various constitutive models, including those proposed by Gent [29] and Fung [55]. Horgan and Murphy [56] further contributed to this field by investigating the problem using the Varga model [57]. Additionally, the superposition of torsion on extension has been studied in the context of transversely isotropic hyperelastic materials and hydrogels, with particular emphasis on the mechanical properties of soft biological tissues [58,59,60]. A comprehensive comparative summary of these studies, highlighting their methodologies, material assumptions, and key findings, is provided in Table A1 in Appendix A.

Numerous studies have extensively explored the linear elasticity of FGMs; however, research on the mechanical behavior of hyperelastic FGMs remains limited. Specifically, investigations addressing the combined extension and torsion loading of FG rubber-like materials are notably scarce. To address this gap, the current study examines the extension and torsion behavior of an FG circular cylinder composed of a nonlinear, incompressible, isotropic material. This analysis is particularly relevant for applications involving soft biological tissues. The study employs the invariant-based exponential–exponential energy density function proposed by Mansouri et al. [36], chosen for its mathematical simplicity, model stability when incorporating the first and second strain invariants, and its demonstrated agreement with experimental data. Additionally, finite element (FE) analysis is performed to validate the proposed analytical solutions. This is achieved by developing a user-defined hyperelastic material subroutine (UHYPER) within the ABAQUS 2024 software environment, allowing for a comprehensive comparison between the analytical and numerical results.

The research utilizes the invariant-based exponential–exponential strain energy density function introduced by Mansouri et al. [36], selected for its mathematical elegance, stability in incorporating the first and second strain invariants, and its strong correlation with experimental observations. To further substantiate the proposed analytical solutions, FE simulations are conducted. This is accomplished by implementing a user-defined material subroutine (UHYPER) within the ABAQUS software framework, enabling a thorough comparison between the analytical and numerical findings. The paper is structured as follows: Section 2 presents the foundational aspects of constitutive modeling and the kinematics of the extension–torsion problem. Subsections within this section introduce a closed-form analytical solution for the invariant-based strain energy function, specifically the exponential–exponential model utilized in this study. Section 3 outlines the finite element modeling approach for the combined extension and torsion of an FG hyperelastic cylinder. Section 4 discusses the results and interpretations derived from the investigation. Finally, Section 5 provides a summary and concludes the study.

2. Extension–Torsion of Cylinder Composed of Hyperelastic Material

In this section, we develop a detailed analytical derivation for the coupled extension–torsion behavior of a functionally graded hyperelastic cylinder. We start by introducing the fundamental concepts of nonlinear elasticity, defining the relevant kinematic quantities and strain energy invariants that govern the material response. To aid the reader in navigating the numerous equations and the overall logical progression of the derivation, a concise flowchart—summarizing the key steps from problem definition (geometry and loading conditions) to the final expressions for stresses and resultant forces—is provided in Appendix B.

In the field of nonlinear elasticity, the mechanical behavior of hyperelastic materials is characterized by the incorporation of strain energy density functions $[eqn]$ . The foundational principles of nonlinear elasticity are established by considering $[eqn]$ as the deformation gradient. The left Cauchy–Green deformation tensor is expressed as $[eqn]$ . By adopting strain energy invariants, a constitutive law for hyperelastic materials can be formulated as follows:

$$[eqn]$$

in which, tensor $[eqn]$ stands for second Piola–Kirchhoff stress. Assuming $[eqn]$ , we have the following:

$$[eqn]$$

in which $[eqn]$ stands for the princial Cauchy stresses. Considering the definition of $[eqn]$ , the strain invariants are as follows:

$$[eqn]$$

where $[eqn]$ , and $[eqn]$ stand for the principal stretches. Regarding the incompressibility constraint, the volume of the cylinder remains constant, which implies $[eqn]$ .

Although soft materials such as hydrogels can absorb large amounts of water and exhibit significant volumetric changes during the shape shifting [61,62], elastomers and rubber-like materials are commonly treated as incompressible and isotropic in the literature. Thus, considering a circular cylinder composed of incompressible isotropic hyperelastic materials subjected to a tensile stretch in the longitudinal direction and a torsional twist, the undeformed and deformed coordinates can be assumed as ( $[eqn]$ ) and ( $[eqn]$ ), respectively. Therefore, we have the following [63]:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the axial stretch and twist per unit stretched length, respectively. It is also notable that $[eqn]$ originated from satisfying the incompressibility constraint. For stretches larger than $[eqn]$ , the cylinder contracts in radial direction, besides the pure torsion assumptions, are obtained when $[eqn]$ .

Considering Equation (4) as a deformation field, the total deformation gradient tensor $[eqn]$ of the problem of extension–torsion is defined as follows:

$$[eqn]$$

Furthermore, left Cauchy–green deformation tensor, $[eqn]$ , and its inverse $[eqn]$ are as follows:

$$[eqn]$$

Accordingly, the strain invariants may be computed as follows:

$$[eqn]$$

The Cauchy stress for an incompressible hyperelastic material is introduced as follows:

$$[eqn]$$

where $[eqn]$ denotes the second order identity tensor multiplied by volumetric pressure $[eqn]$ . Considering the components of Cauchy stress, the resultant axial force and torsional moment can be found as follows:

$$[eqn]$$

In addition, in the current cylindrical configuration, the equilibrium equation in radial direction is defined as follows:

$$[eqn]$$

Assuming an invariant-based strain energy function, in light of Equations (6), (8) and (10) and considering the stress-free boundary condition $[eqn]$ at $[eqn]$ , the Cauchy stresses are obtained as follows:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Moreover, the generated moment and axial force are computed utilizing the stresses introduced in Equations (13) and (14).

$$[eqn]$$

Among the invariant-based constitutive models, an exponential strain energy density developed by Mansouri and Darijani [36] is employed in the present work because of its stability and great correlation with experiments.

$$[eqn]$$

in which $[eqn]$ , and $[eqn]$ are material parameters. The conformity between the results of this model with the various experimental data was studied for hyperelastic materials considering incompressibility and even compressibility assumptions. In this article, a four parameters form of the exp-exp strain energy function is employed to characterize the material behavior of an FG hyperelastic cylinder in the following scheme:

$$[eqn]$$

It should be noted that $[eqn]$ , and $[eqn]$ are parameters specifically modified for FGMs with material variation along the radial direction. It is clear that for the homogeneous materials, these parameters are no longer a function of radius.

As the modified material parameters for FGMs introduced in Equation (17) vary in the radial direction, the material parameters considered in this paper are assumed to alter along the radial direction in an exponential fashion as follows:

$$[eqn]$$

where $[eqn]$ stands for any material parameter in Equation (17) and subscripts $[eqn]$ and $[eqn]$ denote the parameters correspondent with the materials at the center and outer surface of the cylinder, respectively.

Henceforth, we further investigate the influences of torsion on the generated total axial force of a hyperelastic cylinder. Assuming $[eqn]$ stands for the amount of axial force generated by axial stretch $[eqn]$ in the absence of torsional twist ( $[eqn]$ ), $[eqn]$ is the portion of axial force created by the applied torsional twist. Thus, the total axial force generated by the external extension–torsion loading is defined as follows:

$$[eqn]$$

Taking into account the axial load generated by the prescribed deformation, it appears that the component of axial force resulting from the external torsional twist is critical for determining both the direction and magnitude of the total longitudinal force. To explore the contribution of axial force due to torsion, the rate of change in the longitudinal force in a hyperelastic cylinder with respect to the applied torsional twist is analyzed. By employing the exponential–exponential strain energy function and Equation (15), the torsional moment induced within the cylinder, as well as the total axial force, can be determined as follows:

$$[eqn]$$

3. Finite Element Analysis

Considering that finite element (FE) analysis for the extension–torsion problem of a circular cylinder made of functionally graded (FG) materials has not been previously explored, the present study introduces a three-dimensional FE model to validate the proposed analytical solutions. A user-defined UHYPER was developed to define the strain energy function and the corresponding material parameters. The 3D solid cylinder is subdivided into multiple radial layers (e.g., 10, 20, 40, and 80 layers), with each layer assigned a specific set of material properties based on its radial position. In the material module, different properties are defined for each layer and assigned to the corresponding section. The material parameters vary exponentially along the radial direction, as presented in Equation (18).

Figure 1 illustrates the schematic of the 3D FE model and its mesh structure, employing 3D continuum 8-node hexahedral elements (C3D8H). The selection of the optimal number of layers to ensure an accurate and reliable FE analysis of the extension–torsion problem is addressed in the subsequent section. To confirm the convergence of the FE model, mesh independence studies are conducted, and the results are presented.

A prescribed axial displacement and rotation are applied on one end of the cylinder, while the opposite end is fully constrained. Finally, the axial force and torsional moment are computed along the cylinder’s axis using a static general analysis step.

4. Results and Discussion

The analytical solution proposed in this study is applied to the exponential–exponential (exp-exp) constitutive model to assess the stress distribution, as well as the resultant total axial force and moment, in an FG hyperelastic cylinder subjected to a combination of axial extension and torsional loading. Additionally, FE analysis results are employed to validate the accuracy of the analytical solutions. For the present investigation, the hyperelastic materials selected for the core and outer surface of the FG cylinder are VHB 4905 and natural rubber gum, respectively.

These two materials were chosen deliberately due to their markedly different mechanical behaviors, which are ideal for constructing a functionally graded system with a pronounced property gradient. VHB 4905 is a soft elastomer that exhibits notable strain-softening characteristics, while natural rubber gum is stiffer and shows strain-hardening behavior. This contrast enables the model to capture the radial variation in mechanical response more effectively. Furthermore, both materials are well studied in the literature, with available experimental data used for calibrating the parameters of the exp–exp strain energy model. The corresponding material parameters for the exp-exp model are provided in Table 1, and these values are determined using a least squares regression approach [37]. The close match between model predictions and experimental data confirms the reliability of the chosen materials for validating the proposed theoretical framework.

In Equation (18), parameter $[eqn]$ controls the form and intensity of the material variation throughout the radial direction. Various values of $[eqn]$ are considered to study its effect on the stress distribution. For instance, the variation in $[eqn]$ throughout the radial direction is illustrated in Figure 2. It is worth mentioning that the distribution of other material parameters has a similar trend to that of $[eqn]$ .

Figure 3a illustrates the variation in radial stress throughout the radius for $[eqn]$ , and $[eqn]$ considering different numbers of layers in FE modeling of the FG cylinder. Obviously, increasing the number of strips in FGM modeling increases the accuracy and reliability of the results. However, considering the fact that increasing the number of layers results in time-consuming computations, an optimum number of layers that guarantees good agreement with the analytical results should be identified to reduce the computational cost. It is evident in Figure 3 that for 80 layers, an excellent agreement between the analytical solutions and FE analysis is achieved. In addition, it is apparent that the radial stress begins from an extremum value appearing at the center and becomes zero at the outer surface, which reflects the stress-free condition at outer radius $[eqn]$ . The variation in the dimensionless hoop and longitudinal stress for $[eqn]$   $[eqn]$   $[eqn]$ , and various number of layers are plotted in Figure 3b,c, respectively. Similarly to the radial stress, for 80 layers, the results are in great agreement with those of the analytical calculations. Although the longitudinal stress is positive throughout the radial direction, approaching the exterior surface from the central point of the cylinder, the hoop component of the stress coverts from negative to positive.

We further discuss the effect of $[eqn]$ on the stress distribution. Figure 4 indicates the dimensionless radial stress for $[eqn]$ , and different values of $[eqn]$ . It is observed that the highest radial stress occurs at $[eqn]$ and at higher values of $[eqn]$ , the radial stress decreases at the entral point of the cylinder. Moreover, the radial stress vanishes at the exterior surface of the circular cylinder regardless of the value of $[eqn]$ satisfying the stress-free boundary condition at the outer radius $[eqn]$ . Figure 4b demonstrates the distribution of the dimensionless hoop stress along the dimensionless radius for $[eqn]$ , and various values of $[eqn]$ . Regardless of the value of $[eqn]$ , the circumferential stress alters from compressive to tensile. However, the location of the alteration point depends on the value of $[eqn]$ . Figure 4c depicts the variation in longitudinal stress for $[eqn]$ , and various magnitudes of $[eqn]$ . It is perceived that for $[eqn]$ , the longitudinal stress reaches to its maximum approximately at $[eqn]$ while for $[eqn]$ , the point where the maximum longitudinal stress takes place moves toward the outer surface and for $[eqn]$ , it shifts toward the center. This is mainly due to the different distribution of materials through the radial direction for the positive and negative values of $[eqn]$ according to Figure 2. Figure 4d plots the distribution of the dimensionless strain energy density and von Mises stress throughout the radial direction for $[eqn]$ , and various magnitudes of $[eqn]$ . Similarly to the longitudinal stress, the maximum strain energy density and von Mises stress for positive values of $[eqn]$ appears close to the outer surface while for negative $[eqn]$ , it moves toward the center. This similar trend indicates that the contribution of the longitudinal stress to the total energy or von Mises stress is considerably higher than the contribution of radial and hoop stresses.

Considering a constant value for $[eqn]$ , we turn the focus on the effect of torsional twist and longitudinal stretch on the stress distribution. Figure 5 indicates the radial stress distribution for $[eqn]$ , and various twists. It is apparent that increasing the twist escalates the rate of the stress variation and results in higher radial stress. The distribution of tangential and axial components of stress along the radial direction for $[eqn]$ , and different magnitudes of twists are illustrated in Figure 5b,c, respectively. It is observed that regardless of the value of twist, the hoop stress distribution plots intersect approximately at $[eqn]$ and the hoop stress alters from the compressive value to tensile. Referring to Figure 5c, increasing the twist has a more intricate effect on the variation in longitudinal stress. For $[eqn]$ , the maximum longitudinal stress occurs approximately at $[eqn]$ , while for larger twists, two extrema including a maximum and a minimum are identified in the stress variation plot. By increasing the twist, the maximum point moves toward the center whereas the minimum point shifts toward the outer surface. Figure 5d depicts the distribution of the strain energy density function and also the von Mises stress along the radial direction for $[eqn]$ , and various twists. It is revealed that the squared value of longitudinal stress constructs a large portion of von Mises stress and the strain energy density function.

Figure 6 depicts the distribution of radial stress for $[eqn]$ and $[eqn]$ and varying amounts of stretch. Increasing the value of axial stretch has a similar effect as the torsional twist and intensifies the radial stress. The influence of external axial stretch applied on the cylinder on the variation in hoop and longitudinal stress for $[eqn]$ and $[eqn]$ is illustrated in Figure 6b,c, respectively. For the hoop stress, the trend is analogous with the influence of the twist. Regardless of the axial stretch, the hoop stress distribution plots intersect roughly at $[eqn]$ where the tangential component of stress changes from the compressive state to tensile. However, on the contrary to the complicated effect of the twist on the longitudinal stress, the location of the maximum axial stress appears approximately at $[eqn]$ regardless of the magnitude of the axial stretch. The variations in the strain energy density and von Mises stress are displayed in Figure 6d. The maximum value of strain energy function and also von Mises stress is dominated by the longitudinal component of stress, where the form of variations are almost identical.

To further approve the accuracy of the proposed analytical procedure and provide a better visual comprehension, von Mises stress contours obtained from analytical and finite element analysis are compared. Figure 7a illustrates the von Mises stress contours for a sector of the cross-section for different values of $[eqn]$ , which conforms with the results presented in Figure 4d. Moreover, in Figure 7b, the von Mises stress contours for $[eqn]$ , and different values of $[eqn]$ and $[eqn]$ are demonstrated, which recast the results shown in both Figure 5d and Figure 6d.

Figure 8a represents the variation in dimensionless total force along the axial versus the twist for different amounts of $[eqn]$ . The total force is significiantly diminished for larger twists regardless of the value of $[eqn]$ . However, the rate of variation is smooth for a negative $[eqn]$ , while it is intensified for zero or positive values of $[eqn]$ . The variation in the torsional moment versus the twist is demonstrated in Figure 8b. Similarly to the trend of axial force, for negative values of $[eqn]$ , a smooth rate of variation is observed while for zero and positive values of $[eqn]$ , the rate of momentum variation increases.

The variation in dimensionless total axial force and generated moment versus the torsional twist for $[eqn]$ and varying axial stretch is depicted in Figure 9a,b, respectively.

Differentiating from $[eqn]$ in Equation (20), the position where the axial force generated from twist alters from the tensile value to compressive is identified as the alteration point. For $[eqn]$ , $[eqn]$ is always compressive, which means the cylinder tends to elongate upon twisting. However, for $[eqn]$ , a different behavior is observed and for relatively small twists, $[eqn]$ is initially tensile, which temporarily increases the total axial force, while on further twisting, $[eqn]$ becomes compressive and the total axial force decreases. In other words, the cylinder compresses for sufficiently small twists and then elongates on further twisting. Figure 10 indicates the analytical solutions for the variation in dimensionless axial force versus the twist even including large deformations for various values of $[eqn]$ . The circular markers indicate the points where $[eqn]$ and for a constant value of $[eqn]$ , by increasing the axial stretch, this point appears in larger twists.

In addition to the theoretical contributions of this study, our findings offer valuable insights into practical material design and structural applications. The developed analytical model not only deepens our understanding of the complex extension–torsion behavior in functionally graded hyperelastic cylinders but also provides a robust tool for material characterization. By correlating experimental test data with the model’s predictions, engineers and researchers can efficiently extract key constitutive parameters, thereby facilitating the design and optimization of graded elastomeric components. Potential applications include the development of soft actuators, flexible sensors, and biomedical devices—such as vascular grafts and tissue scaffolds—where precise prediction of stress distributions under complex loading is critical. Moreover, integrating our analytical model within numerical simulation workflows can streamline the design process, reduce prototyping cycles, and ultimately lower development costs for advanced soft functional materials and medical devices. This dual utility, for both design and characterization, underscores the real-world impact of our theoretical framework. These results highlight the model’s potential applicability in soft tissue engineering (e.g., vascular grafts, scaffolds) and in the design of soft actuators or sensors in robotics, where accurate modeling of complex, multi-axial deformations is essential.

5. Summary and Conclusions

The current investigation examines the mechanical behavior of an incompressible, isotropic, functionally graded hyperelastic circular cylinder subjected to combined extension and torsion. The material composition of the cylinder was assumed to vary exponentially in the radial direction. Specifically, the core of the cylinder was modeled using VHB 4905, while the outer surface was represented by natural rubber gum. An exponential strain energy density function was employed in this study due to its mathematical simplicity, stability, and excellent agreement with experimental observations under a variety of loading conditions. This model was used to predict the behavior of the FG hyperelastic material constituting the cylinder. To assess the proposed analytical solutions, FEA was carried out using an ABAQUS user-material subroutine (UHYPER v6.6). A parametric study was conducted to determine the optimal number of layers required to accurately model the FG cylinder, followed by a mesh independency analysis to ensure reliable and precise FEA results.

The findings demonstrated that the distribution of stress and the location of maximum stress were significantly influenced by the exponential exponent in the material variation function. Radial, tangential, and longitudinal stress distributions were analyzed for various levels of torsional twist and longitudinal stretch. Radial stress, regardless of the twist or stretch magnitude, exhibited a compressive maximum at the center and diminished to zero at the outer surface. Additionally, the hoop stress distribution, irrespective of the twist or axial stretch, intersected at approximately R/Rout = 0.53, where the hoop stress transitioned from compression to tension. While the axial stretch did not influence the position of maximum longitudinal stress, the effect of torsion was more pronounced. For a torsional twist of τ = 1, the maximum longitudinal stress occurred near R/Rout = 0.5. As the twist increased, two extrema appeared in the longitudinal stress distribution, moving toward both the center and outer surface of the cylinder. The longitudinal stress primarily governed the variation and magnitude of the von Mises stress and strain energy density along the radial direction.

For negative values of the material exponent (k), the variation in the total moment and axial force was relatively smooth. In contrast, for positive values of k, the rate of variation in these quantities increased. Furthermore, the total force in the longitudinal direction transitioned from an initial tensile to a final compressive state due to the torsional effect. A critical alteration point was identified for the axial stretch, where ∂N/∂τ = 0. For strains smaller than this alteration value (γ < γ_alteration_), the axial force resulting from torsion, $[eqn]$ , remained compressive, suggesting that the cylinder tended to elongate under external twisting. However, for strains greater than the alteration value (γ > γ_alteration_), $[eqn]$ transitioned from tensile to compressive, indicating that the cylinder first shortened and then elongated as the torsional twist increased.

In summary, future work may extend the present study in several directions. For instance, dedicated experimental tests using commercially available materials like VHB 4905 could further validate and refine the model. Additionally, incorporating non-axisymmetric loading conditions and geometric imperfections would enhance the practical applicability of the analysis, while extending the formulation to accommodate anisotropic material behavior could better capture the complex responses observed in real-world applications. Furthermore, exploring more complex or non-monotonic gradient functions—beyond the simple exponential variation considered here—and integrating damage or failure criteria into the model represent promising avenues to optimize material performance and assess durability under extreme conditions.