Strain-induced crystallisation in natural rubber: a thermodynamically consistent model of the material behaviour using a serial connection of phases

A thermodynamically consistent concept to model the strain-induced crystallisation phenomenon using a multiphase approach is discussed in Loos et al. (CMAT 32(2):501–526,2020). In this follow-up contribution, the same mechanical framework is used to construct a second model that sets the same three phases in a serial connection, demonstrating an alternative to the former parallel connection of phases. The hybrid free energy is used to derive the constitutive equations. The evaluation of the Clausius–Duhem inequality ensures thermomechanical consistency. The model is based on a one-dimensional derivation that extends with the concept of representative directions to a three-dimensional anisotropic model. After the step-by-step derivation, the performance of the model is analysed in detail, including its comparison to the well-known Flory model, its evaluation for infinite fast and slow excitations, its simulation of uniaxial cycles and its validation via relaxation experiments. Finally, the model is compared comprehensively to the former parallel model showing their equivalent reason for existence.


Introduction and state of the art
High deformations of natural rubber (NR) result in a change of the molecular orientation of its network. Since its discovery in 1925 by Katz [33], the phenomenon of strain-induced crystallisation (SIC) in NR has become the focus of experimental and modelling investigations. Although SIC is a challenging subject due to its complex thermodynamics, polymer physics and sophisticated kinetics, SIC is a rewarding subject for academic research and industrial applications. It is motivated by the beneficial influence of SIC on the mechanical material characteristics such as NR's superior crack growth resistance and excellent tensile properties [2,54]. The main application made from NR is vehicle tires, where 60 -70 % of the world's NR production is utilised. The kinetics of SIC is a highly investigated topic shown with a large number of recent publications and conference contributions. Several studies focus on different aspects during the crystallisation process, such as the nucleation of crystals and their orientation, whereas others focus on different materials and loading conditions [3]. The crystallisation process during the loading of natural rubber is measured in situ with synchrotrons using wideangle X-ray scattering (WAXS) since the beginning of 2000 [60,61]. Huneau published an overview of different l (t) = l a (t) + l n (t) + l c (t) . (1) The existence of a non-crystallisable amorphous phase in polymers was first proved experimentally in 1980 [47], where Menczel and Wunderlich refer to a rigid amorphous phase [46,64]. Since the model uses a serial connection, the total stress is equal to the stresses of the single phases σ = σ a = σ n = σ c . The sketch in Fig.  1 visualises one of these parts in the virtual intermediate configuration, whereas the sketch in Fig. 2 shows the Fig. 1 Phases inside the elastomer in a serial connection in their virtual intermediate configuration: crystallisable amorphous phase with its initial length l a0 (left), non-crystallisable amorphous phase with its initial length l n0 (middle) and crystalline phase with its initial length l c0 (right)

Fig. 2
Phases inside the elastomer in a serial connection: crystallisable amorphous phase with its length l a (left), non-crystallisable amorphous phase with its length l n (middle) and crystalline phase with its length l c (right) elastomer part in the current configuration in stretching direction. In the virtual intermediate configuration, the crystallinity is equal to that of the current configuration, whereas the three phases are undeformed. A virtual experiment helps to envision this configuration: the sample is deformed to a certain state, where it has a defined crystallinity and certain stretches of the individual phases. Then, abruptly, the stress is removed such that the crystallinity remains constant. This virtual experiment follows the method when the viscous strain of a damper is identified in a viscoelastic Maxwell model. The three phases result to their initial lengths l a0 , l n0 , l c0 without any mechanical deformation. The modelling of the microstructure of the elastomer is conducted relative to the unloaded configuration. The crystallinity x is defined by the ratio of the length of the crystalline fraction to the total length in the virtual intermediate configuration Thus, the crystallinity index is related to the volume. The relation between the non-crystallisable amorphous fraction and the crystalline phase is assumed to be directly proportional i.e. each crystalline part possesses a non-crystallisable amorphous part [42]. The mobile amorphous fraction consequently is l a0 l 0 . The straightforward derivation starting from the intermediate configuration in addition to the application of the above definitions shows that the fraction of the crystallisable amorphous phase is expressed as Since the fraction of the crystallisable amorphous phase is non-negative, the crystallinity is limited by a maximum value x 0 . The advantage of this concept is that the maximum crystallinity x 0 is smaller than 1 in comparison with models where the degree of crystallinity can reach 100%. In contrast, the maximum crystallinity x 0 reaches around 20% in reality; see synchrotron measurements published in [10,55]. With these abbreviations, each length fraction simplifies to (stretch of the non-crystallisable amorphous fraction) The following calculation beginning with Eq. (1) derives the expression for the total stretch λ and the stretch of the amorphous phase λ a dependent on crystallinity x and the stretches of the remaining phases λ n , λ c .
Consequently, assuming an amorphous stretch of λ a > 0 and x ≤ x 0 (and setting x = x 0 ), the numerator of the fraction in Eq. (5) results in which is the condition that assures the stretch of the mobile amorphous phase λ a to be non-negative. Furthermore, the well-known model published by Flory in 1947 [20] is a special case included in Eq. (5). With x 0 = 1, the stretch of the amorphous phase λ a reads as where the limit of the stretch for the amorphous phase for x → 1 is lim The stretch of the amorphous phase increases continuously, i.e. it reaches unbounded values. Here, the physical interpretation of Flory's λ c and the λ c of the proposed model is different. Since the proposed model introduces λ c = l c (t) l c0 as an internal variable, which is defined by the corresponding evolution Eq. (45) in Sect. 3.4, the physical interpretation is inherently complicated.
Additionally, under the assumption of high stiffness of the non-crystallisable phase and therefore limited to λ n = 1, the limit for x → x 0 for the stretch of the amorphous phase given in Eq. (5) is unlimited positively as x approaches to x 0 In the following, the stretch of the amorphous phase of Flory's model [20] is compared to the model presented here. In Flory's model, increasing crystallinity x leads to an increase in the stretch of the mobile amorphous phase λ a for a fixed total stretch of λ = 6, as shown in Fig. 3. Moreover, an increasing stretch of the crystalline phase λ c leads to a decrease in the stretch of the mobile amorphous phase λ a visualised in Figs. 3 and 4. For example, with constant values for the maximum crystallinity x 0 = 0.5, the crystallinity x = 0.2, the stretch of the non-crystallisable phase λ n = 1 and the stretch of the crystalline phase λ c = {1, 3, 5, 6}, the here presented model shows the following behaviour: As the total stretch λ increases, the stretch of the mobile amorphous phase λ a increases as visualised in Fig. 5. Moreover, a higher fixed stretch of the crystalline phase also increases the stretch of the amorphous phase, though with less effect as the stretch of crystalline fraction λ c goes to 6. A further point to mention is that an increasing crystallinity x leads to an increase in the stretch of the mobile amorphous phase λ a depicted in Fig. 6 for a fixed stretch of the crystalline phase λ c = 2.

Modelling approach: the structure of the constitutive model
The here presented model is based on the framework published in [42]. With the thermomechanical basics shown in Lion & Johlitz [39], the directional approach applied in Lion et al. [40], based on the works of Miehe & Göktepe et al. [25,49,50], Guilie, Thien-Nga & Le Tallec et al. [27,58], the model is further constructed to focus on the accurate description of SIC with its hysteresis characteristics.

The multiplicative split of deformation gradient
The deformation gradient F is decomposed multiplicatively into the isochoric partF and the volumetric part F, consequently F :=F ·F, first introduced by Flory [21]. Thus, changes in shape and volume are considered independently. The volumetric part is defined asF := J 1 3 1 with the Jacobian of the deformation gradient J = det(F).
The isochoric Green-Lagrange strain tensor is defined asÊ := 1 2 (Ĉ − 1), whereĈ =F T ·F = J − 2 3 C is the isochoric right Cauchy-Green deformation tensor. The Cauchy stress tensor is split into a volumetric and a deviatoric part where p := − 1 3 tr(T) describes the hydrostatic pressure and T D is the deviatoric stress. The second Piola-Kirchhoff stress is defined asT := J F −1 ·T·F −T , and in a similar manner, the deviatoric second Piola-Kirchhoff The latter is only denoted as deviatoric part. A correction term is introduced in Eq.

The concept of representative directions
The concept of representative directions includes an appropriate micro-to-macro transition of micromechanically motivated models. Averaging operations on the unit sphere are carried out with integrals presented in the following. Continuous orientation is replaced by a discrete set of orientations. For the current simulation, the following steps are taken. First, arbitrary states of deformation are projected to one-dimensional deformations along with specific directions. Each directional stretch is used within the one-dimensional material model to calculate directional stresses, crystallinities and other quantities. The directional stresses are then combined via an averaging scheme along a unit sphere to a three-dimensional stress tensor. The major advantage of this approach is that the formulation of constitutive equations becomes accessibly simple since the formulation is executed in one-dimensional states of deformation and stress. Moreover, the micro-sphere concept includes direction-dependent effects like anisotropic behaviour naturally. The one-dimensional quantities are specified per direction e α defined in the spherical coordinate system shown in Fig. 8. The Cartesian unit vectors e 1 , e 2 , e 3 form an orthogonal basis. e α = sin ϑ cos ϕ e 1 + sin ϑ sin ϕ e 2 + cos ϑ e 3 The averaging operator A, which calculates three-dimensional quantities through integration over the sphere of one-dimensional quantities, is defined as where f α is a physical quantity related to the direction e α . Its first time derivative isḟ = A[ḟ α ]. Regarding the operator A, the physical quantities related to the direction e α are introduced: (isochoric part of the specific hybrid free energy) (isochoric part of the specific entropy) The isochoric part of the stress power is formulated as The directional stress tensorT α =σ where is an arbitrary scalar. The thermodynamic consistency is ensured: the additional Ĉ −1 term does not produce stress power in the Clausius-Duhem inequality, Sect. 3.4, due to the orthogonality relation d dt det Ĉ −1 =Ĉ −1 :Ċ = 0 of the unimodular Cauchy-Green deformation tensorĈ. Inserting Eq. (17) intỗ T :Ĉ = 0 and usingĈ −1 :Ĉ = tr(1) = 3, the scalar is given by Consequently, the isochoric second Piola-Kirchhoff stressT has the final form where the relationλ 2 α = e α ⊗ e α :Ĉ is used.

The concept of the hybrid free energy
The constitutive model is based on the formulation of the hybrid free energy [41], which is decomposed additively into two contributions The Gibbs-type energy contribution g represents the volumetric part of the material behaviour and depends on pressure p, temperature θ and possibly on the total crystallinity x.
The energy contribution of Helmholtz-typeψ describes the isochoric part of the material behaviourψ = A ψ α λ α , θ, x α such that the isochoric part of the material behaviour is modelled by using the concept of representative directions explained in Sect. 3.2. Within the current model with its total directional stretch given in Eq. (5), the isochoric part of the specific hybrid free energy depends on three single stretches, i.e. ψ = A ψ α λ a ,λ n ,λ c , θ, x α due to the three individual material's fractions.

Total directional Helmholtz free energy
The specific directional Helmholtz free energies of the single phases are stated aŝ which are all considered to be dependent on their corresponding stretch. The directional total Helmholtz free energy per unit massψ α is the sum of the free energies of the individual phaseŝ where the last term originates from the entropy of mixing, where R is the universal gas constant, M is the molar mass and the function γ (x α ) is the mixing ratio. Its derivation is found in detail in [42]. The empirical exponents m, n, k are introduced to formulate the model with more flexibility to represent experimental data. They do not influence on the thermodynamical consistency of the constitutive model. The mixing ratio of the elastomer part, derived in [42], is For the evolution equation of crystallinity, the first derivative of the mixing ratio with respect to the crystallinity is required and reads as The individual contributions of the isochoric part of the free energyψ α defined in Eq. (21) consist of the directional free energy densities of the crystallisable and non-crystallisable amorphous elastomer, which are both considered to be entropy elastiĉ and the directional free energy density of the crystalline elastomer, which is considered to be energy elastiĉ Here,ψ 0a ,ψ 0n andψ 0c are the initial free energy densities, s 0a , s 0n and s 0c are the initial specific entropies and ϕ αa (λ a ),φ αn (λ n ) andφ αc (λ c ) are the directional strain energies of the amorphous and crystalline phases. The stretch-dependent mechanical contributions to the directional free energies, i.e. the individual strain energieŝ ϕ αa ,φ αn ,φ αc can be chosen arbitrarily, e.g. empirical or physical-based models. The temperature θ G is a reference temperature of the material (e.g. glass transition temperature), and θ B is the basic reference temperature for entropy elasticity, cf. Lion [38, p.29]. θ θ B is the purest form of the temperature dependence of the strain energy functionφ in the context of entropy elasticity [48]. The partial derivative of the isochoric part of the free energy with respect to the crystallinity is Since the Helmholtz-type free energy is dependent on five quantitiesψ = A ψ α λ a ,λ n ,λ c , θ, x α , its time derivative is

Strain energies of elastic materials
Each material fraction of the model, i.e. the crystalline in addition to the crystallisable amorphous and the non-crystallisable amorphous phase, is described with a material model of elasticity. Each phase can be treated separately when it comes to the aspect of finding a suitable material model. The corresponding stress tensors are derived from the strain energy density function ϕ [45, p.837], which is defined per unit volume. The following material models for each phase were chosen after a parameter identification using the multistep experimental data for stretch up to λ ≤ 3.5 by [42, Fig. 11] representing the amorphous material response. Since the presented model uses a serial connection for all phases, the softest material determines the stiffness and the material's stress response. For the amorphous phase (crystallisable and non-crystallisable), the model of extended tube is used, which was introduced by Kaliske & Heinrich [28,32] and is based on physical considerations on the molecular scale. The total elastic free energy change, or elastic potential, of the extended tube model [32, p.606] is divided into two additive terms: the cross-link part and the elastic partφ =φ c +φ e , respectively. The cross-link strain energy for the extended tube model iŝ and for the constraint contribution, the related strain energy readŝ The constants G e , G c , β and δ are material parameters: shear moduli, completeness of cross-link reaction and inextensibility parameter. IĈ =λ 2 1 +λ 2 2 +λ 2 3 =λ 2 + 2λ −2 is the first invariant of the isochoric right Cauchy-Green tensorĈ for uniaxial extension. Here, the strain energies of the amorphous phaseφ αa (λ a ),φ αn (λ n ) in Eqs. (24)(25) are modelled with the extended tube model using the same parameters shown in Tab. 1.
For the crystalline phase, a Yeoh-type model [65] with three terms n = 3 is chosen, such that the model has three material parameters C 1 , C 2 , C 3 . The phenomenological model is using the strain energy density functionφ Here, the strain energy of the crystalline phaseφ αc (λ c ) in Eqs. (26) is modelled with this Yeoh-type model using the parameters shown in Tab. 1 representing a stiff material response.

Volumetric part of the free energy
The model for the volumetric part of the free energy of Gibbs type [42] is It is assumed that the degree of crystallinity does not influence the volume. Thus, the Gibbs-type free energy contribution depends on two quantities g( p, θ). Its time derivative is Thus, the time derivative of the hybrid free energyφ is

Thermodynamical considerations: the Clausius-Duhem inequality
The Clausius-Duhem inequality has to be satisfied by the constitutive law such that the second law of thermodynamics is satisfied. With the stress power given in Eq. (11), the Clausius-Duhem inequality on the reference configuration reads as With the time derivative of the free energy The stress powerT :Ė = A[σ αλα ] given in Eq. (16) reads as after insertion of the time derivative of the total stretch calculated from Eq. (5). The Clausius-Duhem inequality considering the directional quantities introduced in Sect. 3.2, inter alia, and the time derivative of the hybrid free energy Eq. (34) reads as Here, the variablesλ n ,λ c , x α are used as internal state variables, which are defined by the solution of ordinary differential equations (ODEs). To evaluate this inequality, the stretchλ a of the mobile amorphous phase given in Eq. (5) is used as the process variable. The definition of state and process variables is flexible and has to be set once. Requiring that Eq. (38) is non-negative for arbitrary values of the individual variables [12], the following quantities result directly in consequence: the volume strain the specific entropy s is the negative derivative of the free energy with respect to the temperature and the directional heat flux where κ is the thermal conductivity. From the inversion of the linear pressure-dependent volume strain ε vol ∼ p, the pressure p = p (ε vol , θ, x) can be easily computed. The total stress is computed with the usage of the The evolution equation of the directional crystallinity including a prefactored positive function β(θ, ...) ≥ 0 reads asẋ using the partial derivative given in Eq. (27). In the current work, the function β(θ, . . .) in the evolution equation describes the temperature dependence of the crystallisation rate among other optional dependences on internal state variables. It is possible to add dependences on stress, stretch or internal state variables to this function, which remains as a degree of freedom. The temperature dependence is represented by the approach by Williams, Landel, Ferry [63].
with c 1 = 17.44 and c 2 = 51.6 K. Regarding temperature, the function β(θ, . . .) remains constant for the isothermal case. In the current work, the experiment and the simulation are both isothermal at θ = 300 K. The evolution equations for the stretches of the non-crystallisable amorphous phase and the crystalline phase areλ where η c , η n ≥ 0 are the scalar proportional factors of the single phases comparable to the standard relation of linear viscoelasticity for dashpotsλ = 1 η σ .

Model evaluation for special cases: fast and slow excitations
In this section, the index α, which refers to the directional quantity is neglected.

Evaluation for fast deformations
The model can be evaluated for fast and slow thermomechanical excitations. First, the evaluation for infinitely fast deformationsλ → ∞ from the initial stateλ = 1, x = 0 is carried out. For such a high stretch rate, all internal state variables are frozen, i.e. the evolution of crystallinity is not inducedẋ = 0 such that no crystallinity evolves x = 0. Thus, crystallisation remains zero as shown in Fig. 10 by the solid grey line. Likewise, the evolution of the stretches of the non-crystallisable amorphous phase and crystalline phase remains zeroλ n = 0,λ c = 0 such that the stretch of the crystallisable amorphous phase equals the total stretchλ =λ a .
The stress is determined from Eq. (42) by the behaviour of the amorphous phaseσ (λ a , . The sketch of the stress, where no hysteresis is observed, is plotted qualitatively with the solid grey line in Fig. 9. The stress response is higher compared to the equilibrium case.

Evaluation for slow deformations
Secondly, the case of infinitely small stretch rate is analysedλ → 0. Starting from the initial stateλ = 1, x = 0, the equilibrium stress and crystallinity response are reached, plotted in solid black line in Figs. 9 and 10. The following calculations show the absence of hysteresis within the equilibrium response. For the steady state, the time derivativesẋ,λ n ,λ c remain zero. The evolution equation for the degree of crystallinity allows the computation of the equilibrium degree of crystallinity as a function of the individual stretches of each phase. Withẋ = 0, Eq. (43) follows as such that the equilibrium crystallinity x eq can only be calculated with the help of numerical methods due to its implicit expression. The evaluation shows that the stretch of the crystalline phase and the stretch of the non-crystallisable amorphous phase are dependent on the process variablesλ a and x Therefore, the equilibrium degree of crystallinity can be expressed as a function depending on the stretch of the amorphous phase x = x eq λ a , f (λ a , x eq ), g(λ a , x eq ) ⇒x eq (λ a ,x eq ) ⇒x eq (λ a ) .
Thus, the stretch of the amorphous phase is dependent on the total stretcĥ can also be written as a function of the total stretchσ =σ (λ) as well as the crystallinity being dependent on the total stretch x = x(λ). These expressions cannot be resolved analytically due to the implicit formulation. To summarise, in the case of infinitely slow deformations, all quantities of the model can be expressed as a function of the total stretch. Finally, the dependence of stress and the crystallinity response on the global stretch results in the absence of a hysteresis. Figures 11 and 12 show the simulation results for infinitely fast and slow deformations for the here presented model (serial model) and the former presented model (parallel model) [42].

Fig. 11
Simulation of infinitely small and fast excitation: stress over total stretch, parallel and serial model, 3D simulation, experimental data of [10] used for guidance Fig. 12 Simulation of infinitely small and fast excitation: crystallinity over total stretch, parallel and serial model, 3D simulation, experimental data of [10] used for guidance

Macro-micro-macro transition
In the context of a simulation of a homogeneous uniaxial cyclic tension test, the first step within the concept of representative directions is to compute the isochoric directional stretchesλ α from the relationλ α = e α ⊗ e α :Ĉ as presented in Sect. 3.2. The 11-component is consequentlyĈ 11 =λ 2 e 1 ⊗e 1 + 1 λ e 2 ⊗e 2 + 1 λ e 3 ⊗e 3 for uniaxial tension. Next, for each direction, the directional stressσ α is computed by Eq. (42) via the solution of the coupled system of ODEs for x α ,λ n andλ c . Having the directional stressesσ α computed, the second Piola-Kirchhoff stress tensorT is calculated with the relation from Eqs. (10) and (19): To determine the hydrostatic pressure p, the relationT 11 = 3 2T 11 from Lion et al. [40,Eq.(67)] is taken into account. The second Piola-Kirchhoff stress tensorT is related to the first Piola-Kirchhoff stress tensor P by the following operation Moreover, for the micro-to-macro transition, i.e. to computeT from σ α in Eq. (49), the following discrete where N is the number of integration points alias directions and w α is the weight defined per direction for which N α=1 w α = 1 is valid. The widely used integration points and weights of Bazant and Oh [4] are here varied from 21 to 37 directions per half sphere, cf. Figs. 13, 14. Many sets of orientation vectors e α and weights w α are given in this work and also in [50]. For instance, the reader is referred to [31] and [62] for detailed discussions and the comparison of different integration points   16 Crystallinity over stretch of unfilled natural rubber, integration points based on [19] and methods. The usage of a high amount of points does not improve the numerical accuracy significantly. Vice versa, the usage of many integration points leads to a large number of ODEs, which has to be solved at each Gauss point in the context of the finite element method (FEM). As a consequence, in this study, instead of using many integration points, the parameters are optimised for the minimum number of integration points possible. In addition to the mentioned reasons, the change of the simulated stress response is neglectable and only noticeable when the unloading path reaches the loading path at a stretch around λ m ≈ 3. In contrast, the simulated crystallinity response decreases slightly in the unloading path. Furthermore, Miehe et al. stated that the 21-point integration scheme provides sufficient accuracy [50, p.2639, l.14]. Nevertheless, in similar to the work of [31], different amounts of integration points presented by Fliege and Maier [19] are successfully used and presented in Figs. 15, 16 for comparison reasons. In parallel with the expectations, changes in the stress or crystallinity response are not significant when increasing the number of directions, which underlines the statement of Miehe et al. cited above.

Optimisation strategy
In the following subsection, dealing with the parameter identification, the subscript α, which stands for 'directional', is omitted for notation clarity. The first-order ODEs developed in Sect. 3 form a fully coupled system. The evolution equations for the directional crystallinity x, stretch of the crystalline phaseλ c and stretch of the non-crystalline phaseλ n are presented in Eqs. (43), (46), (48), respectively. The stressσ is given in Eq. (42) and the constraint is used to compute the crystallisable amorphous stretchλ a from Eq. (5). This set of equations constitutes a coupled system for each direction. Considering an isothermal state and for the given time interval t ∈ [0, t end ], the coupled system can be written in the clear short forṁ where all evolution equations are highly nonlinear. This system is solved in each direction. Then, the total quantities are acquired by the weighted sum of each direction. The authors made use of a four-step 'hybrid optimisation' approach described in the following: 1. Generate a high number of Monte Carlo samples (10 Mio), uniformly distributed between lower and upper boundaries 2. Sort and pick parameter sets concerning stress and crystallinity 3. Use genetic algorithm with Monte Carlo's result as initial population matrix 4. Select a parameter set from the Pareto front as a compromise between stress and crystallinity response For the Monte Carlo sampling, i = 10 million parameter sets p i = [p 1i , p 2i , . . . , p ni ] T normally distributed in between lower and upper boundaries p l and p u are generated: where R ∈ [0, 1] is a random real number between zero and one. Then, favourites are selected to generate an initial population matrix used in the genetic algorithm gamultiobj in MATLAB's optimisation toolbox.
The aim is to identify the global optimum with these heuristic approaches. They are often used for coupled nonlinear differential equations. In particular, the reader is referred to [24,26,34] for further information about optimisation of ODEs in general and multiobjective optimisation. The total constraint minimisation problem in a permitted set between upper and lower boundaries is stated as where the objective function f(p) = [ f σ , f x ] includes two separated objectives. Here, parameter sets p = [ p 1 , p 2 , . . . , p n ] T laying on the Pareto front are identified with the following weighted objective functions whereσ i andx i are the stress and the crystallinity index obtained from the experiment, and σ i (p) and x i (p) are the stress and crystallinity index computed by the model, respectively. Furthermore, the stress response is divided into s different parts, which are weighted with c j to emphasise the importance of certain parts of the stress response, such as the characteristic plateau exemplary. The advantage of the Pareto front approach is that many parameter sets are available for selection from the Pareto front. This multiobjective optimisation allows treating the objectives separately. On the other hand, the parameter set presented in this work is a compromise between these two objectives: stress and crystallinity. There is no indication that the presented parameter set is the global optimum within the set boundaries. Certainly, the boundaries were selected such Fig. 17 Simulation of a uniaxial tensile test: total stress over total stretch using 21 directions that no parameter never lies on any boundary. The usage of the multiobjective optimisation is also motivated by the available experimental data for the crystallinity index. It is measured as the ratio of intensities of scattered electromagnetic waves in the reciprocal lattice. For further information on WAXS the reader is referred to [9, p.64f] and [10, p.66]. On the contrary, the simulated crystallinity is defined by a ratio of lengths, which is consequently related to the volume, cf. Eq. (2). Thus, the equality of the modelled and the experimental remains unsure.

Evaluation and comparison to experimental data
The simulation results are shown graphically in Figs. 17,18,19,20,21,22,23,24,25,26. The identified material parameters are shown in Table 1, where the set in the column '3D serial model' is valid for this subsection. The abbreviations 1D and 3D refer to the usage of the concept of representative directions. In the former published model (parallel model), directional values were presented (1D). Therefore, four parameter sets are presented in total to make a complete comparison. The experiment and the simulation are both conducted isothermally at θ = 300 K, which results in a constant value for β (θ, . . .) ≥ 0. The simulation results are compared to experimental data of Candau [10].
To begin with, the total stress response over total stretch is shown in Fig. 17, which show high concordance with the experimental data. Also, the crystallinity response in Fig. 18 satisfies the experimental data most satisfactorily. The details will be explained in the latter comparison, although it can already be said at this point that the results obtained are gratifying.
First of all, the analysis is dedicated to the individual directions. The single stretches, stresses and crystallinities of each direction are shown in Figs. 19, 20, 21, 22, 23, 24, 25, 26, where the non-weighted and weighted values are presented. The direct comparison of non-weighted and weighted stretches in Figs. 19 and 20 visualises the influence of the weights within the concept of representative directions. Since the weights sum up to 1, N α=1 w α = 1, it follows that only the non-weighted values represent real values in micro-mechanical manner. The averaging scheme makes the weighted and therefore reduced values only interpretable in a qualitative manner.
Since the number of one direction is not dedicated to the location in the unit sphere, Fig. 21 is key to the three dimensionality and visualises the concept of representative directions. The directional values are shown in the nodes and interpolated in between with a colour bar. The direction parallel to the stretching direction is the one with maximum stretch, which equals the total stretch, whereas the directions orthogonal to the stretching direction experience compression due to transverse deformation. The surrounded directions decrease radially to the stretching direction. The visualisation in Fig. 22 shows the directional stretches using 900 integration points by [19].

Fig. 22
Simulation of a uniaxial tensile test: directional stresses in the unit sphere, 900 integration points by [19] For stress and crystallinity, the same analysis is shown in Figs. 23, 24, 25, 26, 27, 28. Since the third direction is the most stretched one, it shows the highest stress (Fig. 23) and the highest crystallinity evolution (Fig. 25). Other directions show compression in stress, whereas the evolution of crystallinity is not influenced by compression, which follows experimental results. Another point is that no stress hysteresis occurs in directions, which are stretched below the critical stretch for SIC. In these directions, the crystallinity does not evolve either. Figures 28 and 27 show the spheric plots. In terms of the stress presented in Fig. 27, the high nonlinear upturn of the stress is the reason for the point of maximum stress on the tip. The reader is referred to the legend that reaches from 0 to 2.1 MPa. In terms of crystallinity in Fig. 28, it is confirmingly observed that crystallinity only occurs at high deformations.
To analyse mechanic consistency, the free energiesψ a ,ψ n ,ψ a , including the strain energieŝ ϕ a (λ a ),φ n (λ n ),φ c (λ c ) given in Eqs. (24)(25)(26) are presented in Fig. 30. They satisfy the polyconvexity. The polyconvexity of the free energy is one of the physical boundary conditions. The free energy approaches to the infinity asλ decreases to zero: lim

Model validation
The constitutive model is validated via two different deformation processes, where the material's time dependence is analysed. The material's relaxation behaviour in a one-step relaxation and two-step relaxation with decreasing stretch is investigated. For both experiments, videoextensometry is used to identify the exact stretch, which differs from the machine input due to sample clamping. For higher accuracy, the stretch measured via videoextensometry was 2.2% less at the maximum stretch. In the relaxation experiments, the stretch rate during loading and unloading isλ = 0.026 1 s , which is higher in comparison with the stretch rate ofλ = 0.0042 1 s for the cycles shown in the simulations of the former Sect. 5.3. Since the model is rate dependent, deviations in stress and crystallinity are to be expected.

Stress relaxation at a constant stretch, one-step relaxation
The given experiments are published in [42]. First of all, the simulations shown in Figs. 31 and 34 cover the stress relaxation qualitatively. In the following, a few characteristics are accentuated. First, the high concordance of the simulation and the experiment for a maximum stretch of λ max = 3 should be highlighted, where small stress relaxation is observed, i.e. 2.5% stress relaxation from maximum stress. Second, at the stretch of λ max = 6, the simulation shows faster stress relaxation than the experiment. In the simulation, the equilibrium stress is reached earlier. Moreover, the absolute and relative value of the amount of relaxed stress is lower in the simulation. The simulation covers stress relaxation due to SIC.

Stress relaxation at a constant stretch, two-step relaxation
Another experiment is conducted to clearly distinguish between stress relaxation when the sample is stretched below or above the critical stretch of crystallisation onset. A two-step relaxation is conducted starting at a constant stretch of λ I = 6 with a relaxation time of 6 hours, unloading directly to a second step at λ II = 3 for another 6 hours. The simulation of the stress response is shown in Figs. 33 and 34. Starting with Fig.  33, the stress highly agrees with the experiment. Especially at the second step, it excellently agrees with the Simulation of uniaxial tensile test: stress over total stretch for one-dimensional implementation experiment. In contrast, the stress relaxation in the first step shows similar characteristics as discussed in the former subsection: The equilibrium stress is reached earlier in time, and the relative and absolute amount of relaxed stress is smaller in the simulation compared to the experiment. In particular, with the help of the close up of the second step shown in Fig. 34, it is observed that stress reaches the equilibrium directly after unloading to the constant stretch of the second step. This increase in the stress is in parallel with the experimental data, showing that the proposed model captures two-step relaxation experiments. In addition, the crystallinity simulation of the two-step relaxation presented in Fig. 32 is in agreement with the kinetics of SIC. Crystallinity evolves during the first step, where stress relaxation takes place. In contrast, crystallinity vanishes during the second step, where the stretch lies below the critical stretch of crystallisation onset, and stress relaxation is absent. The reader is referred to [10,55] for similar experimental data.

Evaluation and comparison to the parallel model
In [42], a concept to model strain-induced crystallisation was presented using a parallel connection of three phases. Since the simulation results in [42] only covered the 1D approach, both models (parallel and in-series connection) are first compared with each other in 1D. This helps to understand the micro-to-macro transition. Further in this subsection, both models are compared with their 3D results. The micro-sphere approach of the parallel model is successfully implemented, as proposed in the outlook of Loos et al. [42]. All identified parameters are presented in Table 1 for an effortless comparison. The 1D results are shown in Figs. 35 and 36, whereas the 3D results are presented in Figs. 37 and 38. Concerning the simulation results, one detail has to be mentioned: the comparison follows from separately identified parameter sets, which produce unique results as the parallel model and the serial model are different from each other and the parameters result from the optimisation approach explained in Sect. 5.2. The experimental data from Candau [10] are used as guidance for the reader.
The stress stretch and the crystallinity stretch graphs have several characteristics: Capture of the amorphous stress response In the undeformed state, all the material is composed of the amorphous phase. Therefore, the material's stress response in the beginning of a uniaxial tension experiment is totally covered by the amorphous phase. The parameters of the amorphous phase are identified by stretching the material less than the critical stretch of crystallisation onset λ critical ≈ 4.3. This almost linear loading path is captured perfectly for the presented unfilled NR. The main parameters, which influence the behaviour for low stretches, are the shear moduli of the crystallisable amorphous phase G e , G c , excluding the δ since it mostly influences the nonlinearity of stress at high stretches. The reunion of the unloading path with the loading path takes place at λ = 3, which underlines the absence of classical viscoelastic effects in this unfilled natural rubber. Another evidence is that the specimen does not lengthen during the slowly conducted tensile test, i.e.

Fig. 36
Simulation of uniaxial tensile test: crystallinity over total stretch for one-dimensional implementation Fig. 37 Simulation of uniaxial tensile test: stress over total stretch for three-dimensional implementation Fig. 38 Simulation of uniaxial tensile test: crystallinity over total stretch for three-dimensional implementation

Fig. 39
Simulation of uniaxial tensile test: directional stresses in the unit sphere, parallel model no recovery time is needed to return to its initial length. All presented data capture the amorphous response, where the one-dimensional results fit precisely, and the three-dimensional results vary slightly, which is caused by the integration over all directions.
Capture of the stress plateau As the amorphous phase relaxes short after the onset of crystallisation, the onset of crystallisation results in a stress relaxation, which is depicted in the plateau of the stress response close to λ critical ≈ 4.3 in the experimental results in Figs. 35 and 37. The authors see a correlation between the capture of the plateau and the shape of the hysteresis. Therefore, the benefit of the presented models is the ability to capture this stress decrease at crystallisation onset.
Capture of the nonlinear upturn of the stress response The nonlinearity of the stress response is influenced by the nonlinearity of the material models. For the in-series connection, a Yeoh-type model is used for the crystalline phase with its parameters C 1 , C 2 , C 3 . As the amorphous phase modelled by the extended tube model, the parameter δ mostly influences the nonlinearity upturn, which refers to an inextensibility parameter [32]. The higher this parameter, the higher is the nonlinearity upturn.
Capture of the area of the stress hysteresis The model is able to capture the area of hysteresis. The difference between loading and unloading is due to the crystallisation evolution. The area of the hysteresis gives evidence to the work per unit volume done during the cyclic test. In all simulations, the area of the hysteresis of the stress response is captured with high accordance.
Capture of the crystallinity shape First of all, when it comes to the comparison of the simulations with experiments in Figs. 35 and 38, the discrepancy between the crystallinities is qualitatively higher than the small discrepancy between the stresses. This underlines a behaviour observed during the derivation of the model: the concordance between stress and crystallinity response is of minor size, which motivates the introduction of the empirical parameters m, n, k similar to [42], cf. Section 3.3.1 first seen in the Helmholtz free energy in Eq. (21). In the current study, the main focus lies on the capture of the stress response, including its hysteresis. As mentioned in Sect. 5.2, it cannot be ensured that the experimental and the simulated crystallinity indices are identical. Thus, the good concordance in shape and hysteresis of the crystallinity is pleasant. The characteristics of the crystallinity response can be structured with the help of three aspects: first, the stretch at crystallisation onset is higher for the 3D serial model in comparison with the 3D parallel model in Fig. 38. Second, the beginning of the unloading path shows varying rates of increase in the crystallinity, depending on the material parameters. In dependence of the material composition and the kinetics of crystallinity evolution, the stress decreases, but the crystallinity can decrease [10, p.193] or increase [10, p.236]. Third, the qualitative representation of the shape, including the area of the crystallinity hysteresis is of interest. The simulations demonstrate that the 3D simulations provide better results for the crystallinity in comparison with the 1D simulations. Moreover, the 3D parallel model and the 3D serial model vary slightly in terms of crystallinity, but they have a similar behaviour depending on the material parameter. to concentrate on small differences between experiment and simulation. Both models have similar capabilities to simulate the experimental data.

Conclusion and outlook
A thermomechanical approach to model strain-induced crystallisation is presented in continuation of [42]. The here presented formulation of a constitutive model is based on serial connection of three different phases: the crystallisable amorphous phase, the non-crystallisable amorphous phase and the crystalline phase. The stretch relation proposed by Flory [20], which is based on a two-phase approach, is contained as a special case (x 0 = 1). Flory's theory and its in-series connection are taken up by Albouy & Sotta [1,57] exemplary.
The derivation of the serial approach is stated in detail from basic considerations such as multiplicative split of the deformation gradient and the concept of the hybrid free energy. Thermodynamic consistency is used to derive the constitutive equations. The model is evaluated via special cases such as infinitely fast and slow excitations. These evaluations prove physical and mechanical accuracy. Suitable parameter sets are identified using a multistep optimisation method. Within the simulation, ODEs of the first order were solved. The simulations show high agreement with experimental data. To validate the model, two approaches of stress relaxation are simulated, which both show high qualitative concordance with the experimental data. Finally, a complete and comprehensive comparison is made. First, the current model using serial connection is compared for one direction to the former model using parallel connection of phases. Second, the threedimensional simulations of both models are compared since they are equivalent to the experimental data. All in all, the parallel model and the serial model provide similar performance, although the parallel model has one evolution equation compared to three evolution equations of the serial model. All shown simulations cover the experimental data in high concordance in shape and absolute values. The models depict the viscoelastic behaviour of the material depending on the degree of deformation. They also represent cyclic processes precisely and simulate the development of SIC, including its hysteresis behaviour. The models are characterised by their simplicity. They do not make use of artificially introduced case distinctions, which could, for example, distinguish between loading and unloading.
Since all conducted simulations show isothermal results, the implementation of the equation of heat conduction is left for future works, which is of high interest since the phenomenon of strain-induced crystallisation is highly exothermal. Future studies could include filled NR, whereas the current study focuses on unfilled NR.
The evolution equations with internal variables developed here are similar to the theory of viscoelasticity with internal variables, but are not connected.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by Projekt DEAL.