Modeling and Simulation of the Aging Behavior of a Zinc Die Casting Alloy

Abstract

While zinc die-casting alloy Zamak is widely used in vehicles and machines, its solidified state has yet to be thoroughly investigated experimentally or mathematically modeled. The material behavior is characterized by temperature and rate sensitivity, aging, and long-term influences under external loads. Thus, we model the thermo-mechanical behavior of Zamak in the solid state for a temperature range from −40 °C to 85 °C, and the aging state up to one year. The finite strain thermo-viscoplasticity model is derived from an extensive experimental campaign. This campaign involved tension, compression, and torsion tests at various temperatures and aging states. Furthermore, the thermo-physical properties of temperature- and aging-dependent heat capacity and heat conductivity are considered. One significant challenge is related to the multiplicative decompositions of the deformation gradient, which affects strain and stress measures relative to different intermediate configurations. The entire model is implemented into an implicit finite element program and validation examples at more complex parts are provided so that the predicability for complex parts is available, which has not been possible so far. Validation experiments using digital image correlation confirm the accuracy of the thermo-mechanically consistent constitutive equations for complex geometrical shapes. Moroever, validation measures are introduced and applied for a complex geometrical shape of a zinc die casting specimen. This provides a measure of the deformation state for complex components under real operating conditions.

1. Introduction

Zinc die casting alloys are commonly used in the production of components of cars and machines, as well as in the electronic, construction, and appliance industries, see [1,2]. They show an excellent castability along with low production costs. Their suitability for the hot chamber die casting process comes from their low melting point. However, the material implies three main disadvantages: a strong viscoplastic behavior in the temperature range of applications (for the automotive industry this lies approximately between −40 °C and 85 °C), a temperature-dependent mechanical response, see [3] and a change in the thermo-mechanical behavior due to aging of the material, see [4]. This has to be considered in the design process requiring a suitable model.

In comparison to other materials such as steel or aluminum alloys, zinc die casting alloys have not been sufficiently investigated so far. One can find several investigations from the point of view of Material Science, in which the microstructure is characterized for different zinc die casting alloys under different thermal conditions. The microstructural changes of zinc-aluminum alloys became subject of investigations in the last century. Refs. [5,6,7] characterized the micro-structure, phase transformations, and volume changes of the zinc-aluminum and zinc-aluminum-copper alloys over time. Moreover, the reduction of dimensions (shrinkage) in zinc die casting alloys over the course of time is described in [8]. Murray [9] investigated the equilibrium phase diagram of zinc-aluminum alloys as well as their different crystal structures. In [10] the analyses of the influence of the alloying elements on the mechanical properties, studies focusing on other aspects such as the creep behavior, aging, the consequent volume changes, as well as inter-crystalline corrosion can be found. Although in the nineties, the interest in zinc die casting alloys reduced, new investigations in this field started to appear in the past years. At the beginning of the 21st century, phase transformations and precipitation of several Zn-Al alloys were extensively characterized in [11,12,13,14,15,16]. Leis and Kallien [2] show the characterization of the influence of temperature, natural, and artificial aging, the thickness-dependence of the specimen, and die casting conditions on several mechanical parameters such as the Young’s modulus, yield strength, or tensile strength. Further characterizations of the effect of alloying components on the microstructure, phase decomposition, or mechanical characterizations in zinc die casting alloys can be found in the works of Jareno et al. [17], Liu et al. [18], Savaskan and Pasa Hekimoglu [19], Pola et al. [20], and Wu et al. [21].

There are very few proposals for modeling the mechanical behavior of zinc die casting alloys in the literature. Some of them consider the die casting process itself, see [22,23], or investigations from the point of view of Material Sciences, see [2,13,14,15,16,17,18,19,20,21,22,23,24]. A wide amount of models have been proposed for aluminum alloys and steel. One can find several proposals for other alloys, such as titanium, magnesium, nickel or copper alloys, see [25,26,27,28,29]. Most of these alloys are exposed to some kind of mechanical treatment before they are modeled or their behavior differs strongly from zinc alloys. Finally, from the author’s side, a model of thermo-viscoplasticity for small deformations of Zamak 5 is proposed in [3], which was extended to natural aging at room temperature in [4]. Moreover, in [30], the thickness dependence was modeled considering the porosity for a purely elastic material behavior.

One major issue of this work is the modeling of aging effects in zinc die casting alloys and its extension to large strains. This provides a new prediction tool that—for the first time—makes it possible to predict the thermomechanical behavior of a zinc die cast material at different aging states.

In the field of theory of materials in a continuum mechanical framework, Krempl [31] defines aging as a modification of the micro-structure, in such a way that the same loading process, performed in the same environment at different times (or ages of the material), will lead to different mechanical responses. The origin of these modifications lies in the microstructure of the material. The cause of aging is different for polymers and alloys. Chaboche [32] distinguishes three different types of aging in alloys: Dynamic strain aging is caused by “dragging” of the dislocations of the atoms in solution and is associated with the “Portevin-Le Chatelier effect”. Static strain aging implies a growth in material strength with time, which is observed as a reverse effect of static recovery. In the literature, there are several suggestions on how to model these two kinds of aging types in aluminum alloys and steel, but in the present study of Zamak, none of these types of aging was observed. Finally, aging due to phase changes is a consequence of dissolutions, precipitations, and phase transformations in alloys, which fits with the observations made in [33]. In the case of polymers, aging is related to physical and chemical changes, which limits their application period, see [34]. From a macroscopic point of view, the effect of aging can be introduced into the model with an internal variable, see [35,36] and in the area of polymers see [34,37,38,39,40]. In the literature one can find several proposals for modeling aging. Marquis and Lemaitre [35] proposed a model coupling elasto-plasticity, damage and aging for small deformations. Additionally, a model for plasticity and aging is suggested in Marquis and Costa Mattos [41]. In these models, the effect of aging is included with a new internal variable. In [42], chemical aging is included in a model for adhesive bonds in a polymer using an additional field variable to describe the diffusive processes. Moreover, in [38], a model for small deformations with aging as an internal variable is developed. In [34,37,39], a chemo-thermomechanical and a chemical aging model for elastomers for small deformations is suggested, and Dippel et al. [40] offers a solution for polymer bonds at large deformations. A further area of interest of aging effects can be found on non-metallic construction materials such as concrete, see [43,44,45], which are less appropriate for metals. In general, volume changes due to aging are not considered in the mentioned models.

The thermo-mechanical coupled problem has been treated in a large amount of works. Thermodynamical consistent modeling proposals for thermo-viscoplastic material behavior can be found in, for example [36,46,47,48]. Moreover, the work of Chaboche [32] offers an extensive overview of small-strain viscoplasticity models considering the thermo-mechanical coupling. In the case of finite deformations, the introduction of different effects in the model can be performed through a multiplicative decomposition of the deformation gradient. Lion [49] proposes constitutive equations for finite thermo-viscoplasticity deduced from rheological models. In this work, an additional multiplicative split of the deformation gradient is proposed in order to consider each physical effect included in the model: aging, temperature and mechanical effects. The equivalence of the order of decomposition of the deformation gradient for thermoelasticity has been discussed in [50]. The works of [51,52,53] provide investigations in the field of thermo-elasticity and thermo-viscoelasticity. Moreover, in the works of [54,55] one can find an extensive literature review for the coupled thermo-mechanical problem for the case of thermo-viscoelasticity and viscoplasticity for polymers and steel respectively. In the literature one can find many different proposals for aluminum alloys and steel in the field of thermo-viscoplasticity of metals, although, to the author’s best knowledge, there are no proposals for the quasi-static behavior of Zamak apart from the previously made proposals from the authors.

The goal of this work is, on the one hand, to extend the available information about the zinc die casting alloy Zamak 5 providing temperature and aging-dependent experimental results for its thermo-mechanical characterization, and, on the other hand, to offer a constitutive model of finite deformations which is able to reproduce this complex strain rate, temperature and aging-dependent behavior. On the one hand, significant deformations occur in conjunction with substantial rigid body movements, necessitating the objectivity of stress and strain measurements. Conversely, the determination of material parameters does not mitigate the inaccuracies associated with moderate strains. In light of these considerations, it becomes evident that there is a challenge in developing a thermodynamically consistent material model that fulfills the Clausius–Duhem inequality. This resulting model does not only address the interaction between thermal and mechanical influences but also, and more specifically, considers the aging state of the material. Thus, a highly complex coupling is treated. The model is based on several multiplicative decompositions of the deformation gradient, yielding various pure strain measures connected to the individual physical properties. These are considered in the proposal of a specific strain-energy function, where the evolution equations of particular strain quantities are proposed. Moreover, subsets of material parameters can be determined separately by particular experiments. Finally, the entire model, which is implemented into a finite element program, is validated using totally different specimens with complex geometries, where we try to study a validation metric to obtain insight into the predicability of the model. The investigations provided in this work are based on [56].

The following notation is used in this work: geometrical vectors are symbolized by $\overrightarrow{a}$ and second-order tensors $\mathbf{A}$ by bold-faced Roman letters. Matrices at so-called global levels of finite elements are written as bold-faced italic letters $\mathsf{A}$, and matrices on the local (element) level are represented using bold-faced Roman letters $\mathsf{A}$.

2. Experimental Work

In previous investigations published in [3], the thermo-mechanical behavior of Zamak 5 without aging has been characterized. It was found that the material shows a strong dependence on the strain rate and the temperature. Both factors affect the initial slope as well as the height of the stress-strain curve. In [4], the effect of aging at room temperature on the mechanical response of Zamak has been investigated. There, a decrease of the stiffness with the increasing aging time could be observed. Additionally, a shrinkage takes place over the course of time due to microstructural changes in the alloy. Moreover, in [33], it is shown that the thermal conductivity is also influenced by the aging time and the temperature, which is connected to the previously mentioned changes in the microstructure.

In this work, the effect of natural aging on the thermo-mechanical response of Zamak is characterized. To this end, thin-walled cylindrical specimens made of Zamak 5 (ZnAl4Cu1) and having the dimensions $R_i=13 \mathrm{m}\mathrm{m}$ of the inner radius, $d=2 \mathrm{m}\mathrm{m}$ of the wall thickness, and $L=90 \mathrm{m}\mathrm{m}$ of the length (region between the clamping areas) are tested at different natural aging times, see [3,4] as well. Tension and Torsion experiments were performed with specimens with 3, 6 and 12 months of aging at 19 °C (natural aging). Since in [3] it has been observed that the tension-compression asymmetry can be neglected for a typical application of the alloy, only the tension experiments are considered for the additional aging times. In order to consider strain-rate dependence, the experiments are performed for all aging times at four different strain rates each (${\dot{\gamma }}_k=3.2\times {10}^{-(k+3)}{\mathrm{s}}^{-1}$, $k=1,2,3,4$ for torsion and ${\dot{\mathit{\epsilon }}}_k=3.3\times {10}^{-(k+3)}{\mathrm{s}}^{-1}$, $k=1,2,3,4$ for tension). The additional temperatures of $-40 °\mathrm{C}$, $60 °\mathrm{C}$, and $85 °\mathrm{C}$ are used to investigate the temperature dependence since this is the typical application range of the temperature in the automotive industry. The equilibrium stress state—defining that stress state obtained for an infinitely slow loading, unloading, and reloading process—is characterized at the different aging times and temperatures with the help of multi-step relaxation tests. The results are compared to the previous observations made in [3,4], where the results at room temperature and without aging are available. Additionally, shrinkage at three different temperatures is measured. These shrinkage measurements are important for the identification of an internal aging variable in a similar way as in [4], since they directly connect the microstructural changes over the course of time with the aging process.

2.1. Influence of Aging on the Thermo-Mechanical Behavior

The results of the torsion tests for the four temperatures and aging times are shown in Figure 1, and for the case of tension in Figure 2.

Each figure represents the results at one temperature. The measures in the diagrams are defined as follows: $\gamma =\left(\tan \phi R_m\right)/L\approx \left(\phi R_m\right)/L$ defines the shear angle of applied to a thin-walled cylindrical specimen, where $R_m$ and L are the middle radius and the length of the specimen. $\phi$ symbolizes the applied testing machine’s rotation angle. $\tau =M_T/\left(2d{A}_m\right)$ represents the shear stress concerned, with the measured torque $M_T$, the thickness d of the specimen’s wall, and the mean cross-sectional area of the circular region $A_m=\pi R_m^2$. Furthermore, $\mathit{\epsilon }=(L_e-L_0)/L_0$ defines the axial strain, where $L_0=20 \mathrm{m}\mathrm{m}$ and $L_e$ are the initial and current length of an applied extensometer. $\sigma =F/A$ defines the axial stress, with the measured axial force F and the specimens cross-section A. Details can be found in [3,4].

In all figures, different strain rates are represented by different colors and for one strain rate, the four aging times are represented by different line styles. For the chosen aging times, the main influence of aging occurs at room temperature; see Figure 1b and Figure 2b. A reduction of $15\%$ in the stresses over the aging time can be observed at all strain rates. Furthermore, the decrement between the initial state and 3 months of aging is larger than the decrement between 3 and 6 months of aging. This can be explained by the aging kinetics, since the material is in the unaged state far from a microstructural equilibrium, and, thus, the diffusion processes are more pronounced in the beginning of the aging process and decrease over time. One can observe a reduction in the initial slope with increased aging as well, but the effects of the temperature and the strain rate are more pronounced for the initial slope. For higher temperatures and lower strain rates, the influence of aging declines, which indicates that the viscous effects are affected by aging, see Figure 1c,d and Figure 2c,d. At lower strain rates, the aging dependence is reduced until it almost vanishes at $85 °\mathrm{C}$ at the slowest strain rate.

At $-40 °\mathrm{C}$, the rate dependence is very small compared to the other temperatures. The effects of aging and rate dependence are almost overlapping. In the case of tension, see Figure 2a, it was not possible to perform all experiments up to the strain $\mathit{\epsilon }=1\%$ due to the brittleness of the material at this temperature. In relation to the equilibrium stress state, Figure 3 shows multi-step relaxation tests in torsion, and the estimated points of the equilibrium hysteresis for different aging times at four temperatures. Here, a slight dependence of the equilibrium stress state on the aging at room temperature can be observed. Again, this dependence is more pronounced at lower temperature and decreases with increasing temperatures until it becomes hardly perceptible at $85 °\mathrm{C}$. At higher temperatures, the viscous effects are predominant and the influence of aging in the equilibrium stress state is barely present.

The previous observations have to be summarized: aging has the greatest influence on the material behavior at room temperature. Its influence is less pronounced at higher—in the sense of the temperatures applied here—temperatures. It has a stronger effect on the viscous behavior of the material than on the internal equilibrium stress state, see, for the terminology [48] and the literature cited therein. Moreover, the influence of natural aging is moderate compared to the effect of temperature or the strain rate. Nonetheless, the influence due to aging has to be considered for certain application areas. Additionally, it has to be pointed out that the 12 months of natural aging do not represent the final state of aging in the sense that the micro-structure reaches an equilibrium state. The natural aging investigations up to one year were chosen based on the results of the investigations in ([57] p. 62). However, in our studies it has been found that the final state of a fully aged material is approximately five years for applications at room temperature, see Section 2.2. This is also in agreement with the observations made in [8]. Therefore, a further reduction in the initial slope and the maximum stresses is to be expected for longer aging times than 12 months. Furthermore, aging influences other mechanical properties as well, for example those that are not treated in this work, for instance elongation at fracture, see [2].

2.2. Shrinkage

As exposed in [33] and the literature cited therein, the changes in the mechanical response with the course of time originate in microstructural changes such as diffusion, precipitation and decomposition of several phases of the alloy. One of the effects of aging is a reduction in the dimension or shrinkage of the material. In the next section, an internal variable will be used to model the effect of aging on the mechanical response. The measurements of the shrinkage are of particular interest for the definition of the aging variable, since they directly connect the aging process with the measured quantities. The results of the experiments in Figure 1 and Figure 2 show a change in the mechanical response of the alloy in dependence of the aging time, but it is not possible to obtain a unique expression for the development of an aging variable from them. The works [8,57] mention a connection of the shrinkage in Zamak to a change in the crystallization system from fcc (face-centered cubic) to hcp (hexagonal close-packed), as well as to changes in the lattice parameter of some phases of the alloy. The deformation caused by these phase transformations associated to aging is purely volumetric. The measurements for specimens at room temperature from [8] show a maximum value of the shrinkage of about $-0.1\%$. Measurements of the shrinkage are an appropriate option to identify an internal process variable for aging in a unique manner, since this volume change is only connected to the microstructural changes, which occur during the aging process. The process is temperature-activated due to the fact that the microstructural changes during the aging process are related to diffusion and precipitation of several alloying elements. At higher temperatures there is a faster development of the shrinkage. Since it takes more than one year to reach the completely aged stage at room temperature, the measurements of the shrinkage have been performed at higher temperatures, so that the aging takes place in a reasonable period of time.

In the measurements, the specimen is heated up to the temperature ${\Theta }_{\exp }$. The heating rate is chosen to be fast enough so that shrinkage caused by the aging can be neglected during the heating phase and it is equal to $10 \mathrm{K} {\min }^{-1}$. The temperature is then kept constant over the time. The specimen has an initial length $L_0$ and it shrinks over time until the length reaches a saturation value $L_{\infty }$, see Figure 4.

As a measure of the aging development, the current longitudinal deformation ${\mathit{\epsilon }}_{\exp }=(L-L_0)/L_0$ is used.

The measurements of the shrinkage were carried out with a Mettler Toledo (TMA 841^e; performed at the Institute of Mechanics of the Bundeswehr University Munich). Figure 5 displays the results of the measurements, which were carried out at a temperatures of $85 °\mathrm{C}$, $105 °\mathrm{C}$, and $125 °\mathrm{C}$. Figure 5a—measurements over the time in logarithmic scale—is used to compare the shrinkage measurements at $20 °\mathrm{C}$ with results in [8] for the same alloy. We can observe a faster deformation rate at higher temperatures. The measurement at $125 °\mathrm{C}$ reaches a saturation value after approximately $18 \mathrm{h}$. The measurement at $105 °\mathrm{C}$ is still far away from the final saturation state after $24 \mathrm{h}$. The saturation value observed in the measurement at $125 °\mathrm{C}$ is approximately $-0.14\%$, and larger than the one of the measurement from [8], which was reached after approximately 5 years at $20 °\mathrm{C}$.

We notice a deviation between the results from Figure 5 and [57] for Zamak 5 as well. In [57], the alloy reaches the final stage at $120 °\mathrm{C}$ after $0.4$ days. This deviation can be observed as well comparing the measurement from [8] and the results from [57]. In [8] the specimens reach the final state after approximately five years at $20 °\mathrm{C}$ and in [57] only after one year of natural aging. Thus, additional tests of the mechanical behavior were performed in order to determine whether the specimens with 12 months of natural aging represent the completely aged state. A torsion test at room temperature and the fastest strain rate ${\dot{\gamma }}_1$ was performed with further specimens which were artificially aged at temperatures of $80 °\mathrm{C}$ for $24 \mathrm{h}$, $100 °\mathrm{C}$ for $24 \mathrm{h}$, and $120 °\mathrm{C}$ for $44 \mathrm{h}$. According to the investigations in [33] and the results in Figure 5a, the specimen with thermal treatment at $120 °\mathrm{C}$ for $44 \mathrm{h}$ reaches a completely aged state. The results of these new measurements are compared with the natural aged specimens in Figure 5b. Here, one can see that the material with one year of natural aging still does not correspond to the final aged state. The deviation from the results from [57] indicates further dependencies of the development of aging not only on the composition of the alloy, but also on other parameters of the die casting process which were not considered in this work. It is possible as well that the final microstructural state reached after a long heating process does not correspond completely to the same microstructural equilibrium state at lower temperatures. Since a more detailed investigation of this effect is out of the scope of this work, these dependencies are not investigated further. Nevertheless, it is taken into account that 12 months of natural aging does not correspond to the completely aged state.

3. Constitutive Model of Large Deformations for Zamak 5

In the following, we develop a phenomenological constitutive model under the assumption of large strains. Since the possible strains exceed 1% and applications require models for large displacements, a geometrical non-linear theory is necessary. Here, we draw on an overstress-type modeling concept, see, for example [58,59,60,61], where the total stress state is decomposed into a rate-independent equilibrium stress state, and a rate-dependent overstress part. The equilibrium stress part can be either obtained by an infinitely slow loading process, what is not feasible from a practical point of view, or by a multi-step relaxation path, where the termination points of relaxation—or some estimated points—determine the internal equilibrium state. Considering the experimental results, the material shows a strong viscoplastic material response, where the viscous effects are influenced by the temperature and the aging state. The equilibrium stress state is mostly influenced by the temperature. Moreover, a volume change due to aging (shrinkage) and due to temperature changes was observed, see [33]. The proposed constitutive model is based on the experience and knowledge of the small deformations model of [3]. In this case, the model structure is motivated by a rheological model, see Figure 6.

There are three main elements connected in series: an element describing the shrinkage due to aging (I), a thermal element (II) corresponding to the thermal expansion, and a mechanical element (III). In the mechanical element, the total stress $\tilde{\mathbf{T}}$ is decomposed into an equilibrium stress part ${\tilde{\mathbf{T}}}_{\mathrm{eq}}$ and an overstress part ${\tilde{\mathbf{T}}}_{\mathrm{ov}}$. $\tilde{\mathbf{T}}=\left(\det \mathbf{F}\right){\mathbf{F}}^{-1}\mathbf{T}{\mathbf{F}}^{-\mathrm{T}}$ defines the 2nd Piola-Kirchhoff stress tensor, $\mathbf{T}$ the Cauchy stress tensor, and $\mathbf{F}=\mathrm{Grad}{\overrightarrow{\mathit{\chi }}}_{\mathrm{R}}(\overrightarrow{X},t)$ the deformation gradient. $\overrightarrow{x}={\overrightarrow{\mathit{\chi }}}_{\mathrm{R}}(\overrightarrow{X},t)$ represents the motion of the material point $\overrightarrow{X}$ occupying the point $\overrightarrow{x}$ at time t in the current configuration. The equilibrium stress part ${\tilde{\mathbf{T}}}_{\mathrm{eq}}$ represents the equilibrium hysteresis and is composed of an elastic stress part related to the spring ${\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{e}}$ and a plastic component related to the plastic element ${\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{h}}$. The overstress ${\tilde{\mathbf{T}}}_{\mathrm{ov}}$ is modeled with a Maxwell-element, which is responsible for the rate-dependent response. Thus, the total stress is expressed by

$$\tilde{\mathbf{T}}={\tilde{\mathbf{T}}}_{\mathrm{eq}}+{\tilde{\mathbf{T}}}_{\mathrm{ov}}={\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{e}}+{\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{h}}+{\tilde{\mathbf{T}}}_{\mathrm{ov}}.$$

(1)

The decomposition of the stress into equilibrium stresses and overstresses shows advantages in the parameter identification procedure, since the parameters of the equilibrium stress part ${\tilde{\mathbf{T}}}_{\mathrm{eq}}$ can independently be identified of the parameters of the overstress ${\tilde{\mathbf{T}}}_{\mathrm{ov}}$. In the following, the rate-independent hysteretic part ${\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{h}}$ is called plastic part.

In order to consider the different physical effects, the deformation gradient is decomposed into an aging ${\mathbf{F}}_a$, a thermal ${\mathbf{F}}_{\Theta }$, and a mechanical part ${\mathbf{F}}_{\mathrm{M}}$,

$$\begin{array}{c}\hfill \mathbf{F}={\mathbf{F}}_a{\mathbf{F}}_{\Theta }{\mathbf{F}}_{\mathrm{M}}.\end{array}$$

(2)

The thermal and aging parts are assumed to behave purely volumetrical

$$\begin{array}{c}\hfill {\mathbf{F}}_{\Theta }={\phi }_{\Theta }(\Theta )\mathbf{I},{\mathbf{F}}_a={\phi }_a\left(a\right)\mathbf{I},\end{array}$$

(3)

where ${\phi }_{\Theta }(\Theta )$ and ${\phi }_a\left(a\right)$ are scalar functions of the temperature $\Theta$ and of the aging variable a

$$\begin{array}{cc}\hfill {\phi }_{\Theta }(\Theta )& =1+{\alpha }_{\Theta }(\Theta -{\Theta }_0),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill {\phi }_a\left(a\right)& =1-{\alpha }_{a}a,\hfill \end{array}$$

(5)

with ${\alpha }_{\Theta },{\alpha }_a>0$. As suggested by Flory [62], the mechanical deformation gradient ${\mathbf{F}}_{\mathrm{M}}$ is decomposed into a volumetric part ${\widehat{\mathbf{F}}}_{\mathrm{M}}$ and an isochoric part ${\overline{\mathbf{F}}}_{\mathrm{M}}$

$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathrm{M}}={\widehat{\mathbf{F}}}_{\mathrm{M}}{\overline{\mathbf{F}}}_{\mathrm{M}},\end{array}$$

(6)

with

$$\begin{array}{cc}\hfill {\widehat{\mathbf{F}}}_{\mathrm{M}}& ={\left(\det {\mathbf{F}}_{\mathrm{M}}\right)}^{1/3}\mathbf{I},\det {\widehat{\mathbf{F}}}_{\mathrm{M}}=\det {\mathbf{F}}_{\mathrm{M}},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill {\overline{\mathbf{F}}}_{\mathrm{M}}& ={\left(\det {\mathbf{F}}_{\mathrm{M}}\right)}^{-1/3}{\mathbf{F}}_{\mathrm{M}},\det {\overline{\mathbf{F}}}_{\mathrm{M}}=1.\hfill \end{array}$$

(8)

The Maxwell-element in the rheological model is considered with the decomposition of the mechanical part of the deformation gradient into an elastic part ${\mathbf{F}}_{\mathrm{e}}$ and an isochoric viscous part ${\overline{\mathbf{F}}}_{\mathrm{v}}$,

$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathrm{M}}={\mathbf{F}}_{\mathrm{e}}{\overline{\mathbf{F}}}_{\mathrm{v}}={\widehat{\mathbf{F}}}_{\mathrm{e}}{\overline{\mathbf{F}}}_{\mathrm{e}}{\overline{\mathbf{F}}}_{\mathrm{v}}.\end{array}$$

(9)

Moreover, the mechanical deformation gradient is decomposed in a similar way as in [49,63], into an elastic (kinematic hardening part) ${\mathbf{F}}_{\mathrm{k}}$ and a plastic component ${\overline{\mathbf{F}}}_{\mathrm{p}}$ to model the plastic element,

$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathrm{M}}={\mathbf{F}}_{\mathrm{k}}{\overline{\mathbf{F}}}_{\mathrm{p}}={\widehat{\mathbf{F}}}_{\mathrm{k}}{\overline{\mathbf{F}}}_{\mathrm{k}}{\overline{\mathbf{F}}}_{\mathrm{p}}.\end{array}$$

(10)

Due to the fact that there is no distinction between an elastic and a plastic range, such as it is the case in yield function based models, the interpretation of the components differs from [49,63]. The hysteretic response of the equilibrium stress is assigned to the term ${\mathbf{F}}_{\mathrm{k}}$, which corresponds to the strain-like variable $\mathbf{Y}$ in the small deformations case, see [3]. The condition of inelastic incompressibility is considered with the components ${\overline{\mathbf{F}}}_{\mathrm{v}}$ and ${\overline{\mathbf{F}}}_{\mathrm{p}}$ being purely isochoric. The decomposition of the deformation gradient is shown in Figure 7.

There ${\check{\mathit{\chi }}}_t$, ${\overline{\mathit{\chi }}}_t$, ${\widehat{\mathit{\chi }}}_t$, ${\breve{\mathit{\chi }}}_t$, and ${\tilde{\mathit{\chi }}}_t$ represent the intermediate configurations caused by the decomposition of $\mathbf{F}$.

The split of the deformation gradient allows the definition of a number of strain measures, which are formulated relative to particular configurations. Considering the right Cauchy-Green tensor $\mathbf{C}={\mathbf{F}}^{\mathrm{T}}\mathbf{F}$, the Green strain tensor can be decomposed as well,

$$\begin{array}{c}\hfill \mathbf{E}={\displaystyle \frac{1}{2}}(\mathbf{C}-\mathbf{I})={\displaystyle \frac{1}{2}}(\mathbf{C}-{\mathbf{C}}_{\Theta \mathrm{M}})+{\displaystyle \frac{1}{2}}({\mathbf{C}}_{\Theta \mathrm{M}}-{\mathbf{C}}_{\mathrm{M}})+{\displaystyle \frac{1}{2}}({\mathbf{C}}_{\mathrm{M}}-\mathbf{I}).\end{array}$$

(11)

Using ${\mathbf{F}}_{\Theta \mathrm{M}}={\mathbf{F}}_{\Theta }{\mathbf{F}}_{\mathrm{M}}={\phi }_{\Theta }{\mathbf{F}}_{\mathrm{M}}$ and $\mathbf{F}={\phi }_a{\phi }_{\Theta }{\mathbf{F}}_{\mathrm{M}}$, the components right Cauchy-Green type tensors ${\mathbf{C}}_{\mathrm{M}}$, ${\mathbf{C}}_{\Theta \mathrm{M}}$ and $\mathbf{C}$ are defined by

$$\begin{array}{cc}\hfill {\mathbf{C}}_{\mathrm{M}}& ={\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{M}},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\mathbf{C}}_{\Theta \mathrm{M}}& ={\mathbf{F}}_{\Theta \mathrm{M}}^{\mathrm{T}}{\mathbf{F}}_{\Theta \mathrm{M}}={\phi }_{\Theta }^2{\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{M}},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \mathbf{C}& ={\mathbf{F}}^{\mathrm{T}}\mathbf{F}={\phi }_a^2{\phi }_{\Theta }^2{\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{M}}.\hfill \end{array}$$

(14)

Inserting definitions (12)–(14) into Equation (11) results in

$$\begin{array}{c}\hfill \mathbf{E}={\displaystyle \frac{1}{2}}({\phi }_a^2-1){\phi }_{\Theta }^2{\mathbf{C}}_{\mathrm{M}}+{\displaystyle \frac{1}{2}}({\phi }_{\Theta }^2-1){\mathbf{C}}_{\mathrm{M}}+{\displaystyle \frac{1}{2}}({\mathbf{C}}_{\mathrm{M}}-\mathbf{I})={\mathbf{E}}_a+{\mathbf{E}}_{\Theta }+{\mathbf{E}}_{\mathrm{M}}.\end{array}$$

(15)

Thus, the Green strain tensor decomposes into three different components: the aging strain ${\mathbf{E}}_a$, the thermal strain ${\mathbf{E}}_{\Theta }$, and the mechanical strain ${\mathbf{E}}_{\mathrm{M}}$,

$$\begin{array}{c}\hfill {\mathbf{E}}_a={\displaystyle \frac{1}{2}}({\phi }_a^2-1){\phi }_{\Theta }^2{\mathbf{C}}_{\mathrm{M}},{\mathbf{E}}_{\Theta }={\displaystyle \frac{1}{2}}({\phi }_{\Theta }^2-1){\mathbf{C}}_{\mathrm{M}},{\mathbf{E}}_{\mathrm{M}}={\displaystyle \frac{1}{2}}({\mathbf{C}}_{\mathrm{M}}-\mathbf{I}).\end{array}$$

(16)

The mechanical strain tensor ${\mathbf{E}}_{\mathrm{M}}$ is assigned both to the Maxwell-element as well as the rate-independent endochronic element, see Figure 6,

$$\begin{array}{cc}\hfill {\mathbf{E}}_{\mathrm{M}}& ={\displaystyle \frac{1}{2}}({\mathbf{C}}_{\mathrm{M}}-\mathbf{I})={\displaystyle \frac{1}{2}}({\mathbf{C}}_{\mathrm{M}}-{\overline{\mathbf{C}}}_{\mathrm{v}})+{\displaystyle \frac{1}{2}}({\overline{\mathbf{C}}}_{\mathrm{v}}-\mathbf{I})={\mathbf{E}}_{\mathrm{e}}+{\mathbf{E}}_{\mathrm{v}},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\mathbf{E}}_{\mathrm{M}}& ={\displaystyle \frac{1}{2}}({\mathbf{C}}_{\mathrm{M}}-\mathbf{I})={\displaystyle \frac{1}{2}}({\mathbf{C}}_{\mathrm{M}}-{\overline{\mathbf{C}}}_{\mathrm{p}})+{\displaystyle \frac{1}{2}}({\overline{\mathbf{C}}}_{\mathrm{p}}-\mathbf{I})={\mathbf{E}}_{\mathrm{k}}+{\mathbf{E}}_{\mathrm{p}},\hfill \end{array}$$

(18)

with ${\overline{\mathbf{C}}}_{\mathrm{v}}={\overline{\mathbf{F}}}_{\mathrm{v}}^{\mathrm{T}}{\overline{\mathbf{F}}}_{\mathrm{v}}$ and ${\overline{\mathbf{C}}}_{\mathrm{p}}={\overline{\mathbf{F}}}_{\mathrm{p}}^{\mathrm{T}}{\overline{\mathbf{F}}}_{\mathrm{p}}$. ${\overline{\mathbf{C}}}_{\mathrm{v}}={\left(\det {\mathbf{C}}_{\mathrm{v}}\right)}^{-1/3}{\mathbf{C}}_{\mathrm{v}}$ and ${\overline{\mathbf{C}}}_{\mathrm{p}}={\left(\det {\mathbf{C}}_{\mathrm{p}}\right)}^{-1/3}{\mathbf{C}}_{\mathrm{p}}$ are unimodular tensors, which we anticipate here due to the later choice of the evolutionary equations. In other words, the mechanical strain decomposes into an elastic and a viscous part (Maxwell-element) ${\mathbf{E}}_{\mathrm{M}}={\mathbf{E}}_{\mathrm{e}}+{\mathbf{E}}_{\mathrm{v}}$, or an elastic and plastic part ${\mathbf{E}}_{\mathrm{M}}={\mathbf{E}}_{\mathrm{k}}+{\mathbf{E}}_{\mathrm{p}}$ (rate-independent plastic element),

$$\begin{array}{cc}\hfill {\mathbf{E}}_{\mathrm{e}}& ={\displaystyle \frac{1}{2}}({\mathbf{C}}_{\mathrm{M}}-{\overline{\mathbf{C}}}_{\mathrm{v}}),{\mathbf{E}}_{\mathrm{v}}={\displaystyle \frac{1}{2}}({\overline{\mathbf{C}}}_{\mathrm{v}}-\mathbf{I}),\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\mathbf{E}}_{\mathrm{k}}& ={\displaystyle \frac{1}{2}}({\mathbf{C}}_{\mathrm{M}}-{\overline{\mathbf{C}}}_{\mathrm{p}}),{\mathbf{E}}_{\mathrm{p}}={\displaystyle \frac{1}{2}}({\overline{\mathbf{C}}}_{\mathrm{p}}-\mathbf{I}).\hfill \end{array}$$

(20)

These strain measures can be formulated relative to different configurations. For example, the mechanical part of the strain tensor ${\mathit{\epsilon }}_{\mathrm{M}}={\displaystyle \frac{1}{2}}({\mathbf{F}}_{\mathrm{e}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{e}}-{\overline{\mathbf{F}}}_{\mathrm{v}}^{-\mathrm{T}}{\overline{\mathbf{F}}}_{\mathrm{v}}^{-1})$ relative to the elastic intermediate configuration ${\check{\mathit{\chi }}}_t$ can be decomposed in a purely elastic component ${\mathit{\epsilon }}_{\mathrm{e}}$ and a purely viscous part ${\mathit{\epsilon }}_{\mathrm{v}}$,

$$\begin{array}{c}\hfill {\mathit{\epsilon }}_{\mathrm{e}}={\displaystyle \frac{1}{2}}({\mathbf{F}}_{\mathrm{e}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{e}}-\mathbf{I}),{\mathit{\epsilon }}_{\mathrm{v}}={\displaystyle \frac{1}{2}}(\mathbf{I}-{\overline{\mathbf{F}}}_{\mathrm{v}}^{-\mathrm{T}}{\overline{\mathbf{F}}}_{\mathrm{v}}^{-1}),\end{array}$$

(21)

see Figure A1. Furthermore, quantities relative to the plastic configuration ${\tilde{\mathit{\chi }}}_t$ can be defined. Namely the mechanical strain tensor ${\tilde{\mathit{\epsilon }}}_{\mathrm{M}}={\displaystyle \frac{1}{2}}({\mathbf{F}}_{\mathrm{k}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{k}}-{\overline{\mathbf{F}}}_{\mathrm{p}}^{-\mathrm{T}}{\overline{\mathbf{F}}}_{\mathrm{p}}^{-1})$, which can be decomposed into a component depending exclusively on ${\mathbf{F}}_{\mathrm{k}}$ and one depending only on ${\overline{\mathbf{F}}}_{\mathrm{p}}$

$$\begin{array}{c}\hfill {\tilde{\mathit{\epsilon }}}_{\mathrm{k}}={\displaystyle \frac{1}{2}}({\mathbf{F}}_{\mathrm{k}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{k}}-\mathbf{I}),{\tilde{\mathit{\epsilon }}}_{\mathrm{p}}={\displaystyle \frac{1}{2}}(\mathbf{I}-{\overline{\mathbf{F}}}_{\mathrm{p}}^{-\mathrm{T}}{\overline{\mathbf{F}}}_{\mathrm{p}}^{-1}).\end{array}$$

(22)

Figure A1 of the Appendix A shows the transformations of the strain measures relative to the different configurations.

3.1. Resulting Stresses and Evolution Equations

The constitutive model has to satisfy the Clausius-Duhem inequality for fulfilling thermodynamical consistence. In order to minimize the presentation here, the proof is compiled in Appendix C. We assume a specific free energy, which is divided into an elastic ${\psi }_{\mathrm{eq}}^{\mathrm{e}}$ (spring element), a plastic ${\psi }_{\mathrm{eq}}^{\mathrm{h}}$ (rate-independent plastic component), an overstress ${\psi }_{\mathrm{ov}}$ (Maxwell-element), a thermal ${\psi }_{\Theta }$, and an aging part ${\psi }_a$

$$\begin{array}{c}\hfill \psi ({\mathbf{E}}_{\mathrm{M}},{\tilde{\mathit{\epsilon }}}_{\mathrm{k}},{\mathit{\epsilon }}_{\mathrm{e}},k,a,\Theta )={\psi }_{\mathrm{eq}}^{\mathrm{e}}\left({\mathbf{E}}_{\mathrm{M}}\right)+{\psi }_{\mathrm{eq}}^{\mathrm{h}}({\tilde{\mathit{\epsilon }}}_{\mathrm{k}},\Theta )+{\psi }_{\mathrm{ov}}({\mathit{\epsilon }}_{\mathrm{e}},k)+{\psi }_{\Theta }(\Theta )+{\psi }_a\left(a\right).\end{array}$$

(23)

The internal variable k represents a softening variable having the goal of considering a change of the initial slopes in loading and unloading. This was observed in the experiments, see [3]. For the elastic component ${\psi }_{\mathrm{eq}}^{\mathrm{e}}\left({\mathbf{E}}_{\mathrm{M}}\right)$ of the specific free energy, an expression of near-incompressibility is chosen. Here, the specific free energy is decomposed into a volumetric and an isochoric component

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{R}}{\psi }_{\mathrm{eq}}^{\mathrm{e}}\left({\mathbf{E}}_{\mathrm{M}}\right)& =U\left(J_{\mathrm{M}}\right)+\overline{\omega }\left({\overline{\mathbf{C}}}_{\mathrm{M}}\right),\hfill \end{array}$$

(24)

${\rho }_{\mathrm{R}}$ defines the density in the reference configuration. The volume changing part of the energy function $U\left(J_{\mathrm{M}}\right)$ is chosen according to [64] depending on the determinant of the mechanical deformation gradient $J_{\mathrm{M}}:=\det {\mathbf{F}}_{\mathrm{M}}$, see Equation (A47). For the isochoric part component a Neo-Hooke ansatz is chosen, see Equation (A48), depending on the unimodular mechanical right Cauchy-Green tensor ${\overline{\mathbf{C}}}_{\mathrm{M}}$. This term is equal to the total unimodular Cauchy-Green tensor ${\overline{\mathbf{C}}}_{\mathrm{M}}=\overline{\mathbf{C}}={\left(\det {\mathbf{C}}_{\mathrm{M}}\right)}^{-1/3}{\mathbf{C}}_{\mathrm{M}}$ (the aging and the thermal part are purely volumetric).

The plastic element and the Maxwell-element are both able to reproduce a hysteretic behavior, with the difference that the Maxwell-element is rate-dependent. The rate independence is included in the plastic element with definition (38) due to the proportionality of the strain rate ${\stackrel{\Delta }{\tilde{\mathit{\epsilon }}}}_{\mathrm{p}}$ to the rate of the arc-length ${\dot{s}}_{\mathrm{M}}$. For these elements, a Neo-Hookean ansatz is chosen too, in which the parameters ${\mu }_{\mathrm{eq}}(\Theta )$ and ${\mu }_{\mathrm{ov}}\left(k\right)$ are affected by the temperature and the softening respectively

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{R}}{\psi }_{\mathrm{eq}}^{\mathrm{h}}({\tilde{\mathit{\epsilon }}}_{\mathrm{k}},\Theta )& ={\mu }_{\mathrm{eq}}(\Theta )(I_{{\overline{\mathbf{C}}}_{\mathrm{k}}}-3),{\overline{\mathbf{C}}}_{\mathrm{k}}={\left(\det {\mathbf{C}}_{\mathrm{k}}\right)}^{-1/3}{\mathbf{C}}_{\mathrm{k}},\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill {\rho }_{\mathrm{R}}{\psi }_{\mathrm{ov}}({\mathit{\epsilon }}_{\mathrm{e}},k)& ={\mu }_{\mathrm{ov}}\left(k\right)(I_{{\overline{\mathbf{C}}}_{\mathrm{e}}}-3),{\overline{\mathbf{C}}}_{\mathrm{e}}={\left(\det {\mathbf{C}}_{\mathrm{e}}\right)}^{-1/3}{\mathbf{C}}_{\mathrm{e}}.\hfill \end{array}$$

(26)

${\mathrm{I}}_{\mathbf{A}}=\mathrm{tr}\mathbf{A}=\mathbf{A}·\mathbf{I}=a_k^k$ defines the trace operator of a second-order tensor. For the aging part of the free energy, the expression

$$\begin{array}{c}\hfill {\rho }_{\mathrm{R}}{\psi }_a\left(a\right)=C_1-C_{2}a,\end{array}$$

(27)

is chosen and for the thermal part, the expression proposed by [65] is selected

$$\begin{array}{c}\hfill {\rho }_{\mathrm{R}}{\psi }_{\Theta }(\Theta )=c_{\mathrm{p}1}\left((\Theta -{\Theta }_0)-\Theta \ln {\displaystyle \frac{\Theta }{{\Theta }_0}}(1-c_{\mathrm{p}2}\Theta )-{\displaystyle \frac{1}{2}}({\Theta }^2-{\Theta }_0^2)\right).\end{array}$$

(28)

With this function, the specific heat capacity obtains a linear dependence on the temperature, which corresponds to the experimental results shown in [33].

From considering the Clausius-Duhem inequality, see the Appendix C, and the chosen free energy functions, the following expressions arise for the components of the stress tensors

$$\begin{array}{cc}\hfill {{\mathbf{S}}_{\mathrm{M}}^{\mathrm{e}}}_{\mathrm{eq}}& ={\displaystyle \frac{{\rho }_{\mathrm{R}}}{{\phi }_{\Theta }^2{\phi }_a^2}}{\mathbf{F}}_{\mathrm{M}}{\displaystyle \frac{\partial \psi }{\partial {\mathbf{E}}_{\mathrm{M}}}}{\mathbf{F}}_{\mathrm{M}}^{\mathrm{T}}={\left({\phi }_{\Theta }{\phi }_a\right)}^{-2}{\displaystyle \frac{K}{10}}(J_{\mathrm{M}}^5-J_{\mathrm{M}}^{-5})\mathbf{I}+{\left({\phi }_{\Theta }{\phi }_a\right)}^{-2}2c_{10}{\overline{\mathbf{B}}}_{\mathrm{M}}^{\mathrm{D}},\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\check{\mathbf{S}}}_{\mathrm{ov}}& ={\displaystyle \frac{{\rho }_{\mathrm{R}}}{{\phi }_{\Theta }^2{\phi }_a^2}}{\displaystyle \frac{\partial \psi }{\partial {\mathit{\epsilon }}_{\mathrm{e}}}}=2{\mu }_{\mathrm{eq}}{\left({\phi }_{\Theta }{\phi }_a\right)}^{-2}{\mathbf{C}}_{\mathrm{k}}^{-1}{\overline{\mathbf{C}}}_{\mathrm{k}}^{\mathrm{D}},\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{S}}}_{\mathrm{eq}}^{\mathrm{h}}& ={\displaystyle \frac{{\rho }_{\mathrm{R}}}{{\phi }_{\Theta }^2{\phi }_a^2}}{\displaystyle \frac{\partial \psi }{\partial {\tilde{\mathit{\epsilon }}}_{\mathrm{k}}}}=2{\mu }_{\mathrm{ov}}{\left({\phi }_{\Theta }{\phi }_a\right)}^{-2}{\mathbf{C}}_{\mathrm{e}}^{-1}{\overline{\mathbf{C}}}_{\mathrm{e}}^{\mathrm{D}},\hfill \end{array}$$

(31)

see Figure A3 for the different stress tensors. Here, the exponent D represents the deviator of a second-order tensor, ${\mathbf{A}}^D=\mathbf{A}-\left(\mathrm{tr}\mathbf{A}\right)/3\mathbf{I}$. For the implementation of the model, the second Piola-Kirchhoff-type stresses relative to the reference configuration are of interest. Making use of the relations ${\mathbf{F}}_{\mathrm{M}}={\left({\phi }_a{\phi }_{\Theta }\right)}^{-1}\mathbf{F}$, $J_{\mathrm{M}}={\left({\phi }_a{\phi }_{\Theta }\right)}^{-3}J$, $J:=\det \mathbf{F}$, one obtains

$$\begin{array}{cc}\hfill {\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{e}}& ={\mathbf{F}}_{\mathrm{M}}^{-1}{{\mathbf{S}}_{\mathrm{M}}^{\mathrm{e}}}_{\mathrm{eq}}{\mathbf{F}}_{\mathrm{M}}^{-\mathrm{T}}={\left({\phi }_{\Theta }{\phi }_a\right)}^{-2}\left({\displaystyle \frac{K}{10}}(J_{\mathrm{M}}^5-J_{\mathrm{M}}^{-5}){\mathbf{C}}_{\mathrm{M}}^{-1}+2c_{10}{J}_{\mathrm{M}}^{-2/3}{\mathbf{C}}_{\mathrm{M}}^{-1}{\mathbf{C}}_{\mathrm{M}}^{\mathrm{D}}\right)\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{h}}& ={\overline{\mathbf{F}}}_{\mathrm{p}}^{-1}{{\mathbf{S}}_{\mathrm{M}}^{\mathrm{h}}}_{\mathrm{eq}}{\overline{\mathbf{F}}}_{\mathrm{p}}^{-\mathrm{T}}=2{\mu }_{\mathrm{eq}}{J}^{-2/3}{\mathbf{C}}^{-1}{\left(\mathbf{C}{\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\right)}^{\mathrm{D}},\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{T}}}_{\mathrm{ov}}& ={\overline{\mathbf{F}}}_{\mathrm{v}}^{-1}{\check{\mathbf{S}}}_{\mathrm{ov}}{\overline{\mathbf{F}}}_{\mathrm{v}}^{-\mathrm{T}}=2{\mu }_{\mathrm{ov}}{J}^{-2/3}{\mathbf{C}}^{-1}{\left(\mathbf{C}{\overline{\mathbf{C}}}_{\mathrm{v}}^{-1}\right)}^{\mathrm{D}}.\hfill \end{array}$$

(34)

Moreover, after a push-forward operation,

$$\begin{array}{cc}\hfill {\mathbf{S}}_{\mathrm{eq}}^{\mathrm{e}}& =\mathbf{F}{\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{e}}{\mathbf{F}}^{\mathrm{T}}={\mathbf{S}}_{\mathrm{vol}}^{\mathrm{e}}+{\mathbf{S}}_{\mathrm{iso}}^{\mathrm{e}}={\displaystyle \frac{K}{10}}\left(J_{\mathrm{M}}^5-J_{\mathrm{M}}^{-5}\right)\mathbf{I}+2c_{10}{\overline{\mathbf{B}}}_{\mathrm{M}}^{\mathrm{D}}.\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\mathbf{S}}_{\mathrm{eq}}^{\mathrm{h}}& =\mathbf{F}{\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{h}}{\mathbf{F}}^{\mathrm{T}}=2{\mu }_{\mathrm{eq}}{\overline{\mathbf{B}}}_{\mathrm{k}}^{\mathrm{D}},\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill {\mathbf{S}}_{\mathrm{ov}}& =\mathbf{F}{\tilde{\mathbf{T}}}_{\mathrm{ov}}{\mathbf{F}}^{\mathrm{T}}=2{\mu }_{\mathrm{ov}}{\overline{\mathbf{B}}}_{\mathrm{e}}^{\mathrm{D}},\hfill \end{array}$$

(37)

one can see that the plastic stress part ${\mathbf{S}}_{\mathrm{eq}}^{\mathrm{h}}$ and the overstress ${\mathbf{S}}_{\mathrm{ov}}$ in the current configuration are purely deviatoric. ${\mathbf{S}}_{\mathrm{vol}}^{\mathrm{e}}$ is the only volumetric stress contribution and represents the only part connected to elastic volumetric processes, see Equation (A59).

Additionally, expressions for the evolution equations have to be found. The following relations ensure the thermodynamical consistence of the model

$$\begin{array}{cc}\hfill {\stackrel{\Delta }{\tilde{\mathit{\epsilon }}}}_{\mathrm{p}}& ={\beta }_{\mathbf{Y}}{\dot{s}}_{\mathrm{M}}{\mathbf{C}}_{\mathrm{k}}{\tilde{\mathbf{S}}}_{\mathrm{eq}}^{\mathrm{h}},{\beta }_{\mathbf{Y}},{\dot{s}}_{\mathrm{M}}>0,\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill {\stackrel{\Delta }{\mathit{\epsilon }}}_{\mathrm{v}}& ={\displaystyle \frac{1}{\widehat{\eta }}}{\mathbf{C}}_{\mathrm{e}}{\check{\mathbf{S}}}_{\mathrm{ov}},\widehat{\eta }>0,\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill \dot{k}& \ge 0,{\displaystyle \frac{\partial \psi }{\partial k}}\le 0,\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill \dot{a}& ={\beta }_g\left({\phi }_a{\displaystyle \frac{\partial {\phi }_a}{\partial a}}(\mathrm{tr}\breve{\mathbf{S}})-{\displaystyle \frac{\partial \psi }{\partial a}}\right),{\beta }_g>0,\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill {\overrightarrow{q}}_R& =-{\mathit{\kappa }}_{\Theta }^{\mathrm{R}}\mathrm{Grad}\Theta ,\mathrm{and}{\overrightarrow{g}}_R=\mathrm{Grad}\Theta .\hfill \end{array}$$

(42)

Here, ${\stackrel{\Delta }{\tilde{\mathit{\epsilon }}}}_{\mathrm{p}}$ and ${\stackrel{\Delta }{\mathit{\epsilon }}}_{\mathrm{v}}$ are the covariant Oldroyd rates of the plastic and viscous strains, respectively. In the Maxwell-element, the parameter $\widehat{\eta }$ represents the viscosity, which is modeled later on. The plastic element has the goal to describe a rate-independent process. To this end, the term ${\stackrel{\Delta }{\tilde{\mathit{\epsilon }}}}_{\mathrm{p}}$ is defined to be proportional to the rate of the mechanical arc-length ${\dot{s}}_{\mathrm{M}}$

$$\begin{array}{c}\hfill {\dot{s}}_{\mathrm{M}}=\sqrt{{\stackrel{\Delta }{\overline{\tilde{\mathit{\epsilon }}}}}_{\mathrm{M}}·{\stackrel{\Delta }{\overline{\tilde{\mathit{\epsilon }}}}}_{\mathrm{M}}},\end{array}$$

(43)

with ${\tilde{\mathit{\epsilon }}}_{\mathrm{M}}={\displaystyle \frac{1}{2}}({\mathbf{F}}_{\mathrm{k}}^{\mathrm{T}}{\mathbf{F}}_{\mathrm{k}}-{\overline{\mathbf{F}}}_{\mathrm{p}}^{-\mathrm{T}}{\overline{\mathbf{F}}}_{\mathrm{p}}^{-1})$ the mechanical strain expressed with quantities of the intermediate configuration ${\tilde{\mathit{\chi }}}_t$, and ${\overline{\tilde{\mathit{\epsilon }}}}_{\mathrm{M}}$ as its isochoric part. The parameter ${\beta }_{\mathbf{Y}}$ from Equation (38) is defined by

$$\begin{array}{c}\hfill {\beta }_{\mathbf{Y}}={\displaystyle \frac{\widehat{b}}{{\mu }_{\mathrm{eq}}(\Theta )}},\end{array}$$

(44)

where $\widehat{b}$ is a material parameter and ${\mu }_{\mathrm{eq}}(\Theta )$ a function of the temperature. The derived stresses from Equation (31) and Equation (30) are inserted into Equations (38) and (39), which leads to

$$\begin{array}{cc}\hfill {\stackrel{\Delta }{\tilde{\mathit{\epsilon }}}}_{\mathrm{p}}& =2{\mu }_{\mathrm{eq}}{\beta }_{\mathbf{Y}}{\left({\phi }_{\Theta }{\phi }_a\right)}^{-2}{\dot{s}}_{\mathrm{M}}{\overline{\mathbf{C}}}_{\mathrm{k}}^{\mathrm{D}},\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\stackrel{\Delta }{\mathit{\epsilon }}}_{\mathrm{v}}& ={\displaystyle \frac{2{\mu }_{\mathrm{ov}}}{\widehat{\eta }}}{\left({\phi }_{\Theta }{\phi }_a\right)}^{-2}{\overline{\mathbf{C}}}_{\mathrm{e}}^{\mathrm{D}}.\hfill \end{array}$$

(46)

With a pull-back operation and considering the relations ${\mathbf{F}}_{\mathrm{M}}={\mathbf{F}}_{\mathrm{e}}{\overline{\mathbf{F}}}_{\mathrm{v}}={\mathbf{F}}_{\mathrm{k}}{\overline{\mathbf{F}}}_{\mathrm{p}}$ and ${\overline{\mathbf{C}}}_{\mathrm{M}}=\overline{\mathbf{C}}=J^{-2/3}\mathbf{C}$, these equations can be expressed with quantities relative to the reference configuration

$$\begin{array}{cc}\hfill {\dot{\mathbf{C}}}_{\mathrm{p}}& =4{\left({\phi }_a{\phi }_{\Theta }\right)}^{-2}{J}^{-2/3}{\beta }_{\mathbf{Y}}{\mu }_{\mathrm{eq}}{\dot{s}}_{\mathrm{M}}{\overline{\mathbf{C}}}_{\mathrm{p}}{\left({\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\mathbf{C}\right)}^{\mathrm{D}},\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill {\dot{\mathbf{C}}}_{\mathrm{v}}& ={\displaystyle \frac{4{\mu }_{\mathrm{ov}}}{\widehat{\eta }}}{\left({\phi }_a{\phi }_{\Theta }\right)}^{-2}{J}^{-2/3}{\overline{\mathbf{C}}}_{\mathrm{v}}{\left({\overline{\mathbf{C}}}_{\mathrm{v}}^{-1}\mathbf{C}\right)}^{\mathrm{D}}.\hfill \end{array}$$

(48)

The tensors ${\overline{\mathbf{C}}}_{\mathrm{p}}{\left({\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\mathbf{C}\right)}^{\mathrm{D}}$ and ${\overline{\mathbf{C}}}_{\mathrm{v}}{\left({\overline{\mathbf{C}}}_{\mathrm{v}}^{-1}\mathbf{C}\right)}^{\mathrm{D}}$ are symmetric since they can be expressed by the sum of symmetric tensors

$$\begin{array}{cc}\hfill {\overline{\mathbf{C}}}_{\mathrm{p}}{\left({\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\mathbf{C}\right)}^{\mathrm{D}}& ={\overline{\mathbf{C}}}_{\mathrm{p}}\left({\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\mathbf{C}-{\displaystyle \frac{1}{3}}\left(\mathrm{tr}\left({\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\mathbf{C}\right)\mathbf{I}\right)\right)=\mathbf{C}-{\displaystyle \frac{1}{3}}\left(\mathrm{tr}\left({\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\mathbf{C}\right)\right){\overline{\mathbf{C}}}_{\mathrm{p}},\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill {\overline{\mathbf{C}}}_{\mathrm{v}}{\left({\overline{\mathbf{C}}}_{\mathrm{v}}^{-1}\mathbf{C}\right)}^{\mathrm{D}}& ={\overline{\mathbf{C}}}_{\mathrm{v}}\left({\overline{\mathbf{C}}}_{\mathrm{v}}^{-1}\mathbf{C}-{\displaystyle \frac{1}{3}}\left(\mathrm{tr}\left({\overline{\mathbf{C}}}_{\mathrm{v}}^{-1}\mathbf{C}\right)\mathbf{I}\right)\right)=\mathbf{C}-{\displaystyle \frac{1}{3}}\left(\mathrm{tr}\left({\overline{\mathbf{C}}}_{\mathrm{v}}^{-1}\mathbf{C}\right)\right){\overline{\mathbf{C}}}_{\mathrm{v}}.\hfill \end{array}$$

(50)

To determine the rate of the arc-length ${\dot{s}}_{\mathrm{M}}$ in dependence of the variables of the reference configuration, the isochoric mechanical strain rate ${\stackrel{\Delta }{\overline{\tilde{\mathit{\epsilon }}}}}_{\mathrm{M}}$ is considered first

$$\begin{array}{c}\hfill {\stackrel{\Delta }{\overline{\tilde{\mathit{\epsilon }}}}}_{\mathrm{M}}={\overline{\mathbf{F}}}_{\mathrm{p}}^{-\mathrm{T}}{\dot{\overline{\mathbf{E}}}}_{\mathrm{M}}{\overline{\mathbf{F}}}_{\mathrm{p}}^{-1}.\end{array}$$

(51)

The isochoric strain rate is equal to

$$\begin{array}{c}\hfill {\dot{\overline{\mathbf{E}}}}_{\mathrm{M}}={\displaystyle \frac{J^{-2/3}}{2}}\mathbf{C}{\left({\mathbf{C}}^{-1}\dot{\mathbf{C}}\right)}^{\mathrm{D}}.\end{array}$$

(52)

Thus, considering Equation (52), the rate of the arc-length ${\dot{s}}_{\mathrm{M}}$ is given by

$$\begin{array}{c}\hfill {\dot{s}}_{\mathrm{M}}=\sqrt{{\stackrel{\Delta }{\overline{\tilde{\mathit{\epsilon }}}}}_{\mathrm{M}}·{\stackrel{\Delta }{\overline{\tilde{\mathit{\epsilon }}}}}_{\mathrm{M}}}={\displaystyle \frac{J^{-2/3}}{2}}\sqrt{{\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\mathbf{C}{\left({\mathbf{C}}^{-1}\dot{\mathbf{C}}\right)}^{\mathrm{D}}{\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}·\mathbf{C}{\left({\mathbf{C}}^{-1}\dot{\mathbf{C}}\right)}^{\mathrm{D}}}.\end{array}$$

(53)

The following evolution equation is chosen for the softening variable k similarly to the small deformations case, see [3],

$$\begin{array}{c}\hfill \dot{k}={\alpha }_k(1-k){\dot{s}}_{\mathrm{M}}.\end{array}$$

(54)

The softening is a growing rate-independent variable with values between 0 and 1. The value 0 correspond to the initial, undeformed state and 1 the state, in which the maximal softening is reached. Finally, an evolution equation for the aging variable a has to be found. From Equations (41) and (A68), one obtains

$$\begin{array}{c}\hfill \dot{a}={\beta }_a\left(C_2+{\phi }_a{\phi }_{\Theta }^2{\displaystyle \frac{\partial {\phi }_a}{\partial a}}\left(\mathrm{tr}{\mathbf{S}}_{\mathrm{M}}\right)\right).\end{array}$$

(55)

The aging rate is proportional to a parameter $C_2$ and a term depending on the trace of the stress ${\mathbf{S}}_{\mathrm{M}}$. The sign of $\mathrm{tr}{\mathbf{S}}_{\mathrm{M}}$ determines whether the aging process is accelerated or decelerated. A negative aging rate would correspond to a reversible aging process, which is physically unrealistic. Because of that, the term $C_2$ has to be larger than ${\phi }_a{\phi }_{\Theta }^2(\partial {\phi }_a/\partial a)\left(\mathrm{tr}{\mathbf{S}}_{\mathrm{M}}\right)$. At the moment, the literature offers no information that allows to quantify the influence of mechanical deformation or stresses on the aging process of die casting alloys. The microstructure of a ZnAl20.2Cu1.8 alloy after a tensile test was investigated with X-ray diffraction and SEM in [24]. Here, it was observed that the tensile stress was associated to microstructural changes such as phase transformations. Zhu [13] concluded that stresses accelerate the decomposition of metastable phases in tensile tests. Unfortunately, the studies were performed at temperatures over $100 °\mathrm{C}$ and the microstructural changes were not quantified. There is also no information about the compression or torsion behavior. Moreover, Deschamps et al. [66] and Cerri and Leo [67] found an increase in the aging rate with the plastic deformation for aluminum alloys. Since aging takes place in a much larger timescale than the mechanical processes in the application range of the alloy and since there is no information at hand to characterize the influence of the stress on aging, we assume that the term $C_2$ is much larger than ${\phi }_a{\phi }_{\Theta }^2(\partial {\phi }_a/\partial a)\left(\mathrm{tr}{\mathbf{S}}_{\mathrm{M}}\right)$, analogously to the discussion in Lion et al. [68] for the field of curing or, more recently, in [4]. Thus, the following expression is considered

$$\begin{array}{c}\hfill \dot{a}\approx {\beta }_a{C}_2.\end{array}$$

(56)

In order to describe the aging state of the material, the aging variable is defined as a process variable with values between 0 an 1. The value 0 corresponds to the initial unaged stage, while 1 represents the completely aged material. In the measurements of Figure 5 one can see that the isothermal aging process is stronger at the beginning of the process and decelerates when approaching saturation and occurs faster at higher temperatures. In order to fulfill these criteria, the function

$$\begin{array}{c}\hfill {\beta }_a(\Theta ,a)={\displaystyle \frac{{\beta }_{a1}(\Theta )}{C_2}}{\displaystyle \frac{(1-a)}{({\beta }_{a2}+a)}},\end{array}$$

(57)

is proposed, see also [4]. The function ${\beta }_{a1}(\Theta )$ includes the temperature dependence in the development of the aging. An Arrhenius ansatz is chosen for this function, which is a common approach for temperature-driven diffusion processes in metals. Thus, the evolution equation of the aging variable reads

$$\begin{array}{c}\hfill \dot{a}={\beta }_{a1}(\Theta ){\displaystyle \frac{1-a}{{\beta }_{a2}+a}},{\beta }_{a1}(\Theta )={\beta }_{a11}\exp (-{\beta }_{a12}/\Theta ),\end{array}$$

(58)

with ${\beta }_{a12}=E_a/R$, $E_a$ is the activation energy and R the universal gas constant.

3.2. Functions for the Material Parameters

Based on the developed functions for the small deformation case, see [3,4], the functions ${\mu }_{\mathrm{eq}}(\Theta )$, $\eta (\dot{\mathbf{C}},\mathbf{C},\Theta ,{\overline{\mathbf{C}}}_{\mathrm{p}},a)$, and ${\mu }_{\mathrm{ov}}\left(k\right)$ are adapted. The change in the initial slope of the plastic equilibrium stress is modeled by the parameter ${\mu }_{\mathrm{eq}}(\Theta )$, analogously to the parameter $c(\Theta )$ in the small deformation model in [3]

$$\begin{array}{c}\hfill {\mu }_{\mathrm{eq}}(\Theta )={\widehat{c}}_1({\widehat{c}}_2-\tanh \left({\widehat{c}}_3(\Theta -{\widehat{\Theta }}_c)\right)).\end{array}$$

(59)

The following expression is chosen for the shear modulus of the overstress ${\mu }_{\mathrm{ov}}$ depending on the softening variable k

$$\begin{array}{c}\hfill {\mu }_{\mathrm{ov}}\left(k\right)={\widehat{G}}_0({(1-k)}^{{\widehat{n}}_{\mathrm{ov}}}+{\widehat{\alpha }}_{\mathrm{ov}}).\end{array}$$

(60)

The viscosity is modeled by

$$\begin{array}{cc}\hfill \widehat{\eta }(\dot{\mathbf{C}},\mathbf{C},\Theta ,{\overline{\mathbf{C}}}_{\mathrm{p}},a)& ={\widehat{\eta }}_0(\Theta ,a)\left({\left({\displaystyle \frac{1}{s_{\eta }}}{\dot{s}}_{\mathrm{M}}(\dot{\mathbf{C}},\mathbf{C},\Theta ,{\overline{\mathbf{C}}}_{\mathrm{p}},a)+{\widehat{\alpha }}_{\eta }(\Theta )\right)}^{-{\widehat{r}}_{\eta }(\Theta )}+1\right),\hfill \end{array}$$

(61)

and depends on the right Cauchy-Green tensor $\mathbf{C}$, the plastic right Cauchy-Green tensor ${\overline{\mathbf{C}}}_{\mathrm{p}}$, the aging a, and the temperature $\Theta$ considering the dependencies of the variable ${\dot{s}}_{\mathrm{M}}(\dot{\mathbf{C}},\mathbf{C},\Theta ,{\overline{\mathbf{C}}}_{\mathrm{p}},a)$. The temperature and aging-dependent functions were modeled during the identification step

$$\begin{array}{cc}\hfill {\widehat{\eta }}_0(\Theta ,a)& =({\widehat{\eta }}_{01}-{\widehat{\eta }}_{02}\exp \left({\widehat{\eta }}_{03}a\right))(1+\tanh \left({\widehat{n}}_{\eta }(\Theta -{\widehat{\Theta }}_{\eta })\right),\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill {\widehat{r}}_{\eta }(\Theta )& ={\widehat{r}}_{\eta 1}(1-\tanh \left({\widehat{n}}_{\eta }(\Theta -{\widehat{\Theta }}_r)\right),\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill {\widehat{\alpha }}_{\eta }(\Theta )& ={\widehat{\alpha }}_{\eta 1}\exp ({\widehat{\alpha }}_{\eta 2}\Theta ).\hfill \end{array}$$

(64)

A summary of the entire model is presented in

Table 1. Summary of the constitutive model (the equation numbers are related to the specified equations).

Kinematic
$\mathbf{F}={\mathbf{F}}_a{\mathbf{F}}_{\Theta }{\mathbf{F}}_{\mathrm{M}}={\mathbf{F}}_a{\mathbf{F}}_{\Theta }{\mathbf{F}}_{\mathrm{k}}{\overline{\mathbf{F}}}_{\mathrm{p}}={\mathbf{F}}_a{\mathbf{F}}_{\Theta }{\mathbf{F}}_{\mathrm{e}}{\overline{\mathbf{F}}}_{\mathrm{v}}$,    $J=\det \mathbf{F}$	(2), (9), (10)
${\mathbf{F}}_{\Theta }={\phi }_{\Theta }(\Theta )\mathbf{I},{\mathbf{F}}_a={\phi }_a\left(a\right)\mathbf{I}$,       $J_{\mathrm{M}}={\left({\phi }_a\left(a\right){\phi }_{\Theta }(\Theta )\right)}^{-3}J$	(3)
2nd Piola-Kirchhoff stresses in the reference configuration
$\tilde{\mathbf{T}}={\tilde{\mathbf{T}}}_{\mathrm{eq}}+{\tilde{\mathbf{T}}}_{\mathrm{ov}}={\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{e}}+{\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{h}}+{\tilde{\mathbf{T}}}_{\mathrm{ov}}$	(1)
${\displaystyle {\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{e}}=\frac{K}{10}\left(J_{\mathrm{M}}^5-J_{\mathrm{M}}^{-5}\right){\mathbf{C}}^{-1}+2c_{10}{J}^{-2/3}{\mathbf{C}}^{-1}{\mathbf{C}}^{\mathrm{D}}}$	(A57), (A58)
${\tilde{\mathbf{T}}}_{\mathrm{eq}}^{\mathrm{h}}=2{\mu }_{\mathrm{eq}}(\Theta )J^{-2/3}{\mathbf{C}}^{-1}{\left(\mathbf{C}{\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\right)}^{\mathrm{D}}$	(A64)
${\tilde{\mathbf{T}}}_{\mathrm{ov}}=2{\mu }_{\mathrm{ov}}\left(k\right)J^{-2/3}{\mathbf{C}}^{-1}{\left(\mathbf{C}{\overline{\mathbf{C}}}_{\mathrm{v}}^{-1}\right)}^{\mathrm{D}}$	(A65)
Kirchhoff stresses in the current configuration
$\mathbf{S}={\mathbf{S}}_{\mathrm{vol}}^{\mathrm{e}}+{\mathbf{S}}_{\mathrm{iso}}^{\mathrm{e}}+{\mathbf{S}}_{\mathrm{eq}}^{\mathrm{h}}+{\mathbf{S}}_{\mathrm{ov}}$
${\displaystyle {\mathbf{S}}_{\mathrm{vol}}^{\mathrm{e}}=\frac{K}{10}\left(J_{\mathrm{M}}^5-J_{\mathrm{M}}^{-5}\right)\mathbf{I}}$,  ${\mathbf{S}}_{\mathrm{iso}}^{\mathrm{e}}=2c_{10}{\overline{\mathbf{B}}}_{\mathrm{M}}^{\mathrm{D}}$	(A59), (A60)
${\mathbf{S}}_{\mathrm{eq}}^{\mathrm{h}}=2{\mu }_{\mathrm{eq}}(\Theta ){\overline{\mathbf{B}}}_{\mathrm{k}}^{\mathrm{D}}$	(A66)
${\mathbf{S}}_{\mathrm{ov}}=2{\mu }_{\mathrm{ov}}\left(k\right){\overline{\mathbf{B}}}_{\mathrm{e}}^{\mathrm{D}}$	(A67)
Evolution equations in the reference configuration
${\dot{\mathbf{C}}}_{\mathrm{p}}=4{\left({\phi }_a\left(a\right){\phi }_{\Theta }(\Theta )\right)}^{-2}{J}^{-2/3}\widehat{b}{\dot{s}}_{\mathrm{M}}{\overline{\mathbf{C}}}_{\mathrm{p}}{\left({\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\mathbf{C}\right)}^{\mathrm{D}}$	(49)
${\displaystyle {\dot{\mathbf{C}}}_{\mathrm{v}}=\frac{4{\mu }_{\mathrm{ov}}\left(k\right)}{\widehat{\eta }(\dot{\mathbf{C}},\mathbf{C},\Theta ,{\overline{\mathbf{C}}}_{\mathrm{p}},a)}{\left({\phi }_a\left(a\right){\phi }_{\Theta }(\Theta )\right)}^{-2}{J}^{-2/3}{\overline{\mathbf{C}}}_{\mathrm{v}}{\left({\overline{\mathbf{C}}}_{\mathrm{v}}^{-1}\mathbf{C}\right)}^{\mathrm{D}}}$	(50)
${\displaystyle \dot{a}={\beta }_{a1}(\Theta )\frac{1-a}{{\beta }_{a2}+a}}$	(58)
$\dot{k}={\alpha }_k(1-k){\dot{s}}_{\mathrm{M}}$	(54)
${\dot{s}}_{\mathrm{M}}=J^{-2/3}/2\sqrt{{\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}\mathbf{C}{\left({\mathbf{C}}^{-1}\dot{\mathbf{C}}\right)}^{\mathrm{D}}{\overline{\mathbf{C}}}_{\mathrm{p}}^{-1}·\mathbf{C}{\left({\mathbf{C}}^{-1}\dot{\mathbf{C}}\right)}^{\mathrm{D}}}$	(53)
Temperature, softening and aging-dependent functions
${\beta }_{a1}(\Theta )={\beta }_{a11}\exp (-{\beta }_{a12}/\Theta )$	(58)
${\mu }_{\mathrm{eq}}(\Theta )={\widehat{c}}_1({\widehat{c}}_2-\tanh \left({\widehat{c}}_3(\Theta -{\widehat{\Theta }}_c)\right))$	(59)
${\mu }_{\mathrm{ov}}\left(k\right)={\widehat{G}}_0({(1-k)}^{{\widehat{n}}_{\mathrm{ov}}}+{\widehat{\alpha }}_{\mathrm{ov}})$	(60)
$\widehat{\eta }(\dot{\mathbf{C}},\mathbf{C},\Theta ,{\overline{\mathbf{C}}}_{\mathrm{p}},a)={\widehat{\eta }}_0(\Theta ,a)\left({\displaystyle {\left(\frac{{\dot{s}}_{\mathrm{M}}}{s_{\eta }}+{\widehat{\alpha }}_{\eta }(\Theta )\right)}^{-{\widehat{r}}_{\eta }(\Theta )}+1}\right)$	(61)
${\widehat{\eta }}_0(\Theta ,a)=({\widehat{\eta }}_{01}-{\widehat{\eta }}_{02}\exp \left({\widehat{\eta }}_{03}a\right))(1+\tanh \left({\widehat{n}}_{\eta }(\Theta -{\widehat{\Theta }}_{\eta })\right)$	(62)
${\widehat{r}}_{\eta }(\Theta )={\widehat{r}}_{\eta 1}(1-\tanh \left({\widehat{n}}_{\eta }(\Theta -{\widehat{\Theta }}_r)\right)$,  ${\widehat{\alpha }}_{\eta }(\Theta )={\widehat{\alpha }}_{\eta 1}\exp ({\widehat{\alpha }}_{\eta 2}\Theta )$	(63), (64)
${\phi }_{\Theta }(\Theta )=1+{\alpha }_{\Theta }(\Theta -{\Theta }_0)$,    ${\phi }_a\left(a\right)=1-{\alpha }_a$	(4), (5)

.

3.3. Modeling of the Thermo-Physical Properties

For the computation of coupled thermo-mechanical problems, expressions for the specific heat capacity $c_{\mathrm{p}}$ and thermal conductivity ${\mathit{\kappa }}_{\Theta }$ are necessary. For the development of these functions for Zamak 5, the measurements from [33] are taken into account. The linear thermal expansion ${\alpha }_{\Theta }$ was already defined in Equation (4), with the selection of the thermal part of the deformation gradient ${\mathbf{F}}_{\Theta }$.

3.3.1. Specific Heat Capacity

From the considerations in the Clausius-Duhem inequality, see Appendix C, there is a relationship between the specific heat capacity and the specific free energy. Equation (28) results in a linear behavior over the temperature for the specific heat $c_{\mathrm{p}}$

$$\begin{array}{c}\hfill c_{\mathrm{p}}(\Theta )=c_{\mathrm{p}1}(1+c_{\mathrm{p}2}(\Theta -{\Theta }_0)).\end{array}$$

(65)

Here, $c_{\mathrm{p}1}$ and $c_{\mathrm{p}2}$ are material parameters, and ${\Theta }_0$ is the reference temperature.

3.3.2. Thermal Diffusivity

Considering the measurements of the thermal diffusivity in [33], it can be seen that the thermal diffusivity depends on the temperature and the aging variable. In order to model this behavior, a weighting between the initial unaged stage ${\widehat{\lambda }}_0(\Theta )$ and the final stage ${\widehat{\lambda }}_{\infty }(\Theta )$ is carried out with the help of a weighting function $f_{\widehat{\lambda }}(\Theta ,a)$

$$\begin{array}{c}\hfill \widehat{\lambda }(\Theta ,a)={\widehat{\lambda }}_{\infty }(\Theta )f_{\widehat{\lambda }}(\Theta ,a)+{\widehat{\lambda }}_0(\Theta )(1-f_{\widehat{\lambda }}(\Theta ,a)).\end{array}$$

(66)

A linear behavior over the temperature is observed for the completely aged state, see Figure 15,

$$\begin{array}{c}\hfill {\widehat{\lambda }}_{\infty }(\Theta )={\widehat{\lambda }}_{\infty 1}-{\widehat{\lambda }}_{\infty 2}\Theta .\end{array}$$

(67)

As one can see in Figure 15, the initial state shows a constant behavior up to the temperature of ${\Theta }_0=293.15 \mathrm{K}$. After that, an approximately linear reduction is observed. In order to model this behavior, the following function is chosen

$$\begin{array}{c}\hfill {\widehat{\lambda }}_0(\Theta )={\widehat{\lambda }}_{01}-{\widehat{\lambda }}_{02}(1+\tanh (\Theta -{\Theta }_0))(\Theta -{\Theta }_0).\end{array}$$

(68)

For the weighting factor $f_{\widehat{\lambda }}$ from Equation (66) we chose a function of the aging and temperature, which has a growing behavior with the aging. It has values between 0 and 1, where 0 corresponds to the unaged state and 1 to the completely aged state

$$\begin{array}{c}\hfill f_{\widehat{\lambda }}(\Theta ,a)=a-{\alpha }_{\widehat{\lambda }}(\Theta )(1-a)a^2,{\alpha }_{\widehat{\lambda }}(\Theta )={\alpha }_{\widehat{\lambda }1}\exp (-{\alpha }_{\widehat{\lambda }2}\Theta ).\end{array}$$

(69)

This expression was chosen during the identification.

3.3.3. Thermal Conductivity

The thermal conductivity can be obtained from the relation

$$\begin{array}{c}\hfill {\mathit{\kappa }}_{\Theta }(\Theta ,a)=\rho \widehat{\lambda }(\Theta ,a)c_{\mathrm{p}}(\Theta ).\end{array}$$

(70)

With the proposed expressions for the thermal diffusivity and the specific heat capacity, the thermal conductivity is equal to

$$\begin{array}{c}\hfill {\mathit{\kappa }}_{\Theta }(\Theta ,a)=\rho c_{\mathrm{p}1}(1+c_{\mathrm{p}2}(\Theta -{\Theta }_0))({\widehat{\lambda }}_{\infty }(\Theta )f_{\widehat{\lambda }}(\Theta ,a)+{\widehat{\lambda }}_0(\Theta )(1-f_{\widehat{\lambda }}(\Theta ,a))),\end{array}$$

(71)

with the functions ${\widehat{\lambda }}_{\infty }$, ${\widehat{\lambda }}_0$ and $f_{\widehat{\lambda }}$ from Equations (67)–(69).

4. Parameter Identification

The parameter identification in performed in several steps, where each phase considers a different phenomenon. First, the aging variable is considered using the shrinkage measurements. Second, the identification of the mechanical part of the model is performed in a partitioned way. The parameters of the equilibrium stress state are identified first using the termination points of the relaxation. Once these parameters are known, the parameters of the overstress part are identified with the help of the experiments at different strain rates. In the final step, the thermo-physical properties are identified.

4.1. Identification of the Aging Variable

The material parameters related to the aging variable can be identified with the measurements of the shrinkage considering the evolution Equation (58) and the relation of the stretch due to aging ${\lambda }_a={\phi }_a=1-{\alpha }_{a}a=L/L_0$. Here, the parameter ${\alpha }_a$ corresponds to the maximum shrinkage reached by the material. The parameters responsible for the development of the aging over time are ${\beta }_{a11}$, ${\beta }_{a12}$ and ${\beta }_{a2}$, see Equation (58). In a first identification step, the parameters ${\alpha }_a$, ${\beta }_{a2}$ and ${\beta }_{a1}^i$, are identified, with ${\beta }_{a1}^i$, $i=1,2,3$, changing for every temperature ${\Theta }_1=85 °\mathrm{C}$, ${\Theta }_2=105 °\mathrm{C}$, ${\Theta }_3=125 °\mathrm{C}$. The results of this first identification are used in a subsequent second identification step. Here, the parameters of the Arrhenius equation, ${\beta }_{a11}$ and ${\beta }_{a12}$, are identified. To this end, we consider the natural logarithm of Equation (57)

$$\begin{array}{c}\hfill \ln {\beta }_{a1}(\Theta )=\ln \left({\beta }_{a11}\right)-{\beta }_{a12}{\displaystyle \frac{1}{\Theta }}.\end{array}$$

(72)

Equation (72) corresponds to a linear equation $y\left(x\right)=y_0-mx$, with $y\left(x\right)=\ln {\beta }_{a1}$ in dependence of $x=1/\Theta$. The parameters $y_0=\ln \left({\beta }_{a11}\right)$ and $m={\beta }_{a12}$ were found with the help of the points ${\beta }_{a1}^i$ obtained in the previous identification step. The results of the second identification are shown in Figure 8 and

Table 2. Parameter identification of the Arrhenius equation.

Parameter	${\mathit{\alpha }}_{\mathit{a}}$	$\ln \left({\mathit{\beta }}_{\mathit{a}11}\right)$	${\mathit{\beta }}_{\mathit{a}12}$
	−	−	$\mathrm{K}$
initial value	10−3	10−3	1
final value	1.37 × 10−3	$10.24$	8113
conf. interval	10−7	$0.787$	2748

. In Figure 8b one can see the prediction of the model (brown line) at the temperature of $20 °\mathrm{C}$ compared to the measurement from [8] (points). The black line represents the response of the model at $20 °\mathrm{C}$ with the identified parameters but with a change in the saturation value ${\alpha }_a$, where the value identified in [4] is used, since it is clear that the saturation value in the measurement shown in [8] is different. Here, one can see that—although their saturation value ${\alpha }_a$ is different—the process for the evolution of the aging variable represented by the parameters ${\beta }_{a2}$, ${\beta }_{a11}$, and ${\beta }_{a12}$ corresponds to the measurements of [8].

In Figure 8, it is shown that the model can reproduce the main behavior of the experiments. The small values of the function ${\beta }_{a1}$ for temperatures below $80 °\mathrm{C}$ indicate that the aging process happens quite slowly.

4.2. Identification of the Mechanical Response

Contrarily to the small deformations model from [3,4], the case of torsion is not purely deviatoric for finite deformation, since there is a reaction force in longitudinal direction. Thus, the uniaxial tension/compression experiments are used in the identification with the help of the procedure exposed in [69]. The modular structure of the model is drawn on for material parameter identification. First, the parameters of the equilibrium stress state are identified with the equilibrium hysteresis and after that, the overstress with the experiments at different strain rates and multi-step relaxation tests are determined.

4.2.1. Elastic Part of the Equilibrium Stress State

The two parameters K and $c_{10}$ of the elastic equilibrium stress state, see Equations (A59) and (A60), are not influenced by aging or temperature and are taken from the small deformations case in [3] to prevent a correlation between the parameters during the identification process

$$\begin{array}{c}\hfill K=K_{\mathrm{eq}}=49.5 \mathrm{GPa},c_{10}=G_{\mathrm{eq}}/2=0.05 \mathrm{GPa}.\end{array}$$

(73)

The parameter $c_{10}$ corresponds to the shear modulus of the elastic equilibrium stress state and it is proportional to $G_{\mathrm{eq}}$, see [3].

4.2.2. Hysteretic Part of the Equilibrium Stress State

The material parameters of the hysteretic part of the equilibrium stress state $\mathit{\kappa }=\{\widehat{b},{\widehat{c}}_1,{\widehat{c}}_2,{\widehat{c}}_3,{\widehat{\Theta }}_c\}$ from Equation (59), are identified at the termination points of relaxation of multi-step relaxation processes in tension and compression at different temperatures.

Table 3. Identified parameters of the equilibrium stress part in tension/compression.

Parameter	$\widehat{\mathit{b}}$	${\widehat{\mathit{c}}}_1$	${\widehat{\mathit{c}}}_2$	${\widehat{\mathit{c}}}_3$	${\widehat{\mathsf{\Theta }}}_{\mathit{c}}$
	−	MPa	MPa	−	K
initial values	$54.77$	$5793.5$	$1.1$	$0.24$	$283.44$
final values	$58.94$	$5924.8$	$1.12$	$0.0279$	$288.03$
conf. intervals	$4.14$	$541.5$	$0.02$	$0.0019$	$3.77$

and Figure 9 display the results of the identification. Here, the axial component of 1st Piola-Kirchhoff tensor (axial force per cross section in the undeformed state) is represented by ${\sigma }_R$ and the axial stretch by $\lambda$. The initial values for the non-linear optimizer are given by the model for small strains. The confidence intervals are sufficiently small for all material parameters.

4.2.3. Overstress Part

The identification of the overstress part (A67), i.e., the shear modulus (60), is performed in two steps. The parameters of the shear modulus ${\mu }_{\mathrm{ov}}$ and the parameters of the viscosity $\widehat{\eta }$ in Equation (61) are identified first for a constant aging time considering the experiments without aging. The feature of interest here is temperature dependence. In this identification step there are 40 experimental curves in total: four different strain rates plus the multi-step relaxation test at four different temperatures and for tension and compression. In the second step, the parameters considering aging dependence are found with the experiments at different aging times. In this second identification, the temperature-dependent parameters are known.

1st Step: Shear Modulus

There are three parameters that have to be identified for the identification of the shear modulus (60), namely $\mathit{\kappa }=\{{\widehat{G}}_0,{\widehat{\alpha }}_{\mathrm{ov}},{\widehat{n}}_{\mathrm{ov}}\}$. The parameter of the softening is set to ${\alpha }_k=350$, so that it reaches the value 1 at the end of the loading process for a strain equal to 1%. During this identification step, the value of the viscosity is identified independently for each strain rate and temperature for the material without aging.The results of the identification are shown in

Table 4. Results of the parameter identification for ${\mu }_{\mathrm{ov}}$, see Equation (60).

Parameter	${\widehat{\mathit{G}}}_0$	${\widehat{\mathit{\alpha }}}_{\mathbf{ov}}$	${\widehat{\mathit{n}}}_{\mathbf{ov}}$
	MPa	−	−
initial value	1.604 × 104	2 × 10−1	3.5 × 10−1
final value	1.424 × 104	2.531 × 10−1	3.349 × 10−1
conf. interval	$15.3$	5.71 × 10−4	6.8 × 10−4

.

2nd Step: Viscosity

The parameters of the shear modulus are assumed to be fixed for the subsequent identification step, see Table 4. The material parameters $\mathit{\kappa }=\{{\widehat{\Theta }}_{\eta },{\widehat{r}}_{\eta 1},{\widehat{\Theta }}_r,{\widehat{n}}_{\eta },{\widehat{\alpha }}_{\eta 1},{\widehat{\alpha }}_{\eta 2}\}$ in Equations (61) and (62) have to be identified next. They are related to the temperature and rate dependence. A value for the function ${\widehat{\eta }}_0(a=0)={\widehat{\eta }}_{01}-{\widehat{\eta }}_{02}$ of Equation (62) is identified, which is used later on in the identification of the aging dependence. The results of the identification are shown in

Table 5. Results of the parameter identification for the temperature-dependence of the viscosity.

Parameter	${\widehat{\mathsf{\Theta }}}_{\mathit{\eta }}$	${\widehat{\mathit{r}}}_{\mathit{\eta }1}$	${\widehat{\mathsf{\Theta }}}_{\mathit{r}}$	${\widehat{\mathit{n}}}_{\mathit{\eta }}$	${\widehat{\mathit{\alpha }}}_{\mathit{\eta }1}$	${\widehat{\mathit{\alpha }}}_{\mathit{\eta }2}$
	K	−	K	K−1	−	K−1
initial value	$373.6$	$0.480$	$426.3$	$0.0115$	1.26 × 10−9	$0.0135$
final value	$375.2$	$0.479$	$428.2$	$0.0117$	3.47 × 10−10	$0.0173$
conf. interval	$0.837$	2.8 × 10−4	$0.192$	9.7 × 10−6	4.80 × 10−11	4.1 × 10−4

showing that the confidence intervals are less than one order, and appear to be sufficiently accurate considering the complexity of the model.

The correlation matrix is equal to

$$\begin{array}{c}\hfill \mathsf{R}=\left[\begin{array}{ccccccc}& {\widehat{\Theta }}_{\eta }& {\widehat{r}}_{\eta 1}& {\widehat{\Theta }}_r& {\widehat{n}}_{\eta }& {\widehat{\alpha }}_{\eta 1}& {\widehat{\alpha }}_{\eta 2}\\ {\widehat{\Theta }}_{\eta }& 1& 0.916& -0.876& -0.645& 0.394& -0.391\\ {\widehat{r}}_{\eta 1}& & 1& -0.826& -0.539& 0.405& -0.386\\ {\widehat{\Theta }}_r& & & 1& 0.222& -0.452& 0.452\\ {\widehat{n}}_{\eta }& & & & 1& -0.136& 0.134\\ {\widehat{\alpha }}_{\eta 1}& & & & & 1& -0.997\\ {\widehat{\alpha }}_{\eta 2}& & & & & & 1\end{array}\right]\end{array}$$

(74)

showing a strong correlation between the parameters ${\widehat{\Theta }}_{\eta }$, ${\widehat{r}}_{\eta 1}$, and ${\widehat{\Theta }}_r$, which consider the temperature dependence of the functions ${\widehat{\eta }}_0$ and ${\widehat{r}}_{\eta }$, and also between the parameters ${\widehat{\alpha }}_{\eta 1}$ and ${\widehat{\alpha }}_{\eta 2}$ from the function ${\widehat{\alpha }}_{\eta }(\Theta )$. The latter is known for exponential functions, see [70].

The parameters of the aging dependence are identified in the last identification step. Here, the identification for the value ${\widehat{\eta }}_0\left(a\right)$ is performed separately for each aging time. After that, the function of the aging dependence ${\widehat{\eta }}_0\left(a\right)={\widehat{\eta }}_{01}-{\widehat{\eta }}_{02}\exp \left({\widehat{\eta }}_{03}a\right)$ is fitted with the identified points using Matlab.

Table 6. Results of the parameter identification for the viscosity function ${\widehat{\eta }}_0\left(a\right)$ of Equation (62).

Parameter	${\widehat{\mathit{\eta }}}_{01}$	${\widehat{\mathit{\eta }}}_{02}$	${\widehat{\mathit{\eta }}}_{03}$
	MPa s	MPa s	−
initial value	$0.265$	$0.068$	$0.305$
final value	1612	$81.93$	$1.862$
conf. interval	$285.75$	$9.6$	$0.78$

and Figure 10 show the results of this identification step.

The complete results of the identification of the temperature and the aging dependence are shown in Figure 11 and Figure 12.

One can observe a very good agreement between the experiments and the fit for the tensile tests, and a good agreement for compression. The difference originates in the tension/compression asymmetry, which is not considered in the model. In

Table 7. Summary of the identification.

Equilibrium Stress
parameter	K	$c_{10}$	$\widehat{b}$	${\widehat{c}}_1$	${\widehat{c}}_2$	${\widehat{c}}_3$	${\widehat{\Theta }}_c$
value	4.95 × 104	50	$58.94$	5.925 × 103	$1.12$	$0.028$	$288.03$
dimension	$\mathrm{MPa}$	$\mathrm{MPa}$	−	$\mathrm{MPa}$	$\mathrm{MPa}$	${\mathrm{K}}^{-1}$	$\mathrm{K}$
Shear Modulus Overstress
parameter		${\widehat{\alpha }}_k$		${\widehat{G}}_0$		${\widehat{\alpha }}_{\mathrm{ov}}$		${\widehat{n}}_{\mathrm{ov}}$
value		350		1.424 × 104		$0.253$		$0.335$
dimension		−		$\mathrm{MPa}$		−		−
Viscosity
parameter	${\widehat{\eta }}_{01}$	${\widehat{\eta }}_{02}$	${\widehat{\eta }}_{03}$	${\widehat{\Theta }}_{\eta }$	${\widehat{r}}_{\eta 1}$	${\widehat{\Theta }}_r$	${\widehat{n}}_{\eta }$	${\widehat{\alpha }}_{\eta 1}$	${\widehat{\alpha }}_{\eta 2}$
value	1612	$81.93$	$1.862$	$375.2$	$0.479$	$428.2$	$0.012$	3.47 × 10−10	$0.0173$
dimension	$\mathrm{MPa}\mathrm{s}$	$\mathrm{M}\mathrm{Pa}\mathrm{s}$	−	$\mathrm{K}$	−	$\mathrm{K}$	${\mathrm{K}}^{-1}$	−	${\mathrm{K}}^{-1}$

, a summary of the identified parameters is given.

4.3. Identification of the Thermo-Physical Properties

4.3.1. Thermal Expansion

The thermal expansion coefficient was identified in [3] equal to ${\alpha }_{\Theta }=2.8109 \times {10}^{-5} {\mathrm{K}}^{-1}$, with the reference temperature ${\Theta }_0=20 °\mathrm{C}/293.15 \mathrm{K}$, see Figure 13a.

4.3.2. Specific Heat Capacity

The parameters of the specific heat capacity of Equation (65) are identified with the measurements of [33]. The identified parameters are $c_{\mathrm{p}1}=413.724 \mathrm{J}\mathrm{k}{\mathrm{g}}^{-1}{\mathrm{K}}^{-1}$ and $c_{\mathrm{p}2}=5.581 \times {10}^{-4} {\mathrm{K}}^{-1}$, see Figure 13b.

4.3.3. Thermal Diffusivity

The identification of the thermal diffusivity is connected to the modeling and is carried out in two steps. In Equation (66), the thermal diffusivity is expressed as a weighting between the initial unaged stage ${\widehat{\lambda }}_0(\Theta )$ and the final stage ${\widehat{\lambda }}_{\infty }(\Theta )$ by means of a weighting factor $f_{\widehat{\lambda }}(\Theta ,a)$

The first step of the identification considers the initial and final state of the thermal diffusivity, for which the weighting factor is known (for the initial state $f_{\widehat{\lambda }}(\Theta ,a=0)=0$ and for the final state $f_{\widehat{\lambda }}(\Theta ,a=1)=1$). In the second step the weighting factor $f_{\widehat{\lambda }}(\Theta ,a)$ is unknown and has to be determined.

Thus, the functions ${\widehat{\lambda }}_0(\Theta )$ and ${\widehat{\lambda }}_{\infty }(\Theta )$ are considered first, which have a total number of four parameters to be identified ${\mathit{\kappa }}^{\mathrm{T}}=\{{\widehat{\lambda }}_{\infty 1},{\widehat{\lambda }}_{\infty 2},{\widehat{\lambda }}_{01},{\widehat{\lambda }}_{02}\}$. The parameter ${\Theta }_0$ corresponds to the reference temperature and is set to ${\Theta }_0=293.15 \mathrm{K}$.

The identification of these parameter is carried out with the measurements over the temperature for the unaged specimen ($f_{\widehat{\lambda }}(\Theta ,a=0)=0$) and the completely aged specimen ($f_{\widehat{\lambda }}(\Theta ,a=1)=1$).

Table 8. Results of the parameter identification of the functions ${\widehat{\lambda }}_0(\Theta )$ and ${\widehat{\lambda }}_{\infty }(\Theta )$.

Parameter	${\widehat{\mathit{\lambda }}}_{\mathit{\infty }1}$	${\widehat{\mathit{\lambda }}}_{\mathit{\infty }2}$	${\widehat{\mathit{\lambda }}}_{01}$	${\widehat{\mathit{\lambda }}}_{02}$
	mm2 s−1	mm2 s−1 K−1	mm2 s−1	mm2 s−1 K−1
initial value	1	1	42	5 × 10−3
final value	$55.07$	2.92 × 10−2	$41.74$	4.63 × 10−3
conf. interval	$0.345$	0.12 × 10−2	$0.074$	1.2 × 10−3

and Figure 15 display the results of this identification. In Figure 15, the results correspond to the initial and final state (black and blue lines).

In the second identification step, a function for the weighting factor $f_{\widehat{\lambda }}(\Theta ,a)$ is developed. To this end, an experimental value for $f_{\widehat{\lambda }\exp }$ is calculated with the help of Equation (66) using the experimental points by

$$\begin{array}{c}\hfill f_{\widehat{\lambda }\exp }({\Theta }_i,a_j)={\displaystyle \frac{{\widehat{\lambda }}_{\exp }({\Theta }_i,a_j)-{\widehat{\lambda }}_{\exp }({\Theta }_i,0)}{{\widehat{\lambda }}_{\exp }({\Theta }_i,1)-{\widehat{\lambda }}_{\exp }({\Theta }_i,0)}},\end{array}$$

(75)

where ${\widehat{\lambda }}_{\exp }({\Theta }_i,a_j)$ represents the measured point at the temperature ${\Theta }_i$ and the aging time $a_j$, ${\widehat{\lambda }}_{\exp }({\Theta }_i,0)$ denotes the initial state ($a=0$), and ${\widehat{\lambda }}_{\exp }({\Theta }_i,1)$ is the final state ($a=1$). The function $f_{\widehat{\lambda }\exp }$ is determined for the temperatures from $-60 °\mathrm{C}$ to $20 °\mathrm{C}$, since one can assume that the aging does not evolve during the time of the experiments for these temperatures. These experimental points are used to develop the function for the weighting factor.

In the identification, the factor is determined independently for each temperature $f_{\widehat{\lambda }}(\widehat{\Theta }=\mathrm{const}.,a)$. Afterwards, the function ${\alpha }_{\widehat{\lambda }}(\Theta )$ is fitted to the results of this pre-identification, obtaining the values ${\alpha }_{\widehat{\lambda }1}=7.093$ (conf. int. 1) and ${\alpha }_{\widehat{\lambda }2}=4.937 \times {10}^{-3} {\mathrm{K}}^{-1}$ (conf. int. $2.301 \times {10}^{-4} {\mathrm{K}}^{-1}$). The results are shown in Figure 14.

The results of the entire identification of the thermal diffusivity are shown in Figure 15. One can see the small deviation between the simulation and the experiment for the 3 and 6 months old samples. However, the modeled function can reproduce the complex experimental behavior well.

4.3.4. Thermal Conductivity

With the identified values of the thermal diffusivity, the specific heat capacity, and the measured density, the parameter of the thermal conductivity can be determined with Equation (70). The thermal conductivity is shown in Figure 16, dependent on the temperature for different aging times.

5. Stress Computation and Numerical Example

The constitutive model described in Table 1 along with the non-linear kinematics, the balance of linear momentum, the heat equation in connection with the initial and boundary conditions form a coupled system of partial and ordinary differential equations. A common approach is to apply the method of vertical lines, see [71,72]. In a first step, the spatial discretization is carried out yielding a system of differential-algebraic equations. The algebraic part stems from the weak form of the balance of linear momentum, and the ordinary differential parts is result of the weak form of the heat equation and the evolution equations of the constitutive model, see [73,74] for the concept in thermo-mechanical coupled situations, and for treating the non-linear equation [75]. In the subsequent considerations only low order finite elements and a time discretization using the Backward Euler method are applied. This leads in each point in time to a system of non-linear equations

$$\begin{array}{cc}\hfill {\mathsf{G}}_u(\mathsf{u},\Theta ,\mathsf{q})& :={\mathsf{g}}_u(t_{n+1},{\mathsf{u}}_{n+1},{\Theta }_{n+1},{\mathsf{q}}_{n+1})=\mathsf{0},\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill {\mathsf{G}}_{\Theta }(\mathsf{u},\Theta ,\mathsf{q})& :={{\mathsf{C}}_{\mathsf{p}}}_{n+1}{\displaystyle \frac{{\Theta }_{n+1}-{\Theta }_n}{\mathsf{\Delta }t}}-{\mathsf{r}}_{\Theta }(t_{n+1},{\mathsf{u}}_{n+1},{\Theta }_{n+1},{\mathsf{q}}_{n+1})=\mathsf{0},\hfill \end{array}$$

(77)

$$\begin{array}{cc}\hfill {\mathsf{L}}_q(\mathsf{u},\Theta ,\mathsf{q})& :={\displaystyle \frac{{\mathsf{q}}_{n+1}-{\mathsf{q}}_n}{\mathsf{\Delta }t}}-{\mathsf{r}}_q(t_{n+1},{\mathsf{u}}_{n+1},{\Theta }_{n+1},{\mathsf{q}}_{n+1})=\mathsf{0}.\hfill \end{array}$$

(78)

The solution of the non-linear system has to be done at each point in time $t_{n+1}$ to obtain the unknown nodal displacements ${\mathsf{u}}_{n+1}$, the unknown nodal temperatures ${\Theta }_{n+1}$, and internal variables ${\mathsf{q}}_{n+1}$, which are (formally) assembled from all Gauss-points of all elements into a large vector ${\mathsf{q}}_{n+1}$. Furthermore, we have to consider the initial conditions $\mathsf{u}(t=0)={\mathsf{u}}_0$, $\Theta (t=0)={\Theta }_0$, and $\mathsf{q}(t=0)=\mathsf{0}$. Here, ${\mathsf{G}}_u$ represents the discretized form of the balance of linear momentum, ${\mathsf{G}}_{\Theta }$ the integration step of discretized weak form of the heat equation, and ${\mathsf{L}}_q$ the integration step of the discretized evolution equations of the internal variables, see Table 1. The computation of the Cauchy stress $\mathsf{T}{\in \mathbb{R}}^6$ (here, the Voigt notation is chosen instead of the tensor notation), the internal variables $\mathsf{q}{\in \mathbb{R}}^{6+6+3}$, and the consistent tangent ${\mathsf{C}}_t$ required for the global Multilevel-Newton algorithm—see, for details, Refs. [74,76] and the literature cited therein—are computed within the stress algorithm. Afterwards, the stress algorithm of the model is provided and the applicability of the model is demonstrated with two different validation examples.

5.1. Stress Algorithm

In the proposed model, the vector of internal variables ${\mathsf{q}}^{\mathrm{T}}=\{{\mathsf{C}}_{\mathrm{p}}^{\mathrm{T}},{\mathsf{C}}_{\mathrm{v}}^{\mathrm{T}},k,a\}$ has two vectorial components, ${\mathsf{C}}_{\mathrm{p}}{\in \mathbb{R}}^6$ and ${\mathsf{C}}_{\mathrm{v}}{\in \mathbb{R}}^6$, which are expressed in Voigt notation, and two scalar components, k and a, as well as the arc-length $s_{\mathrm{M}}$. By applying a Backward Euler step to the evolution equations of the internal variables at each Gauss point, see Table 1, one obtains

$$\begin{array}{cc}\hfill \mathsf{0}& ={{\mathsf{C}}_{\mathrm{p}}}_{n+1}-{{\mathsf{C}}_{\mathrm{p}}}_n-\mathsf{\Delta }{t}_n\left(4{\left({\phi }_a{\phi }_{\Theta }\right)}^{-2}{J}^{-2/3}b{\displaystyle \frac{{\mathsf{\Delta }{s}_{\mathrm{M}}}_n}{\mathsf{\Delta }{t}_n}}{{\mathsf{C}}_{\mathrm{p}}}_{n+1}\mathsf{D}\left({{\mathsf{C}}_{\mathrm{p}}}_{n+1}^{-1}\mathsf{C}\right)\right),\hfill \end{array}$$

(79)

$$\begin{array}{cc}\hfill \mathsf{0}& ={{\mathsf{C}}_{\mathrm{v}}}_{n+1}-{{\mathsf{C}}_{\mathrm{v}}}_n-\mathsf{\Delta }{t}_n\left({\displaystyle \frac{4{\mu }_{\mathrm{ov}}}{\eta }}{\left({\phi }_a{\phi }_{\Theta }\right)}^{-2}{J}^{-2/3}{{\mathsf{C}}_{\mathrm{v}}}_{n+1}\mathsf{D}\left({{\mathsf{C}}_{\mathrm{v}}}_{n+1}^{-1}\mathsf{C}\right)\right),\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill 0& =k_{n+1}-k_n-\mathsf{\Delta }{t}_n\left({\alpha }_k(1-k_{n+1}){\displaystyle \frac{{\mathsf{\Delta }{s}_{\mathrm{M}}}_n}{\mathsf{\Delta }{t}_n}}\right),\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill 0& =a_{n+1}-a_n-\mathsf{\Delta }{t}_n{\beta }_{a1}\left({\Theta }_{n+1}\right){\displaystyle \frac{1-a_{n+1}}{{\beta }_{a2}+a_{n+1}}},\hfill \end{array}$$

(82)

with the increment of the mechanical arc-length

$$\begin{array}{c}\hfill {\displaystyle {\mathsf{\Delta }{s}_{\mathrm{M}}}_n=s_{n+1}-s_n=\frac{J^{-2/3}}{2}\sqrt{{\left({{\mathsf{C}}_{\mathrm{p}}}_{n+1}^{-1}\mathsf{C}\mathsf{D}\left({\mathsf{C}}^{-1}\mathsf{\Delta }\mathsf{C}\right){{\mathsf{C}}_{\mathrm{p}}}_{n+1}^{-1}\right)}^{\mathrm{T}}{\mathsf{M}}^{-1}\left(\mathsf{C}\mathsf{D}\left({\mathsf{C}}^{-1}\mathsf{\Delta }\mathsf{C}\right)\right)}.}\end{array}$$

(83)

$\mathsf{C}:={\mathsf{C}}_{n+1}{\in \mathbb{R}}^6$ defines the vector containing the components of the right Cauchy-Green tensor in the current iteration of the Multilevel-Newton algorithm, and $\mathsf{D}{\in \mathbb{R}}^{6}6$ represents the deviator operator (in matrix form). In other words, at each Gauss-point the non-linear system Equations (79)–(82) in the form of a system of non-linear equations

$$\begin{array}{cc}\hfill {\mathsf{L}}_{\mathrm{p}}(\mathsf{C},\Theta ,{\mathsf{C}}_{\mathrm{p}},a)& =\mathsf{0},\hfill \end{array}$$

(84)

$$\begin{array}{cc}\hfill {\mathsf{L}}_{\mathrm{v}}(\mathsf{C},\Theta ,{\mathsf{C}}_{\mathrm{p}},{\mathsf{C}}_{\mathrm{v}},k,a)& =\mathsf{0},\hfill \end{array}$$

(85)

$$\begin{array}{cc}\hfill {\mathsf{L}}_k(\mathsf{C},\Theta ,{\mathsf{C}}_{\mathrm{p}},k)& =0,\hfill \end{array}$$

(86)

$$\begin{array}{cc}\hfill {\mathsf{L}}_a(\Theta ,a)& =0\hfill \end{array}$$

(87)

has to be computed. After this step, the second Piola-Kirchhoff stress can be computed for the internal variables ${\mathsf{q}}^{\mathrm{T}}=\{{{\mathsf{C}}_{\mathrm{p}}}_{n+1}^{\mathrm{T}},{{\mathsf{C}}_{\mathrm{v}}}_{n+1}^{\mathrm{T}},k_{n+1},a_{n+1}\}$,

$$\begin{array}{cc}\hfill {\tilde{\mathsf{T}}}_{n+1}& =\tilde{\mathsf{h}}({\mathsf{C}}_{n+1},{\Theta }_{n+1},{{\mathsf{C}}_{\mathrm{p}}}_{n+1},{{\mathsf{C}}_{\mathrm{v}}}_{n+1},k_{n+1},a_{n+1})\hfill \\ \hfill & ={{\tilde{\mathsf{T}}}_{\mathrm{eq}}^{\mathrm{e}}}_{n+1}({\mathsf{C}}_{n+1},{\Theta }_{n+1},a_{n+1})+{{\tilde{\mathsf{T}}}_{\mathrm{eq}}^{\mathrm{h}}}_{n+1}({\mathsf{C}}_{n+1},{\Theta }_{n+1},{{\mathsf{C}}_{\mathrm{p}}}_{n+1})\hfill \\ \hfill & +{{\tilde{\mathsf{T}}}_{\mathrm{ov}}}_{n+1}({\mathsf{C}}_{n+1},{\Theta }_{n+1},{{\mathsf{C}}_{\mathrm{v}}}_{n+1}).\hfill \end{array}$$

(88)

Here, the elastic equilibrium stress state is given by

$$\begin{array}{c}\hfill {\displaystyle {{\tilde{\mathsf{T}}}_{\mathrm{eq}}^{\mathrm{e}}}_{n+1}=\frac{K}{10}\left(J_{\mathrm{M}}^5-J_{\mathrm{M}}^{-5}\right){\mathsf{C}}^{-1}+2c_{10}{J}^{-2/3}{\mathsf{C}}^{-1}\mathsf{D}\mathsf{C},}\end{array}$$

(89)

with $J_{\mathrm{M}}={\left((1-{\alpha }_{a}a)(1+{\alpha }_{\Theta }(\Theta -{\Theta }_0))\right)}^{-3}J$, and the current temperature $\Theta ={\Theta }_{n+1}$. The plastic equilibrium stress state is given by

$$\begin{array}{c}\hfill {{\tilde{\mathsf{T}}}_{\mathrm{eq}}^{\mathrm{h}}}_{n+1}=2{\mu }_{\mathrm{eq}}\left({\Theta }_{n+1}\right)J^{-2/3}{\mathsf{C}}^{-1}\mathsf{D}\left(\mathsf{C}{{\mathsf{C}}_{\mathrm{p}}}_{n+1}^{-1}\right),\end{array}$$

(90)

and the overstress reads

$$\begin{array}{c}\hfill {{\tilde{\mathsf{T}}}_{\mathrm{ov}}}_{n+1}=2{\mu }_{\mathrm{ov}}\left(k_{n+1}\right)J^{-2/3}{\mathsf{C}}^{-1}\mathsf{D}\left(\mathsf{C}{{\mathsf{C}}_{\mathrm{v}}}_{n+1}^{-1}\right).\end{array}$$

(91)

The Multilevel-Newton algorithm, see [76,77], is applied to solve the non-linear system of equations Equations (76)–(78) for which a consistent linearization is required. To this end, the consistent tangent has to be determined

$$\begin{array}{c}\hfill {\mathsf{C}}_t=\left[{\displaystyle \frac{\mathrm{d}\tilde{\mathsf{h}}}{\mathrm{d}\mathsf{C}}}{\displaystyle \frac{\mathrm{d}\tilde{\mathsf{h}}}{\mathrm{d}\Theta }}\right],\end{array}$$

(92)

with

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}\tilde{\mathsf{h}}}{\mathrm{d}\mathsf{C}}}& ={\displaystyle \frac{\partial \tilde{\mathsf{h}}}{\partial \mathsf{C}}}+{\displaystyle \frac{\partial \tilde{\mathsf{h}}}{\partial \mathsf{q}}}{\displaystyle \frac{\mathrm{d}\mathsf{q}}{\mathrm{d}\mathsf{C}}},\hfill \end{array}$$

(93)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}\tilde{\mathsf{h}}}{\mathrm{d}\Theta }}& ={\displaystyle \frac{\partial \tilde{\mathsf{h}}}{\partial \Theta }}+{\displaystyle \frac{\partial \tilde{\mathsf{h}}}{\partial \mathsf{q}}}{\displaystyle \frac{\mathrm{d}\mathsf{q}}{\mathrm{d}\Theta }}.\hfill \end{array}$$

(94)

In the present case, it was computed with the help of the software Acegen, see [78,79], which allows to compute the tangent using Mathematica, combining symbolic and algebraic operations, automatic code generation, automatic differentiation, and simultaneous optimization of expressions.

5.2. Computational Examples

The behavior of the model is shown with the help of two computational examples in which inhomogeneous strain and stress states dominate the structural response. In both cases, experimental results in real components are compared to simulations obtained with the developed material model. The first example includes the effect of the temperature, and the second example treats the effect of aging.

5.2.1. Influence of Temperature

In this first example, the component part shown in Figure 17 is used. The specimens were provided by the manufacturer and were already several years old when the experiments were performed (i.e., a totally different batch is chosen). Thus, the component part has already reached the completely aged state in the simulation, $a\left(t\right)=1$. As shown in Figure 17a, the upper area of the component is fixed with a screw to the testing machine’s fixture and the lower area is pulled with the help of an inserted oval bolt. There are four different isothermal experiments at three different temperatures equal to $20 °\mathrm{C}$, $50 °\mathrm{C}$, and $70 °\mathrm{C}$ in total. Here, the load experienced by the component results from the measured relative displacement between upper region and bolt, which is shown in Figure 18 with the corresponding temperatures for each experiment. Since the experiments were performed at different time scales, the time axis in Figure 18 is normed to the end time of each experiment. For each experiment, the reaction force is measured with a force sensor as well as the surface deformation using a digital image correlation system.

In the following, the experiment is modeled with the commercial finite element program Abaqus as shown in Figure 17b. The component part is meshed with quadratic tetrahedral elements (10 nodes per element) and the bolt with linear, eight-noded hexahedral elements. The problem has totally 54,078 elements and 85,569 nodes. The upper area of the geometry, which would be in contact with the screw, is fixed ($\overrightarrow{u}=\overrightarrow{0}$) and on the bolt a displacement $\overrightarrow{u}=-u_z{\overrightarrow{e}}_z$ is applied, see Figure 17b right. The values of Figure 18 are considered for the displacement $u_z$ and temperature $\Theta =\mathrm{const}.$ for each simulation. The contact formulation “node-to-surface” is chosen for the contact between the bolt and the component part. For the tangential behavior, a friction coefficient of $0.5$ is assumed, see [80]. For the normal and tangential behavior, the standard options of Abaqus for the penalty method are chosen. The bolt is modeled as linear elastic steel (Young’s modulus $210 \mathrm{MPa}$ and Poisson’s ratio $0.3$) and the component part with the proposed material model. The initial value of the aging variable in each simulation is equal to $a(t=0)=1$, i.e., the completely aged state. The initial values for the remaining internal variables correspond to the state without previous deformation ${\overline{\mathbf{C}}}_{\mathrm{p}}\left(0\right)=\mathbf{I}$, ${\overline{\mathbf{C}}}_{\mathrm{v}}\left(0\right)=\mathbf{I}$, and $k\left(0\right)=0$. The simulations are carried out with automatic load step control, which are different for each simulation since the time scale varies as well. The values of the initial minimum and maximum time step for the four simulations are shown in

Table 9. Time step choices.

	Simulation 1	Simulation 2	Simulation 3	Simulation 4
$\mathsf{\Delta }{t}_0$	$0.1$	100	1	1
$\mathsf{\Delta }{t}_{\min }$	${10}^{-3}$	${10}^{-2}$	${10}^{-3}$	${10}^{-3}$
$\mathsf{\Delta }{t}_{\max }$	$0.1$	210	1	1

.

The results of the simulations and experiments are shown in Figure 19. There, the reaction force over the prescribed displacement of the bolt is displayed. In Figure 20 the deformation on the surface of the component part for the third loading case (${\Theta }_3=50 °\mathrm{C}$) is shown as an example, other results are similar here. Figure 20a displays the major strain at the final point in time of the experiment, and Figure 20b the simulation concerned. Values equal or higher than $5\%$ deformation are shown in red. Figure 20c displays the major strain for a chosen point on the surface over time for experiment and simulation. Since the experiments at room temperature ($\Theta =20 °\mathrm{C}$) were performed at two different strain rates, see Figure 19a,b, we can observe that the model is able to reproduce the rate dependence, see Figure 18 with $t_{\mathrm{end}1}=13.4 \mathrm{s}$, $t_{\mathrm{end}2}=9 \mathrm{h}$, although it presents difficulties reproducing completely the relaxation behavior in Figure 19b. In more recent invesigations, see [81,82], considering the Gaussian error propagation concept to estimate the uncertainty due to uncertain input data, it was shown that both the uncertain material parameters as well as the uncertain boundary conditions and geometrical data of the specimens have an essential influence on the simulation of real parts. Particularly, the totally unknown spatial and size distribution of micro-voids, see in this context [30], might have an essential effect for higher loads. This can be observed in all diagrams of Figure 19. However, the exact causes (contact problem with relaxation properties, micro-voids distribution with micro-crack initiation, …) cannot exactly substantiated. The temperature dependence can be considered comparing Figure 19a,c,d. Here, one can see that in these three cases the starting point of the experiment is difficult to reproduce. The reason lies in the inaccuracy of the boundary conditions at these start times, which do not correspond completely to the simulation of the experiment (e.g., contact, position of the bolt). Moreover, the experiments at $50 °\mathrm{C}$ and $70 °\mathrm{C}$ were performed in a thermal chamber which introduces an inaccuracy in the measurement of the digital image correlation system due to the fact that the measurements have to be performed through glass windows of the chamber. Particularly, for the environmental temperature of $70 °\mathrm{C}$, the simulation shows a slightly more stiff response above a certain loading state. This can be seen at the identified state in Figure 19d. In the future, one might consider changing the temperature-dependence in the model.

5.2.2. Influence of Aging

The second example shows the ability of the model to reproduce the aging dependence. To this end, specimens with a known age are necessary. In this case, cylindrical tubes with a hole were used in torsion experiments for the validation as shown in Figure 21a.

One end of the specimens was fixed, and on the other end a rotation was applied as shown in Figure 21c. The experiments were performed at the constant temperature of $27 °\mathrm{C}$ and with specimens with three different aging times: without aging, with twelve months natural aging and a completely (artificially) aged specimen. The reaction torque was measured over time with the help of the testing machine’s torque sensor and the deformation of the specimens was measured on the surface with the digital image correlation system.

The specimen was modeled in Abaqus standard as shown in Figure 21b. It is meshed with 18,608 quadratic, 20-noded hexahedral elements and has 97,759 nodes. For the boundary conditions, the upper clamping area of the specimen is fixed ($\overrightarrow{u}=\overrightarrow{0}$) and on the lower clamping, the rotation $\overrightarrow{\phi }={\overline{\phi }}_y{\overrightarrow{e}}_y$ as shown in Figure 21c is applied. The temperature is constant during the simulation and is equal to $27 °\mathrm{C}$. The proposed constitutive model is used assuming both a “young” specimens (initial aging variable $a(t=0)=0$ for the specimen without aging) as well as $a(t=0)=0.802$ for the specimen with twelve months natural aging (this value was determined solving the evolution equation of the aging variable at the constant temperature of $19 °\mathrm{C}$ for 12 months), and, thirdly, $a(t=0)=1$ for the completely aged specimen. The initial values for the remaining internal variables are ${\overline{\mathbf{C}}}_{\mathrm{p}}\left(0\right)=\mathbf{I}$, ${\overline{\mathbf{C}}}_{\mathrm{v}}\left(0\right)=\mathbf{I}$, and $k\left(0\right)=0$. The results of the simulation and experiment are shown in Figure 22 and Figure 23. Here, the reaction torque over time is displayed. One can observe a very good agreement between experiment and simulation. The completely aged specimen shows the higher difference between experiment and simulation. It is possible that the discrepancy between the observed and predicted results is due to the fact that no experiments with a completely aged material were used for the identification process. This would be a purely predictive outcome of the model. Figure 24 and Figure 25 display the major strain at the final moment of the experiment and the simulation for the specimen with twelve months natural aging as an example result. The other cases are the same. Values larger than $4\%$ are shown in red. Figure 24 shows a very good reproducibility of the experimental results by the simulation.

The major strain over time is displayed in Figure 25 for different points at the surface of the specimen. Here, we can see the good agreement between experiment and simulation as well.

5.2.3. Validation of the Model

The evaluation of the agreement between experimental results and simulation is performed with certain numerical values. Thus, we are interested in the deviation of the simulation results with respect to the experimental measurement. First, the absolute error ${\mathrm{err}}_{F\mathrm{abs}}$ and relative error ${\mathrm{err}}_{F\mathrm{rel}}$ in the force of the evaluations of the complex geometry in Figure 19 are defined by

$$\begin{array}{cc}\hfill {\mathrm{err}}_{F\mathrm{abs}}& =|F_{\exp }-F_{\mathrm{sim}}|,\hfill \end{array}$$

(95)

$$\begin{array}{cc}\hfill {\mathrm{err}}_{F\mathrm{rel}}& =100\times {\displaystyle \frac{{\mathrm{err}}_{F\mathrm{abs}}}{|F_{\exp }|}}.\hfill \end{array}$$

(96)

Second, the torque quantities in the cylindrical tube with the hole of Figure 22 and Figure 23 are defined as

$$\begin{array}{cc}\hfill {\mathrm{err}}_{M\mathrm{abs}}& =|M_{\exp }-M_{\mathrm{sim}}|,\hfill \end{array}$$

(97)

$$\begin{array}{cc}\hfill {\mathrm{err}}_{M\mathrm{rel}}& =100\times {\displaystyle \frac{{\mathrm{err}}_{M\mathrm{abs}}}{|M_{\exp }|}}.\hfill \end{array}$$

(98)

They are displayed in Figure 26 and Figure 27. In Figure 26b and Figure 27b we can see that the relative error is generally larger at the beginning than in the other course of the line. The two main reasons are (1) that the slope in the force-displacement diagram is larger for small displacements, and, because of that, small absolute errors are relatively large. (2) the measured force is approaching zero in certain time intervals at the beginning of the experiments, but not in the simulation, which causes high relative errors. Moreover, one can observe that the relative error is smaller for the torsion experiments at different aging times than for the complex geometry (since the specimens drawn on for the torsion experiments were produced in the same batch as the specimens used for the experimental material characterization). Since this is not the case for the component part, larger deviations are expected. The errors are also larger for the artificially aged specimen, see Figure 27. The model tends to overestimate the torque and force for a fully aged material since these values were not used in the identification process.

In the following, a validation metric is introduced to provide an indicator of the ability of the constitutive model to predict the real behavior of a component part made of Zamak. Since a value of the quantity of interest is measured over time for different processes instead of having one unique value for one experiment that is repeated several times, we decided to adapt the metric proposed in [83] to this case. There are $n_{\exp }$ different experiments with $n_p^i$ data points for the i-experiment. The experimental value of a quantity of interest is represented as $y_e$ and the result of the simulation as $y_m$ (response of the model). The maximum relative error E is used now as an error measure, which is defined as the maximum relative error of all the experiments

$$\begin{array}{c}\hfill E=\max |{ϵ}_i|,{ϵ}_i={\displaystyle \frac{{y_m^i}_j-{y_e^i}_j}{{y_e^i}_{\max }}}.\end{array}$$

(99)

Here, an experiment is denoted by i and a point of the experiment by j. The maximum value of the experiment i is represented by ${y_e^i}_{\max }$. The experiment with the maximum relative error is denoted with the index E. Similarly as proposed in [83], the confidence interval (CI) for this error is a normalized measure by the maximum experimental value evaluated at the point where the maximum relative error metric occurs

$$\begin{array}{c}\hfill \left|{\displaystyle \frac{\mathrm{CI}}{{y_e^E}_{\max }}}\right|={\displaystyle \frac{t_{0.05,\nu }}{\sqrt{n_p^E}}}\left|{\displaystyle \frac{s_{\mathrm{E}}}{{y_e^E}_{\max }}}\right|,\mathrm{with}{s}_E=\left({\displaystyle \frac{1}{n_p^E-1}}\sum _{j=1}^{n_p^E}{({y_m^E}_j-{\overline{y}}_e^E)}^2\right).\end{array}$$

(100)

Here, $t_{0.05,\nu }$ is the t-distribution for $90\%$ confidence, with $\nu =n_p-1$, and the mean value of the experimental points of an experiment i is written as ${\overline{y}}_e^i$. The value ${y_e^E}_{\max }$ has the goal to normalize the confidence interval in the same way as the error E to make them comparable. The model is developed for non-cyclic quasi-static processes up to a maximum strain of approximately $10\%$ with Zamak 5. Processes at different strain rates, temperatures, and aging times are included in the validation with influence of torsion and tension/compression. Thus, by applying the validation metric from Equations (99) and (100) to the measurements, the values of

Table 10. Relative error and confidence interval for every experiment.

	Exp. 1	Exp. 2	Exp. 3	Exp. 4	Exp. 5	Exp. 6	Exp. 7
E in −	$0.2067$	$0.2626$	$0.1816$	$0.3133$	$0.0574$	$0.0364$	$0.1007$
CI in −	$0.0473$	$0.0224$	$0.0053$	$0.0053$	$0.0119$	$0.0113$	$0.0122$

are obtained.

The maximum deviation equal to $E=31.33\%$ with a confidence interval of $0.53\%$ takes place in the fourth experiment and simulation (component part at $70 °\mathrm{C}$). Therefore, the simulation is able to predict the experiments with a maximum relative error of $31.33\%$ for deformations of at least $5\%$. This relative error is, in general, smaller for lower deformations and it is also connected to other factors such as the thickness of the specimen and other microstructural imperfections such as micro-voids, see [30].

6. Conclusions

The provided material model finally provides a comprehensive model for the simulation of components made of the zinc die casting alloymaterial Zamak 5, which is used for large strain applications, temperature ranges in automotive applications, and for the consideration of different aging conditions. Additionally, an implementation in the commercial finite element program Abaqus is presented to interested investigators. The material model implies rate-dependent mechanical processes, temperature, and aging dependence. In addition to the new development of the model, the presentation details how the material parameters occurring in it can be determined as thoroughly as possible from a variety of experiments. The evaluation of the correlation matrix is used for this purpose.

The validation of the model (material model and model for computation) is carried out using a complex component made of Zamak. The force-displacement curves for the fully aged component of the finite element simulations agree very well with the experimental results up to the point of failure of the component. To characterize the validation more precisely, a validation metric that was developed in the context of fluid mechanical simulations is transferred to the solid mechanics problem at hand. Using displacement data derived from full-field data of digital image correlation measurements, the quantified deviation in a tensile-torsion test of 70 $°\mathrm{C}$ leads to approx. 30% with a low confidence interval. This shows that the quantification of simulation results from finite element simulations is very possible. In addition, a validation measure for practical use is presented to enable users to quantify the simulation results in comparison to experimental data.