Role of Fiber Orientations and Invariant Coupling in the Stress–Strain Non-Coaxiality of Bioinspired Fiber-Reinforced Elastomers Open Access

Fiber-reinforced structures are ubiquitously present in nature and strongly influence the mechanical behavior of soft tissues and muscular hydrostats in multiple soft-bodied invertebrates such as octopuses, sea anemones, and worms [1,2]. Such structures can undergo large deformation, exhibit direction-dependent properties, and possess remarkable maneuverability and posture maintenance by leveraging the differential arrangement of muscle fibers that control such responses [3]. Such natural systems have inspired the development of bioinspired fiber-reinforced elastomers (FREs) that possess enhanced functionality and tunability. These systems hold great promise, particularly for soft robotics applications, allowing complex and programmable movements without the need for additional control systems [4,5]. FRE enclosures with asymmetric fiber orientations can control complex movements under pressure. These enclosures exhibit intricate structural changes, making them crucial design elements for soft actuators, grippers, and biomedical devices [6,7].

Bioinspired FREs possess nonlinear mechanical behavior, which depends on multiple factors such as the orientation of the fibers, individual mechanical properties of the fiber and matrix, type of applied loading and deformations, and the interaction between matrix and fibers [8,9]. Such responses are described by using a strain energy density function based on strain invariants in a hyperelastic framework. The number of invariants in the constitutive equation depends on the material symmetries, with the numbers varying from two for isotropic, five for transversely isotropic, and nine for orthotropic materials [10–13]. Formulation of constitutive models for FRE materials undergoing large deformations has been primarily motivated by the seminal works on rubber-like materials by Rivlin and Saunders [14]. These models were later adapted further to incorporate fiber orientation and reinforcement effects, as widely observed in biological tissues like the heart, arteries, and skin [15,16]. Additionally, Dong and Sun [17] have formulated a new hyperelastic model with planar distributed fibers that can incorporate the second type of Poisson effect as part of the material responses. Existing constitutive models for FRE materials consider invariant formulations and homogenization methods to explore failure instabilities and microstructure evolution during loading. However, the validation of these models, especially for FREs with controlled hierarchical microstructures, is currently undetermined by the lack of multiaxial experimental data, an important requirement for their validation. While most existing forms of the strain energy functions often employ an additive decomposition of the matrix and fiber contributions, the role of fiber–matrix interaction, often represented with invariant coupling terms, is crucial in the accurate mathematical modeling of such materials. The role of invariant coupling in understanding the mechanical behavior of anisotropic FREs is vital in the context of anisotropic hyperelastic modeling [16,18]. To effectively simulate the material response under large deformations, the constitutive relations of FREs must exhibit an invariant coupling to capture the anisotropic behavior of these materials in a physically consistent manner [10]. Research by Li and Zhang [19] highlighted how fiber orientation representing the matrix–fiber interaction affects stress–strain responses, integrating both isotropic and anisotropic properties in the models. Sun and Chen [20] underscored the need for comprehensive formulations to capture fiber–matrix interactions. Additionally, Castillo-Méndez and Ortiz-Prado [12] focused on using these invariants to predict failure modes.

Stress–strain coaxiality, a key characteristic of isotropic material class where the principal stress and strain directions align, is not observed in anisotropic materials like FREs, where these directions typically do not coincide. In FREs, the anisotropy introduced by the embedded fibers results in stress and strain tensors with distinct primary directions, which substantially impacts the mechanical behavior of the materials under deformation. In Beatty’s renowned paper [21], it is shown that stress–strain coaxiality is a fundamental mechanism for generating universal relations in materials. Universal relations are fundamental equations describing material-independent stress–strain behavior within a given material class. These universal relations for a given deformation are defined as homogeneous algebraic equations that (i) relate the components of the symmetric Cauchy stress tensor with the strain tensor and (ii) remain independent of the material response functions within the material class. Modeling the universal relations for the mechanical response of incompressible, transversely isotropic, bioinspired FRE has garnered much attention lately. Such biological soft elastomers reinforced with fiber bundles with a single preferred direction are abundant [22]. Researchers [21,23–26] have expanded on the universal relations, especially for isotropic, transversely isotropic, and FREs materials, using similar definitions, some of which apply only to isotropic materials. Beatty [21] highlighted the significance of universal relations for finite deformations in isotropic elastic materials, with the same principles extending to transversely isotropic materials. He derived a class of universal relations for constrained and unconstrained isotropic elastic materials. More recently, other researchers [23–26] have developed a general approach for identifying universal relations in continuum mechanics. However, the role of fiber orientations and invariant coupling in the stress and strain coaxiality/non-coaxiality-driven universal relations for FRE materials requires further exploration. Furthermore, Padovani and Šilhavỳ [27] examined the coaxiality of stress and strain in anisotropic “no-tension” materials. More recently, Melnic et al. [28] investigated the material responses of anisotropic materials by incorporating a coupling invariant term (⁠$I8$⁠) between two fiber families, although they did not derive the universal relations for this specific constitutive model. However, the role of fiber orientations and nonlinear invariant (⁠$I1$ and $I4$⁠) coupling in the stress and strain coaxiality/non-coaxiality-driven universal relations for FRE materials have not been studied in detail and requires further exploration.

In this study, we have investigated the role of invariant coupling in terms of $I1$ representing the matrix response and $I4$ corresponding to the fibers in the universal relations driven by stress–strain coaxiality of FREs with transverse isotropic symmetry. We adopted the form of strain energy function developed using experimental results from one of our earlier studies [16] and have extended it further to investigate how different forms of invariant coupling terms affect the stress–strain coaxiality-driven universal relations for the transversely isotropic constitutive models. Specifically, we calculate the variation in the axial vector terms concerning changes in the applied stretch, fiber orientation, and the degree of nonlinearity of both $I1$ and $I4$ invariants in the invariant coupling terms and compare the results for three distinct classes of deformations, simple shear, uniaxial, and nonequibiaxial tension. Our results show significant changes in the axial vector terms, including invariant coupling terms in the form of strain energy functions for all three deformation modes. Additionally, there are significant dependencies on the degree of nonlinearity of the invariant coupling terms. Together, these results provide some useful insights into the role of fiber–matrix interaction on the constitutive properties of FREs that can help in the development of efficient computational models to predict their behavior under complex loads and deformations.

This section presents an analytical modeling framework to develop the universal relations for incompressible and transversely isotropic nonlinear elastic solids. A second law of thermodynamics-based continuum mechanics approach is followed.

The universal relation (6) is valid for all deformation class $B$ and for all transversely isotropic material class of energy density $Ω(I1,I4)$⁠.

The derived universal relation (9) is valid for all transversely isotropic material classes of energy density $Ω(I1,I4)$ type and all deformation classes $B$⁠. While the relation holds for the chosen form of energy density $Ω(I1,I4)$⁠, its applicability to other strain energy functions $Ω(I1,I2,I3,I4,I5)$ depends on whether they share similar invariant dependencies.

The section addresses the results of the analytical modeling framework established in previous sections for developing universal relations of fiber-reinforced elastomers, with an emphasis on a majorly ignored invariant coupling under simple shear, uniaxial tension, and nonequibiaxial stretch. The selected deformation modes capture fundamental loading conditions in natural and artificial fiber-reinforced systems, allowing assessment of fiber–matrix interactions under different strain states. The selected deformation modes capture fundamental loading conditions in natural and fiber-reinforced systems like soft tissues [35], skeletal muscles [36], and plant cell walls [37] as well as artificial systems like soft actuators and grippers [38,39], allowing assessment of fiber–matrix interactions under different strain states and loading conditions. Illustrations of the effects of various invariant coupling terms on the coaxially driven universal relations between stress and strain are further shown for such material classes. The variations in axial vector terms are further reported for various fiber orientations, stretches, and the level of nonlinearity in the coupling terms with their coupling exponents. The selected fiber orientations (⁠$15deg,30deg,45deg,60deg,75deg$⁠, and $90deg$⁠) are based on multiple previous experimental studies [16,36,37]. These angles represent morphological configurations found in biological tissues such as skeletal muscles, tendons, and arterial walls and engineered elastomeric composites [38,39] in soft robotics and biomedical applications.

We used the material constants reported for the form of energy density function developed using experimental results from uniaxial tensile testing of FRE materials from our earlier study [16] to calculate the nonzero component of the axial vector $A3$ and $γ4$ (Eqs. (12) and (13)), for different applied shear stretches $γ$⁠. In order to examine the effect of invariant coupling terms on the variation of the coaxial term $A3$ and $γ4$ with the applied shear, we calculated these terms for two distinct cases. In the first case, we considered the form of the strain energy function (7) and calculated $γ4$ and $A3$ as a function of the material constants $(C3,C4,C5,α,β)$⁠, fiber orientation (⁠$θ$⁠), and applied shear (⁠$γ$⁠) as shown in Figs. 1(a), 1(c), and 1(e). We herein assumed $α=2$ and $β=1$⁠, representing the exact form of the strain energy function developed previously [16]. In the second case, we have dropped the $I1−I4$ invariant coupling term (associated with the material constant $C3$⁠) and recalculated the variation in $γ4$ and $A3$ for fiber orientations (⁠$θ$⁠) and applied shear (⁠$γ4$⁠) (Figs. 1(b), 1(d), and 1( f)) and compared the results from the two different cases. Our results show significant sensitivity to the inclusion of invariant coupling term in the strain energy function (7) for the variation in $A3$ and $γ4$ with applied shear for different fiber orientations. These differences are pronounced for $θ=45deg$⁠, $60deg$⁠, and $75deg$⁠, as shown in Figs. 1(a) and 1(b) for $γ4$ and Figs. 1(c) and 1(d) for $A3$⁠. Figures 1(e) and 1( f) show 3D surface plots that highlight the combined effect of variation in fiber orientations (⁠$θ$⁠) and applied shear (⁠$γ$⁠) on the coaxial vector term $A3$ with and without the inclusion of the invariant coupling term in the strain energy density (7).

As our calculations highlight that the inclusion of the $I1−I4$ invariant coupling term for multiple fiber orientations significantly influenced the variation in the coaxial vector term $A3$⁠, we wanted to investigate further the effect of nonlinearity in the invariant coupling terms by varying the coupling exponent terms $α$ (for $I1$⁠) and $β$ (for $I4$⁠) to changes in $A3$⁠. To this end, we varied both exponents $α$ and $β$⁠, from 1 to 5 each, and calculated $A3$ for a specific applied shear stretch (⁠$γ=1.2$⁠). Figures 2(a)–2( f) show the 3D surface plots for the variation in $A3$ due to the combined changes in both invariant exponents $α$ and $β$⁠. Interestingly, for all six different fiber orientation angles $θ$⁠, $A3$ shows sensitivity to changes in both $α$ and $β$⁠. These effects are most significant for $θ=45deg,60deg,75deg$⁠, and $90deg$ (Figs. 2(c)–2(f)). While the relative changes to $A3$ are less pronounced for $θ=15deg$ and $30deg$ (Figs. 2(a) and 2(b)), we still see some variations with changes in the exponents. The results suggest that the inclusion of invariant coupling plays a vital role in dictating the material responses for transversely isotropic FRE materials under simple shear.

Next, we used an approach similar to the one described in Sec. 4.1 to compute the coaxial vector term $A3$ and $γ4$ for different applied uniaxial stretches $λ$ (Eqs. (16) and (17)), using the material constants for the strain energy function for different fiber orientations, as obtained from Ref. [16]. We calculated the variation in coaxial term $A3$ and $γ4$ with the applied stretch $λ$⁠, with and without the $I1−I4$ invariant coupling term in the strain energy function, to investigate the effect of invariant coupling in terms of $I1$ and $I4$⁠. As mentioned earlier, in the first case, we assumed a value for the invariant coupling exponent terms, such that $α=2$ and $β=1$⁠, based on our choice of the form of strain energy function (7), then computed $γ4$ and $A3$ as a function of the material constants $(C3,C4,C5,α,β)$⁠, (⁠$θ$⁠), and applied stretch (⁠$λ$⁠) as shown in Figs. 3(a), 3(c), and 3(e). Variations in $A3$ and $γ4$⁠, in the absence of the invariant coupling term (associated with $C3$⁠) with changes in fiber orientations (⁠$θ$⁠) and applied stretch (⁠$λ$⁠), are represented in Figs. 3(b), 3(d), and 3( f). Similar to our previous results in the simple shear case (Fig. 1), we see pronounced differences in $γ4$ and $A3$⁠, with an increase in uniaxial stretch, for the $θ=45deg,60deg$⁠, and $75deg$ fiber orientations, as depicted in Figs. 1(a) and 1(b) for $γ4$ and Figs. 1(c) and 1(d) for $A3$⁠. The 3D surface plots also reveal the sensitivity of the coaxial vector term $A3$ to combined changes in fiber orientations (⁠$θ$⁠) and applied stretch (⁠$λ$⁠) as observed in Figs. 3(e) and 3( f).

Subsequently, we varied the coupling exponent terms $α$ (for $I1$⁠) and $β$ (for $I4$⁠) to investigate the effect of nonlinearity in the invariant coupling terms on the coaxial vector term $A3$⁠. Similarly, as reported in Sec. 4.1, we varied both $α$ and $β$ from 1 to 5 each and calculated $A3$ for a specific applied uniaxial stretch (⁠$λ=1.8$⁠) (Figs. 4(a)–4( f)). Specifically, we see that in the case of fiber orientation angle $θ=15deg$ and $45deg$⁠, the value of $A3$ increases monotonically with an increase in both exponents $α$ (for $I1$⁠) and $β$ (for $I4$⁠) (Figs. 4(a) and 4(c)). For other fiber orientations, such as $θ=15deg$⁠, $60deg$⁠, and $75deg$⁠, we see disparate responses, with the coaxial term $A3$ increasing nonlinearly with an increase in $α$ (⁠$I1$ exponent term) but shows an inverse dependence to changes in $β$ for both $30deg$ and $60deg$ orientations (Figs. 4(b) and 4(d)). There are no observable changes in $A3$ with changes in both $α$ and $β$ for the $90$ deg specimen (Fig. 4( f)), with the value of $A3$ remaining approximately zero throughout. These results provide useful insights into how invariant coupling influences the material responses of bioinspired FRE materials under uniaxial tension, depending on the orientation of fibers with respect to the direction of applied stretch.

Finally, we investigated the effect of invariant coupling on the variation of the coaxial vector term $A3$ and $γ4$ under nonequibiaxial stretch. Specifically, we varied the biaxiality ratio $η=λ2/λ1$⁠, where in $λ1$ and $λ2$ are the stretches in the $X$ and $Y$ directions (Figs. 5 and 6). We varied the value of $η$ from 0.5 to 1.5 at intervals of 0.25 and compared the effect of applied nonequibiaxial stretch to the variations in both $γ4$ (Fig. 5) and the coaxial vector term $A3$ (Fig. 6) with and without the invariant coupling terms, using the material constants for the strain energy function for all six different fiber orientations (⁠$θ$⁠) [16]. Specifically, Figs. 5(a)–5(d) show the variation in $γ4$ for four different values of $η$ (0.5, 0.75, 1.25, and 1.5) in the presence of the $I1−I4$ invariant coupling term in the strain energy function. In contrast, the variation in $γ4$ in the absence of the $I1−I4$ invariant coupling term is shown in Figs. 5(e)–5(h). Similarly, the variations in the coaxial vector term $A3$ subjected to different values of nonequibiaxial stretch (⁠$η=0.5$⁠, 0.75, 1, and 1.5) with considering the $I1−I4$ invariant coupling term in the strain energy function are shown in Figs. 6(a)–6(d) and without the invariant coupling term in Figs. 6(e)–6(h). Additionally, our theoretical investigations revealed that for the equibiaxial deformation case, i.e., $η$=1, both $A3$ and $γ4$ terms become zero, as can be deduced from Eq. (12) and have therefore been excluded from this section. For all these calculations, the invariant coupling exponent terms were fixed at $α=2$ and $β=1$⁠, as mentioned previously, using which $A3$ and $γ4$ were calculated as functions of the material constants (⁠$C3,C4,C5,α,β$⁠), fiber orientation (⁠$θ$⁠), nonequibiaxial ratio (⁠$η$⁠), and applied stretch (⁠$λ1$⁠).

Our results show that for increased values of $η$ (1.25 and 1.5), we observe that the inclusion of the invariant coupling term significantly influences the variation in $γ4$ with an increase in the stretch in the $X$-direction (⁠$λ1$⁠), especially for the fiber orientations $θ=15deg,60deg,$ and $75deg$ (Figs. 5(c), 5(d) and 5(g), 5(h)) for biaxiality ratios $η=1.25$ and 1.5, respectively. Additionally, the change in $γ4$ is significantly higher for all $θ$⁠, at higher values of $η$ (1.25 and 1.5) compared to 0.5 and 0.75. (Fig. 5). Similar trends were observed in the variation of values for the coaxial vector term $A3$ terms (Fig. 6). The variation of $A3$ with applied stretch (⁠$λ1$⁠) was significantly different at higher values of $η$ (1.25 and 1.5), specifically for the $θ=15deg,60deg$⁠, and $75reg$ samples (Figs. 6(c), 6(d) and 6(g), 6(h)). While some changes were also observed at values of $η<1$ (0.5 and 0.75), those changes were not very pronounced. Additionally, the changes in $A3$ were significantly higher in both cases when the value of $η$ was fixed at 1.5 (Fig. 6) compared to the other cases, an observation similar to the value of $γ4$ (Fig. 5). Together, these results highlight how the inclusion of invariant coupling terms can affect the universal relationship in transversely isotropic fiber-reinforced elastomers for different fiber orientations when subjected to nonequibiaxial stretch.

We have explored the role of $I1−I4$ invariant coupling in the stress–strain coaxially driven universal relations of transversely isotropic FRE materials under three different classes of deformations, namely, simple shear, uniaxial, and nonequibiaxial tension. Our results show that including invariant coupling terms significantly altered the coaxial vector terms for fiber orientations $θ=15deg$⁠, $60deg$⁠, and $75deg$ under simple shear deformations. Furthermore, these results were also sensitive to the effect of nonlinearity in invariant coupling terms, with changes in the coupling exponents $α$ and $β$⁠, significantly changing the variation in the coaxial vector term for all three orientations under shear. Similar to the case of simple shear, we also saw that the invariant coupling terms influenced the coaxial vector for the same set of fiber orientations $θ=15deg,60deg,$ and $75deg$⁠. However, the coaxial vector term variation with changes to $α$ and $β$ were different, with the changes being more pronounced for $θ=15deg$ and $45deg$⁠. No changes were observed in the coaxial vector term with respect to the invariant coupling term for $θ=90deg$⁠.

Finally, for the nonequibiaxial case, we again observed that the variation in the coaxial vector terms was regulated by including the invariant coupling terms when the biaxiality ratio $η$ was larger than 1. In all three cases, we observe that the effect of the fiber–matrix coupling term in the strain energy functions of the transversely isotropic FRE materials resulted in significant changes to the stress–strain coaxially driven universal relations for $θ=45deg$⁠, $60deg$⁠, and $75deg$⁠. Furthermore, the findings on $A3$ sensitivity to invariant coupling parameters reveal how fiber–matrix interactions affect stress–strain behavior in FREs. This is a key for optimizing actuator performance, enhancing deformation efficiency, and predicting failure in soft robotics. Bioinspired designs like muscle-mimicking actuators and biomedical implants can benefit by tailoring fiber alignment for desired responses. Such findings provide increased evidence of the importance of fiber–matrix coupling terms in the constitutive properties of FRE materials. Such coupling terms are crucial in developing sophisticated computational tools to predict their behavior under different scenarios with potential applications in soft robotics and biomedical applications. The theoretical framework employed in this study can also be used to investigate the effect of higher order invariants like $I5$ or materials with two families of fibers having orthotropic symmetries to the stress–strain coaxiality responses of FRE materials, which can become a future scope of this work. Another scope of the future work is to extend the experimental validation of invariant coupling effects in fiber-reinforced elastomers to nonequibiaxial tensile tests and to study the inflation–extension responses of hollow cylindrical tubes made up of FRE materials with transverse isotropic symmetry that is particularly relevant in soft robotics and biomedical material-based applications.

AC acknowledges the support from NFSG (New faculty seed grant) at BITS Pilani, Hyderabad.

The authors declare that there are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.