1 Introduction

Compressible rubber components are often needed for sealings in automotive applications. The rubbery basis material gives a high flexibility and the cellular microstructure of foamed rubber provides the compressibility. During their assembly, sealing systems are pre-deformed to satisfy a permanent tightness in use, [33]. Particularly in the industrial development process, simulation tools, like the finite element method (FEM), are used to predict the geometrically complex components. For example, the closing pressure and its relaxation over time are of interest. In order to describe the mechanical properties of rubber foams accurately, suitable material models have to be used. The materials in mind are characterized by the elastic and viscoelastic properties of the rubber itself and by their cellular microscopic structure. In order to describe this interaction between the mechanics of the rubber matrix material and the geometrical porous material structure in a convenient manner, a constitutive model has to account for both.

In the simplest case, rubber-like compressible materials are modeled by the usual hyperelastic models with an additional term for the volumetric deformation, cf. [5, 16]. Then, in contrast to (incompressible) rubbers, the experimental characterization is more complex since no split in an isochoric and volumetric part can be assumed any more. For example, uniaxial or biaxial tension experiments do not longer depend on the isochoric behavior only. In practice, that means the compressibility must be determined in elaborate experiments for every value of rubber porosity. A thermomechanical characterization of foamed rubber by experimental investigations is proposed by [35]. In the complex case, the rubber’s microstructure can be measured by computed tomography and evaluated precisely in representative volume elements. This procedure, which may be formulated in the spirit of multiscale finite elements (FE\(^2\)), leads to an accurate material mapping but a huge computational effort. A compromise is a continuum-level constitutive model for foamed rubber with low to moderate porosity. Such an approach has been presented by Danielsson et al. [10]. This model is derived from a hollow sphere representative volume which is based on earlier works of Abeyaratne [18], Ball [4], Horgan [17] and many others. The homogenized model is able to predict the stress-strain behavior based on the material’s porosity and the parameter of the incompressible matrix material solely. However, it is only valid for moderate compressive deformations because the model has a singularity, when the volume change equals the volume fraction of matrix material in the considered foamed rubber [24]. Despite of these restrictions, we follow here the modeling approach of Danielsson et al. [10], because the experimental characterization requires only the mechanical data of the pore-free matrix material and the porosity of the foam.

The homogenized hollow sphere model was originally derived for Neo–Hookean matrix material [10]. A generalization of the Danielsson ansatz to higher order terms for Mooney–Rivlin series like models was introduced by Ricker et. al. [31, 32]. To avoid numerical integration problems Raghunath et al. [28] used higher order terms of the resulting derived energy density functions directly to cover the high nonlinearity for large strains. Lewis et al. [25] extended the approach to Mooney–Rivlin matrix material and an additional phenomenological ad-hoc volumetric density function to capture the plateau behavior for moderate porosity foamed rubber. They also included material aging and a gas pore effect [25]. A compressible hyperelastic isotropic porosity-based foam (CHIPFoam) model followed in [24]. This model additionally considered the compressibility of the matrix material. Modeling approaches concerning porous materials are manifold and many authors contributed, see e.g.  additionally [12, 20, 22, 27], reference within and others.

Here we present a material model which captures the essential properties of low porosity rubber. The model is based on the homogenized hollow sphere model of [10] in its original representation, and is extended to viscoelastic material properties. An additional numerical hollow sphere model serves as a comparison. Both, the elastic model and the viscoelastic extension, are limited to moderate deformations. The hollow sphere modeling captures the effect of pores and allows for arbitrary deformations as long as limit states like pore collapse or large strain softening are not reached. Correspondingly, the time-depended material response is modeled with a straightforward formulation derived from small-strain viscoelasticity. The approach is justified here because the analyzed rubber mixture only reproduces equal mechanical behavior for moderate deformations up to 40%. Above all, the model allows for a clear interpretation and attributes the experimentally measured properties to parameters like porosity, elastic modulus and relaxation time.

We start with a description of the material and the experiments in Sect. 2. The underlying elastic part for the closed-cell rubber is modeled in Sect. 3 via an homogenization approach and additionally by a packed hollow sphere ansatz. Both versions are compared numerically in Sect. 3.5. The extension to viscoelasticity is described in Sect. 4. The constitutive model is implemented in the finite element software Abaqus. Required details are given for further use in Sects. 3.3 and 4.1. Hereby, the equilibrium state is the basis for the hyperelastic validation in Sect. 5.1. The identification of viscoelastic parameters in Sect. 5.2 concludes the validation. All experimentally determined material parameters are based on the first loading. Consequently, further inelastic effects like Mullins-effect and permanent set as well as material aging or any chemical altering on long time scales are neglected.

2 Experimental investigations

Rubbers typically show a nonlinear mechanical behavior, which result in effects like creep, relaxation, equilibrium hysteresis, permanent set and material softening. These effects are becoming increasingly complex by considering higher filler content, e.g. carbon black or silica, or additional porosity. The modeling approach proposed in this work is motivated by the experimental findings of EPDM rubber and its foam, see Table 1.

Table 1 Recipe of investigated carbon filled EPDM resp. its foam (by ACR)
Full size table

Since in many applications inhomogeneous conditions are expected, the relaxation behavior is covered for different strains in multi-step relaxation experiments. In such experiments ordinary tensile specimens are loaded stepwise with a defined holding time. Caused by the holding time the material is able to relax the instantaneous stress response and reach an equilibrium stress state. Here the investigations are conducted for the pore-free material and the foam with specimens according to DIN 53504 standard [7]. The stretch is increased from \(\lambda _\text {min}=1.1\) gradually by \(\Delta \lambda =0.05\) up to the maximal stretch \(\lambda _\text {max}=1.4\) with a velocity around and an holding time of ten minutes. The pore-free EPDM hereby will be used as input for FEM model parameter identification, and the data of the foam are used to evaluate the model results. The measured stress-time behavior is displayed in Fig. 1. Since the modeling strategy is divided in an elastic and a viscoelastic part, the equilibrium stresses with respect to the initial configuration \(P_\infty \) (end point of each relaxation phase) of each experiment are extracted and used as input and validation set in Sect. 5, cf. Table 2. The direct comparison of the equilibrium stresses shows the effect of geometrical porous material structure.

Fig. 1
figure 1

Multi-step relaxation experiments on pore-free and foamed EPDM rubber

Full size image
Table 2 Equilibrium stress \(P_\infty \) of pore-free and foamed EPDM rubber material, cf. Fig. 1
Full size table

A material property of the foam, which later will serve directly as parameter, is the (initial) porosity of the foam. For comparison with experimental results, the initial porosity \(\phi \) of the examined rubber is determined with the following approximation, using the foamed rubber density \(\rho _\text {f}\) and the matrix density \(\rho _\text {m}\),

$$\begin{aligned} \phi \approx 1-\frac{\rho _\text {f}}{\rho _\text {m}}. \end{aligned}$$
(1)

Our foamed EPDM rubber has an initial porosity of \(\phi =0.34\). From microscopic observations we conclude that the rubber has mainly closed cells, see Fig. 2.

Fig. 2
figure 2

Microscope image of analyzed foamed rubber with porosity of \(\phi =0.34\)

Full size image

3 Elastic modeling of closed-cell rubber

The overall modeling approach is separated like the classical models for finite viscoelasticity in an elastic and viscoelastic part, see Sect. 4 for the latter. In this section we formulate the constitutive model for hyperelastic porous rubbers which is based on a deforming hollow unit sphere. At first, the analytical derivation by Danielsson et al. [10] for a homogenized hollow sphere of incompressible matrix material is presented. At next, the formulation is generalized to a representative hollow sphere volume of arbitrary material. The described approach is motivated by closed-cell foams but is indeed applicable for other foams as well, see e.g. [10, 32].

3.1 Deformation of a hollow sphere

The geometry in the initial configuration \(\mathcal {B}^0\) is completely defined by the porosity \(\phi \), the relation between pore volume \(V_\text {void}\) and total volume \(V_\text {total}\)

$$\begin{aligned} \phi&=\frac{V_\text {void}}{V_\text {total}} , \end{aligned}$$
(2)

respectively by the initial (dimensionless) pore radius , Fig. 3. The overall deformation \(\varvec{\varPhi }\left( \varvec{x}^0,t\right) \) is described by macroscopic stretches in the principal axes, \(\langle \lambda _i\rangle \), \(i=1,2,3\), so that the homogenized macroscopic deformation gradient \(\langle \varvec{F} \rangle \) has the components

$$\begin{aligned}{}[\langle \varvec{F} \rangle ]_{ij}=\langle \lambda _{i}\rangle \, \delta _{ij}, \qquad \left( \text {no sum on}\, i\right) . \end{aligned}$$
(3)
Fig. 3
figure 3

Undeformed and deformed unit hollow sphere

Full size image

The material point in the current configuration \(\mathcal {B}^t\) is at position \(\varvec{x}^t = \varvec{\varPhi }\left( \varvec{x}^0,t\right) \), whereby the deformation \(\varvec{\varPhi }\) has to be kinematically admissible at all points of the hollow sphere but is unknown otherwise. To specify it, we follow the works of Hou and Abeyaratne [18] and Danielsson et al. [10]. A vector-valued dimensionless function of the sphere radius R is introduced, \(\varvec{\varPsi }(R,t)\), and so the material position \(\varvec{x}^t\) can componentwise be written as

$$\begin{aligned} x^t_i&=\varPsi _i(R,t) \, x^0_i, \qquad \left( \text {no sum on}\, i\right) . \end{aligned}$$
(4)

The components of the unknown function \(\varvec{\varPsi }\) satisfy the condition

(5)

With the spherical coordinates in the reference configuration \(x^0_1 = R \, \sin \theta \, \cos \delta \), \(x^0_2 = R \, \sin \theta \, \sin \delta \), \(x^0_3 = R \, \cos \theta \), the components of the microscopic deformation gradient \(\varvec{F}\) follow as

$$\begin{aligned}{}[\varvec{F} ]_{ij}&=\frac{\partial x^t_i}{\partial x^0_j} =\frac{\text {d}\varPsi _i}{\text {d}R}\frac{x^0_i \, x^0_j}{R}+\varPsi _i \, \delta _{ij}, \qquad \left( \text {no sum on}\,i\right) . \end{aligned}$$
(6)

Assuming an incompressible matrix material, i.e. a locally volume preserving deformation with

$$\begin{aligned} \det \varvec{F}&= \varPsi _1 \, \varPsi _2 \, \varPsi _3+R \,\varPsi _1 \, \varPsi _2 \, \varPsi _3' + \frac{\left( x^0_1\right) ^2}{R} \, \varPsi _2\left( \varPsi _1' \, \varPsi _3-\varPsi _1 \, \varPsi _3' \right) + \nonumber \\&\qquad \frac{\left( x^0_2\right) ^2}{R} \, \varPsi _1\left( \varPsi _2' \, \varPsi _3-\varPsi _2 \, \varPsi _3' \right) =1 , \end{aligned}$$
(7)

a system of differential equations for the components of \(\varvec{\varPsi }\) results to

$$\begin{aligned} \varPsi _1' \, \varPsi _3-\varPsi _1 \, \varPsi _3'&=0 ,\nonumber \\ \varPsi _2' \, \varPsi _3-\varPsi _2 \, \varPsi _3'&=0 , \nonumber \\ \varPsi _1 \, \varPsi _2 \, \varPsi _3+R \, \varPsi _1 \, \varPsi _2 \, \varPsi _3'&=1 , \end{aligned}$$
(8)

where \(\left( \cdot \right) '\) denotes the differentiation with respect to R, cf. [18]. To solve the system (8) we make use of the solution for the spherical expansion of a hollow sphere. Here we have a homogeneous mapping \(x^t_i =\varPsi _o(R,t) \, x^0_i\) with

$$\begin{aligned} \varPsi _o =\left( 1+\frac{C^3}{R^3} \right) ^{\nicefrac {1}{3}} , \end{aligned}$$
(9)

and constant \(C=\left( a^3-a_0^3\right) ^{\nicefrac {1}{3}}\) characterizing the pore growth. The asymmetry of deformation is now denoted with positive constants \(\alpha _i\) as

$$\begin{aligned} \varPsi _i(R,t)&=\alpha _i \, \varPsi _o(R,t) . \end{aligned}$$
(10)

with condition \(\alpha _1 \, \alpha _2 \, \alpha _3=1\). With this ansatz the differential equations (8) are solved. The three constants \(\alpha _i\) can be determined from the boundary conditions on the outer surface, \(\langle \lambda _i(t)\rangle =\varPsi _i(1,t)\), which results in the coefficients

$$\begin{aligned} \alpha _i&=\langle \lambda _i\rangle \, \langle J \rangle ^{\nicefrac {-1}{3}}, \qquad \text {with}\quad \langle J \rangle =\langle \lambda _1\rangle \, \langle \lambda _2\rangle \, \langle \lambda _3\rangle , \end{aligned}$$
(11)
$$\begin{aligned} C&=\root 3 \of {\langle J \rangle -1 } , \quad \qquad \text {with}\quad \langle J \rangle > 1-\phi . \end{aligned}$$
(12)

In total, this gives the microscopic deformation

$$\begin{aligned} x_i^t = \langle \lambda _i\rangle \, \langle J \rangle ^{\nicefrac {-1}{3}} \left( 1 + \frac{\langle J \rangle -1}{R^3}\right) ^{{\nicefrac {1}{3}}} x_i^0, \qquad \left( \text {no sum on { i}}\right) \end{aligned}$$
(13)

and the microscopic deformation gradient \(\varvec{F}\) result from (6)

$$\begin{aligned}{}[\varvec{F} ]_{ij}&=\langle \lambda _i\rangle \, \langle J \rangle ^{\nicefrac {-1}{3}} \left( \frac{1-\varPsi _o^3}{\varPsi _o^2 \, R^2} \, x^0_i \, x^0_j+\varPsi _o \, \delta _{ij} \right) , \qquad \left( \text {no sum on { i}}\right) . \end{aligned}$$
(14)

Based on this, the right Cauchy-Green tensor \(\varvec{C}=\varvec{F}^T \varvec{F}\) follows

$$\begin{aligned}{}[\varvec{C} ]_{ij}&= \langle J \rangle ^{\nicefrac {-2}{3}} {\left\{ \begin{array}{ll} \hat{C} \left( x^0_i\right) ^2 + \langle \lambda _i\rangle ^2 \left( 2 \, \frac{1-\varPsi _o^3}{\varPsi _o \, R^2}\left( x^0_i\right) ^2+\varPsi _o^2 \right) , &{} \text{ for } i=j \\ \hat{C} \, x^0_i \, x^0_j + \frac{1-\varPsi _o^3}{\varPsi _o \, R^2} \, x^0_i \, x^0_j \left( \langle \lambda _i\rangle ^2+\langle \lambda _j\rangle ^2 \right) , &{} \text{ otherwise } \end{array}\right. } \nonumber \\ \text {with}\quad \hat{C}&= \left( \frac{1-\varPsi _o^3}{\varPsi _o^2 \, R^2}\right) ^2 \langle \lambda _k\rangle ^2 \left( x^0_k\right) ^2 , \end{aligned}$$
(15)

summation only on k. Now the microscopic first and second principal invariants \(\uppercase {i}\) and \(\uppercase {ii}\) are determined by

$$\begin{aligned} \uppercase {i}&=\langle J \rangle ^{{\nicefrac {-2}{3}}} \left( \langle \uppercase {i} \rangle \, \varPsi _o^2 +\frac{1}{R^2} \, \langle \lambda _i\rangle ^2 \left( x^0_i\right) ^2 \left( \frac{1}{\varPsi _o^4}-\varPsi _o^2\right) \right) , \end{aligned}$$
(16)
$$\begin{aligned} \uppercase {ii}&=\langle J \rangle ^{{\nicefrac {2}{3}}} \left( \frac{\left( \uppercase {ii}\right) }{\langle J \rangle ^2 \, \varPsi _o^2}+\frac{1}{R^2} \, \langle \lambda _i\rangle ^{-2}\left( x^0_i\right) ^2 \left( \varPsi _o^4-\frac{1}{\varPsi _o^2}\right) \right) , \end{aligned}$$
(17)

where \(\langle \uppercase {i} \rangle \) and \(\langle \uppercase {ii} \rangle \) are the corresponding macroscopic first and second principal invariants in principal axes

$$\begin{aligned} \langle \uppercase {i} \rangle&= \langle \lambda _1\rangle ^2+\langle \lambda _2\rangle ^2+\langle \lambda _3\rangle ^2 , \end{aligned}$$
(18)
$$\begin{aligned} \langle \uppercase {ii} \rangle&= \langle \lambda _1\rangle ^2\langle \lambda _2\rangle ^2 + \langle \lambda _2\rangle ^2\langle \lambda _3\rangle ^2 +\langle \lambda _3\rangle ^2\langle \lambda _1\rangle ^2 . \end{aligned}$$
(19)

These expressions for the deformation at every material point of the representative hollow sphere volume can be used to derive a porous rubber model.

We remark that this analytical approach allows for more than a purely spherical expansion of the pores but it is not meaningful for arbitrary deformations. For example, it does not account for pore collapse under compression. Also, the deformation becomes indefinite, if the macroscopic Jacobian \(\langle J \rangle \) reaches the value of porosity, see relations (11) and (12). In consequence, it does not lead to a polyconvex constitutive model.

3.2 Homogenization and macroscopic response

An incompressible elastic material can be modeled with a simple Neo–Hookean energy density

$$\begin{aligned} W_\text {NH} \left( \uppercase {i}\right) =\frac{\mu }{2}\left( \uppercase {i}-3\right) \, \end{aligned}$$
(20)

where \(\mu \) corresponds to the shear modulus at small strains. Assuming such an energy density for the matrix material, gives the homogenized energy density of the spherical shell model as

$$\begin{aligned} W =\frac{1}{V_\text {total}^{{0}}} \int \limits _{V_\text {matrix}^{{0}}} W_\text {NH}\left( \uppercase {i}\right) \,\mathrm {d}V . \end{aligned}$$
(21)

This evaluates for the unit sphere in the reference configuration with volume \(V_\text {total}^{{0}}=\nicefrac {4\pi }{3}\) , first invariant (16) and integral transformation to

$$\begin{aligned} W_\text {NHD} \big (\langle \uppercase {i} \rangle , \langle J \rangle , \phi \big )&= \frac{3\, \mu }{8 \, \pi } \int \limits _{a_0}^1 \int \limits _0^{2\pi } \int \limits _0^\pi \left( \uppercase {i}-3\right) R^2 \, \text {sin}\,\theta \, \,\mathrm {d}\theta \, \,\mathrm {d}\delta \, \,\mathrm {d}R \nonumber \\&=\frac{\mu }{2} \left[ \langle \uppercase {i} \rangle \, \langle J \rangle ^{{\nicefrac {-2}{3}}} \int \limits _{a_0}^1 2\, \varPsi _o^2 \, R^2 + \frac{R^2}{\varPsi _o^4}\,\mathrm {d}R -3 \left( 1-\phi \right) \right] \nonumber \\&= \frac{\mu }{2} \left[ A_1 \big (\langle J \rangle , \phi \big )\, \langle \uppercase {i}\rangle -3\left( 1 -\phi \right) \right] \\ \text {with}\quad A_1 \big (\langle J&\rangle , \phi \big )=2-\frac{1}{\langle J \rangle } - \frac{\phi +2\left( \langle J \rangle -1 \right) }{\langle J \rangle ^{{\nicefrac {2}{3}}} \left( 1 +\phi ^{-1}\left( \langle J \rangle -1\right) \right) ^{\nicefrac {1}{3}}}.\nonumber \end{aligned}$$
(22)

In that way we obtain an analytical expression for the macroscopic energy density of a Neo–Hookean material with spherical pores. For brevity we will name it here Neo–Hookean–Danielsson (NHD) model.

Similarly, when we assume the matrix material of the spherical shell to have a Mooney–Rivlin energy density

$$\begin{aligned} W_\text {MR} \left( \uppercase {i}, \uppercase {ii}\right) =\frac{\mu _1}{2}\left( \uppercase {i}-3\right) +\frac{\mu _2}{2}\left( \uppercase {ii}-3\right) \end{aligned}$$
(23)

then the homogenized energy density of the unit sphere follows, additionally with second invariant (17), as

$$\begin{aligned} W_\text {MRD}\big (\langle \uppercase {i} \rangle , \langle \uppercase {ii} \rangle , \langle J\rangle , \phi \big )&=\frac{\mu _1}{2} \left[ A_1\big (\langle J \rangle , \phi \big )\, \langle \uppercase {i}\rangle -3\left( 1 -\phi \right) \right] \nonumber \\&\qquad +\frac{\mu _2}{2} \left[ \langle \uppercase {ii} \rangle \, \langle J \rangle ^{{\nicefrac {-4}{3}}} \int \limits _{a_0}^1 2\, \frac{R^2}{\varPsi _o^2} + \varPsi _o^4 \, R^2 \,\mathrm {d}R -3 \left( 1-\phi \right) \right] \nonumber \\&=\frac{\mu _1}{2} \left[ A_1\big (\langle J \rangle , \phi \big )\, \langle \uppercase {i}\rangle -3\left( 1 -\phi \right) \right] + \nonumber \\&\qquad \frac{\mu _2}{2} \left[ A_2\big (\langle J \rangle , \phi \big )\, \langle \uppercase {ii}\rangle -3\left( 1 -\phi \right) \right] \nonumber \\ \quad \text {with}\quad A_2 \big (\langle J \rangle , \phi \big )=&-1+\frac{2}{\langle J \rangle }-\left( \frac{\langle J \rangle -1+\phi }{\langle J \rangle \, \phi } \right) ^{\nicefrac {1}{3}}\left( \frac{\phi +1}{\langle J \rangle }-1\right) . \end{aligned}$$
(24)

The corresponding model is named here Mooney–Rivlin–Danielsson (MRD) model. Compared with model (22), we have \(\mu _1+\mu _2=\mu \) here and the two parameters may allow for a better simultaneous fit to experimental data of different loading modes, e.g. uniaxial, planar and equibiaxial tension.

From the macroscopic strain energy density the first Piola–Kirchhoff stress tensor \(\varvec{P}\) follows as the work conjugate,

$$\begin{aligned} \varvec{P} = \frac{\partial W}{\partial \varvec{F}} \end{aligned}$$
(25)

and relates to the Cauchy stress tensor \(\varvec{\sigma }\) by \(\varvec{\sigma }=J^{-1}\varvec{P} \varvec{F}^T\), [16]. For an isotropic tensor function W the latter can directly be expressed as a function of the left Cauchy-Green tensor \(\varvec{B} = \varvec{F} \varvec{F}^T\) and so the Cauchy stress is given by the Rivlin-Ericksen representation theorem

$$\begin{aligned} \varvec{\sigma }=\frac{1}{J} \left[ \gamma _0 \, \varvec{I} + \gamma _1 \, \varvec{B} + \gamma _2 \, \varvec{B}^2\right] \end{aligned}$$
(26)

with the coefficients

$$\begin{aligned} \gamma _0 =J \, \frac{\partial W}{\partial J},\quad \gamma _1 =2 \left( \frac{\partial W}{\partial \uppercase {i}}+\uppercase {i} \, \frac{\partial W}{\partial \uppercase {ii}}\right) ,\quad \gamma _2 =-2 \, \frac{\partial W}{\partial \uppercase {ii}}. \end{aligned}$$

3.3 Finite element implementation

For numerical simulation the material model is implemented as a user material subroutine in the commercial FEM software Abaqus, [11]. In addition to the Cauchy stress tensor \(\varvec{\sigma }\) this requires its linearization for the Newton-Raphson method, i.e. the Eulerian stiffness tensor

$$\begin{aligned} \mathbb {K}=\frac{2}{J} \left( \frac{\partial (J \, \varvec{\sigma })}{\partial \varvec{B}} \, \varvec{B} \right) \cdot \!\!\!\! \cdot \; \mathbb {S}, \end{aligned}$$
(27)

here the symmetrizing fourth-order unit tensor \(\mathbb {S}\) is used, cf. [19]. The tensor product of two second-order tensors \(\varvec{A} \otimes \varvec{B}\) is defined componentwise as \(A_{ij}\otimes B_{kl} = C_{ijkl}\). Its customized part \((\varvec{A} \otimes \varvec{B})^{S_{24}}\) reads in components \((A_{ij}\otimes B_{kl})^{S_{24}}\) and is defined by \(\nicefrac {1}{4} \left( A_{il} \, B_{kj} + A_{ik} \, B_{lj} + A_{jl} \, B_{ki} + A_{jk} \, B_{li} \right) \), we refer to [19] for details of the representation. Then, the Eulerian stiffness tensor has the expression

$$\begin{aligned} \mathbb {K}&= \delta _1 \, \varvec{I} \otimes \varvec{I} + \delta _2 \left( \varvec{I} \otimes \varvec{B} + \varvec{B} \otimes \varvec{I} \right) + \delta _3 \left( \varvec{B}\otimes \varvec{I}\right) ^{\text {S}_{24}} + \delta _4 \, \varvec{B} \otimes \varvec{B} \nonumber \\&\qquad +\delta _5 \left( \varvec{I} \otimes \varvec{B}^2 + \varvec{B}^2 \otimes \varvec{I} \right) + \delta _6 \left( \varvec{B}\otimes \varvec{I}\right) ^{\text {S}_{24}} \varvec{B} + \nonumber \\&\qquad \delta _7 \left( \varvec{B} \otimes \varvec{B}^2 + \varvec{B}^2 \otimes \varvec{B} \right) + \delta _8 \, \varvec{B}^2 \otimes \varvec{B}^2 \end{aligned}$$
(28)

with the coefficients

$$\begin{aligned} \delta _1&= \frac{\partial W}{\partial J} + J \, \frac{\partial ^2 W}{\partial J^2} , \quad \ \delta _2=2\left( \frac{\partial ^2 W}{\partial \uppercase {i} \, \partial J } + \uppercase {i} \, \frac{\partial ^2 W}{\partial \uppercase {ii} \, \partial J } \right) , \quad \ \delta _3=\frac{2}{J} \, \gamma _1 ,\\ \delta _4&=\frac{4}{J} \left( \frac{\partial ^2 W}{\partial \uppercase {i}^2} + \frac{\partial W}{\partial \uppercase {ii}} + 2 \, \uppercase {i} \, \frac{\partial ^2 W}{\partial \uppercase {i} \, \partial \uppercase {ii} } + \uppercase {i}^2 \, \frac{\partial ^2 W}{\partial \uppercase {ii}^2} \right) ,\quad \delta _5 =-2 \, \frac{\partial ^2 W}{\partial \uppercase {ii} \, \partial J } ,\\ \delta _6&= \frac{4}{J} \, \gamma _2 , \quad \ \delta _7 = -\frac{4}{J} \left( \frac{\partial ^2 W}{\partial \uppercase {i} \, \partial \uppercase {ii} } + \uppercase {i} \, \frac{\partial ^2 W}{\partial \uppercase {ii}^2} \right) ,\quad \delta _8 = \frac{4}{J} \frac{\partial ^2 W}{\partial \uppercase {ii}^2} . \end{aligned}$$

3.4 Packed hollow sphere model

Additional to the analytical homogenization we will present here a more versatile model for closed pores in general materials undergoing arbitrary deformations. The model is also based on a hollow sphere, that is, the kinematic of deformation is as above and the major difference lies in the matrix material which now may be compressible, time-dependent and/or elasto-plastic. To describe a closed-cell material with pores of arbitrary size the hollow spheres can be combined to a space-filling packed hollow sphere (PHS) ensemble, cf. [29]. The homogenized porous material behavior follows from numerical evaluation of the PHS model. For computation an adapted Ritz-Galerkin method based on spherical harmonics, specialized quadrature rules, and exact boundary conditions is employed. Thus, on the microscale the stresses, strains and inelastic deformations are not uniform, see Fig. 4. The macroscopic stress state result by applying Hill’s averaging theorem.

Fig. 4
figure 4

Hollow sphere model: a real spherical harmonics \(Y_{25\,10}\) on a unit sphere, b discretization in radial direction and c hollow sphere deformed by a prescribed deformation; the colors indicate tension (red) and compression (blue) (color figure online)

Full size image

The shape functions \(N_a\) of the Ritz-Galerkin method are based on spherical coordinates R, \(\theta \) (polar angle) and \(\delta \) (azimuthal angle). They are calculated by a cartesian products of the piecewise polynomial radial functions \(R_r\) and the real spherical harmonics \(Y_{lm}\). This results for the material point \(x^t_i\) in

$$\begin{aligned} x^t_i\left( R,\theta ,\delta \right)&=\sum _a x^t_{ia} \, N_a\left( R,\theta ,\delta \right) =\sum ^{N_r}_{r=0} \sum ^{N_l}_{l=0} \sum ^{l}_{m=-l} x^t_{irlm} \, R_r\left( R\right) \, Y_{lm}\left( \theta , \delta \right) \end{aligned}$$
(29)

with the coefficients \(x^t_{ia}\) or \(x^t_{irlm}\), the number of radial elements \(N_r\) and the highest degree of spherical harmonics \(N_l\). The real spherical harmonics \(Y_{lm}\) are defined with the associated Legendre functions \(P^m_l\) by

$$\begin{aligned} Y_{lm}\left( \theta , \delta \right) = {\left\{ \begin{array}{ll} N_{l0} \, P^0_l \, \text {cos}\left( \theta \right) , &{} \text{ for } m=0 \\ [1Ex] \sqrt{2} \, N_{lm} \, P^m_l \, \text {cos}\left( \theta \right) \, \text {cos}\left( m \delta \right) , &{} \text{ for } m>0 \\ [1Ex] \sqrt{2} \, N_{l\vert m \vert } \, P^{\vert m \vert }_l \, \text {cos}\left( \theta \right) \, \text {cos}\left( \vert m \vert \delta \right) , &{} \text{ for } m<0 \end{array}\right. } , \end{aligned}$$
(30)

where \(N_{lm}\) is given by [34] \(N_{lm}=\sqrt{\nicefrac {\left( 2l+1\right) \left( l-m\right) !}{4 \pi \left( l+m\right) !}}\). In contrast to the commonly used finite element shape functions the real spherical harmonics represent finite rotations exactly, from which follows that the material symmetry is preserved [8, 29, 36].

The components of the macroscopic deformation gradient \(\langle \varvec{F} \rangle \) define the boundary of the sphere and are directly prescribed as [29]

$$\begin{aligned} x^t_{irlm}\left( \theta , \delta \right) = {\left\{ \begin{array}{ll} -\sqrt{\frac{4 \, \pi }{3}} \, [\langle \varvec{F} \rangle ]_{i2} , &{} \text{ for } \left[ r,l,m\right] =\left[ N_r,1,-1\right] \\ [1Ex] \sqrt{\frac{4 \, \pi }{3}} \, [\langle \varvec{F} \rangle ]_{i3} , &{} \text{ for } \left[ r,l,m\right] =\left[ N_r,1,0\right] \\ [1Ex] -\sqrt{\frac{4 \, \pi }{3}} \, [\langle \varvec{F} \rangle ]_{i1} , &{} \text{ for } \left[ r,l,m\right] =\left[ N_r,1,1\right] \\ [1Ex] \qquad 0 , &{} \text{ otherwise } \end{array}\right. } . \end{aligned}$$
(31)

With these boundary conditions the material’s response is calculated. The homogenized energy density follows from (21). The macroscopic first Piola–Kirchhoff stress tensor \(\langle \varvec{P}\rangle \) is the volume average of the numerically calculated microscopic first Piola–Kirchhoff stress tensor \(\varvec{P}\left( R, \theta , \delta \right) \),

$$\begin{aligned} \langle \varvec{P}\rangle =\frac{1}{V_\text {total}^{{0}}} \int \limits _{V_\text {matrix}^{{0}}} \varvec{P}\left( R, \theta , \delta \right) \, \text {d} V . \end{aligned}$$
(32)

In the following we model the sphere’s matrix material to be Neo–Hookean extended to the compressible range [1],

$$\begin{aligned} W^\text {comp}_\text {NH} \left( \bar{\uppercase {i}},J\right) =\frac{\mu }{2}\left( \bar{\uppercase {i}}-3\right) +\frac{K}{2}\left( J-1\right) ^2. \end{aligned}$$
(33)

Here \(\bar{\uppercase {i}}=J^{\nicefrac {-2}{3}}\uppercase {i}\) denotes to the incompressible part of the first invariant and so the first term is the same as in the Neo–Hookean energy density (20); K is the bulk modulus which controls the compression.

3.5 Comparison of NHD and PHS models

As described in the previous section the approaches for the analytical porous Neo–Hookean material (NHD model) and the numerically evaluated packed hollow sphere material (PHS model) differ slightly. Both presented models are able to describe the homogenized response of the representative volume of a hollow sphere. The main difference between the models is the matrix material. In the analytical derivation of the NHD model a completely incompressible matrix is presumed, whereas the PHS model allows for any matrix material. Here we consider an elastic matrix with energy density (33). For different levels of porosity \(\phi \in \left[ 0.05, \, 0.15,\, 0.3 \right] \) a pure volumetric expansion and a simple shear deformation are applied on the unit spheres, cf. [30]. The shear modulus is set to \(\mu = 2\) MPa. The ratio

$$\begin{aligned} \nu = \frac{3 \, K-2 \, \mu }{6 \, K+2 \, \mu } \end{aligned}$$
(34)

of model (33) is treated here like the small strain Poisson ratio and is varied to study the effect of compressibility for the PHS model.

In Fig. 5 the tensional stress state induced by a pure volumetric expansion up to \(J=4\) is shown. The negative pressure is normalized here with the shear modulus. The mesh parameters for the numerical evaluation of the PHS model result from a convergence analysis as \(N_r=58\) and \(N_l=1\). The curves for different porosities show that the numerical PHS model convergence to the NHD model for \(\nu \rightarrow 0.5\) as expected. We also see the typical non-monotonous stress-strain curves of an expanding spherical structure. The higher the spheres’ porosity, the lower is the deformation induced tension.

Fig. 5
figure 5

Pressure response to a pure volumetric expansion of hollow spheres with different porosities \(\phi \)

Full size image

At next, the hollow sphere is deformed under simple shear up to \(\gamma =\text {tan}\,10^\circ \). As before different values of porosity \(\phi \) and Poisson ratio \(\nu \) are considered. The mesh parameters for the PHS model are in \(N_r=19\) and \(N_l=3\). The normalized shear stresses for both models are displayed in Fig. 6 and show the same tendency as before. For compressive matrix material the induced stresses are smaller and the numerical PHS solution converges to the analytical NHD model for an increasing Poisson ratio (34). We remark that although simple shear is not a deformation in principal axes, the analytical model is still valid here.

Fig. 6
figure 6

Shear stress response to a simple shear deformation of hollow spheres with different porosities \(\phi \)

Full size image

4 Viscoelastic modeling

The time-dependency of the rubber’s mechanical behavior requires a viscoelastic model. Classic linear viscoelasticity uses the framework of rheological models with the small strain deformation \(\varepsilon \) and stress \(\bar{\sigma }\) cf. [2, 14]. The well-known generalized Maxwell model is described by the Young moduli \(E_j\) and the dynamic viscosities \(\eta \), as shown in Fig. 7. Its time dependent normalized relaxation function is defined by a Prony series form

$$\begin{aligned} g\left( t\right) =1+\sum _{j=1}^n g_j \, e^{\nicefrac {-t}{\tau _j}},\qquad \text {with}\quad g_j=\frac{E_j}{E_\infty }\quad \text {and}\quad \tau _j=\frac{\eta _j}{E_j} , \end{aligned}$$
(35)

where \(E_\infty \) represent the equilibrium modulus, in the linear theory constant, and the \(\tau _j\) denote the respective relaxation times. The full stress history \(\bar{\sigma }\left( t\right) \) in dependence of arbitrarily strain history \(\varepsilon \left( t\right) \) is defined by an alternative version of the integral formulation of linear viscoelasticity [13]

$$\begin{aligned} \bar{\sigma }\left( t\right)&= \int \limits _{0}^t g\left( t-s\right) \, \frac{\text {d}\bar{\sigma }_\infty \left( s\right) }{\text {d}s} \, \text {d}s ,\qquad \text {with}\quad \bar{\sigma }_\infty \left( s\right) =E_\infty \, \varepsilon \left( s\right) . \end{aligned}$$
(36)

For further details of the derivation we refer to “Appendix”. Combining this history integral with the generalized Maxwell model (35), the stress \(\bar{\sigma }\left( t\right) \) corresponds to a purely elastic stress response \(\bar{\sigma }_\infty \left( t\right) \) and a viscoelastic component described by history functions \(\bar{h}_j\) [21]

$$\begin{aligned} \bar{\sigma }\left( t\right)&= \bar{\sigma }_\infty \left( t\right) + \sum _{j=1}^n \bar{h}_j\left( t\right) , \nonumber \\ \text {with}\quad \bar{h}_j\left( t\right)&= \int \limits _{0}^t g_j \, e^{\nicefrac {-\left( t-s\right) }{\tau _j}} \, \frac{\text {d}\bar{\sigma }_\infty \left( s\right) }{\text {d}s} \, \text {d}s . \end{aligned}$$
(37)
Fig. 7
figure 7

Generalized Maxwell model for \(n+1\) rheological elements

Full size image

So far the viscoelastic representation is only valid for small and one-dimensional deformations. For our EPDM rubber model it is assumed that the viscoelastic behavior is on the one hand linear and on the other hand small in comparison to the elastic behavior \((g_j \ll 1)\). This corresponds to a normalized relaxation function g(t) which is independent of the deformation level. Thus, it follows for isotropic material the viscoelastic Cauchy stress tensor

$$\begin{aligned} \varvec{\sigma }\left( t\right)&= \varvec{\sigma }_\infty \left( t\right) + \sum _{j=1}^n \varvec{h}_j\left( t\right) , \nonumber \\ \text {with}\quad \varvec{h}_j\left( t\right)&= \int \limits _{0}^t g_j \, e^{\nicefrac {-\left( t-s\right) }{\tau _j}} \, \frac{\text {d} \varvec{\sigma }_\infty \left( s\right) }{\text {d}s} \, \text {d}s \, , \end{aligned}$$
(38)

where the equilibrium Cauchy stresses \(\varvec{\sigma }_\infty \) are defined by (26). We remark that this simplified relation is only applicable in the range of deformations, where the initial tangent modulus is valid.

4.1 Finite element implementation

For the numerical implementation an approximation of the history integrals is needed. Our at first one-dimensional approach follows the derivation of [21] and is originally based on [15, 37]. Considering the time interval \(\left[ t_{n}, t_{n+1}\right] \) with the time step \(\Delta t\) and \(t_{n+1}=t_{n} + \Delta t\) the integral splits into a known and unknown part

$$\begin{aligned} h_j\left( t_{n+1}\right) = e^{\nicefrac {-\Delta t}{\tau _j}} \, h_j\left( t_{n}\right) +g_j \int \limits _{t_n}^{t_{n+1}} e^{\nicefrac {-\left( t_{n+1}-s\right) }{\tau _j}} \, \frac{\text {d}\sigma _\infty \left( s\right) }{\text {d}s} \, \text {d}s. \end{aligned}$$
(39)

The latter can be computed recursively,

$$\begin{aligned} h_j\left( t_{n+1}\right)&=e^{\nicefrac {-\Delta t}{\tau _j}} \, h_j\left( t_{n}\right) + \nonumber \\&\qquad g_j \, \frac{1 - e^{\nicefrac {-\Delta t}{\tau _j}}}{\nicefrac {\Delta t}{\tau _j}} \left( \sigma _\infty \left( t_{n+1} \right) -\sigma _\infty \left( t_{n} \right) \right) \end{aligned}$$
(40)

so that the updating stress \(\sigma \left( t_{n+1}\right) \) is based on the history of the equilibrium response \(\sigma _\infty \left( t\right) \) [21]

$$\begin{aligned} \sigma \left( t_{n+1}\right) =\sigma _\infty \left( t_{n+1}\right) + \sum _{j=1}^n h_j\left( t_{n+1}\right) . \end{aligned}$$
(41)

The viscoelastic tangent modulus \(\nicefrac {\text {d}\sigma \left( t_{n+1}\right) }{\text {d}\varepsilon \left( t_{n+1}\right) }\) follows with Eqs. (40), (41) and the definition \(\sigma _\infty \left( t_{n+1}\right) = E_\infty \left( t_{n+1}\right) \varepsilon \left( t_{n+1}\right) \) to

$$\begin{aligned} E_{\Delta t}&= g_{\Delta t} \, E_\infty \left( t_{n+1}\right) , \qquad \text {with} \quad g_{\Delta t} = 1 + \sum _{j=1}^n g_j \, \frac{1 - e^{\nicefrac {-\Delta t}{\tau _j}}}{\nicefrac {\Delta t}{\tau _j}}, \end{aligned}$$
(42)

where the incremental normalized relaxation function \(g_{\Delta t}\) describes the relation between viscoelastic and pure elastic stiffness. The presented one-dimensional relations can simply be extended to the general three-dimensional stress tensor

$$\begin{aligned} \varvec{\sigma }\left( t_{n+1}\right)&=\varvec{\sigma }_\infty \left( t_{n+1}\right) + \sum _{j=1}^n \varvec{h}_j\left( t_{n+1}\right) ,\qquad \text {with} \nonumber \\ \varvec{h}_j\left( t_{n+1}\right)&=e^{\nicefrac {-\Delta t}{\tau _j}} \, \varvec{h}_j\left( t_{n}\right) + \nonumber \\&\qquad g_j \, \frac{1 - e^{\nicefrac {-\Delta t}{\tau _j}}}{\nicefrac {\Delta t}{\tau _j}} \left( \varvec{\sigma }_\infty \left( t_{n+1} \right) -\varvec{\sigma }_\infty \left( t_{n} \right) \right) , \end{aligned}$$
(43)

and the incremental stiffness tensor

$$\begin{aligned} \mathbb {K}_{\Delta t}= g_{\Delta t} \, \mathbb {K}\left( t_{n+1}\right) , \end{aligned}$$
(44)

with the elastic stiffness tensor \(\mathbb {K}\) from Eq. (27), cf. [21]. We remark that for pore-free rubber we compute corresponding results with standard ABAQUS models, which, however, have their limits too, cf. [3, 9, 38].

5 Parameter identification and model validation

The parameter identification is conducted stepwise. All mechanical model parameters are identified purely on the data of the pore-free EPDM. The porosity enters as a separate parameter from the considered foamed rubber. Here, first the hyperelastic parameters for the equilibrium state are determined and then the viscous parameters on the relaxation data. The hyperelastic and viscoelastic model approaches are adapted to the experimental data of the foam.

5.1 Equilibrium hyperelasticity: NHD and PHS models

As mentioned before, only the elasticity of the rubber matrix material needs to be determined. In order to predict the equilibrium states, the last points of the relaxation phase are considered, see Table  2. In Fig. 8 the technical stress P (the uniaxial first Piola–Kirchhoff stress) is displayed in dependence of the uniaxial stretch up to 40%. For the pure rubber matrix material, a least square fit of the Neo–Hookean model (20) leads to a parameter of \(\mu =3.15\) MPa (blue dashed line).

Fig. 8
figure 8

Parameter identification based on uniaxial data for Neo–Hookean model, and model predictions for the foam

Full size image

The computation in Abaqus is based on a single linear 3d continuum element, whose boundaries are moved and constraint corresponding to a uniaxial deformation. It can be seen that the stress-strain prediction of the analytical constitutive NHD model results in a stress overestimation (solid line). The numerical PHS model with \(\nu =0.17\) and the mesh parameter \(N_r=25\) and \(N_l=7\) fits the experimental data quite well (green dashed line). This underlines the effect of porosity in the matrix material. Clearly, the analytical derivation is idealized and small pores within matrix material, which may result from the vulcanization process, are neglected. Additionally, both models are not able to map the effect of (some) open cells on the overall rubber behavior. With these limits and using only one additional material parameter, namely the porosity \(\phi \), the analytical constitutive model (22) gives a sufficient approximation.

In general, for a reliable parameter identification uniaxial, equibiaxial and pure shear experiments are needed to consider the influence of arbitrary deformations. Therefore additional experiments on solid EPDM rubber have been performed. At this point we refer to the study [31], where the microscopic deformation in a representative foamed rubber element under uniaxial macroscopic tension up to \(\lambda = 2\) is investigated. Because of the cellular structure of porous rubber the deformation and stress state in a uniformly deformed specimen is not homogeneous. The deformation state at numerous microscopic positions was determined and the incompressible principal invariants \(\bar{\uppercase {i}}=J^{\nicefrac {-2}{3}}\uppercase {i}\) and \(\bar{\uppercase {ii}}=J^{\nicefrac {-4}{3}}\uppercase {ii}\) were derived, see left plot in Fig. 9. It was observed, that beside uniaxial stretch, especially pure shear deformations exist and that the local stretch may be up to 185% higher than the macroscopic uniaxial deformation.

Based on these results we have done additional experiments for uniaxial, pure shear and equibiaxial stretching tests to revise the parameter fit, which give similar but not identical values for the elastic parameter. We weight the measured data of the different experiments and obtain a revised parameter of \(\mu =3.19\,\)MPa for the Neo–Hookean material (20). Additionally, the matrix material model is extended by a second principal invariant, Eq. (23). The weighted combined parameter fit based on Mooney–Rivlin material yields \(\mu _1 =3.11\) MPa and \(\mu _2 =0.10\) MPa. Both matrix materials are now used for the analytical constitutive models (22) and (24), see right plot in Fig. 9. In the uniaxial stretch both models coincide and we note that the stress overestimation does not result from the parameter identification of the examined rubber material. The comparison of the quadratic error between simulation and experimental data yield that the simulation based on uniaxial data only, is the best one.

Fig. 9
figure 9

a Incompressible invariant plane of the microscopic state of deformation for macroscopic uniaxial deformation of foamed rubber (up to \(\lambda =2\)) [31]. b Elastic parameter validation based on weighted combined parameter fit

Full size image

5.2 Identification and prediction of viscoelasticity for foamed rubber

The viscoelastic parameters of the solid EPDM rubber material are determined on the basis of a multi-relaxation tensile test, explained in Sect. 2. They are derived from the first measured normalized relaxation function \(g\left( t\right) =\nicefrac {P\left( t\right) }{P_\infty }\), based on an elongation of \(\lambda =1.1\). This is necessary because all following relaxation curves do not refer to an initial zero stress condition. The coefficients \(g_j\) of the prony series in Eq. (35) are calculated with a least square method. A comparison of the error for various numbers of Maxwell elements gives the best result for \(n=4\), see Fig. 10. The corresponding values are presented in Table 3. We observe a small viscoelastic effect.

Fig. 10
figure 10

Viscoelastic parameter identification of pore-free EPDM rubber

Full size image
Table 3 Viscoelastic material parameter \(g_j\) and \(\tau _j\) of pore-free EPDM rubber
Full size table

For validation of our results we simulate the relaxation by a finite element computation in Abaqus. We use the viscoelastic NHD model with a matrix material based on \(\mu =3.15\) MPa and Table 3. In order to compare the measured normalized relaxation function with the viscoelastic model, the stress dependent material response, \(g\left( t\right) =\nicefrac {\sigma \left( t\right) }{\sigma _\infty }\), is calculated from an elongation of \(\lambda =1.1\), a hold time of ten minutes and a deformation rate of \(4\, \nicefrac {1}{\text {min}}\). We remark explicitly that the considered deformation of \(\lambda =1.1\) is in range which can be linearized; please compare with the measured behavior in Fig. 8. Therefore the introduced viscoelastic model is applicable. The results are presented in Fig. 11 and match the experimental data quite well. However for relaxation times lower than one second the stiffness results in an underestimation.

Because typically it is assumed that the foamed rubber has the same time dependence as its matrix material, we also fit the relaxation of the matrix material in Fig. 11. A comparison shows that the matrix material and the foamed rubber have equal viscoelastic properties for \(t>10\,\)s. For lower relaxation times the matrix stiffness is somewhat overestimated. This is caused by non-equal deformation rate during the experiments. In total, these results confirm the assumption that the matrix’ viscoelastic parameters are also valid for the foamed material.

Fig. 11
figure 11

Viscoelastic model prediction for the foam and comparison of pore-free and foamed EPDM rubber

Full size image

6 Summary

The presented viscoelastic model is able to predict the time dependent stress-strain relation for the observed EPDM foamed rubber as a first approximation. We separated the overall model in an hyperelastic and a viscoelastic part. The hyperelastic behavior is described by an homogenized hollow sphere constitutive model and the time-dependency by rheological elements. Both parts are implemented in the finite element software Abaqus.

All mechanical material parameter used here are based on multi-step relaxation experiments of pore-free material. In addition with the porosity of the foam, the model is completely defined. In each case the parameter identification is done with least square fits. The performed equilibrium simulation results in a stress overestimation, also under a optimized parameter identification based on different deformation states. It follows that the strict assumption of incompressible matrix material is not valid. However, we remark that the stiffness and the stress response would reduce, if the effect of possible open cells is considered.

In addition to the analytical hollow sphere model, we employ a numerical hollow sphere model to evaluate the effect of matrix compressibility. In order to compare both model predictions, we subjected the hollow sphere to a pure volumetric expansion and a simple shear deformation. As a result, we obtained that the numerical solution converges to the analytical model for an increasing Poisson ratio. It is noteworthy to notice that the analytical model is also valid for non principal axes deformations, although only the case of principal axes deformations is considered during the derivation. Compared to incompressible analytical model, the numerical model is able to determine the equilibrium stress response for the observed foamed rubber more accurately. This results from the considered matrix compressibility.

The viscoelastic prediction matches the experimental data of the foamed rubber quit well. Only in the case of lower relaxation times the stiffness is underestimated. Moreover a experimental comparison between pore-free and foamed rubber yields a good match. It follows that both materials, pore-free and foam, have nearly equal viscoelastic properties.

In future works the hyperelastic model should consider the effect of an open cell structure and the viscoelastic model should be extended to large deformations. This would allow a more general use of FEM simulations, cf. [22, 24, 39]. Furthermore the integration of additional microscopic parameter, like mean diameter or distribution of cells, may improve the equilibrium prediction for the foamed rubber.