Viscous and Failure Mechanisms in Polymer Networks: A Theoretical Micromechanical Approach

Polymeric materials typically present a complex response to mechanical actions; in fact, their behavior is often characterized by viscous time-dependent phenomena due to the network rearrangement and damage induced by chains’ bond scission, chains sliding, chains uncoiling, etc. A simple yet reliable model—possibly formulated on the basis of few physically-based parameters—accounting for the main micro-scale micromechanisms taking place in such a class of materials is required to properly describe their response. In the present paper, we propose a theoretical micromechanical approach rooted in the network’s chains statistics which allows us to account for the time-dependent response and for the chains failure of polymer networks through a micromechanics formulation. The model is up-scaled to the mesoscale level by integrating the main field quantities over the so-called ‘chains configuration space’. After presenting the relevant theory, its reliability is verified through the analysis of some representative tests, and some final considerations are drawn.


Introduction
Polymers, thanks to their ability to withstand very large deformations prior to failure, in addition to several other useful functionalities, are widely used in various applications including biomechanics, soft actuators, stretchable electronics, and adhesives. However, under mechanical loads, polymers are also prone to failure and show a strain rate dependent response, drawbacks that can limit their applicability. The above mentioned macroscopically observable mechanisms can be conveniently studied and interpreted on the basis of the polymer's underlying molecular network, made of cross-linked chains whose mechanical behavior is well described by the so-called entropic-based elasticity. In fact, the disorder degree of the network chains can be easily related to the material entropy that, upon stretching, reduces because of the greater ordered conformation assumed by the deformed chains [1,2].
The typical polymer's underneath network microstructure is shared also by others classes of polymer-like materials, such as biological tissues and natural matters; natural biopolymers, such as DNA, proteins, cellulose, pectin, and soft tissues, as well as synthetic polymers, are obtained via the polymerization of many small molecules, known as monomers, that usually provides the material characteristic physical properties such as high deformability, good toughness, and a viscoelastic response upon mechanical stress. As an example, pectin gel is made by a chains network constituted by segments join together by crystallization to form a three-dimensional network in which water, sugar, and other materials are held [3].
Very promising theoretical models, suitable to describe the mechanics of polymers based on the so-called statistical mechanics approach, have been proposed so far [4][5][6][7][8]. The above-mentioned models where ψ and r are the free energy of the chain and the unit vector of r, respectively. Since the state of the polymer is completely defined through the conformation of its chains, it is convenient to introduce a scalar function ρ(r) providing the probability density of the chains distribution, i.e., the chain density (number of chains per unit volume of material) having a given end-to-end vector r (Figure 1). The above mentioned function can be identified as the Chain Configuration Density Function (CCDF), whose physical dimension is [L −6 ] (number of chains per unit volume of the configuration space and per unit volume of the physical space). It is convenient to write the CCDF function as: ρ(r) = c a · ϕ(r) where c a is the chain density in the physical space (density of cross-linked chains), while ϕ(r) is the distribution function of r in the configuration space (accordingly, ϕ 0 (r) is the corresponding function in the undeformed state). It is worth noting that, since we are interested in the mechanical response of the polymer, only the bearing (active) chains have to be considered, i.e., those connected at both their extremities to others chains, while the free or dangling chains are not considered; in fact, because of the lack of connection to the network, their contribution to the material's bearing mechanism is negligible. In this perspective, c a represents the density of the active chains only. The integral of ρ(r) over the configuration space provides the actual chains density, i.e., where the last expression corresponds to the integration over the whole configuration space in spherical coordinates (r, α, β). It is worth noting that, being ϕ a probability density function, it must fulfill the condition ϕ = 1 [11].
For sake of exemplification, the simplest 3-chains model assumes the material to be made of several repeated identical unit cells, each one containing three freely jointed chains oriented along the three Cartesian axis; by assuming the affine deformation hypothesis, the free energy per unit volume of the material is given by = ( + + − 3) (with , , being the number of chains per unit volume, the Boltzmann constant and the absolute temperature, respectively, while , = 1,2,3 are the macroscopic stretches existing along the Cartesian axes) that recovers the wellknown neo-Hookean model according to which = ( − 3), where = is the shear moduls of the material, and is the trace of the Cauchy deformation tensor. More complex models (8chains, full network, etc.) are simply based on more general similar assumptions.

The Chain Configuration Space and the Chain Configuration Density Function (CCDF)
By adopting the so-called freely jointed chain model (FJC), the state of a chain depends only on its end-to-end vector, i.e., on the vector linking the chain's ends. In the reference (stress-free) configuration state Ω , the end-to-end chains' vector is identified by , while in a generic deformed state Ω such a vector is stretched to the new conformation given by = , being the deformation gradient tensor describing the macroscopic deformation of the material. The above undeformed-deformed end-to-end vector relationship has been written according to the so-called affine deformation hypothesis [21]; according to such an assumption the deformation of the polymer chains is equal to the deformation of the continuum medium in which the chains are embedded. According to the freely joined chain model, the force experienced by a single chain depends only on its end-to-end vector, so we can associate to each point of the configuration space Ω = { | ∈ ℝ } to the corresponding chain force, i.e., : Ω → ℝ , with where and are the free energy of the chain and the unit vector of , respectively.
Once the chains distribution is known, the energy per unit volume of material can be easily obtained as: The above energy is greater than zero in the initial stress-free state, i.e., Ψ 0 > 0 since it is ρ(r) ≥ 0 ∀r ∈ Ω; the deformation energy Ψ associated to the deformed material has to be written as: According to the above definition, the deformation energy Ψ becomes zero in the undeformed state, i.e., Ψ(F = I) = 0 and ρ(F = I) = ρ 0 .
From the above discussion, it appears that the mechanical state of the material is completely known once the function ρ(r) is provided, so the knowledge of its evolution-accounting for all the phenomena occurring in the material (deformation, chains sliding, etc.)-is sufficient to univocally determine the mechanical state of the material.

Modelling of the Visco-Elastic and the Damage Mechanisms of Polymers
In the present section, we develop a micromechanics-based model aimed at describing the visco-elastic and the damage response of polymers upon mechanical loading. Being the state of the material fully described once the Chain Configuration Density Function is known, its evolution, due to the various micromechanics phenomena involved, allows us to fully determine the mechanical state of the polymer, such as the stress and the damage level.

Microscale Approach: Evolution of the Chain Configuration Density Function
The evolution of the Chain Configuration Density Function (CCDF) in time (hereafter the time can be assumed to represent both the physical time or simply an increasing dimensionless quantity connected to the deformation) provides all the information needed to univocally define the state of the material.
The time derivative of the CCDF function can be decomposed as follows: where . ρ L indicates the variation of ρ induced by the deformation (here L = . F F −1 is the velocity deformation gradient), .
ρ v is the contribution due to the viscous phenomena, and . ρ f provides the contribution of the chains failure to be accounted for if the material's damage is assumed to occur. In the following, the three above contributions to the evolution of the CCDF function are singularly determined.
The first contribution can be evaluated by considering the material behaving elastically; this implies that the number of connected (or active) chains in Ω does not change in time, since no chains are lost (damage) or gained (self-healing) during the deformation process. As a consequence, the time derivative of the chain concentration must vanish at any time, i.e., d dt Ω ρ dΩ = 0; consequently, it can be shown that the CCDF rate, due to the deformation of the material, in the elastic regime is expressed by [11].
where . r i = L ij r j and . r i,i = L ii . It's worth noting that the latter term in (7) (ρL ii ) is zero for an incompressible material since the divergence of . r vanishing for an isochoric deformation (i.e., trL = L ii = 0), so Equation (7) becomes Figure 2).
where ̇ indicates the variation of induced by the deformation (here = ̇− 1 is the velocity deformation gradient), ̇ is the contribution due to the viscous phenomena, and ̇ provides the contribution of the chains failure to be accounted for if the material's damage is assumed to occur. In the following, the three above contributions to the evolution of the CCDF function are singularly determined.
The first contribution can be evaluated by considering the material behaving elastically; this implies that the number of connected (or active) chains in Ω does not change in time, since no chains are lost (damage) or gained (self-healing) during the deformation process. As a consequence, the time derivative of the chain concentration must vanish at any time, i.e., ∫ Ω = 0 Ω ; consequently, it can be shown that the CCDF rate, due to the deformation of the material, in the elastic regime is expressed by [11].
where ̇= and ̇, = . It's worth noting that the latter term in (7) ( ) is zero for an incompressible material since the divergence of ̇ vanishing for an isochoric deformation (i.e., = = 0 ), so Equation (7)   The second term in (6) refers to the time-dependent response of the material, so it must account for the internal microscopic rearrangement responsible for the viscous behavior. A simple model for such an internal material's reorganization can be formulated by considering the network nature of a polymer: due to the fluctuating energy at the material microscale, the chains' bonds can be assumed to detach from their stretched state and to reform again in the stress-free state of the material, i.e., after detaching they start being available again to contribute to the mechanical response [31,32]. From a statistical viewpoint, after the detachment, the chains reattach in the configuration 0 ( ) , corresponding to the initial macroscopically undeformed state of the material. In other words, the current number of cross-linked (i.e., active) chains arises from the dynamic balance between the attachment and the detachment process; this balance can be quantified once the detachment ( ) and attachment ( ) rates of the material are known [33]. These rates, in the simplest case, can be considered to be material's parameter, even if it has been observed that they should depend on the intensity of the chain force; in particular, the force effect on can be described through the Kramer's reaction rate theory, predicting an increasing of the detachment rate with the force increasing [34]. The second term in (6) refers to the time-dependent response of the material, so it must account for the internal microscopic rearrangement responsible for the viscous behavior. A simple model for such an internal material's reorganization can be formulated by considering the network nature of a polymer: due to the fluctuating energy at the material microscale, the chains' bonds can be assumed to detach from their stretched state and to reform again in the stress-free state of the material, i.e., after detaching they start being available again to contribute to the mechanical response [31,32]. From a statistical viewpoint, after the detachment, the chains reattach in the configuration ϕ 0 (r), corresponding to the initial macroscopically undeformed state of the material. In other words, the current number of cross-linked (i.e., active) chains arises from the dynamic balance between the attachment and the detachment process; this balance can be quantified once the detachment (k d ) and attachment (k a ) rates of the material are known [33]. These rates, in the simplest case, can be considered to be material's parameter, even if it has been observed that they should depend on the intensity of the chain force; in particular, the force effect on k d can be described through the Kramer's reaction rate theory, predicting an increasing of the detachment rate with the force increasing [34].
Once the polymerization of the material had occurred, the detachment/reattachment phenomenon reaches the steady state after a sufficiently long time and, in absence of damage or healing, the number of active chains does not change any more. The small deformation elastic shear modulus of the material can thus be associated to the steady state number of active chains. According to the 3-chains model, for instance, the chain concentration-shear modulus relationship is c µ = µ k B T , with c µ < c max , c max being the maximum potential value of active chains, i.e., the sum of all the attached and the detached ones at a given time instant. The evolution equation for c µ can be written as [11,12]: where the difference c max − c a (t) represents the actual concentration of detached chains, while k a and k d are the activation and deactivation cross-link rates of the polymer chains, respectively. It can be easily verified that the steady state solution of Equation (8) is given by: c a ( t → ∞ ) = k a k a +k d c max = c µ , where the last equality is justified because the shear modulus is assumed to be measured in the steady-state condition. By considering that the active chains detach from the actual (deformed) distribution ϕ(r) and the ones that become active begin to be stretched starting from the initial (reference) one, ϕ 0 (r), the density distribution rate due to the material microstructure rearrangement is expressed as: while the corresponding CCDF evolution, evaluated at constant applied deformation, assumes the form: .
where chains conservation at the steady state, dc a /dt = 0, has been leveraged. It is reasonable to assume that, when the polymer is formed, i.e., when the polymerization process is complete, the concentration of active chains corresponds to that in the steady state situation, so c a (t = 0) = c a0 = c µ . Finally, it is worth mentioning that the above mentioned parameters k a and k d can be easily related to the well-known loss and storage modulus, typically used in rheology to characterize the material's viscous behavior [11].
Finally, we consider the CCDF evolution due to the chain scission. The loss of chains entails a modification of the previous chains conservation equation, i.e., now dc a /dt ≤ 0; the general expression of the CCDF evolution provided by Equation (6) must now be updated accordingly: where the second term on the right hand side represents the rate of loss chains in the configuration space. The explicit expression for . ρ f can be obtained by introducing the failure rate function ω f (r) that provides the fraction of broken chains per unit time for a given chain's end-to-end vector r; the CCDF rate due to the chains failure is then This indicates that the number of chains lost per unit time is proportional through the function ω f ; the actual number of active chains for a given end-to-end vector r is quantified by the CCDF ρ(r). The above expression implies a loss of material or, equivalently, of active chains; this entails that the material undergoes a permanent damage because of the chain loss; it is thus possible to define a scalar damage parameter of the material at the current time instant to be defined as: It happens to be D(t) = 0 for an undamaged polymer (i.e., ϕ(t) = 1), while D(t) = 1 for a fully damaged material; the latter case corresponds to c a (t) = c µ ϕ(t) = 0 because no active chains exist anymore if the material is fully broken.
The failure rate function ω f (r) is the key factor to quantify the phenomenon of permanent chain loss; in this context, the chains' bond failure represents the physical phenomenon responsible for the chains break. It is reasonable to describe such an irreversible process by adopting the well-known Eyring's reaction rate theory [35]. For the family of polymer's chains having the end-to-end vector r, the failure rate can be expressed by [36]: In the above expression τ represents a characteristic time, while w b is the chain's segments deformation energy available to induce the chain's bond failure. The function 0 ≤ G(w(r)) ≤ 1, having the unit of t −1 , provides the probability of failure per single chain's segment and per unit time, i.e., the fraction (over the existing chains of a given length) of chain loss (or chain scission) per unit time. By introducing the chain's bond strength energy w, the above equation can be rewritten as follows: where A is a model parameter, w b = 1 2 E b ln 2 λ b is the (enthalpic) deformation energy of one chain's Kuhn segment [37] (E b , λ b being the stiffness and the stretch of the Kuhn segments, respectively, where λ b can be determined from the stationarity condition of the sum of the entropic and the enthalpic parts of the chain energy of the material, [9]), h is the Planck's constant, and γ < 0 is a parameter whose value governs the sharpness of the transition between unfailed (w b < w) and failed state (w b ≥ w) of the chains. The Heaviside step function H (r · . r) has been introduced in (14) in order to distinguish between the state of loading or unloading for a single chain; during loading the chain increases its length (r · . r > 0), and the damage can take place in the material (H = 1). During unloading the chain shortens (r · . r ≤ 0), and its damage cannot increase (H = 0). The graphical representation in the 2-D space r x , r y (the dimensionless quantities r have been used for sake of simplicity) of Equation (15) is provided in Figure 3a,c for a deformation process taking place in plane stress condition characterized by λ y = λ z = λ −1/2 x , where the failure rate ω f is illustrated (Figure 3b unloading the chain shortens ( ⋅̇≤ 0), and its damage cannot increase (ℋ = 0).
The graphical representation in the 2-D space , (the dimensionless quantities ′ = √ , ′ = /( √ ) have been used for sake of simplicity) of Equation (15) is provided in Figure 3a,c for a deformation process taking place in plane stress condition characterized by where the failure rate is illustrated (Figure 3b,d) for a polymer for two assumed values of the bond energy .  The CCDF rate due to the chains failure is depicted in Figure 3b,d; the major portion of the lost chains (see peaks of Figure 3b,d) corresponds to the chains that lie mainly aligned along the stretch direction and has a sufficiently large end-to-end vector length, i.e., for those chains with a deformation energy closed to the bonding one, ≥ ̅. Shorter chains do not fail because they are not enough stretched, while chains longer than the failed ones do not practically exist in the network; it's worth noting that the chains lying along directions mainly elongated in the -direction do not fail either because they are subjected to unloading for the adopted deformation process.
An approach for polymer damage based on the probability of failure in the chain configuration space has been recently proposed in [38]; the present approach differs from the above mentioned one being based on the kinetic of the chain failure through the reaction rate theory, whose main parameter is the chain dissociation energy. On the other hand, the statistical damage approach in [38] assumes an initial failure distribution function in the chains' configuration space and makes it to evolve according to the damage occurred in the material. The problem of polymer's chains failure ahead of a crack in a soft rubbery material has been considered by Hui et al. [39]: they adopted a cohesive model characterized by a failure force whose value depends on the thermal state of the material and on the rate at which the force is transmitted to the bond (thermo-mechanically activated bond dissociation kinetics).

Stress State in the Polymer
As mentioned before, the knowledge of the distribution function enables to completely know the mechanical state of the material; in particular, the stress state can be easily obtained as ( ) = ∂ ( ) + ( ) − , with and being the first Piola stress tensor and the hydrostatic pressure enforcing the isochoric deformation, respectively, and = det being the volumetric deformation.
On the other hand, the Cauchy stress is provided by the relation = −1 and can be demonstrated to be obtainable from the relationship [40]: The CCDF rate ω f ρ due to the chains failure is depicted in Figure 3b,d; the major portion of the lost chains (see peaks of Figure 3b,d) corresponds to the chains that lie mainly aligned along the stretch direction x and has a sufficiently large end-to-end vector length, i.e., for those chains with a deformation energy closed to the bonding one, w b ≥ w. Shorter chains do not fail because they are not enough stretched, while chains longer than the failed ones do not practically exist in the network; it's worth noting that the chains lying along directions mainly elongated in the y-direction do not fail either because they are subjected to unloading for the adopted deformation process.
An approach for polymer damage based on the probability of failure in the chain configuration space has been recently proposed in [38]; the present approach differs from the above mentioned one being based on the kinetic of the chain failure through the reaction rate theory, whose main parameter is the chain dissociation energy. On the other hand, the statistical damage approach in [38] assumes an initial failure distribution function in the chains' configuration space and makes it to evolve according to the damage occurred in the material. The problem of polymer's chains failure ahead of a crack in a soft rubbery material has been considered by Hui et al. [39]: they adopted a cohesive model characterized by a failure force whose value depends on the thermal state of the material and on the rate at which the force is transmitted to the bond (thermo-mechanically activated bond dissociation kinetics).

Stress State in the Polymer
As mentioned before, the knowledge of the distribution function ϕ enables to completely know the mechanical state of the material; in particular, the stress state can be easily obtained as ∂F + p(t)JF −T , with P and p being the first Piola stress tensor and the hydrostatic pressure enforcing the isochoric deformation, respectively, and J = detF being the volumetric deformation. On the other hand, the Cauchy stress is provided by the relation σ = J −1 PF T and can be demonstrated to be obtainable from the relationship [40]: where it has been assumed J = detF = 1 (due to the incompressibility constraint), while f is the force existing in the chains of a given end-to-end distance r (evaluable by a suitable model such as through the well-known Langevin statistics, f = ∂ψ(t) ). It's worth noting that the chain density distribution function in the reference configuration ρ 0 (t) has to be updated in the case the damage occurs in the material, and so the time dependence of ρ 0 has been explicitly indicated in (16). In fact, in this case the stress-free chains' configuration corresponds to the current chain density distribution function ρ(t) brought back to the un-deformed state, i.e., for F = I. In other words, the current distribution function in the reference state, ϕ 0 (t), must correspond to ϕ(t), being ϕ 0 (t) = ϕ(t) ≤ ϕ(t = 0) = 1 due to the chains loss; correspondingly, the density of the actual active chains is now c a (t) ≤ c a (t = 0). The pull-back operation is provided by the following relation i.e., at the current time instant t the value of the reference distribution function ϕ 0 for a given r must correspond to the value assumed by ϕ for the end-to-end vector Fr. The above described pull back operation is unnecessary if the damage is not occurring in the material, and in these cases it is simply ϕ 0 (r, t) = ϕ 0 (r, t = 0).

Polymers with Multiple Networks
The above presented micromechanical model has been developed under the hypothesis that the polymer is made of a single network, i.e., all the chains are made of a fixed number N of Kuhn's segments. However, in real polymers the presence of chains composed by different numbers of segments can be found (these polymers also termed as polydisperse networks or networks with a non-uniform weight distribution) [41,42]. Within this context, the development of polymers with double networks has attracted a lot of efforts in recent years because of their capability to dissipate energy; thanks to the sacrificial role that can be attributed to the more brittle network, the more flexible ones can act as a bridging skeleton, enhancing the fracture resistance of the polymer [42,43]. Theoretical approaches have been also developed to describe the mechanics of polydisperse networks [44,45].
The proposed model can be easily extended to the case of polydisperse networks by introducing the distribution function of the Kuhn's segments number, q(N); accordingly, the energy per unit volume (5) has to be updated as follows: where N min and N max represent the number of segments of the shortest and of the longest network in the polymer, respectively; in the case of a finite number m of different networks, the above expression must be replaced by a summation over all the networks, being N m N i =N 1 q(N i ) = 1, while q(N i ) represents the volume fracion of the i-th network. In this case, the evolution of the CCDF (Equation (6)) has to be evaluated with reference to each single network, i.e., for a specific value of N min ≤ N ≤ N max ,  (N, r), where the contributions of the deformation (7), of the viscous mechanism (10), and of the chains failure (12) have to be evaluated according to the deformation and the stress state acting on each single network.

Thermodynamics of Polymers Undergoing Chains Failure and Bond Exchange
From an energetic viewpoint, the problem under consideration must comply the first principle of thermodynamics; it implies that the material derivative of the internal energy density E n (per unit current volume) can be expressed by: where σ is the Cauchy stress tensor, q the heat flux in a given point, and s is the rate of external heat supply per unit current volume. On the other hand, the Clausius-Duhem inequality provides the second principle of thermodynamics: with ϑ being the entropy per unit current volume and T the absolute temperature. By introducing the Helmholtz free energy per unit current volume, Ξ = E n − Tϑ, the second principle becomes: For sake of convenience, let us assume an incompressible polymer, and postulate that the deformation process takes place isothermally ( . T = 0) and without any external heat supplied (q = 0, s = 0), so σ : L − . Ξ ≥ 0. In this particular case, the rate of change of the Helmholtz free energy per unit current volume becomes: where D represents the dissipated energy. By replacing (7), (10), and (12) in (22), we can recognize the energy loss to be expressed by: i.e., it is provided by the contributions coming from the time-dependent response of the material and the one from the chains failure. Being the viscous and the damage contributions independent of each other, at any time instant t, the positiveness of the dissipated energy is ensured if both the following conditions hold: The first condition is fulfilled if Ω [(ϕ(t) − ϕ 0 )]ψ dΩ ≤ 0; it corresponds to the non-increasing character of the energy stored in the polymer, i.e., the inequality Ω ϕ(t)ψ dΩ ≤ Ω ϕ 0 ψ dΩ guarantees that, because of the network relaxation, the chains detachment is associated with a decrease in the stored elastic energy and that this loss is proportional to the rate of dissociation k d and the stored elastic energy, Ω [c µ (ϕ(t) − ϕ 0 )]ψ dΩ.
The second requirement in (25) is consistent with (15), since ω f , ρ, ψ are all positive functions ∀r, while, according to (14), H (r · . r) is positive when the chain undergoes loading, r · . r > 0 (i.e., increase of its end-to-end vector); meanwhile, it becomes zero in the case of neutral loading, r · . r = 0, or unloading, i.e., r · . r < 0 (no change or decrease of its end-to-end vector).

Application of the Micromechanical Model
In the following some parametric tests as well as simulations of real polymer tests are presented in order to underline the main features of the proposed model for the viscous and damage response of polymers.

Parametric Analyses
In the present section, we present some parametric tests aimed at underlying the role played by the main parameters involved in the above presented time-dependent and damage model to quantify the response of polymers under mechanical actions. Firstly, the rate dependent response is analyzed.

Rate-Dependent Response
In this case, the time-dependent response, in absence of any damage phenomena, is modelled through the chains attachment-detachment mechanism described in Section 3.1, leading to the evolution in time of the network microstructure (see Equation (10)).
In particular, we consider the analysis of a polymer under two stretching cycles (in each one the stretch starts from the un-deformed state, λ = 1 to the deformed one, λ = 2, and, afterwards, the stretch goes back to λ = 1). The material is assumed to be incompressible and characterized by an elastic and a shear modulus, E = 5 MPa, µ = 1.67 MPa, T = 300 K, and N = 50, while the rate parameters k a and k d are made to vary as follows: k a = 1 Hz (kept constant), (i) k d = 0 Hz, (ii) k d = 0.1 Hz, and (iii) k d = 0.3 Hz. Finally, two strain rates have been adopted, i.e., the slowest case with . λ = 0.2 s −1 and the fastest one characterized by . λ = 2 s −1 . In Figure 4, the results obtained for the case of slow strain rate, . λ = 0.2 s −1 , are reported in terms of dimensionless true stress vs deformation (Figure 4a) and vs time (Figure 4b). Greater k d values lead to a greater dissipation (wider hysteresis loops), while, correspondingly, the maximum stress decreases and the stress state after unloading results to be negative; this latter observation indicates that the material has to be compressed in order to be brought to its original un-deformed state, because of the internal resetting occurred during the deformation process. Such an internal resetting is much more pronounced for higher values of k d . In this case, the time-dependent response, in absence of any damage phenomena, is modelled through the chains attachment-detachment mechanism described in section 3.1, leading to the evolution in time of the network microstructure (see Equation (10)).
In particular, we consider the analysis of a polymer under two stretching cycles (in each one the stretch starts from the un-deformed state, = 1 to the deformed one, = 2, and, afterwards, the stretch goes back to = 1). The material is assumed to be incompressible and characterized by an elastic and a shear modulus, = 5 MPa, = 1.67 MPa , = 300 K , and = 50, while the rate parameters and are made to vary as follows: = 1 Hz (kept constant), i) = 0 Hz, ii) = 0.1 Hz, and iii) = 0.3 Hz. Finally, two strain rates have been adopted, i.e., the slowest case with = 0.2s and the fastest one characterized by = 2s . In Figure 4, the results obtained for the case of slow strain rate, = 0.2 s , are reported in terms of dimensionless true stress vs deformation (Figure 4a) and vs time (Figure 4b). Greater values lead to a greater dissipation (wider hysteresis loops), while, correspondingly, the maximum stress decreases and the stress state after unloading results to be negative; this latter observation indicates that the material has to be compressed in order to be brought to its original un-deformed state, because of the internal resetting occurred during the deformation process. Such an internal resetting is much more pronounced for higher values of .
By adopting a faster strain rate value, = 2s , the material responds as indicated in Figure 5; the strain rate is greater than in the first examined case, for the same rate parameters , the  By adopting a faster strain rate value, . λ = 2 s −1 , the material responds as indicated in Figure 5; the strain rate is greater than in the first examined case, for the same rate parameters k a , k d the material has not enough time to reset internally its network, so the behavior is more close to the elastic one. because of the internal resetting occurred during the deformation process. Such an internal resetting is much more pronounced for higher values of .
By adopting a faster strain rate value, = 2s , the material responds as indicated in Figure 5; the strain rate is greater than in the first examined case, for the same rate parameters , the material has not enough time to reset internally its network, so the behavior is more close to the elastic one. In the second parametric test, the material is stretched to a given value, and then the deformation is kept constant (see details in Figure 6); the applied deformation is quantified through the parameter 0 ≤ α(t) ≤ 1.0, being λ(t) = λ 0 · α(t) with λ 0 = 2. It can be observed that the material behaves perfectly linearly when no chains detachment occurs (k d = 0); meanwhile, by increasing the detachment frequency, the stress reached in the material for a given stretch decreases ( Figure 6). Moreover, by keeping the deformation constant, the material relaxes itself more significantly for higher k d values; finally, the elastic case is recovered in the particular case k d = 0. In the second parametric test, the material is stretched to a given value, and then the deformation is kept constant (see details in Figure 6); the applied deformation is quantified through the parameter 0 ≤ ( ) ≤ 1.0, being ( ) = ⋅ ( ) with = 2. It can be observed that the material behaves perfectly linearly when no chains detachment occurs ( = 0 ); meanwhile, by increasing the detachment frequency, the stress reached in the material for a given stretch decreases ( Figure 6). Moreover, by keeping the deformation constant, the material relaxes itself more significantly for higher values; finally, the elastic case is recovered in the particular case = 0. During the loading ramp the stretch rate has been assumed to be equal to 0.2 s -1 (a) and to 2 s -1 (b).

Mechanical Response of a Single Network Polymer in Presence of Damage
In the second case, we analyze the response of a single network polymer by accounting for the occurrence of damage without any internal rearrangement features ( = 0 Hz). The CCDF time evolution, , is here evaluated according to Equation (12), while the failure rate function ( ) is defined according to Equation (15). The material's parameters are the following: Elastic modulus = 20MPa, shear modulus = 6.67MPa, = 50 or = 80, and = 300 K, while the bond energy has been adopted to be = 33.82 (= 1.4 ⋅ 10 J), = 2300 [9], and = 0.012, = −0.8. The material is subjected to two deformation cycles with 1 ≤ ( ) ≤ 2 characterized by strain rates equal to = 0.2s and = 2s . In Figure 7, the stress-deformation response is illustrated for the two different networks considered and for the two strain rates adopted; it can be noted that the polymer with longer chains ( = 80) responds more softly, and that in the high strain rate case the damage is lower because of the short time available to the chains to break during the deformation process (Figure 8b).

Mechanical Response of a Single Network Polymer in Presence of Damage
In the second case, we analyze the response of a single network polymer by accounting for the occurrence of damage without any internal rearrangement features (k d = 0 Hz). The CCDF time evolution, . ρ f , is here evaluated according to Equation (12), while the failure rate function ω f (r) is defined according to Equation (15). The material's parameters are the following: Elastic modulus E = 20 MPa, shear modulus µ = 6.67 MPa, N = 50 or N = 80, and T = 300 K, while the bond energy has been adopted to be w = 33.82 k B T (= 1.4 · 10 −19 J), E b = 2300 k B T [9], and A = 0.012, γ = −0.8. The material is subjected to two deformation cycles with 1 ≤ λ(t) ≤ 2 characterized by strain rates equal to . λ = 0.2 s −1 and . λ = 2 s −1 . In Figure 7, the stress-deformation response is illustrated for the two different networks considered and for the two strain rates adopted; it can be noted that the polymer with longer chains (N = 80) responds more softly, and that in the high strain rate case the damage is lower because of the short time available to the chains to break during the deformation process (Figure 8b). , while the Kuhn segment stiffness parameter has been adopted to be = 2300 .

Mechanical Response of a Polydisperse Polymer in Presence of Damage
In this case, we consider the response of a polydisperse polymer; the material's characteristics have been assumed the same as those presented in the previous section, while now the material is assumed to be made of two networks, one with = 50 and the second with = 80. The volume fractions of the two networks have been assumed to be ( ) = 0.8 , ( ) = 0.2 for the first polydisperse polymer and ( ) = 0.2, ( ) = 0.8 for the second multi-networks polymer. The material is subjected to two deformation cycles, with 1 ≤ ( ) ≤ 2 in the first cycle and 1 ≤ ( ) ≤ 3 in the second one; each of the deformation history has been assumed to be characterized by the strain rates equal = 0.2s or = 2s . has been assumed, while the Kuhn segment stiffness parameter has been adopted to be = 2300 .
In Figure 8, the response of the material is shown for the two double networks polymers and for the two strain rates considered; the case of the polymer richest of the shortest chains ( = 50) shows  (a) and = 2s (b). The bond energy strength has been assumed to be = 33.82 , while the Kuhn segment stiffness parameter has been adopted to be = 2300 .

Mechanical Response of a Polydisperse Polymer in Presence of Damage
In this case, we consider the response of a polydisperse polymer; the material's characteristics have been assumed the same as those presented in the previous section, while now the material is assumed to be made of two networks, one with = 50 and the second with  has been assumed, while the Kuhn segment stiffness parameter has been adopted to be = 2300 .
In Figure 8, the response of the material is shown for the two double networks polymers and for the two strain rates considered; the case of the polymer richest of the shortest chains ( = 50) shows . λ = 2 s −1 (b)) with increasing amplitude. The bond energy strength w = 33.82 k B T has been assumed, while the Kuhn segment stiffness parameter has been adopted to be E b = 2300 k B T.

Mechanical Response of a Polydisperse Polymer in Presence of Damage
In this case, we consider the response of a polydisperse polymer; the material's characteristics have been assumed the same as those presented in the previous section, while now the material is assumed to be made of two networks, one with N 1 = 50 and the second with N 2 = 80. The volume fractions of the two networks have been assumed to be q(N 1 ) = 0.8, q(N 2 ) = 0.2 for the first polydisperse polymer and q(N 1 ) = 0.2, q(N 2 ) = 0.8 for the second multi-networks polymer. The material is subjected to two deformation cycles, with 1 ≤ λ(t) ≤ 2 in the first cycle and 1 ≤ λ(t) ≤ 3 in the second one; each of the deformation history has been assumed to be characterized by the strain rates equal . λ = 0.2 s −1 or . λ = 2 s −1 .
In Figure 8, the response of the material is shown for the two double networks polymers and for the two strain rates considered; the case of the polymer richest of the shortest chains (N 1 = 50) shows a stiffer response with respect to the one with q(N 2 ) = 0.8. As in the previous parametric analyses, a high value of the strain rate implies a lower energy dissipation and the response is very close to the purely elastic case.

Simulations of Experimental Tests
In the present section, we finally consider the simulation of experimental tests [46]. In order to get a satisfactory agreement with the experimental results, we identify that the optimal number of Kuhn segments in the chains of each network is N 1 = 15 (q(N 1 ) = 0.8), N 2 = 80 (q(N 2 ) = 0.1), and N 3 = 160 (q(N 3 ) = 0.1). A deformation cycle with 1 ≤ λ ≤ 2.5 with . λ = 0.025 s −1 has been considered, while E = 5 MPa, µ = 1.67 MPa, and T = 300 K. Moreover, with the same monomer used in each network, the damage parameters have been adopted to be as follows for all the three networks: w = 33.82 k B T (= 1.4 · 10 −19 J), E b = 2300 k B T [9], and A = 0.0019, γ = −0.8. Moreover, the material is supposedly incompressible and to be stretched along the x-direction (parallel to the r x axis). The experimental response and the simulated one, for both the cases of elastic and damage behavior, are illustrated in Figure 9a; the usual stiffening behavior can be appreciated for the elastic case (to this end the damage rate function has been made to vanish, ω f (r) = 0), while a considerable dissipation of energy can be appreciated for the response in presence of damage. a stiffer response with respect to the one with ( ) = 0.8. As in the previous parametric analyses, a high value of the strain rate implies a lower energy dissipation and the response is very close to the purely elastic case.

Simulations of Experimental Tests
In the present section, we finally consider the simulation of experimental tests [46]. In order to get a satisfactory agreement with the experimental results, we identify that the optimal number of Kuhn segments in the chains of each network is = 15 ( ( ) = 0.8), = 80 ( ( ) = 0.1), and = 160 ( ( ) = 0.1 ). A deformation cycle with 1 ≤ ≤ 2.5 with = 0.025s has been considered, while = 5MPa, = 1.67MPa, and = 300 K. Moreover, with the same monomer used in each network, the damage parameters have been adopted to be as follows for all the three networks: = 33.82 (= 1.4 ⋅ 10 J) , = 2300 [9], and = 0.0019 , Moreover, the material is supposedly incompressible and to be stretched along the -direction (parallel to the axis). The experimental response and the simulated one, for both the cases of elastic and damage behavior, are illustrated in Figure 9a; the usual stiffening behavior can be appreciated for the elastic case (to this end the damage rate function has been made to vanish, ( ) = 0), while a considerable dissipation of energy can be appreciated for the response in presence of damage. Further, the cross-section of the distribution function ( ) related to the first network ( = 15, the one present with the greatest volume fraction ( ) = 0.8) along the (Figure 9b,d) and (Figure 9c,e) axis evaluated at various steps of the deformation history has been determined. The deformation states corresponding to the points indicated with A, B, C, and D in Figure 9a have been considered for the elastic as well as for the damage response; it can be appreciated that such a distribution widens in the direction of the applied deformation (Figure 9b,d) and narrows in the normal direction. It can be noticed that the occurrence of the damage is reflected in the decreasing of the area under the distribution function (Figure 9d) due to the chains loss. On the other hand, the chains oriented in the direction normal to the applied stretch do not fail easily, as witnessed by the distribution function's sections along the axis; the sections' pattern related to the elastic response ( Figure 9c) is very similar to that related to the behavior in presence of damage (Figure 9e). . λ = 2 s −1 (b)) with increasing amplitude (a). The bond energy strength w = 33.82 k B T has been assumed, while the Kuhn segment stiffness parameter has been adopted to be E b = 2300 k B T. Cross-section of the distribution function ϕ(r) of the first network along the r x (b,d) and r y (c,e) axis at various steps (see the squares in (a)) of the deformation history. Case of purely elastic response (b,c) and response in presence of damage (d,e). Further, the cross-section of the distribution function ϕ(r) related to the first network (N 1 = 15, the one present with the greatest volume fraction q(N 1 ) = 0.8) along the r x (Figure 9b,d) and r y (Figure 9c,e) axis evaluated at various steps of the deformation history has been determined. The deformation states corresponding to the points indicated with A, B, C, and D in Figure 9a have been considered for the elastic as well as for the damage response; it can be appreciated that such a distribution widens in the direction of the applied deformation (Figure 9b,d) and narrows in the normal direction. It can be noticed that the occurrence of the damage is reflected in the decreasing of the area under the distribution function (Figure 9d) due to the chains loss. On the other hand, the chains oriented in the direction normal to the applied stretch do not fail easily, as witnessed by the distribution function's sections along the r y axis; the sections' pattern related to the elastic response (Figure 9c) is very similar to that related to the behavior in presence of damage (Figure 9e).

Conclusions
The macroscopic mechanical response of polymers, usually characterized by a high nonlinearity degree, damage, and strain rate dependence, depends on many complex mechanisms taking place at the microscale. Among them, it is worth mentioning the existence of chains reversible bonds, bond exchange, chains sliding, chains' bond scission, and chains uncoiling. In the present study, a theoretical micromechanical model enabling to study the time-dependent response of polymers by accounting for the damage mechanism taking place at the microscale has been proposed. The proposed approach has been developed by introducing the so-called Chain Configuration Density Function (CCDF), defined over the chains configuration space; starting from the network's chain statistics, the evolution of such a function has been determined by accounting for the deformation, the rate-dependence and the damage mechanism. The developed approach has been up-scaled to the mesoscale by integrating the main field quantities over the above mentioned chains configuration space. The knowledge of the CCDF at a given instant of the deformation history provides all the necessary information required to fully know the mechanical state of the material.
The thermodynamic consistency of the theoretical framework has been illustrated, and its reliability in modelling the response of polymers has been illustrated through some parametric examples involving the rate dependent response of elastomers and the damage behavior under deformation cycles. An experimental test has also been simulated, and the corresponding CCDF evolution has been shown and discussed. The main advantage of the developed theoretical model can be acknowledged in its physics-based micromechanical nature; the main mechanisms involved at the scale of the polymer's network reflect the macroscopic material's mechanical characteristics, and the need for few parameters having a clear physical meaning (such as the shear modulus, the loss and storage modulus, and the cross-link bond strength) enable the approach to be easily used and implemented in computational code to predict the complex response of polymers. Strain energy of a single Kuhn segment w Bond energy between two segments of a polymer chain γ Parameter of the damage rate function ω f ϕ 0 (r), ϕ(r) Distribution function of the end-to-end vector in the stress-free and in the current state, respectively λ Unidimensional stretch measure ω f (r) Failure rate related to the chains with end-to-end vector r ρ(r) Chain Configuration Density Function (CCDF) . ρ, ρ f Total CCDF rate, CCDF rate due to deformation, viscous effects and failure, respectively ψ Deformation energy for a single chain Ψ 0 , Ψ Network's deformation energy per unit volume in the stress-free and in the current configuration, respectively σ Chauchy stress