Introduction

Elastomer components are used in a wide range of technical applications because, among other things, they offer the following advantages in their use. In addition to the absence of plastic behaviour under high loads, their optimal vibration insulation is highlighted. Due to their chemical structure, elastomers are viscoelastic, so they are mainly used for dynamic loads. But this results in a significant disadvantage, which is the subject of the current work: Elastomers heat up strongly by dissipation if they are periodically mechanically loaded with sufficient amplitudes and frequencies under boundary conditions with insufficient heat removal. The temperature rise leads to a change in viscoelastic material properties and produces thermal strains or thermally induced stresses. Likewise, temperature change leads to ageing phenomena and changes in strength properties.

Concerning the experimental investigations a large number of studies can be found in the literature. Medalia [39] examine the mechanisms of heat production and presented some important effects that can occur during service. Experimental measurements and theoretical analyses of the self-heating process with respect to existence and stability of the steady state temperature distribution were shown by [42], for example. Investigations on the influence of the crosslinking structure and the filler content with respect to the self-heating of carbon black-filled elastomers containing homogeneous initial temperature conditions can be found in [45]. Rennar et al. [48] compare and study the effect of different fillers on dissipative heating. Experimental investigations of the sample surface and core temperatures with respect to varying starting temperatures, amplitudes and frequencies can be found in Refs. [12, 13]. Self-heating experiments using inhomogeneous initial temperature conditions were performed by [31], for example. Behnke et al. [3], on the other hand, considers carbon black-filled EPDM blends in his contribution. Guo et al. [17] examines the influence of the filler content and the initial load on the mechanical dissipation of SBR compounds. Studies on the critical limits of the self-heating temperature, to observe the initiation of damage, can be found in Ref. [28, 29], or for rubber-steel sandwich structures in [37]. Berardi et al. [4] utilise the experimental results to determine the heat source using an inverse method. In contrast, Masquelier et al. [38] develop an experimental procedure to determine dissipation based on the temperature field. In another study, Balandraud et al. [1] use the simultaneous measurement of the thermal and kinematic field for the thermomechanical characterisation. In addition to the experimental activities, a large part of previous contributions includes material modelling to describe the thermomechanical behaviour of elastomers mathematically. A large number of thermomechanical constitutive models of the material classes elasticity, plasticity, viscoelasticity and viscoplasticity have already been developed. These allow the representation of thermomechanical coupling effects by satisfying the dissipation inequality. In the 1970s, Chadwick [7] already formulated constitutive equations which describe the thermomechanical behaviour of solid polymers. An engineering approach for predicting the thermoviscoelastic behaviour of elastomeric dampers under cyclic loading was developed by Hausmann and Gergely [21]. Le Saux et al. [32] also propose an appropriate phenomenological approach. Within this work, three-dimensional models of thermoviscoelasticity and -plasticity are considered, which are in particular suitable for the representation of the thermomechanical behaviour of elastomers. These basically include the mechanical and thermal behaviour and the thermomechanical coupling due to the temperature dependence of parameters and evolution equations and the coupling due to the displacement dependency of dissipative stress power as well as thermoelastic coupling. The models cited can be distinguished on the basis of various criteria. These include, among others, the intermediate configurations introduced in the kinematics as well as the choice of the thermodynamic potential function and determination of the independent variables. Furthermore, the splitting of the potential function e.g. on the basis of kinematic or rheological assumptions as well as their specification by concrete approaches characterise the model. Thus, one obtains individual rules of calculation and material equations for the differently formulated models. Holzapfel and Simo [24] take the temperature-dependent nonlinear elastic behaviour in the form of volumetric-deviatoric finite thermoelasticity into account. For the representation of the linear rate dependence, internal variables of the stress type are used. In contrast, a more suitable concept from a thermodynamic point of view is used by the following authors. This involves the multiplicative decomposition of the deformation gradient into elastic and inelastic components in order to model the nonlinear rate dependence using internal variables of the strain type. Reese and Govindjee [47] develop on this basis for the first time a model of finite thermoviscoelasticity and take into account dissipative effects as well as thermoelastic coupling by means of a deformation- and temperature-dependent heat capacity. In particular, the temperature dependence of the viscoelasticity is modeled in the form of a nonlinear evolution equation. Lion [34] considers the elastic behaviour in terms of entropy elasticity and uses nonlinear temperature-dependent viscosities for the viscous behaviour. In addition, he uses a parallel decomposition to model the static hysteresis behaviour and viscoelasticity. Dippel [11] includes the temperature dependence of the rate-dependent behaviour within the Maxwell springs. In addition, further decompositions of the deformation gradient are used here in order to model temperature strains as well as the incompressible material behaviour. The potential function as well as its specification is adapted to experimental observations. Johlitz [27] further simplifies this model by adding the thermal and mechanical-volumetric part of the free energy, reformulates it and thus omits the thermal intermediate configuration for the representation of thermal strains. Netz and Hartmann [44] show in their model, among other things, the approximation of a highly nonlinear heat capacity and therefore the reduction of the thermoelastic coupling term. These models were further developed by Guo et al. [18], among others, in order to integrate further effects, such as the reorganisation of the network. Meo et al. [40] neglect the thermoelastic coupling effects, whereas Lin et al. [33] use different types of the dissipation inequality. In contrast to the phenomenological material models considered so far, Reese [46] uses a micromechanical approach. Using the transient network theory, the multiplicative decomposition of the deformation gradient and the temperature dependence of the mechanical parameters motivated at this level and formulates on this basis a continuum mechanical model. Behnke et al. [2] introduces an extended tube model for thermoviscoelasticity of rubberlike materials. Khan and Muliana [30], for example, uses the Gibbs potential to derive constitutive relationships and stress and temperature parameters of the material. Miehe [41] presents the basics of discretisation of thermomechanical processes and global algorithms for coupled problems. First, an initial boundary value problem is formulated for which an approximate solution can be calculated on the basis of variational principles, e.g. with the finite element method. These can be subdivided into one-, two- or multi-field variation principles. For this, the balance equations in a suitable representation, usually the quasi-static momentum balance and the first law of thermodynamics in the form of the heat conduction equation, are converted into the weak form and combined to a fully coupled functional. These functionals, which are nonlinear in temperature and displacement, are linearised to be solved with a Newton method. From this, the need for the definition of state vectors and tangent operators becomes apparent. This is followed by the consistent space and time discretisation and the choice of a suitable iteration procedure. The literature presented below provides examples of the steps from consistent material modelling to the implementation in the finite element method for various material models. Among others, [16, 52] dealt with the implementation of finite thermoviscoplasticity. The studies [6, 15, 25] contain various formulations that refer to the requirements for the implementation in commercial FE programs. Implementations for various model formulations of finite thermoviscoelasticity are discussed in the following studies. For example, the derived state vectors and their associated consistent tangent operators can be found in Refs. [2, 22, 47]. Reese and Govindjee [47] for the first time investigate the dissipative effects numerically. By means of several examples, they establish the relations between the influences of the model and the occurring phenomena. In contrast, validations based on experimental studies can be found in Refs. [2, 22] or [4]. The application of such material models in the form of FE-analyses includes not only the implementation but also the investigation and development of solution strategies. All investigations presented so far refer to a linear discretisation scheme, both in space and in time. Netz and Hartmann [44] show the advantages and applicability of higher order approaches with respect to the space and time discretisation of thermoviscoelastic problems. Hamkar [19] applies an iteration-free method of higher order time integration to increase the efficiency in the framework of finite thermoviscoelasticity. There are various possibilities to solve the coupled functionals. A detailed comparison of the different solution methods with respect to thermomechanical couplings can be found in Ref. [5] or [41], for example. As an alternative modelling approach, non-ordinary methods were also used. De Cazenove et al. [9, 10] for instance, use a harmonic damped structure analysis within a sequential iterative scheme to calculate the transient temperature distribution in a three-dimensional structure. This exemplary overview of the state of material modelling, experimental and numerical investigations of thermomechanical elastomer behaviour could by far not cover the relevant literature comprehensively. The publications mentioned represent only an exemplary selection. For further studies, reference is made to the extensive references given in the cited works.

The present contribution, on the other hand, focuses mainly on the model behaviour in relation to initial and boundary conditions as well as the parameterisation of such a model with the aim to provide the user with a deeper understanding of the self-heating effect and of being able to take this into account cost-effectively during the pre-development phase. In the second section, the used material model is presented, the constitutive relationships are introduced and a fully coupled functional is derived. The structure of the computational model, the final parameterisation and the analysis are carried out in the third section. Section four concludes with a discussion of the results and the modelling recommendations.

The mathematical model for thermoviscoelasticity

To formulate a thermomechanically consistent material model, the rheological element of modified finite thermoviscoelasticity illustrated in Fig. 1 is used. Here, the representation of thermal strains is omitted. Lion [34] showed that the thermo-elastic coupling effect under periodic stress only results in an oscillation of the temperature and its influence, therefore, has a subordinate role when considering the temperature evolution through dissipation. The specific free energy \(\varPsi\) stored in the system is additively decomposed into a part which is only depending on the temperature \(\varPsi _{\text {th}}\) with respect to the thermoelement, an equilibrium part \(\varPsi _{\text {eq}}\) and a non-equilibrium part \(\varPsi _{\text {neq}}\), which in turn can be assigned to the springs. Similarly, due to the parallel connection, the stresses of the branches are summed to obtain the total stress \(\mathbf {P}\). In particular, the equilibrium stress \(\mathbf {P}_{\text {eq}}\) and the overstress \(\mathbf {P}_{\text {neq}}\) correspond to the ground elasticity and the Maxwell model, respectively. The strain \(\mathbf {E}\) for the Maxwell model is composed additively into an elastic \(\mathbf {E}_{\text {e}}\) and an inelastic strain \(\mathbf {E}_{\text {i}}\) part. Moreover, the rate dependence of the damper element is described by using the temperature-dependent viscosity \(\breve{\eta }\left( \theta \right)\).

Fig. 1
figure 1

The rheological representation of the modified three-parameter model for the modified finite thermoviscoelasticity

Full size image

Kinematics

At this point, a short introduction to the kinematics is illustrated in Fig. 2. The consideration of large deformations requires the distinction between the reference and current configuration. \(\mathbf {F}\) denotes the deformation gradient.

$$\begin{aligned} \mathbf {F}=\frac{\partial \varvec{\chi }}{\partial \mathbf {X}}=\mathop {\mathrm {Grad}}\nolimits \left( \mathbf {x}\right) \end{aligned}$$
(1)

The motion \(\mathbf {x}=\mathbf {\varvec{\chi }}\left( \mathbf {X},t\right)\) describes the current position \(\mathbf {x}\) of the material point \(\mathbf {X}\) at time t in the current configuration. Furthermore, it is necessary to introduce a volumetric-isochoric intermediate configuration [14] to describe the quasi incompressible material behaviour of elastomers. The deformation gradient \(\mathbf {F}\) is split multiplicatively into a volumetric \(\bar{\mathbf {F}}\) and an isochoric component \(\hat{\mathbf {F}}\) with the aid of the determinant \(\det \left( \mathbf {F}\right) =J\):

$$\begin{aligned} \mathbf {F}=\hat{\mathbf {F}}\cdot \bar{\mathbf {F}}\quad \text{ with }\quad \hat{\mathbf {F}}=J^{-\frac{1}{3}}\mathbf {F}\quad \text{ and }\quad \bar{\mathbf {F}}=J^{\frac{1}{3}}\mathbf {I}. \end{aligned}$$
(2)

The isochoric part is further split into elastic and inelastic parts [36].

$$\begin{aligned} \hat{\mathbf {F}}=\hat{\mathbf {F}}_{\text {e}}\cdot \hat{\mathbf {F}}_{\text {i}} \end{aligned}$$
(3)
Fig. 2
figure 2

The reference configuration (RC), volumetric-isochoric intermediate configuration (VIC), elastic-inelastic intermediate configuration, current configuration (CC)

Full size image

Derivation of the potential expressions

Based on the Clausius–Duhem inequality, the free energy density is introduced in this section as a potential function \(\rho _0\,\varPsi\). The evaluation of the dissipation inequality leads to the potential relations and the thermomechanically consistent representation of occurring terms by satisfying the residual inequality. For this purpose, the right Cauchy–Green tensor \(\mathbf {C}=\mathbf {F}^{\text {T}}\cdot \mathbf {F}\), the isochoric elastic part of the right Cauchy–Green tensor \(\hat{\mathbf {C}}_{\text {e}}=\hat{\mathbf {F}}_{\text {e}}^{\mathrm {T}}\cdot \hat{\mathbf {F}}_{\text {e}}\) and the thermodynamic temperature \(\theta\) are defined as a set of variables.

$$\begin{aligned} \rho _{0}\,\varPsi =\rho _{0}\,\varPsi \big (\mathbf {C},\hat{\mathbf {C}}_\text {e},\theta \big ) \end{aligned}$$
(4)

In the following, the Clausius–Duhem inequality is shown with respect to the reference configuration.

$$\begin{aligned} -\rho _{0}\,\dot{\varPsi }+\mathbf {S}:\dot{\mathbf {E}}-\rho _{0}\,\dot{\theta }\,\eta -\frac{\mathbf {q}_0}{\theta }\mathop {\mathrm {Grad}}\nolimits \theta =\rho _{0}\,\tilde{\delta }\ge 0 \end{aligned}$$
(5)

The first term represents the negative temporal change of the free energy density. \(\rho _0\) is the density related to the reference configuration. \(\mathbf {S}:\dot{\mathbf {E}}\) represents the volume-related stress power where \(\dot{\mathbf {E}}\) is the time derivative of the Green-Lagrange strain tensor and \(\mathbf {S}\) the second Piola-Kirchhof stress tensor. In this context, \(\eta\) indicates the mass-related entropy and \(\mathbf {q}_0\) the Piola heat flux vector, the temperature gradient \(\mathop {\mathrm {Grad}}\nolimits \theta\) refers to the material coordinates \(\mathbf {X}\) and \(\tilde{\delta }\) is the dissipation per unit mass. The insertion of the time derivative of the free energy density (4) into the dissipation inequality (5) leads to:

$$\begin{aligned}&\mathbf {S}:\frac{1}{2}\dot{\mathbf {C}} -\rho _{0}\Bigg (\frac{\partial \varPsi }{\partial \mathbf {C}}:\dot{\mathbf {C}} +\frac{\partial \varPsi }{\partial \hat{\mathbf {C}}_\text {e}}:\dot{\hat{\mathbf {C}}}_{\text {e}} +\frac{\partial \varPsi }{\partial \theta }\dot{\theta }\Bigg )\nonumber \\&\quad -\rho _{0}\,\dot{\theta }\,\eta -\frac{\mathbf {q}_0}{\theta }\mathop {\mathrm {Grad}}\nolimits \theta \ge 0. \end{aligned}$$
(6)

The constitutive relations are obtained in a consistent way by satisfying the inequality. The application of Fourier’s law satisfies the condition of non-negativity for the last term, where \(\lambda\) is the heat conduction coefficient.

$$\begin{aligned} \mathbf {q}_{0}=-\lambda \mathbf {C}^{-1}\mathop {\mathrm {Grad}}\nolimits \theta \qquad \text{ with } \lambda \ge 0 \end{aligned}$$
(7)

Transforming \(\dot{\hat{\mathbf {C}}}_{\text {e}}\) under consideration of kinematic relations leads to the form shown below. \(\hat{\mathbf {L}}_{\text {i}}\) is the inelastic part of the spatial velocity gradient.

$$\begin{aligned} \Bigg [&\mathbf {S}-2\rho _0\Bigg (\frac{\partial \varPsi }{\partial \mathbf {C}}+(\det \mathbf {C})^{-\frac{1}{3}}\hat{\mathbf {F}}_\text {i}\cdot \frac{\partial \varPsi }{\partial \hat{\mathbf {C}}_\text {e}}\cdot \hat{\mathbf {F}}_\text {i}^{-\text {T}}\nonumber \\&-\frac{1}{3}\bigg (\hat{\mathbf {C}}_\text {e}^{\text {T}}\frac{\partial \varPsi }{\partial \hat{\mathbf {C}}_ \text {e}}:\mathbf {I}\bigg )\mathbf {C}^{-1}\Bigg )\Bigg ]:\frac{1}{2}\dot{\mathbf {C}}\nonumber \\&-\rho _0\Bigg [\eta +\frac{\partial \varPsi }{\partial \theta }\Bigg ]\dot{\theta }+2\rho _0 \frac{\partial \varPsi }{\partial \hat{\mathbf {C}}_\text {e}}\cdot \hat{\mathbf {C}}_{\text {e}}^{\text {T}}:\frac{ 1}{2}(\hat{\mathbf {L}}_{\text {i}}^{\text {T}}+\hat{\mathbf {L}}_{\text {i}})\nonumber \\&+\frac{\lambda }{\theta }(\mathop {\mathrm {Grad}}\nolimits \theta )\cdot (\mathbf {C}^{-1}\mathop {\mathrm {Grad}}\nolimits \theta )\ge 0 \end{aligned}$$
(8)

The evaluation of the inequality according to Coleman and Noll [8] yields to the constitutive equations.

$$\begin{aligned} \mathbf {S}&=2\rho _0\Bigg (\frac{\partial \varPsi }{\partial \mathbf {C}}+(\det \mathbf {C})^{-\frac{1}{3}} \hat{\mathbf {F}}_\text {i}^{-1}\cdot \frac{\partial \varPsi }{\partial \hat{\mathbf {C}}_\text {e}}\cdot \hat{\mathbf {F}}_\text {i}^{-\text {T}}\nonumber \\&\quad -\frac{1}{3}\bigg (\hat{\mathbf {C}}_\text {e}^{\text {T}}\frac{\partial \varPsi }{\partial \hat{\mathbf {C}}_\text {e}}: \mathbf {I}\bigg )\cdot \mathbf {C}^{-1}\Bigg ) \end{aligned}$$
(9)
$$\begin{aligned} \eta&=-\frac{\partial \varPsi }{\partial \theta } \end{aligned}$$
(10)

The remaining residual inequality has the form:

$$\begin{aligned} \rho _0\,\tilde{\delta }=&2\rho _0 \frac{\partial \varPsi }{\partial \hat{\mathbf {C}}_\text {e}} \cdot \hat{\mathbf {C}}_{\text {e}}^{\text {T}}:\frac{1}{2}(\hat{\mathbf {L}}_{\text {i}}^{ \text {T}}+\hat{\mathbf {L}}_{\text {i}})\nonumber \\&+\frac{\lambda }{\theta }(\mathop {\mathrm {Grad}}\nolimits \theta )\cdot (\mathbf {C}^{-1}\mathop {\mathrm {Grad}}\nolimits \theta )\ge 0. \end{aligned}$$
(11)

This implies the requirement that the specific mechanical dissipation \(\tilde{\delta }_{\text {M}}\) has to be greater than or equal to zero.

$$\begin{aligned} \rho _0\,\tilde{\delta }_{\text {M}}=2\rho _0 \frac{\partial \varPsi }{\partial \hat{\mathbf {C}}_\text {e}}\cdot \hat{\mathbf {C}}_{\text {e}}^{\text {T}}: \frac{1}{2}(\hat{\mathbf {L}}_{\text {i}}^{\text {T}}+\hat{\mathbf {L}}_{\text {i}})\ge 0 . \end{aligned}$$
(12)

To satisfy the residual inequality, the proportionality relation is introduced and leads to the following evolution equation for the tensor \(\hat{\mathbf {D}}_\text {i}\) and equivalent, for the inelastic right Cauchy–Green tensor \(\hat{\mathbf {C}}_{\text {i}}\):

$$\begin{aligned} \hat{\mathbf {D}}_\text {i}&=\frac{2}{\breve{\eta }(\theta )}\frac{\partial \varPsi }{\partial \hat{\mathbf {C}}_ \text {e}}\cdot \hat{\mathbf {C}}_\text {e}^{\text {T}} \end{aligned}$$
(13)
$$\begin{aligned} \Leftrightarrow \dot{\hat{\mathbf {C}}}_{\text {i}}&=\hat{\mathbf {F}}_{\text {i}}^{\text {T}} \cdot \Bigg [\frac{2}{\breve{\eta }^(\theta )}\frac{\partial \varPsi }{\partial \hat{\mathbf {C}}_\text {e}} \cdot \Bigg [\hat{\mathbf {C}}_\text {e}^{\text {T}}-\frac{1}{3}\text {tr} \Big (\hat{\mathbf {C}}_\text {e}^{\text {T}}\Big )\mathbf {I}\Bigg ]\Bigg ] \cdot \hat{\mathbf {F}}_{\text {i}}, \end{aligned}$$
(14)

where \(\hat{\mathbf {D}}_{\text{ i }}=\frac{1}{2}\left( \hat{\mathbf {L}}_{\text {i}}+\hat{\mathbf {L}}_{\text {i}}^{\text {T}}\right)\) represents the rate of inelastic deformation tensor. Furthermore, \(\breve{\eta }(\theta )\) is interpreted as a temperature-dependent viscosity function, which is expressed by the standard Williams–Landel–Ferry equation [54]:

$$\begin{aligned} \breve{\eta }\left( \theta \right) =\breve{\eta }_0\,\exp \left( \frac{C_1\left( \theta -\theta _0\right) }{C_2+\theta -\theta _0}\right) , \end{aligned}$$
(15)

where \(\breve{\eta }_0\) denotes the viscosity at the reference temperature \({\theta }_0\). Also \(C_1\) and \(C_2\) represent the WLF constants.

Material model

In this section, an overview of the free energy density \(\rho _0\,\varPsi\) is given. It is used for the thermomechanically consistent formulations of the constitutive equations under consideration of the Clausius–Duhem inequality,

$$\begin{aligned} \rho _0\,\varPsi \left( I_{\hat{\mathbf {C}}},I_{\hat{\mathbf {C}}_{\text {e}}},J, \theta \right) =&\rho _0\,\varPsi _{\text {th}}\left( \theta \right) +\rho _0\, \varPsi _{\text {eq}}^{\text {vol}} \left( J\right) \nonumber \\ {}&+\rho _0\,\varPsi _{\text {eq}}^{\text {iso}} \left( I_{\hat{\mathbf {C}}}\right) +\rho _0\,\varPsi _{\text {neq}} \left( I_{\hat{\mathbf {C}}_{\text {e}}}\right) . \end{aligned}$$
(16)

The temperature-dependent part of the free energy density is defined in accordance to [35]:

$$\begin{aligned} \rho _0\,\varPsi _{\text {th}}\left( \theta \right)&=\nonumber \\&\displaystyle {\rho _0\left[ \zeta _0-A\left( \theta \ln \frac{\theta }{\theta _0}-\left( \theta -\theta _0\right) - \frac{B}{2}\left( \theta -\theta _0\right) ^2\right) \right] }. \end{aligned}$$
(17)

The incompressible material behaviour must also be taken into account, which can lead to numerical difficulties [20]. For the approach [53], the Jacobian J is used as an independent variable. Accordingly, the volumetric free energy reads as

$$\begin{aligned} \rho _0\,\varPsi _{\text {eq}}^{\text {vol}}\left( J\right) =\frac{1}{2}\kappa \left[ \left( J-1\right) ^2+\left( \ln J\right) ^2\right] . \end{aligned}$$
(18)

At this point, the formulation of Johlitz [27] et al. should be mentioned, which includes the thermal components of the free energy in its volumetric part to formulate the thermally induced strains. Furthermore, the isochoric free energy density can be found in Ref. [43], for example. Here, a simple approach according to [43] is used depending on the first invariant of the isochoric or isochoricelastic right Cauchy–Green tensor \(I_{\hat{\mathbf {C}}}\) resp. \(I_{\hat{\mathbf {C}}_{\text {e}}}\):

$$\begin{aligned} \rho _0\,\varPsi _{\text {eq}}^{\text {iso}}\left( I_{\hat{\mathbf {C}}}\right)&=c_{10} \left( I_{\hat{\mathbf {C}}}-3\right) , \end{aligned}$$
(19)
$$\begin{aligned} \rho _0\,\varPsi _{\text {neq}}^{\text {iso}}\left( I_{\hat{\mathbf {C}}_{\text {e}}}\right)&=c_{e10} \left( I_{\hat{\mathbf {C}}_{\text {e}}}-3\right) \end{aligned}$$
(20)

The internal variables are described by the evolution equation (14). Using (20), it can be solved numerically according to a method proposed by Shutov et al. [51] and expressed further by the isochoric elastic left Cauchy–Green tensor \(\hat{\mathbf {B}}_e=\hat{\mathbf {F}}\cdot \hat{\mathbf {C}}_{\text {i}}^{-1}\cdot \hat{\mathbf {F}}^{\text {T}}\):

$$\begin{aligned} \dot{\hat{\mathbf {C}}}_{i}=\frac{4 \,c_{e10}}{\breve{\eta }(\theta )}\Bigg [ \hat{\mathbf {C}}-\frac{1}{3}\mathop {\mathrm {tr}}\nolimits \Big (\hat{\mathbf {C}}\cdot \hat{\mathbf {C}}_i^{-1}\Big )\hat{\mathbf {C}}_i\Bigg ]. \end{aligned}$$
(21)

Inserting the free energy densities (18)–(20) into the stress definition (9) and transforming the quantities to the current configuration, the Cauchy stress is obtained:

$$\begin{aligned} \mathbf {\varvec{\sigma }}=\mathbf {\varvec{\sigma }}_{\text {vol}}+\mathbf {\varvec{\sigma }}_{\text {eq}}+\mathbf {\varvec{\sigma }}_{\text {neq}}=\mathbf {\varvec{\sigma }}_{\text {vol}}+\mathbf {\varvec{\sigma }}_{\text {iso}}. \end{aligned}$$
(22)

The volumetric part of the Cauchy stress is:

$$\begin{aligned} \mathbf {\varvec{\sigma }}_{\text {vol}}=J^{-1}\kappa \left( J\left( J-1\right) +\ln J\right) \mathbf {I} . \end{aligned}$$
(23)

The equilibrium part of the Cauchy stress is:

$$\begin{aligned} \mathbf {\varvec{\sigma }}_{\text {eq}}=2J^{-1}c_{10}\left( \hat{\mathbf {B}}-\frac{1}{3}I_{\hat{\mathbf {B}}}\mathbf {I}\right) . \end{aligned}$$
(24)

The non-equilibrium part of the Cauchy stress is:

$$\begin{aligned} \mathbf {\varvec{\sigma }}_{\text {neq}}=2J^{-1}c_{e10}\left( \hat{\mathbf {B}}_{\text {e}}-\frac{1}{3}I_{\hat{\mathbf {B}}_{\text {e}}}\mathbf {I}\right) . \end{aligned}$$
(25)

The inelastic stress power is represented using variables of the current configuration:

$$\begin{aligned} \rho _0\,\tilde{\delta }_{\text {M}}=\frac{1}{\breve{\eta }\left( \theta \right) }\varvec{\sigma }_{\text {neq}}:\varvec{\sigma }_{\text {neq}}. \end{aligned}$$
(26)

The specific entropy is:

$$\begin{aligned} \eta =A\ln \frac{\theta }{\theta _0}+B\left( \theta -\theta _0\right) . \end{aligned}$$
(27)

With these formulations, it is now possible to formulate a fully coupled boundary and initial value problem using the momentum balance and the heat conduction equation, which can be solved with the finite element method.

Thermomechanically coupled functional

In the investigation of thermomechanical problems, the balance equations for mass, momentum, angular momentum and energy as well as the entropy inequality must be fulfilled at all times. The mass balance is satisfied by the relation \(\rho _0=J\,\rho\), the angular momentum balance by the characteristic of the first Piola-Kirchhoff stress tensor \(\mathbf {P}\cdot \mathbf {F}^{\text {T}}=\mathbf {F}\cdot \mathbf {P}^{\text {T}}\). First, the initial boundary value problem is formulated as a set of equations for which an approximate solution can be calculated, for example, with the finite element method based on a variational formulation. This comprises the material-independent basic equations, constitutive equations and initial boundary conditions. The material-independent local equations are the quasi-static balance of linear momentum to describe the mechanical behaviour and the local energy balance in representation of the heat conduction equation.

$$\begin{aligned} \mathop {\mathrm {Div}}\nolimits {\mathbf {P}}+\rho _0\,\mathbf {b}=\mathbf {0} \end{aligned}$$
(28)

Here \(\rho _0\,\mathbf {b}\) describes the body force per unit volume. The divergence operator \(\mathop {\mathrm {Div}}\nolimits\) is referred to the material coordinates. The first law of thermodynamics expressed by the heat conduction equation has the following structure:

$$\begin{aligned} \rho _0\, c\, \dot{\theta }=-\mathop {\mathrm {Div}}\nolimits {\big (\mathbf {q}_0\big )}+\rho _0\, r+\rho _0\, \tilde{\delta }_{\text {M}} \end{aligned}$$
(29)

\(c=A+B\,\theta\) is the linear temperature dependent specific heat capacity and r is the heat source related to the mass. The definition of initial conditions for the temperature field and the internal variable is required.

$$\begin{aligned} \theta \left( \mathbf {X},t_0\right)&= {}_0\theta \left( \mathbf {X}\right) \end{aligned}$$
(30)
$$\begin{aligned} \hat{\mathbf {C}}_{\text {i}}\left( \mathbf {X},t_0\right)&={}_0\hat{\mathbf {C}}_{\text {i}}\left( \mathbf {X}\right) \end{aligned}$$
(31)

Moreover, boundary conditions for displacement \(\mathbf {u}\) and temperature \(\theta\) on the Dirichlet boundary surfaces \(\partial \varOmega _{0u}\) and \(\partial \varOmega _{0\theta }\) as well as Piola heat flow \(\mathbf {q}_{0}\)- and stress vector \(\mathbf {t}_{0}\) on the Neumann boundary surfaces \(\partial \varOmega _{0q}\) and \(\partial \varOmega _{0\sigma }\) are defined with respect to the reference configuration:

$$\begin{aligned}&\mathbf {u}\left( \mathbf {X},t\right) =\bar{\mathbf {u}}\left( \mathbf {X},t\right) \quad \text {on}\quad \partial \varOmega _{0u} \end{aligned}$$
(32)
$$\begin{aligned}&\theta \left( \mathbf {X},t\right) =\bar{\theta }\left( \mathbf {X},t\right) \quad \text {on}\quad \partial \varOmega _{0\theta }\end{aligned}$$
(33)
$$\begin{aligned}&\mathbf {t}_0\left( \mathbf {X},t\right) =\mathbf {P}\,\mathbf {n}_0 =\bar{\mathbf {t}}_0\left( \mathbf {X},t\right) \quad \text {on}\quad \partial \varOmega _{0\sigma }\end{aligned}$$
(34)
$$\begin{aligned}&{q}_0=\mathbf {q}_0\cdot \mathbf {n}_0=\bar{{q}}_0\left( \mathbf {X},t\right) \quad \text {on}\quad \partial \varOmega _{0\sigma }. \end{aligned}$$
(35)

Elastomers are considered to be nearly incompressible, so that the use of the Galerkin approach would lead to pronounced volumetric locking. Therefore, a mixed formulation for nearly incompressible materials is used, which can be derived using the Hu–Washizu functional. The weak form of the equilibrium condition can also be obtained according to Hamilton’s principle as a result of the energy functional variation. However, this is only valid for conservative systems. Even though no stationary potential exists in the general case of inelastic material behaviour, Hartmann [20] as well as Holzapfel et al. [23] show the suitability of the multi-field formulation to be used for viscoelastic materials. Following Hartmann [20], the equilibrium part of the free energy density is taken to construct the functional. In the present study, the four independent fields are the displacement vector \(\mathbf {u}\), an assumed dilation measure \(\tilde{J}\), an assumed hydrostatic pressure p and the temperature \(\theta\). In particular, the Hu–Washizu functional for a general thermoelastic material can be written, by using the volumetric deviatoric split of the free energy density as:

$$\begin{aligned} \Pi _{\text {STP}}^{\text {TE}}\left( \mathbf {u},\theta ,\tilde{J},p\right) =&\int \limits _{\varOmega _0}\rho _0\left( \varPsi _{\text {dev}}\left( \tilde{\mathbf {C}},\theta \right) +\varPsi _{\text {vol}}\left( \tilde{\mathbf {C}},\theta \right) \right) \text {d}V_0\nonumber \\&+\int \limits _{\varOmega _0}p\left( J-\tilde{J}\right) \text {d}V_0\nonumber \\&-\int \limits _{\varOmega _0}\rho _0\,\mathbf {b}\cdot \mathbf {u}\,\text {d}V_0 -\int \limits _{\partial \varOmega _{0\sigma }}\mathbf {t}_{0}\cdot \mathbf {u}\,\text {d}A_0. \end{aligned}$$
(36)

Here, two further unknown scalar variables occur. The quantity p acts as a Lagrange multiplier to satisfy the constraint \(J=\tilde{J}=\det {\mathbf {F}}\) and \(\tilde{J}\) is used as an additional degree of freedom, which results from the introduction of the volumetric-deviatoric split of the modified deformation gradient \(\tilde{\mathbf {F}}\). Hence follows:

$$\begin{aligned} \tilde{\mathbf {F}}=\tilde{J}^{\frac{1}{3}}\mathbf {I}\cdot \check{\mathbf {F}}=\tilde{J}^{\frac{1}{3}}{J}^{-\frac{1}{3}}\mathbf {F}\quad \text {and}\quad \tilde{\mathbf {C}}=\left( \frac{\tilde{J}}{J}\right) ^{\frac{2}{3}}\mathbf {C}. \end{aligned}$$
(37)

The variation \(\delta \left( \circ \right)\) of the functional in terms of the mechanical fields and subsequent formulation in terms of the current configuration, respecting the boundary conditions, yields the following mechanical variation equations:

$$\begin{aligned} {\mathcal {M}}_{\text {u}}\left( \mathbf {u},\theta ,\tilde{J},p,\delta \mathbf {u}\right) =&\int \limits _{\varOmega }\left( \varvec{\sigma }_{\text {dev}}+p\mathbf {I}\right) :\mathop {\mathrm {grad}}\nolimits {\delta \mathbf {u}}\,\text {d}V\nonumber \\&-\int \limits _{\partial \varOmega _{\sigma }}\bar{\mathbf {t}}\cdot \delta \mathbf {u}\,\text {d}A\nonumber \\&-\int \limits _{\varOmega }\rho \,\mathbf {b}\cdot \delta \mathbf {u}\,\text {d}V=0 \end{aligned}$$
(38)
$$\begin{aligned} {\mathcal {M}}_{\text {p}}\left( \mathbf {u},\theta ,\tilde{J},p,\delta p\right) =&\int \limits _{\varOmega }J^{-1}\left( J-\tilde{J}\right) \delta p\,\text {d}V=0 \end{aligned}$$
(39)
$$\begin{aligned} {\mathcal {M}}_{\tilde{\text {J}}}\left( \mathbf {u},\theta ,\tilde{J},p,\delta \tilde{J}\right) =&\int \limits _{\varOmega }\left( \rho \frac{\text {d}\varPsi _{vol}\left( \tilde{J}\right) }{\text {d}\tilde{J}}-p\right) \delta \tilde{J}\,\text {d}V=0 . \end{aligned}$$
(40)

Equation (39) reproduces the constraint \(J=\tilde{J}\) and Eq. (40) gives the constitutive equation for the pressure p. The first equation (38) denotes the weak form of equilibrium with the relations \(\varvec{\sigma }_{\text {dev}}=\varvec{\sigma }_{\text {iso}}\) and \(p\,\mathbf {I}=\varvec{\sigma }_{\text {vol}}\). To obtain the variation formulation of the heat conduction equation, the heat conduction equation (29) is multiplied by the variation of the temperature \(\delta \theta\). The thermal functional \({\mathcal {T}}_{\theta }\) follows from the subsequent reformulation and insertion of the boundary condition (35):

$$\begin{aligned} {\mathcal {T}}_{\theta }(\mathbf {u},\theta ,\delta \theta )=&\int \limits _{\varOmega }\rho \, c\,\dot{\theta }\,\delta \theta \text {d}V -\int \limits _{\partial \varOmega _q}\bar{q}\,\delta \theta \,\text {d}A\nonumber \\&+\int \limits _{\varOmega }\mathbf {q}\cdot \mathop {\mathrm {grad}}\nolimits {\delta \theta }\,\text {d}V\nonumber \\&+\int \limits _{\varOmega }\rho \,(r+\tilde{\delta }_{\text {M}})\,\delta \theta \,\text {d}V=0. \end{aligned}$$
(41)

The thermomechanical problem is therefore completely formulated. All functionals (38)–(41) have to be fulfilled. The equations are combined to get a fully coupled functional \({\mathcal {G}}={\mathcal {M}}+{\mathcal {T}}=0\), which can be solved with the finite element method. The detailed FE implementation shall be dispensed within this paper and the reader is instead referred to [49].

Computational model

Fig. 3
figure 3

The model (side view) showing boundary and initial conditions

Full size image

This section contains the computational structure and the parameterisation of the material model described in Sect. 2. A thermomechanically fully coupled simulation procedure is used to solve simultaneously displacement and temperature fields with Newton’s method. The analysis is transient with respect to the heat conduction equation, whereas inertia effects are neglected. The heat transfer equation is integrated using a back-ward-difference scheme. The geometry of the structure to be simulated represents a cube with an edge length of 10 mm. A simple geometry is chosen so that results can be transferred to more complex geometries. A simple homogeneous deformation that is very well suited for theoretical considerations is the so-called simple shear test shown in Fig. 3. Here, different deformations within the mesh are avoided. In this example, the deformation that enters into the heat conduction equation produces a constant temperature at the centre node, such that dissipation effects can be studied mainly in isolation. The geometry is discretised by four C3D8HT (Three-dimensional 8-node trilinear coupled temperature-dis-placement hybrid elements with constant pressure) elements in the 1 and 2 directions and by twenty elements in the 3 direction. Furthermore, \(\bar{\mathbf {u}}_0\left( \mathbf {X},t\right)\) are displacement boundary conditions. Here \(\mathbf {X}\) and t denote the material position vector and the time, respectively. They are described as follows: \(\bar{\mathbf {u}}_0^{\text {LEF}}=\bar{\mathbf {u}}_0^{\text {RIG}}=\bar{\mathbf {u}}_0^{\text {FRO}}=\bar{\mathbf {u}}_0^{\text {BAC}}=\bar{\mathbf {u}}_0^{\text {TOP}}=\left[ u_1\,0\, 0\right] ^{\text {T}}\) and \(\bar{\mathbf {u}}_0^{\text {BOT}}=\left[ 0\,0\, 0\right] ^{\text {T}}\). The displacement vector \(\mathbf {u}\), on the other hand, serves as an excitation for the experiment. The prescribed displacement depends on the shear angle amplitude \(\tilde{\alpha }\), the excitation frequency \(\omega =2\pi \,f\) and the type of periodic signal. \(_{0}\theta \left( \mathbf {X}\right)\) is the initial temperature field and acts on the whole model. All \(\bar{\theta }_0^{\text {X}}\) are designated as isothermal temperature boundary conditions. Start temperature and boundary conditions are not distinguished within the calculation. If no temperature or special heat flow conditions are applied to a surface, it will automatically be considered as an adiabatic boundary condition (i.e. \(\bar{\mathbf {q}}_0=\mathbf {0}\)). The evaluation of the temperature always takes place at the centre node of the model. The commercial FE software Abaqus is used to perform the calculations. The constitutive relations and the corresponding tangent operators were implemented properly within the thermomechanically fully coupled UMAT interface. For a detailed overview of the implementation of thermomechanically coupled material models in the form of a UMAT, please refer to the article [49]. Table 1 displays the set of parameters that are used. These values, which are typical for a carbon black-filled natural rubber, were taken from [26].

Table 1 The parameter set used for the modified finite strain thermoviscoelastic material model
Full size table

Analysis

In this section, the respective calculations are further specified and evaluated. It also includes the presentation and analysis of these results. First, the influence of the displacement boundary condition is investigated. For this purpose, a sinusoidal displacement signal is applied \(\bar{\mathbf {u}}_0^{\text {TOP}}=\left[ u_0\sin \left( \omega \,t\right) \,0\,0\right] ^{\text {T}}\), where the displacement amplitude \(u_0\) is calculated over the angle \(\alpha\) and the height h of the cube \(u_0=\tan \left( \alpha \right) \,h=\gamma \,h\). The angular velocity is calculated as \(\omega =2\,\pi \,f\) with the frequency f. In the calculation with variable frequencies, the shear amplitude is \(40^{\circ }\), whereas in the variation of the amplitude, the frequency of 1 Hz is kept constant. The starting temperature as well as the isothermal boundary conditions are set to \(353\,\text {K}\). The temperature curve is recorded at the centre node of the cube versus time.

Fig. 4
figure 4

The amplitude variation of the finite thermoviscoelasticity

Full size image
Fig. 5
figure 5

The frequency variation of the finite thermoviscoelasticity

Full size image

Figures 4 and 5 depicts that the temperature in the steadystate equilibrium increases with higher amplitude as well as with higher frequency. The convergence to the steadystate is achieved when the dissipated energy is in a stationary state of balance with the heat transferred to the environment. In addition to the influence of frequency and amplitude, the investigation of critical loads is particularly important for the computations under large deformations. For this purpose, the normalised dissipated energy \(\tilde{\delta }_{\text {N}}=\frac{\tilde{\delta }_{\text {M}}}{\tilde{\delta }_{\text {Mmax}}}\) of a stationary cycle as a function of the shear rate \(|\dot{\gamma }|\) is shown in Fig. 6. The displacement vector is described here by a piecewise constant shear rate signal, known as a sawtooth pattern of various gradients \(\dot{\gamma }\,t\,h=u\left( t\right)\). Both the starting temperature and the isothermal boundary conditions are set to \(293\,\text {K}\) in this evaluation.

Fig. 6
figure 6

The shear rate variation of the finite thermoviscoelasticity

Full size image

This representation is reminiscent of the illustration of the loss modulus, however, it is suitable for large deformations. The greater the dissipation, the steeper the resulting temperature rise and thus the distance from the state of equilibrium. This novel representation is particularly suitable for use in non-uniform deformation fields and non-periodic load signals, as it is also valid beyond small amplitude and frequency. Now, the influence of the temperature boundary conditions is examined. For this purpose, the temperature-time curve is shown as a function of homogeneous starting temperatures in Fig. 7. The displacement boundary condition corresponds to the shear amplitude of \(40^{\circ }\) and a frequency of \(1\,\text {Hz}\). Here, it can also be observed that the temperature increase depends sensitively on the starting temperature. Similarly, the influence of the actual temperature field is observed. If we consider an isothermal characteristic curve, we can see that different states of equilibrium are achieved depending on the existing temperature field. This representation can be used to find the temperature range in which the material has the highest dissipation-sensitivity with regard to the expected loads. This representation concept can also be used for large deformations. In addition, the illustration shows the importance of taking the initial temperature and its distribution into account in the interpretation of test results.

Fig. 7
figure 7

A parameter variation of the finite thermoviscoelasticity with homogenous initial conditions for temperature

Full size image

Furthermore, the influence of inhomogeneous equilibrium temperature fields will be investigated. Basically, four differently stepped frequency sweeps are performed starting from the same homogeneous temperature field. Within each step, the steadystate equilibrium is calculated to compare them with each other. The shear amplitude \(\alpha =40^{\circ }\) and the isothermal boundary conditions on the cube sidewalls of \(293\,\text {K}\) are also applied here. Figure 8 clearly shows that if a general steadystate equilibrium is chosen as the starting condition for the simulation, the result will correspond to that of the respective homogeneous initial state. The proof of such behaviour is of paramount importance for the efficient experimental characterisation. In this way, equilibrium states can be generated and investigated for different loads with respect to a homogeneous initial condition.

Fig. 8
figure 8

The frequency sweep of the finite thermoviscoelasticity

Full size image

Finally, the influence on the steadystate equilibrium of the utilised material parameters is to be investigated. The material parameters will be varied in \(\pm \,10\,\%\) steps and the resulting temperatures will be compared with each other. This allows conclusions to be drawn from the sensitivity of the material parameters on the model. First of all, it is noticeable in Fig. 9 that a variation in density \(\rho _0\) or specific heat capacity c does not influence the final steadystate temperature. It can be demonstrated that the state of the equilibrium under periodic excitation is independent of these material parameters. A detailed proof can be found in the article [50]. The thermal conductivity \(\lambda\), the respective parameters of the Maxwell model \(c_{e10},\,\breve{\eta }\) and the relaxation time \(\tau =\frac{2\breve{\eta }}{c_{e10}}\) exhibit a significant influence concerning the expected state of equilibrium. Therefore, for the parameterisation of the constitutive model, these values should be determined experimentally, whereas for the others the use of literature parameters is sufficient. Thus, calculations can be accelerated enormously by the appropriate choice of density or heat capacity with the same accuracy of the equilibrium state.

Fig. 9
figure 9

Finite thermoviscoelasticity parameter variations and their influences on the steadystate temperature

Full size image

Conclusion

In this work, the influence of boundary and initial conditions as well as the parameterisation of the model on the self-heating effect was investigated using a modified model of finite thermoviscoelasticity. A simple model of modified finite thermoviscoelasticity was formulated. In addition to the introduction of a suitable kinematic description, a thermomechanically consistent derivation of the constitutive relations was presented. As part of the implementation, the unique initial and boundary value problem was formulated using an appropriate representation of the balance equations. A fully coupled functional was derived using a multi-field variational principle. Numerical studies have shown the dependence and characteristics of the temperature profile on frequency and amplitude. Furthermore, a representation of the critical shear velocities in relation to the dissipative characteristic was given, which is also suitable for finite deformations. The special influence of the initial temperature conditions on the self-heating process are also shown, from which recommendations for the experimental recording of the effect are derived. Finally, the influence of the model parameters was investigated and based on this, the experimental programme for characterisation can be adapted. On this basis, it can be concluded that the more efficient calculation of the temperature changes caused by dissipation, in addition to the significant reduction of the experimental effort, allows a wide range of users in the context of industrial applications to take into account and interpret the self-heating effect.