A thermo-mechanical material model for rubber curing and tire manufacturing simulation

Abstract

In this contribution, the phase change of unvulcanized to vulcanized rubber is described by a thermo-mechanical material model within the finite element method (FEM) framework. Before the vulcanization process (curing), rubber exhibits an elasto-visco-plastic behaviour with significant irreversible deformations without a distinct yield surface. After exposing rubber to high temperature, the molecular chains build-up crosslinks among each other and its mechanical behaviour changes to stiffer viscoelastic material. The proposed model assumes, that both phases are present during the vulcanization process. The ratio changes from the uncured phase at the beginning to the cured phase according to the current state of cure. A constitutive curing formulation is introduced into the model, to capture the shape change during the vulcanization and to ensure, that the second law of thermodynamics is fulfilled. A multiplicative split of the deformation gradient is assumed to describe incompressible material. Thermal expansion due to the change of temperature is taken into account in the volumetric part, as well as shrinkage during the vulcanization process. In the isochoric part, the phase change from elasto-visco-plastic to viscoelastic material is described by micro-macro transition based on the micro-sphere model. The consistent formulation of the material model and its tangent are important for a successful implementation into a three-dimensional finite strain FEM framework. The capabilities of the model are shown by the simulation of an axisymmetric tire production process starting at the green tire inserted into the heating press up to a post-cure inflation step.

1 Introduction

Rubber is used in a large variety of applications in industry. With certain fillers like carbon black or sulfur, its properties can be tuned for every specific application. In most cases, rubber is used in its vulcanized state above the glass transition temperature, where rubber has an elastic and incompressible behaviour. By far, the most important application for rubberlike elastomers are car tires. To predict the performance of a tire, engineers have created a lot of material models describing its behaviour at specific loading including temperature- and time-dependencies. However, in recent years, the development of more powerful simulation techniques and computer hardware led to the interest in simulating the production process of the product and the low crosslinked uncured rubber. Before rubber is exposed to the vulcanization temperature and the crosslinks of the molecular chains have not been built up, rubber is in its uncured phase and has a soft and fluid like behaviour. In [1], a strain energy density function for uncrosslinked butadiene rubber is proposed, where a polynomial strain energy density function is fitted to shear and uniaxial experiments. The distinct plastic and viscous flow of the uncured material are not addressed in this formulation. A more advanced step based on various experiments is presented in [2], where uniaxial tension, relaxation and compression tests have been carried out. They showed that plastic flow occured even at low strain without a distinct yield surface. Also, rate-dependent behaviour could be identified by a hysteresis at cyclic loading. They used a generalized Maxwell model to describing the viscoelastic behaviour ignoring the plastic flow. A development for the representation of uncured rubber has been done by [3], where the plastic flow and the viscoelastic part are decoupled from each other and considered separately. The model is capable to represent the short-term moulding effects and they used an endochronic plastic evolution law. For a more complex description, micro-macro transition via the micro-sphere model [4] has been used, which splits the three-dimensional deformation into a finite number of one-dimensional stretches. This material model is used in the contribution at hand to represent the uncured rubber before vulcanization. The phase change during vulcanization is addressed by few researchers. In [5], a small-strain model is given, where the increased stiffness is introduced by adding additional springs to the model. This represents the addition of crosslinks during curing, while the stress response of the newly added springs is only affecting the current and future strain increment and not the previous strains. A continuous model could be achieved by defining the elastic modulus as definite time integrals. They extended the model to finite strains [6] and to viscoelasticity [7]. While curing shrinkage is modelled by an additional element in the volumetric part of the strain energy, the shape change during vulcanization is not addressed in the approach. In [8], a so-called pseudoplastic strain is introduced that takes into account the permanent strain during vulcanization. They introduced that the behaviour of rubber is highly dependent on the amount of crosslinks in the rubber. For predicting the crosslinks, two chemical reactions are proposed, first the raw material is reacting with a vulcanization agent and the resulting product is further reacting into a fully crosslinked material.

In [8], the amount of crosslinks is determined by a moving die rheometer, where the reaction moment of the rubber during vulcanization is measured and related to the state of cure. Various researcher have been focused on describing the amount of cure over time at different temperatures. A widely used model is presented in [9], where the isothermal curing kinetics are described. They evaluated a differential scanning calorimeter to determine the amount of crosslinks over time at constant temperature. They used the measured current and total heat input in the sample to experimentally observe the curing reaction. A so-called induction period has been introduced in [10], which models the temperature-dependent time that elapsed before vulcanization starts. The model in [9] is developed further for non-isothermal cases and with a temperature dependent power term in the equation, see [11, 12]. These models are then applied to a tire cross-section in the mould, where only the heat transfer problem and the vulcanization process are addressed. Therefore, regions of over- or undercure could be identified neglecting the mechanical process of in-moulding and the phase change due to vulcanization. In [13], the stream of hot steam into the bladder has been depicted by computational fluid dynamics (CFD) simulations, resulting in a precise prediction of the temperature of the bladder that heats the tire. However, the fully thermo-mechanical simulation of a tire production process from the uncured green tire that is pressed into a mould, exposed to heat to enforce vulcanization and phase change of the tire to obtain the final cured tire, has not been addressed so far.

In this contribution, a thermo-mechanical consistent framework is presented for the phase change of rubber from the uncured to the cured phase. It is assumed that the behaviour is mainly dependent on the amount of crosslinks inside and that both phases contribute, dependent on the current state of cure, to the overall representation of the material. The phase change and the different geometries of the rubber before and after vulcanization are addressed, while always fulfilling the second law of thermodynamics. To ensure that no additional energy is evolving in the model during vulcanization, a so-called curing strain is introduced. This permanent strain is irreversible and will lead to a new equilibrium configuration after curing.

The paper is structured as follows. In Sect. 2, an evolution law for the state of cure for a non-isothermal curing process is introduced based on the model from [9]. The equations of motion in a continuum mechanical setting are presented in Sect. 3, as well as the fundamental ideas of the material model for the vulcanization process. A micro-sphere based approach and its implementation are discussed in Sect. 4. Its capabilities and the consistent derivation of the evolution laws are presented in Sect. 5, where the model is first applied to simple numerical examples. The production process of a tire from green tire to cured tire is shown and carried out by a finite element simulation.

2 State of curing

Vulcanization of rubber is assumed as crosslinking of the free molecular chains with each other. In this section, a temperature dependent function is introduced, which describes the amount of crosslinks c in the rubber, assuming no crosslinks $c=0$ for completely uncured rubber and a maximum amount of crosslinks $c=1$ for fully cured rubber. Breaking of crosslinks is not considered and only monotonic increasing of the state of cure (SOC) is modelled. Above a certain temperature, at around 403K, the curing process starts and the amount of crosslinks increases. The vulcanization process is highly temperature dependent. For larger temperatures, the process starts earlier and is faster than at lower temperatures.

Experimentally, this phenomenon can be captured via a moving die rheometer (MDR) test, see [8]. A rubber sample is placed in the machine, where a rotational displacement is applied at the top surface and the reaction moment is measured at the bottom of the sample. During this test, the sample is heated and vulcanization takes place. During vulcanization, the material becomes stiffer and the reaction moment increases and reaches a maximum at the end of vulcanization. Via post-processing, the state of cure is computed as ratio of the moment at the current time step divided by the maximum moment. Thus, it is possible to measure the time when vulcanization starts and also the velocity of the process at different temperatures. The rubber sample is small and thin enough to assume constant and homogeneous temperature over the whole body. During the experiment, the temperature of the sample will evolve over a certain time. It is assumed, that due to the small sample the rubber will heat quickly and will reach the temperature of the climate chamber immediately.

In [9], an empirical function for an isothermal curing process, that takes into account the rate- and temperature-dependent behaviour, is presented. This model represents the curing process very well for an isothermal state. During the moulding and heating of a tire in the press via surface heat exchange, the temperature increases slowly and the process cannot be assumed isothermal. Therefore, a model is introduced, where the curing rate is computed by the current state of cure

$$\begin{aligned} \dot{c} = k\, c_{t_n}^m\, \left( 1-c_{t_n}\right) ^n. \end{aligned}$$

(1)

The temperature dependency of the curing evolution is represented by an Arrhenius-type equation

$$\begin{aligned} k = k_0\, \text {exp}\left( -\frac{E}{R\,\theta }\right) . \end{aligned}$$

(2)

The gas konstant R is used and the parameters $k_0$ and E have to be determined via an optimization from experimental results. In [11, 12], a temperature dependent power term is proposed. The different rate of curing at the beginning of the vulcanization process is modelled via a linear temperature dependency of the power term

$$\begin{aligned} m =m_1\, \theta + m_2. \end{aligned}$$

(3)

The new state of cure is then obtained by an update formula with a predefined timestep $\varDelta t$

$$\begin{aligned} c_{t_{n+1}} = c_{t_n} + \dot{c} \, \varDelta t. \end{aligned}$$

(4)

Obviously, the rate of cure, Eq. (1), is zero when the state of cure is equal to one, which means, the rubber is fully cured, but also when it is equal to zero. Therefore, assuming no crosslinks at the beginning would lead to no evolution of the state of cure. It is assumed, that at the beginning of the vulcanization process, a small initial amount of crosslinks is already present

$$\begin{aligned} c_0 = 0.0002. \end{aligned}$$

(5)

This value represents a state of cure of $0.02\%$ at the beginning, which is small enough to ensure, that the unvulcanized rubber model can still be applied and is large enough to start the self-induced curing process. In Fig. 1, the simulated state of cure over time is plotted versus the experimental results from an MDR test for two tire compounds. The curing time and process can be captured nicely at the tested temperature with the proposed function and the behaviour at lower temperature is according to the results from the literature (Tables 1 and 2). To show the capabilities at various temperatures, the model is applied and compared to experimental results from [10]. The parameters are identified via optimization algorithm and are shown in Table 3. The state of cure is plotted over time for 4 different temperatures in Fig. 2. The temperature dependent curing model captures well the different curing kinematics at the tested temperatures.

Fig. 1

Numerically obtained state of cure over time compared to the experimental data determined by a moving die rheometer at $160^\circ \hbox {C}$ for a a tire tread compound and b a rubber containing the tire belts

Table 1 Curing parameters used in the simulation for a tire tread compound

Table 2 Curing parameters used in the simulation for a rubber containing the tire belts

Table 3 Curing parameters used in the simulation for a rubber containing the tire belts

Fig. 2

Numerically obtained state of cure over time compared to experimental data, provided in [10]

3 Thermo-mechanically consistent vulcanization framework

During the production process of a tire, both phases of the rubber are important. The uncured phase applies during moulding in the press and the cured phase in the post curing inflation. Therefore, both phases should be considered adequately. A thermo-mechanical consistent formulation has to be found to assemble both phases into one model in order to be able to describe intermediate transition states as well. The main idea for the proposed material formulation is to use already developed and proven material descriptions and connect them in a way, to fulfil the second law of thermodynamics and that a permanent curing strain is present.

3.1 Fundamental equations of motion

An initially undeformed body $\mathbb {B}_0$ is defined as reference configuration at time $t_0=0.$ The mapping $\varphi (\mathbf {X},t)$ links every material point from its initial position $\mathbf {X}$ to its spatial position $\mathbf {x}$ at time t

$$\begin{aligned} \mathbf {x} = \varphi (\mathbf {X}, t). \end{aligned}$$

(6)

The current deformation of one material point at time $t>t_0$ can be identified as difference of the current position and its reference position

$$\begin{aligned} \mathbf {u} = \mathbf {x} - \mathbf {X}. \end{aligned}$$

(7)

The deformation gradient $\mathbf {F} = \text {Grad} \varphi (\mathbf {X},t)$ maps an infinitesimal line element from the reference configuration to the spatial configuration

$$\begin{aligned} \mathrm {d}\mathbf {x} = \mathbf {F} \, \mathrm {d}\mathbf {X}. \end{aligned}$$

(8)

The gradient operator Grad$[\cdot ]$ denotes the spatial derivatives with respect to the reference configuration and grad$[\cdot ]$ with respect to the current configuration. The volume change $\nicefrac {V}{V_0}$ of an infinitesimal small volume element $\text{ d }V_0$ is represented by the determinant of the deformation gradient

$$\begin{aligned} J=\text {det} \, \mathbf {F} = \frac{V}{V_0}. \end{aligned}$$

(9)

To ensure that no self-penetration results from the deformation $\varphi$, the Jacobian must always be positive $J>0$. The right and left Cauchy Green tensor are derived as

$$\begin{aligned} \mathbf {C} = \mathbf {F}^T \, \mathbf {F}, \quad \mathbf {b} = \mathbf {F}\, \mathbf {F}^T \end{aligned}$$

(10)

by the deformation gradient. In rubberlike materials, incompressible behaviour is found. To enforce this feature, a multiplicative split of the deformation gradient into volumetric and isochoric parts

$$\begin{aligned} \mathbf {F} = \mathbf {F}_\text {vol}\, \overline{\mathbf {F}} \end{aligned}$$

(11)

is achieved by

$$\begin{aligned} \mathbf {F}_\text {vol} = J^{\nicefrac {1}{3}} \mathbf {1} \end{aligned}$$

(12)

and

$$\begin{aligned} \overline{\mathbf {F}} = J^{-\nicefrac {1}{3}}\, \mathbf {F}. \end{aligned}$$

(13)

3.2 Free energy function and Kirchhoff stress

The free energy function of the model is split additively into a volumetric part and an isochoric part

$$\begin{aligned} \varPsi (\mathbf {F},\theta , c, \mathcal {I}) = U(J,\theta ,c) + \overline{\varPsi }(\overline{\mathbf {F}}, \theta , c, \mathcal {I}). \end{aligned}$$

(14)

The volumetric part $U(J,\theta ,c)$ is a penalty term, which enforces the incompressibility of the model. The volumetric temperature expansion is associated with for the volumetric part as well as the curing shrinkage. The volumetric free energy function

$$\begin{aligned} U(J,\theta ,c)&= \kappa (J-\text {ln}(J)-1) \nonumber \\&\quad +\, \kappa \,\alpha _v \, \left( \theta _0-\theta \right) \, \text {ln}(J) + \kappa \,c\, \alpha _c \, \text {ln}(J) \end{aligned}$$

(15)

is represented by the rheological model depicted in Fig. 3. The bulk modulus $\kappa$ should be large enough to ensure incompressible behaviour. The volumetric thermal expansion is characterized by the parameter $\alpha _v$, which defines the volumetric change stemming from a temperature change of 1K. The last term of Eq. (15) describes the shrinkage due to vulcanization. The state of cure is monotonically increasing from $c_0$ to 1, therefore, the maximum curing shrinkage is defined by $\alpha _c$.

Fig. 3

Rheological representation of the volumetric material model: incompressible behaviour, thermal expansion and curing shrinkage are addressed in the volumetric part

It is assumed that during vulcanization some molecular chains are entangled and build connections with others while others can move freely. Therefore, the overall material response is represented by both phases depending on the current state of cure c. The isochoric free energy function is split further

$$\begin{aligned} \overline{\varPsi }(\overline{\mathbf {F}}, c, \theta , \mathcal {I})&= (1-c)\, \overline{\varPsi }_u(\overline{\mathbf {F}},\theta ,\mathcal {I}_u) \nonumber \\&\quad +\, c\, \overline{\varPsi }_c(\overline{\mathbf {F}}_m, \theta , \mathcal {I}_c), \end{aligned}$$

(16)

with a set of internal variables $\mathcal {I}_u$ and $\mathcal {I}_c$, into the uncured and cured rubber part, respectively. In the cured rubber model, the deformation gradient is decomposed further

$$\begin{aligned} \overline{\mathbf {F}} = \overline{\mathbf {F}}_m\, \overline{\mathbf {F}}_c. \end{aligned}$$

(17)

The permanent set due to curing is captured by the curing deformation gradient $\overline{\mathbf {F}}_c$, with the evolution law presented in the following section. In the undeformed and uncured state no curing strain is present, $\overline{\mathbf {F}}_c = \mathbf {1}$.

The associated Kirchhoff stresses are as well additively decomposed into volumetric and deviatoric contributions

$$\begin{aligned} \varvec{\tau }:= \varvec{\tau }_\text {vol} + \mathbb {P}:\overline{\varvec{\tau }}. \end{aligned}$$

(18)

In Eq. (18), the fourth-order deviatoric projection tensor is used

$$\begin{aligned} \mathbb {P}_{cd}^{ab} = \frac{1}{2}\left[ \delta _c^a\,\delta _d^b+\delta _d^a\,\delta _c^b\right] - \frac{1}{3}\,\delta ^{ab}_{cd}. \end{aligned}$$

(19)

The volumetric stresses are the negative pressure derived from the penalty term of the free energy function

$$\begin{aligned} \varvec{\tau }_\text {vol}= J\, \frac{\partial U(J,\theta ,c)}{\partial J} \mathbf {1}. \end{aligned}$$

(20)

The deviatoric stresses are obtained from the isochoric part of the free energy function

$$\begin{aligned} \overline{\varvec{\tau }} = 2\frac{\partial \overline{\varPsi }(\overline{\mathbf {F}},\theta ,c,\mathcal {I})}{\partial \mathbf {g}}. \end{aligned}$$

(21)

The covariant reference metric tensor $\mathbf {G}$ and the current metric tensor $\mathbf {g}$ are introduced. These mappings are identified by the Kronecker delta in the Cartesian coordinate system, $\mathbf {G} = \delta _{AB}$ and $\mathbf {g} = \delta _{ab}$. These metrics are necessary for the identification of the co- and contravariant objects in the Lagrangian and Eulerian manifolds. According to Eq. (16), the isochoric stresses are further decomposed depending on the state of cure

$$\begin{aligned} \overline{\varvec{\tau }} = (1-c)\, \overline{\varvec{\tau }}_u + c\, \overline{\varvec{\tau }}_c. \end{aligned}$$

(22)

From Eq. (22), it can be seen that for a fully cured material, i.e. $c=1$, the isochoric stresses are only derived from the cured free energy function $\overline{\varPsi }_c$.

3.3 Evolution law of the curing strain

Fig. 4

Thermo-mechanical framework of the proposed vulcanization model, where any reasonable material model can be inserted for the uncured and cured rubber

The basic underlying rheology of the framework is shown in Fig. 4, where the uncured and cured material are represented by so-called black boxes. The main feature of the vulcanization is its change of the mechanical behaviour and that a new equilibrium shape is present. As discussed previously, when the state of curing increases, the ratio of the stresses and free energy change from the uncured model to the cured one. For rubber-like materials, the curing process increases the stiffness of the material. For a thermo-mechanical consistent formulation, the free energy should not increase during the vulcanization process. Therefore, a symbolic curing element $\tilde{\mathbf {d}}_c$ is included in the lower branch, which also represents the shape change of the material during curing, see Fig. 4. Instead of defining an evolution law directly for the curing element, an evolution law is derived from the conservation of energy. In other words, the partial derivative of the free energy function with respect to the state of cure always has to be zero $\frac{\partial \overline{\varPsi }}{\partial c} := 0.$ Evaluating this statement with the free energy function Eq. (16) and the chain rule leads to

$$\begin{aligned} \frac{\partial \overline{\varPsi }}{\partial c} = \overline{\varPsi }_c - \overline{\varPsi }_u + c\, \frac{\partial \overline{ \varPsi }_c}{\partial \overline{\mathbf {b}}_m}\frac{\partial \overline{\mathbf {b}}_m}{\partial c}= 0. \end{aligned}$$

(23)

The evolution law for the curing strain ${\tilde{\mathbf {d}}}_c$ can be obtained by rearranging Eq. (23) and multiplying with the curing rate Eq. (1) from the previous section

$$\begin{aligned} {\tilde{\mathbf {d}}}_c =-\frac{\partial \overline{\mathbf {b}}_m}{\partial t} =-\frac{\partial \overline{\mathbf {b}}_m}{\partial c} \frac{\partial c}{\partial t} = -\left( \overline{\varPsi }_u - \overline{\varPsi }_c \right) \frac{\dot{c}}{c} \left( \frac{\partial \overline{\varPsi }_c}{\partial \overline{\mathbf {b}}_m}\right) ^{-1}. \end{aligned}$$

(24)

For a free energy function based on the first invariant, for example Yeoh’s model introduced in [14], the first part is easily obtained as a scalar multiplied by the identity tensor

$$\begin{aligned} \frac{\partial \overline{\varPsi }_c}{\partial \overline{\mathbf {b}}_m} = \frac{\partial \overline{\varPsi }_c}{\partial \overline{I}_{1}} \, \mathbf {1} \end{aligned}$$

(25)

using the relation $\frac{\partial I_A}{\partial \mathbf {A}}=\mathbf {1}$ from [15]. To project the curing strain into the stress direction, the projection tensor

$$\begin{aligned} \mathbf {N}_p = \frac{\varvec{\tau }_c}{\left| \varvec{\tau }_c\right| } \end{aligned}$$

(26)

is introduced with $\left| \varvec{\tau }_c\right| =\sqrt{\varvec{\tau }_c:\varvec{\tau }_c}$. The final constitutive evolution law for the curing strain is

$$\begin{aligned} -\tilde{\mathbf {d}}_c = \frac{1}{2}\mathcal {L}_v(\overline{\mathbf {b}}_m)\overline{\mathbf {b}}_m^{-1} = \left( \overline{\varPsi }_u-\overline{\varPsi }_c\right) \frac{\dot{c}}{c} \left( \frac{\partial \overline{\varPsi }_c}{\partial \overline{I}_{1}}\right) ^{-1}\,\mathbf {N}_p. \end{aligned}$$

(27)

The evolution law is introduced in a similar manner than in [16] using the Lie time derivative. The always positive curing rate and the assumption, that the free energy of the cured rubber is larger than the uncured, ensure, that the curing strain is monotonically increasing. For a fully cured material, the curing rate is zero and the curing strain is not increasing further. If the free energy of the cured rubber is significantly larger than the uncured, almost all of the deformation before the vulcanization process is transferred to the curing strain, resulting in a new equilibrium configuration. Obviously, Eq. (24) is not defined for a state of cure equal to zero. As stated in Eq. (5), a small amount of crosslinks is present even for uncured rubber and an initial state of cure $c_0>0$ is introduced. The proposed framework for the representation of the vulcanization process can be used for any given material model of uncured and cured rubber. Parameters for both phases can be identified independently by experiments that are done before and after the vulcanization process.

3.4 Thermo-mechanical consistency

It is assumed, that all energy during the phase change is used for the evolution of the permanent set. Therefore, the heat generation due to rubber curing, for example investigated in [17], is neglected in this contribution. For a constant temperature, the second axiom of thermodynamics, the so-called dissipation inequality, reads

$$\begin{aligned} \mathcal {D} = \varvec{\tau }:\mathbf {d} - \dot{\varPsi } \ge 0. \end{aligned}$$

(28)

The time derivative of the free energy function can be obtained by the product rule

$$\begin{aligned} \dot{\overline{\varPsi }} = \left( 1-c\right) \,\frac{\partial \overline{\varPsi }_u}{\partial \overline{\mathbf {b}} } \, \overline{\mathbf {b}} : \mathbf {d} + c\, \frac{\partial \overline{\varPsi }_c}{\partial \overline{\mathbf {b}}_m} :\dot{\overline{\mathbf {b}}}_m + \frac{\partial \overline{\varPsi }}{\partial c}: \dot{c}. \end{aligned}$$

(29)

Using the definition of the material time derivative with the found evolution law $\dot{\overline{{\mathbf {b}}}}_m=\mathbf {l}\,\overline{\mathbf {b}}_m+\overline{\mathbf {b}}_m\,\mathbf {l}^T+\overline{\mathcal {L}}_c(\overline{\mathbf {b}}_m)$, Eq. (28) is

$$\begin{aligned} \mathcal {D} =&\left[ \varvec{\tau } - \left( 1-c\right) \,\frac{\partial \overline{\varPsi }_u}{\partial \overline{\mathbf {b}} } \, \overline{\mathbf {b}} - c \frac{\partial \overline{\varPsi }_c}{\partial \overline{\mathbf {b}}_m}\overline{\mathbf {b}}_m \right] : \mathbf {d} \nonumber \\&- c \frac{\partial \overline{\varPsi }_c}{\partial \overline{\mathbf {b}}_m}: \left( \frac{1}{2} \overline{\mathcal {L}}_v(\overline{\mathbf {b}}_m) \overline{\mathbf {b}}_m^{-1} \right) - \frac{\partial \overline{\varPsi }}{\partial c}: \dot{c} \ge 0. \end{aligned}$$

(30)

The equation should be valid for every displacement rate and, therefore, the part in brackets is the definition of the stress, as stated in Eq. (22). The last part is equal to zero via definition of the evolution law. The remaining part is

$$\begin{aligned} \mathcal {D} = -c \frac{\partial \varPsi _c}{\partial I_{1}} \left( \varPsi _u-\varPsi _c\right) \frac{\dot{c}}{c} \left( \frac{\partial \varPsi _c}{\partial I_{1}}\right) ^{-1} \mathbf {1}:\mathbf {N}_p \ge 0 \end{aligned}$$

(31)

with Eq. (25). Due to the definition of the projection tensor and its trace $\mathbf {1}:\mathbf {N}_p=\text {tr}\left( \mathbf {N}_p\right) =1$, the dissipation inequality is fulfilled only if the free energy of the cured rubber is larger than of the uncured rubber, which is one of the main assumptions of the curing model. The model is based on the assumption of not changing free energy during curing, still an always positive internal work is present while the internal variable $\overline{\mathbf {F}}_c$ is evolving. Due to the definition of the monotonic increasing state of cure $\dot{c} > 0$, the dissipation is always positive and the process is irreversible.

4 Vulcanization model based on the micro-sphere model

The linearization of the proposed evolution law is rather difficult and can lead to a complex and numerically expensive implementation. Therefore, a micro-macro transition based on the micro-sphere model is proposed and disscused in this section. A three-dimensional deformation is represented by a number of uniaxial stretches pointing to the surface of a deformed micro-sphere. This approach simplifies the deformation state and even a more complex formulation of the material model can be used.

4.1 Kinematic description

The macroscopic deformation is linked to a single point of the micro-sphere via the Lagrangian unit orientation vector

$$\begin{aligned} \mathbf {r} = \text {cos}\, \varphi \, \text {sin} \, \vartheta \mathbf {e}_1 + \text {sin}\, \varphi \, \text {cos}\, \vartheta \, \mathbf {e}_2 + \text {cos}\, \vartheta \, \mathbf {e}_3. \end{aligned}$$

(32)

This vector with $\left| \mathbf {r}\right| =1$ and the Cartesian basis vectors $\mathbf {e}_1$, $\mathbf {e}_2$ and $\mathbf {e}_3$ points to the surface of the sphere, as illustrated in Fig. 5. The spatial orientation vector of the deformed sphere is introduced as the mapping with the isochoric deformation gradient

$$\begin{aligned} \mathbf {t} = \overline{\mathbf {F}}\,\mathbf {r}. \end{aligned}$$

(33)

The stretch of such a material line element in the orientation direction yields

$$\begin{aligned} \lambda = \sqrt{\mathbf {t} \cdot \mathbf {t}}. \end{aligned}$$

(34)

Homogenization over the unit sphere is done, where the orientation vectors are pointing to an infinitesimal small area $\text {d}A = \text {sin}\,\vartheta \,\mathrm {d}\varphi \, \mathrm {d}\vartheta$. The area of a part of the micro-sphere is expressed in polar coordinates

$$\begin{aligned} A(\varphi ,\vartheta ) = \int _0^\varphi \, \int _0^\vartheta \text {sin}\,\vartheta \,\mathrm {d}\varphi \,\mathrm {d}\vartheta . \end{aligned}$$

(35)

With the polar coordinates in the range of $\varphi \in \left[ 0,2\pi \right]$ and $\vartheta \in \left[ 0,\pi \right]$, the total area of the sphere is $\left| \mathcal {S}\right| = 4\pi .$ The averaging operator is defined as

$$\begin{aligned} \langle \cdot \rangle = \frac{1}{\left| \mathcal {S}\right| } \int _{\mathcal {S}} (\cdot )\, \mathrm{d}A. \end{aligned}$$

(36)

Assuming affinity between the macroscopic and microscopic stretches, $\lambda$ and $\tilde{\lambda }$, respectively, it holds

$$\begin{aligned} \lambda =\tilde{\lambda }. \end{aligned}$$

(37)

Therefore, the macroscopic free energy function can be identified as

$$\begin{aligned} \overline{\varPsi }(\overline{\mathbf {F}},\theta ,c,\mathcal {I}) = \langle \psi ( \lambda ,\theta ,c,\mathcal {P}) \rangle . \end{aligned}$$

(38)

The microscopic free energy is now a function of a one-dimensional stretch and a set of internal variables $\mathcal {P}$. The temperature and the state of cure are constant over the micro-sphere, therefore, for every unit direction the same temperature and state of curing are used for the microscopic free energy function.

Fig. 5

Micro-sphere model in spatial coordinates $\mathbf {e}_1$, $\mathbf {e}_2$ and $\mathbf {e}_3$

A discrete averaging operator is introduced for the numerical implementation of the model. The integral in Eq. (36) is replaced by a number of unit direction vectors and the accompanying weighting factors

$$\begin{aligned} \langle (\cdot )\rangle = \frac{1}{\left| \mathcal {S}\right| } \int _{\mathcal {S}} (\cdot )\, \mathrm{d}A \approx \sum _{i=1}^m (\cdot )^i w^i. \end{aligned}$$

(39)

A set of 21 directional vectors and corresponding weighting factors are taken from [18, 19] using the symmetry of the sphere. These values fulfil the isotropy condition and a stress free state in the reference configuration.

Fig. 6

Combined material model

4.2 Rheological model and free energy function

The proposed combined rheological model is shown in Fig. 6. For uncured rubber, the visco-elasto-plastic material model in [3] is used and, for cured rubber, the viscoelastic generalized Maxwell model is applied. Similar to the framework, see Fig. 4, a curing element is introduced in the cured rubber branch. The curing strain has to be computed in every unit direction of the micro-sphere model and saved as an internal variable $\mathcal {P}_i$. The microscopic free energy function is similar to Eq. (16) represented by a ratio of an uncured and a cured free energy function

$$\begin{aligned} \psi (\lambda ,\theta ,\mathcal {P}) = f_{EQ}(\theta ) \left[ (1-c)\,\psi _u(\lambda ,\mathcal {P}^u) + c\, \psi _c(\lambda ,\mathcal {P}^c)\right] . \end{aligned}$$

(40)

An entropy elastic thermo-mechanical free energy is employed with the temperature scaling function $f_{EQ}$, which is adopted from [16]. During the production process of a tire, the short term moulding process is not much effected by changing temperature in the rubber. The external heat sources from the press will have a larger influence than the internal heat sources from the viscous dissipation of the material. The temperature-dependent mechanical behaviour affects the tire production process, as the hotter material has a lower viscosity and flows into all edges of the mould. However, in this contribution, an isothermal free energy function with $f_{EQ}=1$ is assumed for simplification. Note, that the temperature still has an influence on the evolution of cure and on the volumetric expansion in Eq. (15).

4.2.1 Uncured rubber model

The model used for the uncured rubber is taken from [3] and is presented in short. As shown in [2], rubber before the vulcanization process shows very different behaviour than after curing. Without crosslinks of the molecular chains, the uncured rubber exhibits softer behaviour. The stress-strain response is highly nonlinear with rate-dependency of the initial loading and kinematic hardening. A hysteresis can be observed during cyclic loading and plastic flow is present even at small strains. A distinct yield surface cannot be identified as well as a ground state elastic response. In the used material model, see Fig. 7, the rate-independent elasto–plastic representation is seperated from the rate-dependent viscoelastic part and it is assumed that these two features exist independently from each other. The plastic flow is defined as rate-independent with a so-called endochronic plasticity law. For quasi-static loading, the stress–strain response shows kinematic hardening while plastic strain evolves, which is modelled by an additional spring element in parallel to the modified endochronic dashpot. In the rate-dependent branch, a nonlinear spring is coupled to a single dashpot in parallel to a Maxwell element. The single dashpot describes rate-dependency of the initial loading. The kinematic hardening part shows also rate-dependency, that the stress-strain response yields a larger slope for faster loading rates. This feature is captured by the additional Maxwell element in parallel to the dashpot, consisting of a nonlinear spring and a nonlinear dashpot. The according free energy function of the uncured rubber model is split into rate-independent and rate-dependent parts

$$\begin{aligned} \psi _u(\lambda , \mathcal {P}^u) = \hat{\psi }^p(\lambda , \mathcal {P}^p) + \hat{\psi }^v(\lambda , \mathcal {P}^v). \end{aligned}$$

(41)

A multiplicative split of the total stretch in the orientation direction is introduced for the rate-independent branch

$$\begin{aligned} \lambda = \lambda _{ep} \, \lambda _p. \end{aligned}$$

(42)

The logarithmic strains $\varepsilon := \text {ln}(\lambda )$ are introduced and additively decomposed into

$$\begin{aligned} \varepsilon = \varepsilon _{ep} + \varepsilon _p. \end{aligned}$$

(43)

The rate-independent branch for the uncured rubber model consists of two storage functions for the elastic response and the post-yield hardening part

$$\begin{aligned} \hat{\psi }^p(\lambda ,\mathcal {P}^p) = \psi _{ep}(\lambda _{ep}) + \psi _ {ph}(\lambda _p). \end{aligned}$$

(44)

Generic power type expressions are used for the free energy functions

$$\begin{aligned} \psi _{ep}(\lambda _{ep})= & {} \mu _{ep} \left[ \frac{\lambda _{ep}^{\delta _{ep}}}{\delta _{ep}} - \frac{1}{\delta _{ep}}-\text {ln}(\lambda _{ep})\right] , \end{aligned}$$

(45)

$$\begin{aligned} \psi _{ph}(\lambda _{p})= & {} \frac{\mu _{ph}}{\delta _{ph}} \left[ \lambda _{p}-1\right] ^{\delta _{ph}}, \end{aligned}$$

(46)

where $\mu _{ep}>0,\, \mu _{ph}>0, \, \delta _{ep}\ge 2$ and $\delta _{ph}\ge 2$ are material parameters. The free energy functions are always positive for every stretch and also the condition $\psi (\lambda =1) = 0$ is fulfilled. The stresses are defined as the derivatives with respect to the logaritmic strains and stretches

$$\begin{aligned} \sigma _p=\frac{\partial \psi _{ep}}{\partial \varepsilon } \quad \text {and} \quad \tau ^p = \frac{\partial \psi _{ep}}{\partial \lambda } = \frac{\sigma _p}{\lambda }. \end{aligned}$$

(47)

The evolution of the internal variable for the plastic strain is defined with the power-type evolution equation

$$\begin{aligned} \dot{\varepsilon }_p = \dot{\gamma }_p(\hat{\beta }_p) \hat{\beta }_p, \end{aligned}$$

(48)

with $\dot{\gamma }_p =\frac{\dot{z}}{\eta _p}\left| \hat{\beta }_p\right| ^{m_p}$ and the thermodynamical driving force

$$\begin{aligned} \hat{\beta }_p = -\frac{\partial \hat{\psi }^p}{\partial \varepsilon _p}. \end{aligned}$$

(49)

The evolution of the arclength $\dot{z}=\left| \dot{\varepsilon }\right|$ ensures a smooth hysteresis and an always positive plastic strain, which is not affected by the loading rate. For more information about the specific plasticity approach, the reader is referred to [3, 20].

Similar to the rate-independent branch, the total stretch is decomposed into elastic and viscoelastic stretches

$$\begin{aligned} \lambda = \lambda _{ve}\,\lambda _{v}, \end{aligned}$$

(50)

where the viscoelastic stretch $\lambda _v$ is further decomposed into

$$\begin{aligned} \lambda _v = \lambda _{eh}\,\lambda _h. \end{aligned}$$

(51)

The according logarithmic strains are additively decomposed as

$$\begin{aligned} \varepsilon = \varepsilon _{ve} + \varepsilon _v \end{aligned}$$

(52)

and

$$\begin{aligned} \varepsilon _v = \varepsilon _{eh}+\varepsilon _h. \end{aligned}$$

(53)

The free energy function of the rate-dependent branch is defined in a similar manner with two storage functions

$$\begin{aligned} \hat{\psi }^v(\lambda ,\mathcal {P}^v) = \psi _{ve}(\lambda _{ve}) + \psi _h(\lambda _{eh}). \end{aligned}$$

(54)

Generic power-type expressions, which are always positive and zero for no deformation, are employed

$$\begin{aligned} \psi _{ve}(\lambda _{ve})= & {} \mu _{ve}\left[ \frac{\lambda _{ve}^{\delta _{ve}}}{\delta _{ve}} - \frac{1}{\delta _{ve}}-\text {ln}(\lambda _{ve})\right] , \end{aligned}$$

(55)

$$\begin{aligned} \psi _h(\lambda _{eh})= & {} \frac{\mu _h}{\delta _h} \left[ \lambda _{eh}-1\right] ^{\delta _h}. \end{aligned}$$

(56)

The viscous strains evolve defined by the following evolution laws

$$\begin{aligned} \dot{\varepsilon }_v= & {} \dot{\gamma }_v(\hat{\beta }_v) \hat{\beta }_v \quad \text {with} \quad \dot{\gamma }_v(\hat{\beta }_v) = \frac{1}{\eta _v} \left| \hat{\beta }_v\right| ^{m_v}, \end{aligned}$$

(57)

$$\begin{aligned} \dot{\varepsilon }_h= & {} \dot{\gamma }_h(\beta _h) \beta _h \quad \text {with} \quad \dot{\gamma }_h(\beta _h) = \frac{1}{\eta _h} \left| \beta _h\right| ^{m_h}. \end{aligned}$$

(58)

The thermodynamical driving forces are introduced as the derivatives with respect to the logarithmic strains

$$\begin{aligned} \hat{\beta }_v = - \frac{\partial \hat{\psi }^v}{\partial \varepsilon _v} \quad \text {and} \quad \beta _h = -\frac{\partial \psi _h}{\partial \varepsilon _h}. \end{aligned}$$

(59)

Fig. 7

Uncured rubber material model and definition of the stretches, introduced in [3]

4.2.2 Cured rubber model

For the cured rubber material, the generalized Maxwell model is used, that captures well the viscoelastic properties of the cured rubber. The single nonlinear spring depicts long term elasticity that can be observed if the rubber is stretched and held until relaxation is finished. A finite number of Maxwell branches can be connected in parallel to model a wide range of frequencies. In this contribution, only one branch is considered to show the basic capabilities. However, an extension of the model is possible without much effort. The overall stretch is again decomposed into curing and elastic stretch and their logarithmic strains

$$\begin{aligned} \lambda = \lambda _e\,\lambda _c, \quad \varepsilon = \varepsilon _e+\varepsilon _c. \end{aligned}$$

(60)

For the Maxwell branch, the elastic stretch is further decomposed into elastic and inelastic stretches

$$\begin{aligned} \lambda _e = \lambda _{ce}\, \lambda _{ci}, \quad \varepsilon _e= \varepsilon _{ce} + \varepsilon _{ci}. \end{aligned}$$

(61)

The free energy functions for the cured rubber branch is

$$\begin{aligned} \hat{\psi }^c(\lambda _e, \mathcal {P}^c) = \psi _{ce}(\lambda _e) + \psi _{ci}(\lambda _{ce}). \end{aligned}$$

(62)

For both springs, the same power type function is used multiplied by a scaling factor $\gamma _\infty$ and $\gamma _1$ for the long term response and the viscous spring, respectively,

$$\begin{aligned} \psi _{ce} = \gamma _\infty \, \mu _c \left[ \frac{\lambda _e^{\delta _c}}{\delta _c}-\frac{1}{\delta _c}-\text {ln}(\lambda _e)\right] \end{aligned}$$

(63)

and

$$\begin{aligned} \psi _{ci} = \gamma _1 \, \mu _c \left[ \frac{\lambda _{ce}^{\delta _c}}{\delta _c}-\frac{1}{\delta _c}-\text {ln}(\lambda _{ce})\right] . \end{aligned}$$

(64)

Similar to the viscoelastic branch of the uncured rubber model, the evolution law of the dashpot is defined by the thermodynamically driving force $\hat{\beta }_v = -\frac{\partial \hat{\psi }^c}{\partial \varepsilon _{ci}}$

$$\begin{aligned} \dot{\varepsilon }_{ci} = \dot{\gamma }_{ci}(\hat{\beta }_v) \hat{\beta }_v \quad \text {with} \quad \dot{\gamma }_{ci} = \frac{1}{\eta _{ci}} \left| \hat{\beta }_v\right| ^{m_{ci}}. \end{aligned}$$

(65)

For the evolution of the curing strain, no evolution law is defined directly. As stated in Eq. (23), the free energy should remain constant during vulcanization. In the micro-sphere approach, this remains true for every orientation direction $\mathbf {r}$

$$\begin{aligned} \frac{\partial \hat{\psi }}{\partial c} = 0. \end{aligned}$$

(66)

This leads to the formulation of the evolution law of the curing strain

$$\begin{aligned} \dot{\varepsilon }_c = \frac{\psi _u- \psi _c}{c\, \frac{\partial \psi _c}{\partial \varepsilon _c}}\dot{c}. \end{aligned}$$

(67)

The derivative of the cured free energy with respect to the curing strain is obtained by the chain rule

$$\begin{aligned} \frac{\partial \psi _c}{\partial \varepsilon _c} = -\gamma _\infty \, \mu _c \left[ \lambda _e^{\delta _c}-1\right] - \gamma _1 \, \mu _c \left[ \lambda _{ce}^{\delta _c}-1\right] . \end{aligned}$$

(68)

The free energy of the viscoelastic branch depends on the overall stretch and the curing strain. For a more efficient implementation of Eq. (66), the total derivative $\frac{\mathrm {d} \psi _c}{\mathrm {d} \varepsilon _c}$ has to be used. The total derivative of the viscoelastic stretch with respect to the curing strain reads

$$\begin{aligned} \frac{\mathrm {d} \lambda _{ce}}{\mathrm {d}\varepsilon _c} = \frac{\partial \lambda _{ce}}{\partial \lambda _c}\frac{\partial \lambda _c}{\partial \varepsilon _c} + \frac{\partial \lambda _{ce}}{\partial \lambda _{ci}}\frac{\partial \lambda _{ci}}{\partial \varepsilon _{ci}} \frac{\partial \varepsilon _{ci}}{\partial \hat{\beta }_v}\frac{\partial \hat{\beta }_v}{\partial \lambda _{ce}}\frac{\mathrm {d} \lambda _{ce}}{\mathrm {d} \varepsilon _c}. \end{aligned}$$

(69)

Rearranging Eqs. (68) and (69) leads to the total derivative that is inserted into the evolution law of the curing strain, Eq. (67),

$$\begin{aligned} \frac{\mathrm {d}\psi _c}{\mathrm {d}\varepsilon _c}&= -\gamma _\infty \, \mu _c \left[ \lambda _e^{\delta _c}-1\right] \nonumber \\&\quad +\, \gamma _1\,\mu _c \frac{\lambda _{ce}^{\delta _{ce}}-1}{\lambda _{ce}} \, \frac{\frac{\partial \lambda _{ce}}{\partial \lambda _c}\frac{\partial \lambda _c}{\partial \varepsilon _c}}{1-\frac{\partial \lambda _{ce}}{\partial \lambda _{ci}}\frac{\partial \lambda _{ci}}{\partial \varepsilon _{ci}} \frac{\partial \varepsilon _{ci}}{\partial \hat{\beta }_v}\frac{\partial \hat{\beta }_v}{\partial \lambda _{ce}}}. \end{aligned}$$

(70)

The rather simple chosen material model for the cured rubber, captures its basic features sufficiently. However, there are different models like a non-linear viscoelastic model presented in [21] or [22]. It is necessary to derive the free energy function with respect to the curing stretch and any model can be inserted. Due to the fact that the curing strain is only affecting the cured material model, any model can be implemented for the uncured rubber without further derivations. In [4], a non-affine micro-sphere model for rubber is introduced, where a stretch fluctuation field on the micro-sphere is determined by a principle of minimum averaged free energies. This approach could also be used for modelling the rubber curing process, as the micro-stretches $\lambda$ for cured and uncured rubber are equal. However, even the evolution law of the simple affine visco-elastic model, Eq. (67), is quite complex. Increasing the complexitiy of the cured rubber model would increase the complexitity of the evolution law further.

4.3 Algorithmic stresses and tangents

In a fully coupled thermo-mechanical simulation, the temperature of the new timestep $\varDelta t_{n+1}$ is computed based on the heat transfer equation for every integration point, simultaneously. The current state of cure is obtained by Eq. (4), the state of cure at the previous timestep and the current temperature. The newly obtained curing rate, Eq. (1), will be used to update the curing strain. The new state of cure will lead to a different ratio of uncured to cured rubber, combined with the updated curing strain. The stresses and tangent moduli have to be derived for a new time step $t_{n+1}$, where the stress and internal state variables from the previous time step $t_n$ are known. The history variables are updated with the flow rule in the Eulerian setting.

The volumetric stress is the hydrostatic pressure, derived from the volumetric free energy function, Eq. (15),

$$\begin{aligned} \varvec{\tau }_\text {vol} = p\mathbf {1} \end{aligned}$$

(71)

with

$$\begin{aligned} p=J U'= \kappa \left[ J-1 + \alpha _v (\theta _0-\theta ) + c\, \alpha _c\right] . \end{aligned}$$

(72)

The volumetric tangent terms are derived further and take the form

$$\begin{aligned} \mathbb {C}_\text {vol} = 4\partial ^2_{\mathbf {g}\mathbf {g}} U = (p+s)\mathbf {g}^{-1} \otimes \mathbf {g}^{-1} - 2 p\mathbb {I}, \end{aligned}$$

(73)

using the second order derivative of U, the fourth order identy tensor $\mathbb {I}^{abcd}=\left[ \delta ^{ac}\delta ^{bd}+\delta ^{ad}\delta ^{bc}\right] /2$ and

$$\begin{aligned} s = J^2\,U'' = \kappa \left[ 1 - \alpha _v(\theta _0-\theta ) - c\,\alpha _c\right] . \end{aligned}$$

(74)

From Eq. (21) and the averaging operator, Eq. (36), the isochoric stresses

$$\begin{aligned} \overline{\varvec{\tau }} = \langle \tau \lambda ^{-1} \mathbf {t} \otimes \mathbf {t} \rangle \end{aligned}$$

(75)

are obtained with the identity $2\partial _{\mathbf {g}} \lambda = \lambda ^{-1}\mathbf {t} \otimes \mathbf {t}$. The further derivative of the stresses with respect to the metric and using the expression $\partial _{\mathbf {g}}(\mathbf {t} \otimes \mathbf {t}) = 0$ leads to the algorithmic tangent

$$\begin{aligned} \overline{\mathbb {C}}_\text {algo} = \langle (d-\tau \lambda ^{-1})\lambda ^{-2} \mathbf {t} \otimes \mathbf {t} \otimes \mathbf {t} \otimes \mathbf {t} \rangle . \end{aligned}$$

(76)

The scalar stresses in the microscopic setting are derived from the free energy function, Eq. (40), as

$$\begin{aligned} \tau = \frac{\partial \hat{\psi }}{\partial \lambda } = (1-c) \tau _u + c\, \tau _c = (1-c) (\tau _p+\tau _v) + c\, (\tau _{ce}+\tau _{c1}). \end{aligned}$$

(77)

Analogously, the scalar tangent terms are derived

$$\begin{aligned} d = \frac{\mathrm {d} \tau }{\mathrm {d} \lambda }= (1-c) d_u + c\, d_c = (1-c) (d_p+d_v) + c\, (d_{ce}+d_{c1}). \end{aligned}$$

(78)

The microscopic stresses and tangents of the uncured and cured rubber phase are determined in the next sub-sections. The isochoric tangent term is then obtained as

$$\begin{aligned} \mathbb {C}_\text {iso}&= \mathbb {P} : \left[ \overline{\mathbb {C}}_\text {algo} +\frac{2}{3}\right. (\overline{\varvec{\tau }}:\mathbf {g}) \mathbb {I} \nonumber \\&\quad \left. -\, \frac{2}{3}\left( \overline{\varvec{\tau }}\otimes \mathbf {g}^{-1} + \mathbf {g}^{-1} \otimes \overline{\varvec{\tau }}\right) \right] :\mathbb {P}^{T}. \end{aligned}$$

(79)

The total tangent is the sum of volumetric and isochoric tangent leading to a consistent linearization of the stresses

$$\begin{aligned} \mathbb {C} = \mathbb {C}_\text {vol} + \mathbb {C}_\text {iso}. \end{aligned}$$

(80)

4.3.1 Uncured rubber phase

In this section, the scalar stresses and tangents of the uncured rubber model are presented in short, for more details, the reader is referred to [3]. In direction $\mathbf {r}$, the logarithmic micro-stress of the elasto-plastic branch is derived as

$$\begin{aligned} \sigma _{p} = \frac{\partial \psi _{ep}}{\partial \varepsilon _{ep}} = \mu _{ep} \left( \lambda _{ep}^{\delta _{ep}}-1\right) . \end{aligned}$$

(81)

The scalar stresses are computed by the relation $\frac{\partial \varepsilon _{ep}}{\partial \lambda }=\lambda ^{-1}$ and the chain rule

$$\begin{aligned} \tau _{p} = \frac{\sigma _{p}}{\lambda } = \mu _{ep}\frac{\lambda _{ep}^{\delta _{ep}}}{\lambda }. \end{aligned}$$

(82)

Before computing the stresses, the plastic flow $\varepsilon _p$ in direction $\mathbf {r}$ needs to be known. The updated logarithmic plastic strain is derived from its value at the previous time step $t_n$ and the evolution law of the plastic flow Eq. (48)

$$\begin{aligned} \varepsilon _p^{t_{n+1}} = \varepsilon _p^{t_n} + \dot{\varepsilon }_p|_{t_n}\, \varDelta t. \end{aligned}$$

(83)

This steps leads to the nonlinear residual equation

$$\begin{aligned} r_p=\varepsilon _p^{t_{n+1}} -\varepsilon _p^{t_n} - \dot{\gamma }_p(\hat{\beta }_p) \hat{\beta }_p \varDelta t=0, \end{aligned}$$

(84)

that is solved by a local Newton iteration. In contrast to [3], the scalar tangent term is derived directly and is not obtained with the implicit function theorem

$$\begin{aligned} d_{p} = \frac{\mathrm {d} \tau _{p}}{\mathrm {d} \lambda } = \frac{\mu _{ep}}{\lambda } \frac{\mathrm {d} \lambda _{ep}^{\delta _{ep}}}{\mathrm {d} \lambda } - \frac{\tau _p}{\lambda }. \end{aligned}$$

(85)

The derivative of the elastic stretch is challenging, because it is also dependent on the backstress $\hat{\beta }_p$, which is dependent on the elastic stretch itself. That relation leads to the expression

$$\begin{aligned} \frac{\mathrm {d} \lambda _{ep}^{\delta _{ep}}}{\mathrm {d} \lambda } =&\frac{\partial \lambda _{ep}^{\delta _{ep}}}{\partial \lambda } + \frac{\partial \lambda _{ep}^{\delta _{ep}}}{\partial \varepsilon _p}\frac{\partial \varepsilon _p}{\partial \dot{z}}\frac{\partial \dot{z}}{\partial \lambda } + \frac{\partial \lambda _{ep}^{\delta _{ep}}}{\partial \varepsilon _p}\frac{\partial \varepsilon _p}{\partial \hat{\beta }_p}\nonumber \\&\left[ \frac{\partial \sigma _p}{\partial \lambda _{ep}^{\delta _{ep}}} \frac{\mathrm {d} \lambda _{ep}^{\delta _{ep}}}{\mathrm {d}\lambda } - \frac{\partial \beta _{ph}}{\partial \lambda _p} \left( \frac{\partial \lambda _p}{\partial \lambda } + \frac{\lambda _p}{\partial \lambda _{ep}^{\delta _{ep}}}\frac{\mathrm {d} \lambda _{ep}^{\delta _{ep}}}{\mathrm {d}\lambda } \right) \right] , \end{aligned}$$

(86)

with the definition of the backstress $\hat{\beta }_p = \sigma _p - \beta _{ph}$ and the stress in the kinematic hardening branch $\beta _{ph}=\mu _{ph}\left[ \lambda _p-1\right] ^{\delta _{ph}-1} \lambda _p.$ Due to the self dependency of the elasto-plastic stretch $\lambda _{ep}$, its derivative with respect to the overall stretch is on the left and the right hand side of Eq. (86). Rearranging leads to the expression to calculate the tangent

$$\begin{aligned} \frac{\mathrm {d} \lambda _{ep}^{\delta _{ep}}}{\mathrm {d} \lambda } = \frac{\frac{\partial \lambda _{ep}^{\delta _{ep}}}{\partial \lambda } + \frac{\partial \lambda _{ep}^{\delta _{ep}}}{\partial \varepsilon _p}\frac{\partial \varepsilon _p}{\partial \dot{z}} \frac{\partial \dot{z}}{\partial \lambda } - \frac{\partial \lambda _{ep}^{\delta _{ep}}}{\partial \varepsilon _p}\frac{\partial \varepsilon _p}{\partial \hat{\beta }_p} \frac{\partial \beta _h}{\partial \lambda _p}\frac{\partial \lambda _p}{\partial \lambda } }{1 - \frac{\partial \lambda _{ep}^{\delta _{ep}}}{\partial \varepsilon _p}\frac{\partial \varepsilon _p}{\partial \hat{\beta }_p} \left[ \frac{\partial \sigma _p}{\partial \lambda _{ep}^{\delta _{ep}}} - \frac{\partial \beta _h}{\partial \lambda _p} \frac{\partial \lambda _p}{\partial \lambda _{ep}^{\delta _{ep}}}\right] }. \end{aligned}$$

(87)

The evaluation of all the derivatives and the final form of the derivative can be found in the Appendix. The viscoelastic stresses

$$\begin{aligned} \sigma _v = \frac{\partial \psi _{ve}}{\partial \varepsilon _{ve}} = \mu _{ve}\left( \lambda _{ve}^{\delta _{ve}}-1\right) \end{aligned}$$

(88)

and the stress in the hardening branch

$$\begin{aligned} \beta _h = \mu _h\left[ \lambda _{eh}-1\right] ^{\delta _h-1} \lambda _{eh} \end{aligned}$$

(89)

are derived in a similar way. Due to the viscoelastic kinematic hardening branch, two internal iterations have to be done for the viscous logarithmic strain

$$\begin{aligned} r_v = \varepsilon _v^{t_{n+1}} - \varepsilon _v^{t_n} - \dot{\gamma }_v(\hat{\beta }_v) \hat{\beta }_v \varDelta t = 0 \end{aligned}$$

(90)

and an inner iteration for the hardening strain

$$\begin{aligned} r_h = \varepsilon _h^{t_{n+1}} - \varepsilon _h^{t_n} - \dot{\gamma }_h(\beta _h) {\beta }_h \varDelta t = 0. \end{aligned}$$

(91)

With the two strains at hand, the microscopic stress in direction of $\mathbf {r}$ reads

$$\begin{aligned} \tau _v = \frac{\sigma _v}{\lambda } = \frac{\mu _{ve}\left( \lambda _{ve}^{\delta _{ve}}-1\right) }{\lambda }. \end{aligned}$$

(92)

For the viscoelastic part, the tangent is even more complex, due to the second dashpot in the Maxwell element. Therefore, the total derivative of the hardening stretch $\lambda _h$ also depends on the viscous stretch $\lambda _{ve}$ and both depend on the backstress $\hat{\beta }_v$ and $\beta _h$, respectively. First, the derivative of the elastic stretch of the Maxwell element is computed as

$$\begin{aligned} \frac{\mathrm {d} \lambda _{eh}}{\mathrm {d}\lambda }&= \frac{\partial \lambda _{eh}}{\partial \lambda } + \frac{\partial \lambda _{eh}}{\partial \varepsilon _h}\frac{\partial \varepsilon _h}{\partial \beta _h}\frac{\partial \beta _h}{\partial \lambda _{eh}} \frac{\mathrm {d} \lambda _{eh}}{\mathrm {d}\lambda } \nonumber \\&\quad +\, \frac{\partial \lambda _{eh}}{\partial \lambda _{ve}^{\delta _{ve}}} \frac{\mathrm {d} \lambda _{ve}^{\delta _{ve}}}{\mathrm {d} \lambda }. \end{aligned}$$

(93)

The total derivative in Eq. (93) will now be transferred to the right hand side, and also the total derivative of $\lambda _{ve}^{\delta _{ve}}$ will be separated

$$\begin{aligned} \frac{\mathrm {d} \lambda _{eh}}{\mathrm {d}\lambda }= & {} \frac{\frac{\partial \lambda _{eh}}{\partial \lambda }}{1-\frac{\partial \lambda _{eh}}{\partial \varepsilon _h}\frac{\partial \varepsilon _h}{\partial \beta _h}\frac{\partial \beta _h}{\partial \lambda _{eh}} } + \frac{\frac{\partial \lambda _{eh}}{\partial \lambda _{ve}^{\delta _{ve}}} }{1-\frac{\partial \lambda _{eh}}{\partial \varepsilon _h}\frac{\partial \varepsilon _h}{\partial \beta _h}\frac{\partial \beta _h}{\partial \lambda _{eh}} }\frac{\mathrm {d} \lambda _{ve}^{\delta _{ve}}}{\mathrm {d} \lambda } \end{aligned}$$

(94)

$$\begin{aligned}= & {} a_1 + a_2\frac{\mathrm {d} \lambda _{ve}^{\delta _{ve}}}{\mathrm {d} \lambda }. \end{aligned}$$

(95)

The abbreviations $a_1$ and $a_2$ are introduced for simplification of the notation. The derivative of the stretch of the single spring is computed accordingly

$$\begin{aligned} \frac{\mathrm {d} \lambda _{ve}^{\delta _{ve}}}{\mathrm {d}\lambda } =&\frac{\partial \lambda _{ve}^{\delta _{ve}}}{\partial \lambda }\nonumber \\&+\frac{\partial \lambda _{ve}^{\delta _{ve}}}{\partial \varepsilon _v}\frac{\partial \varepsilon _v}{\partial \hat{\beta }_v} \left[ \frac{\partial \sigma _v}{\partial \lambda _{ve}^{\delta _{ve}}} \frac{\mathrm {d} \lambda _{ve}^{\delta _{ve}}}{\mathrm {d} \lambda } - \frac{\partial \beta _h}{\partial \lambda _{eh}}\frac{\mathrm {d} \lambda _{eh}}{\mathrm {d} \lambda } \right] . \end{aligned}$$

(96)

Inserting Eq. (93) into Eq. (96) leads to the total derivative of the viscoelastic stretch, which is needed for the algorithmic tangent

$$\begin{aligned} \frac{\mathrm {d} \lambda _{ve}^{\delta _{ve}}}{\mathrm {d}\lambda } = \frac{\frac{\partial \lambda _{ve}^{\delta _{ve}}}{\partial \lambda } - \frac{\partial \lambda _{ve}^{\delta _{ve}}}{\partial \varepsilon _v}\frac{\partial \varepsilon _v}{\partial \hat{\beta }_v}\frac{\partial \beta _h}{\partial \lambda _{eh}}\,a_1}{1-\frac{\partial \lambda _{ve}^{\delta _{ve}}}{\partial \varepsilon _v}\frac{\partial \varepsilon _v}{\partial \hat{\beta }_v}\left[ \frac{\partial \sigma _v}{\partial \lambda _{ve}^{\delta _{ve}}} - \frac{\partial \beta _h}{\partial \lambda _{eh}}\,a_2 \right] }. \end{aligned}$$

(97)

The derivation of all parts can be found in the Appendix. Finally, the tangent for the viscoelastic part is obtained by

$$\begin{aligned} d_v=\frac{\mathrm {d}\tau _v}{\mathrm {d} \lambda } = \frac{\mu _{ve}}{\lambda }\frac{\mathrm {d} \lambda _{ve}^{\delta _{ve}}}{\mathrm {d}\lambda } - \frac{\tau _v}{\lambda }. \end{aligned}$$

(98)

4.3.2 Cured rubber phase

For cured rubber, the internal iteration in the Maxwell element is solved in a similar way as Eq. (91). After computing the inelastic strains of the cured and uncured rubber phase, the free energy of all springs is known. After a time step $\varDelta t$ has elapsed and a large enough temperature is present, the state of cure increases according to Eq. (4). Thus, the ratio of uncured to cured free energy changes and to fulfil Eq. (23), an additional Newton iteration for the curing strain has to be carried out. Analogously to Eqs. (83) and (84), the residual equation for the logarithmic curing strain is

$$\begin{aligned} r_c = \varepsilon _c^{t_{n+1}} - \varepsilon _c^{t_n} - \dot{\varepsilon }_c \varDelta t = 0, \end{aligned}$$

(99)

which cannot be solved analytically for $\varepsilon _c^{t_{n+1}}$. Linerization is neccesary to obtain the new curing strain

$$\begin{aligned} \text {Lin }r_c = r_c|_{\varepsilon _c^{t_n}} + \mathcal {K}_c \, \varDelta \varepsilon _c^{k+1} = 0 \end{aligned}$$

(100)

with

$$\begin{aligned} \mathcal {K}_c = \frac{\mathrm {d} r_c}{\mathrm {d}\varepsilon _c}|_{\varepsilon _c^k}. \end{aligned}$$

(101)

The curing strain of the last time step is set as initial value for the update procedure $\varepsilon _c^0 = \varepsilon _c^{t_n}$. The curing strain of the next iteration step $k+1$ is computed as

$$\begin{aligned} \varepsilon _c^{k+1} = \varepsilon _c^k + \varDelta \varepsilon _c^{k+1} \end{aligned}$$

(102)

with

$$\begin{aligned} \varDelta \varepsilon _c^{k+1} = -\mathcal {K}_c^{-1} \, r_c|_{\varepsilon _c^k}. \end{aligned}$$

(103)

If a certain tolerance $\left| r_c\right| <TOL$ is reached, the update procedure is terminated and the curing strain has been found. For evaluating the tangent term $\mathcal {K}_c$, the evolution law of the curing rate is noted in short as

$$\begin{aligned} \dot{\varepsilon }_c = \frac{\psi _u-\psi _c}{c\, \frac{\mathrm {d} \psi _c}{\mathrm {d} \varepsilon _c}}\dot{c} = \frac{e_1}{e_2} \frac{\dot{c}}{c}. \end{aligned}$$

(104)

The tangent term of Eq. (101) using the quotient rule reads

$$\begin{aligned} \mathcal {K}_c = 1 - \frac{\frac{\mathrm {d} e_1}{\mathrm {d}\varepsilon _c}\,e_2-e_1\frac{\mathrm {d} e_2}{\mathrm {d} \varepsilon _c}}{(e_2)^2} \frac{\dot{c}}{c}. \end{aligned}$$

(105)

In the current framework, the free energy of the uncured rubber model does not depend on the curing strain $\varepsilon _c$, therefore, the derivative of the numerator is already known from Eq. (70). The quite complex derivation of the complete tangent term $\mathcal {K}_c$ is shown in the Appendix. In Table 4, the structure of all internal iterations is shown for updating the logarithmic strains. After all the iterations have been done, the internal stretches fulfilling the evolution laws are known and the stress of the cured rubber phase is computed as

$$\begin{aligned} \tau _c = \tau _{ce} + \tau _{ci} \end{aligned}$$

(106)

with

$$\begin{aligned} \tau _{ce} = \gamma _\infty \, \mu _c \frac{\lambda _{e}^{\delta _{c}}-1}{\lambda } \quad \text {and} \quad \tau _{ci} = \gamma _1\, \mu _c \frac{\lambda _{ce}^{\delta _{c}}-1}{\lambda }. \end{aligned}$$

(107)

The scalar tangent terms are derived consequently for the elastic and viscoelastic branch separately

$$\begin{aligned} d_{c} = \frac{\mathrm {d} \tau _c}{\mathrm {d} \lambda } = d_{ce} + d_{ci}. \end{aligned}$$

(108)

The microscopic stress of the elastic branch is derived from Eq. (107a)

$$\begin{aligned} d_{ce} = \frac{\mathrm {d} \tau _{ce}}{\mathrm {d} \lambda } = \gamma _\infty \,\mu _c\,\delta _c\frac{\lambda _e^{\delta _c-1}}{\lambda } \frac{\mathrm {d}\lambda _e}{\mathrm {d}\lambda } - \frac{\tau _{ce}}{\lambda }. \end{aligned}$$

(109)

The elastic stretch depends on the the total stretch, the curing strain and also on the viscoelastic strain

$$\begin{aligned} \frac{\mathrm {d} \lambda _{ce}}{\mathrm {d}\lambda } = \frac{\partial \lambda _{ce}}{\partial \varepsilon _{ci}} \frac{\partial \varepsilon _{ci}}{\partial \beta _i}\frac{\partial \beta _i}{\partial \lambda _{ce}}\frac{\mathrm {d} \lambda _{ce}}{\mathrm {d} \lambda } + \frac{\partial \lambda _{ce}}{\partial \lambda _e}\frac{\mathrm {d} \lambda _e}{\mathrm {d} \lambda }. \end{aligned}$$

(110)

Using the evolution law of the curing strain, Eq. (70), one obtains the final derivative of the elastic strain with respect to the overall stretch

$$\begin{aligned} \frac{\mathrm {d}\lambda _e}{\mathrm {d} \lambda } = \frac{z_1}{z_2}. \end{aligned}$$

(111)

For notation, the numerator $z_1$ and the denominator $z_2$ are written separately,

$$\begin{aligned} z_1 = \frac{\partial \lambda _e}{\partial \lambda }+\frac{\partial \lambda _e}{\partial \varepsilon _c}\frac{\dot{c} \, \varDelta t}{c\,\frac{\partial \psi _c}{\partial \varepsilon _c}}\left[ \frac{\partial \psi _u}{\partial \lambda } - \frac{\partial \psi _{ci}}{\partial \lambda _{ce}}\frac{\frac{\partial \lambda _{ce}}{\partial \lambda }}{1-\frac{\partial \lambda _{ce}}{\partial \varepsilon _{ci}}\frac{\partial \varepsilon _{ci}}{\partial \beta _i}\frac{\partial \beta _i}{\partial \lambda _{ce}}}\right] \end{aligned}$$

(112)

and

$$\begin{aligned} z_2&= 1+ \frac{\dot{c}\,\varDelta t }{\left( c\,\frac{\partial \psi _c}{\partial \varepsilon _c}\right) ^2} \nonumber \\&\quad \left[ \left( \frac{\partial \psi _{ce}}{\partial \lambda _e} + \frac{\partial \psi _{ci}}{\partial \lambda _{ce}} \frac{\frac{\partial \lambda _{ce}}{\partial \lambda _e}}{1-\frac{\partial \lambda _{ce}}{\partial \varepsilon _{ci}}\frac{\partial \varepsilon _{ci}}{\partial \beta _i}\frac{\partial \beta _i}{\partial \lambda _{ce}}} \right) \,c\,\frac{\partial \psi _c}{\partial \varepsilon _c} \right. \nonumber \\&\quad +\, \left. \left( \psi _u-\psi _c\right) \,c\,\frac{\partial ^2\psi _c}{\partial \varepsilon _c\,\partial \varepsilon _c}\frac{\partial \varepsilon _c}{\partial \lambda _e}\right] . \end{aligned}$$

(113)

The tangent term of the inelastic part of the cured rubber model is

$$\begin{aligned} d_{ci} = \frac{\gamma _1\,\mu _c}{\lambda }\, \delta _c\,\lambda _{ce}^{\delta _c-1} \frac{\mathrm {d} \lambda _{ce}}{\mathrm {d} \lambda } - \frac{\tau _{ci}}{\lambda }. \end{aligned}$$

(114)

The total derivative of the elastic stretch of the Maxwell branch is already given in Eq. (110) using the result from Eq. (113). The complete workflow of the material model is shown in Table 4. First, the new state of cure is computed with the current temperature. Then, all internal iterations in the uncured and cured rubber phase are carried out. With the found internal strains, the free energies of the uncured and cured rubber are known and the iteration to find the curing strain is carried out. After every residual reached the tolerance, the stresses and tangent moduli are computed.

Table 4 Simulation workflow for the thermo-mechanical curing model

5 Numerical examples

5.1 Vulcanization of a single material point

Fig. 8

Deformation and temperature driven curing process of a material point. The applied stretch and the prescribed temperature plotted versus time

In this section, the correct derivation of the vulcanization and the conservation of the free energy are addressed. Therefore, the deformation of a single material point is prescribed analytically by

$$\begin{aligned} \mathbf {F} = \begin{bmatrix} \lambda &{} 0 &{} 0 \\ 0 &{} \frac{1}{\sqrt{\lambda }} &{}0 \\ 0 &{} 0 &{} \frac{1}{\sqrt{\lambda }} \end{bmatrix}. \end{aligned}$$

(115)

The material parameters for the uncured rubber are taken from [21] and the material parameters for the cured rubber are arbitrarily chosen. The parameters are given in Tables 1, 6 and 7. Similar to the forming and building process of a tire, the deformation is applied first and kept constant as it is the case in the moulding press. Next, the temperature is increased and held for a longer time, thus, the vulcanization process is started. In Fig. 8, the prescribed stretch and temperature versus time are plotted. In this example, the deformation is kept constant for a long time, thus, the free energy during vulcanization is studied without viscoelastic effects of the uncured rubber. After 200 s, the final deformation of $\lambda =2.5$ is applied and the viscoelastic branch of the uncured rubber relaxes. The nonlinear effects of the elastoplastic and viscoelastic branch of the uncured rubber model cause a nonlinear increase of the free energy function. Due to the relatively slow loading rate, the viscoelastic effects are not significant in this example in contrast to a faster tire moulding process. The stored energy reduces until it is fully relaxed and is only represented by the time-independent elasto-plastic branch. After 1500 s, the temperature is increased and the state of cure evolves over time according to Eq. (4). Due to the evolution of the curing variable, the representation of the free energy changes but its amount remains constant

$$\begin{aligned} \psi = (1-c) \psi _u + c\, \psi _c. \end{aligned}$$

(116)

Cured rubber usually shows much stiffer behaviour than the uncured rubber and, to conserve the free energy during this transformation, curing strain evolves according to Eq. (67). As seen in Fig. 9, the curing strain of the first unit direction vector $\mathbf {r}^1 = \mathbf {e}_1$, which is equal to the uniaxial loading direction, evolves very fast. Due to the exponential free energy function of the cured rubber, the curing strain has to increase rapidly in comparison to the slow evolution of the state of cure to ensure that the free energy is not increasing. The portion of the cured rubber of the overall free energy $c\,\psi _c$ is increasing according to the state of cure, as seen in Fig. 9. The final value of the curing strain in loading direction is $\lambda _c=2.32$, which means that nearly 93% of the applied stretch is transformed into a permanent set, caused by the vulcanization process. In Table 5, the convergence behaviour of the internal iteration for the curing strain is shown. After a maximum number of 2 iteration steps, the residual equation is fulfilled, which verifies a correct derivation and implementation of the tangent term Eq. (101). In Fig. 10, the micro mechanical stress response of the first unit direction vector is shown. At the beginning of the simulation, the overall stress response is represented solely by the uncured rubber model. After 1500 s, the ratio of uncured to cured rubber is changed and the stress is increasing according to the state of cure.

Fig. 9

Resultant free energy and the evolution of the curing strain plotted versus time, the state of cure and the cured part of the free energy function evolves accordingly over time

Table 5 Residuals of the internal iteration for the curing strain

5.2 Rubber block simulation

After demonstrating the correct implementation on an integration point level, the model will be applied to a simple geometry. A fully coupled thermo-mechanical simulation will be carried out with a commercial finite element code, extended by user subroutines to model the material behaviour and evolution of cure. A tension test of a simple rubber block model, see Fig. 11a, is investigated, to show the capabilities of the model on a structural scale. The same material and curing parameters are chosen as in the first example, depicted in Tables 1, 6 and 7. The bottom nodes are fully constrained and the top nodes are driven vertically up to 100% structural deformation while constrained in the other two directions. In Fig. 11b, the deformed configuration is plotted with the displacements in the loading direction. After the total displacement has been applied, the temperature at the top and bottom will be increased up to $160^\circ \text {C}$ and, via conductivity, the temperature of the block will increase non homogeneous. Therefore, the state of cure will evolve non homogeneously and the block will be represented by the uncured and by the cured rubber model, simultaneously. After the block is fully cured, the boundary conditions are removed and the block is deforming to its new equilibrium configuration, Fig. 11c. Approximately 90% of the deformation is transferred to a permanent set and the shape of the block has significantly changed. In Fig. 12 the reaction force in loading direction is plotted versus time. The uncured rubber is softer than the cured rubber and, therefore, the reaction force is increasing while the state of cure is evolving.

Fig. 10

Uncured and cured stress contributions plotted versus time as well as the overall stress in loading direction

Fig. 11

Displacement of the block model in loading direction $u_x$ at the a undeformed configuration, b deformed state and c after vulcanization in the new equilibrium configuration

5.3 Tire production simulation

In tire building, all tread and belt compounds are assembled together on a drum. The carcass, innerliner and sidewall rubber compounds are wrapped around the steel bead and, then, stuck together with the other rubber parts of the belt and tread. The formed green tire on the drum starts to rotate and rollers press all parts together, so that all trapped air within the rubber compounds is squeezed out. At this point, the green tire is measured and a finite element model is created, see Fig. 13.

The fully coupled thermo-mechanical simulation is carried out based on a commercial finite element code, that is extended by user specific subroutines to include the material model and the evolution of cure. The subroutines are written in Fortran. Due to the very stiff behaviour of the mould, it is modelled by rigid surfaces. The simulation is then carried out in an axisymmetric setting, neglecting the tread pattern and only takes into account circumferential grooves. Therefore, the number of degrees of freedom is decreased and the simulation is completed faster in comparison to a full three-dimensional model. Taking into account lateral grooves of the tread pattern leads to large mesh distortions in this area, which have to be addressed by mesh smoothing algorithms like an Arbitrary Lagrangian-Eulerian framework or a meshless discretization often used in impact simulations. The whole production process simulation at hand takes around 3.5 hours on an $\hbox {INTEL}^{\textregistered }$ Core™i7 CPU @ 3.60GHz and 16GB RAM. The results are shown in the following section step by step.

Table 6 Mechanical properties of uncured rubber material used in the simulation

Table 7 Mechanical properties of cured rubber material used in the simulation

Fig. 12

Reaction force in loading direction $\hbox {RF}_x$ versus time

Fig. 13

Discretized green tire and its material sections

5.3.1 Moulding simulation

The short-term in-moulding process is highly dependent on the uncured phase of the rubber. Due to the large inner pressure, that is used to press the green tire into the mould, large strains are present especially in the tread area. Therefore, a correct derivation and implementation of the material is crucial to find a converged state. First, the green tire is placed inside the opened mould, Fig. 14a, contact between the elastic bladder and the green tire is established by applying a low inner pressure. Due to the low stiffness of the uncured nylon cords, the tire will expand until the side of the mould is closed, Fig. 14b. Up to this point, the tire center is held by Dirichlet boundary conditions to keep it inside the mould. After bladder and tire are in contact, this boundary condition is not active any more. The green tire is now only held by contact with bladder and mould. The mould is then closed up to its final position, Fig. 14c, still with the low inner pressure. Establishing contact in the tread area is difficult, therefore, the tire has to be discretiziesed very fine in this region. The tire will change its shape significantly, where the material is flowing into the mould to form the circumferential grooves. The final position of the tire, Fig. 14e, is achieved after the inner pressure of 25 bar is applied. While kept in this position, the viscoelastic branch of the uncured rubber relaxes and the material flows further into the edges of the mould. Validation of the in-moulding process is very difficult, because the press is a closed system and it is not possible to observe the tire during the moulding explicitly. In [23], temperature sensors are placed on the surface of the tire and, when contact is established between the tire and the hot mould, an instant increase of the temperature could be observed. The numerical results correspond to the experimental results, the first part that is in contact with the mould is the bead area. After that the sidewall near the belt edges establishes contact. At last, the area between bead and sidewall is in contact with the mould.

Fig. 14

Single steps of the in-moulding simulation: a inserting of the green tire into the mould, b inflation of the bladder and closing of the sides of the press, c fully closed moulding press with low inner pressure, d increased inner pressure and e finally applied inner pressure of 25 bar

5.3.2 Vulcanization simulation

By the model presented in Sect. 2, the description of unvulcanized rubber, exposed to large temperature, can be associated with the change of its mechanical behaviour. This process will take place at different times and at different rates depending on the mixture of the rubber compound. For heating the tire, the surface temperature of the moulding press and the air inside the bladder are increased and due to heat exchange over the tire surface, it is heated from the outside. For all rubber compounds, the same thermal properties, presented in Table 8, are taken from the literature. It is assumed, that the thermal parameters during vulcanization are constant. The temperature distribution of the cross-section over time is shown in Fig. 15a–h. The surface area in contact with the mould is heated directly and, therefore, faster than the inner surface in contact with the bladder. After 700 s, Fig. 15d, the sidewall has reached the final temperature of $160^\circ C$ while the tread area is significantly colder. This temperature difference in the tire cross-section is another reason for the non-uniform state of cure in the tire, shown in Fig. 16a–h. The parameters for the curing rate, Eq. (4), are fitted for every compound of the tire based on experimental data. The fast heating of the sidewall will cause the material to vulcanize there first while most of the other regions are still uncured, see Fig. 16c. The rubber, that contains the belts, shows vulcanization behaviour that starts early but at slower rate, which is observed in the curing process as well. After 1500 s, the SOC in the belt rubber starts to evolve while the tread material has not started to vulcanize, Fig. 16d. The faster curing process of the tread compound will overtake the belt rubber and after 2100 s, the tread is nearly fully cured, while the belt rubber is only half way cured, Fig. 16g. In the current model, the feature of overcuring is not represented, where the stiffness of the rubber will decrease again due to chemical alterations and breaking of established crosslinks.

In Fig. 17a, the temperature evolution of 3 different points are plotted. The temperature in the relatively small sidewall is increasing very fast in comparison to the thicker parts in the belt and subtread area. Point B at a circumferential groove is heated faster, due to the closer distance to the heating press. In Fig. 17b, the different curing kinematics for different rubber compounds can be seen. As presented in Fig. 1, the belt rubber starts the vulcanization process faster than the tread rubber. But the tread compound will reach the fully cured state faster and so, after 2500 s, both compounds are nearly at the same time fully cured. The curing kinematics of the sidewall and tread compound are similar, but due to the faster temperature evolution, the sidewall area is much faster fully cured.

Table 8 Thermal properties of the rubber material used in the simulation

Fig. 15

Contour plots of the temperature in the tire cross-section at a 100 s, b 300 s, c 500 s, d 700 s, e 900 s, f 1100 s, g 1300 s and h 1500 s

Fig. 16

Contour plots of the current state of cure in the tire cross-section at a 900 s, b 1100 s, c 1300 s, d 1500 s, e 1700 s, f 1900 s, g 2100 s and h 2300 s

Fig. 17

a Temperature and b state of cure of the 3 points of interest are plotted versus time: at the sidewall (Point A), at the belt material under the groove (Point B) and at the subtread (Point C)

5.3.3 Post curing inflation

After the tire has been fully cured, it is taken out of the mould and cooled down to room temperature again. During this step, the different heat expansion parameters of rubber and reinforcements lead to stress inside the tire. To minimize the pre-stress, tire manufacturers add an additional step in the production process, so-called post curing inflation (PCI). The tire is fixed on a rim and an inner pressure is applied while the tire is cooled down. Due to the constraints by the inner pressure, no free thermal shrinkage occurs and the PCI tire has a wider section width than a tire without post curing inflation [24]. Modelling this phenomenon requires a realistic model of the complex mechanics and phase transitions of the nylon cords that undergo large shrinkage during vulcanization and cooling. However, the tire after releasing from the moulding press has now a different shape than the green tire. In Fig. 18, the new equilibrium shape of the cured tire is compared to the stress-free reference configuration of the green tire. The tread grooves of the cured tire are sharper than in the uncured tire and represent the negative geometry of the mould. The width of the tire is decreased after the curing process. The shape of the cured tire could be achieved without changing the reference configuration.

Fig. 18

Comparison of the shape between the uncured green tire and the fully cured tire. (Color figure online)

6 Conclusion

Understanding the production process of a tire is neccessary to optimize the final product. The finite element method is a powerful tool for analysing the process and to identify areas of improvement in the curing process. Therefore, vulcanization of rubber is described by a model valid for rubber compounds of the tire. It could be observed that different rubber mixtures will vulcanize in a different way and, consequently, the tire curing process will be non-uniform over the tire cross-section. These results can then be used to optimize the tire curing process and to obtain a perfectly cured tire in a minimum time.

A thermo-mechanically consistent framework for the phase transition of the vulcanization of rubber has been presented, which includes the evolution of a permanent curing strain, the change of the mechanical behaviour and curing shrinkage. The presented framework is split into an uncured and a cured phase, therefore, the models for both phases can be chosen and identified independently from each other. The material model is implemented in a finite element framework and is used to simulate the tire production process. The procedure can be modelled in an axisymmetric setting from moulding of the green tire, heating of the tire inside the press and non-uniform vulcanization to the final cured tire. The production process is finished after the tire is released from the mould and a changed shape of the cured tire has been achieved. In the current axisymmetric model, lateral grooves are neglected. The shape of these is most important for the behaviour of the tire at wet surfaces and influences the noise of the running tire significantly.

The uncured rubber material flows into the small cavities in the mould due to high inner pressure from the bladder. In a pure Lagrangian setting, this causes large mesh distortions that will lead to numerical difficulties and an unstable simulation. Overcoming these challenges, mesh smoothing techniques, like the Arbitrary Lagrangian-Eulerian (ALE) formulation or a meshless representation have to be applied for modelling the fully three-dimensional production process with lateral grooves. The shrinkage behaviour of the nylon reinforcement cords, due to curing and cooling, has to be addressed to model the influence of the post-cure inflation.

Appendix

In order to provide all information required for the material modelling, additional remarks are given subsequently.

1.1 Appendix A Elastic-plastic Tangent

The total derivative for the rate-independent tangent requires the total derivative of the elasto-plastic stretches with respect to the overall stretch

$$\begin{aligned} \frac{\mathrm {d} \lambda _{ep}^{\delta _{ep}}}{\mathrm {d} \lambda } = \frac{z_1}{z_2}. \end{aligned}$$

(117)

The numerator is

$$\begin{aligned} z_1&= \delta _{ep}\,\lambda _{ep}^{\delta _{ep}-1} \nonumber \\&\quad \left( \frac{1}{\lambda _p} - \frac{\lambda _{ep}}{\eta _p}\left| \hat{\beta }_p\right| ^{m_p}\varDelta t \frac{\partial \dot{z}}{\partial \lambda } + (m_p+1)\dot{\gamma }_p\varDelta t \frac{\partial \beta _{ph}}{\partial \lambda _p} \right) \end{aligned}$$

(118)

and the denominator

$$\begin{aligned} z_2 = 1+\delta _{ep}\,\lambda _{ep}^{\delta _{ep}}(m_p+1)\dot{\gamma }_p\,\varDelta t \left( \mu _{ep}+\frac{\partial \beta _{ph}}{\lambda _p} \frac{\lambda _p}{\delta _{ep}\,\lambda _{ep}^{\delta _{ep}}}\right) , \end{aligned}$$

(119)

with the derivative of the evolution of the arclength

$$\begin{aligned} \frac{\partial \dot{z}}{\partial \lambda } = \frac{\text {ln}(\lambda )-\text {ln}(\lambda ^{t_n})}{\text {ln}(\lambda )-\text {ln}(\lambda ^{t_n})} \frac{1}{\varDelta t \, \lambda } \end{aligned}$$

(120)

and the derivative of the backstress

$$\begin{aligned} \frac{\partial \beta _{ph}}{\partial \lambda _p} = \mu _{ph}\left[ (\lambda _p-1)^{\delta _{ph}-1} + (\delta _{ph}-1)(\lambda _p-1)^{\delta _{ph}-2}\lambda _p\right] . \end{aligned}$$

(121)

1.2 Appendix B Viscoelastic Tangent

The rate-dependent branch of the uncured rubber model is more complex due to the viscoelastic hardening branch and the derivative of the elastic stretch with respect to the overall stretch is

$$\begin{aligned} \frac{\mathrm {d} \lambda _{ve}^{\delta _{ve}}}{\mathrm {d}\lambda } = \frac{ \delta _{ve}\, \frac{\lambda _{ve}^{\delta _{ve}}}{\lambda } + \delta _{ve}\,\lambda _{ve}^{\delta _{ve}}\,(m_v+1)\,\dot{\gamma }_v \, \varDelta t \,\frac{\partial \beta _h}{\partial \lambda _{eh}}a_1 }{ 1+\delta _{ve}\,\lambda _{ve}^{\delta _{ve}}\,(m_v+1)\,\dot{\gamma }_v\, \varDelta t \left[ \mu _{ve}-\frac{\partial \beta _h}{\partial \lambda _{eh}}\,a_2\right] } \end{aligned}$$

(122)

with the help of the substitution

$$\begin{aligned} a_1 = \frac{\frac{1}{\lambda _h\,\lambda _{ve}}}{1+\lambda _{eh}\,(m_h+1)\,\dot{\gamma }_h \, \varDelta t \frac{\partial \beta _h}{\partial \lambda _{eh}}} \end{aligned}$$

(123)

and

$$\begin{aligned} a_2 = \frac{-\frac{\lambda _{eh}}{\delta _{ve}\,\lambda _{ve}^{\delta _{ve}}}}{1+\lambda _{eh}\,(m_h+1)\,\dot{\gamma }_h \, \varDelta t \frac{\partial \beta _h}{\partial \lambda _{eh}}}. \end{aligned}$$

(124)

The used derivative of the backstress is consequently

$$\begin{aligned} \frac{\partial \beta _h}{\partial \lambda _{eh}} = \mu _h \left[ (\lambda _{eh}-1)^{\delta _{h}-1} + (\delta _h-1)\,(\lambda _{eh}-1)^{\delta _h-2}\,\lambda _{eh}\right] . \end{aligned}$$

(125)

1.3 Appendix C Evolution Law for the Curing Strain

The tangent for the local Newton iteration, Eq. (101), is discussed in detail. Substitution of the curing strain rate $\dot{\varepsilon }_c = \frac{e_1}{e_2}\,\frac{\dot{c}}{c}$ and the quotient rule lead to

$$\begin{aligned} \mathcal {K}_c = 1 - \frac{\frac{\mathrm {d} e_1}{\mathrm {d}\varepsilon _c}\,e_2-e_1\frac{\mathrm {d} e_2}{\mathrm {d} \varepsilon _c}}{e_2^2} \frac{\dot{c}}{c}. \end{aligned}$$

(126)

The derivative of the numerator has been already found, Eq (70),

$$\begin{aligned} \frac{\mathrm {d} e_1}{\mathrm {d} \varepsilon _c} = -\frac{\mathrm {d} \psi _c}{\mathrm {d} \varepsilon _c} . \end{aligned}$$

(127)

The derivative of the denominator is obtained by the derivative of $e_2=-\gamma _\infty \,\mu _c\, \left[ \lambda _e^{\delta _c}-1\right] +\frac{b_1}{b_2}$, where the numerator and denominator of the viscoelastic part are derived seperatly by the quotient rule

$$\begin{aligned} \frac{\mathrm {d} e_2}{\mathrm {d} \varepsilon _c} = -\gamma _\infty \, \mu _c\,\delta _c \lambda _e^{\delta _c} \frac{\mathrm {d} \lambda _{ce}}{\mathrm {d}\varepsilon _c} + \frac{\frac{\partial b_1}{\partial \varepsilon _c}b_2 - b_1\frac{\partial b_2}{\partial \varepsilon _c}}{b_2} \end{aligned}$$

(128)

with

$$\begin{aligned} b_1 = \gamma _1\,\mu _c\,\left[ \lambda _{ce}^{\delta _c}-1\right] \quad \text {and} \quad \frac{\partial b_1}{\partial \varepsilon _c} = \gamma _1\,\mu _c\,\delta _c\, \lambda _{ce}^{\delta _c-1}\,\frac{\mathrm {d} \lambda _{ce}}{\mathrm {d}\varepsilon _c}. \end{aligned}$$

(129)

The total derivative of the viscoelastic stretch is already found in Eq. (69), which leads to

$$\begin{aligned} \frac{\mathrm {d} \lambda _{ce}}{\mathrm {d}\varepsilon _c} = - \frac{\lambda _{ce}}{1+(m_v+1)\,\dot{\gamma }_v\,\varDelta t \, \gamma _1\,\mu _c\, \lambda _{ce}^{\delta _c}\,\delta _c}. \end{aligned}$$

(130)

The derivative of the denominator of Eq. (128) is finally

$$\begin{aligned} \frac{\mathrm {d} b_2}{\mathrm {d} \varepsilon _c} =&(m_v+1)\,\varDelta t \gamma _1\,\mu _c\, \delta _c \nonumber \\&\left[ \frac{m_v}{\eta _v}\frac{\left| \hat{\beta }_v\right| ^{m_v}}{\hat{\beta }_v}\gamma _1\,\mu _c\,\delta _c\,\lambda _{ce}^{\delta _c}\frac{\mathrm {d}\lambda _{ce}}{\mathrm {d}\varepsilon _c}\lambda _{ce}^{\delta _c}+\dot{\gamma }_v\,\delta _c\,\lambda _{ce}^{\delta _c-1}\frac{\mathrm {d}\lambda _{ce}}{\mathrm {d} \varepsilon _c}\right] . \end{aligned}$$

(131)