Continuous softening up to the onset of failure: A hyperelastic modelling approach with intrinsic softening for isotropic incompressible soft solids

This paper seeks to put forward an alternative notion for modelling the continuous softening in the finite deformation of incompressible isotropic rubber-like materials, up to the onset of failure. Based on the proposition that ‘there is no fundamental reason as to why the basic constitutive parameters of the (hyper)elastic state of a material cannot determine the natural evolution of softening and failure in that material’, the concept of hyperelasticity with intrinsic softening is posited. In contrast to the currently available theories of continuum damage mechanics and energy limiters, which superpose external parameters onto the basic hyperelastic model, the alternative approach presented here postulates that if a hyperelastic strain energy function accommodates a comprehensive set of constitutive parameters and has an appropriate functional form, it will intrinsically capture the observed continuous softening in the stress–deformation curves of soft solids, up to the onset of failure. Examples of this notion are presented using a specific hyperelastic model previously proposed by the author, applied here to extant datasets of a wide variety of isotropic incompressible soft solids, ranging from soft tissues to protein gels and 3D printed biomaterials, that include the softening behaviour. The favourable modelling results obtained via this approach are demonstrated. It is conferred that these results portend a versatile modelling tool for application to the whole-range deformation of soft solids, up to the onset of failure, and make the case for devising a unified theory for modelling both continuous and discontinuous softening in the finite deformation of rubber-like materials.


Introduction
The softening phenomenon in the finite deformation of rubberlike materials may be classified as continuous and/or discontinuous.Discontinuous softening is referred to the softening observed in the stress-deformation curves of the unloading path, compared with the primary loading path, governed by the maximum deformation experienced by the specimen during the preceding loading phase (see Fig. 1a).This behaviour is akin to the well-known Mullins effect [1,2], and a widely used theoretical framework for modelling this behaviour is the seminal pseudo-elasticity approach [3].Continuous softening, on the other hand, reflects itself progressively in the primary loading stressdeformation curve with increase in deformation, eventually leading to the structural failure of the sample (e.g., Fig. 1b).A versatile modelling approach to capturing this behaviour is perhaps that of Volokh [4,5], using energy limiters in conjunction with the traditional hyperelastic stored energy functions.Another widely used framework in this context is the continuum damage mechanics framework [6], which is also capable of capturing the discontinuous softening effects [7,8].
The foregoing two modelling approaches to continuous softening, however, suffer from unfavourable ramifications in terms of the underlying rationale and modelling applications.Take the continuum damage mechanics approach first.As astutely observed in [9], due to the presence of the so called reduction factor (1 − ) in this approach and that, unlike the pseudo-elasticity framework, this factor is not inactive over any range of deformation, it is not possible to independently obtain the values of the constitutive parameters of the basic (hyper)elastic material by fitting a continuum damage mechanics model to the experimental data.That is, the obtained (hyper)elastic parameters will too be influenced by the presence of the damage factor.Now let us consider the energy limiter approach.Since the softening and failure behaviours in this approach are essentially controlled by only the energy limiter parameters (e.g.,  and ), on using a basic hyperelastic strain energy function  with fixed values of model parameters it is possible to predict and obtain entirely different softening and failure behaviours (i.e., only by varying the values of  and ).See [5] for modelling examples using Yeoh's cubic [10] and Ogden [11] models.In other words, the use of this approach implicitly considers https://doi.org/that the softening/failure behaviour is by and large independent of the constitutive (elastic) properties of the material, and thus the same material (with the same constitutive elastic parameter values) can exhibit different softening and failure behaviours.
In view of the aforementioned disadvantages to the currently existing approaches for modelling continuous softening, I wish to put forward here an alternative proposition which is free from superposing damage parameters or energy limiters onto the basic hyperelastic properties, based on the following question: is there a fundamental reason as to why the basic constitutive parameters of a hyperelastic material cannot determine its softening/failure behaviour?Put differently, why should not the set of the basic constitutive parameters of the (hyper)elastic state of a material determine the natural evolution of softening and failure in that material?I therefore stipulate that if the basic hyperelastic strain energy function  incorporates the appropriate constitutive parameters, i.e., is comprehensive enough and has a rich-enough functional form, it will intrinsically accommodate the softening behaviour of that material.Note that since we now wish to model the softening behaviour as a natural evolution within the (hyper)elastic response, the application of this modelling approach is restricted up to the onset of failure.To cast a precise definition, the failure point is where the maximum stress is reached, and d d = 0 (or, equivalently,   = 0 when dealing with discrete experimental datapoints), where  denotes the stress and  is the stretch.See Fig. 1b for a graphical representation.Note that the behaviour after the failure point cannot be captured by this approach, since the failure process introduces instabilities and involves irreversibilities/dissipative effects that are not included in the framework of hyperelasticity.
To back up this stipulation, here I propose the use of a specific strain energy function  previously presented by the author [12], and demonstrate its application to the softening behaviour, up to the onset of failure, for various isotropic incompressible soft solids ranging from soft tissues to protein gels and 3D-printed biomaterials.Note that the softening behaviour of interest here is distinct from the stress/strain softening observed in natural behaviour of, for example, filled rubbers (e.g., [13]) or liquid crystal elastomers (e.g., [14]).The softening behaviour of interest here is that which leads to (the onset of) failure.The model preliminaries will be presented in Section 2, and the modelling results in application to the considered extant experimental datasets are demonstrated in Section 3. Concluding remarks are conferred in Section 4, discussing the implication of this modelling approach in devising a unified theory which captures both the continuous and discontinuous softening, including the Mullins effect.

Model preliminaries
The strain energy function  considered here for the purpose of showcasing our hyperelasticity with intrinsic softening approach is of the form [12]: , ( where , ,  and  are model parameters with: and parameter  is related to the infinitesimal shear modulus  0 via: Note that for the ln (•) function to be well-defined in a general deformation we must have: subject to the incompressibility constraint: A detailed theoretical background to devising this model has been presented in [12], and its successful application to multiaxial deformation datasets of a variety of isotropic incompressible soft solids has been demonstrated in [12] and [15].For the interested reader, we briefly note that this model is based on a generalisation of the class of limiting chain extensibility models whose response functions were first formulated in [16] based on a [1/1] Padé approximant, and as shown in [12] is the parent to many of the existing landmark models in the literature including the neo-Hookean [17], Mooney-Rivlin [18], Gent [19] and Ogden [11] models (in either the one-term form in Eq. ( 1) or the multi-term form as presented in [12]).In the special case where  = 3, this model recovers that of Anssari-Benam et al. [20] which was devised for application to the deformation of the brain tissue.Other specialisations with  = 2 and { = 3 ,  = 2} recover earlier models proposed by Anssari-Benam and Horgan [21] and Anssari-Benam and Bucchi [22,23], respectively.Due to these features, this model provides a good representative example of a comprehensive model for consideration here, from the range of existing hyperelastic models in literature (see, e.g., reviews [24] and [25]).
From the model in Eq. (1), specific (Cauchy)stress-stretch relationships for each deformation mode may be derived.For the purpose of this study, since most of the datasets characterising continuous softening and failure are presented using uniaxial and (equi-)biaxial deformations, it suffices to note the following two relationships: ) , (5) where  denotes the Cauchy stress and  is the principal stretch.We further note that solving d d = 0 establishes the failure stretch   ; the stretch level at which the failure occurs and  =   .Using the relationships in Eq. ( 5) we find: for the uniaxial deformation, and: for the equi-biaxial deformation.While Eqs. ( 6) and ( 7) are not amenable to finding a closed-form analytical solution, it can be verified numerically that for various combinations of , ,  and , these equations have a positive real root.Numerical examples will be presented in Section 3, demonstrating a close agreement with the experimental data.
Remark.It is instructive to note that d d = 0 does not result in a positive real root in many of the existing models in the literature a priori; rendering those models incapable of predicting the softening behaviour which is of interest here.For example, on using the neo-Hookean model   = 1 For the purpose of modelling the experimental data in the next section, the relationships in Eq. ( 5) are fitted to the relevant dataset(s), where the best fit is obtained by minimising the residual sum of squares (RSS) function defined as: RSS = ∑  (   −   ) 2  , with  being the number of data points.The minimisation is performed via an inhouse developed code in Matlab ® , using the genetic algorithm (GA) function.The coefficient of determination R 2 values are reported as a measure of the goodness of the obtained fits.

Modelling softening up to the onset of failure
Since most studies in the literature utilise uniaxial tensile tests to characterise continuous softening/damage in isotropic soft solids, multiaxial datasets purporting this phenomenon appear scarce, at least to the author's knowledge.However, theoretical predictions of multiaxial deformations using customised models developed for particular applications may be found in the literature.One such study is due to Guo and Zaïri [26], where they present predictions of the softening and failure in silicone specimens using a micro-structurally derived model for uniaxial and equi-biaxial deformations.We start our modelling campaign by presenting their predicted data and the fitting results using the model in Eq. ( 1).The plots in Fig. 2 show the correlation, and Table 1 summarises the model parameter values.The values of   , calculated from Eqs. ( 6) and ( 7) using the obtained model parameters versus those established from the data, have also been presented therein.
Next we consider two datasets pertaining to the uniaxial deformation (up to the onset of failure) of two biomaterials, namely whey protein gels due to Liu and Böl [27] and 3D-printed zigzag composite

Table 2
Model parameter and   values for the softening behaviour of a typical whey protein gel [27] and a 3D-printed composite suture [28] specimens.

Table 3
Model parameter and   values for the softening behaviour of the various soft tissue specimens in Fig. 4. sutures due to Liu and Li [28].These two datasets exhibit a very distinct softening trend, as shown in the plots of Fig. 3.The proposed model captures both datasets favourably.The ensuing model parameters are presented in Table 2.The experimental values of   for both datasets versus those calculated using Eq. ( 6) with the obtained model parameters have also been reported.
Softening is an important biomechanical feature of biological tissues, as it may be indicative of damage and failure, and thus is critical to characterise and/or predict.It is therefore in order here to consider some examples of softening in soft tissues.Accordingly, in the following we present the modelling results of the softening behaviour (up to the onset of failure) of a range of different tissues, including rectus sheath due to Martins et al. [29], liver parenchyma due to Untaroiu and Lu [30], abdominal aortic aneurysm (AAA) due to Volokh [31], and corneal stroma due to Xiao et al. [32], all under uniaxial loading.The plots in Fig. 4 present the modelling results, and Table 3 contains the model parameter values for each fit along with the corresponding values of   calculated using Eq. ( 6) and those of the experimental data.

Some remarks for discussion and conclusion
My intention in this concise manuscript was to put forward the question as to: 'why should we not expect that the basic constitutive parameters of the hyperelastic state of a material could determine the natural evolution of softening and failure in that material'?On this basis, a proposition was posited that a comprehensive hyperelastic model, comprehensive in that it accommodates an appropriate set of constitutive parameters and has a suitable functional form, may be capable of capturing the whole-range deformation of the subject soft solid specimen up to the onset of failure, including its softening behaviour.To this end, a case of hyperelasticity with intrinsic softening modelling was provided using a previously proposed strain energy function  by the author [12].Examples of modelling the continuous softening behaviour for a variety of rubber-like materials were considered and presented, reaffirming the proposed proposition.It is hoped that these results furnish a new venue for constitutive modelling of softening and failure in finite deformation of soft solids.However, there are some points of discussion to confer that may help with further maturation and future development of this concept.1) and those of the model in [26] for: (a) uniaxial; and (b) equi-biaxial deformations.The relationships in Eq. ( 5) were simultaneously fitted to the datasets.Note that the scales of the axes are different for a better clarity of presentation.Fig. 3. Modelling results for the softening behaviour (up to the onset of failure) of: (a) a typical whey protein gel [27]; and (b) a 3D-printed composite suture [28] specimen using the model in Eq. ( 1).The relationship in Eq. ( 5) 1 was fitted to each dataset.Note the different scale of the axes for a better clarity of presentation.
First is perhaps the implication of this modelling approach.As was briefly reviewed in Section 1, the continuum damage mechanics approach currently offers a modelling tool to capture both the continuous and discontinuous softening.However, it suffers from setbacks: in addition to not being able to calibrate the constitutive parameters of the (hyper)elastic state of the subject material in this approach (due to the presence of the reduction factor; see, e.g., [9]), it is well established that the continuum damage mechanics approach encounters difficulties capturing the strain hardening behaviour and also the Mullins effect when the level of softening is high (see, e.g., [33,34]).These drawbacks are not present within the pseudo-elasticity theory; however, that theory does not capture the continuous softening.With the results put forward in this paper, it may now be possible to enhance the framework of pseudo-elasticity and devise a unified theory for capturing both the continuous and discontinuous softening, free from the aforementioned drawbacks.The pseudo-elasticity framework has been extended to account for the permanent set [35], anisotropy [36,37] and recently the induced anisotropy [38] and rate-dependent behaviour [39] too.By incorporating the proposed model here into these extensions, it is plausible to achieve a single theoretical approach that captures features such as the Mullins effect, continuous softening, the permanent set, anisotropy and rate-effects, all in one platform.
Second, an important limitation must be noted.The modelling approach presented here is strictly phenomenological and not predictive.Calibrating the model and obtaining values for the constitutive parameters such that the model produces the softening behaviour requires fitting the model to datasets which already include the softening datapoints in them.If the model is fitted with the same dataset but without the softening part, it will produce a different fit with different parameter values which are not guaranteed to result in the observed softening -see Appendix B of [40] for illustrative examples on this issue using other models.Therefore, it is important that in the fitting process, as was done in the current study, the softening data is considered.For predictive applications, other approaches to modelling the softening effects may be more suitable, including microstructuralbased models (see; e.g., [41,42]).It should also be noted that the energy limiters approach has proved a versatile tool for modelling the failure localisation and crack propagation in elastomers; see, e.g., the works of Volokh and co-workers [43][44][45][46], while the generalisation of the proposed model here for such applications remains untested.In addition, for usages such as modelling the thermoelastic effects on failure of rubber-like materials, where the assumption of separating the material parameters linked with stiffness (e.g., the hyperelastic constitutive parameters) from those connected with strength (e.g., the energy limiter parameters  and ) is justified; for example see [47], the energy limiter approach may be considered an appropriate concept.
A separate approach to softening, based on the spectral formulation for isotropic materials as in the works of Shariff [48][49][50], also exists in the literature, which to the author's knowledge, has been specialised for Mullins-type effects (i.e., discontinuous softening).The hyperelasticity with intrinsic softening model presented here is distinct from that approach, as it captures the continuous softening.I also note a disparate softening effect, namely the ''localised strain softening'', defined, for example, in the works of Haughton and Merodio [51,52] Fig. 4. Modelling results for the softening behaviour (up to the onset of failure) of a verity of soft tissue specimens from: (a) rectus sheath [29]; (b) liver parenchyma [30]; (c) abdominal aortic aneurysm (AAA) [31]; and (d) corneal stroma [32], using the model in Eq. ( 1).The relationship in Eq. ( 5) 1 was fitted to each dataset.Note the different scale of the axes for a better clarity of presentation.
as showing ''a limited amount of strain stiffening followed by strain softening and then by further stiffening''.The focus of their work has been on pathological soft tissues, where the softening behaviour and the resulting instabilities (bifurcations) that ensue were analysed and modelled.The scope of the current work centres on modelling the continuous softening phenomenon; although it would be of interest to investigate the predictions of the proposed model here for such instabilities and softening effects, in order to establish a model for capturing the holistic softening behaviours of various kinds in soft solids.
Lastly, given the promising results here regarding the application of the model to capturing the softening behaviour of soft tissues in the isotropic context, it may be merited to extend the model for the anisotropic case too.By nature, the inclusion of anisotropy results in increased number of model parameters, which combined with the additional damage parameters render a high number of total model parameters for a meaningful calibration, implementation and interpretation.An approach such as the one presented in this study would not need additional damage variables for capturing the continuous softening effect, and as a result would help reducing the total number of parameters for optimisation, calibration and modelling.

Declaration of competing interest
The author declares no conflict of interest.

Fig. 1 .
Fig. 1.Illustrative representation of: (a) discontinuous softening, akin to the Mullins effect, reflected in the unloading paths; and (b) continuous softening which occurs progressively up to the failure point, becoming more visible closer to the failure.

2 ( 2 ,
which is clearly not valid.As another example, for the Ogden model   = ∑ possesses no positive real roots subject to the empirical condition     > 0 etc.

Fig. 2 .
Fig. 2. Comparison between the predictions of the proposed model in Eq. (1) and those of the model in[26] for: (a) uniaxial; and (b) equi-biaxial deformations.The relationships in Eq. (5) were simultaneously fitted to the datasets.Note that the scales of the axes are different for a better clarity of presentation.