Daletal 20 JMPS

See discussions, stats, and author profiles for this publication at: https://www.researchgate.
net/publication/344286291
An extended eight-chain model for hyperelastic and finite viscoelastic

response of rubberlike materials: Theory, experiments and numerical aspects
Article in Journal of the Mechanics and Physics of Solids · September 2020

DOI: 10.1016/j.jmps.2020.104159
CITATIONS READS
45 1,925
3 authors:
Hüsnü Dal Osman Gueltekin

Middle East Technical University Technische Universität Dresden
89 PUBLICATIONS 1,447 CITATIONS 18 PUBLICATIONS 502 CITATIONS
SEE PROFILE SEE PROFILE
Kemal Acikgoz
Middle East Technical University
8 PUBLICATIONS 166 CITATIONS
SEE PROFILE
All content following this page was uploaded by Hüsnü Dal on 17 September 2020.
The user has requested enhancement of the downloaded file.

An extended eight-chain model for hyperelastic and finite viscoelastic
response of rubberlike materials: Theory, experiments and numerical
aspects
Hüsnü Dala,∗, Osman Gültekina , Kemal Açıkgöza

a Department of Mechanical Engineering, Middle East Technical University, Dumlupınar Bulvarı No. 1
Çankaya, 06800, Ankara, Turkey
Paper dedicated to Professor Michael Kaliske on the occasion of his 60th birthday
Abstract
Rubberlike materials exhibit strong rate-dependent mechanical response which manifests itself in creep and relaxation tests as well as in the
hysteresis curves under cyclic loading. Unlike linear viscoelasticity, creep and relaxation response of rubber is nonlinear and amplitude-dependent.
Within this context, the contribution of this work is three fold. (i) On the experimental side, the characterization of equilibrium and non-equilibrium
responses are carried out by means of uniaxial and equibiaxial extension tests. Also performed are the creep and relaxation experiments at various
stress and strain levels. (ii) On the theoretical side, we extend the well-known eight-chain model via incorporating a simple yet instrumental tube-
constraint term composed of the second invariant into the non-affine network contribution reflecting the ground state equilibrium response. For
the non-equilibrium response, we propose a new evolution equation for the creep/relaxation behavior of rubberlike materials based on a relaxation
kinetics of a single polymer chain. The geometric non-linearity is incorporated via the finite deformation kinematics based on the multiplicative
split of the deformation gradient into elastic and viscous parts, whereas the volumetric and isochoric split of the deformation gradient is entirely
discarded. The rheology uses a nonlinear spring responsible for equilibrium elastice response in parallel to n number of Maxwell elements, leading
to a generalized Maxwell-Wiechert viscoelastic solid. (iii) The algorithmic implementation of the model features the spectral decomposition of
the respective terms and is demonstrated within the context of the finite element method. The developed model is validated by fitting both the
elastic and viscoelastic model responses with respect to the experimental data in the sense of uniaxial, (equi)biaxial extensions, and pure shear
tests. Relaxation and creep behavior of the model are thoroughly assessed. Also presented in the manuscript is the capability and the performance
of the model in the face of a relevant non-homogeneous boundary value problem. The quality of the findings earns the model vast utilization areas
from an engineering perspective.
Keywords: Finite viscoelasticity · rubber viscoelasticity · hysteresis · creep · relaxation · thermodynamical consistency · algorithmic setting ·
finite element modeling · experimental analysis
1. Introduction
Rubberlike materials exhibit strong rate-dependent response that earns rubber excellent mechanical properties
highly sought in applications suc as automotive, aerospace and soft robotics, just to name but a few. Viscoelasticity
improves the fracture toughness and energy absorption capacity rendering rubber very suitable as and shock absorbers
and vibration decouplers. The backbone of finite viscoelastictiy of rubberlike materials is an accurate description of
ground state equilibrium elastic response. In this context, the hyperelastic constitutive models for rubber can either be
represented by (i) phenomenological, or (ii) micromechanically based free energy functions. The phenomenological
models are based on invariant or principal stretch-based expressions for the free energy function that successfully
fit the experimental stress-strain curves. The invariant based models usually consist of polynomial combinations
[80, 62, 7, 89] or more sophisticated mathematical functionals [28, 76, 42, 27, 90, 30, 49, 66, 55, 40] of the first and
second invariants of the right/left Cauchy-Green deformation tensor. Various functional forms of free energy have been
proposed in terms of principal stretches e.g. [85, 4, 63, 71, 2, 5, 44]. Some material models have a mixed mathematical
representation in terms of principal stretches and invariants e.g. [32, 37]. Micromechanical models, on the other
hand, are physically-based models, condisdering the material microstructure by linking the macroscopic mechanical
∗ Corresponding author
Email address: [email protected] (Hüsnü Dal)
Preprint submitted to Journal of LATEX Templates September 17, 2020

behaviour to molecular structure of rubberlike materials through statistical mechanics. Herein, the force-displacement
relation of a single free-chain is obtained from a distribution function for the end-to-end distance of the chain by
making use of the concept of entropic elasticity. In this context, a Gaussian distribution function for a free chain with
infinite length was derived by Kuhn [46, 47]. Mooney [62] investigated the form of free energy function in a purely
mathematical context applicable for ideal rubber behaviour. Later on Treloar [80], using the Wall ’s assumptions [86],
and following Kuhn and Grün’s and Mooney’s form, derived the neo-Hooke model. Note that neo-Hooke model is
the simplest constitutive model that extends the concept of linear elasticity to the geometrically nonlinear setting.
Gaussian distribution, lacks the ability to predict stress-strain behaviour at large deformations, especially around the
chain extensibility limit [82]. In order to circumvent these drawbacks, non-Gaussian distribution functions based on
Langevin statistics are utilized in the sense of Kuhn and Grün [48]. Within this context, Wang and Guth [87] proposed
the three-chain model, Flory and Rehner [25] the four-chain model, Treloar and Riding [83] the full network model
and Arruda and Boyce [1] the eight-chain model. All these models make use of a fixed relation between micro-scale
and macro-scale deformations. Miehe et al. [60] proposed the non-affine micro-sphere model which extends the
full-network model of Treloar and Riding [83] by introducing a non-affine average network stretch and an additional
tube-constraint which accounts for the topological constraints. The notion of tube model or tube-constraints is based
on the fact that polymer chains are not free. In three dimensional network chains not lying in the same plane causes
topological constraints on one another. For this reason, the network models that rely on the free motion of the chains
usually underestimate the stress response under biaxial deformation upto moderately large stretches. The tube and
the extended tube model of Kaliske and Heinrich [32, 37, 41], the constraint-junction model of Flory and Erman
[22, 24] and the non-affine micro-sphere model of Miehe et al. [60] take into account the topological constraints
on the free motion of the polymer chains by introducing a tube-like constraint. For an extensive overview of the
aforementioned models, we refer to the review papers [15, 57, 75, 8]. The micro-sphere model shows an excellent
fitting performance to uniaxial and (equi)biaxial tests. However, its computational cost is considerably higher than
classical invariant-based formulations. The eight-chain model succesfully captures the uniaxial tensile behavior but
it lacks the ability of simultanous fitting of uniaxial and (equi)biaxial tests. It is a common practive to improve the
biaxial performance of the first invariant based models by incorporation of a second invariant term [69]. In this respect,
micromechanically motivated models based on the first invariant can be either improved by a second invariant based
term [41] or a principal stretch based term [19, 37, 66]. The extended tube model of Kaliske and Heinrich [37] and
the review paper of Boyce and Arruda [8] utilize a principle stretch based expression as a tube constraint, where the
former uses a Ogden-type term with negative power and the latter uses a constraint term proposed by Flory and Erman
[24]. Moreover, the second invariant based models lack the convenient functional form. A rigorous analysis which
establishes the relation between the tube-constraints and the second invariant has not been treated in the literature
yet. In order to resolve this issue, in this contribution, we propose an extended eight-chain model by reinforcing the
free energy function in terms of a tube constraint part. Within this context, micro-macro transition between the
micro-kinematic tube constraint and the second invariant is established and a simple yet instrumental second invariant
based term for the free energy function, which significantly improves the biaxial response of the eight-chain model, is
derived.
Experimental investigations of viscoelasticity in rubberlike materials are reported in, e.g., references [23, 52, 61,
59]. The fundamentals of linear and nonlinear viscoelasticity theory of polymers are reviewed in, e.g., [9, 65]. The
early treatment of rubber viscoelasticity is based on the finite linear viscoelasticity of Coleman and Noll [11] as an
extension of linear viscoelasticity. Most theories of finite viscoelasticity exploit the general axiom of a fading memory
introduced by Coleman [12], see also Truesdell and Noll [84]. Recent models of finite linear viscoelasticity use a
stress-type internal variables [52, 72, 33, 38]. For an excellent review, we refer to Haupt and Lion [31] and Miehe
and Göktepe [59]. Viscoelasticity is perceptible in the stress-strain response of dielectric and magnetorheological
elastomers [34, 70] as well. Internal variable type formulations of finite viscoelasticity based on the multiplicative
decomposition of deformation gradient are outlined in [56, 53, 6, 68, 43]. An alternative kinematic representation
of finite viscoelasticity based on a Lagrangian viscous metric tensor is proposed by Miehe and Keck [61]. The
algorithmic implementation of internal variable type formulations of finite viscoelasticity is shown in Reese and
Govindjee [68] and Dal and Kaliske [17] in the sense of Simo and Miehe [74] (for finite elastoplasticity) and Simo
[73] (for viscoplasticity). Therein, a predictor-corrector algorithm along with exponential mapping is employed for
the integration and solution of the evolution equation in a time-discrete setting, see [88, 13]. These approaches use the
unconditionally stable backward Euler scheme for the integration of the evolution equation in a time-discrete setting.
2
Eidel et al. [21] reported a considerable speed up in convergence by replacing backward Euler integration scheme
with higher order Runge-Kutta methods. Another approach for finite viscoelasticity of rubberlike materials is the
utilization of the micro-sphere kinematics where the three-dimensional viscoelasticity is obtained by integration of
one-dimensional rheological constructions via numerical quadrature over unit sphere. Within this context, Miehe and
Göktepe [59] proposed a finite viscoelasticity model depending on two micro-kinematic scalar internal variables on
discrete space orientations of the micro-sphere. In fact, two micro-mechanisms for the relaxation process of polymer
chains are introduced: The relaxation of superimposed entanglements and the release of topological constraints are
taken into account by a spectrum of nonlinear evolution laws in the logarithmic space of the discrete space orientations.
Extension of this approach to the viscoplastic response of uncured green rubber is achieved by Dal et al. [18].
Two molecular approaches for the description of the time-dependent behavior of polymer chains exist. Transient
network theory explains the stress relaxation phenomenon as a consequence of breakage and reformation of the cross-
links continuously. The theory was firstly proposed by Green and Tobolsky [29] and further developed by Lodge
[54], Phan-Thien [64] and Tanaka and Edwards [77, 78]. This approach was also adopted to the microsphere model
by Linder et al. [51]. Reptation-type tube models were developed for the definition of the motion of a single chain in
a polymer gel. The constraints on the free motion of a single chain are qualitatively modeled as a tube-like constraint
and the motion of the chain is described as a combination of Brownian motion within and reptational motion along the
tube. The model is proposed by De Gennes [26] and Doi and Edwards [20]. The reptational motion is successfully
adopted to finite viscoelasticity by Bergström and Boyce [6]. Two drawbacks of their approach can be observed
in the evolution equation. That is, the creep rate function γ̇ is singular at the onset of loading leading to an initial
asymptotic behavior in the relaxation curves which is not a characteristic behavior of the rubberlike viscoelastic
behavior. Incorporation of a perturbation parameter to the kinetic term helps to overcome the singularity in the
evolution law. However, the model becomes very sensitive to the perturbation parameter leading to an extra material
parameter in the evolution equation. Moreover, incorporation of a fictitious viscous network stretch λin to the evolution
equation is not physically conceivable and causes problems in the derivation of the consistent tangent necessary for
the finite element implementation. Recently, a rate and amplitude dependent finite viscoelastic model based on the
dynamic flocculation model for filled rubber has been proposed by Raghunath et al. [67].
Afore-mentioned facts motivate the derivation of a simple expression for the creep rate function based on the re-
laxation kinetics of a single chain. Within this context, we will (i) carry out amplitude-dependent creep and relaxation
tests in order to characterize the nonlinear viscous response of rubberlike materials, (ii) propose a micromechanically
motivated creep model based on the relaxation kinetics of a single polymer chain, (iii) put forward a new constitutive
model for the equilibrium and non-equilibrium response of rubberlike materials, and (iv) derive a fast and uncon-
ditionally stable update algorithm in a time-discrete setting along with the algorithmic treatment of the model into
finite element method. On the experimental side, the characterization of equilibrium and non-equilibrium response
of unfilled rubber is carried out. In addition to uniaxial and equibiaxial extension tests, creep and relaxation experi-
ments at various stress and strain levels are performed. On the numerical side, the developed model is validated with
respect to the experimental data presented in the manuscript and the well known experimental data by Treloar [81]
and Kawabata et al. [39]. Moreover, the performance of the new model is put through its paces via three dimensional
boundary value problems.
The paper is organized as follows: In Section 2, kinematics, fundamental maps of continuum mechanics, and
the governing equations of motion for an inelastic solid are briefly outlined. Section 3 presents a finite viscoelastic
constitutive framework. A new hyperelastic constitutive model, the extended eight-chain model is introduced for the
ground-state elastic response of rubberlike materials. A micromechanically motivated evolution equation based on
the relaxation kinetics of a single polymer chain is derived. Section 4 is concerned with the complete algorithmic
setting for the new viscoelastic model suitable for large strain finite element analysis. In Section 5, the experimental
methodology is outlined. The parameter identification and the validation of the constitutive model is demonstrated in
Section 6. In addition, the numerical performance of the finite element implementation is investigated through a three
dimensional representative boundary value problem. The manuscript closes with concluding remarks in Section 7.
2. Fundamental maps for a deformable solid and kinematics

This section introduces the field equations and corresponding state variables of an isotropic viscoelastic solid body
subjected to finite deformations. The kinematics and integrity basis of the deformation are briefly introduced.
3
2.1. Geometric mappings and the field variables
The mechanical response of a viscoelastic body is assumed to be described by two independent set of field vari-
ables, State(X, t) := {ϕ(X, t), I(X, t)}, where the internal field variables I(X, t) denote the micro-motion of the
free dangling chains past the ground elastic network of the polymer solid and the deformation map ϕ(X, t) represents
non-linear defprmation field
that transforms at time t ∈ T ⊂ R+ the material points X ∈ B0 onto spatial points x = ϕt (X) which is also
called the current/Eulerian configuration of the material point.Let TX B0 and Tx B symbolize the tangent spaces in
the reference/Lagrangian and current/Eulerian manifolds. The deformation gradient
F : TX B0 → Tx B; F := ∇X ϕt (X) (1)
maps the unit tangent of the reference/Lagrangian configuration onto its counterpart in the current/Eulerian config-
uration. Subsequently, let TX∗ B0 and Tx∗ B denote the co-tangent spaces in the Lagrangian and Eulerian manifolds.
Weighing cof[F ] by the volume map det[F ], one can describe a normal map per unit reference volume
F −T : TX∗ B0 → Tx∗ B; n = F −T N (2)
which maps the unit normal vector N of the reference/Lagrangian configuration onto its counterpart n in the cur-
rent/Eulerian configuration. The gradient operators ∇X [•] and ∇x [•] indicate the spatial derivative with respect to the
reference X and current x coordinates, respectively. In the sequel, let dX, dA and dV represent the infinitesimal
line, area and volume elements, in the undeformed configuration. Then, the deformation gradient F , its cofactor
cof[F ] = det[F ]F −T and the Jacobian J := det[F ] > 0 characterize the deformation of infinitesimal line, area and
volume elements
dV
dX dA
X X X
(a) (b) (c)
Figure 1: Three fundamental maps of a continuum: (a) The deformation gradient F as a mapping of an infinitesimal line element, (b) its cofactor
cof[F ] as an area map, and (c) and its determinant det[F ] as a volume map.
dx = F dX , da = cof[F ]dA , dv = det[F ]dV , (3)

see also Fig. 1. The condition J := det[F ] > 0 ensures the non-penetrable deformations ϕ. Furthermore, the
reference B0 and the spatial B manifolds are locally furnished with the covariant reference G and current g metric
tensors in the neighborhoods NX of X and Nx of x, respectively. These metric tensors are required for the mapping
between the co- and contra-variant objects in the Lagrangian and Eulerian manifolds [58]. In the sequel, the right
Cauchy Green tensor and the inverse of the left Cauchy Green tensors are defined
C = F T gF and c = F −T GF −1 (4)
as the pull back of the current metric g and push-forward of the Lagrangian metric G, respectively. The left Cauchy
Green tensor or the Finger tensor is denoted by b = c−1 . For a geometric interpretation, we refer to Fig. ??(a-b).
In order to comply with the principle of material objectivity and principle of material frame indifference, the energy
stored in an isotropic material is governed by principle stretches or invariants. Within this context, the eigenvalue
decomposition of the deformation metrics U and C are, respectively,
3
X 3
X 3
X
U := λa N a ⊗ N a , C := λ2a N a ⊗ N a and cof[C] := νa2 N a ⊗ N a , (5)
a=1 a=1 a=1
4
where
νi = J/λi such that ν1 = λ2 λ3 , ν2 = λ3 λ1 , ν3 = λ1 λ2 , (6)
are the areal stretches. Three isotropic invariants of the right Cauchy Green tensor read
I1 := tr[C], I2 := tr[cof[C]], and I3 := det[C] . (7)
The relation between the principal stretches and invariants can be established through
I1 = λ21 + λ22 + λ23 I2 = ν12 + ν22 + ν32 I3 = λ21 λ22 λ23 . (8)

In (7), we have introduced an alternative definition to the second invariant by postulating the equivalence of − tr[C 2 ] = 1
2 I12
tr[cof[C]]. Comparing (3) and (7), one can easily infer that the three isotropic invariants are associated with squares
of linear, areal and volumetric stretches of an infinitesimal cube representing a material point, see Fig. 2. To distin-
guish between volumetric and shear deformations, the deformation gradient F can be decomposed into volumetric
F vol := J 1/3 1 and unimodular F̄ := J −1/3 F parts
F = F vol F̄ . (9)
N2 N2
e3 N3 N3
ν1 a 2
0
2
ν 3a
0
0
0
λ1 a
λ1 a
e2 2
ν 2a
0
a0
a0 a0
a0 λ2 λ2
e1
a0 a0 a0
λ3 λ3
N1 N1
(a) (b) (c)
P
Figure 2: Geometric interpretation of the eigenvalue and eigenvectors of the right stretch tensor U = C 1/2 where U = 3a=1 λa N a ⊗ N a .
2 2 2 2 2 2
The invariants (a) I1 (C) = λ1 + λ2 + λ3 as sum of squares of the principal stretches, (b) I2 (C) = ν1 + ν2 + ν3 as sum of squares of principal
areal stretches where νi = J/λi and (c) I3 = J 2 = λ21 λ22 λ23 as square of volumetric stretch.
2.2. Multiplicative split of the deformation gradient and geometric mappings relative to intermediate configuration
In the sense of Kröner [45] and Lee [50], the deformation gradient can be multiplicatively decomposed into elastic
F e and viscous parts F v , i.e.
F = F eF v . (10)
The multiplicative decomposition of the deformation gradient requires the description of tangent T̄X̄ B and co-tangent
spaces T̄X̄∗ B in a fictitious intermediate configuration, respectively. Moreover, we describe a metric tensor G̃ :
T̄X̄ B → T̄X̄∗ B in the intermediate configuration which maps the unit vector in the tangent space onto its counterpart
in the co-tangent space. In the sequel, we introduce metric and deformation tensors relative to the intermediate
configuration as illustrated in Fig. 3. With the help of the intermediate configuration, one can describe the elastic right
Cauchy Green tensor and the elastic inverse left Cauchy Green tensors
C e = F eT gF e and ce = F e−T G̃F e−1 , (11)
as the pull-back of the current metric g to the intermediate configuration and the push-forward of the intermediate
metric G̃ to the current configuration, respectively. Moreover, it can be shown that the viscous metric C v = F T ce F
is the pull-back of the inverse left Cauchy Green tensor to the reference configuration.
5
2.3. Material and spatial rates of stress and deformation tensors
The material and the spatial velocity gradients are given as
Ḟ := ∇X ẋ and l = ∇x ẋ = Ḟ F −1 , (12)
respectively. By analogy with (10), one can describe the elastic and viscous parts of the spatial velocity gradient as,
le = F˙ e F e−1 and lv = F˙ v F v−1 . (13)
The relationships between l, le and lv can be established via the insertion of (10) into (12)2 yielding
v v
l = le + F e lv F e−1 or l = le + l̃ where l̃ := F e lv F e−1 . (14)
The spatial velocity gradient l can be decomposed into covariant symmetric and covariant skew-symmetric parts
engendering the rate of deformation and the vorticity tensors, i.e.
l =d+w where d = sym gl and w = skew gl , (15)
respectively. In line with the additive decomposition of the spatial velocity gradients given in (14), the symmetric rate
of deformation tensor d can also be additively decomposed into elastic and viscous (inelastic) parts
v
d = de + d̃ . (16)
Next, we emphasize the push-forward of the covariant (•)♭ and the contravariant (•)♯ objects
ϕ ∗ (•)♭ = F −T (•)♭ F −1 and ϕ ∗ (•)♯ = F (•)♯ F T . (17)
By the same token, the pull-back of the covariant (•)♭ and the contravariant (•)♯ objects read
ϕ ∗ (•)♭ = F T (•)♭ F and ϕ ∗ (•)♯ = F −1 (•)♯ F −T . (18)
The Lie derivative of an Eulerian object is the push-forward of the material time derivative of the corresponding
˙
Lagrangian object, i.e. ϕ ∗ {ϕ ∗ (•)}. Accordingly, the objective Lie derivative of the Kirchhoff stress tensor reads
£ν τ = τ̇ − lτ − τ lT , (19)
where τ̇ is the material or the total time derivative of the Kirchhoff stress tensor.
F F
TX B0 T̄X̄ B Tx B TX B0 T̄X̄ B Tx B
Fv Fe Fv Fe
C Ce g Cv G̃ ce = be−1
TX∗ B0 F v−T T̄X̄∗ B F e−T TX∗ B0 F v−T T̄X̄∗ B F e−T

Tx∗ B Tx∗ B
(a) F −T (b) F −T
Figure 3: Push-forward and pull-back relations of the current and intermediate metric tensors. (a) Pull-back of the current metric g to the
Lagrangian C = F T gF and the intermediate C e = F eT gF e configurations. (b) Pull-back of the intermediate metric G̃ to the Lagrangian
C v = F vT G̃F v and the its push-forward to the Eulerian configuration ce = F e−T G̃F −1 .
6
Ψe
ψ1ne γ̇1
σ σ
ψnne γ̇n
F ek F vk
F = F ek F vk
Figure 4: The generalized Maxwell-Wiechert viscoelastic solid with n elements in total for k ∈ {1 . . . n} represented in a geometrically nonlinear
setting. The nonlinear spring Ψe accounts for the ground state elastic response, while the nonlinear spring ψkne and the dashpot with the associated
creep rate γ̇k in each Maxwell element reflect the non-equilibrium (viscous) response of the material considered.
3. A compressible finite viscoelastic constitutive model for rubberlike materials
The emphasis is placed upon the formulation of constitutive relations with regard to the equilibrium and non-
equilibrium response of rubberlike materials in the sense of generalized Maxwell-Wiechert model under finite defor-
mations. Besides, a new evolution equation for the creep/relaxation rate based on the relaxation kinetics of a single
polymer chain is the subject matter of this section.
3.1. Rheology of the constitutive model
Unlike the classical treatment for rubberlike materials, where volumetric and isochoric response are strictly de-
coupled, in what follows, we adopt a specific form of free energy density function
Ψ(g, F , F e ) = Ψe (g, F ) + Ψne (g, F e ) = Ψe (C) + Ψne (C e ) , (20)
which reflects the standard viscoelastic solid, i.e. Maxwell elements are established parallel to the non-linear spring
characterizing the ground state elastic response. Hence, the non-equilibrium part of the free energy in (20) follows
from the rheological structure presented in Fig. 4 that describes the generalized Maxwell-Wiechert model with n
elements,
Xn Xn
Ψne (g, F e ) = ψkne (g, F ek ) or Ψne (C e ) = ψkne (C ek ) , (21)
k=1 k=1
where the index k indicates nothing but the Maxwell elements, i.e. k ∈ {1 . . . n} incorporated into the system.
Recalling (10), the multiplicative split of the deformation gradient for the non-equilibrium part now assumes distinct
deformation process, thereby being recast in the following form
F = F ek F vk . (22)
Accordingly, the total Kirchhoff and second Piola-Kirchhoff stress experienced by the viscoelastic material are given
by
τ := 2∂g Ψ(g, F , F e ) = τ e + τ ne and S := 2∂C Ψ(C, C e ) = S e + S ne , (23)
with the Kirchhoff stress terms τ e and τ ne characterizing the equilibrium (elastic) and the non-equilibrium (viscous)
stresses, respectively. By the same token, the second Piola-Kirchhoff stress terms S e and S ne indicate the equilibrium
(elastic) and the non-equilibrium (viscous) stresses, respectively.
7
3.2. Equilibrium part of the constitutive model: Extended eight-chain model
In view of (23), the equilibrium part of the constitutive response can be defined as follows
τ e := 2∂g Ψe (g, F ) and S e := 2∂C Ψe (C) . (24)
In the following sub-section, we will elaborate on the specific form of the constitutive response in the sense of Kirch-
hoff stress τ e .
3.2.1. Micro-molecular motivation for the extended eight-chain model

The entropy describes the available conformations of a single polymer chain. Following Miehe et al. [60], the
joint probability density of a single chain in a tube-like constraint environment can be multiplicatively split into two
independent events, i.e. p(λ, ν) = pf (λ)pc (ν), where pf and pc represent the probability densities due to free chain
response and tube-constraint parts, respectively. The entropic elasticity postulates
ψ(λ, ν) = −θη(λ, ν) with η = kB ln p(λ, ν) ❀ η(λ, ν) = ηf (λ) + ηc (ν) , (25)
where θ and η stand for the temperature and the entropy, respectively, while kB denotes the Boltzmann constant.
The entropy in (25) additively decomposes into free and constraint parts, i.e. ηf = kB ln pf and ηc = kB ln pc ,
respectively. Consequently, the free energy function of a single chain in a tube-like constraint environment reads
ψ(λ, ν) = ψf (λ) + ψc (ν) . (26)
The non-Gaussian probability density function for the free chain response by Kuhn and Grün [48] reads

−1 L−1 (λr ) λ
pf (λ) = p0 exp −N λr L (λr ) + ln −1
with λr = √ , (27)
sinh L (λr ) N
where L−1 (λr ) is√the inverse Langevin function of the relative stretch λr ∈ [0, 1) normalized with respect to the
extensibility limit N characterized by the number of segments N . The probability density function for the tube-like
constraint in the sense of [60, 20] is given as
" 2 #
r0 −1
pc (ν) = p0 exp −α ν , (28)
d0
√
where r0 = N l is the mean end-to-end distance of the chain, d0 is the tube diameter, α is a shape constant and
ν −1 = (d0 /d)2 is the tube areal contraction and ν = (d/d0 )2 is the tube areal stretch. Incorporation of (27) and (28)
into (25) leads to the free energy functions corresponding to the free-chain

−1 L−1 (λr )
ψf (λ) = N kB θ λr L (λr ) + ln + ψ0 (29)
sinh L−1 (λr )
and the tube-like constraint 2
l
ψc (ν) = αkB θN ν −1 + ψ0 , (30)
d0
respectively. The key issue in the modeling of macroscopic response of a polymer network lies in the non-affine
mean-kinematic variables with the affine macro-kinematic variables. Within this context, the macroscopic free energy
functions can be additively split into
Ψe (g, F ) := Ψef (g, F ) + Ψec (g, F ) (31)
the free-chain and tube constraint parts read

Z
1
Ψef (g, F ) = hnψf (λ̄)i and Ψec (g, F ) = hnψc (ν̄ −1 )i with h(•)i = (•)dV , (32)
|V |
V
8
where n is the volume-specific chain density and the λ̄ and ν̄ −1 are the micro-stretch and the micro-tube areal con-
traction, respectively. The homogenized response of the free energy of the macro-continuum can be assumed to have
the following simple form
hnψf (λ̄)i = nψf (λn ) and hnψc (ν̄ −1 )i = nψc (νn ) (33)
in terms of the average kinematic quantities λn and νn , see Fig. 5. At this point, the average network stretch and the
tube areal stretch are assumed as
r r r r
2 2 2 2 2 2
2 λ1 + λ2 + λ3 2 I1 3 ν1 + ν2 + ν3 3 I2
λn = = and νn = = . (34)
3 3 3 3
Herein, the network stretch is the mean-square root average of the principal stretches λi in the sense of [1]. The mean
non-affine tube areal contraction νn of the polymer network in constraint environment is linked to the macroscopic
areal stretches in terms of (34)2 . As a result, the free-chain and the tube-constraint part of the macroscopic free energy
function for the extended eight-chain model read

e L−1 (λr )
−1
Ψf (g, F ) = µN λr L (λr ) + ln and Ψec (g, F ) = µc (νn − 1) , (35)
sinh L−1 (λr )
√
where λr = λn / N represents the relative average network stretch. The shear modulus µ = nkB θ is associated with
the free-chain response and the µc = α(l/d0 )2 nkB θ is associated with the tube constraint part.
A = νn A 0
e3
A0 = d20 N2
N3
r0
λn
r0
0
λ1 a
a0 e2
λ 2a
0
a0
0
λ3 a
e1 a0 N1
Figure 5: Idealization of the rubber network by eight-chain model and the linkage between macroscopic stretches with mean streches of the
representative chain,Prespectively. N i for i = {1, 2, 3} represent the eigenvectors of the right Cauchy-Green tensor, whereas λi are the principle
stretches, i.e. C = 3i=1 λ2i N i ⊗ N i .
3.2.2. Compressible extended eight-chain model

The specific form of the compressible eight-chain model can be descibed as follows

χ µ 3N − 1 L−1 (λr )
Ψe (J, I1 ) = [(ln J)2 + (J − 1)2 ] − ln J + µN λr L−1 (λr ) + ln , (36)
4 3 N −1 sinh L−1 (λr )
where χ = κ − 23 µ and µ are the first Lamé constant and the shear modulus, respectively. The second expression in
the free energy function is the normalization term that enforces τ e = 0 at undeformed state. The stress expression for
the compressible eight-chain model reads
hχ i
τ e = 2∂g Ψe (J, I1 ) = (ln J + J(J − 1)) − µ̃ g −1 + µ̂(λn )b , (37)
2
9
where the two terms related to the shear modulus are given as
µ 3N − 1 µ 3N − λ2n
µ̃ = and µ̂(λn ) = . (38)
3 N −1 3 N − λ2n
For a Gaussian network where the number of segments approaches infinity, i.e. N → ∞, the parameters simplify to
µ̂ → µ and µ̃ → µ in (38) so that one recovers the compressible neo-Hookean model. Advancing a step forward from
(36), an extension to compressible eight-chain model can be achieved
Ψe (g, F ) = Ψev (g, F ) + Ψef (g, F ) + Ψec (g, F ) = Ψ̃ev (J) + Ψ̃ef (λr ) + Ψ̃ec (ν) , (39)
in conjunction with (35). Therein, the purely volumetric expression for the extended eight-chain model reads
χ µ 3N − 1 4
Ψev (g, F ) = Ψ̃ev (J) = [(ln J)2 + (J − 1)2 ] − ln J − µc ln J , (40)
4 3 N −1 9
in which the second and the third expressions are added in order to normalize the stresses at undeformed configuration.
Then, the stress expression for the compressible extended eight-chain model is recast into the following form,
2 µc
τ e = 2∂g Ψe (g, F ) = κ̂(J)g −1 + µ̂(λn )b + (I1 b − b2 ) , (41)
9 νn2
where
χ µ 3N − 1 4 µ 3N − λ2n
κ̂(J) = (ln J + J(J − 1)) − µ̃, µ̃ = + µc and µ̂(λn ) = . (42)
2 3 N −1 9 3 N − λ2n
The first expression in the stress equation (41) is the purely volumetric part, whereas the latter two expressions are
due to free chain response and the tube constraint, respectively. For various volumetric free energy functions and their
implementation into fnite element method in the quasi-incompressible setting, we refer to Kadapa and Hossain [36].
3.3. Non-equilibrium part of the constitutive model

In view of (23), the non-equilibrium part of the constitutive response can be described as follows
n
X n
X ne
τ ne := 2∂g Ψ(g, F e ) = τ̂ ne
k and S ne := 2∂C e Ψne (C e ) = Ŝ k . (43)
k=1 k=1
Therein, the individual overstress expression for each Maxwell element is given by
ne
τ̂ ne ne e ne e
k = 2∂g ψk (g, F k ) and Ŝ k = 2∂C e ψk (C k ) , (44)
according to the rheology illustrated in Fig 4. In the following sub-section, we will elucidate the individual contribu-
tions arisen from the particular selection of the free energy function ψkne (g, F ek ) in the sense of the Kirchhoff stress
τ̂ ne
k . Before giving an account of the constitutive modeling of the non-equilibrium response, we underline that the
subscript k reflecting the particular Maxwell branch in the rheological structure hereafter drops out of the relevant
equations, e.g., ψkne (g, F ek ), in order not to make the expressions too complicated. In the sequel, the effect of the tube
constraint effect on the non-equilibrium response of the rubber network is neglected. The compressible eight-chain
model for a single nonlinear elastic spring part of the Maxwell branch as follows,

χv µv 3N v − 1 β(λer )
Ψne (J e , I1e ) = ln J e + (J e − 1)2 − ln J e
+ µ v v
N λ e
r β(λe
r ) + ln . (45)
4 3 Nv − 1 sinh β(λer )
p p
Therein, λen = I1e /3 and λer = I1e /3N v denote the network and the relative network stretches for the non-
equilibrium deformations, respectively. β = L−1 (λer ) is the inverse Langevin function, µv is the shear modulus of
the non-equilibrium network, N v is the number of segments of the chains entangled around obstacles. Accordingly,
the Kirchhoff stress expression for the compressible eight-chain model is
v
χ
τ̂ ne = 2∂g ψ ne (J e , I1e ) = (ln J e + J e (J e − 1)) − µ̃v g −1 + µ̂v (λen )be , (46)
2
where the two terms related to the shear modulus read
µv 3N v − 1 µv 3N v − λe2
n
µ̃v = and µ̂v (λen ) = . (47)
3 Nv − 1 3 N v − λe2n
10
3.4. Thermodynamical consistency
The second axiom of thermodynamics restricts the proposed model by the so-called dissipation inequality
1
D := P − Ψ̇ ≥ 0 with P := S : Ċ = τ : d , (48)
2
being the stress power term. The reduced form of the dissipation inequality reads
v
Dred = τ̂ ne : d̃ ≥ 0 . (49)
The thermodynamical consistency is satisfied for the non-negative dissipation Dred . The evolution equation for the
v
viscous rate of deformation tensor d̃ must a priori satisfy the inequality (49).
3.5. Relaxation kinetics of a single chain
r0 r r
(a) (b) (c)
Figure 6: Stretch and relaxation of a single chain entangled around an obstacle: (a) undeformed equilibrium state, (b) non-equilibrium deformed
state after rapid stretch, (c) deformed and fully relaxed state.
In this part, we elucidate the relaxation kinetics of a single free chain [14]. To this end, we revisit the relaxation
kinetics of a single chain entangled around an obstacle at both sides as depicted in Fig. 6. In Fig. 6(a),
√ a single chain
entangled around obstacles at both ends is shown at rest. The end-to-end distance is given by r0 = N0 ˆl where N0
is the number of segments between the two obstacles and l̂ denotes the length of a monomer unit. If the polymer
chain is instantaneously stretched, see Fig. 6(b), it elongates and attains a more ordered conformation decreasing its
entropy and increasing the free energy, while the internal energy of the chain is kept constant, see Treloar [82] (p. 25).
The stretching of the chain triggers the movement of the free ends of the chains and they retract in combination
of reptational and Brownian motion. Hence, the√length of the polymer chain portion lying between the obstacles
elongates until the most favorable state, i.e. r = N∞ l̂, is reached. Fig. 6(c) shows not only the deformed but also
the fully relaxed state in this case with the same end-to-end distance r. Then, the number of segments at the fully
relaxed state reads
r
N∞ = λ̄2 N0 where λ̄ = . (50)
r0
Recall that unlike N0 , which is a material constant, N∞ is variable dependent on the instantaneous stretch level λ̄. Let
us consider an intermediate step between rapid deformation
p as seen in Fig. 6(b) and fully relaxed state as indicated in
Fig. 6(c). We define an end-to-end distance r̄(t) = N (t)l̂ which is probabilistically the most preferred state upon
removal of the load maintaining the end-to-end distance r at time t. In a deformation driven state, λ̄ is given and can
be decomposed into elastic and inelastic (viscous) parts as follows
r r r̄(t) r r̄(t)
λ̄ = = = λe (t)λv (t) where λe (t) := and λv (t) := . (51)
r0 r̄(t) r0 r̄(t) r0
Initially at time t = 0, the end-to-end distance reads r̄(0) = r0 , leading to λe (0) = λ̄ and λv (0) = 1. However, at the
fully relaxed state at time t = ∞, the end-to-end distance reads r̄(∞) = r, leading to λe (∞) = 1 and λv (∞) = λ̄.
For the evolution of the chain length in-between the two obstacles, we adopt a generic ordinary differential equation
(ODE), i.e.
1
Ṅ (t) = [N∞ − N (t)] , (52)
τ
11
where τ designates the relaxation time. Then, the solution of the ODE reads

−t
N (t) = (N∞ − N0 ) 1 − exp + N0 . (53)
τ
Next, we substitute (53) into (51)3 and consider the rate of the viscous stretch λv (t) with respect to time. This, upon
some simple mathematical manipulations, lead to the following
p ! !
v d r̄(t) d N (t)ˆl 1 1 N∞ p
λ̇ (t) = = √ = √ p − N (t) . (54)
dt r0 dt N0 ˆl 2τ N0 N (t)
Substitution of the relations stated below

√
e r N∞ p p
λ = = p and N∞ = λ̄ N0 (55)
r̄(t) N (t)
into (54), we readily obtain the rate equation

v 1 e 1 1 v e2
λ̇ = λ̄ λ − e = λ λ −1 . (56)
2τ λ 2τ
In order to comply with the large-strain kinematics in the Eulerian setting, the one-dimensional logarithmic counter-
part of the rate of deformation tensor is derived as follows,
λ̇v
ǫ̇v = = γ̇0 λe2 − 1 . (57)
λv
3.6. Evolution equation

Consistent with the reduced dissipation postulate (49), we propose an evolution for the inelastic rate of the defor-
mation tensor in the current configuration
v
d̃ := γ̇N , (58)
where γ̇ denotes the effective creep rate. Note that the thermodynamical consistency (49) is satisfied for γ̇ ≥ 0. The
flow direction is determined by the non-equilibrium Kirchhoff stress
τ ne √
N= with kτ ne k := τ ne : τ ne . (59)
kτ ne k
Nonlinear viscoelasticity is an energy activated process. In order to take this fact into account and to reconcile the
finite viscoelasticity with the relaxation kinetics derived above, an additional term (τ ne /τ̂ )m is added to the expression
in (57). Therein, the specific norm of the non-equilibrium stress τ ne is described as follows
kτ ne k
τ ne = √ . (60)
2
p
In what follows, we also replace the λe expression for a single chain by the network stretch λen = I1e /3 of the
eight-chain model where I1e = tr C e . Finally, we end up with the following function for the effective creep rate
m
τ ne
γ̇ := γ̇0 [λe2
n − 1] . (61)
τ̂
In the above expression the relaxation time 1/2τ is replaced by the creep rate constant γ̇0 > 0 for the sake of
convenience. Note that the term m ≥ 1 controls the energy-activated relaxation and leads to viscoplastic material
response for high values. The expression τ̂ is used merely for normalization purposes and assumed to be τ̂ = 1 [MPa]
12
in the subsequent treatments. In the end, we emphasize that the nonlinear viscous flow stated in (61) is governed by
one more additional material parameter in comparison to the that of the finite linear viscoelasticity and reflects from
the continuum point of view the behavior of a single dashpot in the Maxwell element as depicted in Fig. 4. Recall that
replacing the term γ̇0 [λe2
n − 1] → γ̇0 and taking m = 1, we obtain the classical multiplicative finite viscoelasticity
in the sense of Reese and Govindjee [68]. Moreover, taking N → ∞, one obtains the simplest finite viscoelasticity
formulation in the geometrically nonlinear setting based on multiplicative decomposition of the deformation gradient
into elastic and viscous parts.
R EMARK: The proper prediction of viscous flow direction requires γ̇0 [λe2 n − 1] ≥ 0. In incompressible rubber
viscoelasticity, the condition is always satisfied since the network stretch λen is always greater than unity λe2
n ≥ 1.
However, for slightly compressible elastomers, under purely hydrostatic pressure loading, it is possible to have λe2
n ≤
1. To resolve this issue, the term [λe2 e2
n − 1] must be replaced by |λn − 1| in the evolution equation (61) for the analysis
of compressible elastomers.
4. Algorithmic setting for the constitutive model
In this section, the algorithmic setting of the model suitable for a finite element implementation will be discussed.
In particular, the algorithmic framework for the Eulerian stresses and moduli terms of the equilibrium and non-
equilibrium responses will be established. Furthermore, we will propose an efficient implicit update for the principle
elastic stretches of the non-equilibrium Maxwell branches, conceptually following the references [14, 68, 17].
In a time-discrete interval, i.e. ∆t = tn+1 −tn , the current stresses and the tangent expressions need to be obtained
where ∆t stands for the time-increment. The time-discrete counterpart of the continuous expression in (19) reads
∆τ = Calgo : 12 £∆ϕg + ∇x (∆ϕ) τ + τ ∇Tx (∆ϕ) (62)
where ∇x (·) represents the spatial gradient of the term considered. Therein, Calgo is the total algorithmic moduli and
must be consistently derived for the accurate and quadratically convergent nonlinear finite element analysis based on
global iterative solution algorithms, e.g., the Newton-Raphson method. Subsequently, the total algorithmic moduli
term is also additively decomposed into equilibrium and non-equilibrium parts as follows
Calgo := Ce + Cne
algo . (63)
4.1. Stress and moduli terms for the elastic part

The computation of the equilibrium Kirchhoff stress in (41) requires the current value of b at time t = tn+1 which
is readily obtained at each iteration in the global iterative solution algorithm. In the subsequent treatment, we focus
C
on the elasticity moduli e in (63) and define
Ce := 4∂gg2 Ψe(J, I1 , I2 ) (64)
which gives rise to the following explicit form of the equilibrium tangent moduli through the particular form of the
free energy (39),
Ce = (ŝe + p̂e)g−1 ⊗ g−1 − 2p̂eI + n̂e b ⊗ b + fê(I1 b − b2) ⊗ (I1 b − b2) + ĝe(b ⊗ b − Ib ). (65)
Therein, the terms in regard to the Jacobian J reads
∂ 2 Ψe χ ∂ 2 Ψe χ 2
p̂e := J = [ln J + J(J − 1)] − µ̃ and ŝe := J 2 = J − ln J + 1) + µ̃ , (66)
∂J 2 2 ∂J 2 2
while the terms associated with I1 and I2 are written as
4µ 1 8µc 1 4µc 1
n̂e = , fê = − and ĝ e = . (67)
9N (1 − λ2r )2 81 νn5 9 νn2
13
The fourth-order tensor Ib is simply the push-forward of the fourth-order identity tensor I such that
1 ik jl 1 ik jl
Iijkl = δ δ + δ il δ jk and Ib ijkl = b b + bil bjk . (68)
2 2
In the above equation, the inverse Langevin function is computed by Padé approximation L−1 (λr ) ≈ λr (3−λ2r )/(1−
λ2r ) proposed by Cohen [10].
4.2. Stresses and moduli terms for the viscous part

Before we start our discussion, it needs to be mentioned that the non-equilibrium stress and moduli terms will be
calculated for a single Maxwell element, see Fig. 4, in order to rid the mathematical expressions of extra subscript k.
The computation of the non-equilibrium Kirchhoff stress in (46) requires the current value of be at time t = tn+1 . The
computation of be is contingent on the algorithmic treatment of the evolution equation (58) which will be elucidated
in the upcoming sub-sections.
4.2.1. Integration of the evolution equation:

In this sub-section, we begin with the introduction of the identity forming the basis of the subsequent development,
that is
1 v
− £ν be · be−1 = d̃ . (69)
2
The integration of the evolution equation (58) is based on the operator split of the material time derivative of be into
an elastic predictor (E) and an inelastic corrector step (I) such that
e
ḃ := lbe + be lT + £ν be . (70)
| {z } | {z }
E I
Both in (69) and (70) the Lie derivative of the elastic Finger tensor be is given which can be defined as follows
d
£ν be := F C v−1 F T . (71)
dt
During the elastic predictor step E (elastic trial step), d(C v−1 )/dt is equal to zero. Therefore,
Elastic predictor (trial step) E : C v−1

tr = C v−1
tn −→ betr = F C v−1 T
tn F . (72)
It should be highlighted that the subscript denoting the current time step tn+1 hereinafter drops out of the expressions
for convenience. The merit of this proposal is to circumvent the rate of inelastic metric in computation of the consistent
tangent moduli. Afterwards, we focus on the inelastic corrector step I. Here, the spatial velocity gradient l is set to
e
zero yielding ḃ = £ν be . Inserting this result into (69) and substituting the evolution equation proposed in (58) for
v
d̃ , we obtain the following expression,
e
Inelastic corrector step I : ḃ = −2γ̇Nbe tr , (73)
for the evolution equation for the inelastic corrector step. The interested reader is referred to references [14, 17],
regarding integration of the equations [72,73] in a time-discrete setting.
4.2.2. Algorithmic moduli for the non-equilibrium part:

In this part, we will set up a closed-form expression for the consistent tangent moduli for the non-equilibrium
(viscous) part. The spectral decomposition of the trial elastic deformation yields
3
X
F etr = λea, tr na ⊗ N a . (74)
a=1
14
Subsequently, we introduce a fictitious second Piola-Kirchhoff stress tensor
3
X
ne e−1 ne e−T ne
S̃ := F tr τ F tr where S̃ = s̃a N a ⊗ N a , (75)
a=1
in the intermediate configuration. In (75)2, the principal value of the fictitious second Piola-Kirchhoff stress reads
2
s̃a = τa /λea, tr . With the definition (75)1 at hand, the incremental rate equation can be defined as follows
∆S̃
ne
= Cne e
algo : ∆C tr where Cne
algo = 2∂C e
tr
S̃
ne
. (76)
Next, we consider (76)2 via exploiting (75)2 and obtain the moduli expression as follows
3 X
3 3
Cne cab − 2τa δab 2τa
X 2 X
algo = 2 2 N a ⊗ N a ⊗ ∂C etr (λeb, tr ) + 2 ∂C etr (N a ⊗ N a ) , (77)
a=1 b=1 λea, tr λeb, tr e
a=1 λa, tr
along with the respective partial derivatives with respect to the trial elastic right Cauchy-Green tensor,
3
1
Gab + Gba) ,
2 X
∂C etr (λeb, tr ) = N b ⊗ N b and ∂C etr (N a ⊗ N a ) = ( (78)
2(λea, tr
2
− λeb, tr )
2
b6=a
where we have introduced GabIJKL := MIK

a b
MJL a
+ MIL a
+ MJK . This expression in its turn consists of the following
M a := N a ⊗ N a where a
MIJ := NIa NJa . (79)
Upon replacing the derivatives in (77) by (78)1 and (78)2, we finally obtain the compact definition of the moduli
expression in the fictitious intermediate configuration
3 X
3 3 X
3
Cne cab − 2τa δab s̃a − s̃b
Gab + Gba) .
X X
algo = e2 e2
Ma ⊗ Mb + 2 2 ( (80)
a=1 b=1
λa, tr λb, tr a6=b b=1
2(λa, tr − λeb, tr )
e
5. Experimental Analyses
In this section, we will perform uniaxial and equibiaxial extension tests on nitrile butadiene rubber (NBR, see
Table 1 for composition) in order to obtain the hyperelastic and viscoelastic characteristics of the material. In the
Table 1: NBR rubber composition.
Ingredients (phr) Composition

Nitrile 3345
100
(Acrylonitrile 33% at Mooney viscosity, ML4, of about 45 at 373 K)
ZnO 5
Stearic acid 1
TMQ antioxidant 2
Antiozonant 2
Sulfur 1.5
TMTD (tetramethyl thiuram disulfide) accelerator 0.5
TMTM (tetramethyl thiuram monosulfide) accelerator 0.5
sequel, the relaxation and creep behavior of the same NBR material will be assessed by means of uniaxial extension
tests. These experimental studies are conducted on a modular Zwick/Roell Z010 universal testing machine with two
video extensometers dedicated specifically to uniaxial and equibiaxial tests. The TestXpert III software associated with
the test machine is used in implementing experimental protocols, measurements, and control of the testing machine.
Pincer grips, especially designed for rubber specimens, are mounted between the Zwick/Roell XForceP series load cell
with 500 N nominal force capacity.
15
5.1. Uniaxial and equibiaxial extension test on NBR rubber
The specimens are cut out from the same beige colored NBR rubber sheet with 2 mm thickness. The preparation
of the uniaxial test specimen complies with ISO37-Type2 geometry [35]. In order to retrieve the hyperelastic re-
sponse of the material, single cycle loading-unloading tests are performed at various crosshead speeds, i.e {10, 1, 0.1}
[mm/min]. It is observed that at/below 1 [mm/min] constant displacement rate, no hysteresis is observed in the spec-
imen. For the quasi-static investigation, a displacement controlled free test definition with 0.1 [mm/min] constant
displacement rate is chosen and maximum stretch λ = 2 is set via the test configuration options. The result of the
experiment is documented in Fig. 7(a) in the sense of stress (first Piola-Kirchhoff stress) versus stretch values. Con-
0.8 1
0.8
0.6
P11 [MPa]
P11 [MPa]
0.6
0.4
0.4
0.2 0.2 raw data
corrected curve
0 0
1 1.2 1.4 1.6 1.8 2 1 1.1 1.2 1.3 1.4
(a) λ (b) λ
Figure 7: (a) Uniaxial and (b) equibiaxial tension test results in terms of the first Piola-Kirchhoff stress versus stretch values.
cerning the equibiaxial experiments, a custom designed equibiaxial fixture is used. The specimen geometry required is
adopted from Axel Physical Testing Services [3]. The lower plate of the custom equibiaxial test apparatus is mounted
on the Zwick/Roell XForceP series load cell having 10 [kN] nominal force capacity, while its upper plate is fixed to
the cross-head and moves together with it. A virgin specimen is marked with four black dots and placed between 16
slider guided connection points of the upper plate. Steel wires connect the upper and lower plates through pulleys
and allow pulling the sliders as a result of cross-head movement, thereby applying equal forces at the periphery of
the specimen from the connection points. For a quasi-static investigation, a 1 [mm/min] constant displacement rate
is chosen from the test definition settings for both loading and unloading during a three-cycle experiment. The max-
imum stretch value is set to λ = 1.4, whereas the point of load removal at the end of a cycle is chosen to occur at 0
[MPa]. The result of the experiment is illustrated in Fig. 7(b). Within the third cycle, where λ = 1.306, the specimen
ruptured.
5.2. Creep and relaxation experiments on NBR rubber

Although the Mullins’ effect is negligibly small in unfilled rubber, a position controlled cyclic preconditioning
test at 6 [mm/min] constant displacement rate for three cycles is performed. The maximum stretch λmax = 2.25 is
chosen such that no failure occured below this value in the creep and relaxation tests. We wait around 30 minutes after
preconditioning and perform the relaxation tests. In the relaxation tests, the specimen is loaded to a target (ultimate)
stretch value with a crosshead speed 450 [mm/min] and held at this stretch value for 10 minutes. In order to quantify
the relaxation at different stretch levels, we conduct four experiments with four distinct target stretch values. The target
stretch values for the relaxation experiments are λ = 1.05, λ = 1.15, λ = 1.45, and λ = 2.20, respectively. Between
each experiment, we wait for around 30 minutes. The resulting data in terms of the first Piola–Kirchhoff stress versus
stretch curves are provided in Fig. 8. After the final relaxation experiment, we wait another 30 minutes and carry out
the creep tests. In the creep tests, the specimen is loaded to a target (ultimate) stress value as fast as possible and
this stress value is maintained for 10 minutes. Different stress levels become the subject matter of the experiments to
quantify the creep response for a broader spectrum, so we perform four experiments with distinct target stress values.
The maximum possible displacement rate of 450 [mm/min] is chosen to reach the target stress values. The specimen
is held at the target stress values for ten minutes. Then, the specimens are unloaded at 6 [mm/min] displacement rate.
16
1.2 1.2
1 1
λ = 1.05 λ = 1.05
0.8 λ = 1.15 0.8 λ = 1.15
P11 [MPa]
P11 [MPa]
λ = 1.45 λ = 1.45
λ = 2.20 λ = 2.20
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 100 200 300 400 500 600 0 20 40 60 80 100
(a) time [s] (b) time [s]
Figure 8: Relaxation experiments at different stretch levels, i.e. λ = 1.05, λ = 1.15, λ = 1.45, and λ = 2.2 for (a) the entire duration of the
experiments, (b) first 100 [s].
The target stress values for the creep experiments are, respectively, P11 = 0.1, P11 = 0.25, P11 = 0.5, and P11 = 1.
[MPa]. Between successive experiments, a 30 minutes waiting time is applied. The results of the creep experiments
is depicted in Fig. 9.
2.25 2.25
P11 = 0.10 P11 = 0.10

2 2
P11 = 0.25 P11 = 0.25
P11 = 0.50 P11 = 0.50
1.75 P11 = 1.00 1.75 P11 = 1.00
λ
1.5 1.5
1.25 1.25
1 1
0 100 200 300 400 500 600 0 20 40 60 80 100
Figure 9: Creep experiments at P11 = 0.1, P11 = 0.25, P11 = 0.5, and P11 = 1 in [MPa] for (a) the entire duration of the experiments, (b) first
100 [s].
6. Representative Numerical Analyses
The following sub-sections are concerned with some meticulous and highly erudite assessments of the proposed
extended eight-chain model in Sections 2-4 in terms of (i) the validation and the performance of hyperelastic and
viscoelastic frameworks under homogeneous uniaxial, biaxial, equibiaxial extensions and pure shear tests reported in
other studies as well as data obtained from our experiments portrayed in Section 5; (ii) the capability of the model
in describing viscoelastic mechanical response under cyclic loading by either referring to direct experimental mea-
surements conducted in Section 5 or representative geometries of relevant devices such as conical springs. To this
end, the in-house material and element subroutines are coded in the general-purpose finite element analysis program
FEAP 8.5 originally developed by Taylor [79] and is used as a simulation tool. Apart from that, the optimization tools
developed in-house for the nonlinear least squared analysis, is extensively used, see Dal et al. [15].
17
6.1. Comparison of the proposed extended eight-chain model with the eight-chain, Pucci-Saccomandi, and Carroll
models
The performance of the extended eight-chain model is compared with the eight-chain, Pucci-Saccomandi, and
Carroll models, well known 3-parameter models, in view of the classical Treloar [81] and Kawabata et al. [39]
datasets. Treloar data is used for the comparison of the uniaxial (UT), equibiaxial (ET), and pure shear (PS) responses,
whereas the Kawabata data is solely used to examine the biaxial (BT) response of the different models. In line with
Dal et al. [16, 15], the corresponding model fits are obtained with quality of fit values. In part (i) of Fig. 10, the
eight-chain model results are shown for uniaxial extension (UT) only, simultaneous (uniaxial (UT), equibiaxial (ET)
extension and pure shear (PS) tests combined), and biaxial extension (BT) tests with several initial stretch values for
λ1 . Likewise, the extended eight-chain model behavior is presented for the same tests in part (ii) of Fig. 10. The
retrieved material parameters are also listed in Table 2.
Table 2: Eight-chain and extended eight-chain model parameters fitted to the Treloar data and Kawabata (biaxial) data. The parameters in the first
row are obtained by fitting to the UT-data only. The second row parameters are obtained by simultaneous fit to the UT,ET and PS data.
Extended eight-chain model Eight-chain model

Parameter µ [MPa] N [-] µc [MPa] µ [MPa] N [-]
UT-fit (Treloar) 0.2608 25.392 0.2262 0.2672 25.59
Sim. fit (Treloar) 0.2518 25.014 0.4934 0.3280 28.64
Biaxial fit (Kawabata) 0.3123 147.59 0.219 0.4339 100
Experiment, Sect. 5.1) 0.3123 147.59 0.219 0.4339 100
Upon a closer look, the eight-chain model produces quite good results for UT and and PS under uniaxial loading
only. However, ET response of the model greatly underestimates the data points, as visualized in Fig. 10(i)(a). Though
the performance of the model is noticeably improved with slightly underestimated data points for ET in the case of
simultaneous fits, the results this time overestimate the data taken from UT and PS as observed in Fig. 10(i)(b). Such
a response is expected and has been previously observed by others which motivates the development of the extension
of the eight-chain model. Focusing on the data retrieved from Kawabata et al. [39], the predictions of the eight-chain
model agree quite well with the data for relatively small stretch values, i.e. λ1 = {1.04 . . . 1.24}, as depicted in
Fig. 10(i)(c). However, the predictions start to fall short of tracking the data for moderate stretch values and finally
exhibits an asymptotic behavior as the spectrum of the model response in terms of the first Piola-Kirchhoff stress
narrows down to a single response curve for large stretch values, see Fig. 10(i)(d).
Fig. 10(ii)(a) and Fig. 10(ii)(b) signify the substantial enhancement of the modeling quality associated with the
extended eight-chain model for both UT only and simultaneous tests. As matter of fact, the addition of the term
involving the second invariant I2 brings the modeling capability to such an extent that even a UT only tests seem
to be sufficient for the characterization of the homogeneous material behavior. Similarly, the new model are in
excellent conformity with the BT tests irrespective of the degree of the stretch, clearly displayed in Fig. 10(ii)(c) and
Fig. 10(ii)(d). In fact, the hyperelastic response of the material for small, moderate and larger stretches is excellently
traced by the model.
6.2. Assessment of the hyperelastic behavior of the extended eight-chain model

In order to obtain the hyperelastic material parameters of the extended eight-chain model, we perform a simulta-
neous nonlinear fitting procedure of the homogeneous uniaxial (UT) and equibiaxial (ET) extension outlined in Dal et
al. [15]. The UT and ET data are given equal importance during the parameter identification procedure. The attained
material parameters are listed in Table 3 and the corresponding fits against the experimental data is obtained follow-
ing the protocol in Section 5.1, see Fig. 12. The harmony of the model results with the experimental data is salient,
demonstrating the excellent performance of the proposed model. Before closing this part, we highlight that the mate-
rial parameters identified here is used in most of the upcoming numerical investigations in order to reflect the ground
state elastic mechanical response. Though the Lamé’s first parameter χ is not among the parameters fitted through
the nonlinear least-squares analysis, the value used in the forthcoming non-homogeneous finite element analyses is
incorporated into the Table 3.
18
(i) UT only fit Simultaneous UT, ET and PS fits
7 7
6 6
5 5
P11 [MPa]
P11 [MPa]
4 4
3 3
2 2
1 1
0 0
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8
(a) (b)
λ λ
BT fit with λ1 : 1.04 − 1.24 BT fit with λ1 : 1.3 − 3.7
0.4 1.4
0.3 1.05
P22 [MPa]
P22 [MPa]
0.2 λ1 = 1.04 0.7 λ1 = 1.3

λ1 = 1.06 λ1 = 1.6
λ1 = 1.08 λ1 = 1.9
λ1 = 1.10 λ1 = 2.2
λ1 = 1.12 λ1 = 2.5
0.1 λ1 = 1.14 0.35 λ1 = 2.8
λ1 = 1.16 λ1 = 3.1
λ1 = 1.20 λ1 = 3.4
λ1 = 1.24 λ1 = 3.7
BT model fit BT model fit
0 0
0.9 1 1.1 1.2 1.3 0.5 1.25 2 2.75 3.5
(c) λ2 (d) λ2
(ii) UT only fit Simultaneous UT, ET and PS fits
7 7
6 6
5 5
P11 [MPa]
P11 [MPa]
4 4
3 3
2 2
1 1
0 0
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8
(a) λ (b) λ
0.4 1.4
0.3 1.05
P22 [MPa]
P22 [MPa]
0.2 λ1 = 1.04 0.7 λ1 = 1.3

λ1 = 1.06 λ1 = 1.6
λ1 = 1.08 λ1 = 1.9
λ1 = 1.10 λ1 = 2.2
λ1 = 1.12 λ1 = 2.5
0.1 λ1 = 1.14 0.35 λ1 = 2.8
λ1 = 1.16 λ1 = 3.1
λ1 = 1.20 λ1 = 3.4
λ1 = 1.24 λ1 = 3.7
0 0
0.9 1 1.1 1.2 1.3 0.5 1.25 2 2.75 3.5
(c) λ2 (d) λ2
UT data points ET data points PS data points
UT model fit ET model fit PS model fit
Figure 10: (i) Eight-chain, and (ii) extended eight-chain model predictions for (a) uniaxial tension, (b) combination of uniaxial, equibiaxial, and
pure shear loadings using Treloar data (c) biaxial tension loading for λ1 : 1.04 − 1.24, (d) biaxial tension loading for λ1 : 1.3 − 3.7 using
Kawabata data. 19
(i) UT only fit Simultaneous UT, ET and PS fits
7 7
6 6
5 5
P11 [MPa]
P11 [MPa]
4 4
3 3
2 2
1 1
0 0
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8
(a) (b)
λ λ
0.4 1.4
0.3 1.05
P22 [MPa]
P22 [MPa]
0.2 λ1 = 1.04 0.7 λ1 = 1.3

λ1 = 1.06 λ1 = 1.6
λ1 = 1.08 λ1 = 1.9
λ1 = 1.10 λ1 = 2.2
λ1 = 1.12 λ1 = 2.5
0.1 λ1 = 1.14 0.35 λ1 = 2.8
λ1 = 1.16 λ1 = 3.1
λ1 = 1.20 λ1 = 3.4
λ1 = 1.24 λ1 = 3.7
0 0
0.9 1 1.1 1.2 1.3 0.5 1.25 2 2.75 3.5
(c) λ2 (d) λ2
(ii) UT only fit Simultaneous UT, ET and PS fits
7 7
6 6
5 5
P11 [MPa]
P11 [MPa]
4 4
3 3
2 2
1 1
0 0
1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8
(a) λ (b) λ
0.4 1.4
0.3 1.05
P22 [MPa]
P22 [MPa]
0.2 λ1 = 1.04 0.7 λ1 = 1.3

λ1 = 1.06 λ1 = 1.6
λ1 = 1.08 λ1 = 1.9
λ1 = 1.10 λ1 = 2.2
λ1 = 1.12 λ1 = 2.5
0.1 λ1 = 1.14 0.35 λ1 = 2.8
λ1 = 1.16 λ1 = 3.1
λ1 = 1.20 λ1 = 3.4
λ1 = 1.24 λ1 = 3.7
0 0
0.9 1 1.1 1.2 1.3 0.5 1.25 2 2.75 3.5
(c) λ2 (d) λ2
UT data points ET data points PS data points
UT model fit ET model fit PS model fit
Figure 11: (i) Carroll, and (ii) Pucci-Saccomandi’s model predictions for (a) uniaxial tension, (b) combination of uniaxial, equibiaxial, and pure
shear loadings using Treloar data (c) biaxial tension loading for λ1 : 1.04 − 1.24, (d) biaxial tension loading for λ1 : 1.3 − 3.7 using Kawabata
data. 20
Table 3: Ground state elastic material parameters of the extended eight-chain model simultaneously fitting the experimental data reported in
Section 5.1.
Extended eight-chain model parameters

Parameter χ [MPa] µ [MPa] N [-] µc [MPa]
Simultaneous fits (UT and ET) 105 0.3448 20.001 1.3497
0
P11 [MPa]
0
0
UT model fit
0 UT exp
ET model fit
ET exp
0
1 1 1 1 1 2
λ
Figure 12: Extended eight-chain model fits against the data received from UT and ET experiments performed on NBR rubber.
6.3. Assessment of the viscoelastic behavior of extended eight-chain model in terms of relaxation & creep
We now assess how well the model proposed fits the experimental data. For this purpose, the model is simultane-
ously fitted against the data obtained from homogeneous tensile uniaxial tests characterizing relaxation experiments
with ultimate stretches λ = 1.05, λ = 1.15, λ = 1.45, and λ = 2.20 as seen in Fig. 8 and creep experiments with
ultimate stresses P11 = 0.1, P11 = 0.25, P11 = 0.5, and P11 = 1 [MPa], see Fig. 9. For the baseline elastic response,
the parameters identified in Table 3 is used. Firstly, we try to fit the model response against both data concurrently
by means of a single Maxwell branch. However, all the attempts invariably fall flat that we cannot achieve model
fits satisfactorily close to the experimental data. Therefore, one more Maxwell branch is added into the rheological
structure of the model and the simultaneous fitting procedure in the sense of trial and method starts again. As can
be observed in Fig. 13 and Fig. 14, both relaxation and creep response of the viscoelastic model subjected to several
target stretch and stress conditions under uniaxial extension are traced very well. The respective viscous material
parameters are reported in Table 4. Also shown in Fig. 13(b) and Fig. 14(b) are the zoom-in view of both relaxation
and creep responses during the first 100 [s], respectively. Apparently, the slight deviations from data observed in
relaxation, i.e. at λ = 1.45 and λ = 2.20, and the creep curves, i.e. at P11 = 0.5, and P11 = 1 [MPa], vanish as
the relaxation and creep processes continues to evolve in time. It needs to be highlighted that a fully incompressible
viscous free energy, i.e. χv = 0 [MPa], which appears in (45) may yield jumps in the relaxation and creep behavior
upon the descent and ascent of the curves, respectively. For this reason, a sufficiently small value for χv is hereinafter
considered throughout the examples, see Table 4. In addition, the segment number of the non-equilibrium network is
taken equal to the that of the equilibrium network, i.e. N v ≡ N , due to the lack of the physical information.
Table 4: Viscoelastic material parameters in regard to Maxwell branches used.
The eight-chain model (viscous part) & evolution equation

Parameter χv [MPa] µv [MPa] N v [-] γ˙0 /τ̂ m [MPa−m s−1 ] m [-]
st
1 Maxwell branch 1 0.3448 20.001 50 1
2nd Maxwell branch 1 0.0517 20.001 1 3
21
1.2 1.2
1 model fit 1 model fit

λ1 = 1.05 λ1 = 1.05
0.8 λ1 = 1.15 0.8 λ1 = 1.15
P11 [MPa]
P11 [MPa]
λ1 = 1.45 λ1 = 1.45
λ1 = 2.20 λ1 = 2.20
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 100 200 300 400 500 600 0 20 40 60 80 100
Figure 13: Model predictions against the experimental data of relaxation tests at λ1 = 1.05, λ1 = 1.15, λ1 = 1.45, and λ1 = 2.20 for a) the
entire duration of the experiment (600 [s]); b) the first 100 [s].
2.25 2.25
2 model fit 2 model fit

P11 = 0.10 P11 = 0.10
P11 = 0.25 P11 = 0.25
1.75 P11 = 0.50 1.75 P11 = 0.50
λ1 [-]
λ1 [-]
P11 = 1.00 P11 = 1.00

1.5 1.5
1.25 1.25
1 1
0 100 200 300 400 500 600 0 20 40 60 80 100
Figure 14: Model predictions against the experimental data of creep tests at P11 = 0.1, P11 = 0.25, P11 = 0.5, and P11 = 1 [MPa] for a) the
entire duration of the experiment (600 [s]); b) the first 100 [s].
6.4. Non-homogeneous cyclic loading of a conical spring

It is of utmost importance to decipher the range of period of loading as it may well affect the design considerations
such as base insulations of structures or vehicle industry, to name but a few. To this end, we demonstrate the utility of
the constitutive approach by simulating non-homogeneous tension-compression tests on a realistic tailored geometry
of a conical spring made of NBR rubber (light pink color) and stainless steel (dark red color) parts as depicted
in Fig. 15(a). Conical springs are utilized, e.g., in trains etc., alternatively to standard compression springs and
each coil nests inside the one before the other, leading to a cone-shaped geometry. The mesh generated consists of
N UMEL= 5300 eight-noded mixed Q1P0 brick elements connected by N UMNP= 6200 nodes, see Fig. 15(b). The
purely volumetric part of the elastic free energy function stated in (40), χ[(ln J)2 + (J − 1)2 ]/4 is formulated at
element level in the mixed Q1P0 finite element formulation. The tailored geometry basically features the following
radii, such as r1 = 50, r2 = 64, r3 = 80, r4 = 100 and r5 = 72 [mm], while the heights are h1 = 35, h2 = 55,
h3 = 75 and h4 = 125 [mm] as indicated in Fig. 15(c). Also shown is the load acting on the top surface of stainless
steel in accordance with Figs. 16(a), where cyclic displacement is applied.
The elastic and viscous parameters for the NBR rubber are basically obtained from Table 3 and Table 4 with the
Lamé’s fist parameter χ = 104 [MPa] being the sole exception, while the first Lamé parameter and the shear modulus
22
z uz (t) fz (t) z
y r5
x
x
r1 h1 h2 h3 h4
r2
r3
rubber
r4
steel
(a) (b) (c)
Figure 15: (a) Schematic view of a conical spring composed of two materials, namely NBR rubber and stainless steel, (b) Discretization with
N UMEL= 5300 eight-noded mixed Q1P0 brick elements connected by N UMNP = 6200 nodes. (c) The tailored geometry involves the following
radii, i.e. r1 = 50, r2 = 64, r3 = 80, r4 = 100 and r5 = 72 [mm], while the heights are h1 = 35, h2 = 55, h3 = 75 and h4 = 125 [mm].
The domain is subjected to displacement controlled load based on Fig. 16(a).
for the essentially elastic mechanical response of the stainless steel are chosen as χ = 104 [MPa] and µ = 76923
[MPa]. Based on our experience during simulations we note that assigning an identical value for the first Lamé
parameter χ in the NBR rubber and the stainless steel seems indispensable in order to hinder the numerical stability
issues due to the mismatch in the volumetric constraint.
6.4.1. Cyclic loading in a displacement driven context:

Based on the afore-stated geometrical, mesh, and the constitutive parameter data a cyclic sinusoidal displacement
uz (t) is applied on the nodes located at the top of the geometry by using three distinct periods T = 1, T = 10 and
T = 100 [s] along z-direction as depicted in Fig. 16(a). The displacement function uz (t) ranges between 20 and −20
[mm]. Fig. 16(b) addresses a remarkable contrast in the forces Fz experienced by the top surface with respect to the
three distinct periods used during loading cycles. Upon a closer look, the viscous stiffening effect due to the higher
loading rate is palpable which is connected with the high frequency and low period of loading. Besides, the result of
Fig. 16(c) substantiates the findings in Fig. 16(b) with the stiffest response provided by T = 1 [s]. Being an important
20 18
tA tC uz (t)
15 12
10 6
Force Fz [kN]
uz [mm]
5 0
0 −6
T = 1 [s]
−5 −12 T = 10 [s]
−10 −18 T = 100 [s]
−15 −24
tB tD
−20 −30
0 1 2 −20 −10 0 10 20
Number of cycles uz [mm]
(a) (b)
Figure 16: (a) Conical spring undergoes a cyclic sinusoidal loading in a displacement driven context with the displacement uz (t). The corre-
sponding (b) force-displacement curves attest to the fact that the higher the loading rate (small period) is, the stiffer the mechanical response
becomes.
23
aspect of any finite element analysis, we also report in Fig. 17 the Cauchy stress distribution over the conical spring
with four snapshots taken at the instants tA , tB , tC , and tD on the cyclic loading path shown in Fig. 16(a). Aside
from the top portion of the steel which is, in fact, directly subjected to the loading, the absolute maximum stresses are
obviously concentrated in the NBR rubber components. Moreover, the larger areas of viscous overstresses are clearly
demonstrated in the rubber parts which are related to T = 1 [s].
(a) tA tB tC tD
(b) tA tB tC tD
(c) tA tB tC tD
−2.5 σzz [MPa] 2.5

Figure 17: Contour plots of the Cauchy stress σzz over the solid domain at the instants tA , tB , tC , and tD on the cyclic loading path shown in
Fig. 16(a) for (a) T = 1; (b) T = 10; (c) T = 100 [s].
7. Conclusion
In this article, we proposed a new constitutive framework that accounts for the finite viscoelastic mechanical
response of rubber-like materials and touched upon the theoretical, algorithmic, experimental and finally the numerical
aspects thereof. For the ground state elastic response, the eight-chain model, originally proposed by Arruda and Boyce
[1], was extended via the addition of a tube constraint part that is macroscopically linked to the the second invariant
I2 . To this end, the relation between the tube area contraction of the representative chain and the second invariant
is established. Special emphasis was given to the continuum formulation in which the multiplicative split of the
deformation gradient into volumetric and isochoric part was completely dispensed with and all the consequences
thereof drop out of the theoretical formulation. Consequently, compressible counterparts of the eight-chain and the
extended eight-chain models are derived. Such an ansatz recovers the quasi-incompressible behaviour or rubberlike
24
materials as the large κ/µ ratios. Moreover, it allows modeling of slightly compressible filled rubber and description
of relevant degradation mechanisms prior to failure. The proposed extended eight-chain model shows an excellent
fitting performance to Kawabata and Treloar data with only 3-material parameters. The computational efficiency
is considerably better than the well known accurate models in the literature. Algorithmic treatment features the
spectral decomposition of the respective terms and is demonstrated within the context of the finite element method.
Having based our modeling endeavors on micro- and macro-physical grounds, we enriched the study via experimental
observations whose results were directly compared to the model behavior thereby elucidating the performance and the
validity of the model proposed. In fact, the results are in close agreement with the creep and relaxation experiments
carried out at various amplitudes. Later, analyses pertaining to non-homogeneous tests were coducted by means of
a custom designed conical spring undergoing displacement driven cyclic extensions. The finite element simulations
yielded results that are in excellent conformity with the propensities of viscous materials. The proposed numerical
implementation is very robust and quadratically convergent.
Acknowledgment: The authors gratefully acknowledge the financial support from the Scientific Research Projects
Coordination Unit of Middle East Technical University under grant number GAP-302-2020-10203 and T ÜBITAK
Career Award 3501 under grant 315M140.
References
[1] Arruda, E. M., & Boyce, M. C. (1993). A three–dimensional constitutive model for the large stretch behavior of rubber elastic materials.
Journal of the Mechanics and Physics of Solids, 41, 389–412.
[2] Attard, M. M., & Hunt, G. W. (2004). Hyperelastic constitutive modeling under finite strain. International Journal of Solids and Structures,
41, 5327–5350.
[3] Axel Products Physical Testing Services (2019). axelproducts.com. Accessed: 2019-09-28.
[4] Ball, R. C., Doi, M., Edwards, S. F., & Warner, M. (1981). Elasticity of entangled networks. Polymer, 22, 1010–1018.
[5] Bechir, H., Chevalier, L., Chaouche, M., & Boufala, K. (2006). Hyperelastic constitutive model for rubber-like materials based on the first
seth strain measures invariant. European Journal of Mechanics-A/Solids, 25, 110–124.
[6] Bergström, J. S., & Boyce, M. C. (1998). Constitutive modeling of the large strain time–dependent behavior of elastomers. Journal of the
Mechanics and Physics of Solids, 46, 931–954.
[7] Biderman, V. (1958). Calculation of rubber parts. Rascheti na prochnost, 40.
[8] Boyce, M. C., & Arruda, E. M. (2000). Constitutive models of rubber elasticity: A review. Rubber Chemistry and Technology, 73, 504–523.
[9] Christensen, R. (1982). Theory of Viscoelasticity. (2nd ed.). Academic Press.
[10] Cohen, A. (1991). A Padé approximant to the inverse langevin function. Rheological Acta, 30, 270–273.
[11] Coleman, B., & Noll, W. (1961). Foundations of linear viscoelasticity. Quarterly of Applied Mathematics, 33, 239–249.
[12] Coleman, B. D. (1964). Thermodynamics of materials with memory. Archive for Rational Mechanics and Analysis, 17, 1–46.
[13] Cuitino, A., & Ortiz, M. (1992). A material-independent method for extending stress update algorithms from small-strain plasticity to finite
plasticity with multiplicative kinematics. Engineering Computations, 9, 437–451.
[14] Dal, H. (2011). Approaches to the modeling of inelasticity and failure of rubberlike materials: theory and numerics. Ph.D. thesis Technische
Universität Dresden.
[15] Dal, H., Açıkgöz, K., & Badienia, Y. (2021). On the performance of hyperelastic constitutive models for rubber-like materials: A review.
Applied Mechanics Reviews, . To be submitted.
[16] Dal, H., Badienia, Y., Açıkgöz, K., & Denli, F. (2019). A comparative study on hyperelastic constitutive models on rubber: State of the art
after 2006. In Constitutive Models for Rubber XI. Proceedings of the 11th European Conference on Constitutive Models for Rubber (ECCMR
XI), June 25-27, 2019, Nantes, France.
[17] Dal, H., & Kaliske, M. (2009). Bergström–Boyce model for nonlinear finite rubber viscoelasticity: Theoretical aspects and algorithmic
treatment for FE method. Computational Mechanics, 44, 809–823.
[18] Dal, H., Zopf, C., & Kaliske, M. (2018). Micro-sphere based viscoplastic constitutive model for uncured green rubber. International Journal
of Solids and Structures, 132-133, 201–217.
[19] Davidson, J. D., & Goulbourne, N. (2013). A nonaffine network model for elastomers undergoing finite deformations. Journal of the
Mechanics and Physics of Solids, 61, 1784 – 1797.
[20] Doi, M., & Edwards, S. F. (1986). The Theory of Polymer Dynamics. (1st ed.). Clarendon Press, Oxford.
[21] Eidel, B., Stumpf, F., & Schröder, J. (2013). Finite strain viscoelasticity: how to consistently couple discretizations in time and space on
quadrature-point level for full order p ≥ 2 and a considerable speed-up. Computational Mechanics, 52, 463–483.
[22] Erman, B., & Flory, P. J. (1978). Theory of elasticity of polymer networks. ii. the effect of geometric constraints on junctions. The Journal
of Chemical Physics, 68, 5363–5369.
[23] Ferry, J. D. (1980). Viscoelastic Properties of Polymers. (3rd ed.). John Wiley & Sons, Inc., New York.
[24] Flory, P. J., & Erman, B. (1982). Theory of elasticity of polymer networks. 3. Macromolecules, 15, 800–806.
[25] Flory, P. J., & Rehner, J. J. (1943). Statistical mechanics of cross–linked polymer networks: I. rubberlike elasticity. The Journal of Chemical
Physics, 11, 512–520.
[26] de Gennes, P. G. (1971). Reptation of a polymer chain in the presence of fixed obstacles. The Journal of Chemical Physics, 55, 572–579.
[27] Gent, A. N. (1996). A new constitutive relation for rubber. Rubber Chemistry and Technology, 69, 59–61.
25
[28] Gent, A. N., & Thomas, A. (1958). Forms for the stored (strain) energy function for vulcanized rubber. Journal of Polymer Science, 28,
625–628.
[29] Green, M. S., & Tobolsky, A. V. (1946). A new approach to the theory of relaxing polymeric media. The Journal of Chemical Physics, 14,
80–92.
[30] Hart-Smith, L. J. (1966). Elasticity parameters for finite deformations of rubber-like materials. Zeitschrift für angewandte Mathematik und
Physik ZAMP, 17, 608–626.
[31] Haupt, P., & Lion, A. (2002). On finite linear viscoelasticity of incompressible isotropic materials. Acta Mechanica, 159, 87–124.
[32] Heinrich, G., & Kaliske, M. (1997). Theoretical and numerical formulation of a molecular based constitutive tube–model of rubber elasticity.
Computational and Theoretical Polymer Science, 7, 227–241.
[33] Holzapfel, G. A., & Simo, J. C. (1996). A new viscoelastic constitutive model for continuous media at finite thermomechanical changes.
International Journal of Solids and Structures, 33, 3019–3034.
[34] Hong, W. (2011). Modeling viscoelastic dielectrics. Journal of the Mechanics and Physics of Solids, 59, 637–650.
[35] ISO-37:2017(E) (2017). Rubber, vulcanized or thermoplastic – Determination of tensile stress-strain properties. Standard International
Organization for Standardization Geneva, CH.
[36] Kadapa, C., & Hossain, M. (2020). A linearized consistent mixed displacement-pressure formulation for hyperelasticity. Mechanics of
Advanced Materials and Structures, (pp. 1–18). doi:https://doi.org/10.1080/15376494.2020.1762952.
[37] Kaliske, M., & Heinrich, G. (1998). An extended tube–model for rubber elasticity: Statistical–mechanical theory and finite element imple-
mentation. Rubber Chemistry and Technology, 72, 602–632.
[38] Kaliske, M., & Rothert, H. (1997). Formulation and implementation of three–dimensional viscoelasticity at small and finite strains. Compu-
tational Mechanics, 19, 228–239.
[39] Kawabata, S., Matsuda, M., Tei, K., & Kawai, H. (1981). Experimental survey of the strain energy density function of isoprene rubber
vulcanizate. Macromolecules, 14, 154–162.
[40] Khajehsaeid, H., Arghavani, J., & Naghdabadi, R. (2013). A hyperelastic constitutive model for rubber-like materials. European Journal of
Mechanics-A/Solids, 38, 144–151.
[41] Khiêm, V. N., & Itskov, M. (2016). Analytical network-averaging of the tube model:: Rubber elasticity. Journal of the Mechanics and Physics
of Solids, 95, 254–269.
[42] Kilian, H.-G., Enderle, H. F., & Unseld, K. (1986). The use of the van der waals model to elucidate universal aspects of structure-property
relationships in simply extended dry and swollen rubbers. Colloid and Polymer Science, 264, 866–876.
[43] Koprowski-Theiss, N., Johlitz, M., & Diebels, S. (2011). Characterizing the time dependence of filled epdm. Rubber Chemistry and
Technology, 84, 147–165.
[44] Korba, A. G., & Barkey, M. E. (2017). New model for hyper-elastic materials behavior with an application on natural rubber. In ASME 2017
12th International Manufacturing Science and Engineering Conference collocated with the JSME/ASME 2017 6th International Conference
on Materials and Processing (pp. V002T03A020–V002T03A020). American Society of Mechanical Engineers.
[45] Kröner, E. (1960). Allgemeine Kontinuumstheorie der Versetzungen und Eigenspannungen. Archive for Rational Mechanics and Analysis,
4, 273–334.
[46] Kuhn, W. (1934). über die gestalt fadenförmiger moleküle in lösungen. Kolloid-Zeitschrift, 68, 2–15.
[47] Kuhn, W. (1936). Beziehungen zwischen molekülgröße, statistischer molekülgestalt und elastischen eigenschaften hochpolymerer stoffe.
Kolloid-Zeitschrift, 76, 258–271.
[48] Kuhn, W., & Grün, F. (1942). Beziehungen zwischen elastischen konstanten und dehnungsdoppelbrechung hochelastischer stoffe. Kolloid-
Zeitschrift, 101, 248–271.
[49] Lambert-Diani, J., & Rey, C. (1999). New phenomenological behavior laws for rubbers and thermoplastic elastomers. European Journal of
Mechanics-A/Solids, 18, 1027–1043.
[50] Lee, E. H. (1969). Elastic-plastic deformation at finite strains. ASME. Journal for Applied Mechanics, 36, 1–6.
[51] Linder, C., Tkachuk, M., & Miehe, C. (2011). A micromechanically motivated diffusion-based transient network model and its incorporation
into finite rubber viscoelasticity. Journal of the Mechanics and Physics of Solids, 59, 2134–2156.
[52] Lion, A. (1996). A constitutive model for carbon black filled rubber. Experimental investigations and mathematical representations. Contin-
uum Mechanics and Thermodynamics, 8, 153–169.
[53] Lion, A. (1997). On the large deformation behaviour of reinforced rubber at different temperatures. Journal of the Mechanics and Physics of
Solids, 45, 1805–1834.
[54] Lodge, A. S. (1956). A network theory of flow birefringence and stress in concentrated polymer solutions. Transactions of the Faraday
Society, 52, 120–130.
[55] Lopez-Pamies, O. (2010). A new i1-based hyperelastic model for rubber elastic materials. Comptes Rendus Mecanique, 338, 3–11.
[56] Lubliner, J. (1985). A model of rubber viscoelasticity. Mechanics Research Communications, 12, 93–99.
[57] Marckmann, G., & Verron, E. (2006). Comparison of hyperelastic models for rubber-like materials. Rubber Chemistry and Technology, 12,
835–858.
[58] Marsden, J. E., & Hughes, T. J. R. (1983). Mathematical Foundations of Elasticity. Prentice-Hall, Englewood Cliffs, New Jersey.
[59] Miehe, C., & Göktepe, S. (2005). A micro–macro approach to rubber–like materials. Part II: The micro–sphere model of finite rubber
viscoelasticity. Journal of the Mechanics and Physics of Solids, 53, 2231–2258.
[60] Miehe, C., Göktepe, S., & Lulei, F. (2004). A micro–macro approach to rubber–like materials. Part I: The non–affine micro–sphere model of
rubber elasticity. Journal of the Mechanics and Physics of Solids, 52, 2617–2660.
[61] Miehe, C., & Keck, J. (2000). Superimposed finite elastic–viscoelastic–plastoelastic stress response with damage in filled rubbery polymers.
Experiments, modelling and algorithmic implementation. Journal of the Mechanics and Physics of Solids, 48, 323–365.
[62] Mooney, M. (1940). A theory of large elastic deformation. Journal of Applied physics, 11, 582–592.
[63] Ogden, R. (1972). Large deformation isotropic elasticity: on the correlation of theory and experiment for incompressible rubberlike solids.
Proceedings of the Royal Society of London A, 326, 565–584.
26
[64] Phan-Thien, N. (1978). A nonlinear network viscoelastic model. Journal of Rheology, 22, 259–283.
[65] Phan-Thien, P. (2002). Understanding Viscoelasticity. (1st ed.). Springer-Verlag Berlin-Heidelberg.
[66] Pucci, E., & Saccomandi, G. (2002). A note on the gent model for rubber-like materials. Rubber chemistry and technology, 75, 839–852.
[67] Raghunath, R., Juhre, D., & Klüppel, M. (2016). A physically motivated model for filled elastomers including strain rate and amplitude
dependency in finite viscoelasticity. International Journal of Plasticity, 78, 223–241.
[68] Reese, S., & Govindjee, S. (1998). A theory of finite viscoelasticity and numerical aspects. International Journal of Solids and Structures,
35, 3455–3482.
[69] Rivlin, R. S., & Saunders, D. W. (1951). Large elastic deformations of isotropic materials vii. experiments on the deformation of rubber.
Proceedings of the Royal Society of London A, 243, 251–288.
[70] Saxena, P., Hossain, M., & Steinmann, P. (2013). A theory of finite deformation magneto-viscoelasticity. International Journal of Solids and
Structures, 50, 3886 – 3897.
[71] Shariff, M. H. B. M. (2000). Strain energy function for filled and unfilled rubberlike material. Rubber chemistry and technology, 73, 1–18.
[72] Simo, J. C. (1987). On a fully three–dimensional finite–strain viscoelastic damage model: Formulation and computational aspects. Computer
Methods in Applied Mechanics and Engineering, 60, 153–173.
[73] Simo, J. C. (1992). Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the
infinitesimal theory. Computer Methods in Applied Mechanics and Engineering, 99, 61–112.
[74] Simo, J. C., & Miehe, C. (1992). Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation.
ASME Journal of Applied Mechanics, 98, 41–104.
[75] Steinmann, P., Hossain, M., & Possart, G. (2012). Hyperelastic models for rubber-like materials: consistent tangent operators and suitability
for treloar’s data. Archive of Applied Mechanics, 82, 1183–1217.
[76] Swanson, S. R. (1985). A constitutive model for high elongation elastic materials. Journal of engineering materials and technology, 107,
110–114.
[77] Tanaka, F., & Edwards, S. F. (1992). Viscoelastic properties of physically crosslinked networks. I: Non–linear stationary viscoelasticity.
Journal of Non–Newtonian Fluid Mechanics, 43, 247–271.
[78] Tanaka, F., & Edwards, S. F. (1992). Viscoelastic properties of physically crosslinked networks. II: Dynamic mechanical moduli. Journal of
Non–Newtonian Fluid Mechanics, 43, 289–309.
[79] Taylor, R. L. (2017). FEAP – A Finite Element Analysis Program, Version 8.5. University of California at Berkeley, Berkeley, California.
[80] Treloar, L. R. G. (1943). The elasticity of a network of long-chain molecules-ii. Transactions of the Faraday Society, 39, 241–246.
[81] Treloar, L. R. G. (1944). Stress–strain data for vulcanised rubber under various types of deformation. Transactions of the Faraday Society,
40, 59–70.
[82] Treloar, L. R. G. (1975). The Physics of Rubber Elasticity. (3rd ed.). Clarendon Press, Oxford.
[83] Treloar, L. R. G., & Riding, G. (1979). A non–gaussian theory of rubber in biaxial strain. i.mechanical properties. Proceedings of the Royal
Society of London A, 369, 261–280.
[84] Truesdell, C., & Noll, W. (1965). The non–linear field theories of mechanics. In S. Flügge (Ed.), Encyclopedia of Physics. Springer–Verlag,
Berlin volume III/3.
[85] Valanis, K. C., & Landel, R. F. (1967). The strain-energy function of a hyperelastic material in terms of the extension ratios. Journal of
Applied Physics, 38, 2997–3002.
[86] Wall, F. T. (1942). Statistical thermodynamics of rubber. The Journal of Chemical Physics, 10, 132–.
[87] Wang, M. C., & Guth, E. (1952). Statistical theory of networks of non–gaussian flexible chains. The Journal of Chemical Physics, 20,
1144–1157.
[88] Weber, G., & Anand, L. (1990). Finite deformation constitutive equations and a time integration procedure for isotropic hyperelastic-
viscoplastic solids. Computer Methods in Applied Mechanics and Engineering, 79, 173–202.
[89] Yeoh, O. H. (1990). Characterization of elastic properties of carbon-black-filled rubber vulcanizates. Rubber chemistry and technology, 63,
792–805.
[90] Yeoh, O. H., & Fleming, P. D. (1997). A new attempt to reconcile the statistical and phenomenological theories of rubber elasticity. Journal
of Polymer Science Part B: Polymer Physics, 35, 1919–1931.
27
View publication stats

Daletal 20 JMPS

Uploaded by

Copyright:

Available Formats

Daletal 20 JMPS

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Daletal 20 JMPS

Uploaded by

Copyright:

Available Formats

See discussions, stats, and author profiles for this publication at: https://www.researchgate.

An extended eight-chain model for hyperelastic and ﬁnite viscoelastic

Article in Journal of the Mechanics and Physics of Solids · September 2020

Hüsnü Dal Osman Gueltekin

SEE PROFILE SEE PROFILE

The user has requested enhancement of the downloaded file.

Hüsnü Dala,∗, Osman Gültekina , Kemal Açıkgöza

Preprint submitted to Journal of LATEX Templates September 17, 2020

2. Fundamental maps for a deformable solid and kinematics

F −T : TX∗ B0 → Tx∗ B; n = F −T N (2)

dx = F dX , da = cof[F ]dA , dv = det[F ]dV , (3)

le = F˙ e F e−1 and lv = F˙ v F v−1 . (13)

l =d+w where d = sym gl and w = skew gl , (15)

ϕ ∗ (•)♭ = F −T (•)♭ F −1 and ϕ ∗ (•)♯ = F (•)♯ F T . (17)

ϕ ∗ (•)♭ = F T (•)♭ F and ϕ ∗ (•)♯ = F −1 (•)♯ F −T . (18)

TX∗ B0 F v−T T̄X̄∗ B F e−T TX∗ B0 F v−T T̄X̄∗ B F e−T

3. A compressible finite viscoelastic constitutive model for rubberlike materials

3.1. Rheology of the constitutive model

Ψ(g, F , F e ) = Ψe (g, F ) + Ψne (g, F e ) = Ψe (C) + Ψne (C e ) , (20)

τ e := 2∂g Ψe (g, F ) and S e := 2∂C Ψe (C) . (24)

3.2.1. Micro-molecular motivation for the extended eight-chain model

ψ(λ, ν) = −θη(λ, ν) with η = kB ln p(λ, ν) ❀ η(λ, ν) = ηf (λ) + ηc (ν) , (25)

ψ(λ, ν) = ψf (λ) + ψc (ν) . (26)

Ψe (g, F ) := Ψef (g, F ) + Ψec (g, F ) (31)

the free-chain and tube constraint parts read

3.2.2. Compressible extended eight-chain model

3.3. Non-equilibrium part of the constitutive model

3.5. Relaxation kinetics of a single chain

(a) (b) (c)

Substitution of the relations stated below

into (54), we readily obtain the rate equation

3.6. Evolution equation

4. Algorithmic setting for the constitutive model

∆τ = Calgo : 12 £∆ϕg + ∇x (∆ϕ) τ + τ ∇Tx (∆ϕ) (62)

4.1. Stress and moduli terms for the elastic part

Ce := 4∂gg2 Ψe(J, I1 , I2 ) (64)

Therein, the terms in regard to the Jacobian J reads

4.2. Stresses and moduli terms for the viscous part

4.2.1. Integration of the evolution equation:

Elastic predictor (trial step) E : C v−1

4.2.2. Algorithmic moduli for the non-equilibrium part:

where we have introduced GabIJKL := MIK

Ingredients (phr) Composition

5.2. Creep and relaxation experiments on NBR rubber

P11 = 0.10 P11 = 0.10

6. Representative Numerical Analyses

Extended eight-chain model Eight-chain model

6.2. Assessment of the hyperelastic behavior of the extended eight-chain model

0.2 λ1 = 1.04 0.7 λ1 = 1.3

0.2 λ1 = 1.04 0.7 λ1 = 1.3

0.2 λ1 = 1.04 0.7 λ1 = 1.3

0.2 λ1 = 1.04 0.7 λ1 = 1.3

Extended eight-chain model parameters

Table 4: Viscoelastic material parameters in regard to Maxwell branches used.

The eight-chain model (viscous part) & evolution equation

1 model fit 1 model fit

2 model fit 2 model fit

P11 = 1.00 P11 = 1.00

6.4. Non-homogeneous cyclic loading of a conical spring