Thermomechanical Modeling of Microstructure Evolution Caused by Strain-Induced Crystallization

Abstract

The present contribution deals with the thermomechanical modeling of the strain-induced crystallization in unfilled polymers. This phenomenon significantly influences mechanical and thermal properties of polymers and has to be taken into consideration when planning manufacturing processes as well as applications of the final product. In order to simultaneously capture both kinds of effects, the model proposed starts by introducing a triple decomposition of the deformation gradient and furthermore uses thermodynamic framework for material modeling based on the Coleman–Noll procedure and minimum principle of the dissipation potential, which requires suitable assumptions for the Helmholtz free energy and the dissipation potential. The chosen setup yields evolution equations which are able to simulate the formation and the degradation of crystalline regions accompanied by the temperature change during a cyclic tensile test. The boundary value problem corresponding to the described process includes the balance of linear momentum and balance of energy and serves as a basis for the numerical implementation within an FEM code. The paper closes with the numerical examples showing the microstructure evolution and temperature distribution for different material samples.

1. Introduction

The strain-induced crystallization (SIC) is a phenomenon of formation and degradation of tiny crystalline regions within the amorphous polymer structure due to the strain variation. It has a significant influence on the production and application of the affected materials, which strongly motivates its experimental investigation and numerical simulation. The experimental characterization of the process relies on techniques such as volume change measurements [1], electron microscopy [2], small-angle X-ray scattering [3], and terahertz spectroscopy [4]. However, the most frequently used technique is the in situ wide-angle X-ray diffraction. Amongst others, the method has been applied to study crystalline content and orientation [5], the SIC kinetics [6], the effects of dynamic loading [7], and the distribution of the crystallite size [8]. Infrared thermography is also worth mentioning, particularly in the context of investigating thermal effects [9,10,11].

The representative results of a cyclic tensile test performed for the unfilled natural rubber at constant ambient temperature and strain rate are shown in Figure 1. Here, the stress diagram (Figure 1a) forms a hysteresis and thus shows the dissipative nature of the SIC phenomenon. The volume fraction of the crystalline regions, the so-called degree of crystallinity, is also used to show the crystallization process. According to Figure 1a, crystalline regions begin to form after exceeding a stretch of $\lambda =4.3$. The degree of crystallinity rises with increasing strain until saturation at about $\lambda =6.0$ is reached. A further increase in total deformation may lead to inelastic deformations, which is not the subject of this study. The decrease of the degree of crystallinity during the unloading phase is less intensive than its growth by loading. The material becomes completely amorphous at a stretch of $\lambda =3.0$. It is also remarkable that the heat generated by the formation of crystalline regions leads to an increase in temperature, whereas the regression of the crystalline regions during the unloading phase causes a decrease in temperature. Accordingly, the temperature diagram (Figure 1b) has a form analogous to the one shown in Figure 1a.

With regard to the modeling and simulation, the first effort was made by Alfrey and Mark (1942) [13] who investigated the behavior of a single polymer chain. On the basis of the this work, Flory (1947, 1949) [14,15] formulated the classical thermodynamic theory of the SIC. In these works, the degree of crystallization is expressed by the change in crystallization temperature which, in turn, depends on strain. At a later stage, a phenomenological expression to describe the growth of the crystalline phase was proposed by Doufas et al. (1999) [16] using a modified Avrami equation. The work by Ahzi et al. (2003) [17] extends the hyperelastic model of Boyce et al. (1993) [18] by the phenomenological expression of Doufas et al. (1999) in order to explain the growth of the crystalline phase. The works by Negahban (2000) [19] and Rao and Rajagopal (2001) [20] used the large deformation theory to model the amorphous and crystalline phases separately. Afterwards, Tosaka et al. (2004) [21] developed a micromechanical model in which shorter chains are completely stretched under tension and form nucleation sites for crystals. The model by Kroon (2010) [22] deals with the anisotropic nucleation in unfilled rubbers. It predicts both stress–stretch hysteresis and the development of the degree of crystallization, under the assumption that the dissipative process is not caused by crystallization only but also by the viscoelastic behavior of the amorphous phase. An advanced micromechanical continuum model for partially crystallized polymers was developed by Mistry and Govindjee (2014) [23]. In their work, the micromechanical model is connected to the macroscopic scale by using the non-affine microsphere model. The model is able to quantitatively predict the macroscopic behavior of strain-crystallizing rubbers. Recently, Nateghi et al. (2018) [24] similarly proposed a micromechanical model that is incorporated into the affine microsphere model. In a contribution by Dargazany et al. (2014) [25], a micromechanical model for SIC in filled rubbers additionally considers inelastic properties of filled rubbers such as the Mullins effect, the permanent setting effect, and the induced anisotropy. In contrast to the aforementioned models that focus on the modeling of mechanical effects of the SIC, the publication by Behnke et al. (2018) [26] simulates the temperature and time dependency of SIC by taking the induced anisotropy into account. The newer publications by Khiêm and Itskov (2018) [27] and Khiêm et al. (2019) [28] focus on calorimetric effects, such as the temperature change due to SIC and the evolution of heat sources, whereas the work by Felder et al. (2019) [29] proposes a thermomechanically coupled model related to the crystallization by cooling from the melt.

Molecular dynamics is an alternative approach to modeling the microstructural phenomena of SIC, as presented in works by Nie (2015) [30] and Yamamoto (2019) [31]. This approach is suitable for investigations at nanoscale since it directly simulates effects related to interatomic potentials or thermal fluctuations of atoms. However, the approach is time-consuming and subject to a high computational effort. This strongly motivates the continuum mechanical modeling which intrinsically includes the aspects mentioned above and enables the efficient simulation of much larger samples than in the case of the molecular dynamics. Coupled with the multiscale strategies, a continuum mechanical model is even able to simulate practical applications.

As the previous overview shows, the development of measurement techniques has already made a significant contribution to the study of the SIC process. Nevertheless, there are still open issues which have not yet been sufficiently clarified by experimental studies, since phenomena are related to the nanoscale and are thus not accessible by the experimental techniques. Amongst others, the influence of the heat production on the shape, distribution, and mutual interaction of crystalline regions for high strain states can be pointed out. Moreover, already existing models mostly provide data on the effective material behavior without giving insight into the developed microstructure and the related temperature distribution. In contrast to these strategies, the present model treats the microstructural changes of the crystalline regions within the amorphous polymer matrix as well as the heat production of both phases. Our previous work [32] presents a material model which depicts the mechanical characteristics of the SIC phenomenon in agreement with the experimental investigations (Figure 1). Based on this model, the aim of the present contribution is to visualize the microstructural development and the change of temperature distribution within a material sample and to observe their interaction depending on external influences. The approach proposed focuses on the unfilled polymers being the representatives of nearly incompressible materials.

The paper is structured as follows: Section 2 gives an overview of general concepts used for the thermomechanical modeling of the SIC. It introduces the finite deformation kinematics that rely on the multiplicative split of the deformation gradient into an elastic, crystalline, and thermal part and assumes the Arruda–Boyce model for elastically stored energy. The same section investigates the thermodynamic consistency, determines the energetically conjugated pairs, and uses the minimum principle of dissipation potential to derive evolution equations. Subsequently, it also defines the boundary value problem (BVP) including the balance of linear momentum and balance of energy. Whereas Section 2 deals with the general concepts, Section 3 introduces specific assumptions related to the modeling of the SIC. Primarily, it defines the internal variables and proposes expressions for the Helmholtz free energy and dissipation potential. Both are chosen such that the resulting evolution equations simulate the increase and the decrease of the crystalline regions and of the temperature during a cyclic test. Finally, Section 4 provides details on the implementation of the material model into the FEM software FEAP, and Section 5 demonstrates the application. Selected numerical examples pertaining cyclic tensile loads visualize the microstructure evolution and the corresponding temperature distribution. At first, an academic example studies the material response of three crystalline regions with different initial regularity while growing and shrinking. Thereafter, the model is applied to monitor the microstructure evolution and temperature distribution for a sample with a complex initial configuration. The paper finishes with conclusions and an outlook.

2. Thermodynamically Consistent Framework

2.1. Kinematics within the Finite Thermomechanical Framework

The multiplicative decomposition of the deformation gradient $\mathit{F}$ into an elastic part ${\mathit{F}}^{\mathrm{e}}$, an inelastic part ${\mathit{F}}^{\mathrm{i}}$, and a thermal part ${\mathit{F}}^{\mathrm{th}}$

$$\mathit{F}={\mathit{F}}^{\mathrm{e}}·{\mathit{F}}^{\mathrm{i}}·{\mathit{F}}^{\mathrm{th}}$$

(1)

is an appropriate method to simultaneously incorporate inelastic and thermal effects within the finite deformation theory [33,34]. The deformation gradient $\mathit{F}$ is a two-point tensor that maps the initial configuration (${\mathcal{B}}_0$) to the deformed configuration (${\mathcal{B}}_t$). Formula (1), however, implies the existence of two additional intermediate configurations (Figure 2) such that ${\mathit{F}}^{\mathrm{th}}$ maps the initial to the thermal intermediate configuration (${\mathcal{B}}_{\mathrm{th}}$). The latter is then related to the intermediate configuration ${\mathcal{B}}_{\mathrm{i}}$ by the inelastic deformation gradient ${\mathit{F}}^{\mathrm{i}}$. Eventually, the elastic deformation gradient ${\mathit{F}}^{\mathrm{e}}$ maps the inelastic intermediate to the current configuration. It is also worth mentioning that the ordering of intermediate configurations is not uniquely prescribed. Hartmann (2012) [35], for example, defines the multiplicative decomposition into the thermal and mechanical deformation gradient in a reverse order.

Based on assumption (1), the elastic deformation gradient is written as

$${\mathit{F}}^{\mathrm{e}}=\mathit{F}·{{\mathit{F}}^{\mathrm{th}}}^{-1}·{{\mathit{F}}^{\mathrm{i}}}^{-1}$$

(2)

which furthermore yields the relationship for the corresponding rate

$${\dot{\mathit{F}}}^{\mathrm{e}}=\dot{\mathit{F}}·{{\mathit{F}}^{\mathrm{th}}}^{-1}·{{\mathit{F}}^{\mathrm{i}}}^{-1}+\mathit{F}·\dot{\overline{{{\mathit{F}}^{\mathrm{th}}}^{-1}}}·{{\mathit{F}}^{\mathrm{i}}}^{-1}+\mathit{F}·{{\mathit{F}}^{\mathrm{th}}}^{-1}·\dot{\overline{{{\mathit{F}}^{\mathrm{i}}}^{-1}}}.$$

(3)

The relationship above depends on the time derivatives $\dot{\overline{{{\mathit{F}}^{\mathrm{i}}}^{-1}}}$ and $\dot{\overline{{{\mathit{F}}^{\mathrm{th}}}^{-1}}}$, which can be expressed in a more appropriate manner as will be shown by the example of $\dot{\overline{{{\mathit{F}}^{\mathrm{i}}}^{-1}}}$. The procedure starts by taking the time derivative of the identity tensor $\mathit{I}$, which is expressed as a function of ${\mathit{F}}^{\mathrm{i}}$

$${\mathit{F}}^{\mathrm{i}}·{{\mathit{F}}^{\mathrm{i}}}^{-1}=\mathit{I}\Rightarrow \dot{\overline{\left({\mathit{F}}^{\mathrm{i}}·{{\mathit{F}}^{\mathrm{i}}}^{-1}\right)}}={\dot{\mathit{F}}}^{\mathrm{i}}·{{\mathit{F}}^{\mathrm{i}}}^{-1}+{\mathit{F}}^{\mathrm{i}}·\dot{\overline{{{\mathit{F}}^{\mathrm{i}}}^{-1}}}=\mathbf{0}.$$

(4)

Rearranging Equation (4) then yields

$$\dot{\overline{{{\mathit{F}}^{\mathrm{i}}}^{-1}}}=-{{\mathit{F}}^{\mathrm{i}}}^{-1}·{\dot{\mathit{F}}}^{\mathrm{i}}·{{\mathit{F}}^{\mathrm{i}}}^{-1}.$$

(5)

A similar procedure provides an analogous relationship for $\dot{\overline{{{\mathit{F}}^{\mathrm{th}}}^{-1}}}$:

$$\dot{\overline{{{\mathit{F}}^{\mathrm{th}}}^{-1}}}=-{{\mathit{F}}^{\mathrm{th}}}^{-1}·{\dot{\mathit{F}}}^{\mathrm{th}}·{{\mathit{F}}^{\mathrm{rh}}}^{-1},$$

(6)

such that the insertion of Equations (2), (5), and (6) into (3) ends up in

$${\dot{\mathit{F}}}^{\mathrm{e}}=\dot{\mathit{F}}·{{\mathit{F}}^{\mathrm{th}}}^{-1}·{{\mathit{F}}^{\mathrm{i}}}^{-1}-\mathit{F}·{{\mathit{F}}^{\mathrm{th}}}^{-1}·{\dot{\mathit{F}}}^{\mathrm{th}}·{{\mathit{F}}^{\mathrm{th}}}^{-1}·{{\mathit{F}}^{\mathrm{i}}}^{-1}-{\mathit{F}}^{\mathrm{e}}·{\dot{\mathit{F}}}^{\mathrm{i}}·{{\mathit{F}}^{\mathrm{i}}}^{-1},$$

(7)

which is the final expression for the rate elastic deformations only depending on rates and on inverse values of basic quantities $\mathit{F}$, ${\mathit{F}}^{\mathrm{i}}$, and ${\mathit{F}}^{\mathrm{th}}$.

The nearly incompressible elasticity of polymeric materials is considered in works by Flory (1961) [36] and Doll et al. (2000) [37]. These authors introduce the decomposition of the deformation gradient ${\mathit{F}}^{\mathrm{e}}$ into a volmetric part ${\widehat{\mathit{F}}}^{\mathrm{e}}$ and a deviatoric part ${\overline{\mathit{F}}}^{\mathrm{e}}$, both depending on the determinant $J^{\mathrm{e}}$:

$${\mathit{F}}^{\mathrm{e}}={\widehat{\mathit{F}}}^{\mathrm{e}}·{\overline{\mathit{F}}}^{\mathrm{e}},{\widehat{\mathit{F}}}^{\mathrm{e}}:=J^{{\mathrm{e}}^{\frac{1}{3}}}\mathit{I},{\overline{\mathit{F}}}^{\mathrm{e}}:=J^{{\mathrm{e}}^{-\frac{1}{3}}}{\mathit{F}}^{\mathrm{e}},J^{\mathrm{e}}=\det \left({\mathit{F}}^{\mathrm{e}}\right).$$

(8)

For the isochoric deformations, it holds that $J^{\mathrm{e}}=1$ since the determinant $J^{\mathrm{e}}$ represents the measure of the elastic volumetric changes.

2.2. The Arruda–Boyce Elastic Energy

Nowadays, there are plenty of different models for the elastic energy part, as is shown in the overview by Steinmann et al. (2012) [38]. However, the present contribution assumes the Arruda–Boyce model relying on the statistical mechanics of idealized polymer chains. This micromechanical model is suitable for simulating the large elastic deformations and excellently fits to the chosen thermodynamic framework and numerical concept. Moreover, it only depends on two material parameters which are experimentally easily accessible and which are already known for many standard polymers such as natural rubber [39] and unfilled silicone rubber [40]. Since it is well established and explored, the Arruda–Boyce model represents an ideal platform for the present analysis where the focus is set on the investigation of dissipative aspects of the process rather than on the elastic effects.

The Arruda–Boyce energy ${\mathsf{\Psi }}^{\mathrm{AB}}$ originally has the form

$${\mathsf{\Psi }}^{\mathrm{AB}}\left({\overline{I}}_1\right)=\mu {\lambda }_m^2\left(\frac{{\lambda }_{\mathrm{chain}}}{{\lambda }_m}\beta +\ln \frac{\beta }{\sin \left(\beta \right)}\right),\beta =L^{-1}\left(\frac{{\lambda }_{\mathrm{chain}}}{{\lambda }_m}\right),{\lambda }_{\mathrm{chain}}=\sqrt{\frac{{\overline{I}}_1}{3}},$$

(9)

where ${\overline{I}}_1$ is the first invariant of the deviatoric elastic right Cauchy–Green tensor ${\overline{\mathit{C}}}^{\mathrm{e}}$

$${\overline{I}}_1=\mathrm{tr}\left({\overline{\mathit{C}}}^{\mathrm{e}}\right),{\overline{\mathit{C}}}^{\mathrm{e}}={\left({\overline{\mathit{F}}}^{\mathrm{e}}\right)}^T·{\overline{\mathit{F}}}^{\mathrm{e}}=J^{{\mathrm{e}}^{-\frac{2}{3}}}{\left({\mathit{F}}^{\mathrm{e}}\right)}^T·{\mathit{F}}^{\mathrm{e}}.$$

(10)

In Equation (9), $\mu$ denotes the shear modulus, ${\lambda }_m$ is the limiting network stretch, ${\lambda }_{\mathrm{chain}}$ is the chain stretch depending on the deviatoric first invariant, and $\beta$ denotes the inverse Langevin function which is related to the energy of a single random chain. The latter function cannot be expressed explicitly and is usually approximated by the Taylor series truncated up to the certain order [41,42]. The present work assumes that an approximation including three terms of the Taylor series provides sufficient accuracy. In this case, the Arruda–Boyce energy takes the form

$${\mathsf{\Psi }}^{\mathrm{AB}}\left({\overline{I}}_1\right)=\frac{\overline{\mu }}{2}\left[\left({\overline{I}}_1-3\right)+\frac{1}{10{\lambda }_m^2}\left({\overline{I}}_1^2-9\right)+\frac{11}{525{\lambda }_m^4}\left({\overline{I}}_1^3-27\right)\right],$$

(11)

where parameter $\overline{\mu }$ is determined from the consistency condition of model (11) with the linear elasticity theory for small strains. According to [43,44], this condition is expressed as

$${\left.\frac{\partial {\mathsf{\Psi }}^{\mathrm{AB}}}{\partial {\overline{I}}_1}\right|}_{{\overline{I}}_1=3}=\frac{\mu }{2},$$

(12)

and requires the evaluation of the derivative of Equation (11) at ${\overline{I}}_1=3$

$${\left.\frac{\partial {\mathsf{\Psi }}^{\mathrm{AB}}}{\partial {\overline{I}}_1}\right|}_{{\overline{I}}_1=3}=\frac{\overline{\mu }}{2}\left(1+\frac{3}{5{\lambda }_m^2}+\frac{99}{175{\lambda }_m^4}\right).$$

(13)

Finally, the implementation of Equation (13) in Equation (12) yields the expression for $\overline{\mu }$

$$\overline{\mu }=\mu {\left(1+\frac{3}{5{\lambda }_m^2}+\frac{99}{175{\lambda }_m^4}\right)}^{-1}.$$

(14)

Alternatively to the Taylor series, a range of Padé approximations can be applied for the numerical evaluation of the inverse Langevin function. These approximations have different degrees of accuracy and complexity as discussed in the review papers by Jedynak [45] and Carroll [46].

The original form of the Arruda–Boyce model corresponds to the incompressible material behavior; however, it can be easily extended to capture the influence of compressibility by adding a volumetric part ${\mathsf{\Psi }}^{\mathrm{vol}}$. In that case, the total elastic energy ${\mathsf{\Psi }}^{\mathrm{e}}$ turns into

$${\mathsf{\Psi }}^{\mathrm{e}}({\overline{I}}_1,J^{\mathrm{e}})={\mathsf{\Psi }}^{\mathrm{AB}}\left({\overline{I}}_1\right)+{\mathsf{\Psi }}^{\mathrm{vol}}\left(J^{\mathrm{e}}\right).$$

(15)

The chosen expression for the volumetric part of the energy is proposed by Simo and Taylor (1991) [47] and represents a special case of the Ogden model [48]

$$\begin{array}{cc}\hfill & {\mathsf{\Psi }}^{\mathrm{vol}}\left(J^{\mathrm{e}}\right)=\frac{K}{4}\left({\left(J^{\mathrm{e}}\right)}^2-1-2\ln \left(J^{\mathrm{e}}\right)\right),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & \frac{\partial {\mathsf{\Psi }}^{\mathrm{vol}}}{\partial J^{\mathrm{e}}}=\frac{K}{2}\left(J^{\mathrm{e}}-\frac{1}{J^{\mathrm{e}}}\right),\hfill \end{array}$$

(17)

where K denotes the bulk modulus. Assumption (16) is physically motivated with regard to the energy- and stress-free state in the reference configuration since it holds

$${\mathsf{\Psi }}^{\mathrm{vol}}(J^{\mathrm{e}}=1)=0,\frac{\partial {\mathsf{\Psi }}^{\mathrm{vol}}}{\partial J^{\mathrm{e}}}(J^{\mathrm{e}}=1)=0.$$

(18)

Moreover, the energy function (16) is convex and coercive in $J^{\mathrm{e}}$. The latter implies that the energy tends to infinity ${\mathsf{\Psi }}^{\mathrm{vol}}\to \infty$ for high compression and tension modes, i.e., for the cases where $J^{\mathrm{e}}\to 0$ and $J^{\mathrm{e}}\to \infty$. A comprehensive comparison of various models for the volumetric energy with respect to their mathematical and physical properties is shown by Hartmann and Neff (2003) [49], whereas the implementation of these approaches for nearly incompressible materials is discussed by Kadapa and Hossain (2020) [50].

2.3. Thermodynamic Consistency

The second thermodynamic law in the form of the Clausius–Duhem inequality [51] is chosen to investigate the thermodynamic consistency of the model:

$$\mathcal{D}=-\dot{e}+\mathsf{\Theta }\dot{\eta }+\frac{1}{{\rho }_0}\mathit{P}:\dot{\mathit{F}}+{\mathcal{D}}^{\mathrm{cond}}\ge 0,{\mathcal{D}}^{\mathrm{cond}}=-\frac{1}{{\rho }_0\mathsf{\Theta }}{\mathit{q}}_0·{\nabla }_{\mathit{X}}\mathsf{\Theta }.$$

(19)

Here, $\dot{e}$ is the rate of internal energy, $\dot{\eta }$ is the rate of entropy, $\mathcal{\Theta }$ is the temperature, ${\rho }_0$ denotes the density, $\mathit{P}:\dot{\mathit{F}}$ represents the internal power, $\mathit{P}$ is the first Piola-Kirchhoff stress tensor, $\mathit{X}$ are the coordinates in reference configuration, ${\mathcal{D}}^{\mathrm{cond}}$ is the heat conduction contribution to the dissipation $\mathcal{D}$, ${\mathit{q}}_0$ is the heat flux vector, and index “0” refers to the initial configuration. Equation (19) has the standard solution determining the heat transfer through the material. This solution is known as the Fourier law and has the following form in the reference configuration [52]

$${\mathit{q}}_0=-{\lambda }_{\theta }\det \left(\mathit{F}\right){\mathit{C}}^{-1}·{\nabla }_{\mathit{X}}\mathsf{\Theta },\mathit{C}={\mathit{F}}^T·\mathit{F},$$

(20)

where ${\lambda }_{\theta }$ is the thermal conductivity coefficient and $\mathit{C}$ is the right Cauchy–Green deformation tensor.

An alternative formulation of the dissipation inequality (19) is obtained by using the Legendre transformation of the internal energy yielding the Helmholtz free energy $\mathsf{\Psi }$

$$\mathsf{\Psi }=e-\mathsf{\Theta }\eta .$$

(21)

Consequently, the internal energy rate can be expressed as a function of elastic energy rate $\dot{\mathsf{\Psi }}$

$$\dot{e}=\dot{\mathsf{\Psi }}+\dot{\mathsf{\Theta }}\eta +\mathsf{\Theta }\dot{\eta }.$$

(22)

The particular case studied in this contribution assumes the Helmholtz energy $\mathsf{\Psi }$ as a function of three arguments: the elastic deformation gradient ${\mathit{F}}^{\mathrm{e}}$, temperature $\mathsf{\Theta }$, and an additional set of internal variables describing microstructural phenomena $\mathit{\gamma }$, such that its rate turns into

$$\dot{\mathsf{\Psi }}=\frac{\partial \mathsf{\Psi }}{\partial {\mathit{F}}^{\mathrm{e}}}:{\dot{\mathit{F}}}^{\mathrm{e}}+\frac{\partial \mathsf{\Psi }}{\partial \mathsf{\Theta }}\dot{\mathsf{\Theta }}+\frac{\partial \mathsf{\Psi }}{\partial \mathit{\gamma }}:\dot{\mathit{\gamma }}.$$

(23)

Finally, the insertion of Equations (7), (22), and (23) into (19) yields an alternative form of the dissipation inequality

$$\begin{array}{cc}\hfill \mathcal{D}& =\left(\frac{1}{{\rho }_0}\mathit{P}-\frac{\partial \mathsf{\Psi }}{\partial {\mathit{F}}^{\mathrm{e}}}·{\left({{\mathit{F}}^{\mathrm{i}}}^{-1}\right)}^T·{\left({{\mathit{F}}^{\mathrm{th}}}^{-1}\right)}^T\right):\dot{\mathit{F}}-\left(\frac{\partial \mathsf{\Psi }}{\partial \mathsf{\Theta }}+\eta \right)\dot{\mathsf{\Theta }}\hfill \\ \hfill & +\left({\left({{\mathit{F}}^{\mathrm{th}}}^{-1}\right)}^T·{\mathit{F}}^T·\frac{\partial \mathsf{\Psi }}{\partial {\mathit{F}}^{\mathrm{e}}}·{\left({{\mathit{F}}^{\mathrm{i}}}^{-1}\right)}^T\right):\left({\dot{\mathit{F}}}^{\mathrm{th}}·{{\mathit{F}}^{\mathrm{th}}}^{-1}\right)\hfill \\ \hfill & +\left({{\mathit{F}}^{\mathrm{e}}}^T·\frac{\partial \mathsf{\Psi }}{\partial {\mathit{F}}^{\mathrm{e}}}\right):\left({\dot{\mathit{F}}}^{\mathrm{i}}·{{\mathit{F}}^{\mathrm{i}}}^{-1}\right)-\frac{\partial \mathsf{\Psi }}{\partial \mathit{\gamma }}:\dot{\mathit{\gamma }}+{\mathcal{D}}^{\mathrm{cond}}\ge 0,\hfill \end{array}$$

(24)

which is suitable for the application of the Coleman–Noll procedure and for determining the constitutive laws [53]. Following the procedure mentioned, the first two terms in Equation (24) yield definitions for the first Piola–Kirchhoff stress tensor and entropy

$$\begin{array}{cc}\hfill & \mathit{P}={\rho }_0\frac{\partial \mathsf{\Psi }}{\partial {\mathit{F}}^{\mathrm{e}}}·{\left({{\mathit{F}}^{\mathrm{i}}}^{-1}\right)}^T·{\left({{\mathit{F}}^{\mathrm{th}}}^{-1}\right)}^T,\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & \eta =-\frac{\partial \mathsf{\Psi }}{\partial \mathsf{\Theta }},\hfill \end{array}$$

(26)

whereas the remaining terms enable the definition of thermodynamically conjugated pairs. The velocity gradient and Mandel stress tensor related to the thermal intermediate configuration are the first conjugated pair

$$\begin{array}{cc}\hfill & {\mathit{L}}^{\mathrm{th}}:={\dot{\mathit{F}}}^{\mathrm{th}}·{{\mathit{F}}^{\mathrm{th}}}^{-1},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & {\mathit{M}}^{\mathrm{th}}:={\left({{\mathit{F}}^{\mathrm{th}}}^{-1}\right)}^T·{\mathit{F}}^T·\frac{\partial \mathsf{\Psi }}{\partial {\mathit{F}}^{\mathrm{e}}}·{\left({{\mathit{F}}^{\mathrm{i}}}^{-1}\right)}^T\hfill \end{array}$$

(28)

and velocity gradient and Mandel stress tensor related to the inelastic intermediate configuration represent the second conjugated pair

$$\begin{array}{cc}\hfill & {\mathit{L}}^{\mathrm{i}}:={\dot{\mathit{F}}}^{\mathrm{i}}·{{\mathit{F}}^{\mathrm{i}}}^{-1},\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill & {\mathit{M}}^{\mathrm{i}}:={{\mathit{F}}^{\mathrm{e}}}^T·\frac{\partial \mathsf{\Psi }}{\partial {\mathit{F}}^{\mathrm{e}}}.\hfill \end{array}$$

(30)

In addition, both Mandel stress tensors (Equations (28) and (30)) can be related to each other by considering the multiplicative decomposition of the deformation gradient (Equation (1))

$${\mathit{M}}^{\mathrm{th}}={\left({{\mathit{F}}^{\mathrm{th}}}^{-1}\right)}^T·{{\mathit{F}}^{\mathrm{th}}}^T·{{\mathit{F}}^{\mathrm{i}}}^T·{{\mathit{F}}^{\mathrm{e}}}^T·\frac{\partial \mathsf{\Psi }}{\partial {\mathit{F}}^{\mathrm{e}}}·{\left({{\mathit{F}}^{\mathrm{i}}}^{-1}\right)}^T={{\mathit{F}}^{\mathrm{i}}}^T·{\mathit{M}}^{\mathrm{i}}·{\left({{\mathit{F}}^{\mathrm{i}}}^{-1}\right)}^T.$$

(31)

According to the previous relationship, the Mandel stress tensor is mapped from the intermediate inelastic configuration to the intermediate thermal configuration, which corresponds to a pullback operation.

Eventually, the last term in inequality (24) provides the definition for the driving force corresponding to the rate of internal variables $\dot{\mathit{\gamma }}$

$$q_{\gamma }:=-\frac{\partial \mathsf{\Psi }}{\partial \mathit{\gamma }}.$$

(32)

Using the new notation (Equations (27)–(30) and (32)), the dissipation inequality (24) is written in a shorter form

$$\mathcal{D}={\mathcal{D}}^{\mathrm{th}}+{\mathcal{D}}^{\mathrm{i}}+{\mathcal{D}}^{\gamma }+{\mathcal{D}}^{\mathrm{cond}}\ge 0,$$

(33)

which includes the contribution due to thermal deformations ${\mathcal{D}}^{\mathrm{th}}$, the contribution due to the inelastic deformations ${\mathcal{D}}^{\mathrm{i}}$, the contribution due to microstructural changes ${\mathcal{D}}^{\gamma }$, and the contribution due to the heat conduction ${\mathcal{D}}^{\mathrm{cond}}$. The single contributions are defined as follows:

$${\mathcal{D}}^{\mathrm{th}}={\mathit{M}}^{\mathrm{th}}:{\mathit{L}}^{\mathrm{th}}\ge 0,{\mathcal{D}}^{\mathrm{i}}={\mathit{M}}^{\mathrm{i}}:{\mathit{L}}^{\mathrm{i}}\ge 0,{\mathcal{D}}^{\gamma }={\mathit{q}}_{\gamma }:\dot{\mathit{\gamma }}\ge 0,$$

(34)

where the non-negativity of each term is required separately.

2.4. Derivation of Evolution Equations

Evolution equations for internal variables are an essential part of material models for inelastic processes. Their derivation is a challenging task such that different approaches have been established for this purpose [54]. The present contribution follows the minimum principle of dissipation potential [55] which is expressed as follows

$$min\{\mathcal{L}=\dot{\mathsf{\Psi }}+\Delta |\dot{\mathit{\gamma }}\}.$$

(35)

According to this principle, the minimization of the Lagrangian $\mathcal{L}$ including the Helmholtz energy rate $\dot{\mathsf{\Psi }}$ (Equation (23)) and the dissipation potential $\Delta$ leads to the evolution laws for the internal variables $\mathit{\gamma }$. Bearing in mind that the Helmholtz energy $\mathsf{\Psi }$ is a function of three arguments ${\mathit{F}}^{\mathrm{e}}$, $\mathsf{\Theta }$, and $\mathit{\gamma }$, the Lagrange function can be written in an extended form

$$\mathcal{L}=\frac{\partial \mathsf{\Psi }}{\partial {\mathit{F}}^{\mathrm{e}}}:{\dot{\mathit{F}}}^{\mathrm{e}}+\frac{\partial \mathsf{\Psi }}{\partial \mathsf{\Theta }}\dot{\mathsf{\Theta }}+\frac{\partial \mathsf{\Psi }}{\partial \mathit{\gamma }}:\dot{\mathit{\gamma }}+\Delta \left(\dot{\mathit{\gamma }}\right).$$

(36)

Now, the minimization of Equation (36) with respect to $\dot{\gamma }$ reads the conditions

$$\frac{\partial \mathcal{L}}{\partial \dot{\mathit{\gamma }}}=\frac{\partial \mathsf{\Psi }}{\partial \mathit{\gamma }}+\frac{\partial \Delta }{\partial \dot{\mathit{\gamma }}}=\mathbf{0}\Rightarrow \frac{\partial \Delta }{\partial \dot{\mathit{\gamma }}}=-\frac{\partial \mathsf{\Psi }}{\partial \mathit{\gamma }},$$

(37)

and a comparison of Equations (32) and (37) leads to the conclusion that the driving force of an internal variable is equal to the derivative of the dissipation potential with respect to the same quantity

$${\mathit{q}}_{\mathit{\gamma }}=\frac{\partial \Delta }{\partial \dot{\mathit{\gamma }}}.$$

(38)

Equations (35)–(38) define a generic procedure which will be used to derive equations driving the microstructural changes associated with the SIC (Section 3.5).

2.5. Balance Equations

Different from a BVP, which is related to a purely mechanical problem, the BVP corresponding to a thermomechanical problem includes two differential equations, namely the balance of linear momentum and the balance of energy. These differential equations are accompanied by the suitable boundary conditions as follows [56,57]:

$$\begin{array}{cc}\hfill \mathrm{Div}\mathit{P}+{\rho }_0\mathit{b}=\mathbf{0}& \mathrm{in}\mathcal{B},\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill \dot{e}+\frac{1}{{\rho }_0}\mathrm{Div}\left({\mathit{q}}_0\right)=\frac{1}{{\rho }_0}\mathit{P}:\dot{\mathit{F}}+r_{\mathsf{\Theta }}& \mathrm{in}\mathcal{B},\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill \mathit{P}·\mathit{n}=\overline{\mathit{t}}& \mathrm{on}\partial {\mathcal{B}}^t,\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill \mathit{u}=\overline{\mathit{u}}& \mathrm{on}\partial {\mathcal{B}}^u,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill -{\mathit{q}}_0·\mathit{n}={\overline{q}}_{\mathsf{\Theta }}& \mathrm{on}\partial {\mathcal{B}}^q,\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill \mathsf{\Theta }=\overline{\mathsf{\Theta }}& \mathrm{on}\partial {\mathcal{B}}^{\mathsf{\Theta }}.\hfill \end{array}$$

(44)

Here, $\mathit{b}$ is the body force, $r_{\mathsf{\Theta }}$ is the heat source, ${\mathit{q}}_0$ is the heat flux vector defined by the Fourier law (20), and $\mathit{n}$ is the surface normal vector. The balance of linear momentum (39) is supplemented by the Neumann and Dirichlet boundary conditions expressed in terms of prescribed tractions $\overline{\mathit{t}}$, displacements $\mathit{u}$ and prescribed displacements $\overline{\mathit{u}}$. The corresponding boundary parts are denoted by ${\mathcal{B}}^t$ and ${\mathcal{B}}^u$, respectively. Similarly, the balance of energy (40) is accompanied by boundary conditions relating the heat flux ${\mathit{q}}_0$ and temperature $\mathsf{\Theta }$ to the prescribed values ${\overline{q}}_{\mathsf{\Theta }}$ and $\overline{\mathsf{\Theta }}$ acting on the boundary parts ${\mathcal{B}}^q$ and ${\mathcal{B}}^{\mathsf{\Theta }}$.

Internal energy rate $\dot{e}$ in Equation (40) is defined by Equation (22) and thus requires a more precise study of entropy rate $\dot{\eta }$. To this end, definition (26) is used along with the fact that the free energy is a function of three arguments:

$$\dot{\eta }=-\dot{\overline{\left(\frac{\partial \mathsf{\Psi }}{\partial \mathsf{\Theta }}\right)}}=-\left(\frac{{\partial }^2\mathsf{\Psi }}{\partial \mathsf{\Theta }\partial {\mathit{F}}^{\mathrm{e}}}:{\dot{\mathit{F}}}^{\mathrm{e}}+\frac{{\partial }^2\mathsf{\Psi }}{\partial {\mathsf{\Theta }}^2}\dot{\mathsf{\Theta }}+\frac{{\partial }^2\mathsf{\Psi }}{\partial \mathsf{\Theta }\partial \mathit{\gamma }}:\dot{\mathit{\gamma }}\right).$$

(45)

Now, taking into account the definitions introduced in Section 2.3, the internal energy rate turns into

$$\dot{e}=\frac{1}{{\rho }_0}\mathit{P}:\dot{\mathit{F}}-{\mathit{M}}^{\mathrm{th}}:{\mathit{L}}^{\mathrm{th}}-{\mathit{M}}^{\mathrm{i}}:{\mathit{L}}^{\mathrm{i}}-{\mathit{q}}_{\gamma }:\dot{\mathit{\gamma }}+\mathsf{\Theta }\dot{\eta }$$

(46)

such that the final form of the energy balance is obtained by inserting Equations (45) and (46) into (40)

$$\begin{array}{cc}\hfill & c_d\dot{\mathsf{\Theta }}+\frac{1}{{\rho }_0}\mathrm{Div}\left({\mathit{q}}_0\right)=r_{\mathsf{\Theta }}+{\mathit{M}}^{\mathrm{th}}:{\mathit{L}}^{\mathrm{th}}+{\mathit{M}}^{\mathrm{i}}:{\mathit{L}}^{\mathrm{i}}+{\mathit{H}}_{\gamma }:\dot{\mathit{\gamma }}+{\mathit{H}}_F:{\dot{\mathit{F}}}^{\mathrm{e}},\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill & c_d:=-\mathsf{\Theta }\frac{{\partial }^2\mathsf{\Psi }}{\partial {\mathsf{\Theta }}^2},{\mathit{H}}_{\gamma }:={\mathit{q}}_{\gamma }+\mathsf{\Theta }\frac{{\partial }^2\mathsf{\Psi }}{\partial \mathsf{\Theta }\partial \mathit{\gamma }},{\mathit{H}}_F:=\mathsf{\Theta }\frac{{\partial }^2\mathsf{\Psi }}{\partial \mathsf{\Theta }\partial {\mathit{F}}^{\mathrm{e}}}.\hfill \end{array}$$

(48)

The standard notation introduced in (48) includes the heat capacity $c_d$, the latent heat due to internal microstructural processes in the material ${\mathit{H}}_{\gamma }$ and the latent heat due to deformations ${\mathit{H}}_F$. The latter enables, amongst others, the modeling of the Gough–Joule effect [58].

3. Thermomechanical Modeling of the SIC

Whereas Section 2 provides the general framework for modeling thermomechanical processes, Section 3 presents the assumptions specific to the modeling of the SIC. This includes the choice of internal variables, the assumption for the Helmholtz free energy, the definition of coupling conditions, and the formulation of the dissipation potential.

3.1. Definition of Internal Variables

The multiplicative decomposition of the deformation gradient representing the base of kinematics in the finite thermomechanical framework (Equation (1)) already incorporates two internal variables, namely the inelastic deformation gradient ${\mathit{F}}^{\mathrm{i}}$ and its thermal counterpart ${\mathit{F}}^{\mathrm{th}}$. In the case of the SIC, the inelastic deformations are caused by the crystallization, which will be emphasized by introducing index “c” instead of index “i” in the following text (i.e., ${\mathit{F}}^{\mathrm{c}}:={\mathit{F}}^{\mathrm{i}}$, ${\mathit{L}}^{\mathrm{c}}:={\mathit{L}}^{\mathrm{i}}$, ${\mathit{M}}^{\mathrm{c}}:={\mathit{M}}^{\mathrm{i}}$ and ${\mathcal{D}}^{\mathrm{c}}:={\mathcal{D}}^{\mathrm{i}}$ ).

However, the description of microstructural processes related to SIC requires the introduction of two additional internal variables: regularity of polymer chain network ($\chi$) and thermal flexibility of polymer network ($\alpha$). The regularity in this context comprises the information on the ordering of the polymer chains with respect to each other and the degree of order among the polymer atoms. It takes values from the range $[0,1]$, such that values close to zero correspond to an amorphous state, whereas values close to the value of one are classified as crystalline regions. During a tensile test, the regularity evolves, thus simulating the formation/degradation of crystalline regions. The second internal variable captures physical properties related to the molecular mobility in the polymer which is governed by the chain flexibility and the temperature [59]. The flexibility of the polymer chains in the network changes due to warming up and cooling during a cyclic tensile test, where heating induces the chains to become more flexible.

3.2. Assumption for the Helmholtz Free Energy Density

In order to obtain a model suitable for simulating the SIC, the Helmholtz free energy is assumed to consist of four terms combining the benefits of several modeling approaches

$$\mathsf{\Psi }={\mathsf{\Psi }}^{\mathrm{e}}({\overline{I}}_1,J^{\mathrm{e}})+{\mathsf{\Psi }}^{\mathrm{th}}\left(\mathsf{\Theta }\right)+{\mathsf{\Psi }}^{\mathrm{crys}}\left(\chi \right)+{\mathsf{\Psi }}^{\mathrm{coup}}(\alpha ,\mathsf{\Theta }).$$

(49)

Here, the elastically stored energy ${\mathsf{\Psi }}^{\mathrm{e}}$ (Equation (15)) includes the Arruda–Boyce model (Equation (11)) and the volumetric contribution (Equation (16)). A similar combination has been used in works by Gasser and Holzapfel (2002) [60] and by Elguedj and Hughes (2014) [61] when studying nearly incompressible materials behavior. The thermally stored energy part ${\mathsf{\Psi }}^{\mathrm{th}}$ has been proposed by Raniecki and Bruhns (1991) [62] and has recently been used by Mahnken (2013) [52] for thermomechanical modeling of polymers

$${\mathsf{\Psi }}^{\mathrm{th}}\left(\mathsf{\Theta }\right)=c_d\left(\mathsf{\Theta }-{\mathsf{\Theta }}_0\right)-c_d\mathsf{\Theta }ln\left(\frac{\mathsf{\Theta }}{{\mathsf{\Theta }}_0}\right).$$

(50)

The thermal contribution above is a function of reference temperature ${\mathsf{\Theta }}_0$ and holds under the condition that the heat capacity $c_d$ is independent of temperature. The crystalline energy part ${\mathsf{\Psi }}^{\mathrm{crys}}$ is proposed in [32] to model the regularity evolution during the unloading phase. It is a linear function

$${\mathsf{\Psi }}^{\mathrm{crys}}\left(\chi \right)=c_1\chi ,$$

(51)

where parameter $c_1$ is load-dependent and enables to control the evolution direction. The last energy contribution, ${\mathsf{\Psi }}^{\mathrm{coup}}$, is a mixed term in thermal flexibility and temperature, including proportionality constant $c_2$

$${\mathsf{\Psi }}^{\mathrm{coup}}(\alpha ,\mathsf{\Theta })=c_2\alpha \left(\mathsf{\Theta }-{\mathsf{\Theta }}_0\right).$$

(52)

A similar bilinear form is suggested in the work by Hajidehi and Stupkiewicz (2017) [63] for the energy transition in shape memory alloys with the aim to simulate the temperature decrease during unloading.

3.3. Coupling Conditions

The problem formulation presented depends on four internal variables, two of which have a tensorial character. However, some additional information is provided by the physics of the SIC phenomenon which enables the simplification of the underlying setup. The focus is first set on the velocity gradient ${\mathit{L}}^{\mathrm{c}}$ referring to the rate of deformations caused by the crystallization. The closer specification of this quantity relies on the experimental observation, showing that the orientation of crystalline regions is dependent on the external load orientation. This implies that velocity gradient ${\mathit{L}}^{\mathrm{c}}$ is coaxial with the Mandel stress tensor ${\mathit{M}}^{\mathrm{c}}$ since they are a thermodynamically conjugated pair. Furthermore, it can be expected that the intensity of deformations due to the crystallization depends on the regularity degree of network, which leads to the following assumption for the velocity gradient

$${\mathit{L}}^{\mathrm{c}}=k_1\dot{\chi }\frac{{\mathit{M}}^{\mathrm{c},\mathrm{dev}}}{\parallel {\mathit{M}}^{\mathrm{c},\mathrm{dev}}\parallel },{\mathit{M}}^{\mathrm{c},\mathrm{dev}}={\mathit{M}}^{\mathrm{c}}-\frac{\mathrm{tr}\left({\mathit{M}}^{\mathrm{c}}\right)}{3}\mathit{I}.$$

(53)

Here, $k_1$ denotes the proportionality constant and the choice of the deviatoric part of the Mandel stress tensor is substantiated by the fact that the model deals with the unfilled polymers being the representatives of nearly incompressible materials. In a similar way, the thermal velocity gradient and the rate of thermal flexibility are coupled to each other through the linear relationship with the positive proportionality constant $k_2$

$${\mathit{L}}^{\mathrm{th}}=k_2\dot{\alpha }\mathit{I}.$$

(54)

Equation (54) suggests a diagonal form of ${\mathit{L}}^{\mathrm{th}}$ since an isotropic character of thermal effects is expected. After introducing coupling conditions (53) and (54), the problem formulation only depends on the two scalar internal variables $\chi$ and $\alpha$.

3.4. Derivation of Driving Forces

Assumptions made in Section 3.2 and Section 3.3 enable the specification of driving forces corresponding to remaining internal variables $\chi$ and $\alpha$. For this purpose, coupling conditions (53) and (54) are first introduced into the dissipation inequality (33)

$$\begin{array}{cc}\hfill \mathcal{D}& ={\mathit{M}}^{\mathrm{th}}:{\mathit{L}}^{\mathrm{th}}+{\mathit{M}}^{\mathrm{c}}:{\mathit{L}}^{\mathrm{c}}+q_{\chi }\dot{\chi }+q_{\alpha }\dot{\alpha }+{\mathcal{D}}^{\mathrm{cond}}\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill & =k_2\mathrm{tr}\left({\mathit{M}}^{\mathrm{th}}\right)\dot{\alpha }+k_1∥{\mathit{M}}^{\mathrm{c},\mathrm{dev}}∥\dot{\chi }+q_{\chi }\dot{\chi }+q_{\alpha }\dot{\alpha }+{\mathcal{D}}^{\mathrm{cond}}.\hfill \end{array}$$

(56)

Thereafter, $q_{\chi }$ and $q_{\alpha }$ are calculated according to the standard definitions applied to the assumptions on energy contributions (51) and (52)

$$q_{\chi }=-\frac{\partial {\mathsf{\Psi }}^{\mathrm{crys}}}{\partial \chi }=-c_1,q_{\alpha }=-\frac{\partial {\mathsf{\Psi }}^{\mathrm{coup}}}{\partial \alpha }=-c_2\left(\mathsf{\Theta }-{\mathsf{\Theta }}_0\right).$$

(57)

Now, the implementation of Equation (57) in (56) and gathering of terms including the same type of rates yield expressions for sought driving forces

$$\begin{array}{cc}\hfill & {\tilde{q}}_{\chi }:=k_1∥{\mathit{M}}^{\mathrm{c},\mathrm{dev}}∥-c_1,\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill & {\tilde{q}}_{\alpha }:=k_2\mathrm{tr}\left({\mathit{M}}^{\mathrm{th}}\right)-c_2\left(\mathsf{\Theta }-{\mathsf{\Theta }}_0\right),\hfill \end{array}$$

(59)

such that the dissipation inequality turns into

$$\mathcal{D}={\tilde{q}}_{\chi }\dot{\chi }+{\tilde{q}}_{\alpha }\dot{\alpha }+{\mathcal{D}}^{\mathrm{cond}}\ge 0.$$

(60)

The requirement for the non-negativity of each addend yields the conclusion that a driving force must have the same sign as the rate of the corresponding internal variable. This furthermore implies that a decrease of the regularity ($\dot{\chi }<0$) is accompanied by a negative driving force (${\tilde{q}}_{\chi }<0$) during the unloading stage. It is furthermore assumed that the rate of driving force (58)

$${\dot{\tilde{q}}}_{\chi }=k_1\dot{\overline{∥{\mathit{M}}^{\mathrm{c},\mathrm{dev}}∥}}$$

(61)

distinguishes the loading stage (${\dot{\tilde{q}}}_{\chi }\ge 0$) and the unloading stage (${\dot{\tilde{q}}}_{\chi }<0$).

3.5. Assumption for the Dissipation Potential and Derivation of Evolution Equations

The assumed dissipation potential includes two contributions, each of them dependent on the rate of a single internal variable

$$\Delta (\dot{\chi },\dot{\alpha })={\Delta }^{\mathrm{crys}}\left(\dot{\chi }\right)+{\Delta }^{\mathrm{th}}\left(\dot{\alpha }\right).$$

(62)

The choice of the potential part related to the network regularity ${\Delta }^{\mathrm{crys}}$ is motivated by the experimental results shown in Figure 1a. According to this diagram, the crystalline regions start to evolve after exceeding a threshold value and grow up to the end of the loading phase where the limiting stretch of $\approx 600\%$ must not be exceeded. Subsequently, the crystalline regions shrink during the unloading with a lower rate. The dissipation potential suitable for modeling this process has been proposed in [32] and has the following form

$${\Delta }^{\mathrm{crys}}\left(\dot{\chi }\right)=\left(A+B\right)\left|\dot{\chi }\right|.$$

(63)

Here, A is a crystallization limit for the evolution of the regularity and B represents the increment of this limit. The rate of increment B is a function of the regularity rate

$$\dot{B}=\frac{b}{f\left(\chi \right)}\left|\dot{\chi }\right|,f\left(\chi \right)={\beta }_1-{\left(\chi -{\beta }_2\right)}^{{\beta }_3},b=\left\{\begin{array}{cc}{b}_1\hfill & \mathrm{if}{\dot{\tilde{q}}}_{\chi }\ge 0,\hfill \\ b_2\hfill & \mathrm{if}{\dot{\tilde{q}}}_{\chi }<0.\hfill \end{array}\right.$$

(64)

The chosen formulation depends on material constants ${\beta }_1$, ${\beta }_2$, and ${\beta }_3$ and introduces function $f>0$ to favor the evolution of the regularity in regions with higher values. Moreover, the load-dependent parameter b controls the velocity of the evolution. The values for this parameter are different for the loading phase (${\dot{\tilde{q}}}_{\chi }\ge 0$) and unloading phase (${\dot{\tilde{q}}}_{\chi }<0$). Condition $b_2>b_1>0$ indicates a slower decrease of the regularity during the unloading phase than its growth during the loading phase.

The second part of the dissipation potential, ${\Delta }^{\mathrm{th}}$, is defined in terms of the rate of thermal flexibility and depends on material constants $D_1$ and $D_2$

$${\Delta }^{\mathrm{th}}\left(\dot{\alpha }\right)=\frac{D_1}{2}{\chi }^{-D_2}{\dot{\alpha }}^2.$$

(65)

The proposed dissipation potentials now serve as a basis for the application of the minimum principle (Section 2.4). To this end, the notation introduced in Section 3.3 is used to reformulate the rate of Helmholtz energy $\dot{\mathsf{\Psi }}$ (Equation (23)) and the Lagrange function $\mathcal{L}$ (Equation (36))

$$\begin{array}{c}\dot{\mathsf{\Psi }}=\mathit{P}:\dot{\mathit{F}}-\eta \dot{\mathsf{\Theta }}-{\tilde{q}}_{\chi }\dot{\chi }-{\tilde{q}}_{\alpha }\dot{\alpha },\end{array}$$

(66)

$$\begin{array}{c}\mathcal{L}=\mathit{P}:\dot{\mathit{F}}-\eta \dot{\mathsf{\Theta }}-{\tilde{q}}_{\chi }\dot{\chi }-{\tilde{q}}_{\alpha }\dot{\alpha }+\left(A+B\right)\left|\dot{\chi }\right|+\frac{D_1}{2}{\chi }^{-D_2}{\dot{\alpha }}^2.\end{array}$$

(67)

Subsequently, the minimization procedure yields the modified driving forces

$$\begin{array}{cc}\hfill & {\tilde{q}}_{\chi }=\frac{\partial {\Delta }^{\mathrm{crys}}}{\partial \dot{\chi }},\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill & {\tilde{q}}_{\alpha }=\frac{\partial {\Delta }^{\mathrm{th}}}{\partial \dot{\alpha }}.\hfill \end{array}$$

(69)

However, the evaluation of the driving force ${\tilde{q}}_{\chi }$ (Equation (68)) is not straightforward since potential ${\Delta }^{\mathrm{crys}}$ (Equation (63)) depends on the absolute value function which is not differentiable at $\dot{\chi }=0$. For that reason, the subdifferential of the crystalline dissipation potential $\partial {\Delta }^{\mathrm{crys}}\left(\dot{\chi }\right)$ is introduced instead of its derivative

$$\partial {\Delta }^{\mathrm{crys}}\left(\dot{\chi }\right)=\{\begin{array}{}\mathrm{(70)}& {\tilde{q}}_{\chi }= \left(A+B\right)\frac{\dot{\chi }}{\left|\dot{\chi }\right|}\mathrm{for}\dot{\chi }\ne 0,\mathrm{(71)}& \{{\tilde{q}}_{\chi };\left|{\tilde{q}}_{\chi }\right| \le A+B\}\mathrm{for}\dot{\chi }=0.\end{array}$$

Equation (70) represents the derivative of the potential with respect to the regularity rate under the condition that $\dot{\chi }\ne 0$, whereas any value $\left|{\tilde{q}}_{\chi }\right|\le A+B$ can be a solution for $\dot{\chi }=0$. Case (71) applies at the beginning of the tensile test where the material undergoes purely elastic deformations without changing the regularity. After exceeding the crystallization limit, the regularity starts to increase according to the evolution law which is obtained by rearranging Equation (70)

$$\dot{\chi }=\frac{\left|\dot{\chi }\right|}{A+B}{\tilde{q}}_{\chi }=\lambda {\tilde{q}}_{\chi },\lambda \ge 0.$$

(72)

The procedure for determining parameter $\lambda$ has been comprehensively explained in [32]. Within the present contribution, only the main steps are summarized as follows: First, the absolute values are taken from both sides of the Equation (70) and the obtained relationship is squared. In a second step, after taking the time derivative of modified Equation (70) and inserting Equations (58), (61), (64) and (72), the crystallization parameter turns into

$$\lambda =\frac{f{\tilde{q}}_{\chi }{\dot{\tilde{q}}}_{\chi }}{b{\left(A+B\right)}^2}\ge 0.$$

(73)

Equation (70) also holds for the complete unloading stage, where the degradation of crystalline regions occurs. This is achieved by introducing a shift yielding a negative driving force ${\tilde{q}}_{\chi }$ during the unloading (Section 4.2). Simultaneously, the rate ${\dot{\tilde{q}}}_{\chi }$ is negative per definition during the unloading stage (Equation (61)). Since driving force ${\tilde{q}}_{\chi }$ and its rate ${\dot{\tilde{q}}}_{\chi }$ are both negative during the unloading phase, the condition $\lambda \ge 0$ (Equation (73)) holds and relationship (72) applies.

The study of thermal influences accompanying crystallization are based on Equation (69) yielding the following evolution equation

$${\tilde{q}}_{\alpha }=\frac{\partial {\Delta }^{\mathrm{th}}}{\partial \dot{\alpha }}=D_1{\chi }^{-D_2}\dot{\alpha }\Rightarrow \dot{\alpha }=\frac{{\chi }^{D_2}}{D_1}{\tilde{q}}_{\alpha }.$$

(74)

Obviously, the rate $\dot{\alpha }$ is proportional to its driving force which characterizes the viscous type of evolution. In this case, the decrease of the temperature during the unloading phase is controlled by a suitable choice of constant $c_2$ in the driving force (Equation (59)).

Finally, the insertion of evolution Equations (72) and (74) in (60) yields the single dissipation contributions

$${\mathcal{D}}^{\chi }={\tilde{q}}_{\chi }\dot{\chi }=\lambda {\tilde{q}}_{\chi }^2\ge 0,{\mathcal{D}}^{\alpha }={\tilde{q}}_{\alpha }\dot{\alpha }=\frac{{\chi }^{D_2}}{D_1}{\tilde{q}}_{\alpha }^2\ge 0,$$

(75)

where both terms fulfill the non-negativity requirement, which approves the thermodynamical consistency of the model.

4. Numerical Implementation of the SIC Model

4.1. Implementation of the Thermomechanical Coupled Problem into the FEM

Balance Equations (39) and (47) are numerically solved by applying the FEM where a setup appropriate for nonlinear materials and large deformations is developed. This includes the transformation of the energy balance and the development of the weak form of the problem.

The balance of energy (47) is more closely specified by using driving forces (Equations (58) and (59)) along with coupling conditions (53) and (54)

$$c_d\dot{\mathsf{\Theta }}+\frac{1}{{\rho }_0}\mathrm{Div}\left({\mathit{q}}_0\right)=r_{\mathsf{\Theta }}+{\tilde{H}}_{\chi }\dot{\chi }+{\tilde{H}}_{\alpha }\dot{\alpha }+{\mathit{H}}_F:{\dot{\mathit{F}}}^{\mathrm{e}}.$$

(76)

Here, latent heat ${\tilde{H}}_{\chi }$ and ${\tilde{H}}_{\alpha }$ are expressed according to the definition (48)

$${\tilde{H}}_{\chi }:={\tilde{q}}_{\chi }+\mathsf{\Theta }\frac{{\partial }^2\mathsf{\Psi }}{\partial \mathsf{\Theta }\partial \chi },{\tilde{H}}_{\alpha }:={\tilde{q}}_{\alpha }+\mathsf{\Theta }\frac{{\partial }^2\mathsf{\Psi }}{\partial \mathsf{\Theta }\partial \alpha }$$

(77)

whereas ${\mathit{H}}_F$ remains as defined in Equation (48). Furthermore, the energy derivatives in the latent heat contributions are calculated on the basis of assumptions (49)–(52)

$${\tilde{H}}_{\chi }={\tilde{q}}_{\chi },{\tilde{H}}_{\alpha }={\tilde{q}}_{\alpha }+\mathsf{\Theta }{c}_2,{\mathit{H}}_F=\mathbf{0}.$$

(78)

Provided that the maximal cooling of natural rubber due to the Gough–Joule effect is about $0.005$ K [64], the temperature change by latent heat ${\mathit{H}}_F$ is negligible compared to the temperature change due to the SIC. This motivates the choice of energy density $\mathsf{\Psi }$ (Equation (49)), where it holds $\frac{{\partial }^2\mathsf{\Psi }}{\partial \mathsf{\Theta }\partial {\mathit{F}}^{\mathrm{e}}}=\mathbf{0}$, and thus the corresponding latent heat contribution becomes ${\mathit{H}}_F=\mathbf{0}$. Accordingly, balance of energy (76) turns into

$$c_d\dot{\mathsf{\Theta }}+\frac{1}{{\rho }_0}\mathrm{Div}\left({\mathit{q}}_0\right)=r_{\mathsf{\Theta }}+{\tilde{q}}_{\chi }\dot{\chi }+\left({\tilde{q}}_{\alpha }+\mathsf{\Theta }{c}_2\right)\dot{\alpha }.$$

(79)

For the purpose of transforming the strong formulation of the problem into its weak form, two steps are performed: First, the strong form of each differential equation (Equations (39) and (79)) is multiplied by test functions $\delta \mathit{u}$ and $\delta \mathsf{\Theta }$, commonly referred to as virtual displacement and virtual temperature. Thereafter, the equations are integrated over domain $\mathcal{B}$. Finally, the integration by parts and subsequently the Gauss integration theorem are applied to the divergence term in Equation (79) to transform the volume integral over the body $\mathcal{B}$ into an integral over the surface $\partial \mathcal{B}$. The weak form of both balance equations reads

$${\int }_{\mathcal{B}}{\nabla }_{\mathit{X}}\delta \mathit{u}:\mathit{P}\mathrm{d}V-{\int }_{\mathcal{B}}\delta \mathit{u}·{\rho }_0\mathit{b}\mathrm{d}V-{\int }_{\partial {\mathcal{B}}^{\mathit{t}}}\delta \mathit{u}·\mathit{t}\mathrm{d}A=0,$$

(80)

$$\begin{array}{cc}\hfill & {\int }_{\mathcal{B}}\delta \mathsf{\Theta }{c}_d\dot{\mathsf{\Theta }}\mathrm{d}V+{\int }_{\partial \mathcal{B}}\delta \mathsf{\Theta }\frac{1}{{\rho }_0}{\mathit{q}}_0·\mathit{n}\mathrm{d}A-{\int }_{\mathcal{B}}{\nabla }_{\mathit{X}}\delta \mathsf{\Theta }·\frac{1}{{\rho }_0}{\mathit{q}}_0\mathrm{d}V-\hfill \\ \hfill & {\int }_{\mathcal{B}}\delta \mathsf{\Theta }{r}_{\mathsf{\Theta }}\mathrm{d}V-{\int }_{\mathcal{B}}\delta \mathsf{\Theta }{\tilde{q}}_{\chi }\dot{\chi }\mathrm{d}V-{\int }_{\mathcal{B}}\delta \mathsf{\Theta }\left({\tilde{q}}_{\alpha }+\mathsf{\Theta }{c}_2\right)\dot{\alpha }\mathrm{d}V=0.\hfill \end{array}$$

(81)

In a second step, domain $\mathcal{B}\approx {\bigcup }_{e=1}^{n_{\mathrm{el}}}{\mathcal{B}}^e$ is spatially disrcretized into a finite number of elements $n_{\mathrm{el}}$ and integrals are approximated by a sum of integrals over single elements ${\mathcal{B}}^e$. The same applies for the boundary parts:

$$\sum _{e=1}^{n_{\mathrm{el}}}\left[{\int }_{{\mathcal{B}}^e}{\nabla }_{\mathit{X}}\delta {\mathit{u}}^e:\mathit{P}\mathrm{d}V-{\int }_{{\mathcal{B}}^e}\delta {\mathit{u}}^e·{\rho }_0\mathit{b}\mathrm{d}V-{\int }_{\partial {\mathcal{B}}^{e\mathit{t}}}\delta {\mathit{u}}^e·\mathit{t}\mathrm{d}A\right]=0,$$

(82)

$$\begin{array}{cc}\hfill & \sum _{e=1}^{n_{\mathrm{el}}}[{\int }_{{\mathcal{B}}^e}\delta {\mathsf{\Theta }}^e{c}_d{\dot{\mathsf{\Theta }}}^e\mathrm{d}V+{\int }_{\partial {\mathcal{B}}^e}\delta {\mathsf{\Theta }}^e\frac{1}{{\rho }_0}{\mathit{q}}_0·\mathit{n}\mathrm{d}A-{\int }_{{\mathcal{B}}^e}{\nabla }_{\mathit{X}}\delta {\mathsf{\Theta }}^e·\frac{1}{{\rho }_0}{\mathit{q}}_0\mathrm{d}V-\hfill \\ \hfill & \phantom{\sum _{e=1}^{n_{\mathrm{el}}}[}{\int }_{{\mathcal{B}}^e}\delta {\mathsf{\Theta }}^e{r}_{\mathsf{\Theta }}\mathrm{d}V-{\int }_{{\mathcal{B}}^e}\delta {\mathsf{\Theta }}^e{\tilde{q}}_{\chi }\dot{\chi }\mathrm{d}V-{\int }_{{\mathcal{B}}^e}\delta {\mathsf{\Theta }}^e\left({\tilde{q}}_{\alpha }+\mathsf{\Theta }{c}_2\right)\dot{\alpha }\mathrm{d}V]=0.\hfill \end{array}$$

(83)

The field variables $\delta {\mathit{u}}^e$ and $\delta {\mathsf{\Theta }}^e$ and their gradients ${\nabla }_{\mathit{X}}\delta {\mathit{u}}^e$ and ${\nabla }_{\mathit{X}}\delta {\mathsf{\Theta }}^e$ refer to single elements and are interpolated by using an approximation with $n_{\mathrm{en}}$ support points

$$\begin{array}{cc}\hfill & \delta {\mathit{u}}^e\approx \sum _{A=1}^{n_{\mathrm{en}}}\delta {\mathit{u}}^{eA}{N}^A,\delta {\mathsf{\Theta }}^e\approx \sum _{A=1}^{n_{\mathrm{nen}}}\delta {\mathsf{\Theta }}^{eA}{N}^A,\hfill \end{array}$$

(84)

$$\begin{array}{cc}\hfill & {\nabla }_{\mathit{X}}\delta {\mathit{u}}^e\approx \sum _{A=1}^{n_{\mathrm{en}}}\delta {\mathit{u}}^{eA}\otimes {\nabla }_{\mathit{X}}{N}^A,{\nabla }_{\mathit{X}}\delta {\mathsf{\Theta }}^e\approx \sum _{A=1}^{n_{\mathrm{nen}}}\delta {\mathsf{\Theta }}^{eA}{\nabla }_{\mathit{X}}{N}^A,\hfill \end{array}$$

(85)

where $\delta {\mathit{u}}^{eA}$ and $\delta {\mathsf{\Theta }}^e$ are the values at node A of element e. The interpolation is done using C${}^0$-continuous Lagrange basis functions $N^A$. The application of Equations (84) and (85) on the discretized weak forms (Equations (82) and (83)) yields

$$\begin{array}{cc}\hfill & \sum _{e=1}^{n_{\mathrm{el}}}\left\{\sum _{A=1}^{n_{\mathrm{en}}}\delta {\mathit{u}}^{eA}·\left[{\mathit{f}}_{u,\mathrm{int}}^{eA}-{\mathit{f}}_{u,\mathrm{vol}}^{eA}-{\mathit{f}}_{u,\mathrm{sur}}^{eA}\right]\right\}=0,\hfill \end{array}$$

(86)

$$\begin{array}{cc}\hfill & \sum _{e=1}^{n_{\mathrm{el}}}\left\{\sum _{A=1}^{n_{\mathrm{nen}}}\delta {\mathsf{\Theta }}^{eA}\left[f_{\mathsf{\Theta },\mathrm{trans}}^{eA}+f_{\mathsf{\Theta },\mathrm{cond}}^{eA}-f_{\mathsf{\Theta },\mathrm{vol}}^{eA}-f_{\mathsf{\Theta },\mathrm{sur}}^{eA}-f_{\mathsf{\Theta },\mathrm{crys}}^{eA}-f_{\mathsf{\Theta },\mathrm{flex}}^{eA}\right]\right\}=0.\hfill \end{array}$$

(87)

The force contributions in the discretized weak form of the balance of linear momentum (Equation (86)) correspond to the internal forces ${\mathit{f}}_{u,\mathrm{int}}^{eA}$, the volumetric forces ${\mathit{f}}_{u,\mathrm{vol}}^{eA}$, and the surface forces ${\mathit{f}}_{u,\mathrm{sur}}^{eA}$:

$$\begin{array}{cc}\hfill & {\mathit{f}}_{u,\mathrm{int}}^{eA}:={\int }_{{\mathcal{B}}^e}\mathit{P}·{\nabla }_{\mathit{X}}{N}^A\mathrm{d}V,{\mathit{f}}_{u,\mathrm{vol}}^{eA}:={\int }_{{\mathcal{B}}^e}{\rho }_0{N}^A\mathit{b}\mathrm{d}V,\hfill \end{array}$$

(88)

$$\begin{array}{cc}\hfill & {\mathit{f}}_{u,\mathrm{sur}}^{eA}:={\int }_{\partial {\mathcal{B}}^{e\mathit{t}}}{N}^A\mathit{t}\mathrm{d}A,\hfill \end{array}$$

(89)

whereas, the discretized weak form of the balance of energy (Equation (87)) includes fluxes $f_{\mathsf{\Theta },\mathrm{tran}}^{eA}$ describing transient heat transfer, fluxes $f_{\mathsf{\Theta },\mathrm{cond}}^{eA}$ including heat conduction, volumetric fluxes $f_{\mathsf{\Theta },\mathrm{vol}}^{eA}$, surface fluxes $f_{\mathsf{\Theta },\mathrm{sur}}^{eA}$, fluxes $f_{\mathsf{\Theta },\mathrm{crys}}^{eA}$ yielding a temperature change caused by the crystallization, and fluxes $f_{\mathsf{\Theta },\mathrm{flex}}^{eA}$ changing the temperature due to thermal flexibility

$$\begin{array}{cc}\hfill & f_{\mathsf{\Theta },\mathrm{tran}}^{eA}:=\sum _{B=1}^{n_{\mathrm{nen}}}{\dot{\mathsf{\Theta }}}^{eB}{\int }_{{\mathcal{B}}^e}{N}^A{c}_d{N}^B\mathrm{d}V,f_{\mathsf{\Theta },\mathrm{vol}}^{eA}:={\int }_{{\mathcal{B}}^e}{N}^A{r}_{\mathsf{\Theta }}\mathrm{d}V,\hfill \end{array}$$

(90)

$$\begin{array}{cc}\hfill & f_{\mathsf{\Theta },\mathrm{cond}}^{eA}:=-{\int }_{{\mathcal{B}}^e}{\nabla }_{\mathit{X}}{N}^A·\frac{1}{{\rho }_0}{\mathit{q}}_0\mathrm{d}V,f_{\mathsf{\Theta },\mathrm{sur}}^{eA}:=-{\int }_{\partial {\mathcal{B}}^e}{N}^A\frac{1}{{\rho }_0}{\mathit{q}}_0·\mathit{n}\mathrm{d}A,\hfill \end{array}$$

(91)

$$\begin{array}{cc}\hfill & f_{\mathsf{\Theta },\mathrm{crys}}^{eA}:={\int }_{{\mathcal{B}}^e}{N}^A{\tilde{q}}_{\chi }\dot{\chi }\mathrm{d}V,f_{\mathsf{\Theta },\mathrm{flex}}^{eA}:={\int }_{{\mathcal{B}}^e}{N}^A\left({\tilde{q}}_{\alpha }+\mathsf{\Theta }{c}_2\right)\dot{\alpha }\mathrm{d}V.\hfill \end{array}$$

(92)

A special focus lies on the last flux contribution (Equation (92)), which after inserting Equation (59) yields

$$f_{\mathsf{\Theta },\mathrm{flex}}^{eA}:={\int }_{{\mathcal{B}}^e}{N}^A\left(q_{\alpha }\dot{\alpha }+c_2{\mathsf{\Theta }}_0\dot{\alpha }\right)\mathrm{d}V.$$

(93)

The first term in the parenthesis in Equation (93) is non-negative due to dissipation inequality (55). However, the second term represents a product of positive constants $c_2{\mathsf{\Theta }}_0$ with rate $\dot{\alpha }$ which can be negative during the unloading phase. This allows the reduction of temperature in the unloading phase as it is observed in experimental results (Figure 1b). Finally, the element contributions are assembled to a global system of equations under the consideration of kinematic compatibility

$$\begin{array}{cc}\hfill & {\left[\begin{array}{c}\delta \mathit{u}\\ \delta \mathsf{\Theta }\end{array}\right]}^T·\mathit{R}\left(\mathit{u},\mathsf{\Theta }\right)=0,\hfill \end{array}$$

(94)

$$\begin{array}{cc}\hfill & \mathit{R}\left(\mathit{u},\mathsf{\Theta }\right)=\left[\begin{array}{c}{\mathit{f}}_{u,\mathrm{int}}-{\mathit{f}}_{u,\mathrm{vol}}-{\mathit{f}}_{u,\mathrm{sur}}\\ {\mathit{f}}_{\mathsf{\Theta },\mathrm{tran}}+{\mathit{f}}_{\mathsf{\Theta },\mathrm{cond}}-{\mathit{f}}_{\mathsf{\Theta },\mathrm{vol}}-{\mathit{f}}_{\mathsf{\Theta },\mathrm{sur}}-{\mathit{f}}_{\mathsf{\Theta },\mathrm{crys}}-{\mathit{f}}_{\mathsf{\Theta },\mathrm{flex}}\end{array}\right].\hfill \end{array}$$

(95)

Provided that $\mathit{R}=\mathbf{0}$ holds for all admissible $\left[\delta \mathit{u}\delta \mathsf{\Theta }\right]$, problem (94) can be solved by any solver for a nonlinear system of equations. Most commonly, the Newton-Raphson method is applied for this purpose. This procedure relies on the following iteration rule

$$\Delta {\mathit{g}}_i=-{\mathit{K}}^{-1}·\mathit{R},{\mathit{g}}_{i+1}^n={\mathit{g}}_i^n+\Delta {\mathit{g}}_i.$$

(96)

Here, i denotes the Newton iteration counter, n is the current time step and $\Delta {\mathit{g}}_i=[\Delta {\mathit{u}}_i\Delta {\mathsf{\Theta }}_i]$ are the increments of the global variables ${\mathit{g}}_i^n=\left[{\mathit{u}}_i^n{\mathsf{\Theta }}_i^n\right]$ related to the time step n. Finally, $\mathit{K}=\frac{\partial \mathit{R}}{\partial \mathit{g}}$ represents the global stiffness matrix of the system.

4.2. Time Discretization and Simulation of the Unloading Phase

Experimental results show that the loading stage is related to the regularity and temperature increase, whereas the degradation of crystalline regions accompanied by the decrease of the temperature occurs during the unloading stage. In the present model, evolution Equations (72) and (74) control the process described. Their numerical implementation relies on the time discretization of the change of crystallization limit (Equation (64)), of the crystallization parameter (Equation (73)) and of both driving forces (Equations (58) and (59)). The explicit scheme approximating time derivatives by a forward difference quotient is chosen for this purpose:

$$\begin{array}{cc}\hfill & {\chi }_{n+1}={\chi }_n+\Delta \lambda {\tilde{q}}_{\chi n},\hfill \end{array}$$

(97)

$$\begin{array}{cc}\hfill & {\tilde{q}}_{\chi n}=k_1\parallel {\mathit{M}}_n^{\mathrm{c},\mathrm{dev}}\parallel -c_1,\hfill \end{array}$$

(98)

$$\begin{array}{cc}\hfill & \Delta \lambda =\frac{\Delta tf\left({\chi }_n\right){\tilde{q}}_{\chi n}{k}_1\left(\parallel {\mathit{M}}_n^{\mathrm{c},\mathrm{dev}}\parallel -\parallel {\mathit{M}}_{n-1}^{\mathrm{c},\mathrm{dev}}\parallel \right)}{b{\left(A+B_n\right)}^2},\hfill \end{array}$$

(99)

$$\begin{array}{cc}\hfill & B_{n+1}=B_n+\frac{b}{f\left({\chi }_n\right)}\left|{\chi }_{n+1}-{\chi }_n\right|,\hfill \end{array}$$

(100)

$$\begin{array}{cc}\hfill & {\alpha }_{n+1}={\alpha }_n+\Delta t\frac{{\chi }_n^{D_2}}{D_1}{\tilde{q}}_{\alpha n},\hfill \end{array}$$

(101)

$$\begin{array}{cc}\hfill & {\tilde{q}}_{\alpha n}=k_2\mathrm{tr}\left({\mathit{M}}_n^{\mathrm{th}}\right)-c_2\left({\mathsf{\Theta }}_n-{\mathsf{\Theta }}_0\right).\hfill \end{array}$$

(102)

Here, subscript $n+1$ denotes values at the current time step and subscript n refers to the values in the previous time step. It is important to point out that the evaluation of the crystalline and thermal Mandel stresses in the above expressions involve the time integration of tensor valued quantities. For this reason, evolution laws (53) and (54) are numerically solved by applying the exponential map [65] along with definitions (27) and (29)

$$\begin{array}{cc}\hfill & {\mathit{F}}_{n+1}^{\mathrm{c}}=\exp \left(k_1\left({\chi }_{n+1}-{\chi }_n\right)\frac{{\mathit{M}}_n^{\mathrm{c},\mathrm{dev}}}{\parallel {\mathit{M}}_n^{\mathrm{c},\mathrm{dev}}\parallel }\right)·{\mathit{F}}_n^{\mathrm{c}},\hfill \end{array}$$

(103)

$$\begin{array}{cc}\hfill & {\mathit{F}}_{n+1}^{\mathrm{th}}=\exp \left(k_2\left({\alpha }_{n+1}-{\alpha }_n\right)\mathit{I}\right)·{\mathit{F}}_n^{\mathrm{th}}.\hfill \end{array}$$

(104)

The method used in the current approach goes back to the definition of the tensor exponential, where the numerical solution is carried out by calculating a finite truncation of the Taylor series. Equation (104) is an exact solution due to the diagonality of identity $\mathit{I}$. The contribution by Moler and Van Loan (2003) [66] discusses various approaches to compute the exponential of a second order tensor and compares their applicability and efficiency.

Within the framework of the numerical implementation, the simulation of the degradation of crystalline regions deserves special attention. According to Equation (61), the evolution direction of the regularity is controlled by the sign of the driving force due to the non-negativity of $\lambda$. In other words, a negative driving force during the unloading stage leads to a reduction of the regularity which is achieved by choosing a suitable shift $c_1$ in Equation (98). This load-dependent parameter is calculated from the condition for the initial value of driving force ${\tilde{q}}_{\chi }^{\mathrm{un},\mathrm{in}}$ to coincide with the negative crystallization limit if increment B is set to zero

$${\tilde{q}}_{\chi }^{\mathrm{un},\mathrm{in}}=k_1\parallel {\mathit{M}}_{\mathrm{end},\mathrm{ld}}^{\mathrm{c},\mathrm{dev}}\parallel -c_1=-A\Rightarrow c_1=A+k_1\parallel {\mathit{M}}_{\mathrm{end},\mathrm{ld}}^{\mathrm{c},\mathrm{dev}}\parallel .$$

(105)

Here, ${\mathit{M}}_{\mathrm{end},\mathrm{ld}}^{\mathrm{c},\mathrm{dev}}$ is the deviatoric crystalline Mandel stress tensor at the end of the loading stage, subscript “ld” denotes the loading stage, superscript “un” denotes the unloading stage, and “in” is an initial value.

4.3. Algorithmic Aspects of the SIC Model Implementation

The SIC model is implemented in the FEAP-software that is appropriate for an enhanced FE-analysis of the complex material behavior. A new element along with a new material subroutine are developed to this end. Their structure and interconnection with the main program are presented in Figure 3.

The flow chart of the main program (Figure 3a) shows that the element subroutine is used at the beginning of the simulations to allocate material properties to each single element and to generate data on the initial network regularity. Thereafter, the element subroutine is called to calculate the element residual and stiffness matrix which are finally assembled in the global residual $\mathit{R}$ and global stiffness matrix $\mathit{K}$. The latter are used to solve the nonlinear system of Equations (94) and (95) according to the scheme (96). This iterative procedure is performed until the prescribed accuracy is achieved. To this end, a tolerance for the energy norm of displacement increments is set to the value $1\times {10}^{-12}$. The applied monolithic scheme solves the mechanical and thermal balance equations simultaneously with respect to the field variables, displacements and temperature. It provides final values of global external variables ${\mathit{u}}^n$ and ${\mathsf{\Theta }}^n$ for each time step n. The time iteration is performed on the basis of the Newmark method with the default parameters $\beta =0.25$ and $\gamma =0.5$. The incorporation of temporal aspects is needed in order to simulate the external load in an incremental form and, correspondingly, to update the solution of the boundary value problem (39)–(44).

The element subroutine includes two essential tasks, namely the generation of the initial network regularity and the evaluation of the element residual and of the element stiffness matrix (Figure 3b). The generation of the initial regularity is physically motivated by the fact that the natural amorphous entangled microstructure of polymer chains includes places with a different cross-linking degree and polymer chain arrangement. Some of these places are especially suitable for the formation of crystalline regions and behave as nuclei of the crystallization. In order to model such an initial configuration, the network regularity is first set to zero in all elements. In a second step, a uniform distribution is applied to generate initial values of the network regularity in the range $[1\times {10}^{-8},1\times {10}^{-3}]$ as well as to generate the element numbers associated to the particular initial values. Two intrinsic subroutines of the Fortran program language, random_seed and random_number, are used for this purpose. The volume fraction of elements with a higher initial value is varied and fitted with respect to the experimental results shown in Figure 1a. Depending on the initial configuration, the material model determines where the crystalline regions develop and how fast they grow.

The second core assignment of the element subroutine is the definition of the evaluation of the residual and of the element stiffness matrix. The basis for this step is already provided by definitions of forces and fluxes (Equations (88)–(92)). Accordingly, the residual vector corresponding to the element e has the form

$${\mathit{R}}^e=\left[\begin{array}{c}{\mathit{R}}_u^e\\ {\mathit{R}}_{\mathsf{\Theta }}^e\end{array}\right]=\left[\begin{array}{c}{\mathit{f}}_{u,\mathrm{int}}^e-{\mathit{f}}_{u,\mathrm{vol}}^e-{\mathit{f}}_{u,\mathrm{sur}}^e\\ {\mathit{f}}_{\mathsf{\Theta },\mathrm{tran}}^e+{\mathit{f}}_{\mathsf{\Theta },\mathrm{cond}}^e-{\mathit{f}}_{\mathsf{\Theta },\mathrm{vol}}^e-{\mathit{f}}_{\mathsf{\Theta },\mathrm{sur}}^e-{\mathit{f}}_{\mathsf{\Theta },\mathrm{crys}}^e-{\mathit{f}}_{\mathsf{\Theta },\mathrm{flex}}^e\end{array}\right],$$

(106)

and corresponding stiffness matrix ${\mathit{K}}^e$ results from the derivative of the residual (106) with respect to the global field variables $\mathit{g}={\left[\mathit{u}\mathsf{\Theta }\right]}^T$

$${\mathit{K}}^e=\frac{\partial {\mathit{R}}^e}{\partial \mathit{g}}=\left[\begin{array}{cc}{\displaystyle \frac{\partial {\mathit{R}}_u^e}{\partial \mathit{u}}}& {\displaystyle \frac{\partial {\mathit{R}}_u^e}{\partial \mathsf{\Theta }}}\\ {\displaystyle \frac{\partial {\mathit{R}}_{\mathsf{\Theta }}^e}{\partial \mathit{u}}}& {\displaystyle \frac{\partial {\mathit{R}}_{\mathsf{\Theta }}^e}{\partial \mathsf{\Theta }}}\end{array}\right]=\left[\begin{array}{cc}{\displaystyle \frac{\partial {\mathit{f}}_{u,\mathrm{int}}^e}{\partial \mathit{u}}}& \mathbf{0}\\ -{\displaystyle \frac{\partial {\mathit{f}}_{\mathsf{\Theta },\mathrm{crys}}^e}{\partial \mathit{u}}}-{\displaystyle \frac{\partial {\mathit{f}}_{\mathsf{\Theta },\mathrm{flex}}^e}{\partial \mathit{u}}}& {\displaystyle \frac{\partial {\mathit{f}}_{\mathsf{\Theta },\mathrm{tran}}^e}{\partial \mathsf{\Theta }}}+{\displaystyle \frac{\partial {\mathit{f}}_{\mathsf{\Theta },\mathrm{cond}}^e}{\partial \mathsf{\Theta }}}-{\displaystyle \frac{\partial {\mathit{f}}_{\mathsf{\Theta },\mathrm{flex}}^e}{\partial \mathsf{\Theta }}}\end{array}\right].$$

(107)

Amongst others, previous equations depend on the stress tensor and on the latent heats which are calculated in the material model subroutine (Figure 3c). Here, the input variables are temperature and deformation gradient. Criterion (61) distinguishes the loading and the unloading mode. Such a differentiation is needed since the shift $c_1$ (Equation (105)) and the appropriate value for constant b (Equation (64)) are defined at the transition from the loading to the unloading mode. Subsequently, three following tasks are accomplished:

• the update of internal variables according to the Equations (97)–(102), which also requires the calculation of Mandel stress tensors (Equations (30) and (31)),
• the calculations of the crystalline and the thermal deformation gradient are updated by Equations (103) and (104),
• the calculations of the first Piola–Kirchhoff stress tensor (Equation (25)), and latent heats (Equation (78)). These quantities are final results passed to the element subroutine.

5. Numerical Examples

5.1. Simulation of Single Crystalline Regions Embedded in the Matrix Material

The first numerical example deals with the simulation of a tensile test performed on a two-dimensional sample that demonstrates the material behavior under the influence of the SIC. The elastic material constants are chosen for rubber and correspond to the original amorphous phase [67,68]. The crystalline material constants are fitted to the experimental results by Toki et al. (2003) [69] and Candau et al. (2015) [12]. The latter work is also used for fitting thermal material constants (Figure 1b). An overview of all material parameters is given in

Table 1. Elastic, crystalline, and thermal material constants used in simulations.

Elastic Constants
Bulk modulus	K	$5\times {10}^8$	Pa
Shear modulus	$\mu$	$4\times {10}^5$	Pa
Limiting network stretch	${\lambda }_m$	2	-
Crystalline Constants
Coupling constant	$k_1$	$7\times {10}^{-2}$	-
Crystallization limit	A	$1\times {10}^5$	Pa
Hardening constant	${\mathit{b}}_1$	$1.7\times {10}^5$	Pa
Softening constant	${\mathit{b}}_2$	$2\times {10}^5$	Pa
Constant in function f	${\beta }_1$	$0.25$	-
Constant in function f	${\beta }_2$	$0.5$	-
Exponent in function f	${\beta }_3$	2	-
Thermal Constants
Reference temperature	${\mathsf{\Theta }}_0$	300	K
Thermal flexibility modulus	$c_2$	$2\times {10}^5$	Pa/K
Coupling constant	$k_2$	$4\times {10}^{-2}$	-
Thermal flexibility coefficient	$D_1$	$3.3\times {10}^6$	Pa
Exponent	$D_2$	6	-

. The thermal conductivity is neglected (${\lambda }_{\mathsf{\Theta }}=0$), as the focus lies on the heat generated or absorbed by SIC and not on the flow of heat within the material. The conductivity of natural rubber is low compared to the other materials and amounts to 0.15 W/(m K) [70]. Moreover, no external heat sources ($r_{\mathsf{\Theta }}=0$) are applied in the simulations.

The setup corresponding to the tensile test is shown in Figure 4a. The chosen square specimen has the dimensions $100\times 100$ nm and is discretized by $50\times 50$ square elements. For this purpose, four-node quadrilateral elements with bilinear shape functions are assumed. They are fully integrated by using four integration points. In addition, the simulations are performed by applying quadrilateral elements with the reduced one-point integration. However, both types of simulations yield comparable results and no effects of volume locking or of hourglassing are observed. The sample thickness is 1 nm and is significantly smaller than the remaining dimensions, which corresponds to a plane stress state problem. The application of the model to 3D simulations would be straightforward, since the general SIC material model is presented in previous sections. Vertical displacements prescribed at the horizontal boundaries increase linearly up to the maximum value of 250 nm and thereafter decrease linearly to zero (Figure 4b). The displacement increment is set to $\left|\Delta \overline{u}\right|=2.5\times {10}^{-2}$ nm. Here, the bar symbol indicates external influences. The prescribed stretch is calculated according to $\overline{\lambda }=(l+2\overline{u})/l$. The total loading time amounts to 10 s, while the time increment is $\Delta t=1\times {10}^{-3}$ s.

The example has an academic character and investigates the material response under the presence of three separate crystalline regions embedded in the amorphous matrix material (Figure 4). Furthermore, the influence of the initial value of the regularity degree is studied. In the initial configuration (Figure 5a) the regularity is set to a value of ${\chi }_0=1\times {10}^{-3}$ at an element in the bottom part of the sample, to ${\chi }_0=1\times {10}^{-4}$ at an element in the middle of the sample and to ${\chi }_0=1\times {10}^{-5}$ at an element in the upper part of the sample. The evolution of the regularity over the time at the crystalline region in the bottom is shown for the complete load cycle in Figure 5b. As expected, the regularity starts to evolve after exceeding the crystallization limit (Equation (71)) and reaches the value of $\chi \approx 1$ corresponding to full crystallization. During the unloading phase, it decreases with a slower rate (Equation (64)) and ends up at the value of $\chi \approx 0$ complying with the amorphous material. A wide range of experimental results show that crystalline regions completely recede at stretches between 250% and 350% and that the original amorphous phase is completely recovered. Examples giving evidence of this phenomenon are presented in Le Gac et al. (2018) [71] for polychloroprene and in Huneau (2011) [72] for polyisoprene. One exception is shown in the work by Albouy et al. (2005) dealing with the inverse yielding effect. An incomplete recovery only occurs if the sample was previously loaded beyond a certain critical draw ratio. Additional factors that benefit this effect are a low cross-linking density and lower temperatures. The stress–stretch curve (Figure 5c) builds a hysteresis loop which is in agreement with the experimental results (Figure 1a). Here, the ${\mathit{P}}_{22}$ component of the first Piola–Kirchhoff stress tensor is plotted. The temperature curve (Figure 5d) shows similar tendencies as the regularity: it increases during loading and reduces during unloading. The remaining temperature at the end of the unloading phase is due to the produced dissipation (Equation (75)) which contributes to Equation (79) and leads to temperature increase.

In addition, Figure 6a monitors the microstructure development during the loading phase at a stretch of $\overline{\lambda }=4.9$. The resulting temperature distribution is shown in Figure 6b. The crystalline regions build up and grow faster at the areas close to the element with a higher initial value of regularity. The same holds for the temperature which rises more strongly in these areas. Figure 7 monitors the fully developed crystalline regions and the corresponding temperature distribution at the end of the loading phase ($\overline{\lambda }=6.0$). The highest temperature occurs directly in the crystallized region. Interestingly, the area surrounding the crystal has cooled down. Equilibrium is established in the deformed state, with the consequence that heat is absorbed from the surrounding areas.

5.2. Microstructure Evolution for a Sample with the Complex Initial Microconfiguration

A further example simulates the tensile test for a sample with randomly generated initial values of network regularity (Figure 8a), which is the situation to be expected in a real polymer. The initial values are generated in the range [0, $1\times {10}^{-3}$] and randomly distributed over the sample as explained in Section 4.3. The higher values of the regularity represent potential nuclei of crystalline regions. All remaining material and process parameters except the initial regularity are kept as in Section 5.1. Comparable simulations are performed using quadrilateral elements with the reduced integration; however, no evidence of volume locking is noticed.

Three snapshots are chosen to display the microstructure development and the temperature distribution: Figure 9 shows the microstructure corresponding to an applied stretch of $\overline{\lambda }=4.8$ during the loading phase, Figure 10 presents the situation at the end of the loading phase $\overline{\lambda }=6.0$, and Figure 11 shows the microstructure for a stretch of $\overline{\lambda }=3.7$ during unloading. The color scale in Figure 8a differs from the color scale in Figure 9, Figure 10 and Figure 11 in order to visualize the initial microstructure better. Crystalline regions (Figure 9a) form during loading and this process generates heat so that an increase in temperature can be observed (Figure 9b). At the highest stretch state, the degree of crystallinity reaches about 15% (Figure 10a). The highest temperatures are located in regions with a high density of crystals, since this is where the heat development of adjacent crystals accumulates. These temperatures correspond to the red color in Figure 10b. In addition, it can be seen that the immediate surroundings of the crystals cool down (purple areas). As already shown in the example with three crystals, a withdrawal of heat from the amorphous regions also takes place during crystallization. Amorphous regions located at a larger distance from the crystals hardly experience any temperature changes and are shown in blue. During the unloading phase the regression of the crystallization (Figure 11a) absorbs heat, thereby reducing the temperature (Figure 11b). The comparison of the two states from Figure 9 and Figure 11 also demonstrates that the rate of the network regularity is higher during the loading stage than during the unloading stage at the same level of the stretch $\overline{\lambda }$.

The simulations shown also enable an analysis of the effective material behavior, which is illustrated by examples of the crystallinity degree and of the effective temperature. The crystallinity degree in the proposed model is defined as the volume fraction of the crystalline regions with the network regularity over 80%. Its evolution compared to the experimental results is presented in Figure 12a. The crystalline regions start to build at $\overline{\lambda }=4.3$ and their volume fraction gradually increases up to the value of 15% at the end of the loading phase. The crystallinity degree decreases during the unloading stage and crystalline regions completely vanish at $\overline{\lambda }=3.0$. The rate of change during the loading phase is higher than during the unloading stage. A comparison shows that model for a chosen set of parameters predicts a slightly faster formation and degradation of the crystals than is observed in experiments. It also should be noted that data on the start of the crystallization, as well as the value of the crystallinity degree at the peak are used to fit the parameters in Table 1 and to fit the volume fraction of the sample with the initial regularity different than zero.

The effective temperature ${\mathsf{\Theta }}^{\mathrm{ef}}$ is evaluated according to the principle of the volume averaging

$${\mathsf{\Theta }}^{\mathrm{ef}}=\frac{1}{V}{\int }_{\mathcal{B}}\mathsf{\Theta }\mathrm{d}V$$

(108)

and corresponding change with respect to the initial temperature is shown in Figure 12b. The temperature strongly increases after reaching $\overline{\lambda }=4.3$ since the crystallization process starts. At the end of the loading stage, the model predicts a temperature which excellently agrees with the experimental value and corresponds to the temperature change of 6 K. During the unloading phase, the temperature gradually decreases with a moderate slope; however, it does not completely recover its initial value. This discrepancy can be explained by the fact that the present model simulates a closed system, whereas an exchange of the heat with the surroundings is possible in experiments [12].

6. Conclusions and Outlook

This contribution presents a thermodynamically consistent model for the SIC in unfilled polymers under consideration of mechanical and thermal effects. The kinematics is defined by distinguishing configurations due to elastic, crystalline, and thermal deformations, which yields the multiplicative split of the deformation gradient into three contributions. The evolution equations of the internal variables describing the regularity of polymer network and the flexibility of polymer chains as a result of heat generation by crystallization are derived from the dissipation potential that has been established specifically for this type of material. The assumption for Helmholtz energy includes the Arruda–Boyce elastic energy term, the crystalline energy term, the purely thermal energy part, and the mixed energy contribution dependent on temperature and thermal flexibility. Finally, the weak form of the BVP including the balance of the linear momentum and balance of energy along with the suitable boundary conditions is defined for the purpose of the numerical implementation of the model proposed into the FEM-software.

The application of the model is illustrated by examples dealing with the microstructure evolution during a cyclic tensile test. In the first example, the sample contains three nuclei embedded in the amorphous matrix material. The aim of this academic test is to investigate the material response affected by the SIC. The numerical results show the evolution of the network regularity, the stress–stretch curve building a hysteresis, and the temperature change with similar trends as the regularity. In addition, the snapshot of the temperature distribution at the end of the loading path shows that heat is generated in the crystalline region, whereas the adjacent amorphous area cools down. A further example simulates the behavior of a sample with the configuration to be expected in a real amorphous polymer and enables the visualization of the SIC process by means of snapshots in several load stages. In the temperature diagrams, the heat generation and absorption due to the formation or degradation of the crystallization can be observed. It is noticeable that the temperature especially increases at places with clusters of crystalline regions, while the amorphous matrix outside the “active” areas hardly experiences any temperature changes at all. The obtained results also serve as a basis for the evaluation of the crystallinity degree and of the effective temperature which plays an important role in the experimental validation.

Apart from the issues mentioned, the developed model also gives rise to other investigations. In a first step, the assumptions for the Helmholtz energy and the dissipation potential can be extended in order to simulate further effects occurring in filled and unfilled rubbers. Some important topics in this context would be the Mullins effect, a deformation state beyond the elastic limit, induced anisotropy, the Gough–Joule effect, and heat conduction [73,74]. In addition, the model proposed can be coupled to the phase–field approach in order to represent the biphasic nature of material better. Internal variable $\chi$ would correspond to the order parameter in that case, and its evolution could be controlled by the same dissipation potential as proposed in the present work.