On tube models of rubber elasticity: fitting performance in relation to sensitivity to the invariant I2

This study investigates the capability of micromechanical models of rubber elasticity to predict the deformation behaviour of soft materials under various modes of deformation. The free energy of individual chains is decomposed into a freely fluctuating chain contribution and a tube contribution representing topological constraints. Full-network averaging over all chain orientations is considered, along with three-chain and eight-chain approximations. The performance of various tube formulations is analysed in relation to their implicit (or in some cases explicit) dependence on the second invariant I2. We show that micromechanical models that involve the area-stretch of the macroscale continuum lead to I2 dependence when combined with the full-network averaging scheme, whereas micromechanical models that only involve the line-stretch of the continuum show much weaker I2 sensitivity. However, I2 sensitivity can emerge from line-stretch–based micromechanical models when three-chain averaging is used. Comparisons between model predictions and experimental data confirm the direct correlation between strong I2 sensitivity and fitting performance. Overall, our study suggests that micromechanical models of rubber elasticity should involve both the line-stretch and the area-stretch to elicit I2-dependent behaviour and reproduce experimental trends.


Introduction
Models of rubber elasticity can be broadly classified into two main groups: phenomenological models, and micromechanical models.Phenomenological models are obtained by constructing a free energy function written in terms of principal stretches or invariants in order to fit experimental stress-strain curves.In contrast, micromechanical models are obtained by averaging the response of representative chains sampling the orientation space.In principle, micromechanical models can provide insights into structure-property relationships, and thus serve as a guide for the design of rubber-like materials.They also provide a framework to investigate phenomena where directional effects are significant, such as damage-induced anisotropy [13] or strain-induced crystallisation [10].However, they often come at a higher computational cost, while their ability to outperform phenomenological models remains unclear, at least for a purely hyperelastic behaviour [8].The purpose of this paper is to shed some light into the performance of micromechanical models through an analysis of their sensitivity to the second invariant I 2 of the right and left Cauchy-Green deformation tensors.
The formulation of micromechanical models of rubber elasticity usually involves the following steps [5,16,34]: 3. The averaging of the chain energy over representative chain orientations to obtain the macroscopic free energy function.
The free energy of individual chains can be estimated from non-Gaussian statistical mechanics considerations for a freely jointed phantom chain fluctuating between two crosslinking points.The relevant micro-kinematic variable is then the chain stretch λ, which under the assumption of affine motion of the crosslinks, is identified with the line-stretch of the continuum in the chain direction.Averaging the chain energy over a uniform distribution of chain orientations, one obtains the (affine) full-network model [41,44].By replacing full-network averaging with directional sampling in the three principal stretch directions, the threechain model is recovered [24].The Arruda-Boyce eight-chain model [2] is obtained by only considering chains aligned with the diagonals of a cube with edges in the principal stretch directions.The average-stretch full-network model, proposed by Beatty [3], is another simplification of the full-network model, in which the macroscopic energy is obtained from the chain free energy evaluated for the average of the (squared) stretch over all orientations.The final constitutive equations are identical to those of the eight-chain model.In the limit of small stretches and Gaussian statistics, all these models reduce to the neo-Hookean model.
Affine full-network models and their reduced orientation sampling approximations generally show limited capability in reproducing different loading conditions with a single set of material parameters [5,8,22,39].Improving the fitting performance of micromechanical models can however be achieved by including a contribution from topological constraints to the chain free energy.The most popular approach consists in assuming that chain fluctuations are restricted to a confining tube [14].Existing tube formulations differ in their statistical mechanics treatments and in the assumed relationship between the tube dimensions and the macroscopic deformation.Most tube models express the change in tube diameter in terms of the chain stretch λ.Examples include affine tube models [15,19], constant volume tube models [20,33], and non-affine tube models [21,38].A general presentation, which encompasses the previously mentioned models as special cases, has recently been proposed by Darabi and Itskov [9], in combination with three-chain averaging.We refer to tube models that relate the change in tube dimensions to the chain stretch, and hence to the line-stretch of the continuum, as "λ-based" tube models.In contrast, Miehe et al. [34] introduced the tube area contraction ν as a second independent micro-kinematic variable, and further proposed a (non-affine) relationship between the tube area contraction and the area-stretch of the continuum.These authors implemented their tube model within a full-network averaging scheme.This approach was recently simplified by considering eight-chain averaging (or, equivalently, area-stretch averaging in the sense of Beatty [3]); see [1,29,30].We call tube models that relate the tube area contraction to the area-stretch of the continuum as "ν-based" tube models.
In the context of phenomenological approaches, it is well established that comprehensive models of rubber elasticity must depend on both I 1 and I 2 to capture different loading modes [1,8,12].Therefore, the question arises as to whether the improved performance of (some) tube formulations is linked to the introduction of stronger I 2 sensitivity.This is very clear in the case of ν-based eight-chain (or area-stretch average) tube models, where the I 2 dependence appears explicitly [1,29,30].However, for λ-based tube formulations and other orientation averaging schemes, the emergence of the second invariant I 2 in relation to fitting performance remains unclear.Previous comparisons of tube models [8,22,32] with relation to fitting performance do not focus on the systematic differences that arise from different kinematic assumptions, which we find to be central to the problem.Moreover, the difference in performance shown by a given tube model under the three-chain, eight-chain, and full-network schemes has not been systematically investigated before.Finally, the connection between fitting performance on multiaxial data and I 2 dependence has not always been explicitly drawn in previous studies [8,22,32].
This contribution aims to systematically assess the capability of tube models of rubber elasticity to elicit I 2 dependence, and further to relate their I 2 sensitivity to their ability to fit experimental data under different loading modes.To focus on the effect of the tube contribution on the mechanical response, we adopt the freely jointed chain model for the phantom contribution under the affine stretch assumption in all cases.Tube models are assessed using full-network averaging, as well as three-chain and eight-chain approximations.Our results show that λ-based tube models are only weakly sensitive to I 2 in their full-network implementation; however, stronger I 2 dependence can emerge in the three-chain averaging approximation.In contrast, ν-based tube models show strong I 2 sensitivity under all orientation averaging schemes.Our results further confirm that I 2 sensitivity is directly correlated to fitting performance.

Chain free energy
We describe the behaviour of rubber-like materials by considering a distribution of individual chains.The end points of the chains are fixed at crosslinking points, and we write r and r 0 the distances between crosslinking points in the current and reference configurations, respectively.Chain fluctuations are limited to a tube, representing topological constraints due to other chains.We write d and d 0 the diameters of the tube in the current and reference configurations.Following Miehe et al. [34], we introduce two micro-kinematic variables: the chain stretch λ and the tube area contraction ν, defined as: We assume that the free energy ψ of individual chains in the network can be additively decomposed into two contributions: where ψ f represents the free energy of the phantom chain (i.e.free of topological constraints) and ψ c accounts for topological constraints due to the surrounding chains.

Phantom contribution
The free energy of the phantom chain ψ f is obtained from the entropy of conformation for a freely jointed chain consisting of N Kuhn segments of length b; see e.g.[5].The free energy of the freely jointed phantom chain is expressed as: where k B is the Boltzmann constant, T is the absolute temperature, and L −1 (x) is the inverse Langevin function.Here, we use the Padé approximant [7]: Making the assumption that crosslinking points move affinely with the macroscale continuum, the chain stretch is given by: where denotes the line-stretch of the continuum: For a chain aligned with direction n ( | | = 1 ) in the reference configura- tion, the line-stretch is calculated from the macroscopic deformation gradient F as: We note that non-affine localisation rules have also been proposed as an alternative to Eq. ( 5); e.g.[18,34,40].However, these are not considered in this work.
Under the affine assumption (5), and taking r 0 as the most probable random walk end-to-end distance r 0 = √ Nb the free energy (3) of the freely jointed phantom chain can be rewritten as: Equation (7) shows that maximum extensibility of the chain is reached when → √ N , at which point the energy diverges.In the limit of small chain extension, the free energy (7) recovers the expression obtained under the assumption of Gaussian statistics: As a rule of thumb, the Gaussian approximation holds for r ≤ 0.3Nb , or, using r 0 = √ Nb and = , ≤ 0.3 √ N.

Tube contribution
The free energy contribution associated with the reduction in conformational entropy arising from the tube constraint is taken of the following form [14,34]: where α is a geometry factor pertaining to the shape of the tube cross section.
λ-based models Most existing tube models relate the tube area contraction ν to the chain stretch λ, so that the two microkinematic variables are effectively not independent.Under the affine stretch assumption (5), here we consider a model of the general form: The kinematic relation (10) together with the tube free energy expression (9) gives the tube free energy as: where k and β are here treated as fitting parameters.The superscript (λ) has been added to clarify the functional dependence of the tube free energy on the line-stretch only.Expression (11) coincides with the non-affine tube model proposed by Heinrich and coworkers [21].The early affine tube model, which assumes d = d 0 [15, 19], is recovered taking = −2 .The constant volume tube model [20, 33] corresponds to = 1 .Recently, Darabi and Itskov proposed a generalised tube model with free energy of the form: c = k + , with as an additional fitting parameter.The Darabi-Itskov model is obtained as the superposition of two tube contributions of the form (11) with different powers.Various -based tube models of the literature are summarised in Table 1.
ν-based models Miehe et al. [34] proposed to relate the tube area contraction ν to the area-stretch v of the continuum by a non-affine relation: The area-stretch is obtained from the macroscopic deformation gradient F and the chain direction n using Nanson's formula: where J = det( ) .The area-stretch v represents the ratio of deformed area to reference area of a surface element with out- ward unit normal n in the reference configuration.Substituting Eq. ( 12) into expression (9) of the tube free energy, the tube contribution takes the general form: where k and β are fitting parameters and the superscript (ν) indicates the functional dependence of the free energy on the area-stretch.Affine tube [15,19] k Constant volume tube [20,33] k Heinrich et al. [21] k Rubinstein and Panyukov [38] k( + −1

Macroscopic free energy
The macroscopic free energy Ψ of the network is obtained by averaging the chain free energy over representative orientations.We consider three averaging schemes, respectively corresponding to the full-network, three-chain, and eight-chain models.

Full-network averaging
The full-network (FN) approach to rubber elasticity considers a uniform distribution of chain orientations in the reference configuration.The macroscopic free energy is given by: where n is the number density of polymer chains and integration is carried out over the unit sphere S .Using Eq. ( 2), the macroscopic free energy can be additively split into its phantom and tube contributions: where Ψ FN f and Ψ FN c respectively represent the phantom and tube contributions to the macroscopic free energy.The macroscopic phantom contribution is independent of the the specific tube model being considered.The macroscopic tube contribution differs depending on whether a -based or -based tube model is used.For a -based model: where ( )  c is a function of the line-stretch and is given by Eq. ( 11), while for a -based tube model: where ( ) c is a function of the area-stretch and is given by Eq. ( 14).Evaluation of the integral over all the orientations generally requires numerical quadrature.Here, we adopt the 900-point scheme proposed in Ref [17] to minimise integration errors; see also [23,42].

Three-chain averaging
The macroscopic free energy of the three-chain (3C) model is obtained by restricting the orientation averaging to chains that are initially aligned with the principal stretch directions.Let λ i (i = 1, 2, 3) be the principal stretches, and ν i be the area-stretches in the principal directions.According to Eq. ( 13), the latter are given by: ν 1 = λ 2 λ 3 , ν 2 = λ 1 λ 3, and ν 3 = λ 1 λ 2 .The macroscopic free energy of the three-chain model is given by: where the superscript (i) refers to a chain aligned with the i th principal direction.More explicitly, the macroscopic phantom contribution to the free energy is given by: The macroscopic tube contributions for the λ-based and ν-based tube models are respectively given by: and: (15) Eight-chain averaging The eight-chain (8C) model only considers chains aligned with the diagonals of a cube with edges aligned with the principal directions.In a coordinate system aligned with the principal axes, = (±1∕ 3) .Direct application of formulas ( 6) and ( 13) in each of the eight-chain directions gives: and: where the superscript (i) refers to the i th chain ( i = (1, 8)) .Equations ( 23) and (24) show that all the representative chains experience the same line-stretch and the same area-stretch; and therefore, all representative chains have the same energy.Averaging of the chain energy over the eight directions then directly gives: The macroscopic phantom contribution is given by: The macroscopic tube contributions in the λ-based and ν-based schemes are respectively given by: and

Macroscopic stress-strain response
In this work we assume that the material is incompressible (J = 1).The Cauchy stress is given by: where p is an arbitrary scalar to be determined from boundary conditions.The principal stresses are obtained from Eq. ( 29) as: The shear modulus G obtained with the different models can be calculated from the stress-stretch response in a uniaxial tensile test along direction 1 as: The shear modulus can be additively decomposed into a phantom contribution G f and a tube contribution G c .Assuming that the Gaussian limit holds in the reference configuration, which is the case if N ⪆ 11, then the phantom modulus is simply given by G f = nk B T, which holds for the full-network, three-chain, and eight-chain averaging methods.(For a smaller number of Kuhn segments N, chains are pre-stretched beyond their Gaussian regime in the reference configuration, so that finite (23 extensibility effects must be included in the expression of the modulus.)In contrast, the tube contribution to the shear modulus depends on the chosen averaging method.For the full-network, three-chain, and eight-chain models, we respectively find: Note that identical expressions are found for the -based and -based tube models.When the tube energy is quadratic, i.e. = 2 , all three moduli coincide: G c = 2nk 3 .

Evaluation of I 2 dependence
In this work, we aim to relate the fitting performance of various tube model formulations to their dependence on the second invariant I 2 of the left Cauchy-Green deformation tensor = T .Under the assumption of incompressibility, the invari- ants are given by: or, in terms of the principal stretches: For future reference, we also recall the following expression of the Cauchy stress for an incompressible material: and the principal stresses follow as:

Direct calculation of @ @I 2
The most direct way to evaluate the I 2 dependence of a hyperelastic model or experimental data is to calculate the partial derivative ∂Ψ/∂I 2 .For an incompressible, hyperelastic material subjected to plane stress deformation (σ 3 = 0), the partial derivatives of the free energy with respect to the invariants can be expressed in terms of the principal stretches and stresses as [1,37]: Note that these expressions are only valid when all three principal stretches are distinct; i.e. they do not apply to uniaxial tension or equibiaxial tension.Furthermore, their evaluation is problematic in the small deformation regime where all principal stretches approach 1.The stresses 1 and 2 must approach 0 sufficiently rapidly for a finite value of Ψ∕ I 1 and Ψ∕ I 2 to be reached.Numerical or experimental noise can interfere with this limit and renders it challenging to make statements about the I 2 dependence with confidence.For this reason, Rivlin and Saunders in their seminal work do not apply Eqs. ( 39)-( 40) for deformations where I 1 , I 2 < 5.

Universal relations
To circumvent the aforementioned difficulties associated with the direct calculation of the partial derivatives ( 39)-( 40), we also consider universal relations to evaluate the emergence of I 2 dependence in micromechanical models and experimental data.Universal relationships are relationships between principal stresses and stretches that hold true irrespective of the specific form of the free energy function belonging to a certain class.If the stress-stretch response predicted by a specific model violates a universal relationship associated with the class of material behaviour, then the model does not belong to that class.Several universal relations have been proposed for the class of free energy functions defined as Ψ = Ψ(I 1 ) (i.e.where I 2 dependence of the free energy function is excluded) [1,43].In the following, we recall an existing universal relation valid for general biaxial loading (i.e. when all three principal stretches are distinct), and further propose two new universal relations respectively applicable to uniaxial and equibiaxial loading, and to uniaxial tension and compression.

General biaxial tension
Consider first a general biaxial loading with σ 3 = 0. Then the general stress-strain relation (38) gives: If the free energy only depends on I 1 , then the following universal relation for general biaxial test should hold [1,43]: This universal relation is applicable when the three principal stretches are distinct (general biaxial loading).This somewhat limits its applicability, as experimental data for general biaxial loading are not commonly reported in the literature.Therefore, here, we formulate two alternative universal relations which respectively apply to uniaxial and equibiaxial loading, and to uniaxial tension-compression loading.Uniaxial and equibiaxial tension Consider uniaxial (UT) and equibiaxial (EBT) tensile tests with σ 3 = 0.In these two loading cases, the non-zero stresses are respectively given by: If the free energy depends only on I 1 , and if λ U T and λ EBT are evaluated at the same value of I 1 , then the following relation between principal stresses and stretches should hold: The condition that the stretches in the two loading cases correspond to the same I 1 requires that:

Factorising the above relation gives:
Therefore, the tensile stretches in the uniaxial and equibiaxial tension tests which give the same I 1 are related by: (The other admissible root UT = −2 EBT corresponds to a uniaxial compression test and is not considered.)The universal relation (47) can be used to assess I 2 dependence from uniaxal and equibiaxial tensile data in the following way.For each tensile stretch value in the equibiaxial test, calculate the stretch in the uniaxial tension test which gives the same I 1 using relation (49).Next evaluate condition (46).If it is not satisfied, it implies that the free energy function depends on I 2 .Note that the converse is generally not true: satisfaction of condition (46) at a given I 1 does not necessarily guarantee that Ψ is independent of I 2 .

Uniaxial tension and compression Now consider the non-zero stresses in uniaxial tension (UT) and compression (UC) tests:
Evaluating these expressions at the same value of I 1 , and assuming that the free energy is independent of I 2 , we obtain: The axial stretches in uniaxial tension and compression tests which give the same value of I 1 are related by: which is readily obtained from Eq. ( 49) by noting the correspondence between the stretches in a tensile equibiaxial test and in a compressive uniaxial test: UC = −2

EBT
The universal relation (52) can be verified graphically using the Mooney stress, which is defined as: The universal relation (52) means that the Mooney stress M , when represented as a function of I 1 , should be the same in uniaxial tension and uniaxial compression for the material response to be independent of I 2 .

Results
In this section, we use the various tube formulations described in Sect. 2 to reproduce experimental stress-strain curves reported in the literature for several representative rubber-like materials.We also apply the universal relations described in Sect.3.2 to assess I 2 dependence in theoretical and experimental stress-strain curves. (47)

Data of Kawabata et al. (1981)
We first consider the experimental data of Kawabata et al. [26] for an isoprene rubber vulcanisate.The loading conditions correspond to general biaxial tension with σ 3 = 0, where λ 1 is held constant and λ 2 is varied, covering different modes of deformation from uniaxial tension to equibiaxial tension.We fitted λ-based and ν-based micromechanical models using the experimental data for uniaxial and equibiaxial tension, and subsequently kept the parameters constant for predicting the response under general biaxial loading.Fitting was carried out using least squares minimisation as implemented in the Python Scipy optimisation toolbox.The fitting error to be minimised is defined as: where N p represents the number of data points, (i)  e the i th experimental stress data point, and (i) m the corresponding stress calculated using the model.Best-fit parameters and the minimised fitting error for each model are reported in Table 2.
Figure 1 shows the full-network model predictions considering λ-based (Fig. 1a) and ν-based (Fig. 1b) tube models, along with experimental data for uniaxial and equibiaxial tension.Both models predict a stiffer response under equibiaxial tension than under uniaxial tension.However, the λ-based full-network fails to capture the response under equibiaxial tension accurately.In contrast, the ν-based full-network tube model provides an excellent fit.(Note that modelling and experimental curves for equibiaxial tension practically coincide.)Model predictions for general biaxial loading using the same set of parameters are compared to experimental data in Fig. 2. In the figures, each line corresponds to a constant value of the principal stretch λ 1 , as indicated.The λ-based model predictions (Fig. 2a and b) deviate significantly from the experimental data, particularly at the highest values of λ 1 .For a given value of λ 1 , the error is the smallest under uniaxial states of deformation, and increases monotonically with λ 2 .In contrast, the ν-based model captures the general biaxial response well for all modes of deformation (55) 2  Alternatively, we may also consider the universal relation (43), which applies to general biaxial tension with distinct principal stretches.Similarly to the previous case, we define an indicator function: so that G 2 = 0 corresponds to an I 2 independent response.The indicator functions calculated based on the experimental data for general biaxial loading (neglecting data points with repeated principal stretches) and model predictions for the λ-based and ν-based models are shown in Fig. 4. Again, the experimental data show I 2 dependence ( G 2 is non-zero).The ν-based model (Fig. 4b) is able to recover the experimental I 2 dependence at large values of λ 1 and λ 2 , where the λ-based model shows little sign of I 2 dependence (Fig. 4a).For smaller values of λ 1 , both models fail to recover the experimental trend.However, in this regime, the indicator function is extremely sensitive to measurement errors as both denominators in Eq. ( 57) tend to 0.
We next consider the three-chain and eight-chain averaging schemes as a simplification of the full-network model.Model parameters obtained by fitting the experimental data for uniaxial and equibiaxial tension are shown in Table 2. Comparisons between experimental data and model fits are shown in Fig. 5a  (57) Fig. 3 The I 2 indicator function G 1 (I 1 ) for uniaxial and equibiaxial tension, Eq. ( 56), applied to the experimental data of Kawabata et al. [26] and to model predictions using full-network λ-based and ν-based tube models Fig. 4 The I 2 -indicator function G 2 ( 1 , 2 ) for general biaxial loading, Eq. ( 57), applied to the experimental data of Kawabata et al. [26] and model predictions using full-network a λ-based and b ν-based tube models

Data of Dechnarong et al. (2021)
We next consider experimental data reported by Dechnarong et al. [11] for a thermoplastic elastomer (SEBS) subjected to uniaxial and equibiaxial tension.A similar fitting procedure was adopted as in the previous section, and best-fit parameters for the various tube models are reported in Table 3.

Data of Bechir et al. (2006)
Finally, we consider the tension-compression data reported by Bechir et al. [4] for a carbon black-filled natural rubber.In the experiment, the material has been pre-conditioned to render it effectively hyperelastic and eliminates the Mullins effect.The data are shown in Fig. 7 with model predictions obtained with the λ-and ν-based models under different averaging schemes.Corresponding best-fit parameters are shown in Table 4.
All models are able to capture the upturns in the stress-stretch curve under tension and compression and achieve broadly similar fits in the moderate stretch region except for the λ-based eight-chain model, which fails to capture the stress upturn in compression.Close inspection of the phantom and tube contributions for each model reveals that the stress upturn in tension arises from the extensibility limit in the freely jointed chain model (as expected), whereas the stress upturn in compression is due to the tube contribution, which in the λ-based eight-chain model is unable to capture the experimental trend.I 2 dependence is investigated by plotting the Mooney stress (54) against I 1 ; see Fig. 8a and b.Experimental curves show a different response under uniaxial tension and compression, which indicates I 2 dependence according to the universal relation (52).This is well predicted by the λ-based full-network and three-chain models.However, the λ-based eight-chain model shows the same Mooney stress under uniaxial tension and compression, implying that this model is I 2 independent.In contrast, all three ν-based formulations predict significant I 2 dependence in general agreement with the experimental data.
Equivalently, we may define an additional indicator function measuring deviation from the universal relation ( 52): such that I 2 independence implies G 3 = 0 .Figure 8c and d show this indicator for the data of Bechir et al. [4].As expected, Fig. 8c and d convey the same information as Fig. 8a and b: all models are able to recover the experimental I 2 dependence, except the λ-based eight-chain model, which is independent of I 2 .Here, we included both the Mooney-I 1 plot and the indicator function plot for the sake of illustration, as to the best of our knowledge neither of these tests had previously been used to assess I 2 dependence from uniaxial tension and compression data.

Discussion
Main findings from the comparison of model predictions to experimental data can be summarised as follows: 1. ν-Based tube models are able to provide a good fit for a range of experimental data under various loading conditions, regardless of the orientation-averaging scheme used (full-network, three-chain, or eight-chain averaging).2. λ-Based tube models can only provide a good fit under various loading modes when used in conjunction with three-chain averaging.3. Fitting ability of all the models is directly correlated to their I 2 sensitivity, as measured via their departure from universal relations valid for I 2 -independent behaviour.
To further explore the I 2 sensitivity of various tube formulations, we calculated the partial derivatives of the tube contribution to the free energy, ∂Ψ c /∂I 1 and ∂Ψ c /∂I 2 , as described in Sect.as function of the exponent β a state of pure shear with λ 1 = 2, λ 2 = 0.5, and λ 3 = 1 in Fig. 9. Pure shear was considered because the three principal stretches are distinct and formulas ( 39) and ( 40) can be used.The range of exponent values was chosen to be representative of best-fit exponent values obtained in the previous section.For the λ-based full-network model (Fig. 9a and b), both ∂Ψ c /∂I 1 and ∂Ψ c /∂I 2 are small relative to (nk) and have a similar order of magnitude across the considered range of β-values.The λ-based eight-chain model shows zero I 2 dependence.In contrast, the λ-based three-chain model shows significant I 2 dependence at large negative values of β.Regarding ν-based models, all three averaging schemes lead to significant ∂Ψ c /∂I 2 , which increases as β increases (Fig. 9d).The three-chain and full-network models also show significant negative ∂Ψ c /∂I 1 at large positive values of β, while the tube model is independent of I 1 when used in combination with eight-chain averaging (Fig. 9c).These observations for a broader range of β-values are consistent with results of Sect. 4. The I 2 (resp.I 1 ) independence of the -based (resp.-based) eight-chain model is easily understood.Using Expression (36) of the invariants in terms of the principal stretches, Eq. ( 23) directly gives the line-stretch of the representative chains in the -based tube model as 8C = √ I 1 ∕3 , and Eq. ( 24) gives the area-stretch of the representative chains in the -based tube model as 8C = √ I 2 ∕3 .Therefore, Eq. ( 25) provides, for the -based tube model: and for the -based tube model: Thus, λ-based tube models only bring I 1 dependence in the eight-chain averaging approximation, while ν-based tube models only bring I 2 dependence in this averaging scheme.Note that identical expressions are obtained by adopting the stretch-average full-network approximation proposed by Beatty [3].In this approach, the macroscopic tube energy is obtained by evaluating the λ-based tube energy using the full-network average of the squared line-stretch as argument: and, for the ν-based tube energy, using the full-network average of the squared area-stretch as argument: Expressions (61) and (62) are respectively identical to Expressions (59) and (60), recalling the identities [27] (see also Appendix A.2 and A.3): The I 2 dependence of ν-based eight-chain tube formulations has previously been pointed out by several authors [1,16,29,30].Setting β = 2 in Eq. ( 60), one recovers the C 2 -term of the Mooney-Rivlin model, with C 2 = nk/3.For β = 1, one recovers the phenomenological Carroll model [6].In general, the model ( 60) is of the Rivlin-Saunders type [37].
Regarding the three-chain models, one should first recognise that λ-based and ν-based tube formulations are in fact mathematically equivalent.Their contributions to the macroscopic free energy are respectively given by: and: where we have used the condition of incompressibility to express the principal area-stretches in Eq. ( 19) as: −1 , and ν 3 = λ 1 λ 2 = λ 3 −1 .Expressions (64) and (65) are identical provided that β ′ = − β (and that (nk) = (nk) ′ ).
(61) This is consistent with expression (33) the shear modulus, which involves β 2 .In light of this, the same level of performance exhibited by the two three-chain formulations was expected.The best-fit parameters reported in Tables 2, 3, and 4 indeed show the expected correspondence between λ-based and ν-based tube parameters.The expected relations β ′ = − β and (nk) = (nk) ′ are however not exactly satisfied in Tables 2 and 3, which is attributed to the fact that best-fit parameters are non-unique [36].This is confirmed by the fact that the fitting error associated with λ-based and ν-based models is the same in these two examples.Using the set of parameters identified using the λ-based model to predict results with the ν-based model (upon switching the sign of β) does predict an identical response for all datasets, as expected.The mathematical equivalence between λ-based and ν-based three-chain tube models does not however explain their stronger I 2 sensitivity, as compared to their full-network counterpart.We note however that the affine tube model, obtained taking β = − 2 in the λ-based model (i.e.β = 2 in the ν-based model), leads again to the C 2 term of the Mooney-Rivlin model with C 2 = nk/3.It is also worth pointing out the formal similarity between the three-chain tube model and the one-term Ogden model [35], previously noted in [25].Some insights into the I 2 sensitivity of λ-based and ν-based full-network models can be gained by considering Taylor series expansions of the macroscopic tube energy.Detailed derivations are provided in Appendix A.1, along with a numerical verification for small stretches.The λ-based macroscopic tube energy admits the following series expansion about the full-network average of the squared line-stretch, ⟨ only depends on I 1 and coincides with the eightchain model (59).The second invariant I 2 appears in the second-order term, but is coupled with I 1 .Similarly, the series expansion of the ν-based full-network tube energy about the average squared area-stretch ⟨v 2 ⟩ is given by: where we have used the notations Γ ∶= v 2 and w (v) c (Γ) ∶= (v)  c (v) .This expression shows that, to the lowest order, the tube contribution to the macroscopic energy only depends on I 2 and is identical to the eight-chain expression (60).The comparison between Eqs. (66) and (67) suggests that the stronger I 2 sensitivity in the ν-based full-network model arises from its lowest order term.
Similar series expansions can be constructed for the three-chain models, as detailed in Appendix A.1.For the λ-based model, the tube energy can be expressed as: and for the -based three-chain model: Expressions (68) and (69) only differ from their full-network counterparts (66) and (67) by a numerical factor in the second-order term (and by higher order terms).Invoking the same argument as in the full-network case, one would conclude that the three-chain λ-based model is dominated by I 1 , and the three-chain ν-based model is dominated by I 2 .However, as demonstrated above, the two models are formally equivalent.This highlights the importance of higher order terms on the role of invariants in the free energy function.
In addition to fitting performance, the ν-based tube model is also superior in that its best-fit parameters (nk) and β obtained from a given experimental dataset are fairly consistent across all three integration schemes (under the caveat of non-uniqueness of best-fit parameters).This is not the case for the λ-based models, where best-fit parameters show a much larger variation across the different integration schemes for the same experimental dataset; see Tables 2, 3, and 4. Given the molecular basis of the tube model, one would hope for a low sensitivity of the parameters on the orientation averaging scheme used for the parameters to have a physical meaning, considering the three-chain and eight-chain averaging schemes as (66) cost-effective simplifications of the full-network averaging.ν-Based models fulfill this requirement some extent, whereas λ-based models clearly do not.
Best-fit parameters predicted by the ν-based models can also be compared to theoretical values derived from statistical mechanics considerations.The non-affine tube model of the form ( ) c = k B was derived by Heinrich and coworkers [21,25], where three-chain averaging was used.These authors expressed the kinematic relation for the tube deformation as d = d 0 � , where α = 1/2 and β ′ is an empirical fit parameter in the range 0 < β ′ ≤ 1.The original non-affine model is recovered from our formulation by setting β = −2αβ ′ = − β ′ .According to these authors, a value β ′ = 1 (β = −1) is relevant for unswollen, well-connected networks [25].The actual exponent value may however differ due to various constraint release mechanisms during deformation, or due to the presence of rigid fillers leading to strain amplification [25].For the experimental data considered in this work, the three-chain λ-based model provides β = − 0.55 (data of Kawabata et al. [26]), β = −0.98 (data of Dechnarong et al. [11]), and β = − 4.3 (data of Bechir et al. [4].).The first two β-values are consistent with theoretical predictions.The larger exponent value (in magnitude) in the third example might be attributed to strain amplification required to capture the composite effect.

Conclusion
In this work, we have systematically investigated the ability of two classes of tube formulations to capture the response of rubber-like incompressible materials under different modes of deformation.The two classes were defined based on the assumed relationship between the (microscopic) tube-area contraction and the deformation of the macroscopic continuum.In λ-based approaches, the tube contraction is expressed in terms of the line-stretch of the continuum, while in ν-based approaches, it is expressed in terms of the area-stretch.To describe the free energy contribution associated with the tube topological constraints, a power-law expression was adopted, which encompasses a broad range of tube formulations proposed in the literature.In addition, three popular orientation averaging schemes were considered, namely the full-network, three-chain, and eight-chain averaging methods.
Comparisons between model predictions and experimental data suggest that ν-based formulations are in general superior to λ-based formulation, in the sense that (1) they can fit experimental data for any orientation-averaging scheme, and (2) associated best-fit parameters are similar, and thus potentially more amenable to a physical interpretation.In contrast, we found that λ-based formulations can only fit experimental data when used together with the three-chain averaging scheme.We also found that the fitting performance of the various models was directly correlated to their sensitivity to the second invariant I 2 .This confirms previous findings that I 2 dependence is necessary for good fitting performance against experimental data for different modes of deformation.Our study further identifies the combinations of tube kinematic assumption and orientation averaging scheme that lead to I 2 dependence, which had not been investigated before.
In this work, we focused on the mathematical structure of micromechanics-based hyperelastic models in relation to fitting performance and I 2 sensitivity.To this end, the power law description of the tube free energy was effectively treated as a phenomenological model, and we did not attempt to link the parameters k and β to specific molecular processes, nor did we attempt to propose a new tube formulation based on molecular considerations.The present study could however be helpful in guiding the development of improved tube models which combine a good fitting ability with physically based parameters.We also limited our analysis to the classical freely jointed chain model under the assumption of affine chain stretch to describe the elasticity of the phantom network.Other models could be considered, in particular models that relax the assumption of the affinity of the chain stretch.Such considerations are left for future work.

Appendix
Series expansions of the full-network tube free energy λ-Based full-network model Express the tube free energy ( )  c as a function of the squared line-stretch: w ( ) c (Λ) ∶= ( ) c ( ) , with Λ ∶= 2 .The Taylor series expansion of the tube energy about ⟨Λ⟩ writes, to the second order: Note the following results for the average line-stretch and its variance (see Appendix A.2 for derivation Taking the average of Eq. ( 70) over all chain orientations and using the above results then give: Note that the first-order derivative of the free energy plays no role in Eq. (72) since ⟨Λ − ⟨Λ⟩⟩ = 0.

ν-Based full-network model
The series expansion is obtained following similar steps as for the -based model.Define the squared area-stretch as Γ ∶= 2 , and express the tube free energy as w ( ) c (Γ) ∶= ( ) c ( ) .The Taylor series expansion of the tube energy about ⟨Γ⟩ writes, to the second order: The average and variance of the squared area-stretch are given by (see Appendix 7.3 for details): Taking the average of Eq. ( 73) over all chain orientations and using the above results give: Three-chain models Similar series expansion of the average tube energy can be obtained when three-chain averaging is used, instead of full-network averaging.First introduce the following notation: ⟨⋅⟩ 3C = 1 3 ∑ 3 i=1 (⋅) i .We find, for the three-chain average and variance of the squared line-stretch: and, for the three-chain average and variance of the squared area-stretch: See Appendices A.2 and A.3 for details.Series expansions of the average energies can then be expressed as: and: Numerical verification Series expansions (72) and (75) of the full-network λ-based and ν-based tube energies are compared to numerical estimates of the expressions that the series seek to approximate in Fig. 10a and b.The numerical results are exact for the three-chain case and only deviate from the exact solution in the full-network case due to the error in discretisation over the unit sphere.In the figures, the macroscopic energy Ψ = n⟨w c ⟩ has been normalised by (kn) uniaxial tension (UT), pure shear (PS), and tension (EBT).The series expansions are at small stretches, but significantly deviate from the reference numerical solution at moderate to large stretches (not shown).We note that the series expansion of the ν-based model is more accurate than that of the λ-based model.The switch in concavity of the free energy between λ-based and ν-based model is consistent with the sign of the tube modulus (32) when the sign of the exponent β changes.The series expansions of the macroscopic free energy obtained using three-chain averaging, Eqs. ( 78) and (79), are compared to reference numerical results in Fig. 10c and d for the same values of the tube exponent β.Compared to the full-network case, the series expansions are less accurate than their full-network counterparts in this range of deformation.For the considered values of β, the three-chain models are identical in their λ-based and ν-based forms, as confirmed by the numerical results.The series expansions are however not identical, and appear slightly more accurate when the expansion is taken about the squared area-stretch than about the squared line-stretch.

Average and variance of the squared line-stretch in terms of the invariants
We use the notation Λ ∶=

2
, where is the line-stretch of the continuum given by Eq. ( 6).For a chain aligned in direction n in the reference configuration, we have: = ⋅ ⋅ , where C = F T F is the right Cauchy-Green deformation tensor.The average squared line-stretch is given by: where: The variance the squared area-stretch is given by: If three-chain averaging is used, direct calculations give: where 1 = 2 3 , 2 = 1 3 , and 3 = 1 2 .The variance is given by:

Fig. 1 2 UT − − 1 UT 2 EBTFig. 2
Fig. 1 Comparison between model predictions and experimental data of Kawabata et al. [26] for uniaxial (UT) and equibiaxial tension (EBT).Non-affine tube models are considered, with a λ-based and b ν-based tube models.Full-network averaging is used and b for the λ-based and ν-based models, respectively.Corresponding values of the indicator function (56) are shown in Fig. 5c and d.The eight-chain λ-based model, similar to the λ-based full-network model, is unable to predict both responses simultaneously.However, the λ-based three-chain model is able to achieve good agreement with the data for both loading modes.In contrast, the ν-based models are able to predict the experimental response equally well under any averaging scheme.The ability of all of these models to capture simultaneously the uniaxial and equibiaxial responses is reflected in the indicator function.The λ-based three-chain model and all of the ν-based models are able to capture the experimental curve, whereas the λ-based full-network model shows limited I 2 dependence.The eight-chain λ-based model is independent of I 2 .
and b show the model predictions along with the experimental data.The λ-based models show different fitting performances under the different averaging schemes.The full-network model describes the uniaxial response well at small deformations, but departs from the experimental curve at large deformation and underpredicts the equibiaxial response up to moderate deformation.The eight-chain λ-based model insufficiently the different responses and underpredicts the initial slope of the stress-stretch curve.In contrast, the three-chain λ-based model achieves near perfect agreement.The ν-based model gives good agreement with the experimental data under all of the averaging schemes considered.

Fig. 5
Fig.5 Comparison between model predictions and experimental data of Kawabata et al.[26] for uniaxial (UT) and equibiaxial tension (EBT), using a λ-based tube models and b ν-based tube models in combination with the fullnetwork (FN), three-chain (3C), and eight-chain (8C) averaging schemes.I 2 indicator function (56) for c λ-based and d ν-based tube models

Figure 6c and d
Figure 6c and d show the indicator function (56) corresponding to the results in Fig. 6a and b.The full-network λ-based model shows I 2 dependence only at moderate stretches, whereas the three-chain λ-based model shows I 2 dependence consistent with experimental data.The eight-chain λ-based model is independent of I 2 .The indicator function calculated for the ν-based models is very close to the experimental one for all averaging schemes.The degree of I 2 dependence exhibited by all of these schemes correlates with their ability to distinguish between uniaxial and equibiaxial tension.

Fig. 6
Fig.6 Comparison between model predictions and experimental data of Dechnarong et al.[11] for uniaxial (UT) and equibiaxial tension (EBT), using a λ-based tube models and b ν-based tube models in combination with the fullnetwork (FN), three-chain (3C), and eight-chain (8C) averaging schemes.I 2 indicator function (56) for c λ-based and d ν-based tube models

Fig. 8 a
Fig. 8 a, b Comparison of Mooney stress-I 1 responses obtained using various orientation averaging schemes to the experimental data of Bechir et al. [4].c, d Corresponding I 2 indicator functions, Eq. (58)

Fig. 10
Fig. 10 Comparison between the series expansion approximations of the tube contribution to the macroscopic free energy and the reference solution obtained by numerical integration.a λ-Based full-network model, Eq. (72); b ν-based full-network model, Eq. (75); c λ-based threechain model, Eq. (78); and d ν-based three-chain model, Eq. (79).The tube exponent is set as β = − 1 in the λ-based models, and as β = 1 in the ν-based models.Three loading conditions are considered: uniaxial tension (UT), pure shear (PS), and equibiaxial tension (EBT)

Table 1
Overview of λ-based tube models of the literature

Table 2
[26]-fit parameters for models fitted simultaneously to the uniaxial and equibiaxial tensile data of Kawabata et al.[26]under different averaging schemes

Table 3
[11]-fit parameters for models fitted simultaneously to the uniaxial and equibiaxial tensile data of Dechnarong et al.[11]under different averaging schemes

Table 4
[4]t-fit parameters for models fitted simultaneously to the uniaxial and equibiaxial tensile data of Bechir et al.[4]under different averaging schemes