The researchers interested in nonlinear elasticity and in particular in its role in describing phase transition phenomena are no doubt greatly indebted to J.L. Ericksen for his 1975 work Equilibrium of Bars [11], where in a simple and straightforward manner he extended the predictive power of this theory in considering the delicate problem of material instabilities and phase transition phenomena. After many efforts in determining so called constitutive inequalities for constitutive equations, not given by thermodynamics [21], he proposed a new approach in determining the correct form of constitutive assumptions. In his words [12] “my general experience tells me that it is neither wise nor fruitful to impose any restrictions upon the constitutive equations, in general... I merely propose adoption of the mathematician’s criterion: the weaker the hypothesis, the better the theorem”. The main assumption of [11] was the possibility of instability material domains, where the energy density was not rank-one convex, or in other words the possibility of multiwells energy densities. In the conclusion of [11] we read: “Linear elasticity theory is, by its nature, ill-designed to treat such problems, but we should not saddle nonlinear elasticity theory with deficiencies not inherent in it”.

Ericksen’s work has put the basis to enable us to review and understand non-linear elasticity in particular in the direction of prediction of the possible coexistence of equilibrium phases in material science. His seminal paper originated a stream of fundamental results that for more than twenty years reshaped the potential of a mathematical theory fundamental to all continuum mechanics, about many other J. Ball, R.D. James, M. Gurtin, J. Knowles, R. Abeyaratne, E Sternberg results in [1, 2, 4, 10, 20, 23, 24]. Very fruitful has then been the extension to the discrete versions of these models as proposed by I. Muller and P. Villaggio [25] and fully analyzed and adopted in the analysis of shape memory materials by G. Puglisi and L. Truskinovsky [29, 30]. Possible extension to the analysis of non local interactions, mimicking three dimensional effects have been considered in [26]. Recent extension to lower scales thermomechanical problems exhibited the efficacy of this type of approaches also for molecular scale material transition effects [13] and related multiscale modelling [31].

Roger Fosdick was not only one of the pioneers of this important extensions of non-linear elasticity to problems with non (rank-one) convex energy densities, as evidenced by his work with Dunn in 1980 [10], but he established fundamental results in the thermodynamics of continua and nonlinear elasticity theory representing a reference for the entire scientific community. First, the interest in deformation classes that, especially in the case of incompressible materials, reduces the general balance equations through the semi-inverse method to a less complex system of equations (e.g., from a system of PDEs to a system of ODEs) allowing for analytic solutions and deep results also in the complex field of three-dimensional non linear elasticity. Such analytical results allowed for a deeper physical interpretation with both theoretical interest and important applications for real problems. This series of pioneering works, beginning more or less in 1983 and continuing without interruption until 1996 [14–18], has been a landmark for an entire generation interested in the extension of classical methods of non-linear elasticity and thermodynamics.

This work is inspired to the recalled approaches in solving boundary values problems in continuum mechanics describing solid-solid phase transition phenomena. Specifically, we explore the possibility of describing damage phenomena in elastomeric materials modeled as a material hard→soft phase transition. In particular we consider elastomeric materials. These are amorphous materials and therefore their mesoscopic description is highly phenomenological typically based on the kinetic theory of rubber [34]. Despite this, it is widely recognized that approaches based on micro-mechanical interpretations allow for the deduction of very effective material modelling to describe both elastic and inelastic phenomena in rubber-like materials. Amorphous materials can be modelled as a mixture of different materials: in some sense a family of elastic materials [5, 6]. Their damage and fracture behavior can be considered as a homogenized effect of complex phenomena of links breaking and recrosslinking, corresponding to a distribution of activation and fracture thresholds of the different constituent materials at the microstructure scale. At low strains, the material response is mainly regulated by the molecular network properties such as intermolecular or chain-filler interactions. This first stage is characterized by a continuous phenomenon of link scission leading to a softening effect. On the contrary, at larger strain the behavior is mainly independent from the network properties such as density and type of strands, whereas it is regulated by the elastomeric matrix with a hardening effect induced by chains re-orientation.

A number of microstructure inspired or multiscale models have been proposed to describe the recalled phenomenology. We recall the papers by Rajagopal, Wineman, and others based on a continuous transition among different network natural configurations controlled by a scalar function of the strain invariants [22, 33, 35]. Here we consider a different approach based on the assumption of two distinct material phases, undamaged and damaged, characterized by a different energy density. In particular we focus on the simple case of isotropic materials and damage, and consider elastic energy of Neo-Hookean and Gent type. As well known in this case the energy depends on the only first invariant \(I_{1}\) of the left Cauchy Green tensor that gives a measure of the averaged (with respect to orientation) of molecule stretch at the network scale under the classical affinity hypothesis [34]. We then assume the existence of a damage threshold \(I_{1}=I_{cr}\) leading to a decreasing of the shear modulus \(\mu \). Based on the physical interpretation of \(\mu \) as proportional to the number of molecules per unit volume, this assumption corresponds to a transition to a state with a fraction of the molecules in the broken state. We then evaluate the corresponding dissipation. Based then on a Griffith type variational approach, we search for the minima of the total (elastic and damage) energy. As in the problems considered by Fosdick in [14–18], the resulting functional to be minimized is non-convex and for this reason it results in a coexistence of two phases: damaged and virgin material.

While the choice of two states, Neo-Hookean and Gent materials allow us for fully analytical description of the results, with a transparent physical interpretation, several extension of the models can be considered. In particular, we argue that by extending this approach to more phases can be fully descriptive of many damage effects and resulting localization phenomena [6].

1 Constitutive Assumptions

We consider here isotropic damageable incompressible materials. In particular, we analyze the simple case of a two-state (virgin and damaged) material. In the first state the system has a classical Neo-Hookean behavior with an energy density

$$ \bar{\Phi}(\boldsymbol{F})= \Phi _{u}(I_{1})=\frac{\mu}{2}(I_{1}-3) $$

whereas in the second state we consider again a Neo-Hookean behavior

$$ \bar{\Phi}(\boldsymbol{F})= \Phi _{d}(I_{1})= \frac{\alpha \, \mu}{2}(I_{1}-3)+ \frac{\beta \, \mu}{2}, $$

where \(\alpha <1\) measures the damage in the material (damaged vs undamaged modulus) and \(\beta \) measures the dissipation associated with the damage effect. Here \(\boldsymbol{F}\) is the deformation gradient, \(\boldsymbol{B}=\boldsymbol{F}^{T}\boldsymbol{F}\) is the left Cauchy Green tensor and \(I_{1}=\boldsymbol{F}\cdot \boldsymbol{F}\) is its first invariant. Here and in the following we use the index \(u\) and \(d\) to denote the undamaged and damaged material phases, respectively.

While anisotropic effects can be very important in the description of damage and associated residual strains in soft materials [19, 31], here to focus on the main physical effect of describing damage as a material phase transition, we neglect anisotropic effects. As a result, following [9] we assume that damage is regulated by the only first invariant \(I_{1}\). In the case of macromolecular materials it is interesting to observe that \(I_{1}\) measures the average elongation along the different direction in each material point, thus if we assume, according with the classical affinity hypothesis [34], that the molecular deformations coincide with the macromolecular ones and that the damage depends on the macromolecular deformation, then we may assume that the damage is regulated by \(I_{1}\) itself. In particular, under our simplifying two-state hypothesis we assume that there exists a damage threshold \(I^{cr}_{1}\) such that the system is undamaged for \(I_{1}\in (0,I_{1}^{cr})\) and damaged for \(I_{1}\geq I_{1}^{cr}\). Thus the energy function is assigned by (see Fig. 1)

$$ \Phi (I_{1})=\left \{ \textstyle\begin{array}{l@{\quad}l@{\quad}l@{\quad}l} \displaystyle \frac{\mu}{2}(I_{1}-3) & \mbox{ if } & I_{1}^{max}< I_{1}^{cr}, &\mbox{undamaged state}, \vspace{0.2 cm} \\ \displaystyle \frac{\alpha \mu}{2}(I_{1}-3) + \frac{\beta \mu}{2} & \mbox{ if } & I_{1}^{max}\geq I_{1}^{cr}, &\mbox{damaged state}, \end{array}\displaystyle \right . $$
(1.1)

where

$$ \beta =(1-\alpha ) (I_{1}^{cr}-3) $$
(1.2)

measures the dissipation function and \(I_{1}^{max}\) represents the maximum attained value of \(I_{1}\) in the loading history. Correspondingly the Cauchy stress is

$$ \boldsymbol{T}(I_{1})=\left \{ \textstyle\begin{array}{l@{\quad}l@{\quad}l} \displaystyle \frac{\mu}{2}\boldsymbol{B}-p \boldsymbol{I}, & \mbox{ if } & I_{1}^{max}< I_{1}^{cr}, \vspace{0.2 cm} \\ \displaystyle \frac{\alpha \mu}{2}\boldsymbol{B}-p \boldsymbol{I} ,& \mbox{ if } & I_{1}^{max}\geq I_{1}^{cr}. \end{array}\displaystyle \right . $$
(1.3)
Fig. 1
figure 1

Constitutive assumptions. a) Damageable two-phase Neo-Hookean material (1.1): \(\alpha =0.5\), \(I_{1}^{cr}=4\). b) Damageable two-phase Gent material (1.4): \(J_{m}=4\); \(\alpha =2\); \(I_{1}^{cr}=6\)

Full size image

Some comments are in order. By recalling that according with the statistical mechanics interpretation of the shear modulus we have \(\mu =N k_{B} T\), where \(N\) is the number of chains per unit volume, \(k_{B}\) is the Boltzmann constant, and \(T\) is the temperature. Moreover it is easy to interpret \(I_{1}\) as an averaged measure of the elongation in all directions. Thus the underlying physical idea of the model that due to damage a fraction \((1-\alpha )\) chains are broken in the damaged state. Thus under a simple isotropic assumption we simply describe damage as a variation of the parameter \(\mu \). Observe also that based on (1.2) we can give an energetic interpretation to the damage criterion. Indeed if the critical damage threshold \(I_{1}^{cr}\) is assigned, the ‘dissipated’ energy associated to damage is \(\frac{\mu \beta }{2}= \frac{\mu }{2} (1-\alpha ) (I_{1}^{cr}-3)\), i.e. the difference of the elastic energy in the undamaged and damaged state at the critical damage threshold.

It is important to observe that more general assumptions with more phases \(w_{i}(I_{1})=\frac{\alpha _{i} \mu}{2}(I_{1}-3) + \frac{\beta _{i} \mu}{2}\) can be considered, by assuming different limit conditions \(I_{1}^{max}=I_{i}^{cr}\), \(i=1,\ldots,n\). In this perspective (see [34]), we may relate the variation of the moduli to the variation of the number \(N\) of chains per unit volume, i.e. \(\alpha _{i} \sim \Delta N_{i}\), where \(\Delta N_{i}\) is the number of chains broken at the \(i\)-th step and similarly we may assume that the damage energy is proportional to \(\Delta N_{i}\). As a result we can assume a simple relation \(\beta _{i}= k \alpha _{i}\) with \(k\) a single constitutive parameter for the progressive damage. Possible extension to the case of a continuous distribution of damageable materials can be also considered [5, 8].

This simple assumption let us attain a full analytical treatment of the problem with a clear physical interpretation of the results, on the other hand it does not take in consideration important effects such as entropic hardening observed in rubberlike materials [28] and density and volume changes associated with damage induced cavitation effects [3].

Thus, as a second constitutive assumption, in order to take care of the hardening effect, we consider a damageable Gent energy

$$ \Phi (I_{1})=\left \{ \textstyle\begin{array}{l@{\quad}l@{\quad}l} \displaystyle - \frac{\mu}{2}J_{m} \ln \left (1-\frac{I_{1}-3}{J_{m}} \right ), & \mbox{ if } & I_{1}^{max}< I_{1}^{cr}, \vspace{0.2 cm} \\ \displaystyle - \frac{\alpha \mu}{2}J_{m} \ln \left (1- \frac{I_{1}-3}{\alpha J_{m}}\right )+ \frac{\beta \mu}{2}, & \mbox{ if } & I_{1}^{max}\geq I_{1}^{cr}, \end{array}\displaystyle \right . $$
(1.4)

where \(J_{m}\) is a dimensionless parameter defined as \(J_{m} = I_{m} -3\) and \(I_{m}\) indicating the (damage dependent) limiting value for \(I_{1}\). This corresponds to a stress strain relation

$$ \boldsymbol{T}(I_{1})=\left \{ \textstyle\begin{array}{l@{\quad}l@{\quad}l} \displaystyle \frac{\mu J_{m}}{J_{m}-(I_{1}-3)}\boldsymbol{B}-p \boldsymbol{I}, & \mbox{ if } & I_{1}^{max}< I_{1}^{cr}, \vspace{0.2 cm} \\ \displaystyle \frac{\alpha \mu J_{m}}{\alpha J_{m}-(I_{1}-3)} \boldsymbol{B}-p \boldsymbol{I} ,& \mbox{ if } & I_{1}^{max}\geq I_{1}^{cr}. \end{array}\displaystyle \right . $$
(1.5)

In this case we assume \(\alpha >1\) that corresponds to both a decreasing of the tangent modulus in (1.5) and an increase of \(J_{m}\) to \(\alpha J_{m}\) due to damage. This can be again interpreted based on previous microstructure description. As the stretch increases the shortest molecules break. Thus, since \(J_{m}-3\) represents an averaged value of the stretch molecular threshold due to limit chain extensibility [27], after damage we can expect an increase of this value. The explicit dependence of \(\alpha \) and \(\beta \) can be deduced on the basis of microscopic considerations. For example, as described in details in [9] in the case of biopolymer materials or thermoplastic polyurethans [32], damage can be addressed to the progressive unfolding of hard domains into soft unfolded domains. As a result the number \(n\) of monomers (or Kuhn segments) increases. This corresponds to a variation of the contour length \(l_{c}\sim n\) so that \(\alpha \sim n^{3}\) if we consider the case of isotropic damage and the dissipation is \(\beta \sim N n\). In this case then we have \(\alpha \sim \beta ^{3}\). While we may assume that \(J_{m}\) depends only on the elastic part of the deformation, more general assumptions can be considered as in [9] where natural configuration variations and growth mechanics are considered in the constitutive assumptions. Moreover also in this case we may consider progressive damage.

In the simplest framework of a two phase Gent material, here considered, the dissipation energy can be evaluated as

$$ \beta =-J_{m} \ln \left ( 1-\frac{I^{cr}_{1}-3}{J_{m}} \right )+ \alpha J_{m} \ln \left ( 1-\frac{I^{cr}_{1}-3}{\alpha J_{m}} \right ) $$

which again represents the energy difference of the two states at the damage threshold \(I_{1}=I_{1}^{cr}\).

2 Damage as a Phase Transition

We now show through different explicit examples the possibility of describing physically based damage behaviors.

2.1 Piecewise Homogeneous Deformations: Simple Shear

Consider first the simple energetic assumption (1.1) for the rectangular body in Fig. 3 under the hypothesis of assigned displacement \(\delta \) and simple shear deformations \(\boldsymbol{f}\) assigned as

$$ \left \{ \textstyle\begin{array}{l} x=X+kY, \vspace{0.2 cm} \\ y=Y, \vspace{0.2 cm} \\ z=Z, \vspace{0.2 cm} \end{array}\displaystyle \right . $$
(2.1)

where by \(\boldsymbol{X}=(X,Y,Z)\) we denote a material point in the reference configuration and by \(\boldsymbol{x}=\boldsymbol{f}(\boldsymbol{X})=(x,y,z)\) its image under the deformation. The corresponding deformation gradient is then \(\boldsymbol{F}=\boldsymbol{I}+ k \boldsymbol{e_{1}}\otimes \boldsymbol{e_{2}}\), so that the left Cauchy Green tensor is \(\boldsymbol{B}=\boldsymbol{F}\boldsymbol{F}^{T}\). Thus,

$$ \boldsymbol{F}=\left [ \textstyle\begin{array}{c@{\quad}c@{\quad}c} 1 & k & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\displaystyle \right ] \quad \boldsymbol{B}=\left [ \textstyle\begin{array}{c@{\quad}c@{\quad}c} 1+k^{2} & k & 0 \\ k & 1 & 0 \\ 0 & 0 & 1 \end{array}\displaystyle \right ] $$

and \(I_{1}=3+k^{2}\). Here and in the following we indicate with \(\boldsymbol{e}_{i}\), \(i=1,2,3\) the reference frame vectors.

2.1.1 Damageable Neo-Hookean Material

Consider first the case of damageable Neo-Hookean material \(\bar{\Phi}(\boldsymbol{F})=\Phi (I_{1})=\Phi (3+k^{2})=\phi (k)\) with

$$ \phi (k)=\left \{ \textstyle\begin{array}{l@{\quad}l@{\quad}l} \frac{\mu}{2}k^{2} & \mbox{ if } & k^{max}< k^{cr}, \vspace{0.2 cm} \\ \frac{\alpha \mu}{2}k^{2} + \frac{\beta \mu}{2} & \mbox{ if } & k^{max} \geq k^{cr} , \end{array}\displaystyle \right . $$
(2.2)

where \(I_{cr}=3+k_{cr}^{2}\).

The Cauchy stress is then given by \(\boldsymbol{T}=\mu \boldsymbol{B} - p \boldsymbol{I}\) and the Piola stress is given by \(\boldsymbol{S}=\mu \boldsymbol{F} - p \boldsymbol{F}^{-T}\) where \(p\) represents a pressure taking care of the incompressibility constraint. By imposing the boundary conditions \(T_{22}=T_{33}=0\) we obtain \(p=p_{u}= \mu \) for the undamaged phase and \(p=p_{d}=\alpha \mu \) for the damaged phase.

Thus,

$$ \boldsymbol{T}_{u}=\left [ \textstyle\begin{array}{c@{\quad}c@{\quad}c} \mu k^{2}& \mu k & 0 \\ \mu k & 0 & 0 \\ 0 & 0 & 0 \end{array}\displaystyle \right ], \quad \boldsymbol{S}_{u}=\left [ \textstyle\begin{array}{c@{\quad}c@{\quad}c} 0& \mu k & 0 \\ \mu k & 0 & 0 \\ 0 & 0 & 0 \end{array}\displaystyle \right ], $$

with a shear stress

$$ \tau _{u}(k)=\mu k. $$

And for the damaged phase,

$$ \boldsymbol{T}_{d}=\left [ \textstyle\begin{array}{c@{\quad}c@{\quad}c} \alpha \mu k^{2}&\alpha \mu k & 0 \\ \alpha \mu k & 0 & 0 \\ 0 & 0 & 0 \end{array}\displaystyle \right ], \quad \boldsymbol{S}_{d}=\left [ \textstyle\begin{array}{c@{\quad}c@{\quad}c} 0& \alpha \mu k & 0 \\ \alpha \mu k & 0 & 0 \\ 0 & 0 & 0 \end{array}\displaystyle \right ], $$

with a shear stress

$$ \tau _{d}(k)=\alpha \mu k. $$

According with previous constitutive hypothesis, under the normalizing assumption of unitary section in the \(\pi _{\boldsymbol{e}_{1},\boldsymbol{e}_{3}}\) plane, we may consider single phase solutions as in (2.1) with energy

$$ \Psi = \int _{0}^{L} \frac{\mu}{2} k^{2}dY=\frac{\mu}{2 L} \delta ^{2}, $$
(2.3)

where we imposed the kinematic constrain

$$ \int _{0}^{L} k(Y) dY=k L=\delta . $$
(2.4)

We may then, by following classical approaches in the phase transition theory [20] and our energetic assumption, describe the evolution of damage as a phase transition phenomenon. While in the tridimensional case the analysis of continuum bodies with non (quasi)-convex energy densities constitute a very complex (typically not analytically solvable) mathematical problem, since we are here in the special cases of general Neo-Hookean material and under the special assumption on piecewise homogeneous shear deformations, as in the case of simple shear, with interfaces having normal \(\boldsymbol{n}\equiv \boldsymbol{e}_{2}\), we may use the necessary conditions to obtain local minimizers as in [20].

Thus if \(\boldsymbol{F}_{1}\) and \(\boldsymbol{F}_{2}\) are the deformation gradients corresponding to shears \(k_{1}\) and \(k_{2}\) respectively, the solutions must respect

$$ \textstyle\begin{array}{l@{\quad}l} \tau (k_{2})=\tau (k_{1})=:\tau _{Max}& \mbox{Interface equilibrium,} \vspace{0.4 cm} \\ \phi (k_{2})-\phi (k_{1})=\tau _{Max}(k_{2}-k_{1}) & \mbox{Maxwell condition.} \end{array}\displaystyle . \vspace{0.2 cm} $$
(2.5)

Observe that the considered two phase solutions verify the kinematic compatibility Hadamard condition (\(\boldsymbol{n}\) being the normal to the interface)

$$ \boldsymbol{F}_{2}-\boldsymbol{F}_{1}=\boldsymbol{a}\otimes \boldsymbol{n} =(k_{2}-k_{1}) \boldsymbol{e}_{1}\otimes \boldsymbol{e}_{2} $$
(2.6)

and the equilibrium condition at the interface, namely

$$ \boldsymbol{S}(\boldsymbol{F}_{2})\boldsymbol{n}=\boldsymbol{S}( \boldsymbol{F}_{1})\boldsymbol{n}=\tau _{Max}\boldsymbol{e}_{1}. $$

Thus, in the present case, we may reduce the three-dimensional problem within a one dimensional setting. Indeed we have to compare the single phase solution with energy (2.3) with two-phase solutions. We then search for the deformation minimizers \(k=\hat{k}(Y)\) of the total potential energy

$$ \min _{k}\left [ \int _{0}^{L} w(k(Y)) \, dY - \tau \left (\int _{0}^{L} k(Y) \, dY -\Delta \right )\right ], $$
(2.7)

with \(\tau \) representing a Lagrange multiplier taking care of the kinematical constraint (2.4) physically interpreted as the reaction shear applied by the constraint on the edge \(Y=L\). We are then lead to the classical Maxwell construction represented in Fig. 2. The two phases solutions according with (2.2) and (2.5) the Maxwell stresses and shears are:

$$ \textstyle\begin{array}{lll} k_{a}&=&\displaystyle \vspace{0.2 cm} \frac{\tau _{Max}}{\mu}= \frac{\sqrt{\alpha \, \beta }}{\sqrt{1-\alpha }}; \\ k_{b}&=&\displaystyle \vspace{0.2 cm} \frac{\tau _{Max}}{\alpha \mu}= \frac{\sqrt{\beta }}{\sqrt{(1-\alpha ) \alpha }}; \\ \tau _{Max}&=&\displaystyle \frac{\sqrt{\alpha \beta }}{\sqrt{1-\alpha }} \mu .\end{array} $$
(2.8)
Fig. 2
figure 2

Scheme for the energy minimization in (2.7) for a damageable Neo-Hookean material. Here \(\alpha = 0.2; \beta = 0.5, \mu = 2\)

Full size image

We restrict here to the case of irreversible damage. Two points are worth of noticing; we are here neglecting all rate and interfacial energy effects. In other words time represents just an order parameter so that the results are rate-independent (invariant with respect to loading time rescaling) and all the solutions with fixed phase fractions are energetically equivalent. As a result, in the following we focus on the two phase solutions (see Fig. 3) with the body with \(0\leq Y\leq s\) in the undamaged regime and with \(s \leq Y\leq L\) in the damaged regime. This solution (or its symmetric) is chosen because it is the one with one single interface, so that it would be preferred if interfacial energy terms were considered.

Fig. 3
figure 3

Scheme of a two-phase damage evolution in a rectangular body under simple shear deformations. Here \(\alpha = 0.2; \beta = 0.5, \mu = 2\)

Full size image

We may then measure the damage by introducing the overall stiffness function

$$ \hat{\mu}(k_{max})=\left (1 +(\alpha -1)\frac{s (k_{max})}{L}\right ) \mu \in (\alpha \mu , \mu ), $$

with the system in the virgin undamaged configuration \(\hat{\mu}=\mu \) for \(k_{max}\leq k_{a}\) with homogeneous elastic modulus \(\mu \) and in the fully damage configuration (damage saturation limit) if \(k_{max}\geq k_{b}\) with homogeneous damaged elastic modulus \(\hat{\mu}=\alpha \mu \).

Consider first the case of monotonic loading starting from the virgin configuration \(\hat{\mu}=\mu \). According with previous description and our hypothesis of total (elastic plus damage dissipated) energy minimization (corresponding to the convex envelope of the energy \(w\)), it is easy to show that, if we define the averaged assigned shear \(\bar{k}:=\frac{\delta}{L}\), the system remains in the virgin configuration \(s=L\) and

$$ \tau (\bar{k})=\mu \bar{k} $$

until

$$ \bar{k}_{max}\leq k_{a}, $$

where by \(\bar{k}_{max}\) we indicate the maximum value of the past attained averaged shear (for monotonic loading \(\bar{k}=\bar{k}_{max}\)).

For \(\bar{k}\in (k_{a},k_{b})\) we have a two-phase configuration with

$$ s(\bar{k}_{max})= \frac{\beta \mu \sqrt{\frac{\alpha }{\beta -\alpha \beta }}}{\bar{k}_{max}} $$

and a variable overall system stiffness

$$ \hat{\mu}(\bar{k}_{max})= \sqrt{\frac{\alpha \beta}{1 -\alpha}} \frac{ \mu}{\bar{k}_{max}}.$$
(2.9)

The system follows the transition plateaux \(\tau =\tau _{Max}\) until \(\bar{k}=k_{b}\) and the system is in the fully damaged configuration. For higher values of assigned displacement the system follows the fully damaged branch

$$ \tau (\delta )=\frac{\alpha}{\mu }\bar{k}. $$

Observe that we can extend easily previous considerations also when cyclic loading paths are considered. Thus if after monotonic loading with \(\bar{k}_{max}\in (k_{a}, k_{b})\) (path O-A-B in Fig. 3) we unload, the system behaves elastically, within a damaged state and a homogenized elastic modulus (2.9) (path B-O in the figure). If reloaded the system follows again the same path until \(\tau =\tau _{Max}\) and \(\bar{k}=\bar{k}_{max}\) (point B) and the system follows again the stress plateaux. The damage is complete and \(\hat{\mu}=\alpha \mu \) when \(\bar{k}_{max}=k_{b}\). Such conclusions are easily attained if one considers that at fixed \(k_{max}\) there is a fraction \(L-s(\bar{k}_{max})\) irreversibly damaged. Thus we can minimize the energy using again the same (local) conditions in (2.5) for the only undamaged fraction.

Possible healing effects, due to molecule recrosslinking, can be considered by following the approach in [7].

2.1.2 Damageable Gent Material

Consider now the constitutive assumptions (1.4), (1.5). In this case for the considered deformation class we have an energy \(\bar{\Phi}(\boldsymbol{F})=\Phi (I_{1})=\Phi (3+k^{2})=\phi (k)\) with

$$ \phi (k)=\left \{ \textstyle\begin{array}{l@{\quad}l@{\quad}l} -\frac{\mu}{2}J_{m} \ln \left ( 1-\frac{k^{2}}{J_{m}} \right ) & \mbox{ if } & k^{max}< k^{cr}, \vspace{0.2 cm} \\ -\frac{\alpha \mu}{2}J_{m} \ln \left ( 1-\frac{k^{2}}{\alpha J_{m}} \right )+\beta \frac{\mu}{2} & \mbox{ if } & k^{max}\geq k^{cr} , \end{array}\displaystyle \right . $$
(2.10)

where the damage shear threshold is obtained by solving \(I_{cr}=3+k_{cr}^{2}\).

After imposing the boundary conditions for the one phase undamaged solution we have \(p=p_{u}=\frac{\mu J_{m}k^{2}}{J_{m}-k^{2}}\) corresponding to

$$ \boldsymbol{T}_{u}= \begin{bmatrix} \frac{\mu J_{m}k^{2}}{J_{m}-k^{2}} &\frac{\mu J_{m}k}{J_{m}-k^{2}} &0 \vspace{0.2 cm} \\ \frac{\mu J_{m}k}{J_{m}-k^{2}} &0 &0 \vspace{0.2 cm} \\ 0 &0 &0 \end{bmatrix}, \hspace{1 cm} \boldsymbol{S}_{u}= \begin{bmatrix} 0 &\frac{\mu J_{m}k}{ J_{m}-k^{2}} &0 \vspace{0.2 cm} \\ \frac{\mu J_{m}k}{ J_{m}-k^{2}} &0 &0 \vspace{0.2 cm} \\ 0 &0 &0 \end{bmatrix} . $$

When we consider the two phase solutions, in the damaged fraction we have

$$ \boldsymbol{T}_{d}= \begin{bmatrix} \frac{\mu \alpha J_{m} k^{2}}{\alpha J_{m}-k^{2}} & \frac{\mu \alpha J_{m}k}{\alpha J_{m}-k^{2}} &0 \vspace{0.2 cm} \\ \frac{\mu \alpha J_{m}k}{\alpha J_{m}-k^{2}} &0 &0 \vspace{0.2 cm} \\ 0 &0 &0 \end{bmatrix}, \hspace{1 cm} \boldsymbol{S}_{d}= \begin{bmatrix} 0 &\frac{\mu \alpha J_{m}k}{\alpha J_{m}-k^{2}} &0 \vspace{0.2 cm} \\ \frac{\mu \alpha J_{m}k}{\alpha J_{m}-k^{2}} &0 &0 \vspace{0.2 cm} \\ 0 &0 &0 \end{bmatrix} . $$

By following the same theoretical considerations of the previous case we are then lead to the variational problem in (2.7). Specifically the Maxwell rule leads to equations

$$ \textstyle\begin{array}{l} \displaystyle \int _{k_{a}}^{k_{cr}}\frac{\mu J_{m}k}{J_{m} - k^{2}}dk - \int _{k_{cr}}^{k_{b}}\frac{\alpha \mu J_{m}k}{\alpha J_{m}-k^{2}}dk= \tau _{Max} \left ( k_{b} - k_{a} \right ), \vspace{0.2 cm} \\ \displaystyle \frac{\mu J_{m} k_{a}}{J_{m}-k_{a}^{2}}= \frac{\alpha \mu J_{m} k_{b}}{\alpha J_{m}-k_{b}^{2}}=\tau _{Max}, \end{array} $$
(2.11)

that can be solved numerically to determine the (unique) solution \(k_{a}\), \(k_{b}\), \(\tau _{Max}\) (see Fig. 4).Footnote 1

Fig. 4
figure 4

Scheme for the energy minimization in (2.7) for a damageable Gent material. Here \(\alpha =2\); \(\mu = 1\), \(J_{m}=4\), \(I_{cr}=6\) so that \(\beta =8\ln \frac{5}{4}\sim 1.78\), \(\tau _{max}\sim 4.068\), \(k_{a}=\sim 1.57\), \(k_{b}=\sim 2.01\)

Full size image

Once determined these solutions we may obtain the overall behavior by following the same reasoning of the previous case with both monotonically loading (path O-A-C-D) and unloading paths (path O-A-B-O) as exhibited in Fig. 5 where dashed lines represent the (fixed damaged fraction) unloading force-displacement curves.

Fig. 5
figure 5

Overall shear stress-strain relation for a damageable Gent material under the hypothesis of homogeneous simple-shear deformations. Here \(\alpha =2; \mu = 1\), \(J_{m}=4\), \(I_{cr}=6\) so that \(\beta =8 \ln \frac{5}{4}\)

Full size image

2.2 Inhomogeneous Deformations: Antiplane Shear for a Circular Cylinder

To show the effectiveness of the proposed model also for more complex non homogeneous deformations, we consider now a long hollow cylinder \(\mathcal {C}\) (see Fig. 6), bonded to a fixed rigid support at its outer radius \(R=R_{e}\) and to a rigid sleeve at its inner radius \(R=R_{i}=\bar{\rho}R_{e}\). We assume that a force per unit cylinder length \(F=2\pi R_{i} \tau _{i} \) is applied to the inner sleeve. We introduce the cylindrical coordinates \(\Theta \in (0, 2 \pi )\), \(R\in (R_{i},R_{e})\), and \(Z\in (0,L)\) and the corresponding local frame references \(\boldsymbol{e}_{\Theta}\), \(\boldsymbol{e}_{R}\), and \(\boldsymbol{e}_{Z}\). We assume the antiplane shear class [16]

$$ \left \{ \textstyle\begin{array}{l@{\quad}l@{\quad}l} r = R=\rho R_{e}, & \rho \in (\bar{\rho}, 1), \vspace{0.2 cm} \\ \theta =\Theta , & \Theta \in (0,2 \pi ),& \vspace{0.2 cm} \\ z=Z+h(R), & Z \in (0,1), \end{array}\displaystyle \right . $$
(2.12)

where \(\rho :=\frac{R}{R_{e}}\) is the non dimensional radius and

$$ w(\rho )=\frac{h(\rho R_{e})}{R_{e}} $$
(2.13)

represents an unknown, radius dependent, non dimensional displacement function along the cylinder axes direction. We finally denote by \(\delta =\frac{d}{R_{e}}\) the non dimensional displacement of the inner sleeve, so that we have

$$ w(\bar{\rho})=\delta , \,\,\, w(1)=0. $$
(2.14)
Fig. 6
figure 6

Scheme of the case of inhomogeneous deformations for a damageable Neo-Hookean material. Here we assumed \(\alpha =0.2\), \(\beta =0.5\); \(\mu =2\), \(\bar{\rho}=0.5\). In A the system is an undamaged state, in B-C the system is in a partially damaged state (two phase solutions), in D the system is in the fully damaged state. Dashed lines correspond to unloading, where the interface between damaged and undamaged states is fixed as well as the phase fractions

Full size image

The deformation gradient is then

$$ \boldsymbol{F}=\boldsymbol{e}_{R} \otimes \boldsymbol{e}_{R}+ \boldsymbol{e}_{\Theta }\otimes \boldsymbol{e}_{\Theta }+k \boldsymbol{e}_{Z}\otimes \boldsymbol{e}_{R}+ \boldsymbol{e}_{Z} \otimes \boldsymbol{e}_{Z}, $$
(2.15)

where

$$ k=h'(R)=\dot{w} (\rho ) $$
(2.16)

is the shear. Here we denote by \((\, . \,)'=\frac{d (\, . \,)}{dR}\) and by \(\dot{(\, . \,)}= \frac{d (\, . \,)}{d\rho}\). Thus we have the right Cauchy Green tensor is

$$ \boldsymbol{B}=\boldsymbol{e}_{R} \otimes \boldsymbol{e}_{R}+ \boldsymbol{e}_{\Theta }\otimes \boldsymbol{e}_{\Theta}+ (1+k^{2}) \boldsymbol{e}_{Z} \otimes \boldsymbol{e}_{Z} +k (\boldsymbol{e}_{R} \otimes \boldsymbol{e}_{Z}+\boldsymbol{e}_{Z}\otimes \boldsymbol{e}_{R}), $$
(2.17)

with

$$ I_{1}=I_{2}=3+k^{2}. $$
(2.18)

2.2.1 Damageable Neo-Hookean Material

Under the Neo-Hookean assumption we obtain

$$ \boldsymbol{T}_{u}=\mu \boldsymbol{B}_{u}- p_{u}(r) \boldsymbol{I}, \hspace{0.5 cm} \boldsymbol{T}_{d}=\alpha \mu \boldsymbol{B}_{d}- p_{d}(r) \boldsymbol{I}, $$

where \(p_{u}\) and \(p_{d}\) are pressures taking care of the incompressibility constraint.

By imposing the boundary conditions

$$ \boldsymbol{T}_{\Theta \Theta}=\boldsymbol{T}_{RR}=0, $$
(2.19)

we obtain \(p_{u}=\mu \) and \(p_{d}=\alpha \mu \). As a result the first Piola stress tensor has the only non zero components

$$ \textstyle\begin{array}{l} \tau _{u}(\rho ):=\boldsymbol{S}_{R Z}=\boldsymbol{S}_{Z R}=\mu \, k_{u}= \mu \, \dot{w}_{u} (\rho ), \vspace{0.2 cm} \\ \tau _{d}(\rho ):=\boldsymbol{S}_{R Z}=\boldsymbol{S}_{Z R}=\alpha \mu \, k_{d}= -\alpha \mu \, \dot{w}_{d} (\rho ) \end{array} $$
(2.20)

that for a fully undamaged state delivers

$$ \frac{\partial \boldsymbol{S}_{RZ}}{\partial R}+\frac{1}{R} \frac{\partial \boldsymbol{S}_{\Theta Z}}{\partial \Theta}+\frac{\partial \boldsymbol{S}_{ZZ}}{\partial Z}+ \frac{1}{R}\boldsymbol{S}_{RZ}=0. $$
(2.21)

Consider first the one phase, fully undamaged solution. We obtain the local equilibrium equation

$$ \ddot{w}_{u}+\frac{\dot{w}_{u}}{\rho }=0. $$
(2.22)

Observe that the same equation can be obtained by considering the minimization of the Gibbs energy density

$$ \textstyle\begin{array}{lll} g&=& \displaystyle \frac{G}{2\pi R_{e}^{2}}=\int _{\bar{\rho}}^{1} \left [\frac{\mu}{2}\left ( \dot{w}_{u}(\rho ) \right )^{2} \rho \right ] \, d \rho - \tau _{i} \delta \vspace{0.2 cm} \\ &=&\displaystyle \int _{\bar{\rho}}^{1} \left [ \frac{\mu}{2}\left ( \dot{w}_{u}(\rho ) \right )^{2}\rho + \tau _{i} \dot{w}(\rho ) \right ] d \rho , \end{array} $$
(2.23)

corresponding to the Euler–Lagrange equation

$$ \dot{\overline{\rho \,\,\, {\dot{\!}\!\! w}_{u}}}=0 $$
(2.24)

and the natural boundary condition

$$ -\dot{w}_{u}(\bar{\rho})=\frac{t}{\bar{\rho}}, $$
(2.25)

where we introduced the non dimensional assigned axial force

$$ t=\frac{\tau _{i} \bar{\rho}}{\mu}. $$
(2.26)

This equation can be solved easily as

$$ w_{u}(\rho )=A-t \ln \rho . $$
(2.27)

Thus by imposing \(w_{u}(1)=0\), we obtain

$$ w_{u}(\rho )=-t \ln \rho. $$
(2.28)

We then deduce the equilibrium shear strain and stress

$$ k_{u}(\rho )=\frac{t}{ \rho}, \hspace{0.5 cm} \tau _{u}(\rho )=\frac{\mu \, t}{\rho}. $$
(2.29)

Based on the previous assumption of a two-wells energy, we consider the possibility of two phase solutions. Since the shear decreases as \(\rho \) increases, we assume that there exists \(\rho _{Max} \in (\bar{\rho},1)\) such that the material is undamaged for \(\rho \in (\rho _{Max}, 1)\) and damaged for \(\rho \in (\bar{\rho}, \rho _{Max})\). Using (2.17) and (1.1) the non-dimensional total potential energy can written as

$$ \textstyle\begin{array}{lll} g(w)&=& \displaystyle \int _{\bar{\rho}}^{\rho _{Max}} \left [\alpha \frac{\mu}{2} \dot{w}_{d}^{2}(\rho ) \rho +\tau _{i} \dot{w}_{d}(\rho ) \right ] d\rho \vspace{0.2 cm} \\ &+& \displaystyle \int _{\rho _{Max}}^{1} \left [ \frac{\mu}{2} \dot{w}_{u}^{2}(\rho ) +\tau _{i} \dot{w}_{u}(\rho )\rho \right ] d \rho , \end{array} $$
(2.30)

where in this case both the function \(w_{d}(\rho )\) in the damaged region and \(w_{u}(\rho )\) in the undamaged cylinder portion and \(\rho _{Max}\) have to be determined by energy minimization. The Euler Lagrange equations in both regions are again given by (2.24) as

$$ \dot{\overline{\rho \dot{w_{d}}}}=0, \hspace{1 cm} \dot{\overline{\rho \dot{w_{u}}}}=0, $$
(2.31)

whereas the natural conditions are

$$ \textstyle\begin{array}{l} \alpha \mu \dot{w}_{d}(\bar{\rho})=\frac{t}{\bar{\rho}} , \vspace{0.2 cm} \\ \alpha \dot{w}_{d}(\rho _{Max})=\dot{w}_{u}(\rho _{Max}), \end{array} $$
(2.32)

to which we add the continuity of the displacement at the interface and the assigned displacement at the external cylinder

$$ \textstyle\begin{array}{l} w_{d}(\rho _{Max})=w_{u}(\rho _{Max}), \vspace{0.2 cm} \\ w_{u}(1)=0. \end{array} $$
(2.33)

Observe that in this way we respect the interface condition (\(\boldsymbol{n}_{R}\) being the normal to the interface)

$$ \boldsymbol{F}_{2}-\boldsymbol{F}_{1}=\boldsymbol{a}\otimes \boldsymbol{n}_{R} =(\dot{w}_{u}(\rho _{Max})-\dot{w}_{d}(\rho _{Max})) \boldsymbol{e}_{z}\otimes \boldsymbol{e}_{R} $$
(2.34)

and the equilibrium condition at the interface

$$ \boldsymbol{S}(\boldsymbol{F}_{u})\boldsymbol{n}_{R}-\boldsymbol{S}( \boldsymbol{F}_{d})\boldsymbol{n}_{R}=(\mu \dot{w}_{u}(\rho _{Max})- \alpha \mu \dot{w}_{d}(\rho _{Max}))\boldsymbol{n}_{R}=\boldsymbol{0}. $$

The equilibrium configurations are then obtained by imposing the boundary condition (2.32) and (2.33) to the equilibrium equations (2.31)

$$ \left \{ \textstyle\begin{array}{l@{\quad}l@{\quad}l} w_{d}= \displaystyle -\frac{t}{\alpha} \ln (\rho ^{\alpha }\rho _{Max}^{ \alpha -1}),& \rho \in (\bar{\rho}, \rho _{Max}),& \mbox{damaged region,} \vspace{0.2 cm} \\ w_{u}=-\displaystyle t \ln (\rho ), & \rho \in (\rho _{Max}, 1),& \mbox{undamaged region}. \end{array}\displaystyle \right . $$
(2.35)

By following the same consideration in (2.5), that can be considered as classical Weierstrass Erdman corner conditions, corresponding to determine the convex envelope of the energy assigned by (2.8). In this case, based on the described monotonicity results of the stress, we have that the two phase solutions are geometrically characterized by a stress-dependent interface \(\rho _{Max}(t)\) such that the system is undamaged for \(\rho \in (\bar{\rho}, \rho _{Max}(t))\) and damaged for \(\rho \in (\rho _{Max}(t), 1)\). The solution can be determined by the following equilibrium condition at the interface: \(\tau _{d}(\rho _{Max})=\tau _{u}(\rho _{Max})=\tau _{max}=\sqrt{ \frac{\alpha \beta}{1-\alpha}}\mu \) that gives

$$ \rho _{Max} = \frac{\mu \, t}{\tau _{max}}. $$
(2.36)

We remark that in both undamaged and damaged phases the system respects the same variation of the shear stress

$$ \tau _{u}(\rho )=\frac{\mu \, t}{\rho}, \hspace{0.5 cm} \tau _{d}(\rho )=\frac{\alpha \mu \, t}{\rho}, $$

and shear strains

$$ \left \{ \textstyle\begin{array}{l@{\quad}l@{\quad}l} k_{d}=\displaystyle \frac{t}{\alpha \rho} & \rho \in (\bar{\rho}, \frac{\mu t}{\tau _{Max}})& \mbox{damaged region} \vspace{0.2 cm} \\ k_{u}=\displaystyle \frac{t}{\rho}, & \rho \in ( \frac{\mu \, t}{\tau _{Max}}, 1)& \mbox{undamaged region}. \end{array}\displaystyle \right . $$
(2.37)

Observe that the system is in the fully undamaged configuration if \(\rho _{Max}(t)<\bar{\rho}\), i.e. for \(t<\frac{\bar{\rho}\tau _{Max}}{\mu }\). On the other hand the system is in the fully damaged configuration when if \(\rho _{Max}(t)>1\), i.e. for \(t>\frac{\tau _{Max}}{\mu }\). Moreover in the two phase regime, when the cylinder is partially damaged, the global (average) force-displacement relation is given by

$$ \delta =\delta _{d}=-\frac{t}{\alpha } \ln (\bar{\rho}^{\alpha }\rho _{Max}^{ \alpha -1}). $$
(2.38)

By using (2.27) we may synthesize the obtained behavior by considering a monotonically increasing loading \(\tau _{i}\).

The system behaves as follows

$$ \textstyle\begin{array}{l@{\quad}l@{\quad}l} \delta _{u}=-t \ln \bar{\rho}, & t\leq \frac{\bar{\rho}}{\mu}\, \tau _{Max}& \mbox{undamaged} \vspace{0.2 cm} \\ \delta _{d}=-\frac{t}{\alpha} \ln (\bar{\rho}^{\alpha }\rho _{Max}^{ \alpha -1}), &\frac{\bar{\rho}\, \tau _{Max}}{\mu}< t< \frac{\tau _{Max}}{\mu}& \mbox{partially damaged} \vspace{0.2 cm} \\ \delta _{fd}=-\frac{t}{\alpha} \ln \bar{\rho}, & t\geq \frac{\tau _{max}}{\mu}& \mbox{fully damaged} \end{array} $$
(2.39)

The global shear stress displacement curve is represented in Fig. 6 and the corresponding elastica in Fig. 7.

Fig. 7
figure 7

Displacement fields for the same system in Fig. 6 (damageable Neo-Hookean material) for different values of the assigned force \(t\). In A the system is in the undamaged state, in B-C the system is in a partially damaged state (two phase solutions), in D the system is in the fully damaged state

Full size image

2.2.2 Damageable Gent Material

When a damageable Gent material is considered for the hollow cylindrical shape, by using (1.5), (2.19) and (2.15), (2.16) we obtain the stresses

$$ \boldsymbol{S}_{u}= \begin{bmatrix} 0 &0 &\frac{k\mu J_{m}}{ J_{m}-k^{2}} \vspace{0.2 cm} \\ 0 &0 &0 \vspace{0.2 cm} \\ \frac{k\mu J_{m}}{ J_{m}-k^{2}} &0 &0 \end{bmatrix} , \hspace{1 cm} \boldsymbol{S}_{d}= \begin{bmatrix} 0 &0 &\frac{k\mu \, \alpha \, J_{m}}{\alpha \, J_{m}-k^{2}} \vspace{0.2 cm} \\ 0 &0 &0 \vspace{0.2 cm} \\ \frac{k\mu \, \alpha \, J_{m}}{\alpha \, J_{m}-k^{2}} &0 &0 \end{bmatrix}. $$
(2.40)

The only non trivial local equilibrium equation (2.21) for the fully undamaged configuration

$$ \dot{\overline{\rho \, \tau _{u}(\dot{\!\!\!w})}}=0, \mbox{ with } \tau _{u}= \frac{-\dot{w}_{u}(\rho ) \mu J_{m}}{ J_{m}-{\dot{w}_{u}(\rho )}^{2}}, $$
(2.41)

with the boundary conditions

$$ \tau _{u}(\bar{\rho})=\tau _{i}, \hspace{1 cm} w(1)=0. $$

Thus we have

$$ \rho \, \tau _{u}(\rho )=\rho \frac{-\dot{w}_{u}(\rho ) \mu J_{m}}{ J_{m}-{\dot{w}_{u}(\rho )}^{2}}= \bar{\rho}\tau _{i}. $$
(2.42)

We then obtain the displacement field

$$ w_{u}(\rho )= \left [\rho \left (-\sqrt{ \left (\rho ^{2}+t^{2} \right )}\right )+ \left (\rho ^{2}+t^{2} \log \left ( \frac{\sqrt{t^{2}+1}+1}{\rho +\sqrt{\rho ^{2}+t^{2}}}\right )-1 \right )+ \left (t^{2}+1\right )\right ]\frac{\sqrt{J_{m}}}{2 t}, $$

where we introduced the non-dimensional applied force per unit cylinder length

$$ \bar{T}=\frac{2 \bar{\rho}\tau _{i}}{\sqrt{J_{m}}\, \mu}. $$

Correspondingly we find the shear stretch and stress

$$ k_{u}(\rho )= \frac{\sqrt{J_{m}} \left ( \sqrt{\rho ^{2}+t^{2}}-\rho \right )}{t}, \tau _{u}(\rho )=\sqrt{J_{m}} \, \mu \frac{ t}{\rho }.$$
(2.43)

This solution represents the global minimum of the energy until \(\tau (\bar{\rho})=\tau _{Max}\), where \(\tau _{Max}\) is the solution of (2.11), i.e. for

$$ t_{Max}\leq \frac{2 \bar{\rho}\tau _{max}}{\sqrt{J_{m}} \mu }, $$

where \(t_{Max}\) represents the maximum attained value of the non dimensional force \(t\).

If \(t\) is increased more, following the same discussion of previous sections, two phase solutions correspond to the global energy minimization, with an internal damaged portion of the hollow cylinder assigned by

$$ \rho \in (\bar{\rho}, \rho _{Max}), $$

where, by using (2.43)b, we obtain

$$ \rho _{Max}=\sqrt{J_{m}}\, \mu \frac{ t}{\tau _{Max} }. $$

In this case the displacement \(w_{d}(\rho )\) of the damaged fraction can be again determined by the Euler-Lagrange equilibrium condition (2.41) with the continuity boundary condition \(w_{d}(\rho _{Max})=w_{u}(\rho _{Max})\),

$$ \textstyle\begin{array}{lll} w_{d}(\rho )&=&\left [(\alpha -1) \left (-J_{m}^{3/2}\right ) \mu ^{2} t^{2}+4 \sqrt{J_{m}} \tau _{Max}^{2} \left (\alpha \rho ^{2}-1\right ) \right . \vspace{0.2 cm} \\ &+& \left . 4 \tau _{Max}^{2} \left (\rho ^{3} \left (-\sqrt{ \frac{\alpha ^{3} J_{m}}{\alpha \rho ^{2}+t^{2}}}\right )-\rho t^{2} \sqrt{ \frac{\alpha J_{m}}{\alpha \rho ^{2}+t^{2}}}+\sqrt{J_{m} \left (t^{2}+1 \right )}\right ) \right . \vspace{0.2 cm} \\ &+& \left . 4 \sqrt{J_{m}} t^{2}\tau _{Max}^{2} \log \frac{\left (\sqrt{t^{2}+1}+1\right ) \left (\sqrt{\alpha \left (\alpha J_{m} \mu ^{2}+4 \tau _{Max}^{2}\right )}+\alpha \sqrt{J_{m}} \mu \right )}{\left (\sqrt{J_{m} \mu ^{2}+4 \tau _{Max}^{2}}+\sqrt{J_{m}} \mu \right ) \left (\alpha \rho +\sqrt{\alpha \left (\alpha \rho ^{2}+t^{2}\right )}\right )} \right . \vspace{0.2 cm} \\ &+& \left . J_{m} \mu t^{2} \left (\sqrt{\alpha \left (\alpha J_{m} \mu ^{2}+4 \tau _{Max}^{2}\right )}-\sqrt{J_{m} \mu ^{2}+4\tau _{Max}^{2}} \right )\right ]\frac{1}{8 t\tau _{Max}^{2}}.\end{array} $$

Correspondingly the shear stress and strain are given by

$$ \textstyle\begin{array}{l} k_{d}(\rho )= \frac{3 \rho ^{2} t^{2} \sqrt{\frac{\alpha ^{3} J_{m}}{\alpha \rho ^{2}+t^{2}}}+2 \alpha \rho ^{4} \sqrt{\frac{\alpha ^{3} J_{m}}{\alpha \rho ^{2}+t^{2}}}+\sqrt{J_{m}} \left (t^{2} \left (\sqrt{\alpha \left (\alpha \rho ^{2}+t^{2}\right )}-2 \alpha \rho \right )-2 \alpha ^{2} \rho ^{3}\right )+t^{4} \sqrt{\frac{\alpha J_{m}}{\alpha \rho ^{2}+t^{2}}}}{2 \left (t^{3}+\alpha \rho ^{2} t\right )}, \vspace{0.2 cm} \\ \tau _{d}(\rho )=\sqrt{J_{m}}\, \mu \frac{ t}{\rho }. \end{array} $$
(2.44)

Two phase solutions with \(w=w_{d}(\rho )\), \(\rho \in (\bar{\rho}, \rho _{Max})\) and \(w=w_{u}(\rho )\), \(\rho \in (\rho _{Max},1)\) represent the global minima until the interface attains the external boundary \(\rho _{Max}=1\). Thus in this case two-phase solutions are global minimizers for

$$ \frac{2 \bar{\rho}\, \tau _{max}}{\sqrt{J_{m}} \mu }\leq t_{Max}\leq \frac{2 \tau _{max}}{\sqrt{J_{m}}\mu}. $$

Finally for larger maximum attained forces the system is fully damaged with a displacement field

$$ \textstyle\begin{array}{lll} w_{fd}(\rho )&=&\left [ \frac{\rho \left (\rho \sqrt{\alpha ^{3} J_{m} \left (\alpha \rho ^{2}+t^{2}\right )}-\alpha \sqrt{J_{m}} \left (\alpha \rho ^{2}+t^{2}\right )\right )-t^{2} \sqrt{\alpha J_{m} \left (\alpha \rho ^{2}+t^{2}\right )} \log \left (\alpha \rho +\sqrt{\alpha \left (\alpha \rho ^{2}+t^{2}\right )}\right )}{\sqrt{\alpha \left (\alpha \rho ^{2}+t^{2}\right )}} \right . \vspace{0.1 cm} \\ &+& \left . \sqrt{J_{m}} \left (-\alpha +\sqrt{\alpha \left (\alpha +t^{2} \right )}+t^{2} \log \left (\alpha +\sqrt{\alpha \left (\alpha +t^{2} \right )}\right )\right ) \right ] \frac{1}{2 t}\end{array} $$

corresponding to shear stress and strain fields

$$ \textstyle\begin{array}{lll} k_{fd}(\rho )&=& \frac{\alpha ^{2} J_{m}^{3/2} \left (\alpha \rho ^{2}+t^{2}\right )^{2}-\rho \left (\rho ^{2} \sqrt{\alpha ^{7} J_{m}^{3} \left (\alpha \rho ^{2}+t^{2}\right )}+t^{2} \sqrt{\alpha ^{5} J_{m}^{3} \left (\alpha \rho ^{2}+t^{2}\right )}\right )}{J_{m} t \left (\alpha \left (\alpha \rho ^{2}+t^{2}\right )\right )^{3/2}} \vspace{0.2 cm} \\ \tau _{fd}(\rho )&=&\sqrt{J_{m}}\, \mu \frac{ t}{\rho }. \end{array} $$
(2.45)

The damage evolution, the elastica, and the force-displacement diagrams are reported in Fig. 8.

Fig. 8
figure 8

Overall force-displacement curves for inhomogeneous deformations with a damageable Gent material. Here we assumed \(\alpha =2\), \(I_{cr}=6\); \(\mu =1\), \(\bar{\rho}=0.5\), \(J_{m}=4\). Left figure force-displacement curves (dashed lines correspond to unloading). Right figure displacement fields for A an undamaged state, B-C partially damaged states (two phase solutions), and a fully damaged state D

Full size image

3 Concluding Remarks

We have shown the possibility of describing damage effects as a material phase transition phenomenon. Specifically we considered damageable versions of the Neo-Hookean and Gent constitutive assumptions for rubberlike incompressible materials. Thus the material is characterized by two states, damaged and undamaged. The damaged state is characterized by a lower shear modulus, interpreted as a fraction of broken molecules per unit volume at the network scale. We then assign, under an isotropic assumption of damage, a critical damage threshold, \(I_{1}=I_{1}^{cr}\), measuring a critical value of the averaged molecular stretch at each material point, inducing material damage. The hard→soft (undamaged→damaged) transition corresponds then to a finite energy dissipation, proportional to the broken chains’ fraction, that we deduce as a function of the threshold \(I_{1}^{cr}\). Many possible extension to more complex and realistic constitutive laws, possibly with more than one damaged state with progressive damage, can be considered. The choice of this paper has been guided from the aim of attaining fully analytical solutions, allowing for a clear physical interpretation of the results. Specifically we considered a Griffith type approach minimizing the global (elastic plus dissipated) energy, thus ending up to a variational problem with non (rank-one) convex energy that can be solved following the approaches suggested in many papers of Roger Fosdick to whom this paper is dedicated. Despite the simplicity of our considered constitutive assumptions, the model qualitatively reproduce well known effects in 3-D non linear damage mechanics that we show by considering the simple case of homogeneous simple shear and inhomogeneous antiplane shear deformation classes. We believe this paper can represent a ‘proof of concept’ on the possibility of efficiently study damage problems in the field of Griffith type total energy minimization models with non-convex energies.