Thermodynamically consistent modeling of viscoelastic solids under finite strain

Abstract

The present article is concerned with modeling the viscoelastic behavior of Polydimethylsiloxane (PDMS) in large-strain regime. Starting from the basic principles of thermodynamics, an incremental variational formulation is derived. Within this model, the free energy density and dissipation function determine elastic and viscous properties of the solid. The main contribution of this paper is the estimation of the parameters in the proposed phenomenological model from measurements conducted on Sylgard184 samples. This non-linear material model simplifies to a Prony-series representation in frequency domain in case of small deformations. The coefficients of this Prony-series are detected from dynamical temperature mechanical analysis measurements. Time-temperature superposition allows to combine measurements at different temperatures, such that a sufficiently large frequency range is available for subsequent fitting of Prony-parameters. A set of material parameters is thus provided. The incremental variational formulation directly lends itself to finite element discretization, where an efficient and stable choice of elements is proposed for radially symmetric problems. This formulation allows to verify the proposed model against experimental data gained in ball-drop experiments.

1. Introduction

Silicone-based polymers, like Polydimethylsiloxane (PDMS), are known for their high flexibility and excellent damping properties. Thermal and chemical stability, corrosion resistance, and optical transparency add to the features of these materials. Not least their easy accessiblity at low cost and simple generation of samples allow for numerous applications in science and engineering. We cite the recent review papers by Ariati et al. (2021) and Zaman et al. (2019), both of which provide a great overview over the properties and applications of PDMS.

Typical applications of PDMS include, for example, microfluidics, where channels are required for filtering, separating and mixing fluids (cf. Borók et al., 2021), the insulation of electronic components, or coating of sensitive materials to protect them from external effects (cf. Eduok et al., 2017). In these applications, the damping nature of PDMS is of high advantage. An increase of the damping behavior could be realized through micro-structure design of sponges in Zhu et al. (2017). Treatment of surface layers was used to generate custom wrinkling patterns in specimen in Knapp et al. (2021).

The reliable realization of many of these applications requires accurate and efficient simulation techniques to predict the behavior of the system under consideration. Due to the high flexibility of PDMS, such a model must be fit to represent large deformations, and the damping properties have to be characterized properly.

Although frequently used in biomechanical or medical applications, the mechanical characterization of PDMS is not yet sufficiently accurate. Concerning the hyperelastic behavior, Nunes (2011) proposed a simple shear strain test to obtain material parameters for Mooney-Rivlin type materials. They proposed to use a variant to the classical Mooney-Rivlin material law to better fit their measurements. Later, Cardoso et al. (2018) fitted different hyperelastic material laws, including Mooney-Rivlin, Ogden, and Yeoh models. The influence of surface patterns on the elastic response was analyzed by Lee et al. (2021). Characterizations of the damping properties are scarce, Lin et al. (2009) experimentally detected Prony-parameters from punch tests and a specially developed dynamical mechanical analysis (DMA); however, the frequency range of this DMA is limited and the influence of temperature is not treated. Deguchi et al. (2015) measured storage and loss modulus at a single frequency of 1 Hz for a range of temperatures and several PDMS polymers. In their technical note, Long and Brown (2017) treat the applicability of time-temperature-superposition to Sylgard 184.

We propose an approach for modeling and subsequent simulation of viscoelastic solids, which we apply to Sylgard 184 ¹ in particular. Our approach ensures thermodynamic consistency in large deformation regime, as the incremental formulation is based on the Clausius-Duhem inequality and is defined not only through hyperelastic energy densities, but also through a dissipation function to model viscoelastic losses. In Kunzemann et al. (2023), a similar characterization was extended to include electromechanical coupling due to dielectric effects. To determine the material parameters used in these functions, a specimen of Sylgard 184 is subjected to dynamical temperature mechanical analysis (DTMA) in a limited range of frequencies. Through time-temperature superposition, a sufficiently large range of frequencies can be covered. From this data, a set of material parameters for a Prony-series expansion is found through a non-negative, non-linear least squares fit, similar to Kraus et al. (2017). These Prony-parameter directly correspond to the parameters of the proposed non-linear model.

In order to verify the proposed constitutive model, a finite element model of a ball-drop experiment is set up. The ball size and drop height are chosen in such a manner that large deformations will occur. The incremental variational formulation can directly be transferred to the context of finite element simulation. A classical stable pair of displacement-pressure elements is proposed to account for the incompressible nature of PDMS. Additional degrees of freedom for viscous strains satisfy the incompressibility condition exactly for radially symmetric problems. Time integration using the Newmark- $\beta$ scheme in variational context avoids any numerical damping.

The remainder of this paper is organized as follows: After the introduction, the constitutive model based on a hyperelastic potential, viscous contributions, and a dissipation function is discussed. Next, the correspondence of this non-linear model to a Prony-series expansion is discussed. The estimation of Prony-parameters from DTMA measurements is conducted using time-temperature superposition as well as non-negative least squares approximation. An incremental variational formulation, which replaces D’Alembert’s principle in dissipative processes, is derived and shown to be equivalent to classical viscoelastic models. This directly leads to the definition of a finite element model for radially symmetric problems in the following section. Finally, computational and experimental results are compared for a ball-drop experiment.

2. Viscoelastic modeling

2.1. Kinematics

We consider a body $\mathcal{B}$ consisting of a viscoelastic material. Each material point $\mathcal{P}$ is identified by its position vector $\mathit{X}$ in reference configuration. The motion of the body is described by a time-dependent mapping $\chi$ , which maps each material point to its current spatial position at time $t$ , $\mathit{x}=\chi (\mathit{X},t)$ . We use standard notations for all relevant kinematic quantities, starting with the displacement field $\mathit{u}=\mathit{x}-\mathit{X}$ . Throughout the following, we identify derivatives of fields with respect to reference coordinates $\mathit{X}$ by using capital operators. Derivatives with respect to spatial coordinates $\mathit{x}$ are identified using lower case operators.

The deformation gradient is then defined as $\mathbf{F}=\mathrm{Grad}\mathbf{\chi }=\mathbf{I}+\mathbf{Grad}\mathit{u}$ . For modeling large deformation viscoelasticity, an extension of the one-dimensional rheological standard solid model Haupt (2013) is proposed. In the standard solid model, also known as Wiechert model, a single spring element and $N$ spring-dashpot combinations in series, also known as Maxwell element, are connected in parallel, see also Figure 1.

Figure 1.

Linear, one-dimensional rheological Wiechert model for viscoelastic solids.

For the representation of incompressible elastomers, a splitting of the deformation gradient into its volumetric and isochoric components is adopted, with a hat denoting the isochoric part $\mathbf{F}=J^{1/3}\hat{\mathbf{F}}$ with $J=det\mathbf{F}$ . Based on the above definitions, the right Cauchy-Green tensor $\mathbf{C}$ and isochoric component of the right Cauchy-Green tensor $\hat{\mathbf{C}}$ are given by,

$$\mathbf{C}={\mathbf{F}}^{\mathbf{T}}·\mathbf{F},$$

(1)

$$\hat{\mathbf{C}}={\hat{\mathbf{F}}}^{\mathbf{T}}·\hat{\mathbf{F}}=J^{-2/3}\mathbf{C}.$$

(2)

For incompressible materials, its valid to say that we have $J\equiv 1$ which leads to $\hat{\mathbf{F}}\equiv \mathbf{F}$ and therefore leads to $\hat{\mathbf{C}}\equiv \mathbf{C}$ .

At large deformations, as an extension of the linear model, $N$ independent stress-free, volume-preserving intermediate configurations are postulated, to which the respective elastic material laws can be applied. For each of these configurations, the isochoric part of the deformation gradient decomposes multiplicatively into an elastic part ${\hat{\mathbf{F}}}_{e,i}$ and a viscoelastic part ${\mathbf{F}}_{v,i}$ , see Figure 2, or see for example, Ask et al. (2012),

$$\hat{\mathbf{F}}={\hat{\mathbf{F}}}_{e,i}·{\mathbf{F}}_{v,i},\mathrm{for}i\in \{1,\dots N\}.$$

(3)

Figure 2.

Representation of the multiplicative decomposition.

Above, the viscoelastic parts ${\mathbf{F}}_{v,i}$ of the deformation resemble the transformation to the stress-free $i^{\mathrm{th}}$ intermediate configuration. The tensor ${\mathbf{F}}_{v,i}$ is not necessarily compatible, which means there is no continuous displacement field generating a global intermediate configuration; the decomposition is purely local. However, it is assumed that the viscoelastic contributions are modeled in an isochoric way, such that $J_{v,i}=det{\mathbf{F}}_{v,i}\equiv 1$ for all $i=1\dots N$ . Based on the $i^{\mathrm{th}}$ elastic isochoric part of the deformation gradient ${\hat{\mathbf{F}}}_{e,i}=\hat{\mathbf{F}}·{\mathbf{F}}_{v,i}^{-1}$ , the corresponding elastic and viscoelastic Cauchy-Green tensors can be defined via,

$${\hat{\mathbf{C}}}_{e,i}={\hat{\mathbf{F}}}_{e,i}^{\mathbf{T}}·{\hat{\mathbf{F}}}_{e,i}={\mathbf{F}}_{v,i}^{-\mathbf{T}}·\hat{\mathbf{C}}·{\mathbf{F}}_{v,i}^{-1},$$

(4)

$${\mathbf{C}}_{v,i}={\mathbf{F}}_{v,i}^{\mathbf{T}}·{\mathbf{F}}_{v,i}.$$

(5)

2.2. Constitutive equations

In the following, a thermodynamically consistent model for a viscoelastic material shall be developed. To this end, the free energy density $\mathrm{\Psi }$ shall be defined for independent Cauchy-Green tensor $\mathbf{C}$ or deformation gradient $\mathbf{F}$ . The viscoelastic Cauchy-Green tensors ${\mathbf{C}}_{v,i}$ will serve as internal variables, and pressure $p$ will act as a Lagrangian multiplier for the incompressibility constraint.

The free energy is assumed to be additively composed of several independent potentials, compare Ask et al. (2012). We consider a decomposition into a hyperelastic potential ${\mathrm{\Psi }}_{\infty }$ and a viscoelastic energy densities ${\mathrm{\Psi }}_{v,i}$ , which shall be defined in the following,

$$\mathrm{\Psi }(\mathbf{C},{\mathbf{C}}_{v,i},p)={\mathrm{\Psi }}_{\infty }(\mathbf{C},p)+\sum _{i=1}^N{\mathrm{\Psi }}_{v,i}(\hat{\mathbf{C}},{\mathbf{C}}_{v,i}).$$

(6)

For the formulation of the constitutive equations we use the framework of the Clausius-Duhem inequality. This inequality states essentially that the dissipation $\mathcal{D}$ is non-negative (cf. Doghri, 2013), which can be written as,

$$\mathcal{D}=\frac{1}{2}\mathbf{T}:\stackrel{·}{\mathbf{C}}-\stackrel{·}{\mathrm{\Psi }}-s\stackrel{·}{T}-\mathrm{grad}T·\frac{\mathbf{q}}{T}\ge 0.$$

(7)

Above, $\mathbf{T}$ is the second Piola-Kirchhoff stress, $T$ is the absolute temperature, $s$ is the entropy, and $\mathit{q}$ is the heat flux vector. Since we assume isothermal processes in the following, with $\stackrel{·}{T}=0$ and no external heat flux $\mathit{q}=0$ , the expression simplifies to,

$$\mathcal{D}=\frac{1}{2}\mathbf{T}:\stackrel{·}{\mathbf{C}}-\stackrel{·}{\mathrm{\Psi }}\ge 0.$$

(8)

Computing the rate $\stackrel{·}{\mathrm{\Psi }}$ explicitly, equation (8) becomes,

$$\mathcal{D}=-(\frac{\partial \mathrm{\Psi }}{\partial \mathbf{C}}-\frac{1}{2}\mathbf{T}):\stackrel{·}{\mathbf{C}}-\sum _{i=1}^N\frac{\partial \mathrm{\Psi }}{\partial {\mathbf{C}}_{v,i}}:{\stackrel{·}{\mathbf{C}}}_{v,i}-\frac{\partial \mathrm{\Psi }}{\partial p}\stackrel{·}{p}\ge 0.$$

(9)

Equation (9) implies the classical constitutive equations for the second Piola-Kirchhoff stress $\mathbf{T}$ as,

$$\mathbf{T}=2\frac{\partial \mathrm{\Psi }}{\partial \mathbf{C}}.$$

(10)

Moreover, pressure $p$ acting as a Lagrangian multiplier ensures the condition of incompressibility as $\partial \mathrm{\Psi }/\partial p=0$ . The reduced dissipation inequality states that the negative inner product of generalized driving forces $\partial \mathrm{\Psi }/\partial {\mathbf{C}}_{v,i}$ and viscous strain rates ${\stackrel{·}{\mathbf{C}}}_{v,i}$ is non-negative,

$$\mathcal{D}=-\sum _{i=1}^N\frac{\partial \mathrm{\Psi }}{\partial {\mathbf{C}}_{v,i}}:{\stackrel{·}{\mathbf{C}}}_{v,i}\ge 0.$$

(11)

Let us collect specific representations for the energy densities constituting the augmented free energy $\mathrm{\Psi }$ . The elastic part of the free energy ${\mathrm{\Psi }}_{\infty }$ will be governed by an incompressible Neo-Hookean hyperelastic potential. We use ${\mu }_{\infty }=G_{\infty }$ for the long-term shear modulus,

$${\mathrm{\Psi }}_{\infty }(\mathbf{C},p)=\frac{{\mu }_{\infty }}{2}(\mathrm{tr}\mathbf{C}-3)-{\mu }_{\infty }\log J+p\log J.$$

(12)

For small strains, the Neo-Hookean material model reduces to a linear correlation, accurately described by Hooke’s law. Pressure $p$ serves as a Lagrangian multiplier to ensure the incompressibility constraint,

$$\frac{\partial {\mathrm{\Psi }}_{\infty }}{\partial p}=\log J=0,$$

such that $J\equiv 1$ .

Concerning the $i^{\mathrm{th}}$ multiplicative decomposition, a similar incompressible Neo-Hookean hyperelastic potential models the elastic energy for ${\hat{\mathbf{C}}}_{e,i}$ with shear modulus ${\mu }_{v,i}$ . As $\mathrm{tr}{\hat{\mathbf{C}}}_{e,i}=\hat{\mathbf{C}}:{\mathbf{C}}_{v,i}^{-1}$ , we use the representation,

$${\mathrm{\Psi }}_{v,i}(\hat{\mathbf{C}},{\mathbf{C}}_{v,i})=\frac{{\mu }_{v,i}}{2}(\hat{\mathbf{C}}:{\mathbf{C}}_{v,i}^{-1}-3).$$

(13)

Viscoelastic evolution is defined implicitly through a dissipation function $\mathrm{\Phi }$ which satisfies,

$$\sum _{i=1}^N\frac{\partial \mathrm{\Phi }}{\partial {\stackrel{·}{\mathbf{C}}}_{v,i}}:{\stackrel{·}{\mathbf{C}}}_{v,i}\ge 0.$$

(14)

As for the free energy densities, we also assume the dissipation function to be composed additively from contributions due to the $N$ internal variables ${\mathbf{C}}_{v,i}$ . Each of the terms is governed by a viscosity parameter ${\eta }_i$ . To be specific, we use,

$$\mathrm{\Phi }=\sum _{i=1}^N{\mathrm{\Phi }}_i({\mathbf{C}}_{v,i},{\stackrel{·}{\mathbf{C}}}_{v,i}),$$

(15)

with,

$${\mathrm{\Phi }}_i({\mathbf{C}}_{v,i},{\stackrel{·}{\mathbf{C}}}_{v,i})=\frac{1}{12}\sum _{i=1}^N{\eta }_i({\mathbf{C}}_{v,i}^{-1}·{\stackrel{·}{\mathbf{C}}}_{v,i}):({\mathbf{C}}_{v,i}^{-1}·{\stackrel{·}{\mathbf{C}}}_{v,i}).$$

(16)

With the choice of equations (16), (14) is satisfied. We justify this choice in the third section Incremental minimization principle by deriving a variational formulation and deducing the according flow-rule, equation (47). To analyze the implications of free energy $\mathrm{\Psi }$ and dissipation function $\mathrm{\Phi }$ on the evolution of ${\mathbf{C}}_{v,i}$ , an incremental variational formulation based on these quantities is posed later. There, we see that classical viscoelastic material models are reproduced for the above choices. Additionally, the incremental variational formulation can directly be transferred to a finite element discretization, which we do for the special case of rotationally symmetric problems.

3. Temperature-time shifting and parameter fitting of measured data

3.1. Linear viscoelasticity

The viscoelastic parameters such as the long-term modulus ${\mu }_{\infty }$ and the viscous properties ${\mu }_{v,i}$ and ${\eta }_i$ of our non-linear constitutive model are characterized from small-strain measurements, such that an identification using the linear Wiechert model is valid. Within this model, each spring is characterized through its elastic modulus $E_i$ . The long-term modulus of the system equals $E_{\infty }$ , that is, the modulus of the single spring. Its instantaneous modulus $E_0$ can be found as,

$$E_0=\sum _{i=1}^N{E}_i+E_{\infty },$$

(17)

for a thorough discussion we refer to Haupt (2013). For incompressible materials with a Poisson ratio of $\nu =0.5$ , shear moduli ${\mu }_i$ and Young’s moduli $E_i$ are connected via,

$${\mu }_{v,i}=\frac{E_i}{1+2\nu }=\frac{E_i}{3},i\in \{1,\dots N,\infty \}.$$

(18)

Additionally, ${\tau }_i$ defines the relaxation time of the $i^{\mathrm{th}}$ Maxwell element by,

$${\eta }_i={\tau }_i{E}_i.$$

(19)

At time $t>0$ , the relaxation modulus is found as, see Brinson and Brinson (2015),

$$E(t)=\sum _{i=1}^N{E}_i{e}^{-t/{\tau }_i}+E_{\infty }.$$

(20)

Representation (20) is often referred to as Prony-series. After Fourier transformation to frequency domain, one obtains the more common series expansion,

$$E^{*}(\omega )=\underset{\begin{array}{c}E\prime (\omega )\hfill \end{array}}{\underbrace{E_{\infty }+\sum _{i=1}^N\frac{E_i{\tau }_i^2{\omega }^2}{1+{\tau }_i^2{\omega }^2}}}+j\underset{\begin{array}{c}E″(\omega )\hfill \end{array}}{\underbrace{\sum _{i=1}^N\frac{E_i{\tau }_i\omega }{1+{\tau }_i^2{\omega }^2}}},$$

(21)

where $E^{*}(\omega )$ is termed the complex modulus of the material. In this context, the real part $\Re (E^{*}(\omega ))$ is called storage modulus $E\prime (\omega )$ and the imaginary part $\Im (E^{*}(\omega ))$ is known as the loss modulus $E″(\omega )$ , see Tschoegl (2012).

Furthermore we define $e_i=E_i/E_0$ with the instantaneous modulus $E_0$ as defined in (17). Following our model, the material is now fully characterized when estimating the $2N$ parameters $e_i$ and ${\tau }_i$ as well as the long-term modulus $E_{\infty }$ . We often refer to $N$ as the order of the Prony-series.

3.2. Dynamical temperature mechanical analysis

To determine these Prony-parameters from experiments, dynamical temperature mechanical analysis (DTMA) measurements have been performed. These DTMA measurements provide storage and loss modulus over a given frequency range at different temperatures. The measurements were carried out on a Netzsch DMA GABO Eplexor^®500N HT.

The specimen material, Sylgard 184, was prepared using a $10:1$ (base:crosslinker) mixing ratio by weight. The mixture was stirred for 5 min and then cast into a mold and cured at $75°\mathrm{C}$ for 24 h. To create the DTMA samples, a square mold was used, and then cut to size. The cylindrical samples used for the ball-drop tests were produced using a corresponding circular mold. The production of samples is described in detail in Preuer (2018). The dimensions of the sample were $10\mathrm{mm}\times 40\mathrm{mm}\times 3.2\mathrm{mm}$ $(W\times H\times D)$ and a clamping length of 20 mm was chosen. The static strain was 3%, the dynamic strain 2% with a control range of 0.0045% and the pre-stress was 0.05 MPa. The pre-cooling time in the temperature chamber was at least 2 min. Measurements were taken in the range $0°\mathrm{C}$ to $40°\mathrm{C}$ in $2.5°\mathrm{C}$ steps between 0 and 100 Hz. Each measurement contains 12 cycles, from which 10 are used to compute the average. Due to resonance peaks for all curves at around 72 Hz we did not use the entire measurement curve, but only from 0 Hz to 60 Hz, as shown in Figure 3(a) and (b).

Figure 3.

DTMA measured data: (a) measured storage modulus $E\prime$ , (b) measured loss modulus $E″$ .

This range is too small to describe the whole range of frequencies excited in applications. The time-temperature superposition principle (TTSP) allows to generate a so-called master curve that can cover a much larger frequency spectrum from measurements at different temperatures. Conversely, predictions can also be made about the long-term behavior, see also Brinson and Brinson (2015), Ferry (1980).

3.3. Time-temperature superposition

Concerning time-temperature superposition, Sylgard 184 is treated as a thermorheologically complex material, cf. Klompen and Govaert (1999). In the modeled area, the cross-linked material is well above the glass transition temperature $T_g$ and the melting temperature $T_m$ , cf. Dewimille et al. (2005). This results in a melt-like mobility, which can be described by means of thermorheologically complex modeling. The selected approach combines the classic TTSP for solid polymers (complex Young’s modulus) and the TTSP for polymer melts (complex viscosity). Thus, horizontal frequency shifts as well as vertical shifts in the modulus are applied. In Tschoegl et al. (2002), the following general form of time-temperature shifts was identified as,

$$E_T^{*}(\omega )=b_T{E}_{T_0}^{*}(a_T\omega ).$$

(22)

Shifting parameters $a_T$ and $b_T$ allow to transfer the material’s modulus from reference temperature $T_0$ to any other temperature $T$ , at least within a certain temperature range. A simple approach for shifting the modulus is given by,

$$b_T=\frac{T_0{\rho }_0}{T\rho },$$

(23)

where $\rho$ and ${\rho }_0$ are the mass density of the polymer at $T$ and $T_0$ . Since we assume no relevant density differences within the sampled temperature range of $0°\mathrm{C}$ – $40°\mathrm{C}$ , this density ratio is neglected in the following.

For the horizontal frequency shift a large variety of models can be found. Probably the one most commonly used is the Williams-Landel-Ferry model (WLF). However, a main criterion for this is free-volume availability, which is despite the amorphous morphology at ambient temperature limited by the intermolecular bonds created in cross linking. As shown in Lomellini (1992), the WLF model can be used very well for a temperature range of $T_g<T<T_g+100°\mathrm{C}$ .

Clarson et al. (1985) and also Li et al. (2006) found that the glass transition temperature of Sylgard 184 is in the range of $-150°\mathrm{C}<T_g<-120°\mathrm{C}$ and the melting temperature $T_m$ lies between $-50°\mathrm{C}$ and $-70°\mathrm{C}$ . This implies that the WLF model is not applicable to Sylgard 184 at ambient temperatures. Lomellini (1992) states that in this relatively high temperature range of $T>T_g+100$ °C the so-called Arrhenius equation is applicable,

$$a_T=e^{\alpha (\frac{1}{T}-\frac{1}{T_0})}.$$

(24)

The parameter $\alpha =E_a/\mathcal{R}$ is independent of temperature, but can be computed from the activation energy $E_0$ and the universal gas constant $\mathcal{R}$ . We use the DTMA measurement data to compute $\alpha$ , and thereby the activation energy $E_a$ , using a least-squares method.

In order to estimate $\alpha$ , at reference temperature of $T_0=20°\mathrm{C}$ measurements within an extended frequency range were conducted. While it is of course not possible to extend this range towards higher frequencies, measurements for lower frequencies approaching zero could easily be conducted; the results are displayed in Figure 4. Then, data sets for temperatures from $20°\mathrm{C}$ to $40°\mathrm{C}$ were used to determine the unknown shift parameter $\alpha$ , as higher temperatures correspond to lower frequencies. A least squares method is used to find the optimal parameter $\alpha$ in order to minimize the squared vertical distance of the respective data curves when shifted by $b_T$ and $a_T$ as in (23) and (24), respectively. As the shifted frequencies for the data points do not coincide in general, a linear interpolation of the reference data $E_{T_0}^{*}$ in between data points is used instead. For ${\omega }_i,i=1\dots n_{\omega }$ the frequencies at which DTMA data is available, the least-squares problem finally reads,

$$\sum _{T>T_0}\sum _{i=1}^{n_{\omega }}{|b_T{E}_{T_0}^{*}(a_T{\omega }_i)-E_T^{*}({\omega }_i)|}^2\to \underset{\alpha >0}{min}.$$

(25)

Figure 4.

Master curves for storage and loss modulus at $T_0=20°\mathrm{C}$ gained through time-temperature superposition: (a) storage modulus $E_{T_0}\prime$ at $T_0=20°\mathrm{C}$ , (b) loss modulus $E_{T_0}″$ at $T_0=20°\mathrm{C}$ .

For the available DTMA data, this procedure results in $\alpha =7389.124K$ . Then, both the storage and loss moduli from all DTMA data temperatures can be shifted according to (22) to get a master curve at temperature $T_0=20°\mathrm{C}$ as shown in Figure 4. These shifted data sets are used in the subsequent detection of Prony-parameters; the additional data in the extended frequency range for temperature $T_0$ are not included in this process anymore.

3.4. Prony-parameters

Once data in an extended frequency range has been generated using the TTSP, Prony-parameters $E_{\infty }$ as well as $e_i$ and ${\tau }_i$ can be detected. To this end, a non-linear least-squares optimization problem for the $2N+1$ unknown Prony-parameters needs to be solved. For the presented data, the non-negative least squares based routine curve_fit provided within python’s scipy package could be applied successfully.

This function works with complex data sets, and thereby ensures that the parameters determined are accurate for both the storage and loss modulus. Moreover, range conditions have to be met for the different parameters. A minimal set of these conditions is to seek for

• a positive long-term modulus $E_{\infty }>0$ ,

• positive relaxation times ${\tau }_i>0$ ,

• relative viscoelastic moduli $0<e_i<1$ ,

• satisfying additionally $\sum _{i=1}^N{e}_i\le 1$ .

Simple range conditions like the former three items can be included in the parameters passed to the least-squares fitting routine. The latter condition has to be verified manually, once the parameters are computed. As one often observes in non-linear problems, choosing a feasible starting value is of utmost importance. We observed that setting all values to one worked well when seeking the long-term modulus equivalent to 1 MPa.

Finally, the order $N$ of the Wiechert model has to be chosen in accordance with the available data set. As a rule of thumb, good approximation is obtained when using one Maxwell element per measured frequency decade. Since after TTSP, measurements are available for frequencies from 0.1 Hz up to almost 1000 Hz (see Figure 4), fitting Prony-parameters for models of third and fourth order seems appropriate. The respective results are depicted in Figure 5. Tables 1 and 2 show the corresponding Prony-parameters $E_{\infty }$ , $E_i=e_i{E}_0$ and ${\tau }_i$ as well as the shear moduli ${\mu }_i$ and viscosities ${\eta }_i$ computed from these values via (18) and (19).

Figure 5.

Fitted Prony-series, comparison $N=3$ and $N=4$ .

Table 1.

Prony-parameters for $N=3$ and $T_0=20°\mathrm{C}$ .

$i$	$E_i$ in MPa	${\tau }_i$ in s	${\mu }_i$ in MPa	${\eta }_i$ in MPa s
$\infty$	2.11904	–	0.70635	–
1	0.93789	0.00037	0.31263	0.00035
2	0.34398	0.00424	0.11466	0.00146
3	0.24635	0.04684	0.08212	0.01154

Table 2.

Prony-parameters for $N=4$ and $T_0=20°\mathrm{C}$ .

$i$	$E_i$ in MPa	${\tau }_i$ in s	${\mu }_{v,i}$ in MPa	${\eta }_i$ in MPa s
$\infty$	2.07639	–	0.69213	–
1	1.06895	0.00023	0.35632	0.00024
2	0.34258	0.00198	0.11419	0.00068
3	0.24093	0.01323	0.08031	0.00319
4	0.17836	0.1386	0.05945	0.02472

In both Tables 1 and 2, we note that the relaxation times ${\tau }_i$ are distributed more or less evenly across the measured decades, as is to be expected due to Tschoegl (2012).

As parameters at reference temperature $T_0=20°\mathrm{C}$ are determined, according parameters for different ambient temperatures $T$ can be computed through the time-temperature superposition principle (22), leading to

$$\begin{array}{cc}\multicolumn{1}{c}{{\mu }_{v,i}(T)}& =b_T{\mu }_{v,i}(T_0),\hfill \\ \multicolumn{1}{c}{{\tau }_i(T)}& =a_T{\tau }_i(T_0),\hfill \\ \multicolumn{1}{c}{{\eta }_i(T)}& =3{\mu }_{v,i}(T){\tau }_i(T).\hfill \end{array}$$

(26)

Parameters $a_T$ and $b_T$ are defined in (24) and (23), values for specific temperatures are listed in Table 4.

4. Incremental minimization principle

Having determined the relevant constitutive parameters, we now seek a finite element formulation for the problem at hand. As a starting point, a variational formulation for viscoelastic solids is derived from an incremental minimization problem. A similar procedure has been presented for modeling various dissipative processes by the group around Miehe, we cite non-exhaustively Miehe (2002); Miehe et al. (2002, 2011). For reasons of simplicity of presentation, we restrict ourselves to the case of a conservative volume force density $\mathit{f}$ . With ${\rho }_0$ the mass density, the kinetic energy density $T$ and local power of external forces ${\mathcal{P}}_{\mathrm{ext}}$ read,

$$T=\frac{1}{2}{\rho }_0\stackrel{·}{\mathit{u}}·\stackrel{·}{\mathit{u}},{\mathcal{P}}_{\mathrm{ext}}=\mathit{f}·\stackrel{·}{\mathit{u}}.$$

(27)

Consider a finite time step $[t_0,t_0+\mathrm{\Delta }t]$ with known initial conditions at time $t_0$ . We introduce the incremental potential for this time step as,

$$\begin{array}{c}\multicolumn{1}{c}{{\mathrm{\Pi }}_t^{t+\mathrm{\Delta }t}={\int }_{t_0}^{t_0+\mathrm{\Delta }t}(\stackrel{·}{T}+\stackrel{·}{\mathrm{\Psi }}+\mathrm{\Phi }-{\mathcal{P}}_{\mathrm{ext}})\mathrm{dt}}\\ =\underset{\mathrm{kinetic}}{\underbrace{T_{t_0+\mathrm{\Delta }t}-T_{t_0}}}+\underset{\mathrm{stored}}{\underbrace{{\mathrm{\Psi }}_{t_0+\mathrm{\Delta }t}-{\mathrm{\Psi }}_{t_0}}}+\underset{\mathrm{dissipated}}{\underbrace{{\int }_{t_0}^{t_0+\mathrm{\Delta }t}\mathrm{\Phi }\mathrm{dt}}}\hfill \\ -\underset{\mathrm{external}}{\underbrace{{\int }_{t_0}^{t_0+\mathrm{\Delta }t}{\mathcal{P}}_{\mathrm{ext}}\mathrm{dt}}}.\hfill \end{array}$$

(28)

The incremental potential ${\mathrm{\Pi }}_{t_0}^{t_0+\mathrm{\Delta }t}$ is then minimized with respect to all admissible displacements and viscoelastic internal parameters, and maximized with respect to hydrostatic pressure. For time step size $\mathrm{\Delta }t$ approaching zero, we arrive at,

$$\stackrel{·}{T}+\stackrel{·}{\mathrm{\Psi }}+\mathrm{\Phi }-{\mathcal{P}}_{\mathrm{ext}}\to \underset{\stackrel{·}{\mathbf{u}}}{min}\underset{{\stackrel{·}{\mathbf{C}}}_{v,i}}{min}\underset{\stackrel{·}{p}}{max}.$$

(29)

We can now rewrite equation (29) for the specific choice of free energy $\mathrm{\Psi }$ and get,

$$\begin{array}{c}\rho \stackrel{··}{\mathit{u}}·\stackrel{·}{\mathit{u}}+\frac{\partial \mathrm{\Psi }}{\partial \mathbf{C}}:\stackrel{·}{\mathbf{C}}+\sum _{i=1}^N\frac{\partial \mathrm{\Psi }}{\partial {\mathbf{C}}_{v,i}}:{\stackrel{·}{\mathbf{C}}}_{v,i}+\frac{\partial \mathrm{\Psi }}{\partial p}\stackrel{·}{p}\hfill \\ +\mathrm{\Phi }-\mathit{f}·\stackrel{·}{\mathit{u}}\to \underset{\stackrel{·}{\mathit{u}}}{min}\underset{{\stackrel{·}{\mathbf{C}}}_{v,i}}{min}\underset{\stackrel{·}{p}}{max}.\hfill \end{array}$$

(30)

Under common assumptions on convexity and smoothness of the involved energies and dissipation functions, the optimization problem is equivalent to finding critical points as solutions to the corresponding variational equation, also known as D’Alembert’s principle,

$$\begin{array}{c}{\int }_{\mathcal{B}}\rho \stackrel{··}{\mathit{u}}·\delta \mathit{u}\mathrm{dV}\hfill \\ +{\int }_{\mathcal{B}}(\frac{\partial \mathrm{\Psi }}{\partial \mathbf{C}}:\delta \mathbf{C}+\sum _{i=1}^N\frac{\partial \mathrm{\Psi }}{\partial {\mathbf{C}}_{v,i}}:\delta {\mathbf{C}}_{v,i}+\frac{\partial \mathrm{\Psi }}{\partial p}\delta p)\mathrm{dV}\hfill \\ +{\int }_{\mathcal{B}}\sum _{i=1}^N\frac{\partial \mathrm{\Phi }}{\partial {\stackrel{·}{\mathbf{C}}}_{v,i}}:\delta {\mathbf{C}}_{v,i}\mathrm{dV}-{\int }_{\mathcal{B}}\mathit{f}·\delta \mathit{u}\mathrm{dV}=0.\hfill \end{array}$$

(31)

For a simplified presentation, we rewrite D’Alembert’s principle in more compact form, although abusing notation concerning the variation of the dissipation function,

$${\int }_{\mathcal{B}}(\rho \stackrel{··}{\mathit{u}}·\delta \mathit{u}+\delta \mathrm{\Psi }+\delta \mathrm{\Phi }-\mathit{f}·\delta \mathit{u})\mathrm{dV}=0,$$

(32)

with,

$$\delta \mathrm{\Psi }:=\frac{\partial \mathrm{\Psi }}{\partial \mathbf{C}}:\delta \mathbf{C}+\sum _{i=1}^N\frac{\partial \mathrm{\Psi }}{\partial {\mathbf{C}}_{v,i}}:\delta {\mathbf{C}}_{v,i}+\frac{\partial \mathrm{\Psi }}{\partial p}\delta p,$$

(33)

and

$$\delta \mathrm{\Phi }:=\sum _{i=1}^N\frac{\partial \mathrm{\Phi }}{\partial {\stackrel{·}{\mathbf{C}}}_{v,i}}:\delta {\mathbf{C}}_{v,i}.$$

(34)

The variational equation (31) implicitly defines the evolution law for ${\stackrel{·}{\mathbf{C}}}_{v,i}$ for arbitrary but fixed index $i\in \{1,\dots N\}$ . We show that, for the proposed choice of dissipation function from equation (16), the standard model as for example, proposed by Haupt (2013) is re-established. We collect all terms containing $\delta {\mathbf{C}}_{v,i}$ from equation (31), and obtain,

$${\int }_{\mathcal{B}}(\frac{\partial {\mathrm{\Psi }}_{v,i}}{\partial {\mathbf{C}}_{v,i}}+\frac{1}{6}{\eta }_i({\mathbf{C}}_{v,i}^{-1}·{\stackrel{·}{\mathbf{C}}}_{v,i}·{\mathbf{C}}_{v,i}^{-1})):\delta {\mathbf{C}}_{v,i}\mathrm{dV}=0.$$

(35)

Since $det{\mathbf{C}}_{v,i}=1$ , we observe for its rate and variation,

$$\frac{d}{\mathrm{dt}}(det{\mathbf{C}}_{v,i})={\mathbf{C}}_{v,i}^{-1}:{\stackrel{·}{\mathbf{C}}}_{v,i}=0,$$

(36)

$$\delta (det{\mathbf{C}}_{v,i})={\mathbf{C}}_{v,i}^{-1}:\delta {\mathbf{C}}_{v,i}=0.$$

(37)

Thus, a contribution of the form $\gamma {{\mathbf{C}}_{v,i}}^{-1}:\delta {\mathbf{C}}_{v,i}=0$ can be added without changing the equation (35), which leads to,

$${\int }_{\mathcal{B}}(\frac{\partial {\mathrm{\Psi }}_{v,i}}{\partial {\mathbf{C}}_{v,i}}+\gamma {\mathbf{C}}_{v,i}^{-1}+\frac{1}{6}{\eta }_i{\mathbf{C}}_{v,i}^{-1}·{\stackrel{·}{\mathbf{C}}}_{v,i}·{\mathbf{C}}_{v,i}^{-1}):\delta {\mathbf{C}}_{v,i}\mathrm{dV}=0.$$

(38)

Here, $\gamma$ is an arbitrary constant that still has to be determined in the following. As $\delta {\mathbf{C}}_{v,i}$ is otherwise arbitrary, the equality holds in classical sense,

$$\frac{\partial {\mathrm{\Psi }}_{v,i}}{\partial {\mathbf{C}}_{v,i}}+\gamma {\mathbf{C}}_{v,i}^{-1}+\frac{1}{6}{\eta }_i{\mathbf{C}}_{v,i}^{-1}·{\stackrel{·}{\mathbf{C}}}_{v,i}·{\mathbf{C}}_{v,i}^{-1}=0,$$

(39)

which directly implies,

$${\stackrel{·}{\mathbf{C}}}_{v,i}=-\frac{6}{{\eta }_i}{\mathbf{C}}_{v,i}·(\frac{\partial {\mathrm{\Psi }}_{v,i}}{\partial {\mathbf{C}}_{v,i}}+\gamma {\mathbf{C}}_{v,i}^{-1})·{\mathbf{C}}_{v,i}.$$

(40)

Contraction with ${\mathbf{C}}_{v,i}^{-1}$ leads through equation (36) to,

$$\begin{array}{c}\multicolumn{1}{c}{0={\mathbf{C}}_{v,i}^{-1}:[{\mathbf{C}}_{v,i}·(\frac{\partial {\mathrm{\Psi }}_{v,i}}{\partial {\mathbf{C}}_{v,i}}+\gamma {\mathbf{C}}_{v,i}^{-1})·{\mathbf{C}}_{v,i}]}\\ =\mathbf{I}:[({\mathbf{C}}_{v,i}·\frac{\partial {\mathrm{\Psi }}_{v,i}}{\partial {\mathbf{C}}_{v,i}}+\gamma \mathbf{I})·{\mathbf{C}}_{v,i}·{\mathbf{C}}_{v,i}^{-1}].\hfill \end{array}$$

(41)

Since ${\mathbf{C}}_{v,i}·{\mathbf{C}}_{v,i}^{-1}=\mathbf{I}$ , we can simplify even further,

$$\begin{array}{c}\multicolumn{1}{c}{0={\mathbf{C}}_{v,i}·\frac{\partial {\mathrm{\Psi }}_{v,i}}{\partial {\mathbf{C}}_{v,i}}:\mathbf{I}+\gamma \mathbf{I}:\mathbf{I}}\\ =-{\mathbf{M}}_{v,i}:\mathbf{I}+3\gamma .\hfill \end{array}$$

(42)

Above, we introduce ${\mathbf{M}}_{v,i}$ as the $i^{\mathrm{th}}$ Mandel-type stress tensor,

$${\mathbf{M}}_{v,i}=-{\mathbf{C}}_{v,i}·\frac{\partial {\mathrm{\Psi }}_{v,i}}{\partial {\mathbf{C}}_{v,i}}.$$

(43)

With this, we finally have an expression to determine the unknown parameter $\gamma$ ,

$$\gamma =\frac{1}{3}\mathrm{tr}{\mathbf{M}}_{v,i}.$$

(44)

We can now substitute $\gamma$ back into equation (40) and after a few algebraic manipulations we get,

$${\stackrel{·}{\mathbf{C}}}_{v,i}=\frac{6}{{\eta }_i}\mathrm{dev}({\mathbf{M}}_{v,i})·{\mathbf{C}}_{v,i},$$

(45)

which is stated in similar form by Ask et al. (2012) for viscoelastic electrostrictive polymers.

To compare the final evolution law to standard theory, we explicitly compute the Mandel-type stresses for the viscoelastic potential ${\mathrm{\Psi }}_{v,i}$ from equation (13) and obtain,

$${\mathbf{M}}_{v,i}=-\frac{{\mu }_{v,i}}{2}(\hat{\mathbf{C}}·{\mathbf{C}}_{v,i}^{-1}).$$

(46)

After combining the two equations (45) and (46) we finally obtain the evolution law for ${\mathbf{C}}_{v,i}$ as stated in Haupt (2013),

$$\frac{2}{3}\eta {\stackrel{·}{\mathbf{C}}}_{v,i}=2{\mu }_{v,i}(\hat{\mathbf{C}}-\frac{1}{3}\mathrm{tr}(\hat{\mathbf{C}}·{\mathbf{C}}_{v,i}^{-1}){\mathbf{C}}_{v,i}).$$

(47)

Note that, in the finite element approach proposed below, this evolution equation is not integrated explicitely, but is implied implicitely through the incremental variational formulation including variations $\delta {\mathbf{C}}_{v,i}$ .

5. Finite elements for radially symmetric problems

For radially symmetric problems, a finite element model respecting incompressibility of the viscous strains explicitely can be devised as follows. This setup is designed such that it allows for verification against various measurement setups, like the ball-drop experiment in the last section to determine rebound resilience.

5.1. Spatial discretization for radially symmetric problems

Let the body $\mathcal{B}$ as well as all external forces and kinematic constraints be in accordance with radial symmetry. Without restriction of generality, we consider ${\mathbf{e}}_z$ as axis of rotation. Deformation gradient $\mathbf{F}$ and right Cauchy-Green $\mathbf{C}$ are defined in the standard way, assuming that the solution is independent of the angular coordinate $\phi$ .

Assume $\mathcal{T}=\{T\}$ to be a triangular finite element mesh of the section of $\mathcal{B}$ under consideration. Concerning the elastic problem, we propose to use the pair of displacement-pressure elements as proposed by Taylor and Hood (cf. Taylor and Hood 1973), where both $\mathit{u}$ and $p$ are represented by continuous basis functions ${\mathbf{N}}_i^{(u)}$ and $N_i^{(p)}$ respectively, and the approximation order for $\mathit{u}$ is one higher than for $p$ . Additionally, kinematic boundary conditions $\mathit{u}=0$ on the boundary part ${\Gamma }_{\mathrm{fix}}$ have to be satisfied. Summing up, finite element approximations for $\mathit{u}$ and $p$ satisfy,

$$\mathit{u}\in \{\mathit{v}:\mathit{v}{|}_T\in (P^k(T){)}^2,\mathit{v}\mathrm{cont}.,\mathit{v}{|}_{{\Gamma }_{\mathrm{fix}}}=0\},$$

(48)

$$\mathit{u}=\sum _{i=1}^{n_u}{q}_i^{(u)}{\mathbf{N}}_i^{(u)},$$

(49)

$$p\in \{q:q{|}_T\in P^{k-1}(T),q\mathrm{cont}.\},$$

(50)

$$p=\sum _{i=1}^{n_p}{q}_i^{(p)}{N}_i^{(p)},$$

(51)

with $P^k(T)$ meaning the space of bivariate polynomials of order $k$ on triangular elements $T\in \mathcal{T}$ .

In addition, it is necessary to provide a finite-dimensional approximation for the viscoelastic internal variables ${\mathbf{C}}_{v,i}$ . Here, it is important that the kinematic constraint of incompressibility is satisfied by construction. For a rotationally symmetric setting, this can be achieved by a non-linear representation of ${\mathbf{C}}_{v,i}$ in terms of the independent viscoelastic strains $s_{i,\mathrm{rr}}$ , $s_{i,\mathrm{zz}}$ and $s_{i,\mathrm{rz}}$ . The $i^{\mathrm{th}}$ viscoelastic Cauchy-Green tensor respects the axisymmetric structure,

$${\mathbf{C}}_{v,i}=\left[\begin{array}{ccc}2s_{i,\mathrm{rr}}+1\hfill & 0\hfill & 2s_{i,\mathrm{rz}}\hfill \\ \multicolumn{1}{c}{0}& c_{i,\phi \phi }\hfill & 0\hfill \\ \multicolumn{1}{c}{2s_{i,\mathrm{rz}}}& 0\hfill & 2s_{i,\mathrm{zz}}+1\hfill \end{array}\right].$$

(52)

Component $c_{i,\phi \phi }$ is chosen such that the incompressibility constraint $det{\mathbf{C}}_{v,i}\equiv 1$ is satisfied a-priorily, which implies,

$$c_{i,\phi \phi }={((2s_{i,\mathrm{rr}}+1)(1+2s_{i,\mathrm{zz}})-4s_{i,\mathrm{rz}})}^{-1}.$$

(53)

The viscoelastic strains $s_{i,\mathrm{rr}},s_{i,\mathrm{zz}}$ and $s_{i,\mathrm{rz}}$ are defined element-wise polynomial, without continuity, one polynomial order lower than the displacements,

$$\begin{array}{l}{s}_{i,rr},s_{i,zz},s_{i,rz}\in \{s:s{|}_T\in P^{k-1}(T)\},\\ s_{i,\alpha \beta }=\sum _{i=1}^{n_s}{q}_i^{(s,\alpha \beta )}{N}_i^{(s)}.\end{array}$$

(54)

5.2. Time integration

Introducing the finite element discretizations equations (50), (51), and (54) in D’Alembert’s principle (31) leads to a system of ordinary differential equations for the degree of freedom vector $\mathit{q}$ . As the physical problem includes accelerations $\stackrel{··}{\mathit{u}}$ and viscoelastic strain rates ${\stackrel{·}{\mathbf{C}}}_{v,i}$ as well as the Lagrangian multiplier $p$ , a differential-algebraic system of equations is generated after spatial discretization. It is well known that the Newmark- $\beta$ scheme is well-suited to solve differential-algebraic equations of index 2, where the choice of parameter $\beta$ further determines its numerical damping properties.

We shortly present the Newmark- $\beta$ scheme adapted for the viscoelastic problem at hand. To this end, we assume that velocities $\mathit{v}=\stackrel{·}{\mathit{u}}$ are approximated by an independent linear combination of displacement basis functions ${\mathbf{N}}_i^{(u)}$ , via additional degrees of freedom $q_i^{(v)}$ . For any time step $[t_n,t_{n+1}=t_n+\mathrm{\Delta }t]$ , the quantities ${\mathit{u}}_n,{\mathit{v}}_n,p_n$ as well as ${\mathbf{C}}_{v,i,n}$ , ${\stackrel{·}{\mathbf{C}}}_{v,i,n}$ (or rather the viscoelastic strains $s_{i,\alpha \beta ,n}$ and strain rates ${\stackrel{·}{s}}_{i,\alpha \beta ,n}$ ) at time $t_n$ are assumed to be available from the previous step. According updates $\mathit{u},\mathit{v},p,s_{i,\alpha \beta },$ and ${\stackrel{·}{s}}_{i,\alpha \beta }$ at time $t_{n+1}=t_n+\mathrm{\Delta }t$ are to be computed, where we omit $n+1$ as an index in favor of brevity of presentation. The Newmark- $\beta$ scheme interpreted in terms of the underlying variational equations leads to the coupled non-linear system, where we use $\delta \mathit{v}$ as the (independent) variation of the velocity $\mathit{v}$ ,

$$\begin{array}{c}{\int }_{\mathcal{B}}\rho (\mathit{u}-{\mathit{u}}_n-\mathrm{\Delta }t{\mathit{v}}_n)·\delta \mathit{u}\mathrm{dV}\hfill \\ +{\left(\mathrm{\Delta }t\right)}^2{\int }_{\mathcal{B}}\beta (\delta \mathrm{\Psi }+\delta \mathrm{\Phi }-\mathit{f}·\delta \mathit{u})\mathrm{dV}\hfill \\ +{\left(\mathrm{\Delta }t\right)}^2{\int }_{\mathcal{B}}(\frac{1}{2}-\beta )(\delta {\mathrm{\Psi }}_n+\delta {\mathrm{\Phi }}_n-{\mathit{f}}_n·\delta \mathit{u})\mathrm{dV}=0,\hfill \end{array}$$

(55)

$$\begin{array}{c}{\int }_{\mathcal{B}}\rho (\mathit{v}-{\mathit{v}}_n)·\delta \mathit{v}\mathrm{dV}\hfill \\ +\mathrm{\Delta }t{\int }_{\mathcal{B}}\gamma ({\frac{\partial \mathrm{\Psi }}{\partial \mathit{u}}|}_{t_{n+1}}·\delta \mathit{v}-\mathit{f}·\delta \mathit{v})\mathrm{dV}\hfill \\ +\mathrm{\Delta }t{\int }_{\mathcal{B}}(1-\gamma )({\frac{\partial \mathrm{\Psi }}{\partial \mathit{u}}|}_{t_n}·\delta \mathit{v}-{\mathit{f}}_n·\delta \mathit{v})\mathrm{dV}=0.\hfill \end{array}$$

(56)

Note that equation (55) does not only implement the Newmark integration scheme for the displacements, but also the incompressibility constraint via $\partial \mathrm{\Psi }/\partial p·\delta p=0$ and the evolution law for the viscoelastic strains as indicated by the dissipation function. The variation of the dissipation function $\delta \mathrm{\Phi }$ is to be understood in the same manner as in equation (32).

6. Validation

To validate the material characterization developed above, we present a comparison of computational results to measurements for a ball-drop test setup. Ball-drop tests are frequently used to determine the rebound resilience of materials by comparing the original drop height to the rebound height of the ball. When varying temperature, drop height as well as size and mass of the ball, changes in the rebound resilience are of interest.

For the purpose of computational realization we use the open-source software package Netgen/NGSolve ² . Via a python interface, the developed variational time stepping scheme (55)–(56) can be defined symbolically; all proposed finite element spaces are available at arbitrary order. Automatic differentiation takes care of computing correct tangent stiffnesses in a straightforward manner.

As sketched in Figure 6, we consider a cylindric specimen made from Sylgard 184 with a radius of $r_S=30\mathrm{mm}$ and a height of $h_S=30\mathrm{mm}$ . A steel ball of radius $r_B=15\mathrm{mm}$ and mass $m_B=109\mathrm{g}$ is dropped from an initial height of $h_0=450\mathrm{mm}$ above the cylinder’s surface. With these parameters, large local deformations exceeding 30% can be seen in the later results. All parameters are once again summarized in Table 3.

Figure 6.

Computational setup of the ball-drop test.

Table 3.

Parameters of the ball-drop test.

Specimen radius	$r_S$	$30\mathrm{mm}$
Specimen height	$h_S$	$30\mathrm{mm}$
Specimen density	${\rho }_{\mathrm{Syl}}$	$965\mathrm{kg}{\mathrm{m}}^{-3}$
Drop height	$h_0$	$450\mathrm{mm}$
Ball radius	$r_B$	$15\mathrm{mm}$
Ball mass	$m_B$	$109\mathrm{g}$

A finite element model exploiting the model’s radial symmetry is generated, where the finite element mesh consists of $223$ triangular elements. Third order elements are used in simulations, which results in a total of $2212$ displacement degrees of freedom. While the proposed material characterization is employed to describe the cylindric specimen, the ball is modeled as linear elastic body with Young’s modulus $E=210\mathrm{GPa}$ and Poisson ratio $\nu =0.3$ . We choose a penalty-formulation for frictionless contact to describe the interaction of ball and PDMS. In experiments, the polymer’s surface is powdered to reduce friction and adhesion effects to a minimum.

Assume $t_0=0\mathrm{s}$ to denote the time at which the ball is released. The actual simulation is restricted to the time interval around the ball’s first impact. To be concise, time stepping is initialized for the ball at height $h_1=0.02\mathrm{m}$ , which is close to the specimen’s surface. This results in an initial time $t_1$ an initial velocity ${\mathit{v}}_1=v_1{\mathit{e}}_z$ set to,

$$t_1=\sqrt{\frac{2}{g}(h_0-h_1)},v_1=-g{t}_1.$$

(57)

In Figure 7, these time instants as well as the corresponding velocities and observed heights are visualized.

Figure 7.

Idealized trajectory of the ball.

In simulations, the Newmark parameters are chosen as $\beta =0.25$ and $\gamma =0.5$ to avoid numerical damping. The time step size is set to to $\mathrm{\Delta }t=1\times {10}^{-4}\mathrm{s}$ to generate accurate results when using the second-order accurate Newmark- $\beta$ method. The simulation is stopped once the ball exceeds the initial distance of $h_1$ , let this time instant be denoted as $t_2$ . The first rebound height $h_{r,1}$ can then be computed from height $h_2$ and vertical velocity $v_2$ at time $t_2$ ,

$$h_{r,1}=h_2+\frac{g}{2}{v}_2^2.$$

(58)

In Figure 8, the ball’s computed trajectory at different temperatures is presented, where the third order material characterization from Table 1 has been used. Five different ambient temperatures were chosen such that they cover the range of DTMA measurements. We use $10,20,30,40\mathrm{and}50°\mathrm{C}$ respectively, leading to the shifting parameters provided in Table 4.

Figure 8.

Ball-drop trajectory, $h_0=0.45\mathrm{m},r_B=15\mathrm{mm}$ , using material data for $N=3$ from Table 1.

Table 4.

Shifting parameters for $T_0=20°\mathrm{C}$ and different ambient temperatures $T$ .

$T$ in $°\mathrm{C}$	10	20	30	40	50
$a_T$	2.4388	1.0	0.4387	0.1994	0.0959
$b_T$	1.0353	1.0	0.967	0.9361	0.9071

Free fall and rebound trajectories above height $h_1$ and $h_2$ respectively were generated analytically and are identified as dashed lines. Results of the actual simulation are highlighted in blue, and presented in detail in Figure 9. In addition to the trajectories, the ball’s vertical velocity is displayed for all temperatures. An increase in the rebound height for increasing temperatures can be clearly observed.

Figure 9.

Impact trajectory (cmp. Figure 8), solid lines indicate the ball’s trajectory, dots mark its velocity.

A closer look at the ball trajectory in Figure 9 shows that even for $T=10°\mathrm{C}$ representing the stiffest sample, the indentation amounts to approximately 4.71 mm. The resulting compressive strain then corresponds to mean value of about 15%, and a peak value of around 35%. This clearly indicates large deformations, while an approximation by using the Neo-Hookean law is still justifiable. Figure 10 shows the according strain distribution at the time of maximum penetration of the ball. We visualize $E_{\mathrm{zz}}$ , being the compressive component of the strain tensor $\mathbf{E}=\frac{1}{2}(\mathbf{C}-\mathbf{I})$ .

Figure 10.

Strain distribution $E_{\mathrm{zz}}$ at maximum penetration state for $T=10°\mathrm{C}$ , $h_0=0.45\mathrm{m},r_B=15\mathrm{mm}$ .

Damping properties of the material are often characterized through the rebound resilience or loss factor in the ball drop test, see Emminger et al. (2021). The loss factor is given by $\tan \delta =h_{r,1}/h_0$ , often also referred to as $R$ . Both in the experiment and in the simulation, a number of parameters can be changed in order to better study the behavior. For example, the variation of the original height $h_0$ , the ball’s radius $r_B$ , the actual temperature as well as the order of the material model can provide suitable information.

Figure 11 shows the variation of the rebound resilience for different temperatures. From these results, we see that the third and fourth order material model lead to similar results; thus, a restriction to third order is viable to save computational effort. On the other hand, we observe that by decreasing the original drop height to 0.25 m, an increase in the rebound resilience is achieved. For increasing temperatures, the resilience increases independently of the drop height. At higher temperature, the material acts less viscous.

Figure 11.

Rebound resilience comparison.

The computational results are compared to measurements taken in an according setup. For this purpose, an experimental setup has been designed, where the ball is released from an electromagnetic mount. A high-speed optical camera was used to track the trajectory of the ball, as shown in Figure 12. In this figure the six frames show the trajectory of the ball drop experiment. Video frames have been selected in such a manner that impact as well as maximum rebound height are visible. Indicated times do not correspond to simulation times but are computed according to the camera frame rate, being 1600 fps. More information on the setup and procedure of these measurements can be found in Preuer (2018). There, further samples were measured and compared using the same method.

Figure 12.

Video frames of ball trajectory and impact.

Care was taken to ensure that the surface of the sample was neatly powdered before each measurement in order to exclude possible friction and adhesion effects. Several of these measurements were carried out and are shown as dotted lines in Figure 13. The trajectories obtained from simulations, as shown in Figure 9, are compared to these measurement data curves around the time of impact. Again, third and fourth order show no difference for the height course and also the velocity.

Figure 13.

Ball-drop measurements for $h_0=0.45\mathrm{m}$ , $r_B=15\mathrm{mm}$ at $T=30°\mathrm{C}$ .

It shows that the depth of penetration observed in simulation and experiment fit closely. In simulations, however, a higher rebound resilience of $R=72.6\%$ is observed as compared to $R=52.9\%$ in the experiments. Consistently, the velocity of the ball after impact is higher. It is not clear where the discrepancy originates from; certainly the friction-less contact model, the assumption of perfect rotational symmetry in displacements and also the rigid bedding of the PDMS specimen add to the computed rebound height, as dissipation is inhibited.

In a second set of experiments, a ball of radius $r_B=5\mathrm{mm}$ was dropped again from the original drop height of $h_0=0.45\mathrm{m}$ , the ambient temperature is at $T=30°\mathrm{C}$ . In Figure 14, corresponding measurements and simulation results are compared. Again, correspondence is observed for the indentation depth. Concerning rebound resilience and velocity after impact, the results fit better for this smaller ball. The rebound resilience in simulation is found at $R=50.38\%$ , whereas the measured rebound resilience was determined as $R=42.8\%$

Figure 14.

Ball-drop measurements for $h_0=0.45\mathrm{m}$ , $r_B=5\mathrm{mm}$ at $T=30°\mathrm{C}$ .

7. Conclusion

In the present article, we have introduced a thermodynamically consistent characterization for viscoelastic materials for modeling the behavior of the PDMS Sylgard 184. To determine the model’s parameters, a dynamical temperature mechanical analysis was conducted in a limited frequency range. The time-temperature superposition principle was used successfully to generate an extended frequency response. For this data, Prony-parameters based on the linear Wiechert model were identified using a non-negative least squares fitting algorithm. To verify this constitutive model, a mixed finite element model for rotationally symmetric problems was created. A ball-drop test setup was designed in order to conduct measurements on the rebound resilience for different ball sizes and drop heights. Data generated in the finite element simulation were then compared to optical tracking data from the experiment. Our numerical results match the measured data almost perfectly as far as the indentation depth is concerned. This also confirms the accuracy of linear parameters in a non-linear framework Rebound resiliences differ, however, as we observe less dissipation in computations than in experiments. To analyse these deviations in detail, the effects of friction-less contact and the imposed radial symmetry in computations need to be quantified in the future. Additionally, DTMA measurements in an even larger range of temperature may lead to a better characterization of the viscous behavior at high frequencies, which are certainly excited in a ball-drop test.