Investigating the Triaxial Mechanical Behaviour of Silicone Rubber Material

TL;DR

This paper introduces a new model to better understand how silicone rubber behaves under complex mechanical stresses.

Contribution

A novel constitutive model is proposed to capture the nonlinear triaxial mechanical behavior of silicone rubber.

Findings

The model accurately represents nonlinear large-deformation behavior and sensitivity to biaxiality ratio.

Simulation results align well with experimental data in deformation patterns and stress–strain curves.

The model balances parameter count and fitting accuracy better than classical hyperelastic models.

Abstract

Silicone rubber is extensively used in engineering applications due to its toughness and impact resistance; however, traditional characterisation methods fail to capture its nonlinear deformation characterisation and triaxial mechanical behaviour. To address this, we derived a constitutive model within the framework of continuum mechanics that assumes a condition of near incompressibility and conducted uniaxial, planar, and equibiaxial tension tests to fit the model parameters. Through systematic analysis of triaxial mechanical responses under these three loading modes, we determined the material’s nonlinear large-deformation behaviour and sensitivity to the biaxiality ratio. Comparative analyses with classical hyperelastic models show that the proposed model achieves a good balance between the number of parameters and fitting accuracy. After the parameter-fitting process, we performed finite element simulations of the three loading modes. The simulation results show good agreement with experimental data in terms of deformation patterns and stress–strain curves. This study provides a novel theoretical tool for evaluating the mechanical properties and structural designs of soft materials.

Figures

Full text

1. Introduction

Elastomeric polymers exhibit extensive application prospects in the fields of transportation, security, medical health, sports, and many others due to their exceptional mechanical properties, including high elasticity (up to 100% or even 1000%), toughness, and impact resistance [1,2]. Of these, silicone rubber has also become an ideal material for numerous practical applications, such as soft robotics [3] and flexible energy harvesters [4,5], due to its high extensibility and durability, low dissipation behaviour, and biocompatibility. These practical applications require the material to withstand deformations under complex mechanical states, highlighting the need to predict silicone rubber’s complex mechanical behaviour. Studies have demonstrated that silicone materials exhibit significant sensitivity to the biaxiality ratio $[eqn]$  [6]. This variable is defined as $[eqn]$ , where $[eqn]$ and $[eqn]$ represent the maximum and minimum principal stretches, respectively, and serves as a parameter characterising loading conditions and reflecting the degree of inhomogeneity in the deformation states experienced by materials. Biaxial mechanical behaviour is an important research direction in materials science and engineering, as it focuses on investigating material responses under different deformation modes.

In recent years, many studies have been conducted on the mechanical behaviour of silicone rubber materials. Liao et al. [7] performed systematic experiments to decouple the material’s mechanical response characteristics, revealing that it exhibits weak viscoelastic characteristics but significant stress softening and recovery behaviours. Kumar et al. [8] investigated the static and dynamic mechanical properties of polydimethylsiloxane (PDMS) under uniaxial tensile conditions, examining its viscoelastic behaviour under varying strain rates and force-controlled cyclic loading. However, the existing studies provide limited descriptions of mechanical behaviour under varying biaxiality-ratio conditions. Some researchers have also developed biaxial tensile testing systems for biomaterials, with PDMS materials used to validate the equipment’s accuracy and stability, e.g., Roth et al. [9] and Corti et al. [10]. Gao et al. [11] studied the mechanical behaviour of knee joint cartilage under biaxial cyclic loading conditions and characterised the associated ultrastructural changes. Their findings indicated that constraints in the Y-direction (orthogonal direction) could effectively reduce the maximum tensile strain in the X-direction. Liao et al. [6] performed cyclic loading–unloading tests on commercial silicone rubber under various deformation modes (uniaxial, planar, and equibiaxial tension) across multiple strain levels, revealing the material’s mechanical behaviour and stress recovery characteristics under repeated loading. Luo et al. [12] employed finite element analysis to evaluate the accuracy and comparability of three equibiaxial stretching methods: inflation stretching, equibiaxial plane stretching, and radial stretching.

Research on the mechanical properties of rubber-like materials continues to develop in tandem with theoretical frameworks of their physical constitutive relations. A central focus has developed: quantitatively characterising the intrinsic correlation between strain and stress in these materials while fully accounting for their critical mechanical characteristics, such as biaxiality-ratio sensitivity. Hyperelastic constitutive relations are a crucial theoretical framework for describing their nonlinear large-deformation mechanical behaviour. This theory posits the existence of a Helmholtz free energy function $[eqn]$ , which is a continuous scalar function dependent on the deformation gradient tensor $[eqn]$ , capable of characterising the material’s energetic state during deformation processes. According to Holzapfel et al. [13], the strain energy function $[eqn]$ can be expressed through three independent invariants. Beyond uniaxial deformation, scholars have investigated constitutive models under various deformation modes. Ogden [14] demonstrated that the Neo-Hookean model [15] can only describe shear deformation in rubber within small-to-moderate strain ranges, whereas the Mooney–Rivlin model [16,17] incorporating the second invariant demonstrates broader applicability across diverse deformation modes. Xiao et al. [18] found that conventional isotropic damage models fail to capture the stress response of tough gels under multiaxial loading; they, therefore, developed a non-affine theory based on the microsphere model, which successfully predicted experimental results from pure shear and unequal biaxial tests in addition to damage cross-effects. Ostadrahimi et al. [19] proposed a physics-informed machine learning framework to simulate nonlinear, history-dependent viscoelastic mechanical behaviour under multiaxial cyclic loading conditions.

Based on these findings, we can summarise the current limitations in the field using the following categories: most experiments focus on uniaxial loading, with the use of systematic mechanical testing and characterisation techniques for complex triaxial stress states being relatively scarce; and existing models are mainly based on parameter fitting from uniaxial tests, neglecting the complex effects of biaxiality-ratio sensitivity inherent to elastometric polymers.

In-depth investigations, theoretical quantitative descriptions, and numerical implementations of the mechanical behaviour of silicone rubber, therefore, hold fundamental significance for effectively predicting the deformation, optimising the performance, and promoting the engineering applications of related materials.

This study focuses on the mechanical behaviour of silicone rubber under different deformation modes, aiming to thoroughly investigate the biaxiality-ratio-sensitive characteristics of related materials by developing a constitutive model. To achieve this objective, we designed and conducted uniaxial, planar, and biaxial tension tests; parameter fitting; and comparative analyses to evaluate the model’s accuracy regarding the nonlinear large deformation and biaxiality-ratio sensitivity of the material. Our results hold significance for enhancing precise simulation analyses of elastometric polymers under various deformation modes.

2. Constitutive Modelling

In this section, based on continuum mechanics, we develop a constitutive model to characterise the nonlinear large deformations and sensitivity-to-biaxiality ratio of the elastometric material.

2.1. Kinematics

First, consider a continuum body $[eqn]$ in three-dimensional Euclidean space, and introduce an initial coordinate system with the origin O and basis vectors $[eqn]$ , as shown in Figure 1. The continuum $[eqn]$ undergoes motion from its initial configuration (reference configuration) $[eqn]$ to the current coordinate system $[eqn]$ . When any material particle $[eqn]$ moves from position $[eqn]$ to position $[eqn]$ , the deformation gradient is given as

$$[eqn]$$

The determinant of the deformation gradient,

$$[eqn]$$

describes the volume change of the differential element, where V is the initial volume of the differential element, and v is the current volume. The right Cauchy–Green strain tensor is defined as

$$[eqn]$$

2.2. Clausius–Duhem Inequality

According to the Clausius–Duhem inequality [13], we have

$$[eqn]$$

where $[eqn]$ denotes the second Piola–Kirchhoff stress tensor; $[eqn]$ , $[eqn]$ , and $[eqn]$ represent the rate forms of the right Cauchy–Green strain tensor, strain energy function, and temperature, respectively; e signifies internal energy; $[eqn]$ represents the heat flux vector; and $[eqn]$ indicates the temperature gradient. The above equation establishes the thermodynamic relationship governing the mechanical response between stress and strain. Substituting the rate form of the strain energy function into Equation (4) yields

$$[eqn]$$

Considering the arbitrariness of $[eqn]$ , $[eqn]$ can be written as

$$[eqn]$$

Similarly, the relationship between the first Piola–Kirchhoff stress tensor $[eqn]$ and the strain energy function $[eqn]$ can be written as

$$[eqn]$$

2.3. Nearly Incompressible Framework

Further studies have shown that elastomeric polymers exhibit a certain degree of compressibility. For example, in uniaxial tensile tests, vulcanised and pure rubber exhibit volume changes of 1% and 2%, respectively, at deformations of up to 700% [20,21]; in other deformations, the change is even lower [22]. Some softer materials, such as polyurethane and hydrogels, exhibit greater compressibility [23,24]. It is necessary to consider the nearly incompressible nature of elastomeric polymers; not only can this more accurately describe their mechanical behaviour (even slight volume changes can have a significant effect), but it can also alleviate numerical difficulties and instabilities caused by “incompressibility” in finite element simulations [25].

For nearly incompressible models, the strictly incompressible model is modified into an isochoric strain energy function with an added volumetric deformation part, known as the “penalty function” method [25]. In this framework, the deformation gradient is multiplicatively decomposed into volume-changing (dilational) and volume-preserving (distortional) parts, i.e.,

$$[eqn]$$

where $[eqn]$ is the modified deformation gradient; its determinant can be derived as $[eqn]$ , which demonstrates the volume-preserving characteristic of this tensor. Under this decomposition approach, subsequent deformation analysis requires the use of the isochoric form of the right Cauchy–Green strain tensor, i.e.,

$$[eqn]$$

In a nearly incompressible framework, the strain energy function of elastomeric polymers can be decomposed into isochoric and volumetric parts, i.e.,

$$[eqn]$$

where $[eqn]$ is the isochoric strain energy function, and

$$[eqn]$$

represents the isochoric principal invariants. Note that $[eqn]$ in Equation (11) denotes the trace operation. In Equation (10), the volumetric strain energy function

$$[eqn]$$

is commonly referred to as the penalty function [25], which helps improve computational stability in numerical simulations. The parameter $[eqn]$ is related to the bulk modulus and is typically large for elastomeric polymers. Determining the volumetric parameter usually requires volumetric deformation experiments to be performed; however, research in this area is relatively limited. Ogden [22] suggested introducing a coupling strain energy term between the isochoric and volumetric parts to further improve the accuracy of the problem description. However, doing so removes the advantageous separability of the stress and elasticity tensors in such constitutive relations and increases the difficulty of determining the coupling term parameters. Therefore, for simplicity in problem handling, we adopt an additive decoupling approach for the strain energy function.

In Equation (6), the second Piola–Kirchhoff stress tensor can be similarly decomposed into isochoric and volumetric parts, i.e.,

$$[eqn]$$

Isochoric stress can be expressed as

$$[eqn]$$

and volumetric stress can be expressed as

$$[eqn]$$

To facilitate understanding of Equation (14), we explain the differential process of $[eqn]$ below. Considering

$$[eqn]$$

we have

$$[eqn]$$

where $[eqn]$ denotes the fourth-order identity tensor

$$[eqn]$$

and $[eqn]$ is the Kronecker delta function

$$[eqn]$$

Note that the symbol ⊗ in Equation (14) and Equation (18) denotes a tensor product operator. The transpose of the fourth-order tensor $[eqn]$ in Equation (17) is the fourth-order projection tensor

$$[eqn]$$

In Equation (14), $[eqn]$ is the inverse of $[eqn]$ , and

$$[eqn]$$

is the modified second Piola–Kirchhoff stress tensor for nearly incompressible problems. The tensor $[eqn]$ in the above equation is the second-order identity tensor defined as

$$[eqn]$$

Equation (14) is the complete expression of the isochoric second Piola–Kirchhoff stress tensor when the strain energy function $[eqn]$ depends on the isochoric principal invariants $[eqn]$ . Note that in Equation (15),

$$[eqn]$$

is the hydrostatic pressure. Note that the scalar p may only be determined through the equilibrium equations and the boundary conditions [13].

According to the above decomposition process, it can be seen that different forms of the strain energy function $[eqn]$ only change $[eqn]$ and p. These stresses are key tensors in the subsequent process of deriving constitutive models using specific strain energy functions. Of course, the premise of determining these stress tensor expressions is to identify each specific strain energy function and then differentiate between them to determine the stress expression. At this point, a three-dimensional large deformation constitutive model framework has been constructed. However, since direct fitting in three-dimensional form is impossible during the parameter-fitting process, the constitutive model needs to be simplified to a one-dimensional form. Therefore, in the next section, we will derive the one-dimensional general solution form of the three-dimensional constitutive model. Since it is difficult to solve volumetric stress in the one-dimensional form directly, it can only be temporarily obtained using a non-direct method under a nearly incompressible framework. Previous studies [26] have shown that this simplification has little effect on the parameter–fitting accuracy.

2.4. Constitutive Relationship Under Three Deformation Modes

In this section, we derive stress expressions under three typical deformation modes: uniaxial (UT), planar (PT), and equibiaxial tension (ET). According to the definition of the biaxiality ratio $[eqn]$ , the three deformation modes respectively correspond to $[eqn]$ , $[eqn]$ , and $[eqn]$ [6]. The superscripts represent different deformation modes. The derivation in this section gives the one-dimensional stress expression in parametric form, which provides a theoretical basis for the subsequent parameter fitting.

2.4.1. General Solution Form

Since the stresses in the three deformation modes are usually given in the form of the first Piola–Kirchhoff stress tensor P (nominal stress) during experiments, the general solution is given in the form of P below. Considering the principal direction deformation, the componential formulation of P can be expressed as

$$[eqn]$$

According to Equation (24), general expressions under the three deformation modes are

$$[eqn]$$

where subscripts 1 and 2 represent the loading directions. The methods for deriving the general expression under each of the three deformation modes are presented below.

Uniaxial Tension, UTThe specimen is subjected to tensile loading exclusively along the 1-direction. Therefore, under the uniaxial tension mode, the deformation gradient $[eqn]$ and the right Cauchy–Green strain tensor $[eqn]$ are given by

$$[eqn]$$

where $[eqn]$ represents the stretch ratio in the 1-loading direction. The invariants of $[eqn]$ are defined as

$$[eqn]$$

Since the lateral contraction is free during uniaxial deformation, the first Piola-Kirchhoff stress tensor can be defined as

$$[eqn]$$

Then, the hydrostatic pressure can be obtained as

$$[eqn]$$

Under uniaxial tension, the general solution for the stress is derived as follows:

$$[eqn]$$

Planar Tension, PTThe specimen is subjected to tensile loading exclusively along the 1-direction. Therefore, under the planar tension mode, the deformation gradient $[eqn]$ and the right Cauchy–Green strain tensor $[eqn]$ are given by

$$[eqn]$$

The invariants of $[eqn]$ are defined as

$$[eqn]$$

Since the contraction in the 3-direction is free during planar tension, the first Piola–Kirchhoff stress matrix can be written as

$$[eqn]$$

It should be noted that, considering the particularity of the planar tension deformation mode, there is no deformation in the 2-direction (approximately), and the constraint of the fixture causes $[eqn]$ . Moreover, $[eqn]$ is difficult to obtain directly through experiments. The hydrostatic pressure can, therefore, be obtained as

$$[eqn]$$

Under planar tension, the general solution for the stress can be derived as

$$[eqn]$$

Equibiaxial Tension, ETThe specimen is subjected to tensile loading exclusively along the 1- and 2-directions, with equal stretch ratios. Therefore, under the equibiaxial tension mode, the deformation gradient $[eqn]$ and the right Cauchy–Green strain tensor $[eqn]$ can be obtained as follows:

$$[eqn]$$

The invariants of $[eqn]$ are defined as

$$[eqn]$$

Since the contraction in the 3-direction is free during equibiaxial tension, the first Piola–Kirchhoff stress matrix can be written as

$$[eqn]$$

Then, the hydrostatic pressure can be obtained as

$$[eqn]$$

Under equibiaxial tension, the general solution for the stress can be derived as

$$[eqn]$$

2.4.2. Specific Model

When developing a constitutive model, the key lies in determining its isochoric strain energy function $[eqn]$ , while the volumetric strain energy function $[eqn]$ can be determined by setting a large bulk modulus to define its contribution. According to research, in order to determine the nonlinear large-deformation characteristics of elastometric materials, $[eqn]$ needs to be a function of strain invariants as independent variables; to consider the biaxiality-ratio sensitivity, a second strain invariant must be introduced [27]. Therefore, we construct the following Carroll-like form [28] of the isochoric strain energy function

$$[eqn]$$

where a, b, and c are material parameters. Based on the previous derivations, the general solutions for the stress expressions of Equation (41) under the three different deformation modes—uniaxial, planar, and equibiaxial tension—are, respectively, given as

$$[eqn]$$

Moreover, for the purpose of our subsequent comparative analysis, we also introduce several classic models below—the Neo-Hookean [15], Mooney–Rivlin model [16,17], Yeoh [29], and Ogden models [30]—and modify them based on our nearly incompressible framework.

The isochoric strain energy function of the Neo-Hookean model is given by

$$[eqn]$$

where $[eqn]$ is the material shear modulus. Based on the previous derivations, the stress expressions under the three deformation modes can be derived as

$$[eqn]$$

The isochoric strain energy function of the Mooney–Rivlin model is given by

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are material parameters. Based on the previous derivations, the stress expressions under the three deformation modes can be derived as

$$[eqn]$$

The isochoric strain energy function of the Yeoh model is given by

$$[eqn]$$

where $[eqn]$ is a material parameter. Based on the previous derivations, the stress expressions under the three deformation modes can be derived as

$$[eqn]$$

The isochoric strain energy function of the Ogden model is given by

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are material parameters, N denotes the number of terms, and $[eqn]$ represents the equivalent principal stretch. Based on the derivation presented earlier, the stress expressions under the three deformation modes can be derived as

$$[eqn]$$

3. Experiments and Parameter Fitting

3.1. Experiments

The dimensions of the specimens are given in Figure 2a–c, each with a thickness of 1.5 mm. Ecoflex^TM^ silicone rubber was used in our experiments (00-20, 00-30, and 00-50); the material consists of two components, which were weighed to ensure a 1:1 mass ratio and manually mixed in a clean container until a homogeneous state was achieved. Subsequently, the mixture was poured into custom-made moulds (fabricated through resin polymer printing; see Figure 2d). Given that stirring introduces air bubbles, the mixture was then degassed in a vacuum chamber, typically for 5–10 min, until all visible bubbles were eliminated. Then, the specimens were left to cure at room temperature for at least 24 h to ensure complete crosslinking. A sample specimen is shown in Figure 2e.

All tests were conducted using a biaxial testing machine (see Figure 3, load cell with a capacity of 500 N). Each specimen was loaded at a quasi-static strain rate of 0.01/s under ambient laboratory temperature conditions. Note that the velocity dependence of the material was not taken into account in this study, as no significant disparity was observed in the stress values during pre-trial tests under different strain rates (0.001/s, 0.01/s, and 0.1/s). To prevent slippage during testing, aluminium tabs were bonded at the gripping locations. Four parallel trials were conducted for each test type, with specimens from the same batch used for each set of replicates to eliminate potential variations caused by the fabrication process. The nominal stress and strain of each specimen were calculated based on the recorded load-displacement history from the load cells. The testing machine was programmed to stop when specimen fracture occurs to thoroughly characterise failure behaviour. Prior to the test, speckled patterns were affixed to the surfaces of specimens to facilitate high-precision strain measurement using Digital Image Correlation (DIC). The DIC system utilised a camera sampling frequency of 1 Hz; the clamping photos are depicted in Figure 4.

The engineering stress–strain curves are shown in Figure 5, where error bands representing the standard deviation are shown for each curve. It can be seen that these error bands are relatively small, indicating a high degree of consistency in the test results. The higher the biaxiality ratio (UT → PT → ET), the higher the stress–strain curve of the material. This suggests that the mechanical behaviour of silicone rubber is sensitive to biaxiality. The results presented in Figure 5 fully reflect this sensitivity, which has often been overlooked in previous studies. From a micro perspective, during deformation, the continuous internal rotation of the main chain and end sliding occur within the polymer, causing the molecular chains to change from a coiled to an oriented state, macroscopically corresponding to a large degree of deformation. Moreover, this orientation behaviour leads to hardening behaviour under loading, particularly in the equibiaxial tension curve. However, the molecular chain orientation patterns of each specimen under different biaxiality ratios are not consistent. As the ratio increases, the number of molecular-chain orientation directions increases. Since the chain segments in each orientation are not independent of each other, they synergistically strengthen the molecular chain network to resist external deformation, manifesting macroscopically as a higher stress response under higher deformation and biaxiality ratios.

Under various biaxiality ratios, differences in specimen failure can also be observed. Table 1 summarises the failure strains under different deformation modes at the 00-30 Shore hardness. The results show that as the biaxiality ratio increases, the failure strain decreases, indicating a relationship between the two. As previously discussed with regard to the microscopic mechanism, an increase in biaxiality ratio leads to orientation occurring in more directions, increasing the probability of disentanglement in the molecular chain network during deformation. When large-scale chain separation or covalent bond breakage occur, the polymer chain network separates, corresponding macroscopically to the material’s failure behaviour.

3.2. Parameter Fitting

The parameters for each model are fitted simultaneously by using the test data from all three deformation modes without applying weighting. Through the application of the least-squares method, the objective of the fitting procedure is to minimise the difference between the experimental stress–strain data and the predictions from the constitutive models via parameter adjustment. To evaluate the goodness of fit for each model, the root mean square error ( $[eqn]$ ) is calculated, where $[eqn]$ represents the target values, $[eqn]$ represents the predicted values, and n is the total number of data points. The fitting results are presented in Table 2.

The fitting results in Figure 6 demonstrate that different constitutive models describes mechanical behaviour differently under different stress–biaxiality ratios:

1. The Neo-Hookean (Figure 6a) and Mooney–Rivlin models (Figure 6b) both fail to accurately describe the mechanical behaviour under different biaxiality ratios, especially under equibiaxial tension conditions. This is because neither considers the material’s biaxiality-ratio sensitivity, and the Neo-Hookean model relies only on the first strain invariant, limiting its ability to describe complex deformation behaviours.
2. The Yeoh (Figure 6c), Ogden (Figure 6d), and proposed models (Figure 6e–g) demonstrate relatively good fitting performance. Both the Yeoh model and the proposed model contain three parameters; however, the proposed model achieves better fitting results with lower root mean square error (RMSE) values. The Ogden model, despite its higher fitting precision (lower RMSE), requires twice the number of parameters compared to the proposed model. Overall, the proposed model achieves a good balance between parameter number and fitting accuracy, demonstrating better engineering practicality.

Based on our fitting results, we discuss further the strategy for developing the proposed model. The isochoric energy is dependent on the sequential expansion of free energy to minimise the errors remaining in the stress responses of previous terms. The first term is constructed to identify stress–strain data associated with uniaxial tension. In accordance with the Neo-Hookean model, the term $[eqn]$ is employed. However, according to the fitting results shown in Figure 6a, the Neo-Hookean model fails to capture the hardening behaviour of the material under uniaxial tension, especially in the large-deformation stage. Therefore, inspired by the Yeoh model in Figure 6c, a higher-order term $[eqn]$ is introduced to capture this hardening behaviour. Furthermore, the fitting results of the Neo-Hookean model and the Yeoh model indicate that merely incorporating $[eqn]$ is insufficient to describe the hardening behaviour of the ET curve. Thus, the second strain invariant term $[eqn]$ is introduced to capture this residual stress. This phenomenological strategy was initially proposed by Carroll [28], which we use to effectively determine the biaxiality-ratio sensitivity of elastomeric materials.

After parameter fitting, to further verify the accuracy of the proposed model, we conducted finite element simulations of uniaxial, planar, and equibiaxial tension tests with a mesh size of 1 mm. The specimen size, boundary conditions, and loading methods were consistent with our experiments, while the near-incompressible status of the material was reflected by setting a large bulk modulus of 1000 MPa. The deformation maps of our simulations and experiments are compared for the three modes in Figure 7a, Figure 7b, and Figure 7c, respectively, while the stress–strain comparison is shown in Figure 8. Overall, the constitutive model proposed in this paper can effectively describe the nonlinear large deformation and biaxiality-ratio sensitivity of silicone rubber materials. It should also be noted that no standardised stress calculation method currently exists for biaxial deformation, since it does not signify a uniform deformation condition. As such, our study integrates experiments, parameter fitting, and finite element simulations to iterate a correction factor for calculating biaxial stress at the specimen’s centre [31] in order to guarantee the accuracy of the calculation.

4. Conclusions

Silicone rubber materials have been widely applied in fields of engineering due to their excellent mechanical properties and chemical stability under complex deformation modes, which have attracted interest from both the academic and engineering communities.

In this study, we focused on developing a constitutive model for silicone rubber within the framework of continuum mechanics and an assumption of near incompressibility. Subsequently, uniaxial, planar, and biaxial tension experiments were systematically conducted at three Shore hardness levels; as the stress–biaxiality ratio and Shore hardness increased, the stress values of the material rose. Our experimental results exhibit high reproducibility (with relatively small error bands), indicating consistent and dependable experimental data. Our subsequent parameter calibration results indicate that this model can accurately and effectively describe the biaxiality-ratio sensitivity of silicone rubber materials (under the condition of 00-30 Shore hardness, the RMSE values for UT, PT, and BT were 0.0030 MPa, 0.0036 MPa, and 0.0013 MPa, respectively). Comparative analyses with four other classical hyperelastic constitutive models were also conducted, demonstrating the effectiveness and reliability of the proposed model. The Yeoh, Ogden, and proposed models all display good fitting performance, while the proposed model achieves a good balance between the number of parameters and fitting accuracy. After obtaining the model parameters, we performed finite element simulations of the uniaxial, planar, and biaxial tension tests. Our simulation results show good agreement with the experimental data in terms of deformation patterns and stress–strain curves. The maximum RMSE values between the simulated and fitted stress values were less than 0.0015 MPa across all deformation modes.

Overall, this research systematically integrates theoretical modelling, experimental testing, parameter fitting, and numerical simulation to explore the mechanical behaviour of silicone rubber material, especially its sensitivity to different deformation modes. Our approach and findings provide a reference for practical applications in sealing systems (e.g., optimising O-ring compression set performance under pressure loading) and soft robotics (e.g., enhancing the actuation precision of dielectric elastomer actuators). The limitations of this study are that the proposed model does not take into account the temperature effect, multiaxial fatigue behaviour, or anisotropic Mullins softening phenomena. These issues will be tackled in future research.