# Thermodynamically consistent modelling of viscoelastic solids under finite
strain

Mario Kunzemann, Leonhard K. Doppelbauer, Rene Preuer and Astrid Pechstein

###### Abstract

The present article is concerned with modelling 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 PDMS 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 [1] and [30], 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. [3]),
the insulation of electronic components, or coating of sensitive materials to
protect them from external effects (cf. [10]). 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 [31].
Treatment of surface layers was used to generate custom wrinkling patterns in
specimen in [15].

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, [26] 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, [5] fitted different hyperelastic material laws,
including Mooney-Rivlin, Ogden and Yeoh models. The influence of surface
patterns on the elastic response was analyzed by [18]. Characterizations of
the damping properties are scarce, [20] 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. [7] 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, [22] treat the applicability of time-
temperature-superposition to Sylgard 184.
We propose an approach that ensures thermodynamic consistency when dealing
with large deformations in viscoelastic solids. The ensuing 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 [17], 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 [16]. 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 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 modelling

###  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
$\boldsymbol{X}$ in reference configuration. The motion of the body is
described by a time-dependent mapping ${\boldsymbol{\chi}}$, which maps each
material point to its current spatial position at time $t$,
${\boldsymbol{x}}={\boldsymbol{\chi}}({\boldsymbol{X}},t)$. We use standard
notations for all relevant kinematic quantities, starting with the
displacement field ${\boldsymbol{u}}={\boldsymbol{x}}-{\boldsymbol{X}}$.
Throughout the following, we identify derivatives of fields with respect to
reference coordinates ${\boldsymbol{X}}$ by using capital operators.
Derivatives with respect to spatial coordinates ${\boldsymbol{x}}$ are
identified using lower case operators.

The deformation gradient is then defined as
$\boldsymbol{\mathrm{F}}=\operatorname{Grad}{\boldsymbol{\chi}}=\boldsymbol{\mathrm{I}}+\operatorname{Grad}{\boldsymbol{u}}$.
For modeling large deformation viscoelasticity, an extension of the one-
dimensional rheological standard solid model [13] 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.

$\sigma$$\sigma$${E_{\infty}}$${E_{N}}$${\eta_{N}}$${E_{1}}$${\eta_{1}}$${\boldsymbol{\cdot}}$${\boldsymbol{\cdot}}$${\boldsymbol{\cdot}}$

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
$\boldsymbol{\mathrm{F}}=J^{1/3}\hat{\boldsymbol{\mathrm{F}}}$ with
$J=\det\boldsymbol{\mathrm{F}}$. Based on the above definitions, the right
Cauchy-Green tensor $\boldsymbol{\mathrm{C}}$ and isochoric component of the
right Cauchy-Green tensor $\hat{\boldsymbol{\mathrm{C}}}$ are given by,

$\displaystyle\boldsymbol{\mathrm{C}}$
$\displaystyle=\boldsymbol{\mathrm{F}}^{\top}\cdot\boldsymbol{\mathrm{F}},$
(1) $\displaystyle\hat{\boldsymbol{\mathrm{C}}}$
$\displaystyle=\hat{\boldsymbol{\mathrm{F}}}^{\top}\cdot\hat{\boldsymbol{\mathrm{F}}}=J^{-2/3}\boldsymbol{\mathrm{C}}.$
(2)

For incompressible materials, its valid to say that we have $J\equiv 1$ which
leads to $\hat{\boldsymbol{\mathrm{F}}}\equiv\boldsymbol{\mathrm{F}}$ and
therefore leads to
$\hat{\boldsymbol{\mathrm{C}}}\equiv\boldsymbol{\mathrm{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{\boldsymbol{\mathrm{F}}}_{e,i}$ and a viscoelastic part
$\boldsymbol{\mathrm{F}}_{v,i}$, see Figure 2, or see e.g., [2],

$\hat{\boldsymbol{\mathrm{F}}}=\hat{\boldsymbol{\mathrm{F}}}_{e,i}\cdot\boldsymbol{\mathrm{F}}_{v,i},\text{
for }i\in\\{1,\dots N\\}.$ (3)

Above, the viscoelastic parts $\boldsymbol{\mathrm{F}}_{v,i}$ of the
deformation resemble the transformation to the stress-free $i^{th}$
intermediate configuration, where $\boldsymbol{\mathrm{F}}_{v,i}$ is not
necessarily compatible. However, it is assumed that the viscoelastic
contributions are modeled in an isochoric way, such that
$J_{v,i}=\det\boldsymbol{\mathrm{F}}_{v,i}\equiv 1$ for all $i=1\dots N$.
Based on the $i^{th}$ elastic isochoric part of the deformation gradient
$\hat{\boldsymbol{\mathrm{F}}}_{e,i}=\hat{\boldsymbol{\mathrm{F}}}\cdot\boldsymbol{\mathrm{F}}_{v,i}^{-1}$,
the corresponding elastic and viscoelastic Cauchy-Green tensors can be defined
via,

$\displaystyle\hat{\boldsymbol{\mathrm{C}}}_{e,i}$
$\displaystyle=\hat{\boldsymbol{\mathrm{F}}}_{e,i}^{\top}\cdot\hat{\boldsymbol{\mathrm{F}}}_{e,i}={\boldsymbol{\mathrm{F}}}^{-\top}_{v,i}\cdot\hat{\boldsymbol{\mathrm{C}}}\cdot{\boldsymbol{\mathrm{F}}}_{v,i}^{-1},$
(4) $\displaystyle\boldsymbol{\mathrm{C}}_{v,i}$
$\displaystyle=\boldsymbol{\mathrm{F}}_{v,i}^{\top}\cdot\boldsymbol{\mathrm{F}}_{v,i}.$
(5)

…$\boldsymbol{\mathrm{F}}_{v,N}$$\boldsymbol{\mathrm{F}}_{v,1}$$\hat{\boldsymbol{\mathrm{F}}}$$\boldsymbol{\mathrm{F}}_{e,1}$$J^{1/3}\boldsymbol{\mathrm{I}}$$\boldsymbol{\mathrm{F}}$$\boldsymbol{\mathrm{F}}_{e,N}$

Figure 2: Representation of the multiplicative decomposition.

###  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 $\Psi$ shall
be defined for independent Cauchy-Green tensor $\boldsymbol{\mathrm{C}}$ or
deformation gradient $\boldsymbol{\mathrm{F}}$. The viscoelastic Cauchy-Green
tensors $\boldsymbol{\mathrm{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 [2]. We consider a decomposition into a hyperelastic
potential $\Psi_{\infty}$ and a viscoelastic energy densities $\Psi_{v,i}$,
which shall be defined in the following,

$\Psi(\boldsymbol{\mathrm{C}},\boldsymbol{\mathrm{C}}_{v,i},p)=\Psi_{\infty}({\boldsymbol{\mathrm{C}}},p)+\sum_{i=1}^{N}\Psi_{v,i}(\hat{\boldsymbol{\mathrm{C}}},\boldsymbol{\mathrm{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. [9]), which can be written as,

$\displaystyle\mathcal{D}$
$\displaystyle=\frac{1}{2}\boldsymbol{\mathrm{T}}:\dot{\boldsymbol{\mathrm{C}}}-\dot{\Psi}\geq
0.$ (7)

Above, $\boldsymbol{\mathrm{T}}$ is the second Piola-Kirchhoff stress.
Computing the rate $\dot{\Psi}$ explicitly, equation (7) becomes,

$\mathcal{D}=-\left(\frac{\partial\Psi}{\partial\boldsymbol{\mathrm{C}}}-\frac{1}{2}\boldsymbol{\mathrm{T}}\right):\dot{\boldsymbol{\mathrm{C}}}-\sum_{i=1}^{N}\frac{\partial\Psi}{\partial\boldsymbol{\mathrm{C}}_{v,i}}:\dot{\boldsymbol{\mathrm{C}}}_{v,i}-\frac{\partial\Psi}{\partial
p}\dot{p}\geq 0.$ (8)

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

$\displaystyle\boldsymbol{\mathrm{T}}=2\frac{\partial\Psi}{\partial\boldsymbol{\mathrm{C}}}.$
(9)

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

$\displaystyle\mathcal{D}$
$\displaystyle=-\sum_{i=1}^{N}\frac{\partial\Psi}{\partial\boldsymbol{\mathrm{C}}_{v,i}}:\dot{\boldsymbol{\mathrm{C}}}_{v,i}\geq
0.$ (10)

#### Free energy

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

${\Psi}_{\infty}({\boldsymbol{\mathrm{C}}},p)=\frac{\mu_{\infty}}{2}\left(\operatorname{tr}{\mathbf{C}}-3\right)-\mu_{\infty}\log
J+p\log J.$ (11)

Pressure $p$ serves as a Lagrangian multiplier to ensure the incompressibility
constraint,

$\frac{\partial\Psi_{\infty}}{\partial p}=\log J=0,$

such that $J\equiv 1$.

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

$\Psi_{v,i}({\hat{\boldsymbol{\mathrm{C}}},\mathbf{C}_{{v},i}})=\frac{\mu_{v,i}}{2}\left(\hat{\boldsymbol{\mathrm{C}}}:\boldsymbol{\mathrm{C}}_{v,i}^{-1}-3\right).$
(12)

#### Dissipation

Viscoelastic evolution is defined implicitly through a _dissipation function_
$\Phi$ such that,

$\mathcal{D}=\sum_{i=1}^{N}\frac{\partial\Phi}{\partial\dot{\boldsymbol{\mathrm{C}}}_{v,i}}:\dot{\boldsymbol{\mathrm{C}}}_{v,i}.$
(13)

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

$\displaystyle\Phi=\sum_{i=1}^{N}\Phi_{i}(\boldsymbol{\mathrm{C}}_{v,i},\dot{\boldsymbol{\mathrm{C}}}_{v,i}),$
(14)

with,

$\displaystyle\Phi_{i}(\boldsymbol{\mathrm{C}}_{v,i},\dot{\boldsymbol{\mathrm{C}}}_{v,i})=\frac{1}{12}\sum_{i=1}^{N}\eta_{i}({\boldsymbol{\mathrm{C}}}_{v,i}^{-1}\cdot\dot{\boldsymbol{\mathrm{C}}}_{v,i}):({\boldsymbol{\mathrm{C}}}_{v,i}^{-1}\cdot\dot{\boldsymbol{\mathrm{C}}}_{v,i}).$
(15)

To analyze the implications of free energy $\Psi$ and dissipation function
$\Phi$ on the evolution of $\boldsymbol{\mathrm{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_{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}$, i.e. the modulus of the single
spring. Its instantaneous modulus $E_{0}$ can be found as,

$\displaystyle E_{0}=\sum_{i=1}^{N}E_{i}+E_{\infty},$ (16)

for a thorough discussion we refer to [13]. For incompressible materials with
a Poisson ratio of $\nu=0.5$, shear moduli $\mu_{i}$ and Young’s moduli
$E_{i}$ are connected via,

$\displaystyle\mu_{i}$ $\displaystyle=\frac{E_{i}}{1+2\nu}=\frac{E_{i}}{3},$
$\displaystyle i\in\\{1,\dots N,\infty\\}.$ (17)

Additionally, $\tau_{i}$ defines the relaxation time of the $i^{\text{th}}$
Maxwell element by,

$\displaystyle\eta_{i}=\tau_{i}E_{i}.$ (18)

At time $t>0$, the relaxation modulus is found as, see [4],

$\displaystyle E(t)=\sum_{i=1}^{N}E_{i}e^{-t/\tau_{i}}+E_{\infty}.$ (19)

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

$\displaystyle
E^{*}(\omega)=\underbrace{E_{\infty}+\sum_{i=1}^{N}\frac{E_{i}\tau_{i}^{2}\omega^{2}}{1+\tau_{i}^{2}\omega^{2}}}_{\begin{subarray}{c}E^{\prime}(\omega)\end{subarray}}+j\underbrace{\sum_{i=1}^{N}\frac{E_{i}\tau_{i}\omega}{1+\tau_{i}^{2}\omega^{2}}}_{\begin{subarray}{c}E^{\prime\prime}(\omega)\end{subarray}},$
(20)

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^{\prime\prime}(\omega)$, see [28].

Furthermore we define $e_{i}=E_{i}/E_{0}$ with the instantaneous modulus
$E_{0}$ as defined in (16). 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 dimensions of the Sylgard 184 sample were
$$10\text{\,}\mathrm{mm}$\times$40\text{\,}\mathrm{mm}$\times$3.2\text{\,}\mathrm{mm}$$
($W\times H\times D$) and a clamping length of $20\text{\,}\mathrm{mm}$ was
chosen. The static strain was $3\text{\,}\mathrm{\char 37\relax}$, the dynamic
strain $2\text{\,}\mathrm{\char 37\relax}$ with a control range of
$0.0045\text{\,}\mathrm{\char 37\relax}$ and the pre-stress was
$0.05\text{\,}\mathrm{MPa}$. The pre-cooling time in the temperature chamber
was at least $2$ minutes. Measurements were taken in the range $0$ to
$40\text{\,}\mathrm{\SIUnitSymbolCelsius}$ in
$2.5\text{\,}\mathrm{\SIUnitSymbolCelsius}$ steps between $0$ to
$100\text{\,}\mathrm{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\text{\,}\mathrm{Hz}$, we did not use the entire measurement curve,
but only from $0\text{\,}\mathrm{Hz}$ to $60\text{\,}\mathrm{Hz}$, as shown in
Figure 3(a) and Figure 3(b).

(a) Measured storage modulus $E^{\prime}$.

(b) Measured loss modulus $E^{\prime\prime}$.

Figure 3: DTMA measured data

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 behaviour, see also [4], [12].

###  3.3 Time-temperature superposition

Concerning time-temperature superposition, Sylgard 184 is treated as a
thermorheologically complex material, cf. [14]. In the modeled area, the
cross-linked material is well above the glass transition temperature $T_{g}$
and the melting temperature $T_{m}$, cf. [8]. 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 [29], the following general form of time-temperature shifts
was identified as,

$\displaystyle E_{T}^{*}(\omega)=b_{T}E_{T_{0}}^{*}(a_{T}\omega).$ (21)

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,

$\displaystyle b_{T}=\frac{T_{0}\rho_{0}}{T\rho},$ (22)

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\text{\,}\mathrm{\SIUnitSymbolCelsius}$ to
$40\text{\,}\mathrm{\SIUnitSymbolCelsius}$, 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 [21], the WLF model
can be used very well for a temperature range of
$T_{g}<T<T_{g}+$100\text{\,}\mathrm{\SIUnitSymbolCelsius}$$.

[6] and also [19] found that the glass transition temperature of Sylgard 184
is in the range of $-150\text{\,}\mathrm{\SIUnitSymbolCelsius}$
$<T_{g}<$-120\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ and the melting
temperature $T_{m}$ lies between $-50\text{\,}\mathrm{\SIUnitSymbolCelsius}$
and $-70\text{\,}\mathrm{\SIUnitSymbolCelsius}$. This implies that the WLF
model is not applicable to Sylgard 184 at ambient temperatures. [21] states
that in this relatively high temperature range of
$T>T_{g}+$100\text{\,}\mathrm{\SIUnitSymbolCelsius}$$, the so-called Arrhenius
equation is applicable,

$\displaystyle a_{T}=e^{\alpha\left(\frac{1}{T}-\frac{1}{T_{0}}\right)}.$ (23)

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\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ 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\text{\,}\mathrm{\SIUnitSymbolCelsius}$ to
$40\text{\,}\mathrm{\SIUnitSymbolCelsius}$ 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 minimise the squared vertical distance of the respective
data curves when shifted by $b_{T}$ and $a_{T}$ as in (22) and (23),
respectively. As the shifted frequencies for the data points do no 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,

$\displaystyle\sum_{T>T_{0}}\sum_{i=1}^{n_{\omega}}\left|b_{T}E^{*}_{T_{0}}(a_{T}\omega_{i})-E^{*}_{T}(\omega_{i})\right|^{2}\to\min_{\alpha>0}.$
(24)

For the available DTMA data, this procedure results in
$\alpha=$7389.124\text{\,}\mathrm{K}$$. Then, both the storage and loss moduli
from all DTMA data temperatures can be shifted according to (21) to get a
master curve at temperature $T_{0}=$20\text{\,}\mathrm{\SIUnitSymbolCelsius}$$
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.

(a) Storage modulus $E_{T_{0}}^{\prime}$ at
$T_{0}=$20\text{\,}\mathrm{\SIUnitSymbolCelsius}$$.

(b) Loss modulus $E_{T_{0}}^{\prime\prime}$ at
$T_{0}=$20\text{\,}\mathrm{\SIUnitSymbolCelsius}$$.

Figure 4: Master curves for storage and loss modulus at
$T_{0}=$20\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ gained through time-
temperature superposition.

###  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}\leq 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\text{\,}\mathrm{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\text{\,}\mathrm{Hz}$ up to almost $1000\text{\,}\mathrm{Hz}$ (see Figure
4), fitting Prony-parameters for models of third and fourth order seems
appropriate. The respective results are depicted in Figure 5. Table 1 and
Table 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 (17) and (18).

Figure 5: Fitted Prony-series, comparison $N=3$ and $N=4$. Table 1: Prony-parameters for $N=3$ and $T_{0}=$20\text{\,}\mathrm{\SIUnitSymbolCelsius}$$.  $i$ |  $E_{i}$ in $\mathrm{MPa}$ |  $\tau_{i}$ in $\mathrm{s}$ |  $\mu_{i}$ in $\mathrm{MPa}$ |  $\eta_{i}$ in $\mathrm{MPa}\text{\,}\mathrm{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\text{\,}\mathrm{\SIUnitSymbolCelsius}$$. $i$ |  $E_{i}$ in $\mathrm{MPa}$ |  $\tau_{i}$ in $\mathrm{s}$ |  $\mu_{i}$ in $\mathrm{MPa}$ |  $\eta_{i}$ in $\mathrm{MPa}\text{\,}\mathrm{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 Table 1 and Table 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 [28].

##  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 [24, 25, 23]. For reasons of simplicity of presentation, we
restrict ourselves to the case of a conservative volume force density
${\boldsymbol{f}}$. With $\rho_{0}$ the mass density, the kinetic energy
density $T$ and local power of external forces $\mathcal{P}_{ext}$ read,

$\displaystyle T$
$\displaystyle=\frac{1}{2}\rho_{0}\dot{{\boldsymbol{u}}}\cdot\dot{{\boldsymbol{u}}},$
$\displaystyle\mathcal{P}_{ext}$
$\displaystyle={\boldsymbol{f}}\cdot\dot{{\boldsymbol{u}}}.$ (25)

Consider a finite time step $[t_{0},t_{0}+\Delta t]$ with known initial
conditions at time $t_{0}$. We introduce the _incremental potential_ for this
time step as,

$\displaystyle\begin{split}\Pi_{t}^{t+\Delta t}&=\int_{t_{0}}^{t_{0}+\Delta
t}(\dot{T}+\dot{\Psi}+\Phi-\mathcal{P}_{ext})\,dt\\\
&=\underbrace{T_{t_{0}+\Delta
t}-T_{t_{0}}}_{\text{kinetic}}+\underbrace{\Psi_{t_{0}+\Delta
t}-\Psi_{t_{0}}}_{\text{stored}}+\underbrace{\int_{t_{0}}^{t_{0}+\Delta t}\Phi
dt}_{\text{dissipated}}-\underbrace{\int_{t_{0}}^{t_{0}+\Delta
t}\mathcal{P}_{ext}\,dt}_{\text{external}}.\end{split}$ (26)

The incremental potential $\Pi_{t_{0}}^{t_{0}+\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 $\Delta t$ approaching zero, we arrive at,

$\displaystyle\dot{T}+\dot{\Psi}+\Phi-\mathcal{P}_{ext}\to\min_{\dot{{\boldsymbol{u}}}}\min_{\dot{\boldsymbol{\mathrm{C}}}_{v,i}}\max_{\dot{p}}.$
(27)

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

$\begin{split}\rho\ddot{\boldsymbol{u}}\cdot\dot{\boldsymbol{u}}+\frac{\partial\Psi}{\partial\boldsymbol{\mathrm{C}}}\colon\dot{\boldsymbol{\mathrm{C}}}+\sum_{i=1}^{N}\frac{\partial\Psi}{\partial\boldsymbol{\mathrm{C}}_{v,i}}\colon\dot{\boldsymbol{\mathrm{C}}}_{v,i}+\frac{\partial\Psi}{\partial{p}}\dot{p}+\Phi-{\boldsymbol{f}}\cdot{\boldsymbol{\dot{u}}}\to\min_{\dot{{\boldsymbol{u}}}}\min_{\dot{\boldsymbol{\mathrm{C}}}_{v,i}}\max_{\dot{p}}.\end{split}$
(28)

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{split}&\int_{\mathcal{B}}\rho\ddot{{\boldsymbol{u}}}\cdot\delta{\boldsymbol{u}}\,dV+\int_{\mathcal{B}}\left(\frac{\partial\Psi}{\partial\boldsymbol{\mathrm{C}}}\colon\delta\boldsymbol{\mathrm{C}}+\sum_{i=1}^{N}\frac{\partial\Psi}{\partial\boldsymbol{\mathrm{C}}_{v,i}}\colon\delta{\boldsymbol{\mathrm{C}}}_{v,i}+\frac{\partial\Psi}{\partial{p}}\delta{{p}}\right)\,dV\\\
&+\int_{\mathcal{B}}\sum_{i=1}^{N}\frac{\partial\Phi}{\partial\dot{\boldsymbol{\mathrm{C}}}_{v,i}}\colon\delta\boldsymbol{\mathrm{C}}_{v,i}\,dV-\int_{\mathcal{B}}{\boldsymbol{f}}\cdot\delta{\boldsymbol{u}}\,dV=0.\end{split}$
(29)

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

$\displaystyle\int_{\mathcal{B}}\left(\rho\ddot{{\boldsymbol{u}}}\cdot\delta{\boldsymbol{u}}+\delta\Psi+\delta\Phi-{\boldsymbol{f}}\cdot\delta{\boldsymbol{u}}\right)\,dV$
$\displaystyle=0,$ (30)

with,

$\displaystyle\delta\Psi:=\frac{\partial\Psi}{\partial\boldsymbol{\mathrm{C}}}\colon\delta\boldsymbol{\mathrm{C}}+\sum_{i=1}^{N}\frac{\partial\Psi}{\partial\boldsymbol{\mathrm{C}}_{v,i}}\colon\delta{\boldsymbol{\mathrm{C}}}_{v,i}+\frac{\partial\Psi}{\partial{p}}\delta{{p}},$
(32)

and

$\displaystyle\delta\Phi:=\sum_{i=1}^{N}\frac{\partial\Phi}{\partial\dot{\boldsymbol{\mathrm{C}}}_{v,i}}\colon\delta\boldsymbol{\mathrm{C}}_{v,i}.$
(33)

The variational equation (29) implicitly defines the evolution law for
$\dot{\boldsymbol{\mathrm{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 (15), the standard model as e.g. proposed by [13] is
re-established. We collect all terms containing
$\delta\boldsymbol{\mathrm{C}}_{v,i}$ from equation (29), and obtain,

$\int_{\mathcal{B}}\left(\frac{\partial\Psi_{v,i}}{\partial\boldsymbol{\mathrm{C}}_{v,i}}+\frac{1}{6}\eta_{i}\left(\boldsymbol{\mathrm{C}}_{v,i}^{-1}\cdot\dot{\boldsymbol{\mathrm{C}}}_{v,i}\cdot\boldsymbol{\mathrm{C}}_{v,i}^{-1}\right)\right)\colon\delta\boldsymbol{\mathrm{C}}_{v,i}\,dV=0.$
(34)

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

$\displaystyle\frac{d}{dt}(\det\boldsymbol{\mathrm{C}}_{v,i})=\boldsymbol{\mathrm{C}}_{v,i}^{-1}\colon\dot{\boldsymbol{\mathrm{C}}}_{v,i}=0,$
(35)
$\displaystyle\delta(\det\boldsymbol{\mathrm{C}}_{v,i})=\boldsymbol{\mathrm{C}}_{v,i}^{-1}\colon\delta{\boldsymbol{\mathrm{C}}}_{v,i}=0.$
(36)

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

$\displaystyle\int_{\mathcal{B}}{\left(\frac{\partial\Psi_{v,i}}{\partial\boldsymbol{\mathrm{C}}_{v,i}}+\gamma\boldsymbol{\mathrm{C}}_{v,i}^{-1}+\frac{1}{6}\eta_{i}\boldsymbol{\mathrm{C}}_{v,i}^{-1}\cdot\dot{\boldsymbol{\mathrm{C}}}_{v,i}\cdot\boldsymbol{\mathrm{C}}_{v,i}^{-1}\right)}\colon\delta\boldsymbol{\mathrm{C}}_{v,i}\,dV={0}.$

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

$\frac{\partial\Psi_{v,i}}{\partial\boldsymbol{\mathrm{C}}_{v,i}}+\gamma\boldsymbol{\mathrm{C}}_{v,i}^{-1}+\frac{1}{6}\eta_{i}\boldsymbol{\mathrm{C}}_{v,i}^{-1}\cdot\dot{\boldsymbol{\mathrm{C}}}_{v,i}\cdot\boldsymbol{\mathrm{C}}_{v,i}^{-1}=\boldsymbol{\mathrm{0}},$
(37)

which directly implies,

$\dot{\boldsymbol{\mathrm{C}}}_{v,i}=-\frac{6}{\eta_{i}}\boldsymbol{\mathrm{C}}_{v,i}\cdot\left(\frac{\partial\Psi_{v,i}}{\partial\boldsymbol{\mathrm{C}}_{v,i}}+\gamma\boldsymbol{\mathrm{C}}_{v,i}^{-1}\right)\cdot\boldsymbol{\mathrm{C}}_{v,i}.$
(38)

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

$\displaystyle
0=\boldsymbol{\mathrm{C}}_{v,i}^{-1}:\left[\boldsymbol{\mathrm{C}}_{v,i}\cdot\left(\frac{\partial\Psi_{v,i}}{\partial\boldsymbol{\mathrm{C}}_{v,i}}+\gamma\boldsymbol{\mathrm{C}}_{v,i}^{-1}\right)\cdot\boldsymbol{\mathrm{C}}_{v,i}\right]=\boldsymbol{\mathrm{I}}:\left[\left(\boldsymbol{\mathrm{C}}_{v,i}\cdot\frac{\partial\Psi_{v,i}}{\partial\boldsymbol{\mathrm{C}}_{v,i}}+\gamma\boldsymbol{\mathrm{I}}\right)\cdot\boldsymbol{\mathrm{C}}_{v,i}\cdot\boldsymbol{\mathrm{C}}_{v,i}^{-1}\right].$

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

$\displaystyle
0=\boldsymbol{\mathrm{C}}_{v,i}\cdot\frac{\partial\Psi_{v,i}}{\partial\boldsymbol{\mathrm{C}}_{v,i}}:\boldsymbol{\mathrm{I}}+\gamma\boldsymbol{\mathrm{I}}:\boldsymbol{\mathrm{I}}=-\boldsymbol{\mathrm{M}}_{v,i}:\boldsymbol{\mathrm{I}}+3\gamma.$

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

$\displaystyle\boldsymbol{\mathrm{M}}_{v,i}=-\boldsymbol{\mathrm{C}}_{v,i}\cdot\frac{\partial\Psi_{v,i}}{\partial\boldsymbol{\mathrm{C}}_{v,i}}.$
(39)

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

$\displaystyle\gamma=\frac{1}{3}\operatorname{tr}\boldsymbol{\mathrm{M}}_{v,i}.$
(40)

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

$\displaystyle\dot{\boldsymbol{\mathrm{C}}}_{v,i}$
$\displaystyle=\frac{6}{\eta_{i}}\operatorname{dev}(\boldsymbol{\mathrm{M}}_{v,i})\cdot\boldsymbol{\mathrm{C}}_{v,i},$
(41)

which is stated in similar form by [2] for viscoelastic electrostrictive
polymers. To compare the final evolution law to standard theory, we explicitly
compute the Mandel-type stresses for the viscoelastic potential $\Psi_{v,i}$
from equation (12) and obtain,

$\displaystyle\boldsymbol{\mathrm{M}}_{v,i}$
$\displaystyle=-\frac{\mu_{v,i}}{2}\left(\hat{\boldsymbol{\mathrm{C}}}\cdot\boldsymbol{\mathrm{C}}_{v,i}^{-1}\right).$
(42)

After combining the two equations equation (41) and equation (42) we finally
obtain the evolution law for $\boldsymbol{\mathrm{C}}_{v,i}$ as stated in
[13],

$\displaystyle\frac{2}{3}\eta\dot{\boldsymbol{\mathrm{C}}}_{v,i}=2\mu_{v,i}\left(\hat{\boldsymbol{\mathrm{C}}}-\frac{1}{3}\operatorname{tr}\left(\hat{\boldsymbol{\mathrm{C}}}\cdot\boldsymbol{\mathrm{C}}_{v,i}^{-1}\right)\boldsymbol{\mathrm{C}}_{v,i}\right).$
(43)

##  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 ${\boldsymbol{e}}_{z}$ as axis of rotation.
Deformation gradient $\boldsymbol{\mathrm{F}}$ and right Cauchy-Green
$\boldsymbol{\mathrm{C}}$ are defined in the standard way, assuming that the
solution is independent of the angular coordinate $\varphi$.

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. [27]), where both ${\boldsymbol{u}}$ and $p$ are
represented by continuous basis functions ${\boldsymbol{N}}^{(u)}_{i}$ and
$N^{(p)}_{i}$ respectively, and the approximation order for ${\boldsymbol{u}}$
is one higher than for $p$. Additionally, kinematic boundary conditions
${\boldsymbol{u}}={\boldsymbol{0}}$ on the boundary part $\Gamma_{fix}$ have
to be satisfied. Summing up, finite element approximations for
${\boldsymbol{u}}$ and $p$ satisfy,

$\displaystyle{\boldsymbol{u}}$
$\displaystyle\in\\{{\boldsymbol{v}}:{\boldsymbol{v}}|_{T}\in(P^{k}(T))^{2},{\boldsymbol{v}}\text{
cont.},{\boldsymbol{v}}|_{\Gamma_{fix}}={\boldsymbol{0}}\\},$ (44)
$\displaystyle{\boldsymbol{u}}$
$\displaystyle=\sum_{i=1}^{n_{u}}q_{i}^{(u)}{\boldsymbol{N}}^{(u)}_{i},$ (45)
$\displaystyle p$ $\displaystyle\in\\{q:q|_{T}\in P^{k-1}(T),q\text{
cont.}\\},$ (46) $\displaystyle p$
$\displaystyle=\sum_{i=1}^{n_{p}}q_{i}^{(p)}N^{(p)}_{i}.$ (47)

In addition, it is necessary to provide a finite-dimensional approximation for
the viscoelastic internal variables $\boldsymbol{\mathrm{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 $\boldsymbol{\mathrm{C}}_{v,i}$ in terms of the
independent viscoelastic strains $s_{i,rr}$, $s_{i,zz}$ and $s_{i,rz}$. The
$i^{th}$ viscoelastic Cauchy-Green tensor respects the axisymmetric structure,

$\boldsymbol{\mathrm{C}}_{v,i}=\begin{bmatrix}2s_{i,rr}+1&0&2s_{i,rz}\\\
0&c_{i,\varphi\varphi}&0\\\ 2s_{i,rz}&0&2s_{i,zz}+1\end{bmatrix}.$ (48)

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

$\displaystyle
c_{i,\varphi\varphi}=\left((2s_{i,rr}+1)(1+2s_{i,zz})-4s_{i,rz}\right)^{-1}.$
(49)

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

$\displaystyle\begin{split}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^{(s)}_{i}.\end{split}$
(50)

###  5.2 Time integration

Introducing the finite element discretizations equations (46), (47) and (50)
in D’Alembert’s principle (29) leads to a system of ordinary differential
equations for the degree of freedom vector ${\boldsymbol{q}}$. As the physical
problem includes accelerations $\ddot{{\boldsymbol{u}}}$ and viscoelastic
strain rates $\dot{\boldsymbol{\mathrm{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
${\boldsymbol{v}}=\dot{{\boldsymbol{u}}}$ are approximated by an independent
linear combination of displacement basis functions
${\boldsymbol{N}}^{(u)}_{i}$, via additional degrees of freedom $q^{(v)}_{i}$.
For any time step $[t_{n},t_{n+1}=t_{n}+\Delta t]$, the quantities
${\boldsymbol{u}}_{n}$, ${\boldsymbol{v}}_{n}$, $p_{n}$ as well as
$\boldsymbol{\mathrm{C}}_{v,i,n}$, $\dot{\boldsymbol{\mathrm{C}}}_{v,i,n}$ (or
rather the viscoelastic strains $s_{i,\alpha\beta,n}$ and strain rates
$\dot{s}_{i,\alpha\beta,n}$) at time $t_{n}$ are assumed to be available from
the previous step. According updates ${\boldsymbol{u}}$, ${\boldsymbol{v}}$,
$p$, $s_{i,\alpha\beta}$, and $\dot{s}_{i,\alpha\beta}$ at time
$t_{n+1}=t_{n}+\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{\boldsymbol{v}}$ as the (independent) variation
of the velocity ${\boldsymbol{v}}$,

$\displaystyle\int_{\mathcal{B}}$
$\displaystyle\rho\left({\boldsymbol{u}}-{\boldsymbol{u}}_{n}-\Delta
t\,{\boldsymbol{v}}_{n}\right)\cdot\delta{\boldsymbol{u}}\,dV+\left(\Delta
t\right)^{2}\int_{\mathcal{B}}\beta(\delta\Psi+\delta\Phi-{\boldsymbol{f}}\cdot\delta{\boldsymbol{u}})dV$
$\displaystyle+\left(\Delta
t\right)^{2}\int_{\mathcal{B}}\left(\tfrac{1}{2}-\beta\right)(\delta\Psi_{n}+\delta\Phi_{n}-{\boldsymbol{f}}_{n}\cdot\delta{\boldsymbol{u}})\,dV=0,$
(51) $\displaystyle\int_{\mathcal{B}}$
$\displaystyle\rho\left({\boldsymbol{v}}-{\boldsymbol{v}}_{n}\right)\cdot\delta{\boldsymbol{v}}\,dV+\Delta
t\int_{\mathcal{B}}\gamma\left(\left.\frac{\partial\Psi}{\partial{\boldsymbol{u}}}\right|_{t_{n+1}}\cdot\delta{\boldsymbol{v}}-{\boldsymbol{f}}\cdot\delta{\boldsymbol{v}}\right)dV$
$\displaystyle+\Delta
t\int_{\mathcal{B}}\left(1-\gamma\right)\left(\left.\frac{\partial\Psi}{\partial{\boldsymbol{u}}}\right|_{t_{n}}\cdot\delta{\boldsymbol{v}}-{\boldsymbol{f}}_{n}\cdot\delta{\boldsymbol{v}}\right)\,dV=0.$
(52)

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

##  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/NGSolve111https://ngsolve.org. Via a python interface, the
developed variational time stepping scheme (51)-(52) 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\text{\,}\mathrm{mm}$$ and a height of
$h_{S}=$30\text{\,}\mathrm{mm}$$. A steel ball of radius
$r_{B}=$15\text{\,}\mathrm{mm}$$ and mass $m_{B}=$109\text{\,}\mathrm{g}$$ is
dropped from an initial height of $h_{0}=$450\text{\,}\mathrm{mm}$$ above the
cylinder’s surface. All parameters are once again summarized in Table 3.

Table 3: Parameters of the ball-drop test specimen radius | $r_{S}$ | $30\text{\,}\mathrm{mm}$
---|---|---
specimen height | $h_{S}$ | $30\text{\,}\mathrm{mm}$
specimen density | $\rho_{Syl}$ | $965\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$
drop height | $h_{0}$ | $450\text{\,}\mathrm{mm}$
ball radius | $r_{B}$ | $15\text{\,}\mathrm{mm}$
ball mass | $m_{B}$ | $109\text{\,}\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 modelled as linear
elastic body with Young’s modulus $E=$210\text{\,}\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.

$r_{B}$$r_{S}$$h_{S}$$h_{0}$$m_{B}$$\rho_{Syl}$$z$$x$$h_{1}$

Figure 6: Computational setup of the ball-drop test.

Assume $t_{0}=$0\text{\,}\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\text{\,}\mathrm{m}$$, which is close to the specimen’s
surface. This results in an initial time $t_{1}$ an initial velocity
${\boldsymbol{v}}_{1}=v_{1}{\boldsymbol{e}}_{z}$ set to,

$\displaystyle t_{1}$ $\displaystyle=\sqrt{\frac{2}{g}(h_{0}-h_{1})},$
$\displaystyle v_{1}$ $\displaystyle=-gt_{1}.$ (53)

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

$t$$h$$t_{0}$$t_{1}$$t_{2}$$h_{0}$$h_{1},v_{1}$$h_{2},v_{2}$$h_{r,1}$

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
$\Delta t=$1\text{\times}{10}^{-4}\text{\,}\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}.$ (54)

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 temperatures were chosen such that they cover the
range of DTMA measurements, we use $10,20,30,40\text{ and
}$50\text{\,}\mathrm{\SIUnitSymbolCelsius}$$ respectively. 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 8: Ball-drop trajectory,
$h_{0}=$0.45\text{\,}\mathrm{m}$,r_{B}=$15\text{\,}\mathrm{mm}$$, using
material data for $N=3$ from Table 1. Figure 9: Impact trajectory (cmp. Figure
8), solid lines indicate the ball’s trajectory, dots mark its velocity.

Damping properties of the material are often characterized through the rebound
resilience or loss factor in the ball drop test, see [11]. 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 behaviour. 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 10 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\text{\,}\mathrm{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 10: 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. 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 11. 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 11: Ball-drop measurements for
$h_{0}=$0.45\text{\,}\mathrm{m}$,r_{B}=$15\text{\,}\mathrm{mm}$$ at
$T=$30\text{\,}\mathrm{\SIUnitSymbolCelsius}$$.

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\text{\,}\mathrm{\char 37\relax}$$ is observed as compared to
$R=$52.9\text{\,}\mathrm{\char 37\relax}$$ 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.

Figure 12: Ball-drop measurements for
$h_{0}=$0.45\text{\,}\mathrm{m}$,r_{B}=$5\text{\,}\mathrm{mm}$$ at
$T=$30\text{\,}\mathrm{\SIUnitSymbolCelsius}$$.

In a second set of experiments, a ball of radius
$r_{B}=$5\text{\,}\mathrm{mm}$$ was dropped again from the original drop
height of $h_{0}=$0.45\text{\,}\mathrm{m}$$, the ambient temperature is at
$T=$30\text{\,}\mathrm{\SIUnitSymbolCelsius}$$. In Figure 12, 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\text{\,}\mathrm{\char
37\relax}$$, whereas the measured rebound resilience was determined as
$R=$42.8\text{\,}\mathrm{\char 37\relax}$$

##  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
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. 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 behaviour at high
frequencies, which are certainly excited in a ball-drop test.

## References

  * [1] Ronaldo Ariati, Flaminio Sales, Andrews Souza, Rui A Lima, and João Ribeiro.  Polydimethylsiloxane composites characterization and its applications: a review.  Polymers, 13(23):4258, 2021.
  * [2] Anna Ask, Andreas Menzel, and Matti Ristinmaa.  Electrostriction in electro-viscoelastic polymers.  Mechanics of Materials, 50:9–21, 2012.
  * [3] Alexandra Borók, Kristóf Laboda, and Attila Bonyár.  Pdms bonding technologies for microfluidic applications: A review.  Biosensors, 11(8):292, 2021.
  * [4] Hal F. Brinson and L. Catherine Brinson.  Polymer Engineering Science and Viscoelasticity - An Introduction.  Springer, Berlin, Heidelberg, 2015.
  * [5] Cátia Cardoso, Carla S. Fernandes, Rui Lima, and João Ribeiro.  Biomechanical analysis of pdms channels using different hyperelastic numerical constitutive models.  Mechanics Research Communications, 90:26–33, 2018.
  * [6] S.J Clarson, K Dodgson, and J.A Semlyen.  Studies of cyclic and linear poly(dimethylsiloxanes): 19. glass transition temperatures and crystallization behaviour.  Polymer, 26(6):930–934, 1985.
  * [7] Shinji Deguchi, Junya Hotta, Sho Yokoyama, and Tsubasa S Matsui.  Viscoelastic and optical properties of four different pdms polymers.  Journal of Micromechanics and Microengineering, 25(9):097002, 2015\.
  * [8] Ludivine Dewimille, Bruno Bresson, and Liliane Bokobza.  Synthesis, structure and morphology of poly (dimethylsiloxane) networks filled with in situ generated silica particles.  Polymer, 46(12):4135–4143, 2005.
  * [9] Issam Doghri.  Mechanics of Deformable Solids - Linear, Nonlinear, Analytical and Computational Aspects.  Springer Science and Business Media, Berlin Heidelberg, 2013.
  * [10] Ubong Eduok, Omar Faye, and Jerzy Szpunar.  Recent developments and applications of protective silicone coatings: A review of pdms functional materials.  Progress in Organic Coatings, 111:124–163, 2017.
  * [11] Carina Emminger, Umut D Çakmak, Rene Preuer, Ingrid Graz, and Zoltán Major.  Hyperelastic material parameter determination and numerical study of TPU and PDMS dampers.  Materials, 14(24):7639, 2021.
  * [12] John D. Ferry.  Viscoelastic Properties of Polymers -.  John Wiley & Sons, New York, 1980.
  * [13] Peter Haupt.  Continuum Mechanics and Theory of Materials.  Springer Science and Business Media, Berlin Heidelberg, 2013.
  * [14] ETJ Klompen and LE Govaert.  Nonlinear viscoelastic behaviour of thermorheologically complex materials: A modelling approach.  Mechanics of Time-dependent Materials - MECH TIME-DEPEND MATER, 3, 01 1999.
  * [15] André Knapp, Lisa Julia Nebel, Mirko Nitschke, Oliver Sander, and Andreas Fery.  Controlling line defects in wrinkling: a pathway towards hierarchical wrinkling structures.  Soft Matter, 17(21):5384–5392, 2021.
  * [16] Michael A Kraus, Miriam Schuster, Johannes Kuntsche, Geralt Siebert, and Jens Schneider.  Parameter identification methods for visco-and hyperelastic material models.  Glass Structures & Engineering, 2(2):147–167, 2017.
  * [17] Mario Kunzemann, Astrid Pechstein, and Alexander Humer.  Thermodynamically consistent modelling of dielectric viscoelastic solids under finite strain.  In D.A. Saravanos, A. Benjeddou, N. Chrysochoidis, and T. Theodosiou, editors, Proceedings of the X ECCOMAS Thematic Conference on Smart Structures and Materials, 2023.
  * [18] Hoo Min Lee, Jaebum Sung, Byeongjo Ko, Heewon Lee, Sangyeun Park, Hongyun So, and Gil Ho Yoon.  Modeling and application of anisotropic hyperelasticity of pdms polymers with surface patterns obtained by additive manufacturing technology.  Journal of the Mechanical Behavior of Biomedical Materials, 118:104412, 2021.
  * [19] Wanwan Li, Feng Liu, Liuhe Wei, and Tong Zhao.  Synthesis, morphology and properties of polydimethylsiloxane-modified allylated novolac/4,4-bismaleimidodiphenylmethane.  European Polymer Journal, 42(3):580–592, 2006.
  * [20] I-Kuan Lin, Kuang-Shun Ou, Yen-Ming Liao, Yan Liu, Kuo-Shen Chen, and Xin Zhang.  Viscoelastic characterization and modeling of polymer transducers for biological applications.  Journal of Microelectromechanical systems, 18(5):1087–1099, 2009\.
  * [21] Paolo Lomellini.  Williams-landel-ferry versus arrhenius behaviour: polystyrene melt viscoelasticity revised.  Polymer, 33(23):4983–4989, 1992.
  * [22] Kevin Nicholas Long and Judith Alice Brown.  A linear viscoelastic model calibration of sylgard 184.  Technical report, Sandia National Lab.(SNL-NM), Albuquerque, NM (United States), 2017.
  * [23] C Miehe, D Rosato, and B Kiefer.  Variational principles in dissipative electro-magneto-mechanics: a framework for the macro-modeling of functional materials.  International Journal for Numerical Methods in Engineering, 86(10):1225–1276, 2011.
  * [24] Christian Miehe.  Strain-driven homogenization of inelastic microstructures and composites based on an incremental variational formulation.  International Journal for numerical methods in engineering, 55(11):1285–1322, 2002.
  * [25] Christian Miehe, Nikolas Apel, and Matthias Lambrecht.  Anisotropic additive plasticity in the logarithmic strain space: modular kinematic formulation and implementation based on incremental minimization principles for standard materials.  Computer methods in applied mechanics and engineering, 191(47-48):5383–5425, 2002.
  * [26] L.C.S. Nunes.  Mechanical characterization of hyperelastic polydimethylsiloxane by simple shear test.  Materials Science and Engineering: A, 528(3):1799–1804, 2011.
  * [27] Cedric Taylor and Paul Hood.  A numerical solution of the Navier-Stokes equations using the finite element technique.  Computers & Fluids, 1(1):73–100, 1973.
  * [28] Nicholas W. Tschoegl.  The Phenomenological Theory of Linear Viscoelastic Behavior - An Introduction.  Springer Science & Business Media, Berlin Heidelberg, 2012.
  * [29] Nicholas W. Tschoegl, Wolfgang Knauss, and Igor Emri.  The effect of temperature and pressure on the mechanical properties of thermo- and/or piezorheologically simple polymeric materials in thermodynamic equilibrium – a critical review.  Mechanics of Time-Dependent Materials, 6:53–99, 2002.
  * [30] Qaiser Zaman, Khalid Mahmood Zia, Mohammad Zuber, Yahia Nasser Mabkhot, Faisal Almalki, and Taibi Ben Hadda.  A comprehensive review on synthesis, characterization, and applications of polydimethylsiloxane and copolymers.  International Journal of Plastics Technology, 23:261–282, 2019\.
  * [31] Deyong Zhu, Stephan Handschuh-Wang, and Xuechang Zhou.  Recent progress in fabrication and application of polydimethylsiloxane sponges.  Journal of Materials Chemistry A, 5(32):16467–16497, 2017.