CONSTITUTIVE MODELING OF CROSS-LINKED POLYMERS DURING AGING

By

Seyedamir Bahrololoumi Tabatabei

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

Civil Engineering - Doctor of Philosophy
Mechanical Engineering - Dual Major

2024

ABSTRACT

Due to unique features such as high stretchability, chemical toleration, wear resistance, anti-

slipping, and durability; industrial applications of cross-linked polymers have been growing through

the past decades. Their practical function can cover a broad range of applications from insensitive

to sensitive industries such as automotive, medical, mining, and petroleum. In most such cases, the

polymeric components are simultaneously exposed to mechanical loads as well as environmental

elements like humidity and oxygen. This affair can considerably reduce the nominal service life

of the polymers and cause irreversible damage to their functionality. In the literature, this chem-

ical/physical phenomenon is generally regarded as the aging mechanism, which is classified into

different subcategories based on the source of aging. Thermo-oxidation, hydrolysis, and hygrother-

mal are some of the main types of chemical aging in which oxygen, temperature, and water play

major roles in the deterioration of polymer matrix.

Such classes can be either the result of individual aging elements or due to combined synergized

effects of the above-mentioned environmental factors. Any type of aging can induce irreversible

damage to the polymer matrix, which results in the reduction of materials’ nominal service life.

While predicting any of these types of aging is of great interest to researchers in the last few years,

there is still a long way to find a comprehensive model which can consider the synergized effects

of chemical aging.

Here, multi-physics models are devised to describe the solo and combined long-term effects

of different environmental elements on the constitutive behavior of cross-linked polymers. These

models are relevant for all decay mechanisms formed by the occurrence of two simultaneous micro-

mechanisms; (i) formation/reduction of the cross-links, and (ii) chain scission. We modeled each

sub-mechanism through a linear kinetic equation which is then coupled into the modular network

concept to allow the calculation of synergies and parallel damage accumulation through multiple

sources. In this respect, hydrolytic aging is considered to be induced by chain scission and reduc-

tion of the cross-link, while thermal aging is considered a superposition of two competing effects of

cross-link formation and chain scission. Besides, hygrothermal damage can be regarded as a race

among three micro-structural phenomena; i) chain scission due to the presence of temperature, ii)

reduction of cross-links attributed to the attendance of water, and iii) formation of cross-links as a

consequence of oxygen reactions.

These proposed models are micro-mechanically based and are mainly relevant to thin samples

due to our underlying assumption of homogeneous diffusion of oxygen and water throughout the

matrix. These models have been validated against extensive data sets obtained from experiments

we specifically designed for concept validation. In view of their interpretation, precision, and deep

insight, it provides insight into the nature of damage accumulation. The model is a good choice

for advanced implementation in FE applications.

Copyright by
SEYEDAMIR BAHROLOLOUMI TABATABAEI
2024

"In respectful dedication to the cherished memory of my esteemed family members: my beloved
wife, Gelareh, whose unwavering love has consistently fortified me; my father, Majid, a guiding
light and and my mother, Saeedeh, whose unwavering support has shaped my path. It remains
deeply regrettable that my father, who selflessly invested in my future, couldn’t witness its
fulfillment due to the tragic impact of COVID-19."

v

ACKNOWLEDGEMENTS

As I began penning this acknowledgment, I encountered the immense challenge of expressing ade-

quate gratitude for the myriad forms of assistance, support, and understanding offered by numerous

individuals throughout the journey of crafting this dissertation. The fear of inadvertently overlook-

ing or inadequately representing someone’s invaluable contributions loomed large. It is undeniable

that an extensive community has played an integral role in guiding me toward completing this

work, and the individuals mentioned herein are just a fraction of this supportive network.

First and foremost, I extend my heartfelt gratitude to my esteemed advisers, Dr. Roozbeh

Dargazany and Dr. Thomas Pence. Your guidance and belief in my capabilities have been instru-

mental in nurturing my growth as an independent researcher.

I wish to express profound appreciation to my esteemed committee members, Dr. Nizar Lajnef,

Dr. Sara Roccabianca, and Prof. Rafael Auras. Their wealth of expertise, combined with the priv-

ilege of learning from their classes, significantly contributed to shaping my intellectual evolution

throughout my graduate studies.

I am sincerely thankful to my invaluable colleagues and friends within the High-Performance

Materials (HPM) group, notably Vahid Morovati, Mamoon Shaafaey, Hamid Mohammadi, Aref

Ghaderi, Yang Chen, Sharif Alazhary, and Ramin Akbari. Their collaborative spirit and insightful

discussions have been invaluable contributions to this endeavor.

I am profoundly indebted to Dr. Manoochehr Najmi for his unwavering guidance and invalu-

able advice, which have been pivotal throughout my journey. Additionally, I extend heartfelt grati-

tude to my mother-in-law, Fariba. While the list of friends is extensive, I am immensely grateful to

each of you for transforming my journey into a joyous and supportive experience, akin to a second

family away from home. I extend my sincerest wishes for success in all your future endeavors.

The successful completion of this dissertation owes immeasurably to my beloved wife, Gelareh

Najmi, my best friend. Her unwavering support, enduring patience, and selfless sacrifices during

my trials, absences, moments of frustration, and impatience are beyond commendation.

vi

CHAPTER 1 Introduction .

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

TABLE OF CONTENTS

CHAPTER 2 Fundamentals of Continuum Mechanics and Non-linear Finite Element
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Method .

.

.

.

.

.

.

.

1

3

CHAPTER 3 A Multi physics Constitutive Model to Predict Hydrolytic Aging in Quasi

static Behaviour of Thin Cross linked Polymers . . . . . . . . . . . . . . . 21

CHAPTER 4 A Physically Based Model for Thermo-Oxidative and Hydrolytic Aging

of Elastomers .

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

CHAPTER 5 A Multi Physics Approach on Modeling of Hygrothermal Aging and its

Effects on Constitutive Behavior of Cross linked Polymers . . . . . . . . . 79

CHAPTER 6 Thermal Aging Coupled with Cyclic Fatigue in Cross-linked Polymers:

Constitutive Modeling & FE Implementation . . . . . . . . . . . . . . . . 113

CHAPTER 7 Summary and Future Works . . . . . . . . . . . . . . . . . . . . . . . . . 139

BIBLIOGRAPHY .

.

.

.

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

vii

CHAPTER 1

Introduction

1.1 Motivation

Cross-linked polymers have become ubiquitous throughout the world due to their excellent prop-

erties such as high extensibility, strength, abrasion and chemical resistance, and durability [1, 2].

They are widely used in different industries such as automotive, transportation, power distribution,

and aerospace. In view of their demanding applications, rubber-like components should be reli-

ably designed to be subjected to a high number of cyclic loading in harsh environments. However,

their reliability considerably decreases over time due to a process known as aging. This process

usually causes a significant change in the properties of rubber-like materials over their service life

[3, 4]. To this end, reliable lifetime estimation of products made by cross-linked polymers can be

extremely challenging, especially for products that are exposed to multiple damages.

The aging process can happen due to an individual or combined effects of different environ-

mental factors such as humidity, temperature, and radiation [5, 6]. These environmental factors

facilitate degradation mechanisms. These mechanisms can be categorized into those with major

impacts such as hygro-thermal [7, 8], hydrolysis [9, 10], thermo-oxidation [11, 6, 12], and those

with the less severe impact such as chemical corrosion [13], and ozone cracking [14]. While these

mechanisms potentially can significantly reduce the service life of an elastomeric component, the

latter group is not that prevalent in nature. Although predicting the individual contribution of dif-

ferent mechanisms has been the center of attention in the literature, there is still a long way to go

in modeling aging mechanisms.

Usually, polymers are used in applications where in situ testing is expensive due to inaccessi-

bility. Therefore, predictive models are needed to prevent the catastrophic failure of the structures

1

because of aging. Subsequently, a computationally simple yet precise model that can predict the

mechanical properties of cross-linked polymers during aging is in high demand. Albeit their popu-

larity, theoretical efforts for both thermo-oxidative aging and hydrolysis are limited. These efforts

can be categorized into two different approaches, phenomenological and micro-mechanical. Phe-

nomenological approaches use a mathematical model, usually without a physical meaning, to de-

scribe the changes in material properties based on a thermodynamic framework [15, 16]. Though

phenomenological models can accurately model the behavior, they do not offer any insight into

the relationship between model parameters and real material ones. On the other hand, micro-

mechanical models are based on physical motivation. They developed from statistical mechanics

of polymer structure and physics of polymer chain network. Consequently, their parameters do

have physical meaning [17]. The lack of a proper micro-mechanical model for modeling elastomer

constitutive behavior prompted us to begin this research. The main goal of this research is to de-

velop a micro-mechanical model that can accurately explain the elastic and inelastic responses of

elastomers in various aging environments.

2

CHAPTER 2

Fundamentals of Continuum Mechanics and Non-linear Finite Element

Method

The constitutive modeling of elastomers involves nonlinear continuum mechanical quantities, which

describe the behavior of materials undergoing a large deformation. Thus, a basic foundation on

continuum mechanics, including fundamental geometric mappings and basic stress measures of a

solid body undergoing large deformations, is briefly discussed in this chapter. For a comprehensive

description and understanding of the topic, the reader can refer to continuum mechanics reference

materials such as [2].

2.1 Some Notes on Continuum Mechanics

In continuum mechanics theory, it is assumed that an object fully occupies the space with its sub-

stances. Physical properties of a solid or fluid medium are related through mathematical tensor

measures. These measures are independent of their observed coordinate system in general. How-

ever, they can be represented in different coordinate systems. In this section, the relation of various

mechanical measurements is briefly reviewed.

2.1.1 Deformation gradient

A continuum body Bt0 in a 3D Euclidean space at time t0, in which any arbitrary point P0 with

respect to an arbitrary basis can be represented by X ∈ E3. As the body deforms in the space, the

medium occupies its current configuration at Bt. Given that the geometrical mapping of the body

from Bt0 to Bt is one-to-one, any point from the initial configuration, P0, uniquely maps to its new

configuration coordinate, P , with an arbitrary basis of x ∈ E3.

3

The geometrical mapping of the regions of the body from current configuration Bt0 to the

reference configurations Bt in any arbitrary coordinate system, ei (i = 1, 2, 3), can represent the

position of point P and P0 by

x = ˆx (θ1, θ2, θ3, t), X = ˆX (θ1, θ2, θ3) = ˆx (θ1, θ2, θ3, t0).

i = 1, 2, 3 ,

(2.1)

Therefore, the displacement vector, u, of point P can be written as

u = ˆu(θ1, θ2, θ3, t) = x − X

(2.2)

Given a Euclidean space with a set of orthonormal (Cartesian) basis vectors, say e i

(i =

1, 2, 3), one can express each point of the body as

X = X ie i, X j = X · e j,

u = uie i,

uj = u · e j,

j = 1, 2, 3,

j = 1, 2, 3,

(2.3)

x = xie i,

xj = x · e j = X j + uj,

j = 1, 2, 3,

where the Einstein notation, summation over repeated indices, is applied. Deformation can be

related to the tangent vectors of the coordinate lines in each of the configurations. Tangent vectors

of sufficiently differentiable X and x can be written as

G i =

∂X
∂θi =,

g i =

∂x
∂θi ,

i = 1, 2, 3 ,

(2.4)

The relative motion of an arbitrary point with respect to its adjacent point can be calculated in a

direction a through the directional derivative as

d
ds

(cid:12)
(cid:12)
ˆx (cid:0)X + sa, t(cid:1)
(cid:12)
(cid:12)s=0

= lim
∆s→0

ˆx (cid:0)X + ∆sa, t(cid:1) − ˆx (cid:0)X , t(cid:1)
∆s

= (cid:0)Grad x (cid:1)a.

(2.5)

Therefore, second-order deformation gradient tensor, (cid:0)Grad x (cid:1), can be written by a linear mapping

4

of vector a into vector d

ds ˆx (cid:0)X + sa, t(cid:1)

(cid:12)
(cid:12)
(cid:12)
(cid:12)s=0

as

F = Grad x =

dˆx
dX

.

(2.6)

where, dX is a infinitesimal element before deformation and dˆx is the same element on the body

after deformation. This tensor plays a pivotal role in deformation kinematics, which can describe

the relative motion of material elements during deformation. Therefore, one has the relation of

infinitesimal elements in reference and current configurations as

dx = FdX ,

dX = F−1dx .

(2.7)

Moreover, the change of volume and surface elements can be related to the deformation gradient

tensor through these relations. To this end, a volume in the reference configuration of three non-

co-planar vectors, dX 1, dX 2 and dX 3, can be defined as

dV0 = [dX 1dX 2dX 3] = (cid:0)dX 1 × dX 2

(cid:1) · dX 3.

(2.8)

Using Eq. 2.7, each of these vectors deforms in the current configuration to

dx 1 = FdX 1,

dx 2 = FdX 2,

dx 3 = FdX 3.

(2.9)

Thus, the volume of the element in the current configuration can be written as

dV = (cid:2)dx 1dx 2dx 3

(cid:3) = (cid:0)dx 1 × dx 2

(cid:1) · dx 3 = JdV0,

where

J =

dV
dV0

= (cid:12)

(cid:12)Fi
.j

(cid:12)
(cid:12) = detF > 0.

(2.10)

In addition, by defining the surface element of reference and current state as dA0 = dX 1 × dX 2

5

and dA = dx 1 × dx 2, the relation of the surface areas can be calculated by substituting Eq. 2.9

into dV0 = JdV as

dA = JF−T dA0.

(2.11)

where dA = |dA| and dA0 = |dA0| are the surface areas in the current and reference configuration,

respectively. An element length in reference and current states can be calculated similarly as

∥dx ∥2 = dx · dx = dX (cid:0)FTF(cid:1)dX = dX CdX ,

∥dX ∥2 = dX · dX = dx (cid:0)F−TF−1(cid:1)dx = dx b−1dx ,

(2.12)

where C = FTF and b = FFTare the right and left Cauchy-Green tensors, respectively. The

change in the length of a linear element, the stretch of a material element, is defined as the ratio of

the deformed to the referenced length of the material element. Deformation changes the length of

an element dX in direction N in initial state to dx in direction n in the current state as

λ(cid:0)N (cid:1) =

dx
dX

=

(cid:115)

(cid:114)

∥dx ∥2
∥dX ∥2 =

dXN CN dX
dX 2

= (cid:0)N CN (cid:1) 1
2 ,

and

λ(cid:0)n(cid:1) = (cid:0)nb−1n(cid:1)− 1
2 .

Another measure of the change in element length during deformation can be written as

(2.13)

(2.14)

∥dx ∥2 − ∥dX ∥2 = 2dX EdX

= 2dx edx

(2.15)

(cid:0)C − I(cid:1) = 1

2

(cid:0)FTF − I(cid:1) called the Green-Lagrange strain tensor and e = 1

2

(cid:0)I −

where E = 1
2
b−1(cid:1) = 1

2

(cid:0)I − F−TF−1(cid:1) is Almansi strain tensor.

6

2.1.2 Deformation rate

The material velocity gradient can be defined as similar to the deformation gradient by

L = Grad ˙x =

(cid:34)

∂x (X , t)
∂t

(cid:33)

(cid:35)

=

∂
∂t

(cid:32)

∂x
∂X

∂
∂X

= ˙F.

(2.16)

The spatial velocity gradient, the derivative of a spatial velocity field v with respect to the current

configuration, can be written as

l = grad ˙x =

∂v
∂X

∂X
∂x

= ˙FF−1.

The spatial gradient of velocity can be decomposed to a symmetric part d = 1
2

(2.17)

(cid:16)

l + lT(cid:17)

and a

skew-symmetric part w = −wT = 1

2(l − lT). The symmetric part is called the deformation rate

d =

1
2

l + lT(cid:17)
(cid:16)

=

1
2

(cid:16) ˙FF−1 + F−T ˙FT(cid:17)

=

F−T ˙CF−1.

1
2

In addition, the skew+symmetric part is referred to spin (vorticity) tensor

w = −wT =

1
2

(l − lT) =

1
2

(cid:16) ˙FF−1 − F−T ˙FT(cid:17)

.

The rate of Green-Lagrange strain tensor ˙E can be written as

˙E =

˙C =

1
2

1
2

FTdF.

(2.18)

(2.19)

(2.20)

The rate of the volume change can be calculated by the time derivative of the determinant of

the deformation gradient as

˙J =

∂detF
∂t

= J trd

(2.21)

7

2.1.3 Stress measures

The stress is defined by traction force vector dFs per unit of area, which is defined as the neigh-

boring continuum points of the body dA. The stress in each arbitrary point P0 and P in reference

and its counterpart in the current configuration is defined based on infinitesimal areas dA and da

around them with the unit vectors of N and n normal to them, respectively. Thus, the relation of

the traction force and stress can be written as

dFs = T dA = tda,

(2.22)

where the vectors T and t are the traction force in reference and deformed configurations, respec-

tively. Therefore, second-order tensors σσσ and PPP can be expressed based on the Cauchy’s stress

theorem as

t = σσσn T = PPP N ,

(2.23)

where σσσ represents the Cauchy stress tensor related to force on the body surface in the current

configuration and PPP is the first Piola-Kirchhoff stress second-order tensor defined based on the

surface area in the initial configuration. The first Piola-Kirchhoff stress and Cauchy stress tensor

can be related through Nanson’s formula 2.11 and 2.23 as

dF = tdA = σσσdA = JσσσF−Tda.

(2.24)

The second Piola-Kirchhoff stress tensors S can be calculated by pulling back the force vector dfs

in the current state to the reference configuration. Similarly, another spatial stress measure, the

Kirchhoff stress tensor τ , can be calculated and converted to other measures as the following set

of equations

P = JσσσF−T = τ F−T = FS, S = JF−1σσσF−T,

τττ = Jσσσ

(2.25)

8

Note that nominal/first Piola-Kirchhoff stress is not a symmetric tensor due to its two-point coor-

dinate systems.

2.1.4 Balance principles

For a body B with volume V , mass M , and boundary surface A in the current state, the rate of

change of the linear momentum is directly proportional to the force applied on the body, which

can be defined by

(cid:90)

M

v dM =

(cid:90)

V

ρ ˙x dV,

where ρ is the density and v = ˙x is the velocity vector of a particle. The momentum can be

calculated by summation of a body force and a surface force as

d
dt

(cid:90)

V

ρ ˙x dV =

(cid:90)

V

f dV +

(cid:90)

A

tdA.

(2.26)

The surface forces can be obtained through the Cauchy theorem (2.23) and the divergence theorem

(cid:90)

A

tdA =

(cid:90)

A

σσσndA =

(cid:90)

V

div σσσdV.

By substituting Eq. 2.27 into Eq. 2.26, the balance equation is then

(cid:90)

V

(cid:0)divσσσ + f − ρ¨x (cid:1)dV = 0.

For any arbitrary point of the body, one can rewrite Eq. 2.28 as

(2.27)

(2.28)

div σσσ + f = ρ¨x .

(2.29)

9

In order to calculate the balance of mechanical energy, one can multiply Eq. 2.29 with the velocity

vector v as

v · div σσσ + v · f = ρv · ¨x .

(2.30)

Eq. 2.30 can be further simplified considering the symmetry of the Cauchy stress tensor as

v · divσσσ = div (vσσσ) − σσσ : grad v = div (vσσσ) − σσσ : d,

(2.31)

which gives

div(vσσσ) − σσσ : d + v · f = ρ

d
dt

(cid:18) 1
2

(cid:19)

v · v

.

(2.32)

Integrating Eq. 2.32 over the volume of the body and considering Eq. 2.27 yield to

d
dt

(cid:90)

M

(cid:18) 1
2

(cid:19)

v · v

dM +

(cid:90)

V

(σσσ : d) dV =

(cid:90)

A

(v · t) dA +

(cid:90)

V

(v · f ) dV.

(2.33)

The balance of mechanical energy can be rewritten by taking into account the physical meaning of

each term in the Eq. 2.33 as

˙K + W = P,

(2.34)

where P,

˙K, and W refer to the power of external forces, the stress power, and the changes in the

kinetic energy, respectively. The power of external forces can be written as

P =

(cid:90)

A

(v · t) dA +

(cid:90)

V

(v · f ) dV.

The kinetic energy and the stress power of the system is then formulated by

K =

W =

(cid:90)

V
(cid:90)

V

(cid:19)

v · v

dV,

(cid:18)

ρ

1
2

(σσσ : d) dV.

10

(2.35)

(2.36)

Note that the stress power can be calculated through other strain-stress measures and each of the

stress measures σσσ, P, and S are conjugate with deformation measures d, ˙F, and ˙E, respectively.

For instance, the stress power can be calculated using the first Piola-Kirchhoff stress and deforma-

tion gradient as

W =

(cid:90)

V0

(cid:16)

(cid:17)

P : ˙F

dV0 =

(cid:90)

V

J −1 (cid:16)

P : ˙F

(cid:17)

dV.

Thus, the stress power can also be obtained from

(cid:16)

(cid:17)

P : ˙F

= (σσσ : d) =

(cid:16)
S : ˙E

(cid:17)

=

(cid:18)

S :

(cid:19)

˙C

.

1
2

In addition, the total energy dissipation due to the deformation D can be obtained

D = W −

(cid:90)

V0

˙Ψ dV0

where ˙Ψ refers to the rate of stored energy per unit volume in the reference state.

(2.37)

(2.38)

(2.39)

2.1.5 Thermo-elasticity

The relations between mechanical and thermal energy is called thermo-elasticity, which can be

derived using the first and second laws of thermodynamics. The 1st law of thermodynamics is the

conservation of energy, which is formulated as

dU = dQ + dW.

(2.40)

where dU , dQ, and dW are the changes in the internal energy of the system, the heat absorbed by

the system, and the work done on the system. The U and W can be expressed as

U =

(cid:90)

V

ρU dV, W =

dW
dt

(2.41)

11

where U is the internal energy density. According to the second law, the changes in the entropy of

a reversible process (elastic deformation), dS at absolute temperature T is formulated as

dS =

dQ
T

.

(2.42)

Moreover, the change in the Helmholtz energy A is defined as

dA = dU − SdT − T dS.

(2.43)

At a constant temperature, the change in Helmholtz energy is equal to the change in the work done

on the system dW = dA . In the small displacement dl of solid-like structures, dW due to force f

can be written by

dW = f dl − phdV,

(2.44)

where ph and dV denote the hydro-static pressure and the volume change. In the case of incom-

pressible material like elastomers, dV can be neglected. So, the force of an isothermal process can

be obtained using Eq. 2.43 and Eq. 2.44 as

f =

dU
dl

− T

dS
dl

.

(2.45)

The first part of Eq. 2.45 represents the energetic interactions of a single molecule or the volume

change of the body [18]. However, experimental evidences suggest that the contribution of ener-

getic force in the total force is negligible in moderate and large deformations of elastomers. Thus,

the force of the system can be approximated only by entropy component as

f = −T

dS
dl

.

12

(2.46)

For an entropic material, the Helmholtz free energy required to perturb the entropy of the system

is

dA = −T dS.

(2.47)

2.2 Thermodynamic Consistency

According to the second law of thermodynamic that is used in continuum mechanics, elastic energy

of the system mainly increases due to the decrease in entropy, which is known as Clausius-Duhem

inequality

D ≥ 0.

(2.48)

The constitutive model is thermodynamically consistent if and only if Clausius-Duhem inequal-

ity is satisfied on all points of body at any time. In other words, the energy balance should be

satisfied during a thermodynamical deformation process. The energy balance can be rewritten by

substituting Eq. 2.40 to Eq. 2.34 as

˙K − P =

d
dt

(U − Q) = ˙U − ˙Q,

where the thermal power ˙Q is defined as

˙Q =

(cid:90)

V

ρQ dV −

(cid:90)

A

ρq · n dA.

(2.49)

(2.50)

In Eq. 2.50 ,Q and q are the transferred heat and the surface heat flux per unit of mass, respectively.

Moreover, Eq. 2.49 can be rewritten as

σσσ : l = ρ ˙U − ρQ + div q .

(2.51)

13

Therefore, Clausius-Duhem inequality can be obtained by

ρD = ρT ˙S − (ρQ − div q ) ≥ 0,

(2.52)

where S is the entropy density. Clausius-Duhem inequality can be further simplified by substitution

of the energy balance law (2.51) and mass conservation (ρ0 = Jρ) as

D = T ˙S −

(cid:18)

˙U −

1
ρ0

(cid:19)

PPP : ˙F

≥ 0.

(2.53)

considering the stored energy density per unit of reference volume Ψ = ˆΨ (F, T, ΩΩΩ) as a function

of set of internal variables ΩΩΩ. In most of the constitutive models, the internal variables have been

adopted to describe the history dependent dissipative effects. Note that all the quantities of the

Clausius-Duhem inequality, such as S, q , and P should be functions of internal variables. For a

more comprehensive review of thermo-elasticity, the reader is referred to [19, 20]. Substituting the

Helmholtz free energy Ψ = ρ0 (U − T S) into Eq. 2.53 results in

ρJD = − ˙Ψ − ρJS ˙T + PPP : ˙F ≥ 0

(2.54)

where ˙Ψ (F, T, ΩΩΩ) can be obtained as

˙Ψ =

∂Ψ
∂F

: ˙F +

∂Ψ
∂T

˙T +

∂Ψ
∂Ω

· ˙ΩΩΩ.

Eq. 2.54 can be rewritten considering Eq. 2.55 as

ρJD = −

∂Ψ
∂T

˙T −

∂Ψ
∂Ω

· ˙ΩΩΩ − ρJS ˙T +

(cid:18)

PPP −

(cid:19)

∂Ψ
∂F

: ˙F ≥ 0.

(2.55)

(2.56)

Eq. 2.56 can be further simplified for the case of isothermal process and considering a physical

expression for the first Piola-Kirchhoff stress PPP = ∂Ψ

∂F , which leads to an inequality of internal

14

energy dissipation D as

−

∂Ψ
∂Ω

· ˙ΩΩΩ ≥ 0.

(2.57)

This inequality must always be satisfied over the body during deformation.

2.3

Incompressible Materials

Resistance of most of elastomers to change volume known as incompressibility condition should

be considered in the constitutive modelling. Therefore, the internal energy function of an incom-

pressible material is modified by

Ψ = ˆΨ (F, T, ΩΩΩ) − p(J − 1),

(2.58)

where J = detF = dV
dV0

= 1 is the incompressibility constraint and p is a Lagrange multiplier to

satisfy the boundary conditions. The rate of the strain energy function Ψ is

˙Ψ =

(cid:18) ∂Ψ (F, T, ΩΩΩ)
∂F

− pF−T

(cid:19)

: ˙F,

(2.59)

where the incompressibility condition is considered as the volume remains constant, ˙V = ˙J = 0.

Substituting Eq. 2.59 to Eq. 2.39, one has

(cid:18)(cid:90)

V0

PdV0 −

(cid:90)

V0

(cid:16)

∂F

ˆΨ (F, T, ΩΩΩ) − pF−T (cid:17)

dV0

(cid:19)

: ˙F ≥ 0.

(2.60)

Eq. 2.61 should be held for every deformation rate, ˙F, which leads to

P =

∂ ˆΨ (F, T, ΩΩΩ)
∂F

− pF−T ,

(2.61)

This equation is valid for the case of

• Incompressible hyperelastic material

15

Figure 2.1: The schematic figure of continuum domain in the different configuration; initial γ0 ,
current γn, and incremental configurationγn+1.

• Negligible contribution of internal energy in deformation

• Isothermal deformation

• Moderate and large range of deformation

2.4 Nonlinear Finite Element Analysis

The macro-scale in the multi-scale framework is assumed as a continuum body ( see Fig. 2.1 ),

which is analyzed using the nonlinear finite element method. Accordingly, by considering the ref-

erence configuration, the governing equilibrium equation of the continuum body can be expressed

as

∇.P + b = 0,

(2.62)

where the P is used to denote the first Piola-Kirchhoff stress tensor, and b represents the body

force in the reference configuration. In this study, the following fonts are used, X for scalar, X

represents vector, X for second-order tensor, and X for fourth-order tensor.

16

x, Xy, YXn+1XnΓnΓ0X0Γn+1(γn+1)Ωn+1ΩnΩ0(γn)(γ0)The weak form of Eqn. 2.62 can be achieved by multiplying a virtual displacement δU to

its both side and then integrating over the domain in the initial configuration Ω0 which could be

written as

(cid:90)

Ω0

δU T (∇.P + b) dΩ0 = 0.

oBy applying the Divergence theorem, Eqn. 2.63 can be rewritten as

(cid:90)

Ω0

∂δU T
∂x 0

PdΩ0 =

(cid:90)

Γt

δU T ¯tdΓ +

(cid:90)

Ω0

δU T bdΩ0,

(2.63)

(2.64)

where ¯t denotes the the applied traction on the surface Γt, and x 0 is the displacement vector in

reference configuration. Generally, the total or updated Lagrangian description can be utilized

to solve the Eqn. 2.64. Defining the state variable at the known step is n, the variable at step

n + 1 is what needed to be determined. Different from the updated Lagrangian description, where

the state variable at n + 1th step is resolved by the previous step n, the state variable in the total

Lagrangian description are obtained at step n + 1 with respect to the initial configuration. The

schematic description of the initial configuration n, and incremental configuration of step n + 1 is

represented in Fig. 2.1. Here, the Lagrangian description is adapted to the analysis of the macro

domain. Hence, Eqn.2.64 can be expressed as

(cid:90)

Ωn

∂δU T
∂xn

PdΩn =

(cid:90)

Γt

δU T ¯tdΓ +

(cid:90)

Ωn

δU T bdΩn,

(2.65)

where, the term ∂δU T /∂xn is the deformation gradient in the updated Lagrangian description

δFn . Hence, the Eqn .2.65 could be rewritten as

(cid:90)

Ωn

δFn : PdΩn =

(cid:90)

Γt

δU T ¯tdΓ +

(cid:90)

Ωn

δU T bdΩn.

(2.66)

Due to the non-symmetric nature of first PK stress P, the nodal forces are generally expressed in

terms of second PK stress S = Fn

−1P. Having the second PK stress, and its work conjugate strain

17

E = 1

2(Fn

T Fn − 1), one can rewrite Eqn. 2.66 as

(cid:90)

Ωn

δEn : SdΩn =

(cid:90)

Γt

δU T ¯tdΓ +

(cid:90)

Ωn

δU T bdΩn.

(2.67)

Furthermore, it should be noted that the Green strain tensor E is the summation of linear and

nonlinear terms, which can be written as

E = EL + ENL = EL +

1
2

Aθ.θ,

(2.68)

where EL and ENL are used to define two-dimensional plane stress/strain problems as

EL =









∂u
∂xn

∂ν
∂n
+ ∂ν
∂xn

∂u
∂yn









,









=

1
2









∂u
∂xn

0

∂u
∂yn

∂ν
∂xn

0

∂ν
∂yn

0

∂u
∂yn

∂u
∂xn































∂u
∂xn

∂ν
∂xn

∂u
∂yn

∂ν
∂yn

0

∂ν
∂xn

∂ν
∂xn

≡

1
2

Aθ.θ

(2.69)









ENL =

1

2( ∂u
∂xn
2( ∂u
∂yn

1

)2 + 1

)2 + 1

2( ∂ν
∂xn
2( ∂ν
∂yn

)2

)2

∂u
∂xn

∂u
∂yn

+ ∂ν
∂xn

∂ν
∂yn

Implementing the standard finite element method in line with the Galerkin discretization approach,

the shape functions N can be applied to approximate the displacement field within each element

as U = N ¯U . Taking variation from Eqn. 2.68 results in

δE = BLδ ¯U + AθGδ ¯U = (BL + BNL) δ ¯U ≡ Bδ ¯U ,

(2.70)

where BL and G are matrices that contain the derivatives of shape functions which can be written

18

in terms of the node j as

BL

j =









∂Nj
∂xn

0

∂Nj
∂yn









0

∂Nj
∂yn

∂Nj
∂xn

and Gj =












∂Nj
∂xn

0

∂Nj
∂yn

0












0

∂Nj
∂xn

0

∂Nj
∂yn

(2.71)

Substituting Eqn. 2.70 into Eqn. 2.67, one can obtain the governing equation of the finite element

analysis by

Ψ (cid:0) ¯U (cid:1) = f int − f ext =

(cid:90)

Ωn

BT SdΩ −

(cid:18)(cid:90)

Γt

NT¯tdΓ +

NT bdΩ

(cid:19)

= 0,

(cid:90)

Ω

(2.72)

where the term Ψ (cid:0) ¯U (cid:1) is the residual force vector which can be divided into two parts namely;

the internal force f int and external force f ext. Besides, due to the nonlinear nature of Eqn. 2.72, it

should be solved iteratively. Hence, the Newton-Raphson scheme is utilized as

Ψ i+1
n+1

(cid:0) ¯U (cid:1) = Ψ i

n+1

(cid:0) ¯U (cid:1) + dΨ ≃ 0,

(2.73)

where, the ith iteration of Newton-Raphson method is denoted using the superscript i, and the

subscript n represents the nth step of simulation process. Besides, the residual force term dΨ can

be delivered with

dΨ =

(cid:90)

Ωn

BT SdΩ +

(cid:90)

Ωn

BT dSdΩ.

(2.74)

Having the material property tensor CSE aligned with the 2nd Piola-Kirchhoff stress and Green

strain and constitutive relation as dS = CSE : dE, one can rewrite the residual force term dΨ in

terms of nodal displacement. Accordingly, by substituting Eqn. 2.70 into Eqn. 2.74, and having

CSE in plane stress/strain condition as C , dΨ can be written as

dΨ =

(cid:18)(cid:90)

Ωn

GT MSGdΩ +

(cid:19)

BT C BdΩ

(cid:90)

Ωn

d ¯U = (Kgeo + Kmat) d ¯U ≡ KT d ¯U ,

(2.75)

19

where KT is the tangential stiffness matrix which can be decomposed into two distinct matrices

namely; the geometric stiffness matrix Kgeo and the material stiffness matrix Kmat. The symmetric

matrix MS comprising the 2nd PK stresses in Eqn. 2.75 can be composed as

MS =






S11I2x2 S12I2x2

S21I2x2 S22I2x2




 ,

(2.76)

where I is the identity matrix. The linearized system of equations can be achieved by substituting

Eqn.2.75 into Eqn.2.73 which results in

KT d ¯U = −Ψ i

n+1

(cid:0) ¯U (cid:1) .

(2.77)

It should be noted that the tangential stiffness matrix is a function of the 2nd PK stress tensor

and material property matrix C , which is obtained by a multi-scale approach from the statistical

mechanics of polymer chains.

20

CHAPTER 3

A Multi physics Constitutive Model to Predict Hydrolytic Aging in Quasi

static Behaviour of Thin Cross linked Polymers

3.1

Introduction

In view of the high resistance to abrasion, corrosion, and chemical degradation, elastomers are

widely used in many sensitive applications to transfer load in high deformation regimes or to dis-

sipate kinetic energy in high frequency. In many industrial applications, it is necessary to provide

a realistic approximation of the nominal service life of a sample subjected to a specific loading

profile. In many cases, elastomers are expected to sustain millions of load cycles or sustain ex-

treme environmental loads for years, if not decades. However, environmental factors along with the

mechanical loads will contribute to chemo-physical damages that will eventually lead to material

failure [21]. Understanding the relation between chemical damages and decay in mechanical per-

formance is critical for the prediction of elastomer failure in extreme applications, such as offshore

platforms. Such decay can be formed due to an individual or combined effects of different degrada-

tion processes [22]. Those mechanisms can be categorized into those with the major impacts such

as thermo-oxidation [23, 24], hygro-thermal [25], hydrolysis [26, 27], deformation-induced aging

[28, 29, 30] , and those with the less severe impact such as chemical corrosion [13], and ozone

cracking [14]. While each of those mechanisms can significantly reduce the life of elastomeric

material in a certain condition, the latter group is not as prevalent in nature as those of the first

group are. The aging process can be formed due to an individual or combined effects of differ-

ent degradation processes such as thermos-oxidation, photo-oxidation, and hydrolysis [31]. While

predicting the individual contribution of those mechanisms has been the center of attention in the

last few years, there is still a long way to go in modeling the damage accumulation in materials

21

exposed to multiple damages simultaneously.

Elastomeric systems are used in a variety of applications ranging from tires to adhesives,

sealants, and thin films [32]. In view of their excellent properties, e.g. abrasion resistance and

durability, elastomeric components can be reliably designed to be subjected to a high number of

cyclic loading in harsh environments. However, their reliability considerably decreases over time

due to the so-called aging process. The aging process, often described by a gradual decay of

the mechanical performance, is a menace to the polymeric components such as joints, adhesives,

sealant,s and fillers that are used in different components in multiple applications.

In the last 15 years, our predictive abilities have been significantly improved, and many theories

have been developed to describe different aging mechanisms, most of which were focused on

thermo-oxidative aging. To this end, Loeffel and Anand [33] investigated a multi-physics model for

the thermal oxidation aging of coatings in severe operating conditions. Some approaches described

the thermo-oxidative aging of a polymer matrix with respect to the chemical reactions, diffusion

and mechanical coupling [34, 35, 36], and some described it using the mechanism of curing [37,

38]. In all the aforementioned studies, the damage has been micro-mechanically modeled by a

coupled approach to link chemical/physical aging to the thermo-mechanical state of the network

[39, 40].

The effect of time on the constitutive behaviour of elastomers is often described by visco-elastic

behavior, which is more of the short-term effects of time [41, 42, 43]. Visco-elastic properties

have been studied extensively both theoretically [44, 45, 46] and experimentally [47]. Ayoub et

al. [48] proposed a constitutive model using a Zener-type framework. Their model predicts the

behavior by integrating the physics of polymer chains and their alteration under cyclic loading.

Furthermore, Khan et al. [49] proposed a phenomenological model to characterize the thermo-

mechanical behavior of visco-elastic polymers.

Here, we mainly focus on the Quasi-static behaviour of samples submerged in water for a long

time (hydrolytic aging) [50], e.g. those used in off-shore plants, biological systems or within the

cooling systems [51]. At this stage, the model cannot take into account the effects of loading

22

frequency as well as visco-elastic effects. In such applications, the predictive models can prevent

catastrophic failure of the systems due to the aging of components. To this end, a model for

hydrolytic aging and its effect on the mechanical properties of rubber-like material can be subject to

broader interest in our community. Hydrolytic aging can be described in two categories: physical

or chemical [52]. Physical aging has no effect on the chemical structure and is a result of the

movement of polymer chains. Chemical aging, on the other hand, is an irreversible process that

changes the morphology of the matrix and directly affects the mechanical performance of the

sample and the probability of premature failure.

Hydrolytic aging has been mostly studied through experimental studies. So far, only a few phe-

nomenological models, and no micro-mechanical model, are available that can successfully take

into account the effect of hydrolytic aging on the behavior of elastomeric materials. Vieira et al.

[53] proposed a model for hydrolytic aging of PLA-PCL fibers used in surgeries, and Breche et al.

[54] who used a non-linear visco-elastic model to predict the response of biodegradable materials

throughout hydrolytic aging. From the experimental side, several studies were conducted to ex-

plore the relation between hydrolytic aging and loss of mechanical performance in different mate-

rials, such as polyester urethane [55], and thermoplastic polyurethanes [56]. On another approach,

many studies were focused on describing the correlation between aging and changes of the matrix

morphology, such as molecular weight [57], cross-link density [58] and formation/detachment of

cross-links [59].

In this work, a new micro-mechanical model is developed to predict the non-linear behavior

of rubber-like materials such as the Mullins effect and Permanent set during hydrolytic aging.

"Mullins effect" refers specifically to the softening observed between the first loading and the sub-

sequent reloadings. However, for the sake of simplicity, many models [60, 61, 62] exclude the

difference between unload and reload, and thus come up with a stabilized stress-strain curve that

has considerably less noise for fitting. The stabilized softening is often referred to as the Idealized

Mullins effect [63], and our model is also developed to predict this behaviour. The fundamentals

of hydrolysis kinetics and assumptions are first discussed in Section 3.2. Then, the experimental

23

study is presented in Section 5.2. Next, in section 5.3, the constitutive model is presented. Af-

ter a parameter study in section 3.5, an evaluation of the proposed model in comparison to our

experimental data and other sets available in the literature is presented (section 6.4). Finally, the

summary and discussion are provided in section 3.7.

3.2 Hydrolysis Kinetics and Assumptions

Hydrolytic aging results from the interaction of elastomer matrix with the hydroxyl(OH−) or

hydrogen(H+) ions in water [64]. In general, this process can occur in any media with OH− or

H+ e.g. acids or alkaline, while the hydrolysis reaction rate is influenced mainly by ions con-

centration. Esters, amide, imide, and carbonates are particularly vulnerable to hydrolysis which

occurs during water attack. Regardless of the location of those groups in macro-molecules, hy-

drolytic attack of water molecules to those groups can cause chain scission which consequently

leads to the production of carboxylic and alcohol end-groups (see Eq.3.1). Thus aging results in

increased concentration of carboxylic end-groups. On the other hand, hydrolytic attack reduces

the molecular weight of the matrix, which consequently, influences the mechanical behaviour.

Kinetic equations are mathematical models to predict and describe a chemical reaction. In this

respect, the hydrolysis process can be best modeled via a first-order kinetic equation with respect

to the carboxylic end groups’ formation rate[57]

(3.1)

d[COOH]
dt

= ζ[Ester][W ater][COOH] = K[COOH],

(3.2)

where

[COOH] = [COOH]0 exp(Kt). Here, the [Ester],[Water], and [COOH] are the concen-

trations of ester, diffused water, and carboxyl end-groups in the polymer matrix, respectively. The

24

H2O+ROCOHOROHCO+parameter ζ is the hydrolysis constant which defines the rate of hydrolytic attack K. Here, t is

the time variable, and subsequently, we define the subscript •0 to represent the magnitude of a

non-kinematic parameter at t = 0, e.g. [COOH]0 represents the initial concentration of carboxylic

group. Since carboxyl end group concentration is inversely related to number-average molecular

weight Mn, one can determine K by measuring Mn during hydrolytic aging using

[COOH] =

1
Mn

⇒ Mn = Mn,0 exp(−Kt).

(3.3)

In view of Eq.3.3, the factors that influence the hydrolysis rate K are the [Ester], [Water], [COOH]

and the hydrolysis constant ζ which itself is a function of temperature, and media’s pH, where

• temperature increases the polymer chain mobility and thus increases the ζ [53]

• PH of the degradation environment affects reaction rates through catalysis. and thus in-

creases the ζ [65, 66, 67]

In this study, ζ is assumed to be constant since the temperature, and the PH were kept constant.

Accordingly, the proposed model in this work is based on the following assumptions

• The process of water diffusion into the sample is considered to be complete and thus d[W ater]

dt =

0. In fact, the proposed model is based on the assumption that the water diffusion time

tdif f inside the sample is significantly smaller than the hydrolytic aging time taging ( i.e.

tdif f << taging) [68]. The proposed model is mainly relevant for super-thin samples where

the assumption of homogeneous diffusion is relevant. Here, we will not investigate in-

homogeneous water diffusion due to the sample thickness. Accordingly, all experiments

used for validation were performed on super-thin samples (thickness< 2mm ) that are stored

for a duration much longer than a few minutes [69].

• Following completion of water diffusion, the water concentration all over the sample is iden-

tical and equal to the maximum saturation level possible [W ater]|S∈V = where V represents

the sample volume.

25

• The model is not valid for short time intervals or ultra thick samples in which ∇[W ater] ̸= 0

[70].

• Likewise, the concentration of the ester groups [Ester] located at the backbone chains are

considered to be almost constant, despite their random scission during the early stages of the

aging process [71].

Satisfying the conditions used for the assumption of constant ζ, [Ester] and [water], one can derive

the setting at which we can assume K to be constant. Obviously, this assumption is not relevant

in cases where ζ, [Ester] and [water] are non-constant, such as non-uniform water diffusion or

varying temperature, and those conditions are not the subject of interest in this work.

3.3 Experiment

Material: Styrene-butadiene rubber (SBR) procured in rectangular sheets of 24” × 12” × 0.125”

were used for validation of the model. The SBR sheets were processed from one batch and pro-

vided by one supplier. Standardized dumbbell shape samples were cut from the sheets with a punch

die (Die C from the ASTM-D412).

Accelerated aging: Samples were fully immersed in several sealed containers filled with dis-

tilled water at temperatures 60°C and 80°C under constant pressure. After a specific aging time,

samples were removed and dried at room temperature (i.e. 23°C) for ten days before characteriza-

tion.

Gravimetric measurement: To measure water uptake, each sample was individually weighed

before and after aging/drying with an accuracy of 0.01mg. Samples were aged in two individual

ovens (i.e. 60°C & 80°C) for 1, 10, and 30 days each. Immediately after aging, the samples were

carefully dried with tissue paper and weighed on the same electronic balance. In this respect, the

percentage of the water uptake WC can be written as

Wc =

Wt − W0
W0

× 100

,

26

(3.4)

a)

b)

Figure 3.1: Water content of samples (a) extent of water absorption versus the square root of
hydrolytic aging time for 60°C aging temperature, and (b) the loss of absorbed water versus drying
time for aged samples.

where Wt and W0 are the sample weights at time t and virgin state, respectively. Fig. 3.1a is the

plot of the water absorption content versus the square root of hydrolytic aging time for the aging

temperature of 60°C. Our data shows that the water uptake is linearly correlated with the square

root of aging time (i.e. Wc ∝

√

t), representing that it follows Fick’s second law [72].

In order to investigate the efficacy of the drying process, the aged samples were kept at ambient

temperature (i.e. 23°C) and their weights Wt∗ were recorded at three different drying times t∗ (

i.e. 3, 5, and 10 days). For each case, all the weights were then compared with the weight of the

virgin sample by substituting Wt with Wt∗ in Eq.3.4. Fig. 3.1b shows the loss of absorbed water

along the drying period for two different accelerated aging conditions. The water maintained in

the sample reduced exponentially with time and after 10 days, the water content nearly reached its

equilibrium value.

Mechanical test: United testing machine SFM-20 with a load cell of 1000 lb. was used for

quasi-static tensile tests. All tests were displacement controlled with the strain rates of 43.29

%
min . The distance between the machine grips was set to 2.71 inches and all the experiments

were performed at room temperature. In monotonic failure tests, the samples were stretched until

breakage while in the cyclic test, the samples were stretched to preset amplitudes of 1.3, 1.6, 1.9,

and 2.1. Each sample is subjected to a non-relaxing cyclic test with increasing amplitude, where

27

02.557.51007001400Water uptake (%)Aging time1/2(s1/2)T = 60 C02.550510Water uptake (%)Drying time (day)Aged at T = 60 C, t = 10 dAged at T = 60 C, t = 1 da)

b)

Figure 3.2: Constitutive behaviour of the dog-bone SBR sample at various aging times shown in
(a) failure, and (b) cyclic tests that were performed at 60°C.

the samples are loaded two times at each amplitude. In the course of deformation, the elongation

of the central zone has been measured by an external extensometer (see Fig. 3.2). The cyclic

experiment is designed to illustrate the evolution of the permanent set and stress softening during

the primary loading for both unaged and aged samples.

Data treatment: The aim of the present work is the modeling of the nonlinear inelastic phe-

nomena such as permanent set, idealized Mullins effect during hydrolytic aging. To this end, the

experimental data must be treated to highlight the characteristics of these inelastic features. To de-

rive the idealized Mullins effect, the following data treatment will be applied to the experimental

data

• The difference between the reloading and unloading curves will be skipped (see Fig. 3.3a).

• The unloading responses are extended to meet the primary loading that is obtained from the

monotonic tensile test. The maximum stress is set to that of the maximum stretch in loading

history.

The modified experimental results for unaged and aged SBR at 60°C are presented in Fig. 3.3b.

In this figure, the dashed lines refer to the result of monotonic tensile tests and are considered as the

primary loading curves. The solid lines are referred to as the secondary curves which correspond

to the unloading curves. Each sample experiences a continuous cyclic test regime with increasing

28

0250500750123Stress (psi)Stretch χUnaged6day aged10day agedAgingχ010020030040050011.52Stress (psi)Stretch χUnaged6day aged10day agedTimeχa)

b)

Figure 3.3: (A) a schematic diagram of idealized Mullins effect, and (b) the treated experimental
results of SBR at 60°C for modeling purposes. dashed lines refer to the result of the monotonic
tensile tests, while solid lines refer to the result of cyclic tests.

stretches in each cycle. The experimental data after the above-mentioned treatment will be used

as the basis for the discussion on the modeling of inelastic phenomena such as idealized Mullins

effect and permanent set in rubber-like material during hydrolytic aging.

3.4 Constitutive Model

Damage in the polymer matrix should be described with respect to three independent variables

temperature T , time t, and deformation which is given by deformation gradient tensor F. To

describe the state of damage with respect to time, one can use the status of the matrix at time 0 and

∞ as the reference points and then define a shape function to interpolate the damage status at time

t. Representing the energy of the matrix in the virgin state by Ψ0 and after indefinite aging time by

Ψ∞, one can approximate the strain energy of the matrix ΨM at time t as

ΨM (t, T, F) = N (t, T )Ψ0(F) + N ′(t, T )Ψ∞(F),

(3.5)

where N (t, T ) and N ′(t, T ) = 1−N (t, T ) represent the weight of each state. To define the changes

of the matrix in the course of hydrolytic aging, one needs to describe the matrix degradation with

respect to the shape function N and the energy states Ψ0, and Ψ∞. By representing the evolution

29

Stress Stretch χabcdeχTimeabcde010020030040050011.52Stress (psi)Stretch χUnaged6day aged10day agedMonotonic testTimeχ(a)

(b)

(c)

Figure 3.4: (a) Shape functions in one dimension using Arrhenius function, and (b) Relaxation
tests designed to capture the changes in the stress during a deformation controlled aging scenario
[24] , and (c) the reconstructed surface of N (t, T ).

of the shape function by the Arrhenius function [73], one can write N (t, T )

N (t, T ) = exp

−γ exp(−

(cid:18)

(cid:19)

)t

,

Ea
RT

(3.6)

where γ is the degradation constant that defines the rate of hydrolytic aging, Ea the activation

energy, and R = 8.314[J]/[mol][K] the gas constant (see Fig. 3.4a).

Experimental Derivation of N (t, T ) can be carried out through specific relaxation tests, in

which samples were aged at constant stretch amplitude for an indefinite time. The changes of the

mechanical response from time 0 to ∞ can be associated with the energy states of the material,

namely Ψ0 and Ψ∞. Having the test results at multiple temperatures, the resulting curves can be

used to derive the 3D surface of the N (t, T ) directly from the experimental data. Similar approach

can be used in other types of aging including, thermo-oxidative or photo-oxidative. For example,

by considering similar tests that were performed by Johlitz et al [24] on thermo-oxidative aging

(see Fig. 3.4b), one can derive the 3D surface of N (t, T ) directly from the experiments (see

Fig. 3.4c). However, in view of the excessive costs associated with relaxation tests, especially in

hydrolysis conditions, we introduced another approach to approximate N (t, T ) based on simpler

experiments and a fitting procedure.

In the matrix of rubber-like materials, the polymer macro-molecules are often cross-linked

or bonded to aggregate surfaces at many points along their length. Thus, each macro-molecule

30

00.250.50.7510∞N(t,T)N'(t,T)Time11.520Stress (Mpa)60 C80 C100 CTime (s)2*105105∞TemperatureN(t, T)01Timehas been often bonded at many locations and forms multiple polymer strands which are named

"chains" throughout this manuscript [74]. Here, the definition of a chain is considered as the

whole or a part of a polymer molecule restricted between two constrained segments, which are

either cross-linked or entangled with other chains or alternatively bonded to the aggregate surface.

Considering that the length of an active chain is measured between two cross-linked regions along

a polymer molecule, the less the cross-link density within a polymeric material, the lengthier the

active chains will be. Accordingly, once water attacks the matrix, two parallel phenomena occur

[58] (i) reduction of the cross-links which increases the average chain length in the matrix, and (ii)

Energy dissipation due to the reduction of polymer active agents (see Fig. 3.5). To this end, a fully

attacked matrix Ψ∞ can be decomposed into two independent networks, namely a newly morphed

network Ψm and a deactivated network Ψd as given below

Ψ∞ = αΨm + (1 − α)Ψd,

(3.7)

where Ψm and Ψd represent the energies of the morphed and deactivated networks, respectively.

The morphed network results from the reduction of the unattacked network cross-links, and con-

sequently, it has chains that are longer than those of the virgin matrix. Moreover, the deactivated

network results from a water attack on polymer chains. Hence, in this network, the chain length

distribution is similar to that of the original polymer matrix. The parameter α defines the con-

tribution of each network. Substituting Eq.5.6 into Eq.5.3, one gets the total strain energy of the

polymer matrix as

ΨM = N (t, T )Ψ0 + α N ′(t, T )Ψm + (1 − α)N ′(t, T )Ψd.

(3.8)

Fig. 3.6 shows the evolution of the polymer matrix with respect to time. At time t = 0, the polymer

matrix can be represented by Ψ0 state only. As time passes, the contribution of Ψ∞ grows which

itself can be decomposed into two networks, m and d. This process continues until the contribution

of Ψ0 becomes negligible and ΨM ≈ Ψ∞.

31

Figure 3.5: Development of longer chains due to chain scission.

Figure 3.6: Evolution of the networks throughout hydrolytic aging (cube dimensions represent the
number of their chains).

32

nP0(n)cross-linkchain scissionProbabilityPm(n)ProbabilitynR0RmR0>t = t0t = t1t = t2t = tN0NFull networkdUnattacked NetworkMorphed Network Deactivated NetworkWe assume that all chains that are exposed to water attack within the deactivated network will

form free-end chains with no contribution to entropic energy, so Ψd = 0. So in Eq.3.8, ΨM is

mainly determined by the energies of the other two networks, namely Ψ0 and Ψm.

Networks and Sub-networks Each network is considered to have a unique composition and de-

scribes a specific energy-dissipating damage mechanism. Using the concept of micro-sphere,[75,

76] each network is considered as a 3D composition of infinite 1D sub-networks that are distributed

in all spatial directions. Sub-networks can only be subjected to uni-axial deformation λ, and thus

will experience different deformations based on their directions (see Fig. 3.6)[77, 78]. To develop

a model for the sub-network in direction d , only a simplified form of entropic energy is needed,

with respect to uni-axial deformation, λd . The parameter ψd represents the energy of the sub-

network in an arbitrary direction, d . Integrating a sub-network in all directions, the consequent

network is a representation of that concept in a 3D configuration. Here, by assuming the isotropic

spatial distribution of the polymer matrix in reference configuration, the macroscopic energy of an

arbitrary network, Ψ•, can be written as

Ψ• =

1
As

(cid:90)

S

d
ψ•d

d
u ∼=

k
(cid:88)

i=1

d i
ψ• wi, ⇒

d
ψM = N (t, T )

d
ψ0 + αN ′(t, T )

d
ψm,

(3.9)

where As is the surface area of a micro-sphere S, and ud the unit area of the surface with the normal

direction d . The parameters wi are weight factors corresponding to the collocation directions

d i (i = 1, 2, . . . , k). A set of k = 45 integration points on the half sphere was found to yield

the best optimization between computational expenses and the resulting error [79, 80]. Fig. 3.6

depicts the composition of the polymer matrix into the three networks 0, m and d to describe

hydrolysis-induced damages.

Probability Distribution Function of a Polymer Chain Assuming all sub-networks of each

network to have the same composition of chains, one can derive a Gaussian probability function

like P•(L|σ, µ•, R•) which holds for both network [81]. Here, L = nl and R• represent the contour

33

length and the end-to-end distance of a polymer chain, while µ• and σ denote the mean and the

can normalize L and R by segment length to achieve normalized length n = L

standard deviation. Considering n to be the number and l to be the length of the segments, one
l , and ¯R• = R•
l .
Moreover, considering ¯R• to be constant in each sub-network, the distribution function can be

abstractly written as P•(n) := P•(L|σ, µ•, R•) where

P• (n) =

1
√
πσ2
2

exp(

(n − µ•)2
−2σ2

),

(3.10)

where the subscript X• represents the network the parameter X is associated with.

Strain energy of a single chain at networks 0, and m Based on the non-Gaussian theory of

rubber elasticity and its approximation through Kuhn-Grün (KG) model, the entropic energy of a

single freely jointed chain is given by

ψc (n, ¯r•) = nKbT

φβ + ln

(cid:18)

(cid:19)

β
sinh β

= nKbT

(cid:90) φ

0

β dφ,

β = L−1(φ) where

coth(β) −

1
β

= L(φ),

(3.11)

where β = L−1(φ) is the inverse Langevin function, φ = r•

L = ¯r•

n the extensibility ratio ( ¯r•

n ), and

¯r• is end-to-end distance of the chain in the deformed configuration . Recent studies show that the

relative error of Kuhn-Grün is significantly large in shorter length (n < 40) [82], and can cause

significant overs estimation of force in highly cross-linked materials. Therefore, the KG is not a

proper option for estimating the entropic energy of a matrix with short chains. To this end, an

enhanced version of the non-Gaussian distribution function of short chains will be used, namely

ˆψc (n, ¯r•) = nKbT

(cid:90) φ

0

ˆβdφ,

ˆβ =

(cid:20)

1 −

1 + φ2
n

(cid:21)

β.

(3.12)

Having the same computational cost as the KG model, the enhanced model can provide us with a

significantly more accurate description of force in shorter chains.

34

ILF Approximation The majority of the micromechanical models either calculate the ILF im-

plicitly or approximate it with rational functions, as it cannot be explicitly derived. In the whole

range of extensibility of polymer chains, it is preferable to approximate the ILF with a highly accu-

rate simple approximation. Thus, a first-order fractional approximation with two polynomial terms

(relative error of 1.0%), is used to approximate ILF [83, 84]

L−1 (x) ∼=

x
1 − x

+ 2x −

8
9

x2.

(3.13)

3.4.1 Damage

Damage in the polymer matrix will take place with respect to two main factors: time t and defor-

mation F. Here, we consider temperature to only maximize the rate of damage rather than being

an initiating factor in the damage. The time and deformation-induced damages are considered to

have taken place in the following forms in networks 0 and m

• Time-induced damage is mainly associated with the dissolution of the cross-links in the

course of aging which results in homogeneous detachment of polymer chains regardless of

their length (random chain scission). Such damage can be best modeled by the changes in

the end-to-end distance ¯R• of a polymer chain (see Fig. 3.7) and the peak location of the

distribution function P• (n) (see Fig. 3.5).

• Deformation-induced damage is mainly associated with the detachment of shorter chains

due to entropic force that exceeds the strength of the weakest link of the chain F f b. Such

damage results in the deactivation of shorter polymer chains and does not have any effects

on other longer chains and thus does not change the distribution profile (see Fig. 3.8 ).

Hence, the shortest available chain in an arbitrary direction of a deformed network is given

by

F

(cid:17)

(cid:16) ¯r•
n

≤ F f b −→ n•min = ν

d
λM

¯R•,

(3.14)

35

Figure 3.7: Illustration of time-induced damage due to change of the network morphology.

(a)

(b)

Figure 3.8: Illustration of deformation-induced Damage. (a) Detachment of shorter chains during
primary loading, and (b) Distribution function and the force developed by chains with an equal
initial relative end-to-end distance ¯R• and pre-stretched by λM .

36

Rut = t0RuRmt = t>AgingUnattacked networkMorphed networkDeactivated networkλ=1λ=λ1λ=λ2>λ1R.λ1R.λ2R.R.MR.MR.where λM

d is the maximal micro-stretch previously reached in the loading history, while

ν = L( F f b
KbT )

−1

> 1 refers to a material parameter representing the sliding ratio of a polymer

chain (see [85] for details). Considering n•M = µ• + 5σ to be the length of the longest avail-

able chain in each network and use it as a cut-off length of P• (n), can significantly reduce

the computational cost of the model and eliminate the negligible contribution of super-long

chains. Consequently, the length of available chains in the direction d of an arbitrary net-

work is given through the following set

(cid:18)d

λM

(cid:19)

(cid:26)

=

n

D•

(cid:19)

(cid:18)d

λM

n•min

(cid:12)
(cid:12)
(cid:12)
(cid:12)

(cid:27)

≤ n ≤ n•M

.

(3.15)

3.4.2 Network Rearrangement

Time-induced Rearrangement

The chain length distribution function P• (n) in each network

is not a function of time and only its multiplicative amplitude, N (t, T ) is evolving over time.

Thus, one can show that the total number of bonds at time 0 and ∞ are identical as a result of the

conservation of the number of active chains in the matrix. However, the chains of virgin network

transformed into the new chains of attacked networks, which can be summarized as follows

(cid:90)

N (0, T )N0

D0(1)

P0 (n) ndn = N ′(∞, T )αNm

(cid:90)

Pm (n) ndn

Dm(1)

+ N ′(∞, T )(1 − α)Nd

(cid:90)

Pd (n) ndn,

(3.16)

Dd(1)

where N• represents the number of chains in an arbitrary network. Here, one should note that the

process of hydrolysis can be defined as the conversion of chains of the virgin network N0 to the

chains of the deactivated and morphed networks αNm + (1 − α)Nd. In other words, initial polymer

chains transform into two states in the course of hydrolytic aging; (i) longer chains formed by the

merger of original detached chains and (ii) detached chains that could not stay active in the system.

Here, we assumed water attacked the part of chains in the virgin network and deactivated them.

37

Figure 3.9: Alteration of the morphed network chain length distribution with respect to the µm.

Thus, the chain length distribution, Pd (n) and number of chains, Nd, in the deactivated networks

are considered to be the same as the ones of the virgin network, Pd (n) ≈ P0 (n) and Nd = N0.

Thus, Eq.3.16 can be rewritten as

(cid:90)

N0

D0(1)

P0 (n) ndn = Nm

(cid:90)

Pm (n) ndn

(3.17)

D0(1)
(cid:82)

P0 (n) n dn

Pm (n) n dn

,

⇒ ˆN =

Nm
N0

=

D0(1)
(cid:82)

Dm(1)

where ˆN is a normalization constant, which satisfies the mass conservation in the networks during

aging time.

In this respect, Fig. 3.9 shows the alteration of the chain distribution function due to the change

of Nm with respect to µm in the morphed network. By increasing µm, Nm decreases and moves

the distribution function to the right and down (See Fig. 3.9).

Deformation-induced Rearrangement Considering that the detached chains still remain part

of longer macro-molecules, two presumptive conclusions can be made; (i) the detachment of the

chains does not lead to full dissipation of their energy, and (ii) the number of active segments within

a network remains constant in the course of deformation (see e.g. To implement this assumption,

38

00.010.020.03050100150200Chain Length DistributionNumber of segments nμm= μ0μm= 2μ0μm= 3μ0Figure 3.10: Effect of chain detachment on the number of active segments.

we consider the detached chains to be evenly distributed among available chains and thus increase

the number of active chains, though with the same probability function P• (n). Accordingly, in the

course of primary loading shorter chains detach and become part of the longer chains ( see Fig.

3.10). In the unloading and reloading, no further damage occurs, as long as λM remains the same.

[85, 86] for details). The network rearrangement can be best quantified through an amplification

factor Φ•(λM ) which amplifies the distribution profiles P• (n) as more chains are detached from

the matrix [87]. This process can take place in m or 0 networks depending on the λM on that

direction. Assuming the number of segments to be identical before and after deformation, one can

write

(cid:90)

N•Φ•(1)

D•(1)

P• (n) ndn = N•Φ•(λM )

(cid:90)

P• (n) ndn.

(3.18)

D•(λM )

The detachment process is irreversible; we consider no new bonds to be formed during reloading

and unloading, and thus one can rewrite Eq. 5.14 by considering Φ•(1) = 1 as

Φ• (λM ) =

(cid:82)

D•(1)
(cid:82)

P• (n) ndn

P• (n) ndn

.

D•(λM )

39

(3.19)

Deactivation of some segmentsNo effect on number of active segmentsActivation of some segmentsχ3.4.3 Strain Energy of Networks

The free energy of networks 0 and m in an arbitrary direction d can be written with respect to Eqs.

5.10, 3.17, and 5.15

d
ψ• =

(cid:90)

(cid:19)

(cid:18) d
λM

Φ•

D•(λM

d)

ψc (n, ¯r•) P• (n) dn,

(3.20)

where D• represents the available chains’ length in •th network.

Micro-Macro Scale Transition is used to define the relationship between the micro and macro

deformation within the polymer matrix. Here, the non-affine directional model is used to relate the

macro-stretch χ, and the micro-stretch λ. Based on this concept, a non-homogeneous distribution

of the stretch within the network domain is assumed. Here, a classical strain amplification function

which is often used in rubber-like materials [88] is used, however the model can be changed for

polymers with more elaborate behaviour.

d
λ =

χd − C p
1 − C p

(3.21)

where C ∈ [0, 1] is the volume fraction of hard segments per unit volume of the rubber matrix, and

p depends on the structure of the network (p ≈ 1/3 for details see [89]).

3.4.4 Constitutive formulation

For an incompressible polymer matrix

det F = 1

(3.22)

where F stands for the macro-scale deformation gradient. The first Piola-Kirchhoff stress tensor P

can be written as

P =

∂ΨM
∂F

− pF−T =

∂Ψ0
∂F

+

∂Ψm
∂F

− pF−T ,

(3.23)

where p denotes an arbitrary scalar parameter that can be defined to assure incompressibility, and

40

Table 3.1: The reference set of parameters of the proposed model.

N0kbT [M pa]
40

¯R0
3

¯Rm µ0 µm σ
18

13

15

4.5

α

ν
1.0065 0.02

γ exp(− Ea
RT )
0.1

∂Ψ0
∂F

∂Ψm
∂F

= N (t, T )

k
(cid:88)

j=1

wj

d
ψ0
d j
λ

∂

∂

d j
λ
d jχ

∂

∂

= αN ′(t, T )

k
(cid:88)

j=1

∂

wj

d
ψm
d j
λ

∂

∂

∂

2
d j
λ
d jχ

1
d jχ

∂d j

¯Cd j
∂ ¯F

∂ ¯F
∂F

,

:

∂d j

¯Cd j
∂ ¯F

∂ ¯F
∂F

,

:

1
d jχ

2

(3.24)

where C denotes the right Cauchy–Green tensor, J 2 = det C, and ¯C = J −2/3C . Detailed

description of the derivation procedure of each of the terms of Eq. 5.3.4 has been provided in

Appendix A.

3.5 Parameters Study

Parameter study is carried out to evaluate the contribution of the physical parameters and validate

the properties associated with them. First, a reference set of values is selected for all parameters

which is given in Table 3.1. The values are specifically chosen to provide high resolution in ma-

terial behavior during the aging time (See Fig. 3.11). To investigate the effects of each parameter

on model predictions, we investigate the increased and decreased values of that parameter while

all others are kept constant. In view of different effects of parameters at different aging times, we

investigated each parameter at three different times, namely t = 0, t = t1, and t = ∞.

The proposed model has ten material parameters, four represents the behaviour of unchanged

network ( ¯R0, µ0, σ, and ν), five those of morphed network ( ¯Rm, µm, α, E and γ), and one pa-

rameter namely N0 just serves as a multiplication factor to describe the response of the polymer

41

Figure 3.11: Model prediction with the selected material parameters during the aging time.

matrix. Due to similarity of the effect of Ea

RT and γ on the response of the material in a constant

aging temperature, here, we only investigate the variation of γ exp(− Ea

RT ) as a rate of aging. In

order to analyze the parameters, first, the effect of unchanged network parameters are investigated

along the aging trajectory. Fig. 3.12 shows a summary of the parametric analysis of ( ¯R0, µ0, σ, ν,

α, and γ exp(− Ea

RT )) during aging time.

• Parameter ¯R0, which represents the end-to-end distance of unchanged network, governs
¯R0 ) in the unattacked network in any
the minimum available chain length (n0min = νλM
arbitrary direction. Thus by increasing ¯R0, n0min increases while the mean length of polymer

chains in the network, which is governed by µ0, remains constant. As shown in Fig. 3.12a,
decreasing ¯R0 decreases the inelastic effects such as the hysteresis ratio. As expected, the
effects of ¯R0 on material response diminishes along aging time since contribution of Ψ0

becomes negligible .

• Parameter µ0 defines the peak of chain distribution, P0(n), and can be represented as a
multiplicative factor of ¯R0. In virgin material, µ0 plays a significant role in defining the
evolution of damage. Accordingly, the damage will increase when ¯R0 reaches the peak of

probability, µ0 (see Fig. 3.12b) as the network loses more chain due to deformation. In

addition, this parameter control the number of chain in the morphed network through the

normalization constant, Nm. Thus, response of the material will change along the aging

42

0501001501.21.8Stress (Mpa)Stretch χt = 0t = t1t = ∞Agingtime. It is shown (Fig. 3.12b) that by increasing the µ0 the softening of material response

decreases at t = 0 and increases at t = ∞.

• Parameters σ defines the variance of of chain distribution, P•(n), and thus have a direct

effect on the amount of damage the material experience during deformation. Regardless of

aging time, higher values of σ represents wider distribution and consequently lower damage

as shown in figures see Fig. 3.12c.

• Parameter v governs the magnitude of energy dissipation in detachment of a chain, and thus

the hysteresis ratio and also the profile of the unloading curves. Regardless of aging time,

lower values of v represents higher energy dissipation per chain and consequently results in

stronger and yet more damageable networks as shown in Figures 3.12d.

• Parameter α defines the composition of the transformed matrix, in particular it is correlated

with the percentage of the polymers that transform into m and inversely correlated with those

that transform into d networks. Considering the zero contribution of the d network in energy

of the matrix, parameter α can be considered as a multiplicative scaling factor for Ψm that

can linearly scale the contribution of morphed network. As expected, the contribution of α

in the final response increases over time as Ψm grows with aging as shown in Fig. 3.12e.

• Parameter γ exp(− Ea

RT ), or the hydrolysis constant, mainly governs the decay rate of the 0

network and formation rate of the m network. Higher values of γ exp(− Ea

RT ), thus, represents

a faster transformation from 0 to m networks (See Fig. 3.12f).

3.6 Validation and Results

To validate the proposed model, its predictions were bench-marked against a set of new experi-

mental data specifically designed to capture the effects of (1) deformation χ, (2) deformation

history χM , (3) aging time t and (4) aging temperature T . Hereafter, t refer to the time that

43

Figure 3.12: Sensitivity analysis of the material parameters included in the model. Each plot
consists of three independent lines. Solid lines represent the reference curve, and dashed lines
depict model predictions due to the variation of the control parameters. In some cases, all three
lines are on top of each other.

44

t = 0t = t1t = 02001.4χR0= 1R0= 3R0= 501001.4χR0= 1R0= 3R= 5001601.4P(Mpa)0= 0.5R0μ0= 2R0μ0= 4R0μ01001.40= 0.5R0μ0= 2R0μ0= 4R0μμ001001.40= 0.5R0μ0= 2R0μ0= 4R0μ01001.4= 3σ= 4.5σ= 9σ10001.4= 3σ= 4.5σ= 9σσ02001.4P(Mpa)= 3σ= 4.5σ= 9σ02001.4P(Mpa)= 1.0005ν=1.0065ν= 1.0013ν02001.4= 1.0005ν= 1.0065ν= 1.0013ν02001.4= 1.0005ν= 1.0065ν= 1.0013ν(a)(b)(c)(d)02001.4P(Mpa)= 0.01= 0.1= 0.5(f)02001.401001.4χ02001.4P(Mpa)= 0.01α= 0.02α= 0.0401001.4= 0.01α= 0.02α= 0.04αα01001.4= 0.01α= 0.02α= 0.04αα(e)γexp(E/RT)γexp(E/RT)γexp(E/RT)γexp(E/RT)= 0.01= 0.1= 0.5γexp(E/RT)γexp(E/RT)γexp(E/RT)= 0.01= 0.1= 0.5γexp(E/RT)γexp(E/RT)γexp(E/RT)Table 3.2: Material parameters of the proposed model for SBR.

N0kbT [psi]
95.4

¯R0
3.0547

¯Rm
3.3849

µ0
7.04

µm
7.24

σ
7.28

ν
1.0055

α
0.001

γ[1/day]
1459

Ea
R [K]
3700

samples were aged before inducing the deformation χ. The model was optimized by fitting the

aforementioned ten material parameters to the following set of load-unload curves

• Primary loading curve of the unaged sample, to study χ

• One unload cycle of the 1.3 stretch amplitude of the unaged sample, to study χM

• Only loading curve of 10 days aged sample at 80°C, to study χ, χM , and t

• One point: stress at 1.5 stretch for 6 days aged sample 60°C, to study t, and T

To this end, the least square error function was minimized with the aid of the Levenberg–Marquardt

algorithm. Once the material parameters are obtained, the model can predict the other loading-

unloading curves of materials for a different amount of aging times and strains. In this regard, the

model predictions are validated against our own experimental tests for different aging times and

temperatures in Fig. 3.13, while the obtained values of the material parameters are given in Table

6.2.

Here, we validate the model predictions on different types of damages as given below

• Fig. 3.13a: Damages induced by deformation, and deformation history in multiple cycles

• Fig. 3.13b: Damages induced by time and temperature in one cycle

• Fig. 3.13c and d: Effect of time on damages induced by deformation, and deformation

history in multiple cycles

• Fig. 3.13e and f: Effect of temperature on damages induced by deformation, and deformation

history in multiple cycles

45

(a)

(c)

(e)

(b)

(d)

(f)

Figure 3.13: Validation of model predictions for SBR in multiple conditions; (a) unaged, (b) Com-
parison of one cycle load on samples stored at 60°C and 80°C for 6 and 10 days, (c) constitutive
behaviour for 6 days aged sample at temperature 60°C, (d) constitutive behaviour for 10 days aged
sample at temperature 60°C, (e) constitutive behaviour for 6 days aged sample at temperature 80°C,
(f) constitutive behaviour for 10 days aged sample at temperature 80°C.

46

010020030040011.52Stress (psi)Stretch χExperimental dataFitPrediction010020030040050011.52Stress (psi)Stretch χUnagedT = 60 C t = 6 dT = 80 Ct = 10 dT = 80 C t = 6 dT = 60 C t = 10 d010020030040011.52Stress (psi)Stretch χExperimental dataFitted pointPrediction010020030040011.52Stress (psi)Stretch χExperimental dataPrediction010020030040011.52Stress (psi)Stretch χExperimental dataPrediction010020030040011.52Stress (psi)Stretch χExperimental dataPredictionFitMoreover, to quantify the deformation-induced damages of polymers in the course of hy-

drolytic aging, the dissipated energy in each cycle in the model is predicted and compared to

those of the experimental data. Throughout this contribution, the term ”Hysteresis” will be used

to refer to "quasi-static" energy dissipation observed in one load-unload cycle in uni-axial tensile

tests. In elastomers, while hysteresis is strongly dependent on deformation rate, even at close to

zero rates, certain hysteresis can be observed which is referred to as rate-independent or "quasi-

static" hysteresis. We generally refer to the hysteresis of the first cycle as idealized Mullins effect,

while those of the second and next cycles are referred to as hysteresis. Please note that this "Hys-

teresis" is different than frequency dependent visco-elastic hysteresis since this one is resulted by

a combination of reversible and irreversible energies observed in quasi-static loading as the result

of sub-structural damages such as bond detachments. At any aging time t and temperature T , the

dissipated energy per unit of volume during the first cycle can be calculated from the experiments

as follows

t,T

Uhys

=

(cid:90)

t,T

P

X ⊗ X dχ −

(cid:90)

t,T

P

X ⊗ X dχ,

(3.25)

loading

unloading

where X is the deformation vector applied on the matrix. Fig. 3.14 show the evolution of dissi-

pated energy against χM for different aging times at 60°C, and 80°C, respectively. As expected,

the dissipated energy increases exponentially with χM while the growth rate is being inversely

correlated to aging time t.

Considering the micro-mechanical nature of the mode, the concept can be further generalized

and implemented for other types of cross-linked polymers with complex inelastic patterns such as

hardening behavior and alteration of curvature. Here, we generalized the model to predict the

hydrolytic aging behavior of three other compounds form experimental data that are available in

literature

• Polyurethane [90]

• Suture fibers of PLA–PCL[53]

• Natural rubber[91]

47

(a)

(b)

Figure 3.14: Evolution of dissipated energy during aging at different temperatures; (a) 60°C, (b)
80°C.

Table 3.3: Material parameters of the proposed model for different Set of rubber-like materials.

Ref.

Gac et al.[90]
Vieira et al.[53]
Stevenson.[91]

N0kbT
[M pa]
0.32
4.88
0.13

¯R0

¯Rm

µ0

µm

σ

ν

α

13.374 17.535
0.3
10.567 21.795 13.04 15.5 4.28 1.01 0.128
6.027 12.582 17.04 19.04 8.28 1.055 0.36

5.1 1.005

5.5

7.5

γ exp(− Ea
RT )
[1/day]
0.0209
0.0443
0.001

To fit the parameters for each material, two loading curves of unaged and aged samples, as well

as one point of an additional loading curve from another sample with a different aging period

were considered for fitting. The parameters for each of the above mentioned materials are derived

and summarized in Table 3.3. In this table, the decay constant γ exp(− Ea

RT ) presented as single

parameter as material responses captured only in one temperature. Fig. 3.15 shows the predictive

capabilities of the presented model against the experiments for all three materials, and the good

agreement for rest of curves were obtained automatically.

Depending on the polymer matrix, we assume that all chains which are exposed to water attack

within the deactivated network will form free-end chains, which will have no further contribution

to entropic energy. Accordingly, the ratio of morphed network to deactivated network α describes

how the matrix elasticity is perceived at time ∞. In this work, we assume that α has been inversely

correlated to Ψ∞, so for α = 0 we expect the material to show no retention force at time ∞, and

48

08016011.52UhysχMPredictionUnaged6 days10 days08016011.52UhysχMPredictionUnaged6 days10 daysFigure 3.15: Validation of model predictions against aging of three different materials; (a)
PolyUrethane (PU) submerged in seawater at 100°C[90], (b) PLA–PCL fiber submerged in phos-
phate buffer solution at 37°C[53], (c) Natural rubber (vulcanized A) submerged in seawater at 40°C
[91].

49

010203023Stress (Mpa)Stretch χUnaged40 days80 days120 days200 daysFitted point PredictionFit(a)02004006008001.21.5Stress (Mpa)Stretch χUnaged2 weeksFitted point8 weeks4 weeks PredictionFit(b)00.81.62.43.21.21.62Stress (Mpa)Stretch χ588 daysUnagedFitted point168 days410 days824 days PredictionFit(c)for α = 1 we expect the system to remain fully elastic. Since SBR loses most of its elasticity due

to hydrolysis, the response can be simply described by transformation of the virgin network to the

deactivated network with almost negligible contribution of morphed network α ≈ 0. Hence, the

energy at infinity state can be considered zero (i.e. Ψ∞ = 0) which shows that the material will

be annihilated through aging. However, for other type of the material such as polyurethane which

retain elasticity at time ∞, α will remain relatively high and Ψ∞ = αΨm.

Using micro-sphere as the basis of deactivated and morphed networks in our model, the pro-

posed model can be independently used to predict aging in different mechanical loading scenarios.

Once the model has been fitted to uni-axial data, its predictions were validated against shear and

uni-axial tensile loads in virgin-state (see Fig. 3.16). While we did not find data on the multi-axial

performance of elastomeric samples under hydrolytic aging, model predictions on the aging of

samples under different loading were predicted but only validated with respect to uni-axial load-

ing. Using the data of [92] on uni-axial and pure shear of filled silicone rubber, the proposed model

validated for different loading scenarios in Fig.3.16. For the fitting, we used one loading-unloading

cycle of the 1.413 stretch amplitude in the uni-axial tensile direction. The good agreement with

all loading-unloading curves in pure shear, as well as other loading-unloading curves of uni-axial

tensile, was attained automatically. The obtained values of the material parameters are given in

Table 3.4.

Table 3.4: Material parameters of the proposed model for filled silicone rubber.

Ref.

N0kbT ¯R0
[M pa]

¯Rm µ0

µm

σ

ν

α

γ
[1/day]

Machado et al.[92] 0.2858 2.0640 4.25 13.04 15.04 6.28 1.004 0.04 1759

Ea
R
[K]
3700

3.7 Summary and Discussion

We have developed a large strain three-dimensional micro-mechanical constitutive model to de-

scribe the hydrolytic aging of a cross-linked polymers matrix with permanent chemical cross-links

50

a)

c)

b)

d)

Figure 3.16: Validation of model predictions for silicone rubber under different loading scenarios
against experimental data from [92]; (a) unaged sample under uni-axial tension, (b) unaged sample
under pure shear, (c) Comparison of one cycle of uni-axial tension load on samples stored at 60°C
and 80°C for 5 and 10 days, (d) Comparison of one cycle of pure shear load on samples stored at
60°C and 80°C for 5 and 10 days.

that can break and reform the matrix with respect to four external load factors of (i) deformation,

(ii) deformation history, (iii) aging time and (iv) aging temperature. Our model show excellent

agreements with the experimental data on five different types of loading and on four different

materials. To the best of our knowledge, this is the first micro-mechanical model of hydrolytic

aging.

The model has been based on the concept that the chemical cross-links can break, can conse-

quently transform the network in the course of aging. We hypothesized two different types of net-

work transformation depending on the source of cross-link breakage, namely time or deformation.

During a test with continuous loading, the deformation-induced damage will lead to detachment

51

00.350.712Stress (MPa)Stretch χPredictionFitExperimental data00.3511.5Stress (MPa)StretchχPredictionExperimental data00.350.712Stress (Mpa)Stretch χUnagedT = 60 C,t = 5 dT = 60 C,t = 10 dT = 80 C,t = 5 dT = 80 C,t = 10 d00.3511.5Stress (Mpa)Stretch χUnagedT = 60 C,t = 5 dT = 60 C,t = 10 dT = 80 C,t = 5 dT = 80 C,t = 10 dof shorter chains and thus results in a rearranged network with slightly longer polymer chains. In

the course of aging and in the absence of deformation, we assume that the chain detachment due

to water attack on polymer active agents occurs randomly and thus the reformed network will have

a different composition from the original network with fewer chains. Rearrangement of the chains

contributes significantly to the overall stress of the cross-linked polymer matrix and was accurately

represented by our model. Although our model has ten independent parameters, two of them can

be immediately eliminated if deformation controlled relaxation aging tests are available. Here, Our

formulation neglects rate dependency and visco-elasticity while the major focus of this work is on

relatively long aging times and low deformation rates where the effects of visco-elasticity is as-

sumed to disappear. Thus, if the deformation rate is sufficiently low one can assume the tests, pre-

and post-aging, as quasi-static tests which makes the unloading curves insensitive to the loading

rate and mainly a function of deformation history.

The proposed model captures the basic physical laws governing the hydrolytic aging of cross-

linked polymers as frequently encountered in nature and thus is relevant to two other types of

chemically cross-linked polymers such as rubber-like materials, adhesives and sealants. For poly-

meric systems with more complex morphology, e.g. those with two or more types of cross-links

or polymers, while the proposed concepts are still relevant, the models describing the dynamics of

breaking of cross-links and network rearrangement should be updated and probably, would be way

more complex than those used in this work. The proposed model has ten material parameters, all

of which have a clear physical meaning. The excellent performance of the proposed method was

proven by validating against different experimental data on different materials that are particularly

selected to reveal the evolution of inelastic behaviour during hydrolytic aging. The model pre-

dictions was further validated by comparing the evolution of hysteresis against experiments. Our

model and experiments illustrate how polymer networks with their mechanical behaviour domi-

nated by entropic chains, can respond to time, temperature and deformation at the micro-structural

level and further explore the relation to describe macroscopic response with respect to micro-

structural changes.

52

CHAPTER 4

A Physically Based Model for Thermo-Oxidative and Hydrolytic Aging of

Elastomers

4.1

Introduction

There is a high demand for using cross-linked polymers in demanding applications due to their ex-

cellent properties, such as high extensibility, strength, abrasion and chemical resistance, and dura-

bility [32, 21, 93, 94, 95]. They are extensively used in a wide range of industrial applications, such

as automotive, transportation, power distribution, and aerospace. Due to their usage in sensitive

applications, they are constantly exposed to high numbers of cyclic loading in unfavorable environ-

ments. However, these extreme environmental conditions can severely degrade their mechanical

performance and hence their reliability [96]. In the literature, this process is often known as aging.

Individual or combined environmental factors such as high temperature, radiation, and humidity

can cause the aging process [5]. These environmental factors facilitate degradation mechanisms.

Based on the prevalence of their environmental factors, aging mechanism can be categorized as

ones with severe impact such as hygro-thermal [97, 98], hydrolysis [99, 58], thermo-oxidation

[100, 101, 102, 103], photo-oxidation [104, 105, 106] and those that are not that prevalent in na-

ture such as ozone cracking [14], and chemical corrosion [13]. Although there is a vast literature

on proposed aging constitutive models to predict the material behavior in one aging scenario, there

is still a long way to go for simple and efficient models that can be used alternately for different

aging scenarios at the same time. Therefore, in this research, we focused on proposing a com-

putationally simple model that can be used for thermo-oxidative aging or hydrolysis separately.

Regarding thermo-oxidative aging, there are extensive experimental studies on cross-linked poly-

mers using accelerated testing. Celina et al. [107, 108] thoroughly reviewed such studies. Several

53

potential problems of accelerated testing methods were discussed in detail in [109, 110, 111, 112].

Moreover, Shaw et al. [113] investigated the effects of thermo-oxidative aging on oxidative scis-

sion, viscoelastic relaxation, and thermal expansion of elastomers. Rabanizada et al. [114] used

air, distilled water, seawater, and freshwater as different mediums to perform accelerated aging

tests on natural rubber. Furthermore, they studied the changes in visco-elastic behavior through

dynamic mechanical analysis (DMA) tests. In a similar fashion, the thermal aging of peroxide-

cured nitrile butadiene rubber (NBR) has been examined by Pazur et al. [115]. Then, they used

their experimental results to study the effect of thermal aging on the bulk properties of the rubber.

They suggested that to understand the decay of material one should study the cross-link distribu-

tion within the bulk specimen in detail. The experimental approaches to understanding hydrolysis

are not as detailed as thermo-oxidative aging, but there is still a rich literature on this subject. In

this regard, experimental research on the mechanical properties of thermoplastic polyurethanes in

various thermal, mechanical, and humidity conditions has been conducted by Slater et al. [56]. The

results revealed that the water acts as a plasticizer and produces a greater magnitude of compres-

sion set. A series of experiments on polyglyconate B for distinct initial polymer molecular weights

have been done by Farrar and Gillson [57] to show that the degradation rate depends on polymer

morphology. Besides, they suggest an empirical relation between tensile strength and molecular

weight. Later on, similar decreasing trends for tensile strength and molecular weight were found

experimentally by Vieira et al. [116], during hydrolytic aging. Compared to the experimental ef-

forts in this field, there are fewer theoretical approaches in the literature. These theoretical efforts

can be mainly categorized into two groups of phenomenological and micro-mechanical models.

Phenomenological approaches describe the degradation of elastomers based on a thermodynamic

framework, along with a mathematical model without a physical meaning [15, 16]. On the other

hand, micro-mechanical models are based on the statistical mechanics of polymer structure. There-

fore, their parameters do have physical meaning [117, 17]. Recently, with the growth of machine

learning techniques in engineering fields, a new data-driven approach to model behavior is coming

into the picture as well. In this regard, Ghaderi et al. [118, 119] proposed infusing machine learn-

54

ing algorithms and polymer physics to model inelastic behavior of rubber-like materials and their

degradation. Phenomenological models are more prevalent than the micro-mechanical ones. In

thermo-oxidative aging, Ha-Anh and Vu-Khanh [120] used the Mooney-Rivlin model along with

an Arrhenius decay to analyze hyperelastic behavior of polychloprene during aging. Furthermore,

Lion and Johlitz [121] introduced a 3-D phenomenological model for the chemical behavior of

rubbers. Helmholtz free energy has been differentiated into three parts, temperature-dependent

hyper-elasticity, volumetric material behavior, and a function of deformation history. Then, these

parts are added to the decomposition of the deformation gradient. Their model was successful in

simulating the experimental data. In addition, Johlitz [122] proposed a phenomenological model

which was rheologically motivated. Later on, he used finite strain theory and simulated the behav-

ior of aged automobile bearings [24]. In order to study the aging behavior of an adhesive, Dipple

et al. [123] used a coupled chemo-mechanical modeling approach. Using finite element simula-

tion, they investigated the dependency of aging to the geometry between substrate and adhesive.

Recently, Konica and Sain [124, 125] numerically studied the coupled mechanical behavior and

diffusion-reaction of polymers during aging. Thermal aging can also be modeled using molecular

dynamics simulations [126], although it requires higher computational costs, implementation of

coarse-graining methods could accelerate these simulations [127]. For hydrolysis, so far, aside

from our previous model, only a few phenomenological models have been proposed [26]. To

model the constitutive behavior of PLA-PCL fibers during hydrolysis, Vieira et al. [53] devel-

oped a phenomenological model and validated it with their experimental data. To this end, they

used Bergström and Boyce’s constitutive model, which decomposes the mechanical behavior into

two separate parts; i) time-independent and ii) time-dependent response. Recently, Breche et al.

[54] predicted the response of biodegradable materials throughout hydrolytic aging by using a

linear-quasi viscoelastic model. Recently, the demand for physically-based models is increasing.

Schlomka et al. [128] used the Neo-Hookean model to model the polymer’s stiffness during aging.

They used the same staggered solution presented in [129] to solve the problem of the time scale

difference between chemical reactions and mechanical deformations. Beurle et al. [130] proposed

55

a micro-mechanical model that uses a network dynamics model to capture the changes in polymer

network topology due to cross-link formation and chain scission. Mohammadi et al. [131, 12] pre-

dicted the constitutive behavior of elastomers during thermo-oxidative aging. Their model used the

dual network hypothesis of Tobolsky for decay modeling alongside the network evolution model

for deformation-induced damage. Their model has the ability to model idealized Mullins effect as

well as follow complex loading, aging, and unloading profiles. In a similar fashion, Bahrololoumi

et al. [9, 132] proposed a physically-based model for polymer behavior undergoing hydrolytic ag-

ing. Although this model is pretty accurate and versatile in predicting different inelastic features, it

is computationally expensive. The main objective of this work is to provide a generalized method

to reduce complexity and computational cost associated with training of our recent aging models

[9, 133, 12]. The proposed model will use a generic description of two competing phenomena

that are generally observed in aging, namely the formation and dissolution of cross-links, which

results in varying densities of active chains throughout the process. The model simplifies the ag-

ing by assuming the mechanical and environmental damages to be independent and thus can be

modeled as multiplicative damage components. Here, we validate the model with respect to two

aging processes, namely hydrolysis, and thermo-oxidation, and we show that it can capture both

mechanisms with relatively good accuracy.

4.2 Kinetics and Assumptions

In order to develop a micro-mechanical model, one should first understand the kinetics of chemical

reactions during thermal and hydrolytic degradation.

4.2.1 Thermo-Oxidative Kinetics

The rate of chemical oxidation during thermo-oxidative aging can be defined as [134]

56

−

d [P ]
dt

= k [P ]q ,

(4.1)

where [P ] is the chemical compound concentration of P, k the coefficient of reaction rate, and q

the reaction order. Note that Eq. (4.1) is only meaningful when q has a discrete value [135]. First-

order kinetic equations ( q = 1) are usually used for describing the chemical reaction governing

a degradation process[136]. Moreover, in the presence of homogeneous conditions, low stretches

during aging, and assuming there is no diffusion-limited oxidation in place, k would follow the

famous Arrhenius function and therefore is only a function of temperature.

k = τ0exp(−

Ea
RT

),

(4.2)

where R = 8.314[JK −1mol−1] is the ideal gas constant, Ea[Jmol−1] the activation energy of

the chemical reaction, and τ0[s−1] a pre-exponential factor. By solving Eq.(4.1) when q = 1 and

substituting it into Eq.(4.2), following equation for describing [P ] can be derived

(cid:18)

[P ] = A0exp

−τ0exp(−

(cid:19)

)t

.

Ea
RT

(4.3)

Where A0 and τ0, in general, are functions of temperature to satisfy Eq. (4.1).

Now, that we understand the kinetics of reactions, we can further expand it to derive our micro-

mechanical model foundations. Chemical reactions between the polymer backbone and oxygen

during thermo-oxidative aging are the main reason for changes in polymer matrix [137]. Two

types of chemical reactions occur which have been represented in Fig. 4.1 ; (1) cross-linking,

which often leads to material embrittlement [138], and (2) polymer scission. Therefore, in view of

Eq. (4.3), one can define the new peroxide cross-link density during thermo-oxidative aging as

57

Figure 4.1: A schematic diagram of chemical aging (i.e. thermal and hydrolytic aging). Thermal
aging results from two simultaneous phenomena; cross-linking, and chain scission while hydrolytic
aging is due to the reduction of cross-links and chain scission.

[CRP ] = CRP∞ − CRP1exp

−τ0exp(−

(cid:18)

(cid:19)

)t

.

Ea
RT

(4.4)

Where [CRP ] is the peroxide cross-link concentration, CRP∞ is its concentration in a totally aged

stage and CRP1 is the difference between the initial and final concentration of cross-links.

4.2.2 Hydrolysis Kinetics

Similar to thermo-oxidative aging, hydrolysis is generally modeled via first-order kinetics equa-

tions, but with respect to the carboxylic end groups’ formation rate

58

FillerCross-linkActive siteChainsCross-linkingScissionThermal AgingHydrolytic AgingScissionReduction of cross-linksd[C]
dt

= ζ[C],

(4.5)

where [C] is the concentration of carboxyl end groups in the polymer matrix. The parameter ζ

is the hydrolytic degradation rate, which mainly defines the speed of aging. In this respect, me-

chanical stress, the ester group condensation, polymer structure, water diffusion, temperature, and

the degradation environment’s pH are several factors that affect the degradation rate. Accordingly,

implementing mechanical stress during hydrolytic degradation can raise the probability of chain

cleavage and subsequently increase the hydrolysis rate constant ζ. Moreover, temperature esca-

lates the ζ as it increases the polymer chain mobility. Similarly, the degradation environment’s pH

can affect the hydrolysis rate constant through catalysis. In this study, the hydrolysis rate ζ re-

mains constant during the degradation process due to the absence of mechanical loading, constant

temperature, and constant domain pH. Besides, the time for diffusing the water into the sample’s

volume is negligible in comparison to the hydrolytic aging. In fact, the water spread out in the

sample’s volume homogeneously. Hence, the concentration of the water molecules is assumed to

be constant over the degradation time. In this work, to simplify the formulation, the degradation

rate is assumed to be constant in the course of aging for a specified temperature T . Hence the

hydrolytic degradation rate can be rewritten as

ζ = γ exp(−

Ea
RT

),

(4.6)

where γ is the aging constant that describes the rate of hydrolytic aging. Since carboxyl end

groups’ concentration is inversely related to number-average molecular weight Mn, one can write

the evaluation of the molecular weight along the hydrolytic aging trajectory as

Mn = Mn0 exp(−γ exp(−

Ea
RT

)t) + Mn∞,

(4.7)

where Mn0 is the initial value of the number-average molecular weight and Mn∞ is its final value.

59

Figure 4.2: A schematic diagram of idealized Mullins effect (considerable softening between first
loading and unloading). The idealized Mullins effect, which is generally regarded as stabilized
softening, excludes the difference between unloading and reloading.

4.3 Constitutive Model

Generally, constitutive models of rubber-like materials are divided into two groups, namely phe-

nomenological and micro-mechanical. The difference between these two groups lies in the source

of construction of their strain energy function W (F). Regarding the phenomenological models,

W (F) is often expressed as a function of strain invariant terms (i.e. I1, and I2). On the other hand,

the micro-mechanical models developed constitutive behavior of the material based on the poly-

mer physics [74, 139]. The aim of this section is to build a physically-based constitutive model for

cross-linked polymers, which is able to take into account the nonlinear inelastic mechanical behav-

ior such as idealized Mullins effect along with each type of chemical aging (i.e. thermo-oxidative,

and hydrolytic aging). The idealized Mullins effect, which is often referred to as stabilized soften-

ing, neglects the difference between unloading and reloading (see Fig. 4.2).

4.3.1 Statistical mechanics of polymers

Experimental observations revealed that elastomers constructed out of the long polymer chains

which are randomly scattered in different directions [140, 141]. A single polymer chain consists of

60

Stress Stretch χabcdeχTimeabcdea)

b)

Figure 4.3: Schematic figures of a) a single polymer chain consisting of n segments with length l,
and contour length nl and b) the unit micro-sphere with the orientation vector d = d1e 1 + d2e 2 +
d3e 3.

several monomers which are constrained by two linkages at the end (see Fig. 6.6a) each of which

might be resulted either from chemical or physical cross-linking. In order to duplicate mathemat-

ically, the morphology of the polymer network is assumed as a 3D composition of unlimited 1D

chains that are distributed in all available spatial directions (for details, see [117]). As a result, the

macroscopic free energy of polymer network W (F) is expressed as the sum of microscopic strain

energies of all active chains available within the network, which can be numerically calculated by

integration over the unit sphere (see Fig. 6.6b). Within this concept, a freely rotating chain (FRC)

consists of a certain number n, and length l of Kuhn segments with an end-to-end distance of

√

nl

and a density of ρd

r . Due to the one-dimensional nature of sub-networks, they can only experience

uni-axial stretches λd =

d TFTFd and accordingly, they will be subjected to unique deforma-

√

tions due to their different directions [77]. Hereafter, the following font styles are used: for a scalar

X, a vector X , a second-order tensor X, and a fourth-order tensor X. Moreover, the bar sign over

the parameters ¯X = ¯X

l refers to their normalized value with respect to the Kuhn length l.

3-D generalization Integrating the entropic energy of a sub-network in all directions, the macro-

scopic energy of an arbitrary three-dimensional matrix can be written as

61

lrL = nlddAa)

b)

c)

Figure 4.4: Polymer network models employed for modeling the chemical aging. a) the three-
chain model, b) the eight-chain model, and c) the micro-sphere model.

W (F) =

1
4π

(cid:90)

S

ρd
r

d
ψd

d

A ∼=

ρd
r wi

d i
ψ,

k
(cid:88)

i=1

(4.8)

where the parameter
d
A is the unit area of the surface, normal to the direction d . The parameters

d
ψ serves as the entropic energy of a single chain in an optional direction,

d . The parameter

wi are weight factors corresponding to the collocation directions d i (i = 1, 2, . . . , k). A different

number of integration points on the half sphere yields different polymer network models with

distinct properties. By choosing three points on the half sphere, the resulted network will be a

three-chain model (see Fig. 4.4a) while by integrating on a four-point, the eight-chain model will

be produced (see Fig. 4.4b). The eight and three-chain models can only model the idealized

Mullins effect among all the inelastic phenomena. However, a set of k = 21 integration points on

the half sphere can be used to capture more complex inelastic phenomena such as deformation

induced anisotropy (see Fig. 4.4c).

Energy of a single chain Based on an enhanced version of non-Gaussian statistics (see [142] for

details), the strain energy density of a freely rotating chain is written as

ˆψc (¯n) = ¯nKbT

(cid:90) T

0

ˆβd ˆT + ψ0,

(cid:20)

1 −

ˆβ =

1 + T 2
¯n

(cid:21)

β.

(4.9)

62

diwhere β = L−1(T ) is the inverse Langevin function, T = λd
√

¯n the extensibility ratio in the de-

formed configuration, Kb and T denote Boltzmann’s constant and the absolute temperature, and

ψ0 is the entropy of the chain in reference state. An appropriate approximation approach for the

inverse langevin function can be selected depending on the elongation range of polymer chains.

Here, a simple and accurate approximation with first-order rational function (relative error of 0.9%)

is used (see [84] for details)

4.3.2 Damage

L−1 (x) ∼=

x
1 − x

+ 2x −

8
9

x2.

(4.10)

The main sources of damage in the polymer matrix are, namely: i) mechanical, and ii) environ-

mental. Here the main assumption of our recent model is that the nature of these two phenomena

is independent of each other, and as a result, their effects on polymer matrix can superpose on

top of each other. Hence, in order to take into account the chemical degradation along with me-

chanical loading, the deterioration of the polymer matrix should be characterized by three external

variables; (i) storage temperature T , (ii) storage time t, and (iii) deformation F. All the aforemen-

tioned parameters are similar for thermo-oxidative and hydrolytic aging.

• In the course of mechanical, the induced damage results from the debonding of shorter chains

from cross-links or aggregates surface during the primary loading. By stretching the mate-

rial, the junction between cross-links and polymer chains will be broken. Accordingly, the

Mullins effect can be regarded as the evolution of the polymer macro-molecules through the

maximum micro-stretch in an arbitrary direction d previously reached in the loading his-

tory λd

m = max

τ ∈(−∞,t]

λd (τ ) . In this respect, Marckmann et al.[143] considered the evolution

of the polymer chains as a result of two phenomena; (i) active chains’ length ¯nd increase

which are involved in polymer network elasticity, and (ii) number of active chains per unit

of volume ρd

r decrease. Moreover, in the course of deformation-induced matrix transition,

they assumed that the number of active polymer monomers per unit of volume would remain

almost constant throughout the deformation. Based on the network alteration theory, the

63

mechanical damage in each direction d can be best modeled via an exponential function of

the maximum deformation endured previously by the material λd
m:

¯nd = n0 exp (cid:0)µλd

m

(cid:1)

and ρd

r = ρr0 exp (cid:0)−µλd

m

(cid:1) ,

(4.11)

where n0, ρr0, and µ are the material parameters that will be determined through the fitting

process. Due to the assumption that the number of active segments remains constant during

the deformation, the rate of the exponential functions µ in both of the equations represented

in Eq.4.11 are equal.

• During the chemical aging without the presence of mechanical loading, damage is mainly

related to the alteration of the number of active polymer chains resulting in a change in the

morphology of the polymer matrix. Accordingly, in order to see the effects of the storage

time and temperature on the number of available active chains in each direction ρd

r , one can

assume that the alteration of the chain density along the aging trajectory is equal to the perox-

ide cross-link density for thermo-oxidation and the change of the average molecular weight

Mn for hydrolysis. Hence, the modified chain density function in the case of hydrolytic or

thermal oxidation aging is given by

r = ρr∞ exp (cid:0)−µλd
ρd

m

(cid:1)






+

ρr0 exp (cid:0)− exp(− Ea

RT )t − µλd

m

(cid:1) ,

ρr0 > ρr∞

if hydrolytic Aging

−ρr0 exp (cid:0)− exp(− Ea

RT )t − µλd

m

(cid:1) , ρr0 < ρr∞

if thermo-oxidative Aging

(4.12)

In view of the reduced computational cost of the proposed model in comparison to [9, 133, 12],

the proposed model has certain limitations

• Model needs to be fitted to multiple loading types; otherwise the predictions obtained by

fitting to uni-axial tension only, may not be accurate for bi-axial loading.

64

• The model is relevant only for compounds in which mechanical and environmental damages

are entirely independent of each other.

• Model is built on the presumption of homogeneous and consistent oxygen/water absorption

throughout the material, which is only the case in thin samples that are aged for a sufficiently

long time.

4.3.3 Constitutive formulation

For an incompressible polymer matrix

where F stands for the macro-scale deformation gradient. The first Piola-Kirchhoff stress tensor P

det F = 1

(4.13)

can be written as

P =

∂W (F)
∂F

− pF−T ,

∂W (F)
∂F

=

k
(cid:88)

j=1

∂

wjρd j
r

d j
ψ c
d j
λ

∂

∂d j

¯Cd j
∂ ¯F

∂ ¯F
∂F

:

1
d j
λ

2

(4.14)

where parameter p guarantees the incompressibility, C stand for the right Cauchy–Green tensor,

J 2 = det C, and ¯C = J −2/3C. The Eq.4.14 can be even more simplified by using the following

procedure

∂ψc

(cid:0)¯n, λd (cid:1)
∂λd

√

=

¯nKBT ˆβ

(cid:18) λd
√
¯n

(cid:19)

, ¯n

,

∂d ¯Cd
∂ ¯F

∂ ¯F
∂F

:

= 2 ¯F(d ⊗ d ) : J − 1

3 I = 2J − 1

3 ¯F(d ⊗ d ).

In the Eq. 4.15, ˆβ is the modified version of nonlinear Langevin force which is suitable for

(4.15)

65

polymers including short chains. Thus, Eq.4.14 yields

P = KBT

√

ρd i
r

k
(cid:88)

i=1

¯nd iL−1





d i
λ
¯nd i

√





wi
d i
λ

J − 1

3 ¯F(d i ⊗ d i) − pF−T .

(4.16)

4.4

Finite Element Linearization

In what follows, the algorithmic framework of the devised model which is suitable for an FE imple-

mentation is presented. To this end, first, the second Piola-Kirchhoff stress S and consistent tangent

modulus C are derived. Then, The in-compressible plane stress condition in total Lagrangian finite

strain is formulated. Since the system of the equations is nonlinear, the Newton-Raphson method

is utilized to iteratively solve the problem [144]. Hence, the second Piola-Kirchhoff stress should

be linearized, and it could be written as

δS = C : δE,

(4.17)

where E is the Green-Lagrange strain tensor (i.e. E = 1

2(C − I)). The second Piola-Kirchhoff

stress S is the first derivatives of the strain energy W (F) with respect to Green-Lagrange strain

tensor E. By considering incompressibility condition, the second Piola-Kirchhoff stress S could

be written as

S =

∂W (F)
∂E

− pJC−1 = 2

∂W (F)
∂ ¯C

:

∂ ¯C
∂C

− pJC−1,

(4.18)

where ∂ ¯C

∂C is a fourth order deviatoric projection tensor in Lagrangian description and could be

written as

∂ ¯C
∂C

= J

−2
3

(cid:18)

I −

C ⊗ C−1

(cid:19)

,

1
3

(4.19)

where I is the fourth order identity tensor. By substituting Eq. 4.19 in Eq. 4.18, one can rewrite

66

the second Piola-Kirchhoff stress as

S = 2J

−2
3

∂W (F)
∂ ¯C

−

2
3

−2
3

J

∂W (F)
∂ ¯C

: (cid:0)C ⊗ C−1(cid:1) − pJC−1.

(4.20)

Here, the plane stress assumption is used to determine Lagrangian multiplier. In this respect, the

Cauchy stress σ in the out-of-plane direction is set to zero. Accordingly, by defining Sn = 2 ∂W (F)
∂ ¯C ,

the hydrostatic pressure p can be obtained as

p = J

−5

3 Sn

33F2

33 −

−5

3 Sn : C.

J

1
3

(4.21)

By substituting hydrostatic pressure p in the Eq. 4.21, one can write the second Piola-Kirchhoff

stress tensor in the case of plane stress as

S = J

−2

3 (cid:0)Sn − Sn

33F2

33C−1(cid:1) .

(4.22)

The Lagrangian elasticity tensor C is derived thanks to the differentiation of the second Piola–Kirchhoff

stress tensor S with respect to the Green-Lagrange strain tensor E :

C =

∂S
∂E

= 2

∂S
∂C

= 2

∂S
∂ ¯C

:

∂ ¯C
∂C

.






0 0

0 0

0

0

0 0 F2
33






, and Z′ =






0 0

0 0

0

0

0 0 2F33






(4.23)

the Lagrangian tangent

By defining Cn = 4 ∂2W (F)
∂ ¯C2

, Z =

modulus is calculated:

67

C =

−2

3 (cid:0)Sn − Sn

33F2

33C−1(cid:1) ⊗ C−1

−2
3
−4
3

J
(cid:18)

+ J

Cn −

Cn : C ⊗ C−1 − C−1 ⊗ Cn : Z +

1
3

(C : Cn : Z) C−1 ⊗ C−1

(cid:19)

(4.24)

1
3
(cid:19)

− 2Sn : ZC−1 ⊗ C−1

,

− J

−2
3

(cid:18)

2Sn : Z′C−1 ⊗

∂F33
∂C

where

∂F33
∂C

=

∂F33
∂F

∂F
∂C

=

1
2

F−1F′, where F′ =

∂F33
∂F

= F2
33






−F22 F21

0

0


0




0

F12 −F11 0

.

(4.25)

4.4.1 The Plane-Stress Condition

The general 3-D incremental stress-strain relation in Eq. 4.17 can be simplified for a certain special

case, namely plane stress. In the plane stress case, all stresses are parallel to the xy-plane and

normal to the z-axis and therefore all the out of plane shear component of the stress tensor are

equal to zero ( i.e. δS13 = δS31 = δS23 = δS32 = 0) . By using Voigt-Mandel notation Eq. 4.17

can be rewritten as






δS11

δS22

δS12

δS33











=

C1111 C1122 C1112 C1133

C2211 C2222 C2212 C2233

C1211 C1222 C1212 C1233

C3311 C3322 C3312 C3333











δE11

δE22

2δE12

δE33






.

(4.26)

In Eq. 4.26, the strain in the out-of-plane direction δE33 is not an independent variable. As

δS33 = 0 and C3333 = 0, the static condensation cannot be applied. Hence in order to remove δE33

from the system of equation, it should be written in terms of in-plane strains ( i.e. δE11, δE22, and

68

δE12). By knowing

δE33 could be written as

E33 =

1
2

(C33 − 1) =

1
2

(cid:0)F2

33 − 1(cid:1) ,

(4.27)

δE33 =

δE33
δE

: δE =

δE33
δF33

δF33
δC

δC
δE

: δE =

1
2

F33F−1F′ : δE.

(4.28)

By defining A = 1

2F33F−1F, the plane stress Lagrangian elasticity tensor C could be written as

C =






C1111 C1122 C1112

C2211 C2222 C2212

C1211 C1222 C1212






+



A11C1133 A22C1133 2A12C1133

A11C2233 A22C2233 2A12C2233

A11C1233 A22C1233 2A12C1233






.

(4.29)

Finally, the detailed derivation of Cn and Sn are given in the appendix.

4.5 Validation and Results

4.5.1 Comparisons with our experimental data

In order to evaluate the capability of the devised model, its predictions were bench-marked against

our experimental results, which particularly designed to show the effect of chemical aging on the

constitutive response of polymers. In this respect, a universal testing machine (TestResources 311

Series Frame) was utilized for uni-axial tensile tests. All the experiments were displacement con-

trolled under the quasi-static condition with a cross-head speed of 50 mm

min . The vertical gap between

the extensometer grips was fixed to 25.4 mm and, all the tests were operated at room conditions

(i.e. 22 ± 2°C, 50 ± 3% RH). Each sample was stretched until breakage, and all the aging condition

has been repeated with five samples for reliability control. In this study, two different cross-linked

polymers, namely PDMS (Polydimethylsiloxane)-based adhesives, and Styrene-butadiene rubber

(SBR) was used to evaluate the generality of the model for different polymer types. Hence, in the

69

case of the hydrolytic aging, the PDMS was used due to its high sensitivity to the water attack

while in the thermo-oxidative aging, the SBR was used owing to its susceptibility to the thermal

degradation. The experimental procedure was as follows. First, all the samples were matched the

standard dimensions mentioned in ASTM D412 (Die C). Next, in the case of hydrolytic aging,

all the samples were entirely submerged in several closed buckets filled with de-ionized water at

temperatures 60°C, 80°C, and 95°C, while in the thermo-oxidative aging samples were put in the

ovens with the relative humidity of zero (i.e. RH= 0%) at temperatures 45°C, 60°C, and 80°C.

All the specimens were aged with a constant pressure, and after a certain amount of aging time,

specimens were removed from containers and immediately dried with tissue paper. Using the six

material parameters, the model was fitted for the following three cases:

• The primary loading curve of the virgin specimen to investigate the effects of λ and λm

• one point of loading curve of 10 days aged specimen to study the effects of λm, and t

• one point of loading curve of 30 days aged specimen to investigate the effects of t, and T

To this end, Levenberg–Marquardt algorithm was used to minimize the least square error function.

The so-obtained parameters of the devised model for silicon adhesive during hydrolytic aging

are given in Table 4.2 while the material parameters of SBR during the thermo-oxidative aging

are given in Table 4.1. The predictions of the proposed model against the experimental data for

different types of aging and different amounts of aging times, temperatures, and deformations

are plotted in Figs. 5.15a-e and Figs. 4.6a-e. It should be noted that for the case of unloading

experiments, we used the concept of data treatment (see [61] for details). Moreover, the changes

in the loading-unloading response of the silicon adhesive aged at 80°C along with the hydrolytic

aging are compared in Fig. 5.15f, while the alteration of the cyclic response of the SBR during

thermal degradation at 60°C is plotted in Fig. 4.6f. Due to the simplicity of the current model,

it cannot capture complex behaviors such as curvature change, so this is the best fit that can be

achieved. In fact, it is a trade-off between the accuracy and simplicity of the model. Furthermore,

during aging experiments, there are different sources of error, such as diffusion-limited oxidation

70

Table 4.1: Material parameters of the proposed model for SBR subjected to the thermo-oxidative
aging.

ρr0KbT [M pa] ρr∞KbT [M pa] n0

µ

1.4416

2.2072

2.497 0.52

γ[1/d]
4.09875e+18

Ea
R [K]
1.6148e+04

Table 4.2: Material parameters of the proposed model for silicon adhesive subjected to the hy-
drolytic aging.

ρr0KbT [M pa] ρr∞KbT [M pa] n0

0.1594

0.0106

2.497

µ
0.52

γ[1/d]
4.387e+14 1.323e+04

Ea
R [K]

(DLO), incomplete water diffusion, and over-curing. Since considering all of these phenomena in

the model will result in a more complicated model, which would not be computationally efficient,

we neglect their effects in our model. Therefore, facing errors as high as 10% (for instance in the

case of 10 days at 95oC) is still acceptable.

Compared with our previous models, our recent model reduces the number of parameters sig-

nificantly which is mainly due to the change of the model’s kernel and the correlation of mechanical

and environmental damages. In our previous models, the network evolution model was utilized as

the kernel of the models while in this paper, the network alteration model is used to model the

rubber inelasticity. In addition, here we assumed that the mechanical and environmental damages

are independent while in our previous models, these damages were coupled to each other. To

demonstrate this reduction quantificationally, four material parameters were reduced in the case of

our hydrolytic aging model while four and six material parameters were reduced for our first and

second models in the thermo-oxidative aging, respectively.

To demonstrate quantificationally the computational efficiency of the proposed model, its per-

formance was compared to our recent models. In this respect, Table 4.3 shows the computational

cost of each model for calculating the second Piola-Kirchhoff stress during one loading-unloading

cycle. The computational cost of our devised model generally depends on the number of inte-

gration points on the micro-sphere. By choosing the four-point on the half-sphere, the model can

71

a)

c)

e)

b)

d)

f)

Figure 4.5: Validation of model predictions against experimental results for silicon adhesive in
unaged state, and submerged in distilled water for 10 and 30 days in multiple temperature condi-
tions; (a) 60°C , (b) 80°C, (c) 95°C, (d) constitutive behavior for unaged sample, (e) constitutive
behavior for 30 days aged sample at temperature 80°C, and (f) Comparison of one cycle of loading-
unloading on samples stored at 80°C for 10 and 30 days. As it can be seen, the devised constitutive
model is able to predict the softening caused by hydrolytic aging. Solid lines represent model pre-
dictions while dot points depict the experimental data. The arrows demonstrate the curves and the
points which are used for fitting.

72

00.511.522.5123456Stress (Mpa)Stretch χFitted curveFitted pointExperimental DataUnagedAged at T = 60C, t = 10dAged at T = 60C, t = 30d00.511.522.5123456Stress (Mpa)Stretch χExperimental DataUnagedAged at T = 80C, t = 10dAged at T = 80C, t = 30d00.511.522.5123456Stress (Mpa)Stretch χFitted pointExperimental DataUnagedAgedat T = 95C, t = 10dAged at T = 95C, t = 30d00.511.51234Stress (Mpa)Stretch χExperimental DataUnaged00.5123Stress (Mpa)Stretch χExperimental DataAged at T = 80C, t = 30d00.511.5212345Stress (Mpa)StretchχUnagedAged at T = 80C, t = 10dAged at T = 80C, t = 30da)

c)

e)

b)

d)

f)

Figure 4.6: Validation of model predictions against our experimental data for SBR in unaged state
and during thermo-oxidative aging for 10 and 30 days in multiple temperature conditions; (a) 45°C
, (b) 60°C, (c) 80°C, (d) constitutive behavior for unaged sample, (e) constitutive behavior for 30
days aged sample at temperature 60°C, and (f) Comparison of one cycle of loading-unloading on
samples stored at 60°C for 10 and 30 days. As it can be seen the devised constitutive model can
predict the hardening effect caused by thermo-oxidation. Solid lines represent model predictions
while dot points depict the experimental data. The arrows demonstrate the curves and the points
which are used for fitting.

73

01.534.5123Stress (Mpa)Stretch χFitted pointFitted curveExperimental DataUnagedAged at T = 45C, t = 10dAged at T = 45C, t = 30d01.534.5123Stress (Mpa)Stretch χFitted pointExperimental DataUnagedAged at T = 60C, t = 10dAged at T = 60C, t = 30d01.534.5123Stress (Mpa)Stretch χExperimental DataUnagedAged at T = 80C, t = 10dAged at T = 80C, t = 30d0123412Stress (Mpa)Stretch χExperimental DataUnaged00.511.522.512Stress (Mpa)Stretch χExperimental DataAged at T = 60C, t = 30d01.534.567.5924Stress (Mpa)Stretch χUnagedAged at T = 60C, t = 10dAged at T = 60C, t = 30dreduce the run time of simulation by almost 400 times in the case of hydrolytic aging and 600 times

in thermo-oxidative aging.

Table 4.3: Comparison of computational costs of our recent models. In the table below, characters
PM4, PM21, and PM45 represent the proposed model with 4, 21, and 45 Gaussian points on the
half sphere respectively. Moreover, the character HM stands for our hydrolytic aging model [9]
while TM1 and TM2 depict our first [133], and second [12] models for thermo-oxidative aging
receptively.

Hydrolytic aging
Model
Run time (sec)
Thermo-Oxidative aging
Model
Run time (sec)

PM4
0.209

PM21 PM45 HM
0.287
0.262

84.87

PM4
0.209

PM21 PM45 TM1
0.287
0.262

121.834

TM2
236.401

4.5.2 Numerical simulation results

In order to demonstrate the capability of the proposed method in modeling the mechanical response

of cross-linked polymers in large deformation along with aging, two numerical examples are pre-

sented. Indeed, the objective of this simulation is to assess the accuracy of the proposed model

in predicting complex strain fields combined with either hydrolytic or thermo-oxidative aging. In

both examples, the parameter values of the previous section are utilized.

• In the first example, a classical dog-bone specimen is modeled under prescribed displace-

ments along the direction. Its length is equal to 115 mm, thickness is 3.75 mm, and the

extensometer’s gauge length is 25.4 mm (see Fig. 4.7a for details). It is meshed with linear

elements in displacement. Due to its thickness, the plane stress assumption is adopted. The

strain–stress response predicted by the simulation in the central section is identical to the an-

alytical response for the different amounts of deformations, aging times, and temperatures.

Next, in order to see the effects of aging on the 3D state of stress, the contours of the second-

Piola Kirchhoff stress S11 are plotted. In hydrolysis, the second-Piola Kirchhoff stress is

74

reduced by the aging time (see Fig. 4.9) while in thermal degradation, it is vice-versa (see

Fig. 4.8). However, in both aging mechanisms, the induced deformation (i.e. mechanical

damage) causes increases in the stress.

• In the last example, the capability of the proposed constitutive model is presented by mod-

eling of a polymer sheet with large in-homogeneous deformations. A polymer sheet with

the dimensions 62 × 82.5 × 1.75 mm and five 20 mm diameter holes C1, C2, C3, C4, and

C5 is modeled under prescribed displacements along-y direction. The detailed positions of

the holes on the sheet are given in Fig. 4.7b. Here, the eight-chain model is utilized to plot

the contour of chain density ρd

r . In fact, due to the symmetry of the eight-chain model, all

the directions d experience the same deformation (i.e. λd =

(cid:113) I1

3 ), and as a result, the

directional chain density ρd

r is the same for all directions and it can be a plottable quantity.

Fig.4.10 shows the evolution of chain density, which is proportional to the number of active

monomers during the hydrolytic and thermo-oxidative aging. Accordingly, the smaller the

chain density in the network, the greater the stress-softening effect. Here, the chain density

is a function of the first invariant of the deformation and aging time. Accordingly, due to

the homogeneous aging assumption, the chain density will be the same over the entire sheet

due to the effects of aging time while deformation causes the in-homogeneous distribution

of chain density especially near the edges of the holes in the rubber sheet.

These two numerical results validate the implementation of the presented constitutive model and

illustrate the independency of mechanical and environmental damage. Such simulations are crucial

to show the capacities of structures in which the level of stress-softening through environmental

and mechanical damages is noticeable.

4.6 Conclusion

The focus of this study was to reduce the number of model parameters without sacrificing the

precision of the model to remedy the complexity of previous models. Here, a large strain three-

75

a)

b)

Figure 4.7: Detailed geometry of the samples used for FE (a) dog-bone sample and (b) rubber
sheet with appropriate coordinates for its holes position.

Figure 4.8: Evolution of the second-Piola Kirchhoff stress S11 of the SBR through thermo-
oxidative aging and mechanical loading. As it can be seen, higher aging durations caused higher
stress values due to hardening effect of thermo-oxidative aging.

76

25mmR=14mm115mm6mml=33mmR=25mml0=25.4mmA62 mm82.5 mmOC1C5C4C2C3OC1C2C3C4C5X (mm)Y (mm)001447.531.514.547.52321.240.55859Figure 4.9: Evolution of the second-Piola Kirchhoff stress S11 of the silicon adhesive through
hydrolytic aging and mechanical loading. As it can be seen, higher aging durations caused lower
stress values due to softening effect of hydrolytic aging.

a)

b)

Figure 4.10: Evolution of the chain density ρr of (a) SBR during thermo-oxidative and (b) silicon
adhesive during hydrolytic aging. Based on stress evolution, we expected that chain density in-
crease during thermo-oxidation due to cross-linking phenomenon and the opposite for hydrolysis
which both can be seen here.

77

dimensional network alteration theory was adapted to understand and predict the inelastic behavior

of elastomers during hydrolytic and thermo-oxidative aging. Despite significant efforts in model-

ing the aging and its influence on the mechanical performance of cross-linked matrix materials,

there have been very few models that can reliably predict the aging with reasonable computational

efforts. Assuming the alteration of the chain density along the aging trajectory is identical to the

peroxide cross-link density for thermo-oxidation, and the change of the average molecular weight

for hydrolysis, the strain energy of polymer matrix can be modified to predict chemical aging. The

proposed model includes only six physically inspired material parameters. Thus, while it is com-

putationally efficient, it shows good agreement with the experimental data. The proposed model

breaks the polymer matrix corresponding to four external load scenarios of (i) deformation, (ii)

deformation history, (iii) storage time, and (iv) aging temperature. In this respect, in the course of

deformation-induced damage, the continuous loading can detach chains from either cross-linked or

aggregate surfaces and reform the composition of the polymer matrix with longer chains. However,

during the immersion of polymer matrix in the water and the absence of the induced deformation,

the aging can reduce the polymer active agents, which are involved in the rubber elasticity. In

contrast, in thermal degradation, the aging mechanism is vice versa and increases polymer chains’

density. The excellent performance of the proposed method was proven by validating against a

new set of own experimental data, particularly designed to reveal the evolution of constitutive be-

havior during hydrolytic and thermo-oxidative aging. Moreover, the proposed model implemented

for FE simulation and the tangent moduli for large deformation regimes is formulated in terms of

internal variables. The model also suffers from some limitations such as inability to predict in-

homogeneous aging, which could be addressed in the near future. Here, our formulation neglects

the rate dependency so-called visco-elastic behavior of the polymer matrix, while the primary focus

of this work is on the quasi-static behavior. With respect to its computational efficiency, simplicity,

accuracy, and interpret-ability, the model is the right choice for advanced implementations in FE

programs.

78

CHAPTER 5

A Multi Physics Approach on Modeling of Hygrothermal Aging and its

Effects on Constitutive Behavior of Cross linked Polymers

5.1

Introduction

Due to the unique features such as high stretchability, chemical toleration, wear resistance, anti-

slipping, and durability; industrial applications of the cross-linked polymers have been growing

through past decades [93]. Their practical function can cover a broad range of applications from

insensitive to sensitive industries such as automotive, medical, mining, petroleum, and energy

storage [145, 146]. In most such cases, the polymeric components are simultaneously exposed

to mechanical loads as well as environmental elements like humidity and oxygen. This affair can

considerably reduce the nominal service life of the polymers and cause irreversible damage to their

functionality [96, 147, 148]. In the literature, this phenomenon is generally regarded as the aging

mechanism, which is classified into different subcategories based on the chemical source of aging.

Thermo-oxidation [101, 102, 103, 149], hydrolysis [99, 58], hygrothermal [97, 98] and photo-

oxidation [104, 105, 106] are some of the main types of aging in which oxygen, water, and UV

play major roles in polymer degradation, respectively (see Fig. 5.1). Such classes can be either the

result of individual aging elements or due to combined synergized effects of the above-mentioned

environmental factors [5]. For instance, hygrothermal aging is one of the combined mechanisms in

which oxygen and water molecules both are the leading actors in polymer deterioration. Generally,

since these factors coexist in nature, natural aging is seen as a combined effect of all conditions

together. Lately, solo environmental damage models have been of significant interest in the literary

and research communities. However, there are still wide gaps in the modeling of coupled-aging

phenomena since finding the mutual effect and correlation of each environmental element is chal-

79

lenging and needs a huge set of experimental observations. Accordingly, realizing the correlation

of the mechanical and coupled-environmental elements (damage accumulation) is crucial for pro-

viding a realistic prediction of such materials’ behaviour. There currently exists no constitutive

model of polymeric materials that can simultaneously consider two damage mechanisms. There-

fore, in this paper, we mainly focus on providing a multi-Physic model to predict the constitutive

response of cross-linked polymers along with hygrothermal aging, which is a coupled concurrent

effect of thermo-oxidative and hydrolysis aging.

So far, we have studied the individual effect of hydrolytic [9, 132] and thermo-oxidative aging

[11] and the contribution of their decay functions on the constitutive behavior of cross-linked

polymers in our recent works. Following the fundamental formulations presented in those studies,

our goal here is to provide a newly physically-based model to predict the constitutive behavior

of cross-linked polymers exposed to hygrothermal aging which is assumed as the combination

of thermal and hydrolytic aging. Due to the synergistic effects of combined hydrolytic, thermo-

oxidative, and hygrothermal aging, the devised model should be able to capture the constitutive

behavior of polymers in all three mentioned aging scenarios simultaneously.

Moreover, our model also should be able to predict the inelastic properties of cross-linked

polymers as well, in particular, the idealized Mullins effect (stress softening after the first loading

cycle) and permanent set (time-independent residual strain) which are necessary to be considered

for reliable prediction of material response. Many of these inelastic features have been modeled

so far in the context of micro-mechanical and phenomenological approaches [150? , 151]. While

micro-mechanical approaches are based on the physically inspired models developed from statis-

tical mechanics of polymer structure and physics of polymer chain network [152, 153, 154, 155],

the phenomenological ones use a mathematical model, usually without physical meaning [156].

Subsequently, a growing interest in micro-mechanical models is developing recently. For instance,

Guo and Zaïri [157] developed a micro-mechanical model to predict failure in elastomers while

Ben Hassine et al. combined the fracture mechanics and the intrinsic defect concept to establish

a failure criterion accounting for elastomers during thermal aging [158]. Accordingly, we chose

80

Figure 5.1: Different categories of aging.

micro-mechanical approaches for this study.

Aging is a complex continuous phenomenon that involves several chemo-mechanical aspects

of the material [? ]. While testing in the real-time aging trajectory can be expensive and chal-

lenging, the experimental investigation of this phenomenon is among key interests of the mechan-

ical engineering community. Therefore, several experimental studies were conducted to investi-

gate the behavior of the polymeric components along with the solo chemical aging mechanisms

[159, 160, 161]. In this respect, some of them used the concept of accelerated aging to keep the

practical relevance of experiments alive besides reducing the experimental costs [162]. Within this

concept, the elevated temperature has been utilized to reduce the actual aging time and to perform

the test at a reasonable cost. Accordingly, Celina’s [107, 108] literature review on this concept

shed light on its various related aspects. Although this approach significantly cuts down the test

81

Polymer materialOther sensitivitiesBacteria (Biodegradation),Metal contaminants, etc Creep,HTemperatureHydrolysisOHOxygenTemperatureThermo-Oxidative agingHygrothermal agingOxygenTemperatureMoistureOxygenTemperatureUVPhoto-Oxidative agingexpenses, it is not always the best scenario for the investigation of aging phenomena [109, 110]. In

fact, it can impose some additional chemical anomalies, such as diffusion-limited oxygen (DLO)

on the material. DLO, generally, causes heterogeneous aging and occurs when the rate of input

oxygen is lower than oxygen consumption in the polymer matrix [163]. These potential problems

are widely investigated by several researchers [111, 112]. The accelerated aging method can also

be conducted using various mediums such as air, distilled water, and salinated water [? 164]. Ac-

cordingly, Rabanizada et al. [114] investigated the effects of different aging media on the natural

rubber (NR). They performed DMA tests to investigate the change in the viscoelastic behavior of

NR due to aging. Moreover, thermo-oxidative aging has effects on viscoelastic relaxation and ther-

mal expansion of elastomers, which is shown by Shaw et al. [113]. Similarly, Pazur et al. [115]

studied the bulk properties of rubber during thermal degradation.

In comparison to thermal degradation, the experimental efforts for hydrolytic aging are not

as abounding. However, there are still some extensive studies that exclusively investigate the

effects of water aging on the elastomers. In this respect, Slater et al. [56] conducted experimental

research to find the effects of different conditions (i.e.

thermal, mechanical, and humidity) on

the mechanical performance of thermoplastic polyurethanes. They observed that the presence of

water causes plasticizing of the polymer matrix and results in a greater amount of compression set.

Ensuing similar trend, Farrar and Gillson [57] performed a series of experiments on polyglyconate-

B to find the relation between polymer morphology and degradation rate. Later on, Vieira et al.

[116] showed that during hydrolysis, the alteration of molecular weight average has a similar trend

with the change of tensile strength.

From the modeling side, so far, aside from our recent models [9, 11, 12], both phenomeno-

logical and micro-mechanical models are presented to describe the constitutive behavior of the

polymer matrix for solo chemical aging. In the case of thermo-oxidative aging, Ha-Anh and Vu-

Khanh [120] established a phenomenological model based on the Arrhenius function to predict

the behavior of polychloroprene rubber during thermo-oxidative aging. Moreover, Johlitz et al.

[24] presented a physically-based model to simulate the mechanical response of automobile bear-

82

ings during thermo-oxidative aging. In addition, Dippel et al. [123] introduced a material model

that describes material aging as a function of position from the surface to the center of the joints.

Recently, Konica and Sain [124] established a chemo-physical model that can predict the hetero-

geneous oxidation profile in the polymeric material due to high temperature. However, only a

few studies were conducted on the modeling of hydrolytic aging. In this regard, Vieira et al.[53]

adapted the Bergström and Boyce model in a way that it can predict the mechanical response of

biodegradable polymers during the aging trajectory. Later on, they presented a four-dimensional

computation approach to modeling biodegradable fibrous material during hydrolysis [26]. Besides,

the aging phenomenon can also be modeled using molecular dynamics simulations, although it re-

quires higher computational costs, implementation of either coarse-graining schemes [127? ] or

multi-scaling approaches [144] could accelerate these simulations.

The main objective of this work is to provide a physically-inspired constitutive model to predict

the inelastic response of cross-linked polymers during hygrothermal aging. The proposed model is

developed based on the concepts of network evolution [85] and network decomposition [17]. The

devised model takes into account the damage accumulation of coupled-aging phenomena by us-

ing a generic description of two competing solo-aging phenomena, namely thermo-oxidative, and

hydrolytic aging. The model is based on the assumption that each of the aging phenomena works

independently, and as a result, they can superpose each other. The model includes relatively fewer

material parameters and is validated with respect to new sets of our experimental data. The paper

is outlined as follows. The experimental study is first discussed in section 5.2. Then, the devised

constitutive model is presented in section 5.3, while the sensitivity analysis of model parameters is

discussed in section 5.4. Next, the parameter identification and model verification are presented in

section 6.4 and finally, the conclusion is provided in section 5.6.

83

5.2 Experimental Characterization

5.2.1 Material

In order to investigate the behavior of cross-linked polymers during hygrothermal aging and demon-

strate the generality of the devised model, the behavior of two different elastomers was experimen-

tally investigated. The neoprene rubber and Ethylene Propylene Diene Monomer (EPDM) rubber

in this study were procured as rectangular sheets of nearly 0.0625” thickness from one supplier.

From the same batch, each sample was cut with a punch die in the form of the standardized dumb-

bell shape (Die C from the ASTM-D412). These two rubbers were selected due to the fact that

their constitutive response was affected by hygrothermal aging. In fact, the material was chosen in

a way that all involving environmental factors (i.e. oxygen, temperature, and humidity) were able

to degrade them at the same time.

5.2.2 Hydrolytic aging

To simulate the hydrolytic aging condition, samples were fully immersed in sealed containers

filled with distilled water ( pH = 7 ) under constant pressure. Grit guards were placed under the

specimens to ensure that maximum surface area of samples was exposed to water (see Fig. 5.2a).

Elevated temperatures of 80◦ C, and 95◦ C were utilized to simulate accelerated aging conditions.

After certain amounts of storage times (i.e. 10, 20, and 30 days), samples were removed from the

containers and dried at room condition (i.e. 22 ± 2◦C, 50 ± 3% RH).

5.2.3 Hygrothermal aging

To achieve hygrothermal aging, saturated saltwater solutions were utilized to maintain particular

values of relative humidities inside sealed containers (see Fig. 5.2b). Salt solutions were replaced

on a weekly basis to avoid drying-up and to maintain constant relative humidity conditions inside

the containers. Two different humidities of 0%, and 80%, two temperatures 80◦C, and 95◦C and

84

three aging duration 10, 20, and 30 days were adopted for this aging experiment.

In fact, the

temperature setting for this experiment is kept in the range of 80◦C and 95◦C because further

increase in the temperatures would make the materials susceptible to diffusion-limited oxidation

(DLO) while lower temperatures would significantly slow down the kinetics of accelerated aging

of materials under study [165, 166]. Potassium chloride KCl is utilized as a chemically pure salt to

obtain the relative humidity of 80%. Whereas, for 0% RH (i.e. thermo-oxidative aging), molecular

sieves were used in containers.

(a)

(b)

Figure 5.2: Schematic sketches of the storage conditions of specimens in a) hydrolytic, and b)
hygrothermal aging.

5.2.4 Mechanical test

Mechanical tests were performed by a uni-axial universal testing machine (TestResources 311

Series Frame). Tests were done in a displacement control manner with a very slow strain rate

(i.e. 43.37 %

min ) to assure the quasi-static condition. The distance between the extensometer grips
was set to 25.4mm and all the experiments were performed at room condition (i.e. 22 ± 2◦C,

50 ± 3% RH). Each test has been repeated with 3 samples for reliability control. The error bar

is added to the experimental curves to show the variation of tests. Accordingly, sample reliability

was also not increased since the variation of mechanical tests is less than 5% ( see Fig 5.3 ). In

monotonic cyclic test, the samples were stretched to preset amplitudes of 1.25, 1.5, 1.75, and 2.

At each stretch amplitude, four cycles were applied. While in the failure tests, the samples were

subjected to an increasing monotonic tension load until they break (see Figs. 5.3 and 5.4 ). In the

course of deformation, elongation of the central zone of 25.4mm has been measured by an external

85

SpecimensDesalinated waterGrit guard Salt SolutionSpecimensGrit guardextensometer. The cyclic experiment is designed to illustrate the evolution of the permanent set

and stress softening during the primary loading for both virgin and aged samples.

a)

c)

b)

d)

Figure 5.3: The constitutive response of rubber specimens aged at multiple conditions; (a) Neo-
prene at T= 80◦C, RH= 80%, (b) Neoprene at T= 95◦C, (c) EPDM at T= 80◦C, RH= 80%, and
(d) EPDM at T= 95◦C, RH= 80%.

5.2.5 Water uptake

To measure water-uptake (see Fig. 5.5), unaged specimens were weighed on an electronic balance

(0.01mg accuracy) before being exposed to aging conditions. Specimens were then placed in

designated humidity and temperature conditions for planned aging time spans. After aging, the

specimens were carefully dried with a tissue paper and weights were recorded once again. The

water content Wc in the sample was measured as the % weight increase in the sample at completion

86

1234560481216Stress (MPa)Stretch (c) Unaged RH=80%, T=80°C, t=20d RH=80%, T=80°C, t=30d1234560481216Stress (MPa)Stretch (c) Unaged RH=80%, T=95°C, t=20d RH=80%, T=95°C, t=30d123450481216Stress (MPa)Stretch (c) Unaged RH=80%, T=80°C, t=20d RH=80%, T=80°C, t=30d123450481216Stress (MPa)Stretch (c) Unaged RH=80%, T=95°C, t=20d RH=80%, T=95°C, t=30da)

b)

Figure 5.4: Comparison of tensile stress at 100% strain influenced by hygrothermal aging versus
hydrolytic and thermo-oxidative aging at temperatures of 80◦C a) Neoprene, and b) EPDM.

of aging. The water uptake percent was calculated as

Wc(%) =

(cid:21)

(cid:20)Wt − W0
W0

× 100,

(5.1)

where Wt is the weight of aged sample at time t and W0 is the weight of unaged sample.

a)

b)

Figure 5.5: Comparison of water uptake specimens aged in submerged condition and in relative
humidity of 80%; (a) Neoprene, (b) EPDM.

87

01020303456Stress at 100% Strain (MPa)Time (days) RH=0%, T=80°C RH=80%, T=80°C Submerged, T=80°C01020303456Stress at 100% Strain (MPa)Time (days) RH=0%, T=80°C RH=80%, T=80°C Submerged, T=80°C010203002468Water uptake (%)Time (days) Submerged, T=80°C RH=80%, T=80°C010203002468Water uptake (%)Time (days) Submerged, T=80°C RH=80%, T=80°C5.2.6 Cross-link density

The cross-link density analysis for Neoprene and EPDM rubber was done using equilibrium swelling

method [167]. Both unaged and aged specimens were immersed in toluene (molar volume 106.3 mL
mol )

and the specimen weight changes were recorded till a constant weight signifying maximum sol-

vent absorbance was reached. The cross-link density and average molecular weights between the

chains can be calculated [168] as

ne
V0

=

ρ
Mx

(1 −

2Mx
M

),

(5.2)

where ne is the molar cross-linking, V0 Is the specimen volume before swelling, ρ is the density

of the solvent, Mx is solvent molar mass and M is the molar volume of the polymer sample.

It is observed that cross-linking agents in Neoprene rubber are relatively more active and sus-

ceptible to temperature changes as compared to EPDM ( see Fig. 5.6). Especially, this phenomenon

is highlighted in environments favorable for oxidation reactions such as thermo-oxidative and hy-

grothermal ( 80%RH) aging conditions. But regardless of any behavioral variations in both mate-

rials, they responded to the aging environment in a similar fashion.

5.2.7 The competitive aging scenario

The results of cross-link density changes of samples aged in different aging environments ( see

Fig. 5.6), water uptake ( see Fig. 5.5), and mechanical tests ( see Fig. 5.4 ) show that the ef-

fects of hygrothermal aging are always between hydrolytic and thermo-oxidative aging. Indeed,

thermo-oxidative aging is induced by two competing sub-mechanisms of cross-link formation Cf

( see Fig. 5.6), and chain scission Sc which results in the variation of cross-link density CD as

δCDT ∝ C T

f − ST

c . Hydrolytic aging is mainly induced by chain scission which results in the

reduction of cross-link density δCDHd ∝ −SHd

c

and massive deactivation of polymer chains.

In thermo-oxidative aging, for most materials cross-link formation is prevalent thus C T

f > ST
c .

Hence, hygrothermal aging has two effects: i) addition of C T

f , and ii) reduction of SHd

c

due to

88

a)

c)

b)

d)

Figure 5.6: The cross-link density of rubber specimens during thermo-oxidation, hydrolytic and
hygrothermal environments aged at multiple conditions; (a) EPDM at temperature 80◦C, (b) EPDM
at temperature 95◦C, (c) Neoprene at temperature 80◦C, and (d) Neoprene at temperature 95◦C.

the reduced water uptake in hygro in comparison to hydro (see Fig. 5.5). Accordingly, one can

conclude that hygrothermal aging is induced by competitive aging of two concurrent phenomena

namely, hydrolytic and thermo-oxidative aging.

5.3 Constitutive Model

Mechanical loads and environmental conditions play critical role in the damage accumulation

within the polymer matrix. Damages induced by environmental loads are often irreversible due

to the chemical reaction of the polymer components with environmental factors such as oxygen

and moisture. Damages induced by mechanical loads are composed of both permanent and re-

89

01020300.00.10.20.3Cross-Link Density (mol/cm3)Time (days) RH=0%, T=80(cid:176)C  RH=80%, T=80(cid:176)C Submerged, T=80(cid:176)C01020300.00.10.20.3Cross-Link Density (mol/cm3)Time (days) RH=0%, T=95(cid:176)C  RH=80%, T=95(cid:176)C Submerged, T=95(cid:176)C01020300.00.10.20.3Cross-Link Density (mol/cm3)Time (days) RH=0%, T=80(cid:176)C  RH=80%, T=80(cid:176)C Submerged, T=80(cid:176)C01020300.00.10.20.3Cross-Link Density (mol/cm3)Time (days) RH=0%, T=95(cid:176)C  RH=80%, T=95(cid:176)C Submerged, T=95(cid:176)Cversible damages. Permanent damages are often induced by chemical changes in the network,

such as detachment of cross-links, and chain scission, while reversible damages are often induced

by changes in matrix topology like chain slippage, detachment of physical bonds, or chain ad-

sorption/detachment from particles surfaces. Hygrothermal aging is a kind of environmental dam-

age that results from two simultaneous sub-aging mechanisms, namely, thermo-oxidative and hy-

drolytic aging (see Fig.5.7). We consider thermo-oxidative aging to be governed by two competing

micro-mechanisms [137]: chain scission and cross-link formation (see Fig. 5.9), while hydrolytic

aging is considered as the result of parallel reduction in the active chains and cross-links [58] (see

Fig. 5.10). Accordingly, one should separately formulate the effect of each aging mechanism and

couple them to develop a combined aging constitutive model.

Hereafter, the following notations are employed: for a scalar X, a vector X , a second-order

tensor X, and a fourth-order tensor X. The subscript X• depicts the network which the parameter
X is associated with and the bar sign over the parameters ¯X = ¯X

l refers to their normalized value

with respect to the Kuhn length l.

5.3.1 Damage Accumulation

Combining hydrolytic and thermo-oxidative in the form of hygrothermal aging leads to strain en-

ergy of the polymer matrix as a function of three parallel damage sub-mechanisms, C T

f , ST

c , and

SHd
c

. To capture the mutual effects and the correlation of thermo-oxidative and hydrolytic aging,

we assume that each of the aging phenomena works individually, and consequently, they super-

pose each other due to their independent relation. Accordingly, for each damage sub-mechanism,

a separate degradation rate has been considered, which eventually defines the contribution of each

mechanism in the overall damage imposed on the matrix.

To characterize the effects of aging time on the constitute behavior of polymer matrix, we

defined two end-states of the material as the state of polymer matrix at initial state Ψ0 and fully

aged state at time infinity Ψ∞. The end-states are considered as extreme points of aging, and

accordingly, strain energy of the material in all other states of aging can be calculated through

90

Figure 5.7: Evolution of the sub-aging networks throughout hygrothermal aging. At unaged state
(i.e. t = 0), there is only unchanged network. As time passes, the oxygen and water molecules
attack the polymer matrix and lead to two sub-aging networks namely, thermo-oxidative and hy-
drolytic networks. The unchanged network continues its shrinkage till it disappeared and got
replaced by two new networks with different morphology at the fully aged state (i.e. t = ∞).

predefined shape functions, N (t, T ) and N ′(t, T ), as

ΨM (t, T, RH, F) = N (t, T )Ψ0 + N ′(t, T )Ψ∞;

N ′(t, T ) = 1 − N (t, T ),

(5.3)

where N (t, T ) is the shape function that represents the contribution of each initial state in the

current state. In most solo aging of the polymeric systems, the evolution of the shape function can

91

t = t0t = t1t = tt = t2Unaged stateFully aged stateUnchanged networkHydrolytic aging networkThermo-oxidative aging networkΨ0htΨΨhΨtΨΨ0Figure 5.8: Decomposition of the hygrothermal aging into two aging networks namely, thermo-
oxidative and hydrolytic networks. The hydrolytic network, itself, decomposes to the two morphed
and deactivated networks.

be well captured by the Arrhenius decay function [73] as

N (t, T ) = exp

−γ exp(−

(cid:18)

(cid:19)

)t

.

Ea
RT

(5.4)

Here, γ is the underlying degradation rate, Ea the activation energy, and R = 8.314[J]/[mol][K]

the universal gas constant. The main challenge is to define the status of the polymer matrix at the

fully aged state Ψ∞. Here, based on the theory of network decomposition, we assume that fully

aged polymer matrix decomposes into two different networks (see Fig. 5.7); one of them resulted

from the diffusion of oxygen to the polymer matrix Ψ∞

t , and the second one modeled the effect of

water on the polymer matrix Ψ∞

h . Hence, the infinity state of the polymer matrix can be written as

Ψ∞ = (1 − β)Ψ∞

t + βΨ∞
h ,

(0 ≤ β ≤ 1),

(5.5)

where β(t, T, RH) defines the portion of each network, which contributes to the complete aging

of polymer matrix for a certain amount of aging temperature and humidity (see Fig. 5.8). The

92

N(t,T)Ψ0β(t,T,RH)αmΨdΨhΨtΨFigure 5.9: A schematic diagram of thermo-oxidative aging. Thermal aging results from two
simultaneous phenomena; cross-linking, and chain scission and it can be best modeled via a new
network with shorter chains [12].

hydrolysis network Ψ∞

h , itself, decomposes into two independent sub-networks namely morphed

Ψ∞

m and deactivated networks Ψ∞

d (see [9] for details) and can be written as

Ψ∞

h = αΨ∞

m + (1 − α)Ψ∞
d ,

(0 ≤ α ≤ 1),

(5.6)

where the parameter α depicts the contribution of each sub-microstructural phenomenon and is

considered as a fitting parameter. The next step is to find the portion of each aging mechanism (i.e.

the parameter β ), which is involved in the construction of the polymer elasticity. In this respect,

experimental studies reveal that the more water in the polymer matrix, the lower Tg it has and, as a

result, the more plasticizing damage [169]. Hence, the parameter β should be a function of water

absorption. Here, the main assumption of the model is that the contribution of the water network

93

t = tnFillerCross-linkActive siteUnattacked networkThermo-oxidative networkCross-linkingScission<ProbabilityThermal agingt= t0nProbabilityR0R0RtP0P0Ptin the total entropy of polymer matrix is linearly related to the amount of water in the matrix. Base

on this assumption, the decay function can be assumed similar to the water uptake function. It

is clear that the rate of damage in the polymer matrix is different from the rate of water uptake

in the samples since the water first penetrates into the samples, then it starts imparting attrition

to the polymer matrix. Within this concept, experimental observations confirm that [? 170] the

water uptake is a nonlinear function of time t, temperature T , and humidity RH and could be best

modeled via Arrhenius function [169]. Hence one can define the function β(t, T, RH) as

β(t, T, RH) = erf

(cid:18) RH
Q

(cid:19)(cid:115)

(cid:20)

θ exp

−

(cid:21)

t,

Eb
RT

(5.7)

where the parameter Eb is the activation energy which is attributed to the water diffusion, and

parameters Q and θ are adjusting parameters that keep the range of β(t, T, RH) between 0, and 1

(i.e. 0 ≤ β ≤ 1). The lower bound of the β (i.e. β = 0) indicates that there is no humidity effect on

mechanical response and only thermal degradation Ψb causes the matrix deterioration. In contrast,

the upper bound (i.e. β = 1) shows that the hydrolytic aging wins the competition between the two

sub-aging phenomena, and the effects of oxygen are negligible regarding polymer degradation.

5.3.2 3D to 1D transition

Based on the experimental observation, it is assumed that the micro-structure of the polymer net-

work is embodied by the isotropic spatial distribution of polymer chains, which are scattered in all

directions (see Fig. 6.6a). Accordingly, the macroscopic energy of polymer network ΨM can be

carried out by integrating overall available configurations at micro-scale [171, 172]. This process

can be done by numerical integration scheme by utilizing the concept of the unit micro-sphere (see

Fig. 6.6b, and see [85] for details). Hence, the isochoric strain energy of the polymer network can

94

Figure 5.10: A schematic diagram of hydrolytic aging. Hydrolytic aging is due to the reduction of
cross-links and chain scission and it can be modeled with a new network with longer chains [9].

be written as

ΨM =

1
As

(cid:90)

S

d
ψM d

d

A ∼=

d i
ψM wi ⇒

k
(cid:88)

i=1

d
ψM = N (t, T )

d
ψ0 + (1 − β(t, T, RH))N ′(t, T )

d
ψt + αβ(t, T, RH)N ′(t, T )

d
ψm,

(5.8)

where the parameter As is the surface area of unit sphere, Ad the unit area normal to the direc-

tion d , wi Gaussian weight coefficients according to the direction d i (i = 1, 2, . . . , 45) , and ψ•

represents the total energy of each network in direction d .

Chain length distribution function Each network is made by many sub-networks, each having

the same set of chains with a different relative length (number of segments n) in an arbitrary

direction d . The polymer chains in each network are bonded to cross-link only at their first and

95

Hydrolytic AgingScissionReduction of cross-linkst = tProbabilitynFillerCross-linkActive siteUnattacked networkMorphed networkDeactivated networkt = t0nProbabilityRm>R0R0P0P0Pma)

b)

Figure 5.11: Schematic figures of a) a freely rotating chain consisting of n monomers with Kuhn
length of l, and b) the unit micro-sphere representation with the unit basis vectors of e 1, e 2, and
e 3 and normal direction d = d1e 1 + d2e 2 + d3e 3.

last segments, giving specific end-to-end distances of ¯R•. Here, the chains’ distribution function
at a given distance ¯R• is considered as normal distribution and can be written as

P• (n) =

1
√
πσ2
2

exp(

(n − µ•)2
−2σ2

),

(5.9)

where, σ is the standard deviation, and µ• the mean value. In order to summarize our framework,

the damage accumulation is modeled according to the following functions:

• N (t, T ) is a temporal shape function that is mainly influenced by t.

• β(t, T, RH) is a nonlinear function of the RH, t, and T .

• α is a compound based function and that is why it is constant for each material.

• Ψ0 is the strain energy of the polymer matrix in the virgin state.

• Ψ∞
t

is the strain energy of the thermo-oxidative network in the infinity state. This network

is constructed with higher cross-link density and shorter chain (see [12] for details) which

can be best modeled by changing the mean value of chain length distribution and end-to-end

distance of chain (i.e. ¯Rt < ¯Ro, and µt < µo).

96

lrL = nld• Ψ∞

m represent the stain energy of the morphed network (see [9] for details) which consists of

longer chains with respect to the virgin network (i.e. ¯Rm > ¯Ro, and µm > µo).

Energy of a single chain The strain energy density of a freely jointed chain (FJC) were obtained

by utilizing an enhanced version of non-Gaussian statistics (see [142] for details) and can be written

as

ˆψc (n, ¯r•) = nKbT

(cid:90) Ω

0

ˆβdΩ,

ˆβ =

(cid:20)

1 −

1 + Ω2
n

(cid:21)

β,

(5.10)

L = ¯r•

where Ω = r•

n ) of a single polymer chain, while ¯r• is end-to-end
distance in the deformed configuration. The parameter β = L−1(Ω) denotes the inverse Langevin

n is the stretchability ratio ( ¯r•

function (ILF), which can not be explicitly derived and should be implicitly determined. There are

several accurate rational functions that are developed in the literature [173] in which choosing the

suitable approximation function generally depends on the maximum stretchability of the polymer

chain. Here, a simple yet accurate (maximum error less than 1%) first-order rational function is

being utilized as ILF [84], and can be written as

L−1 (x) ∼=

x
1 − x

+ 2x −

8
9

x2.

(5.11)

Polymer chain detachment

In all networks, permanent damage results from the debonding of

shorter chains from cross-links or aggregates’ surface during the primary loading. Detachment

or sliding occurs if the chain force fp goes beyond the finite breakage force of a polymer chain

denoted by fc. Obviously, this force approaches infinity fp → ∞ as the end to end distance of

the chain ¯r• becomes close to its contour length ¯r• → n. This detachment starts with the shortest

chain and gradually affects the longer chains. Considering the parameter ν as the sliding ratio

of a polymer chain (see [85] for details), the length of the shortest available chain in a deformed

network is

n•min = νλd

max

¯R•,

(5.12)

97

while λd

max is the maximal micro-stretch previously reached in the loading history in an arbitrary

direction d . On the other hand, the length of longest available chain in each network n•max is

considered as a cut-off length of probability distribution function of polymer chains P• (n). Ac-

cordingly, the set of available relative chain lengths in the direction d of networks are given by

(cid:18)d

λmax

(cid:19)

(cid:26)

=

n

D•

(cid:19)

(cid:18)d

λmax

n•min

(cid:12)
(cid:12)
(cid:12)
(cid:12)

(cid:27)

≤ n ≤ n•max

,

where

n•max = µ• + 5σ.

(5.13)

5.3.3 Network rearrangement

The network rearrangement is categorized into two independent groups namely mechanical-induced

and environmental-induced rearrangement.

Mechanical-induced rearrangement

. The main assumption here is that the detachment does

not lead to the complete energy dissipation of the chains since the chains will still serve as part

of a larger macro-molecule. As such, the detachment will lead to partial dissipation of the energy

only, and part of it will remain as kinetic energy in the new, longer chains that are formed during

this rearrangement process (see [85] for details). All in all, the number of active segments remain

constant during deformation and given by the following expression

(cid:90)

D•(1)

N0•P• (n) ndn =

(cid:90)

D•(λmax)

˜N0•(λmax)P• (n) ndn.

(5.14)

The debonding process is irreversible; no new bonds will appear because of loading and unloading

the specimen. For this reason, the number of active chains can be defined as

˜N0•(λmax) = Φ• (λmax) N0•,

where

Φ• (λmax) =

(cid:82)

D•(1)
(cid:82)

P• (n) ndn

P• (n) ndn

.

(5.15)

D•(λmax)

Environmental-induced rearrangement As the probability of the chain length changes over

time, conservation of active segment numbers should be considered in the number of chains in

98

each network (see [9, 12] for details). Accordingly, by keeping the number of segments constant

at any time between zero and infinity, we have

N0• = ϕ•N0,

where

(cid:82)

P0 (n) ndn

ϕ• =

D•(1)
(cid:82)

D•(∞)

P• (n) ndn

.

(5.16)

Evolution of sub-network energy Average strain energy of an arbitrary sub-network in direction

d , namely ψd

• can be written with respect to chain’s length distribution in that direction. Thus,

considering Eqs. 5.10, 5.15, and 5.16, ψd

• is given by

d
ψ• = N0

(cid:90)

(cid:19)

(cid:18) d

λmax

ϕ•Φ•

ˆψc (n, ¯r•) P• (n) dn

(5.17)

D•(λmax

d)

Accordingly, the strain energy for each network can be written as

d
ψ0 = N0

(cid:90)

(cid:19)

(cid:18) d

λmax

Φ0

ˆψc (n, ¯r0) P0 (n) dn

D0(λmax

d)

d
ψt = N0

(cid:90)

(cid:19)

(cid:18) d

λmax

ϕtΦt

ˆψc (n, ¯rt) Pt (n) dn

Dt(λmax

d)

d
ψm = N0

(cid:90)

(cid:19)

(cid:18) d

λmax

ϕmΦm

ˆψc (n, ¯rm) Pm (n) dn.

(5.18)

Dm(λmax

d)

5.3.4 Final formulation

The total strain energy of polymer matrix can be written in terms of deformation gradient F,

storing temperature T , time t, and relative humidity RH by substituting Eqs. 5.4, 5.5, 5.6, 5.7,

and 5.18 in Eq. 5.3. Hereafter, F stands for the deformation gradient while C represents the right

Cauchy–Green tensor. By assuming incompressibility condition (i.e. J 2 = det C = 1), the first

99

Piola-Kirchhoff stress tensor T can be obtained as

T =

∂ΨM
∂F

−pF−T = N (t, T )

(cid:20)

+N ′(t, T )

∂Ψ0
∂F

β(t, T, RH)

∂Ψ∞
t
∂F

+ (1 − β(t, T, RH))

(cid:21)

∂Ψ∞
h
∂F

−pF−T ,

(5.19)

where p is a scalar parameter to guaranty the incompressibility. In view of the Eq. 5.19, the

first Piola–Kirchhoff stress tensor T is a nonlinear tensor valued function of the state variable F.

By utilizing the chain rule, it can be obtained by as

∂Ψ0
∂F

=

k
(cid:88)

j=1

wj

∂

∂

∂d j

¯Cd j
∂ ¯F

∂ ¯F
∂F

,

:

1
d j
λ

2

d
ψ0
d j
λ
d
ψm
d j
λ

∂

∂
d
ψt
d j
λ

∂

∂

∂Ψ∞
h
∂F

= α

k
(cid:88)

j=1

wj

∂Ψ∞
t
∂F

=

k
(cid:88)

j=1

wj

∂d j

¯Cd j
∂ ¯F

∂ ¯F
∂F

,

:

1
d j
λ

2

∂d j

¯Cd j
∂ ¯F

∂ ¯F
∂F

,

:

1
d j
λ

2

(5.20)

The detailed derivation of the stress tensor is provided in appendix A.

As a summary, this model should be considered an extension of our recent models [9, 132, 11]

which predict the polymer constitutive behavior during dual aging mechanisms. Similar to our

recent works, network evolution theory [85] considered as a kernel of our model, which is a tool

to predict deformation-induced damages. The term deformation-induced damage is used to refer

to any types of damages which are induced by the deformation such as strain-induced anisotropy,

and idealized Mullins effect. This kernel could be replaced by simpler theories such as Analytical

network-averaging [174] to improve the computational efficiency of the model. However, the

reason that why we choose a numerical scheme method as the kernel of our model is that it be

able to model the inelastic directional-dependent phenomena such as the idealized Mullins effect.

In fact, it is a trade-off between the accuracy and simplicity of the model. In order to simplify

our model kernel, we assumed the homogeneous distribution of the micro-stretch λd within the

network. Accordingly, the relationship between the micro-stretch λd and macro-stretch χd in the

100

direction d is assumed as χd = λd . Here, the main novelty of our current model lies in the

extension in modeling of environmental damages (see Fig. 5.12). Our recent works can only

consider one type of single aging mechanism such as hydrolytic or thermo-oxidative aging. In this

paper, we offered a new generic model to predict the behavior of cross-linked polymers during

hygro-thermal aging which is a dual aging mechanism. Hence, hydrolytic and thermo-oxidative

aging considered as subsets of hygro-thermal aging since it is due to the synergized effects of

oxygen, water, and temperature. From the mathematical point of view, with one set of material

parameters, our generic model is able to cover these two types of single aging mechanisms as well

as their nonlinear interaction.

Figure 5.12: Schematic figure of different damage mechanisms’ modeling.

101

λ=1λ=λ1λ1λ=λ2>λ1λ2RRRDeformation induced damage<R0RtThermo-oxidativeRm>R0HydrolysisR0Virgin materialHygrothermal Nonlinear interactionEnvironmental damagesMechanical damagesMaterial Parameter Description

Table 5.1: Summary of material parameters.

N0
¯R0

µ0
σ
ν
¯Rm

µm

α

¯Rt

µt

γ
Q, θ
Ea, Eb

The number of chains of virgin network per unit volume
Normalized end-to-end distance of the chains in the unchanged
network
Mean value of unchanged network chain length distribution
Standard deviation of chain length distribution in all networks
Sliding ratio denoting the strength of bonds
The normalized end-to-end distance of chains in the morphed net-
work which is bigger than the unchanged one ( ¯Rm > ¯Ro)
The mean value of morphed network chain length distribution
which is bigger than the unchanged one (µm > µo)
The percentage of active chains in the hydrolytic network after
aging
The normalized end-to-end distance of chains in the thermo net-
work which is smaller than the unchanged one ( ¯Rt < ¯Ro)
The mean value of thermo network chain length distribution
which is smaller than the unchanged one (µt < µo)
Arrhenius rate factor
Adjusting parameters to keep β between zero and one
Activation energies

5.4 Parameter sensitivity analysis

The developed model has fifteen material parameters, five of which belong to the unchanged net-

work (N0, ¯R0, µ0, σ, and ν), three to the morphed ( ¯Rm, µm, and α), two to the thermo-oxidative
network ( ¯Rt, and µt), and the rest are for the decaying functions (γ, Q, θ, Ea, and Eb). The material

parameters for the unchanged, morphed, and thermo-oxidative networks along the aging trajectory

are comprehensively studied in our previous work (for detail see [9, 12]). As a summary, each of

these material parameters is described in Table 5.1.

In order to better understand the effect of the damage accumulation parameters (Q, θ and Eb)

on the inelastic response of the material, a parametric sensitivity analysis is utilized to capture the

variation of each of these parameters. The reference sets of material parameters are presented in

Table 5.2. To highlight the interaction of the deformation-induced damage and aging degradation,

the analysis of the parameters is carried out in two different aging times of 10, and 100 days.

102

Table 5.2: The model’s reference material parameters for parametric analysis.

N0kbT ¯R0
[M P a]
3

¯Rm

¯RT µ0 µm µT σ ν

α Q γ

[1/day]

θ
[1/day] [K]

Ea
R

Eb
R
[K]

5 7

2

10 15 5 4 1.01 0.1 100 5.42e+12 1e+9

1.17e+04 0.965e+04

In fact, the parameter Ea

R and γ are calibrated in a way that temporal shape function N ′(t, T ) is

approached to its upper limit at aging time 100 days and temperature 95◦C (i.e. N ′(t, T ) t=100d−−−−−→
T =95◦C

1

) to represent the infinity state of the material. Accordingly, the effect of material parameters

involved in β(t, T, RH) function on the inelastic response of the reference sample is investigated

and described below.

• Parameter Q, controls the effect of humidity on the response of the material. Indeed, for

a certain value of humidity, having the higher value of Q, means the contribution of the

hydrolytic network in the infinity state of the material is lower, and as a result thermo network

will be dominant ( see Fig. 5.13a).

• Parameter θ, the water uptake rate factor, controls the effect of aging time t on the contri-

bution of the hydrolytic network in the infinity state. Here, the aging time is involved in

N (t, T ), and B(t, T, RH) functions and as a result, it is a very complex variable. Indeed,

aging time t can decrease the contribution of the unchanged network as well as increase the

role of the hydro network in the response of the material. Hence, a higher value of the θ

results in a higher effect of water on the response of the material ( see Fig. 5.13b).

• Parameter Eb

R , the water uptake activation energy, control the effect of aging temperature T

on the contribution of the hydro network. Similar to the aging time t, aging temperature T is

playing role in both temporal and humidity functions. Accordingly, the higher value of Eb
R ,

the less contribution of water in the response of material ( see Fig. 5.13c).

103

Figure 5.13: Sensitivity analysis of aging trajectory material parameters in the proposed model.

104

a) Variation of QStretch χStretch χb) Variation of θStretch χStretch χt = 10 days, T = 95 C, RH = 80%t = 10 days, T = 95 C, RH = 80%t = 10 days, T = 95 C, RH = 80%t = 100 days, T = 95 C, RH = 80%t = 100 days, T = 95 C, RH = 80%t = 100 days, T = 95 C, RH = 80%c) Variation ofEb/R061218123Stress (MPa)Stretch χQ = 100Q = 1000Q = 10Stretch χQ = 100Q = 1000Q = 10061218123Stress (MPa)θ= 1e7θ= 1e9θ= 4e9061218123Stress (MPa)θ= 1e7θ= 1e9θ= 4e9061218123Stress (MPa)b/R = 9.15e3Eb/R = 9.65e3Eb/R = 1e4E061218123Stress (MPa)b/R = 9.15e3Eb/R = 9.65e3Eb/R = 1e4E5.5 Parameter identification and model verification

5.5.1 Parameter identification

To validate the proposed constitutive model, its responses were compared to our experimental

data. Depending on the specific aging phenomena to be modeled, different portions of the strain

energy function can be employed, as depicted in Fig. 5.14. The model is a kind of modular

Figure 5.14: Aging phenomena to be modeled and the corresponding network decompositions.

platform to describe damage accumulation in the course of hygrothermal corrosion. Each network

is responsible to model the energy reduction of sub-aging mechanisms,i.e. i) hydrolytic, and ii)

thermo-oxidative. The correlation between these two independent aging mechanisms is held by

a nonlinear function of time, temperature, and relative humidity. Therefore, each network should

be fitted to the experiments separately. Accordingly, by setting β to its upper limit (i.e. β = 1 ),

the thermo-oxidative network will be eliminated and therefore the three material parameters of the

morphed network, as well as the unchanged one, can be achieved by fitting to the following three

105

Unchanged networkHydrolytic networkThermo networkStrain energy functionMaterial parametersUnchanged networkHydrolytic networkThermo networkCoupling parameterscases;

• One loading-unloading cycle of the virgin sample to find the effects of λ and λm

• Only, one point of loading curve of 20 days hydrolytically aged sample at 80◦C (i.e. sub-

merged condition) to find the effects of t, and T

• The primary loading curve of 30 days hydrolytically aged sample at 95◦C (i.e. submerged

condition) to find the effects of λ, t, and T

On the other hand, the material parameters related to the thermo-oxidative network can be obtained

by adjusting B = 0 and fitting them to the following cases:

• The primary loading curve of 30 days thermally aged sample at 95◦C (i.e. RH = 0) to find

the effects of λ, t, and T

• Only one point of loading curve of 20 days thermally aged sample at 80◦C (i.e. RH = 0) to

find the effects of t, and T

Lastly, the parameters in function β(t, T, RH) can be obtained by fitting to the one point of

loading curve of 20 days hygrothermally aged sample at 80◦C (i.e. RH ̸= 0) to find the effects of

RH. Once all the material parameters are achieved, the model can predict other loading-unloading

curves for the two thermo-oxidative and hydrolytic aging and their collaborative damages.

5.5.2 Model verification

Figs. 5.15, 5.16, 5.17, and 5.18 shows the prediction of the devised model which were bench-

marked against our experimental results, while the Tables 6.2, and 5.4 depicts the so-obtained

material parameters.

In these figures, solid lines represent model predictions while dot points

depict the experimental data. The arrows demonstrate the curves and the points which are used for

fitting.

To assess the goodness of fit and robustness of the model, the coefficient of determination R2

is calculated in each aging temperature (See Figs. 5.15, 5.16, 5.17, 5.18). In this respect, the

106

a)

c)

b)

d)

Figure 5.15: Validation of the proposed model for neoprene rubber in the unaged state, thermo-
oxidative (i.e. RH = 0%), hygrothermal (i.e. RH = 80%), and hydrolytic aging at T = 80◦C , in
multiple aging conditions; (a) t = 20d, (b) t = 30d, (c) constitutive behavior of virgin sample, (d)
Comparison of one cycle load on samples aged at different relative humidities and aging times.

107

051015123456Stress (MPa)Stretch χExperimental DataUnagedRH = 0%RH = 80%SubmergedT = 80 C, t = 20dFitted pointsR2 = 0.91Fitted curve051015123456Stress (MPa)Stretch χT = 80 C, t = 30dUnagedRH = 0%RH = 80%SubmergedExperimental DataFitted pointR2  = 0.93012311.52Stress (MPa)Stretch χExperimental DataUnagedFitted pointR2  = 0.95048123Stress (MPa)Stretch χRH = 90%, t = 20dT = 80 CRH = 40%, t = 30dRH = 90%, t = 30dRH = 40%, t = 20da)

c)

b)

d)

Figure 5.16: Validation of the proposed model for neoprene rubber in the unaged state, thermo-
oxidative (i.e. RH = 0%), hygrothermal (i.e. RH = 80%), and hydrolytic aging at T = 95◦C , in
multiple aging conditions; (a) t = 20d, (b) t = 30d, (c) constitutive behavior of aged sample at t =
30d, RH = 0% , (d) Comparison of one cycle load on samples aged at different relative humidities
and aging times.

108

061218123456Stress (Mpa)Stretch χUnagedRH = 0%RH = 80%SubmergedExperimental DataT = 95 C, t = 20dR2= 0.904061218123456Stress (MPa)Stretch χUnagedRH = 0%RH = 80%SubmergedExperimental DataR2= 0.91T = 95 C, t = 30dFitted curves0481211.52Stress (MPa)Stretch χExperimental DataT = 95 C, t = 30d, RH = 0%R2= 0.91051015123Stress (MPa)Stretch χT = 95 CRH = 90%, t = 20dRH = 40%, t = 30dRH = 40%, t = 20dRH = 90%, t = 30da)

c)

b)

d)

Figure 5.17: Validation of the proposed model for EPDM rubber in the unaged state, thermo-
oxidative (i.e. RH = 0%), hygrothermal (i.e. RH = 80%), and hydrolytic aging at T = 80◦C , in
multiple aging conditions; (a) t = 20d, (b) t = 30d, (c) constitutive behavior of virgin sample, (d)
Comparison of one cycle load on samples aged at different relative humidities and aging times.

109

05101512345Stress (MPa)Stretch χUnagedRH = 0%RH = 80%SubmergedExperimental DataT = 80 C, t = 20dR2  = 0.91Fitted curveFitted points05101512345Stress (MPa)Stretch χT = 80 C, t = 30dUnagedRH = 0%RH = 80%SubmergedExperimental DataR2= 0.92Fitted point012311.52Stress (MPa)Stretch χUnagedExperimental DataFitted pointR2  = 0.950510123Stress (MPa)Stretch χT = 80 CRH = 90%, t = 20dRH = 40%, t = 30dRH = 90%, t = 30dRH = 40%, t = 20da)

c)

b)

d)

Figure 5.18: Validation of the proposed model for EPDM rubber in the unaged state, thermo-
oxidative (i.e. RH = 0%), hygrothermal (i.e. RH = 80%), and hydrolytic aging at T = 95◦C ,
in multiple aging conditions; (a) t = 20d, (b) t = 30d, (c) constitutive behavior of aged sample at
t = 30d and submerged condition, (d) Comparison of one cycle load on samples aged at different
relative humidities and aging times.

110

05101512345Stress (MPa)Stretch χUnagedRH = 0%RH = 80%SubmergedT = 95 C, t = 20dExperimental DataR2= 0.9105101512345Stress (MPa)Stretch χUnagedRH = 0%RH = 80%SubmergedExperimental DataR2= 0.93T = 95 C, t = 30dFitted curves04811.52Stress (MPa)Stretch χExperimental DataSubmerged, t = 30dT = 95 CR2 = 0.8970510123Stress (Mpa)Stretch χRH = 90%, t = 20dRH = 40%, t = 30dRH = 90%, t = 30dRH = 40%, t = 20dT = 95 CTable 5.3: Material parameters of the devised model for neoprene rubber.

¯Rm

N0kbT ¯R0
[M P a]
1.663 2.65 3.88 2.54 8.04 9.06 1.9 3.28 1.01 0.685 100 5.42e+12 3.50e+9 1.17e+04 0.965e+04

¯RT µ0 µm µT σ

θ
[1/day]

Eb
R
[K]

Ea
R
[K]

[1/day]

Q γ

α

ν

Table 5.4: Material parameters of the devised model for EPDM rubber.

¯Rm

N0kbT ¯R0
[M P a]
1.927 2.65 3.86 2.34 8.04 8.06 5.63 3.27 1.009 0.7639 80 5.42e+12 3.2e-1 1.17e+04 0.74e+03

¯RT µ0 µm µT σ

θ
[1/day] [K]

Eb
R
[K]

[1/day]

Q γ

Ea
R

α

ν

devised model is capable of predicting the constitutive behavior of aged samples With a marginal

error. The major source of error, here, is ignoring the effects of nonhomogenous aging. In fact,

hygrothermal aging ( i.e. water penetration) starts from the surface and is highly dependent on

diffusion factors that have not been considered in our recent model. The model is then used to

describe different cross-linked polymer damage behaviors as they age once it has been validated.

To this end, the changes in the loading-unloading response of the material aged at various aging

temperatures and humidity levels are compared (See Figs. 5.15d, 5.16d, 5.17d, 5.18d). Finally, it

should be noted that the model can be used for any loading scenario; all that needs to be done is to

change the deformation gradient F to match the desired loading. Accordingly, we mainly focused

on uni-axial tensile loading due to the ease of experiments.

5.6 Conclusion

We propose the first physically based model to account for multi-factor damage accumulation

of hygrothermal aging. Our model can take into account three damage mechanisms namely hy-

drolytic aging, thermo-oxidative aging, and deformation-induced damage at the same time. Each

of those damage mechanism is induced by sub-mechanisms that are extensively studied before but

111

has not been modeled in parallel so far. We modeled each sub-mechanism through a linear kinetic

equation which is then coupled into the modular network concept to allow calculation of synergies

and parallel damage accumulation through multiple sources. In this respect, hydrolytic aging is

considered to be induced by chain scission and reduction of the cross-link, while thermo-oxidative

aging is considered as a superposition of two competing effects of cross-link formation and chain

scission as previously discussed in our recent models [9, 11, 12]. Deformation-induced damage is

modeled by implementing the concept of network evolution [85]. Combining these three damage

mechanisms can define the damage accumulation in the case of hygrothermal aging. In order to

capture the mutual effects of thermo-oxidative and hydrolytic aging, the assumption has been made

that each of the aging phenomena can be superposed upon each other. From the micro-mechanical

perspective, hygrothermal damage can be regarded as a race among three micro-structural phe-

nomena; i) chain scission due to the presence of temperature, ii) reduction of cross-links attributed

to the attendance of water and iii) formation of cross-links as a consequence of oxygen reactions.

Utilizing the theory of network decomposition, all phenomena and their correlation were modeled.

The proposed model is micro-mechanically based and is mainly relevant on thin samples due to

our underlying assumption of homogeneous diffusion of oxygen and water throughout the matrix.

The model has been validated against extensive data-sets obtained from experiments we specif-

ically designed for concept validation. In view of its interpretation, precision, and deep insight,

it provides an insight into the nature of damage accumulation. The model is a good choice for

advance implementation in FE applications.

112

CHAPTER 6

Thermal Aging Coupled with Cyclic Fatigue in Cross-linked Polymers:

Constitutive Modeling & FE Implementation

6.1

Introduction

In view of polymeric materials’ prevalence in a myriad number of industries such as aerospace

where a component failure could potentially lead to disastrous outcomes; having a reliable model

that can predict their mechanical properties degradation induced by different aging conditions is a

necessary process [7, 11, 175, 9, 176, 177, 176]. The resistance of polymeric materials against ex-

ternal factors such as mechanical fatigue and environmental elements is the center of attention for

decades [178, 179, 180, 181, 182, 149, 145]. These external factors induce irreversible damages to

the polymer matrix, resulting in reduction of material’s nominal service life [183, 184, 147]. Ac-

cordingly, damage accumulation in the polymeric matrix is mostly due to two sources; mechanical

damages and environmental aging.

- Environmental aging: Environmental aging is an irreversible process that occurs due to

changes in the molecular structure of the polymer matrix [185, 186, 103, 157]. Environmental

aging can be induced either by a single environmental element such as oxygen [187, 8], water

[188, 189], UV [190, 191] or a synergized effect of multi-factor environmental elements such as

hygrothermal which is due to mutual effects of oxygen and water. Extensive efforts have been

made to describe the effects of environmental aging on the mechanical behaviour of polymers, and

many successful models were developed [114, 113]. However, when it comes to synergized effects

of chemical and physical aging [192], there is still a large gap in our understanding. Although the

popularity of the models to predict the behavior of materials during aging is increasing, aside from

our previous models [11, 7, 9], the theoretical efforts for that are limited. In this respect, Konica

113

and Sain [124] developed a 3D large deformation model for polymers that simultaneously consid-

ers the diffusion of chemical reaction to the material during high deformation regime. Besides,

they implemented their model in commercial finite-element software to numerically investigate

the coupled diffusion-reaction and constitutive response of polymers during oxidation aging. In a

similar fashion, Dippel et al. [193] provided a constitutive model to predict the behavior of ad-

hesives in the joints during chemical aging. Their model can consider the internal structure of an

adhesive by introducing an interface zone where the material properties alter from the surface to

the center of the joint. Moreover, in the same direction, Lamnii et al. [194] investigated the effects

of photo-oxidation on the fatigue life of low-density polyethylene. They used dissipated energy

in each cycle as a damage indicator and tried to model the evolution of damage with respect to

physical and chemical aging.

- Mechanical damages: Generally, mechanical damages retain the morphology of the polymer

matrix intact, and the decay is mainly attributed to the movement of polymer chains. Due to the

high cost of in-situ testing of elastomers, a reliable predicting tool becomes a dire need. In this

respect, mechanical damage models are mostly focused on fatigue analysis and there are very few

that describe the gradual loss of performance of the material in multi-cycle loading [195, 196].

In this regard, recently, we developed a 3D micro-mechanical model for elastomers Shakedown

where the samples are subjected to a high number of cyclic deformations [175]. Moreover, Wang

et al. [197] developed a phenomenological model to predict the shakedown of tough hydrogels

under cyclic loading. Accordingly, they provide a new evolution of damage variables for multiple

cycles, which is achieved by the experimental observations.

The concurrent exposure of rubber to fatigue and environmental aging has been mostly studied

through experimental efforts [198, 199, 200, 201, 202, 203, 204, 205]. So far, only a few phe-

nomenological models, and no micro-mechanical models, are available. In this respect, Moon et

al. [206] predicted the changes in the fatigue life of tires with respect to aging by implementing

the Arrhenius equation. Moreover, Neuhaus et al. [207] investigated the direct and indirect effects

of temperature on the response of the rubber-like material. They also proposed a new model which

114

can take into account the effects of thermal aging as well as the instantaneous effects of tempera-

ture on the material lifetime. Jouan and Constantinescu [208] conducted an experimental study to

investigate the effects of thermal aging on the fatigue life of silicone adhesive joints. Accordingly,

they proposed a law based on two mechanical predictors including the dissipated energy for adhe-

sive joints exposed to thermal aging. Besides, Wang et al. [209] proposed a relationship between

thermal aging conditions and the aging lifetime of carbon-black-filled styrene-isoprene-butadiene

rubber (SIBR-CB). To the best of the authors’ knowledge, currently, there are no physically based

models that can predict synergized effects of mechanical and environmental damages on the con-

stitutive response of polymeric matrix. Despite few successful models that can predict the consti-

tutive behavior of polymeric components exposed to either mechanical or environmental damages,

describing the material behavior in coupled aging conditions remains to be understood. Here, We

first formulate the constitutive behavior of virgin material using the network alteration concept and

then couple it with environmental aging and mechanical damages which are thermal aging and

low cycle fatigue. This paper is outlined as follows, first, the experimental characterization and

observations are discussed in section 6.2. Then the constitutive model is presented in 6.3. Finally,

validation and results are presented in section 6.4.

6.2 Experimental observations

6.2.1 Material and procedure

For this study, we used sulfur-vulcanized SBR sheets supplied by a local retailer. The SBR spec-

imens were cut into dog-bone samples using ASTM D412-C size die with its geometry given in

Fig. 6.1a. The experimental characterizations are briefly summarized here. In this respect, a uni-

axial universal testing machine (Test Resources 311 series frame) was used for quasi-static tensile

tests. All tests were displacement controlled and the extensometer grips holding each specimen

were 25.4 mm apart and all the experiments were performed at normal room temperature. To

induce thermal aging damage: Thermal aging progressively alters the rubber materials’ chemi-

115

a)

b)

c)

d)

Figure 6.1: (a) Detailed geometry of dog-bone sample. Constitutive behavior of a dog-bone SBR
sample damaged by exposure to (b) only thermal aging, (c) only low cycle fatigue, and (d) fatigue-
induced damage on thermally aged specimen stored at 80°C for 20 days.

116

25mmR=14mm115mm6mml=33mmR=25mml0=25.4mmA02.55123Stress (MPa)Stretch χUnagedT = 60oC, t = 10dT = 60oC, t = 20dT = 80oC, t = 10dT = 80oC, t = 20d11.11.200.511.522.5True Stress (MPa)Stretch Low-cycle fatigue11.11.200.511.522.5True Stress (MPa)Stretch (%)Thermal Aging (80�C 20days)+ cyclic fatiguecal composition, which subsequently modifies the structure of their macro-molecular network. In

the last decade, a number of investigations worked on identifying the different aging mechanisms

induced by thermo-oxidation [210, 103, 211, 212, 213, 214]. In contrast with the thermal treat-

ment (which usually should be done in an oxygen-free environment), thermal aging involves the

creation of free radicals P °. The free radicals P ° are extremely reactive with oxygen and con-

sequently favor a chain reaction of chemical processes in the polymer, denoted as propagation.

Mechanical characterization of material during aging could be categorized into intermittent and

contentious investigations. In intermittent tests for thermal aging, material samples are kept in an

air convection oven for several weeks or months at a constant temperature and with ambient oxy-

gen [192, 24, 215]. The aged samples are removed from the medium at predefined intervals and

submitted to a short-term test at room temperature[183]. Here, the intermittent test was utilized.

Accordingly, the aged samples are removed from the medium at predefined intervals and submitted

to a short-term test at room temperature [183]. The specimens were stored in stress-free state at

60°C or 80°C with 0%RH for 10 or 20 days (see Fig. 6.1b).

To induce low cycle fatigue damage: Aged and virgin samples were exposed to uni-axial

tensile cycles under displacement control settings, with a certain amplitude. In this respect, low

cycle fatigue is induced by non-relaxing cyclic loads with different amplitudes of 10% or 20% with

a cross-head speed of 50 mm/min. Here, we characterized accumulated damage after a certain

number of cycles namely 5, 50, 100, and 500 cycles by measuring the changes in constitutive

behavior under uni-axial loading (see Figs. 6.1c and 6.1d).

Table 6.1 shows the summary of all test conditions.

Table 6.1: Summary of coupled fatigue+thermal aging tests performed on SBR specimens, and the
changes in their constitutive behaviour.

Strain(%) Number of cycles (j) Aging time (days) Aging temperature (°C)

10%
20%

5, 50, 100, 500
5, 50, 100, 500

Unaged, 10, 20
Unaged, 10, 20

60, 80
60, 80

117

a)

b)

Figure 6.2: Schematic diagrams of evolution of constitutive response of elastomers during (a)
thermal aging, and (b) low cycle fatigue.

6.2.2 Accumulative damage

• Thermal damage: In view of thermal aging, two sub-structural phenomena namely chain

session and cross-link formation occur that make the material more brittle and reduce its

toughness (see Fig. 6.2a). In most elastomers, this phenomenon can be simply modeled via

first-order kinetics and it can be best approximated by the famous Arrhenius function [11].

• Fatigue damage: Due to fatigue-induced damage, the stress–stretch curves of elastomers

soften in each cycle toward a steady-state value often known as fatigue limit. As shown in

Fig. 6.2b, the peak stress level σ1

max is at its maximum in the first cycle, and continuously

decrease with number of cycles toward a steady-state σ∞

max. Since experimental characteri-

zation of σ∞

max is extremely costly, many studies are focused on derivation of σ∞

max through

a limited number of cycles. By normalizing the peak stress value to that of the first cycle

σ1
max, one can define a dimensionless, strain independent function

ˆσj =

σj
max(λi)
σ1
max(λi)

≈ j−ζ. ∀j,

(6.1)

where ζ is a material constant. In fact, this is an important function for the concept of intrin-

sic fracture toughness, which defines the minimum energy required for a crack to propagate.

118

σλThermal agingλf0J0J1J2λf1λf2J0>J1>J2Low-cycle fatigueσ 1 maxσ  maxσ j maxInterestingly, experimental data show that ˆσj is independent of the stretch amplitudes and

is an ever-decreasing function that depends only on j (see Fig. 6.3). It should be noted

that the extent of the normalized softening (beginning of the softening until reaching the

steady zone) depends on how far is this zone from the fracture. Accordingly, the normalized

softening plots till fracture have three-zone [216] :

i) In the first zone, commencing with the second cycle, a stress softening occurs until the

steady zone is reached.

ii) In the second zone, there is a minor fluctuation of the normalized maximum stress lying

up to achieve the critical damage amount.

iii) The third and final zone corresponds to a rapid decline in the normalized maximum stress

until the sample fails.

Here, since we mostly focused on the synergized effect of thermal aging and fatigue we

narrowed our model to the first zone and for that reason we choose a relatively small fatigue

amplitude. While the fatigue of metals has been extensively studied, there are very few

studies capable of describing fatigue phenomena in soft materials, especially with respect to

mechanical modeling [175, 197].

• Thermal + fatigue damage: Different sub-structural phenomena such as chain scission, or

cross-link formation are responsible for the damage accumulation in various conditions of

mechanical (chain breakage) and environmental degradation (chain scission + cross-link for-

mation). Thermal aging results in two sub-structural phenomena namely chain scission and

cross-link formation which made the material more brittle[217]. Accordingly, an increase

in apparent modulus during thermal aging is mainly attributed to the cross-link formation,

while the physical explanations for stress softening in cyclic fatigue could be bond rupture

[218], molecules slipping [219], filler rupture [220], and disentanglement [221, 63]. Sum-

marizing our experimental results in Fig.6.4, data suggests that the evolution of σj

max, in

quasi-static cyclic tests, depends mainly on the number of cycles j and status of thermal

119

a)

c)

b)

d)

Figure 6.3: Maximum Stress against number of cycles for different stretch amplitude: (a) σj
versus cycle number j for unaged material, (b) σj
°C for 20 days, (c) ˆσj
and (d) ˆσj
20 days.

max(λi)
max versus cycle number j for material aged at 80
max normalized maximum stress versus cycle number j for unaged material,
max normalized maximum stress versus cycle number j for for material aged at 80 °C for

120

010020030040050000.0050.010.0150.02Stress (MPa)Cycles cmax=1.1 cmax=1.2Unaged010020030040050000.0050.010.0150.02Stress (MPa)Cycles cmax=1.1 cmax=1.280(cid:176)C 20days010020030040050000.380.761.14Normalized StressCycles cmax=1.1 cmax=1.2Unaged010020030040050000.380.761.14Normalized StressCycles cmax=1.1 cmax=1.280(cid:176)C 20daysa)

c)

b)

d)

Figure 6.4: Maximum Stress against number of cycles for different aging conditions: (a) σj
sus cycle number j for λ = 1.1, (b) σj
maximum stress versus cycle number j for λ = 1.1, and (d) ˆσj
versus cycle number j for λ = 1.2.

max ver-
max normalized
max normalized maximum stress

max versus cycle number j for λ = 1.2, (c) ˆσj

damage (aging time and temperature). Normalizing the peak stress value to that of the first

cycle σ1

max for each aging condition, suggests that we can decouple thermal and fatigue dam-

age and represent fatigue damage as a function of the number of cycles (see Figs. 6.4c and

6.4d).

Following such experimental conclusion, our goal here is to predict the constitutive response of

material that are concurrently exposed to fatigue and thermal aging for different amount of time.

121

01002003004005000.00340.00510.00680.0085Stress (MPa)Cycles Unaged 60°C, t = 10d 60°C, t =  20d 80°C, t =  10d 80°C, t =  20dcmax=1.101002003004005000.00660.00990.01320.0165Stress (MPa)Cycles Unaged 60°C, t = 10d 60°C, t =  20d 80°C, t =  10d 80°C, t =  20dcmax=1.2010020030040050000.380.761.14Normalized StressCycles Unaged 60°C, t = 10d 60°C, t =  20d 80°C, t =  10d 80°C, t =  20dcmax=1.1010020030040050000.380.761.14Normalized StressCyclescmax=1.2 Unaged 60°C, t = 10d 60°C, t =  20d 80°C, t =  10d 80°C, t =  20d6.3 Constitutive model

The experimental results confirm that the behavior of an elastomeric sample exposed to concurrent

damage can be described as a hyper-elastic behavior with superimposed mechanical and environ-

mental damages. Accordingly, the damage accumulation in the elastomers is a contribution of two

phenomena; (i) environmental damage, and (ii) mechanical damage which is induced by fatigue.

The proposed approach allows us to integrate damage within any constitutive law that is based on

chain density N. Here for the sake of simplicity, we choose the full-network model to describe

the mechanical behavior of the SBR network [222], while other models like network alteration

[223, 224], or network evolution can be also used [7]. As shown in Fig. 6.5, the modular nature of

the proposed model can be best described by the decomposition of the behavior into three parts

1. Full-network model to describe the constitutive behaviour of unaged matrix [224].

2. Thermal aging model to describe the damage induced by thermal oxidation [11].

3. Fatigue model to describe the damage induced by fatigue during multiple cyclic loading

[197].

Figure 6.5: Schematic figure of modular nature of the proposed model.

Hereafter, the following notations are utilized: y as a scalar, y as a vector, Y as a second-order

tensor, and Y as a fourth-order tensor. Furthermore, we use the overhead index

d i
X to describe

122

OxygenTemperatureThermal agingLow-cycle fatigueUnaged polymeric networka)

b)

Figure 6.6: Schematic figures of a) the unit micro-sphere having a orientation vector d = d1e 1 +
d2e 2 + d3e 3, and b) the schematic representation of polymer network system.

certain entity X along direction d i.Note that we have reserved the superscript indices for power-

law relations and subscript indices for the state of the networks/matrix.

6.3.1 Full-network model

Assuming the elastomers to be isotropic and in-compressible, a representative volume element

(RVE) can be represented by a micro-sphere made by polymer chains homogeneously and equiv-

alently distributed in all spatial directions within the sphere. In each direction, one can assume a

similar set of chains with density per direction given as

d
N . Integrating the energy of chains in all

directions (see Fig. 6.6), one can write the volumetric strain energy of the RVE as

Ψ (F) =

d
N

1
4π

(cid:90)

S

d
ψd

d

A ∼=

1
4π

k
(cid:88)

i=1

d i
N wi

d i
ψ,

(6.2)

where

d
ψ is the strain energy of a single polymer chain in direction d , and F stands for the macro-

scale deformation gradient. This integration could be numerically achieved by the Gaussian inte-

gration scheme [145, 225, 226] (see Fig. 6.6b). Here, wi is the weight factor corresponding to each

integration directions d i (i = 1, 2, . . . , k). Since chains in different directions are subjected to dif-

123

ferent micro-stretches
d
λ as

chain with respect to

d
λ =

√

d TFTFd , it is important to rewrite the energy of a single polymer

d
ψ := ψ

(cid:18)d
λ,

(cid:19)

d
n

,

where

ψ (λ, n) = nKbT

d
n = n0 exp

(cid:19)

(cid:18)

d
λm

µ

and

d
N = N0 exp

(cid:18)

−µ

d
λm

(cid:18) λ
√
n
(cid:19)

β + ln

(cid:19)

β
sinh β

+ c0,

(6.3)

where Kb is Boltzmann’s constant, T temperature, and c0 the entropy of a chain in the reference

state. Here, n = L/l represents the normalized chain of the length with respect to the segment

Kuhn length l. It should be noted that, for the sake of simplicity, we assumed the homogeneous

distribution of the micro-stretch

d
λ within the network. Accordingly, the relationship between the

micro-stretch

d
λ and macro-stretch

d
χ in the direction d is assumed as

d
χ =

d
λ. The evolution of

damage in matrix can be considered to be a function of maximum stretches in each direction, e.g.
d
λm = max

d
λ(τ ). Deformation-induced irreversible damage appears only when the material is

τ ∈(−∞,t]

stretched to a level never experience before, or in other words when

and reloading, irreversible damage does not appear anymore and we have

d
λ. During unloading

d
λm =
d
λm = Const. Thus, the

evolution of matrix can be best modeled by representing parameters n and N as exponential func-

tions of

d
λm. Moreover n0, and µ are material parameters that will be determined through the fitting

procedure. Using the concept of Network alteration theory [224, 216], the full network model can

further be advanced to include the irreversible deformation-induced damage, by assuming the total

number of active polymer segments remain unchanged. To this end, one has

d
n

d
N = n0N0 = Const.

(6.4)

where n0 and N0 are both material parameters. Considering chains to have an original end-to-

end distance of R0 =

√

nl, and maximum end-to-end distance equal to the chain contour length

124

Rmax = L = nl, one can write

R =

√

d
λ

nl where R ≤ Rmax

and

d
λ ≤

√

n

(6.5)

The parameter β = L−1

(cid:19)

(cid:18) d
A

stands for inverse Langevin function, while

d
A =

d
λ R0
nl =

d
λ√
n

denotes the extensibility ratio in direction d . The inverse Langevin function L−1(x) could not be

explicitly derived and it should be approximated by a rational function. Here, a simple yet accurate

approximation is utilized [11] and it could be written as

L−1 (x) ∼=

x
1 − x

+ 2x −

8
9

x2.

(6.6)

6.3.2 Thermal model

In the course of thermal aging, we assumed that the changes in the chain density Nev follow the

same trend as the alteration of peroxide cross-link density ( see [11] for details). Hence, in the

absence of mechanical damages, the Ψaged can be written as

Ψaged =

1
4π

(cid:90)

S

d
N ev

d
ψd

d
A where

d
N ev = ρ∞ − ρ0 exp

(cid:18)

(cid:18)

−γ exp

−

(cid:19)

(cid:19)
t

,

Ea
RT

ρ0 < ρ∞

(6.7)

where ρ0, and ρ∞ are material parameters which describe the concentration of polymer chains

in two states of virgin and infinity aged matrix. Moreover, R = 8.314 [J.K −1.mol−1] is the

ideal gas constant, Ea [J.mol−1] the activation energy of the chemical reaction, and γ [(day)−1] a

degradation constant which is considered as a material constant. Figure 6.7a shows the evolution

of the normalized value of Nev through aging time and temperatures.

125

a)

b)

Figure 6.7: Evolution of the model parameters for multiple conditions; (a) Nev for a different aging
times, and (b) Nvib, nvib for a different number of cycles.

6.3.3 Fatigue model

In fatigue-induced damage, the irreversible damage sustained by matrix is mainly associated to

detachments/re-attachments of the physical bonds, the total number of which reduces over time

[227]. Gradual reduction of physical cross-links leads to increase of the average length of polymer

chain n and consequently decrease in the number of polymer chains N , since n.N = Const. After

the irreversible damage occurred in the first cycle, the evolution of the macro-molecular network is

mainly governed by the gradual loss of physical bonds which is also irreversible and its associated

changes in number N and length of polymer chains n. It is important to note that this assumption

is valid as long as experimental data confirms that the evolution of irreversible damage which is

governed through n and N , is fully independent of the fatigue stretch amplitude λi(see Figs. 6.3c

and 6.3d).

To describe the contribution of the irreversible damage in the jth cycle, we represent the com-

position of the damaged matrix through nvib and Nvib using power-law functions [197] as

Nvib
N

= j−ζ

considering

n.N = nvib.Nvib

⇒

nvib
n

= jζ.

(6.8)

It should be noted that the parameter ζ in Eqs. 6.8 and 6.1 are the same. Figure 6.7b depicts the

126

11.520255075100Normalized NevAging time (t) T = 80 °CT = 60 °CAging 0.80.911.11.20150300450Normalized NvibNormalized nvibNumber of Cycles (j)evolution of the normalized values of Nvib, and nvib through fatigue loading. Note that here, the

total number of segments is altered during thermal aging while the length of the chains remain

unchanged. In fact, keeping the chains length unchanged during thermal aging is a simplification

assumption to keep the number of material parameters low. In cyclic fatigue, the length of the

polymer chain N increases while the number of the active chain n decreases which makes the

total number of segments remain unchanged. This means that during mechanical loading conser-

vation of the total number of segments is held while along the aging trajectory the total number of

segments is altered due to the change of the polymer matrix [59]..

6.3.4 Damage superposition

The main assumption here, is that mechanical and environmental aging collaboratively induce

damage and thus, they can be modeled as two parallel damage mechanisms. Here, damage accu-

mulation is given in terms of multiplicative decomposition of damages induced by thermal aging

and fatigue-induced damages ( see Fig. 6.5 ). Accordingly, accumulated damage within the matrix

can be described by inserting Eqs. 6.3, 6.7, and Eq. 6.8 into 6.2 which yields

ΨD =

d i
˜N =

1
4π

k
(cid:88)

d i
˜N ψ

wi

i=1

(cid:18)d i

λ, n0jζ exp

(cid:19)(cid:19)

(cid:18)

d i
λ m

µ

, where

1
(cid:18)

µ

(cid:19)

d i
λ m

jζ exp

(cid:18)

˜N∞ − ˜N0 exp

(cid:18)

(cid:18)

−γ exp

−

(cid:19)

(cid:19)(cid:19)
t

.

Ea
RT

(6.9)

where the two material parameters ˜N∞ and ˜N0 are defined by multiplication of other constants,

namely

˜N∞ = ρ∞N0

˜N0 = ρ0N0,

and thus both can be directly derived from fitting.

127

6.3.5 Constitutive formulation

To derive a constitutive equation for an in-compressible polymer matrix (i.e. J = det F = 1), the

strain energy function can be postulated by

ΨM = ΨD − p(J − 1),

(6.10)

where the scalar p can be regarded as an indeterminate Lagrange multiplier, which is gener-

ally known as hydro-static pressure and it can be determined from boundary conditions ( detailed

derivation of p for the plane stress condition is discussed in [11]). The first Piola-Kirchhoff ten-

sor P can be achieved by differentiating the Eq. 6.10 with respect to the deformation gradient F

[228, 229, 230] and using the identity of ∂J

∂F = JF−T . Hence, P can be written as

P =

∂ΨM (F)
∂F

− pF−T ,

∂ΨM (F)
∂F

k
(cid:88)

d i
˜N

wi

=

i=1

d i
ψ
d i
λ

∂

∂

1
d i
λ

2

∂d i

¯Cd i
∂ ¯F

∂ ¯F
∂F

,

:

(6.11)

where C denotes for the right Cauchy–Green tensor, ¯F = J −1/3F, and ¯C = J −2/3C. Besides, the

Eq. 6.11 could be even more simplified by using the following identities

d i
ψ

∂

(cid:18)

n,

(cid:19)

d i
λ

d i
λ

∂

√

=

nKbT β





d i
λ
√
n



, n

 ,

∂d ¯Cd
∂ ¯F

∂ ¯F
∂F

:

= 2 ¯F(d ⊗ d ) : J − 1

3 I = 2J − 1

3 ¯F(d ⊗ d ),

P = KbT

√

nd iL−1

k
(cid:88)

d i
˜N

i=1





d i
λ
nd i

√





wi
d i
λ

J − 1

3 ¯F(d i ⊗ d i) − pF−T .

(6.12)

6.3.6 FE linearization

To analyze a boundary value problem with the proposed model, the macro-scale domain is as-

sumed as a homogeneous continua body analyzed in the framework of the nonlinear finite element

128

method [231], as shown in Fig. 6.8. In what follows, we present the consistent tangent modu-

Figure 6.8: A schematic figure of the hierarchical micro-macro multi-scale nature of the model.

lus C of the proposed constitutive model for an FE implementation. Accordingly, one can derive

the Lagrangian elasticity tensor C by differentiation of the second Piola–Kirchhoff stress tensor

S = F−1P with respect to the Green-Lagrange strain tensor E = 1

). Having the plane stress condition and defining 2 ∂Ψ(F)

∂C , Cn = 4 ∂2W (F)
∂ ¯C2



2(FT F − I) ( see [11] for details






0 0 F2
33

, A =

0 0

0 0

0

0

,

129

ΓuΓttuFC,SxxxxxxxxxxxxxxxxGaussIntegrationpointDiscretized domainInitial configurationDeformed configurationB =






0 0

0 0

0

0

0 0 2F33






, one can derive the Lagrangian tangent modulus as:

C =

−2

3 (cid:0)Sn − Sn

33F2

33C−1(cid:1) ⊗ C−1

−2
3
−4
3

J
(cid:18)

+ J

Cn −

Cn : C ⊗ C−1 − C−1 ⊗ Cn : A +

1
3

(C : Cn : A) C−1 ⊗ C−1

(cid:19)

(6.13)

1
3
(cid:19)

− 2Sn : AC−1 ⊗ C−1

,

− J

−2
3

(cid:18)

2Sn : BC−1 ⊗

∂F33
∂C

where

∂F33
∂C

=

∂F33
∂F

∂F
∂C

=

1
2

F−1F′, where F′ =

∂F33
∂F

= F2
33






In order to determine Cn and Sn, we have

−F22 F21

F12 −F11 0

0

0


0




0

Snd ⊗ d = 2

∂Ψ (F)
∂C

k
(cid:88)

d i
˜N wi

∂

∼= 2

i=1

d i
λ)

d i
ψ(
∂C

Cnd ⊗ d ⊗ d ⊗ d = 4

∂2Ψ (F)
∂C2

k
(cid:88)

d i
˜N wi

∼= 4

i=1

∂2

d i
λ)

d i
ψ(
∂C2

.

.

The derivatives used in the above equations, could be derived as follows

.

(6.14)

(6.15)

(6.16)

(cid:19)

(cid:18)d i
λ

∂2ψ

√

= KbT

∂C2

n
2

d i
4
λ

(cid:19)

(cid:18)d i
λ

∂ψ

= KbT

∂C

(cid:33)

β

d ⊗ d .

(cid:32)√

n
d i
2
λ

d ⊗ d ⊗ d ⊗ d .

(6.17)

(cid:33)

(cid:32)

∂β
d i
λ

∂

−

β
d i
λ

130

Table 6.2: Material parameters of the proposed model for SBR.

˜N0[M P a] ˜N∞[M P a] n
1.3700

2.5007

µ

Ea[j][mol]−1
2.687 0.52 0.06218 8.5860e+10 9.7207e+03

γ[day]−1

ζ

=

∂β
d i
λ

∂

−16A3 + 50A2 − 52A + 27

√
9

n (A − 1)2

where A =

d i
λ
√
n

.

(6.18)

6.4 Validation and results

To assess the proposed model’s capability, we compared its predictions to our experimental results,

which were specifically designed to demonstrate the negative impact of chemical and physical

aging on the constitutive response and energy absorption of polymer matrix. The model was fitted

for the following three cases using the 7 material parameters:

• The loading-unloading curve of the virgin specimen to investigate the effects of λ and λm.

• Two points of loading curves of 10 days chemically aged specimens at temperatures 60°C,

and 80°C to study the effects of T , and t

• One points of loading curve of physically aged specimen at cycles j = 200 to study the

effects of fatigue.

The least square error function was minimized using the Levenberg–Marquardt algorithm. Ac-

cordingly, Table 6.2 lists the so-obtained material parameters of the proposed model for the SBR.

Besides, Figs. 6.9, and 6.10 show the prediction of the devised model against experimental results

for various aging times, and the number of cycles for temperatures 60 and 80 °C, respectively.

The devised model can predict the hardening response of material due to thermal aging as well

as softening due to fatigue damage. Besides, it can also predict the behavior when both damage

mechanisms are accumulated on top of each other. One should notice that the current terminol-

ogy used to associate the damage to the alteration of chain density and chain length is based on

131

a)

b)

c)

Figure 6.9: Validation of the model predictions against our experimental data for SBR in the
unaged state, and during thermal aging for 10, and 20 days in the temperature 60°C in multiple
conditions; (a) the constitutive response of samples including their corresponding failure points,
(b) maximum stress versus the number of cycles with maximum stretch of χmax = 1.1, and (c)
maximum stress as versus the number of cycles with a maximum stretch of χmax = 1.2.

a statistical representation of the matrix and has not been experimentally validated. In this study,

the proposed evolution describing the alteration of the macro-molecular network through n and

N are mainly derived based on our macro-scale observations of the mechanical behavior and our

familiarity with the aging mechanisms.

To illustrate the performance and applicability of the derived model in modeling the mechan-

ical behavior of polymers during concurrent exposure to thermal aging and fatigue, a numerical

example is solved and presented here. The elasticity tensor, which is generated from the strain en-

ergy function mentioned in previous sections, is used to model the materials’ nonlinear behavior.

132

02.55123Stress (MPa)Stretch χUnagedT = 60oC, t = 10dT = 60oC, t = 20dExperimental DataFitted curveFitted point010100200300400500Stress (MPa)Cycle (j)T = 60oC, t = 10dT = 60oC, t = 20dExperimental DataUnagedχmax=1.1Fitted point0120100200300400500Stress (MPa)Cycle (j)T = 60oC, t = 10dT = 60oC, t = 20dExperimental DataUnagedχmax=1.2a)

b)

c)

Figure 6.10: Validation of the model predictions against our experimental data for SBR in the
unaged state, and during thermal aging for 10, and 20 days in the temperature 80°C in multiple
conditions; (a) the constitutive response of samples including their corresponding failure points,
(b) maximum stress versus the number of cycles with maximum stretch of χmax = 1.1, and (c)
maximum stress as versus the number of cycles with a maximum stretch of χmax = 1.2.

133

02.55123Stress (MPa)Stretch χUnagedT = 80oC, t = 10dT = 80oC, t = 20dExperimental DataFitted point020100200300400500Stress (MPa)Cycle (j)T = 80oC, t = 10dT = 80oC, t = 20dExperimentalDataUnagedχmax = 1.1020100200300400500Stress (MPa)Cycle (j)T = 80oC, t = 10dT = 80oC, t = 20dExperimentalDataUnagedχmax=1.2Accordingly, the free edge of a perforated square plate with a length of 50 mm is subjected to a

uni-axial deformation. Three evenly spaced holes with radii of 5, 6, and 8 mm are included on the

plate ( the detailed geometry is provided in Fig. 6.11a ). Plane stress conditions are assumed, and

quadratic elements with four Gauss points are used ( see Fig. 6.11b ). The right edge of the plate

Figure 6.11: Finite element simulation of a perforated plate under uni-axial deformation: a) de-
tailed geometry and boundary conditions; b) finite-element discretization.

is subjected to tensile deformation of 25 mm over 50 load stages, while the left edge remains sta-

tionary. The boundary value problem was solved with the different numbers of elements (i.e. Ne =

292, 891, 1195, and 1796). As can be seen in Fig. 6.12, the the maximum von Mises stress of the

boundary value problem converged at a reasonable limit. Accordingly, for optimal computational

efficiency, we solved our boundary value problem with 1195 elements.

Figures 6.13, 6.14, and 6.15 illustrate the evolution of second-Piola Kirchhoff stress S contours

of the perforated plate under uni-axial deformation during thermal aging and fatigue using the

suggested model. The results show that increasing the number of cycles reduces stress in the

specimen while increasing the aging period makes the material stiffer. Accordingly, the largest

value of Sxx occurs at the top and bottom margins of circles, particularly the bigger circle for

the environmentally aged specimen. Similarly, the Sxx takes its minimum at the right and bottom

edges of circles for the mechanically damaged sample.

134

Figure 6.12: Number of elements versus the maximum von Mises stress of the boundary value
problem.

Figure 6.13: Evolution of the second-Piola Kirchhoff stress Sxx of the SBR through thermal aging
and fatigue. As it can be seen, higher aging duration caused higher stress values due to the hard-
ening effect of thermal aging while more number of cycles cause the stress softening.

135

Figure 6.14: Evolution of the second-Piola Kirchhoff stress Sxy of the SBR through thermal aging
and fatigue. As it can be seen, higher aging duration caused higher stress values due to the hard-
ening effect of thermal aging while more number of cycles cause the stress softening.

136

Figure 6.15: Evolution of the second-Piola Kirchhoff stress Syy of the SBR through thermal aging
and fatigue. As it can be seen, higher aging duration caused higher stress values due to the hard-
ening effect of thermal aging while more number of cycles cause the stress softening.

137

6.5 Conclusion

In this paper, we proposed a micro-mechanical constitutive model based on the concepts of dam-

age accumulation to describe the damage induced by temperature and low cycle fatigue on the

constitutive response of cross-linked polymers. Each of those damage mechanisms is induced by

mechanisms that have been extensively studied before but have not been modeled in parallel so

far. We modeled each sub-mechanism through a linear kinetic equation which is then coupled into

the modular network concept to allow calculation of synergies and parallel damage accumulation

through multiple sources. This model could be regarded as an extension of our recent models

[7, 11, 175]. Furthermore, the devised model is applied for FE simulation and the tangent modulus

for large deformation regimes are expressed in terms of internal variables. A full set of experi-

mental data was carried out to validate the model. With reasonable accuracy, the proposed model

shows promising results. In addition, with meaningful parameters, the model is relatively simple.

138

CHAPTER 7

Summary and Future Works

The main goal of this research was to create a constitutive model for cross-linked polymers under

various aging scenarios. Each section of the dissertation is briefly summarized in this chapter. An

introduction to cross-linked polymers and their applications was presented in the first part of this

dissertation.

7.1 General Remarks

• In chapter 3, effects of Hydrolytic aging on cross-linked polymer behavior were investi-

gated. Following that, a comprehensive study of hydrolysis kinetics was presented. A micro-

mechanical model to accurately predict the inelastic behavior of elastomers during aging has

been presented later, combining the concepts of network evolution and the dual network hy-

pothesis. Finally, the model’s capabilities were tested by comparing model predictions to

our experimental data.

• In chapter 4, a computationally efficient model is proposed to capture the loss of mechanical

performance due to chemical aging that is formed as the competition of chain scission and

cross-link formation/dissolution, such as thermo-oxidative aging or hydrolytic aging. Later

on, the model is linearized for implementation into Finite Element (FE) simulations. The

proposed model includes only six physically inspired material parameters. Thus, while it is

computationally efficient, it shows good agreement with its own experimental data, which

was performed on various ranges of accelerated aging temperatures and times.

• In chapter 5, a multi-physics model is proposed to predict the changes in the constitutive

139

behavior of cross-linked polymeric systems due to damages induced by deformation, oxida-

tion, moisture, and temperature. Coupling the concept of network evolution, hydrolysis, and

thermo-oxidation aging, we hypothesize that the synergized effects of deformation-induced

damage, as well as different environmental factors such as humidity, temperature, and oxy-

gen, can be taken into account through a generic model. Finally, the model is validated with

respect to extensive sets of our experimental data. In view of its interoperability, precision,

and deep insight it provides into the nature of the aging phenomenon, the model is a good

choice for advanced implementation in FE applications.

• In chapter 6, a micro-mechanical constitutive model is presented to predict the concurrent

effects of thermal aging and cyclic fatigue on the constitutive behavior of cross-linked poly-

mers. The damage associated with each of those aging conditions is induced by mechanisms

that have been extensively studied individually. Here, the main goal was to model the damage

accumulated when those mechanisms work in parallel. Later on, the model is implemented

into finite element simulations. For validation, the devised model is bench-marked against a

comprehensive set of experimental data. The proposed model shows promising results with

reasonable precision.

7.2 Future Research

Aging can be either the result of individual aging elements or due to combined synergized effects

of the different environmental factors. Since these factors coexist in nature, natural aging is seen as

a combined effect of all conditions together. Lately, solo environmental damage models have been

of significant interest in the literary and research communities. However, there are still wide gaps in

the modeling of coupled-aging phenomena since finding the mutual effect and correlation of each

environmental element is challenging and needs a huge set of experimental observations. Accord-

ingly, realizing the correlation of the mechanical and coupled-environmental elements (damage

accumulation) is crucial for providing a realistic prediction of such materials’ behavior.

140

BIBLIOGRAPHY

[1] Lallit Anand and Sanjay Govindjee. Continuum mechanics of solids. Oxford University

Press, 2020.

[2] Gerhard A. Holzapfel. Nonlinear Solid Mechanics: A Continuum Approach for Engineer-

ing. John weily & Sons, 2005.

[3] Scott W Case, Arit Das, Michael J Bortner, David A Dillard, and Christopher C White.
Durability and accelerated characterization of adhesive bonds. In Advances in Structural
Adhesive Bonding, pages 675–710. Elsevier, 2023.

[4] Amir Bahrololoumi, Mamoon Shaafaey, Georges Ayoub, and Roozbeh Dargazany. A failure
model for damage accumulation of cross-linked polymers during parallel exposure to ther-
mal aging & fatigue. International Journal of Non-Linear Mechanics, 146:104142, 2022.

[5] Rod Martin. Ageing of composites. Elsevier, 2008.

[6] Amir Bahrololoumi, Mamoon Shaafaey, Georges Ayoub, and Roozbeh Dargazany. Thermal
aging coupled with cyclic fatigue in cross-linked polymers: Constitutive modeling & fe
implementation. International Journal of Solids and Structures, 252:111800, 2022.

[7] Amir Bahrololoumi, Vahid Morovati, Mamoon Shaafaey, and Roozbeh Dargazany. A multi-
physics approach on modeling of hygrothermal aging and its effects on constitutive behavior
of cross-linked polymers. Journal of the Mechanics and Physics of Solids, page 104614,
2021.

[8] Mamoon Shaafaey, Amir Bahrololoumi, Hamid Mohammadi, Sharif Alazhary, and Roozbeh
Dargazany. Investigation of hygrothermal aging on the polyurethane-based (pub) adhesive:
substantiating competition scenario between sub-aging thermo-oxidation and hydrolytic
phenomena. Journal of Polymer Research, 28(12):1–25, 2021.

[9] Amir Bahrololoumi, Vahid Morovati, Emad A Poshtan, and Roozbeh Dargazany. A multi-
physics constitutive model to predict quasi-static behaviour: Hydrolytic aging in thin cross-
linked polymers. International Journal of Plasticity, page 102676, 2020.

[10] Wanru Miao, Amir Bahrololoumi, Shaafaey Mamoon, Mohammadi Hamid, Omid Nabine-
jad, and Roozbeh Dargazany. Effects of hydrolytic aging on constitutive behavior of silicone
adhesives in saline and desalinated water. Submitted to Polymer Degradation and Stability,
2020.

[11] Amir Bahrololoumi, Hamid Mohammadi, Vahid Moravati, and Roozbeh Dargazany. A
physically-based model for thermo-oxidative and hydrolytic aging of elastomers. Interna-
tional Journal of Mechanical Sciences, 194:106193, 2021.

[12] Hamid Mohammadi, Vahid Morovati, Emad Poshtan, and Roozbeh Dargazany. Under-
standing decay functions and their contribution in modeling of thermal-induced aging of
cross-linked polymers. Polymer Degradation and Stability, page 109108, 2020.

141

[13] Philip A Schweitzer et al. Fundamentals of corrosion: mechanisms, causes, and preventa-

tive methods. CRC press, 2009.

[14] Ján Kruželák, Ivan Hudec, and Rastislav Dosoudil. Influence of thermo-oxidative and ozone
ageing on the properties of elastomeric magnetic composites. Polymer degradation and
stability, 97(6):921–928, 2012.

[15] M. Gigliotti, M. Minervino, M.C. Lafarie-Frenot, and J.C. Grandidier. Effect of thermo-
oxidation on the local mechanical behaviour of epoxy polymer materials for high tempera-
ture applications. Mechanics of Materials, 101:118 – 135, 2016.

[16] R. Behnke and M. Kaliske. Numerical modeling of thermal ageing in steady state rolling

tires. International Journal of Non-Linear Mechanics, 2018.

[17] Roozbeh Dargazany, Vu Ngoc Khiêm, and Mikhail Itskov. A generalized network decom-
International

position model for the quasi-static inelastic behavior of filled elastomers.
Journal of Plasticity, 63:94, 2015.

[18] LR G Treloar. The physics of rubber elasticity. 1975.

[19] B. D. Coleman and W. Noll. The thermodynamics of elastic materials with heat conduction

and viscosity. Archive for Rational Mechanics and Analysis, 13:167, 1963.

[20] B. Coleman and M. E. Gurtin. Thermodynamics with internal state variables. The Journal

of Chemical Physics, 47:593, 1967.

[21] Sima Kashi, Russell Varley, Mandy De Souza, Salwan Al-Assafi, Adriano Di Pietro, Chris-
telle de Lavigne, and Bronwyn Fox. Mechanical, thermal, and morphological behavior of
silicone rubber during accelerated aging. Polymer-Plastics Technology and Engineering,
2018.

[22] Michael Johlitz. On the representation of ageing phenomena. The Journal of Adhesion,

88(7):620–648, 2012.

[23] Tianying Chang, Xiansheng Zhang, and Hong-Liang Cui. Evolution of terahertz dielectric
permittivity of rubber during thermo-oxidative aging. Microwave and Optical Technology
Letters, 60, 12 2017.

[24] Michael Johlitz, Nico Diercks, and Alexander Lion. Thermo-oxidative ageing of elastomers:
A modelling approach based on a finite strain theory. International Journal of Plasticity,
63:138–151, 2014.

[25] Chao Yu, Guozheng Kang, and Kaijuan Chen. A hygro-thermo-mechanical coupled cyclic
constitutive model for polymers with considering glass transition. International Journal of
Plasticity, 89:29–65, 2017.

[26] André C Vieira, Rui M Guedes, and Volnei Tita. Damage-induced hydrolyses modelling
of biodegradable polymers for tendons and ligaments repair. Journal of biomechanics,
48(12):3478–3485, 2015.

142

[27] Amir Bahrololoumi and Roozbeh Dargazany. Hydrolytic aging in rubber-like materials: a
In American Society of Mechanical Engineers

micro-mechanical approach to modeling.
(ASME) 2019 International Mechanical Engineering Congress and Exposition, 2019.

[28] Marcos Lugo, Jason E Fountain, Justin M. Hughes, Jean-Luc Bouvard, and Mark F Horste-
meyer. Microstructure-based fatigue modeling of an acrylonitrile butadiene styrene (abs)
copolymer. Journal of Applied Polymer Science, 131(20), 2014.

[29] Antoine Launay, Yann Marco, MH Maitournam, Ida Raoult, and Fabien Szmytka. Cyclic
behavior of short glass fiber reinforced polyamide for fatigue life prediction of automotive
components. Procedia Engineering, 2(1):901–910, 2010.

[30] Amir K Shojaei and Pieter Volgers. Fatigue damage assessment of unfilled polymers in-

cluding self-heating effects. International Journal of Fatigue, 100:367–376, 2017.

[31] Wulff Possart and Markus Brede. Adhesive Joints: Ageing and Durability of Epoxies and

Polyurethanes. John Wiley & Sons, 2019.

[32] Mayank Sinha, Alborz Izadi, Rebecca Anthony, and Sara Roccabianca. A novel approach to
finding mechanical properties of nanocrystal layers. Nanoscale, 11(15):7520–7526, 2019.

[33] Kaspar Loeffel and Lallit Anand. A chemo-thermo-mechanically coupled theory for elastic–
viscoplastic deformation, diffusion, and volumetric swelling due to a chemical reaction.
International Journal of Plasticity, 27(9):1409–1431, 2011.

[34] Marco Gigliotti and Jean-Claude Grandidier. Chemo-mechanics couplings in polymer ma-
trix materials exposed to thermo-oxidative environments. Comptes Rendus Mécanique,
338(3):164–175, 2010.

[35] Marco Gigliotti, Jean-Claude Grandidier, and Marie Christine Lafarie-Frenot. Assessment
of chemo-mechanical couplings in polymer matrix materials exposed to thermo-oxidative
environments at high temperatures and under tensile loadings. Mechanics of Materials,
43(8):431–443, 2011.

[36] Michael Johlitz and Alexander Lion. Chemo-thermomechanical ageing of elastomers
based on multiphase continuum mechanics. Continuum Mechanics and Thermodynamics,
25(5):605–624, 2013.

[37] Alexander Lion and PHILIPP Höfer. On the phenomenological representation of curing

phenomena in continuum mechanics. Archives of Mechanics, 59(1):59–89, 2007.

[38] K Kannan and KR Rajagopal. A thermodynamical framework for chemically reacting sys-

tems. Zeitschrift für angewandte Mathematik und Physik, 62(2):331–363, 2011.

[39] Hamid Mohammadi and Roozbeh Dargazany. A micro-mechanical approach to model ther-
mal induced aging in elastomers. International Journal of Plasticity, 118:1–16, 2019.

143

[40] H Mohammadi, A Bahrololoumi, Y Chen, and R Dargazany. A micro-mechanical model for
constitutive behavior of elastomers during thermo-oxidative aging. In Constitutive Models
for Rubber XI: Proceedings of the 11th European Conference on Constitutive Models for
Rubber (ECCMR 2019), June 25-27, 2019, Nantes, France, page 542. CRC Press, 2019.

[41] Prashant Saxena, Mokarram Hossain, and Paul Steinmann. A theory of finite deformation
magneto-viscoelasticity. International Journal of Solids and Structures, 50(24):3886–3897,
2013.

[42] Jianyou Zhou, Liying Jiang, and Roger E Khayat. A micro–macro constitutive model for
finite-deformation viscoelasticity of elastomers with nonlinear viscosity. Journal of the
Mechanics and Physics of Solids, 110:137–154, 2018.

[43] GP Potirniche, A Pascu, N Shoemaker, PT Wang, MF Horstemeyer, D Stillman, and T-L
Lin. A visco-hyperelastic model for the thermo-mechanical behavior of polymer fibers.
International Journal of Damage Mechanics, 20(7):1002–1020, 2011.

[44] Aref Ghaderi and Roozbeh Dargazany. A data-driven model to predict constitutive and fail-
ure behavior of elastomers considering the strain rate, temperature, and filler ratio. Journal
of Applied Mechanics, 90(5):051010, 2023.

[45] H Khajehsaeid, M Baghani, and R Naghdabadi. Finite strain numerical analysis of elas-
tomeric bushings under multi-axial loadings: a compressible visco-hyperelastic approach.
International Journal of Mechanics and Materials in Design, 9(4):385–399, 2013.

[46] Aurelien Maurel-Pantel, Erwan Baquet, Jérôme Bikard, Jean-Luc Bouvard, and Noelle Bil-
lon. A thermo-mechanical large deformation constitutive model for polymers based on
material network description: Application to a semi-crystalline polyamide 66. International
Journal of Plasticity, 67:102–126, 2015.

[47] Babak Farrokh and Akhtar S Khan. A strain rate dependent yield criterion for isotropic poly-
mers: low to high rates of loading. European Journal of Mechanics-A/Solids, 29(2):274–
282, 2010.

[48] Georges Ayoub, Fahmi Zaïri, Moussa Naït-Abdelaziz, Jean Michel Gloaguen, and G Kridli.
A visco-hyperelastic damage model for cyclic stress-softening, hysteresis and permanent set
in rubber using the network alteration theory. International Journal of Plasticity, 54:19–33,
2014.

[49] Akhtar S Khan, Oscar Lopez-Pamies, and Rehan Kazmi. Thermo-mechanical large defor-
mation response and constitutive modeling of viscoelastic polymers over a wide range of
strain rates and temperatures. International Journal of Plasticity, 22(4):581–601, 2006.

[50] Jingzhe Pan. Modelling degradation of bioresorbable polymeric medical devices. Elsevier,

2014.

[51] Yanqun Yu, Zongyu Chang, Yaoguang Qi, and Dehua Feng. Research on new technology
for offshore heavy oil thermal recovery with rod pumping. Journal of Petroleum Exploration
and Production Technology, 8(3):947–955, 2018.

144

[52] Leif Kari. Dynamic stiffness of chemically and physically ageing rubber vibration isolators
in the audible frequency range. Continuum Mechanics and Thermodynamics, 29(5):1027–
1046, 2017.

[53] André C Vieira, Rui M Guedes, and Volnei Tita. Constitutive modeling of biodegradable
polymers: Hydrolytic degradation and time-dependent behavior. International Journal of
Solids and Structures, 51(5):1164–1174, 2014.

[54] Quentin Breche, Grégory Chagnon, Guilherme Machado, Benjamin Nottelet, Xavier Garric,
Edouard Girard, and Denis Favier. A non-linear viscoelastic model to describe the mechani-
cal behavior’s evolution of biodegradable polymers during hydrolytic degradation. Polymer
degradation and stability, 131:145–156, 2016.

[55] Thorsten Pretsch, Ines Jakob, and Werner Müller. Hydrolytic degradation and functional
stability of a segmented shape memory poly (ester urethane). Polymer Degradation and
Stability, 94(1):61–73, 2009.

[56] C Slater, C Davis, and M Strangwood. Compression set of thermoplastic polyurethane un-
der different thermal–mechanical-moisture conditions. Polymer degradation and stability,
96(12):2139–2144, 2011.

[57] DF Farrar and RK Gillson. Hydrolytic degradation of polyglyconate b: the relationship
between degradation time, strength and molecular weight. Biomaterials, 23(18):3905–3912,
2002.

[58] Emmanuel Richaud, Pierre Gilormini, Marie Coquillat, and Jacques Verdu. Crosslink den-
sity changes during the hydrolysis of tridimensional polyesters. Macromolecular Theory
and Simulations, 23(5):320–330, 2014.

[59] A.V. Tobolsky, I.B. Prettyman, and J.H. Dillon. Stress relaxation of natural and synthetic

rubber stocks. Journal of Applied Physics, 15(4):380–395, 1944. cited By 87.

[60] Sabine Cantournet, Rodrigue Desmorat, and Jacques Besson. Mullins effect and cyclic
stress softening of filled elastomers by internal sliding and friction thermodynamics model.
International Journal of Solids and Structures, 46(11-12):2255–2264, 2009.

[61] Andri Andriyana, Mei Sze Loo, Grégory Chagnon, Erwan Verron, and Shiau Ying Ch’ng.
Modeling the mullins effect in elastomers swollen by palm biodiesel. International Journal
of Engineering Science, 95:1–22, 2015.

[62] RW Ogden and DG Roxburgh. A pseudo–elastic model for the mullins effect in filled
rubber. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and
Engineering Sciences, 455(1988):2861–2877, 1999.

[63] Julie Diani, Bruno Fayolle, and Pierre Gilormini. A review on the mullins effect. European

Polymer Journal, 45(3):601–612, 2009.

145

[64] Morgan Deroiné, Antoine Le Duigou, Yves-Marie Corre, Pierre-Yves Le Gac, Peter Davies,
Guy César, and Stéphane Bruzaud. Accelerated ageing of polylactide in aqueous environ-
ments: comparative study between distilled water and seawater. Polymer degradation and
stability, 108:319–329, 2014.

[65] Hideto Tsuji and Yoshito Ikada. Properties and morphology of poly (l-lactide) 4. effects
of structural parameters on long-term hydrolysis of poly (l-lactide) in phosphate-buffered
solution. Polymer Degradation and Stability, 67(1):179–189, 2000.

[66] Hideto Tsuji and Yoshito Ikada. Properties and morphology of poly (l-lactide). ii. hydrolysis
in alkaline solution. Journal of Polymer Science Part A: Polymer Chemistry, 36(1):59–66,
1998.

[67] Jacob Song Kiat Lim, Chee Lip Gan, and Xiao Matthew Hu. Unraveling the mechanistic

origins of epoxy degradation in acids. ACS Omega, 4(6):10799–10808, 2019.

[68] Suming Li. Hydrolytic degradation characteristics of aliphatic polyesters derived from lactic
and glycolic acids. Journal of Biomedical Materials Research: An Official Journal of The
Society for Biomaterials, The Japanese Society for Biomaterials, and The Australian Society
for Biomaterials, 48(3):342–353, 1999.

[69] SB Harogoppad and TM Aminabhavi. Diffusion and sorption of organic liquids through
polymer membranes. ii. neoprene, sbr, epdm, nbr, and natural rubber versus n-alkanes. Jour-
nal of applied polymer science, 42(8):2329–2336, 1991.

[70] Suming Li. Hydrolytic degradation characteristics of aliphatic polyesters derived from lactic
and glycolic acids. Journal of Biomedical Materials Research: An Official Journal of The
Society for Biomaterials, The Japanese Society for Biomaterials, and The Australian Society
for Biomaterials, 48(3):342–353, 1999.

[71] Achim Gopferich and Robert Langer. Modeling of polymer erosion. Macromolecules,

26(16):4105–4112, 1993.

[72] Man Wang, Xiaowei Xu, Jin Ji, Yang Yang, Jianfeng Shen, and Mingxin Ye. The hy-
grothermal aging process and mechanism of the novolac epoxy resin. Composites Part B:
Engineering, 107:1–8, 2016.

[73] Rong Long, H Jerry Qi, and Martin L Dunn. Modeling the mechanics of covalently adapt-
able polymer networks with temperature-dependent bond exchange reactions. Soft Matter,
9(15):4083–4096, 2013.

[74] L.R.G. Treloar. The Physics of Rubber Elasticity. Oxford University Press, 2005.

[75] C. Miehe, S. Göktepe, and F. Lulei. 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, 2004.

146

[76] Reza Rastak and Christian Linder.

A non-affine micro-macro approach to strain-
crystallizing rubber-like materials. Journal of the Mechanics and Physics of Solids, 111:67–
99, 2018.

[77] Hüsnü Dal, Christoph Zopf, and Michael Kaliske. Micro-sphere based viscoplastic con-
stitutive model for uncured green rubber. International Journal of Solids and Structures,
132:201–217, 2018.

[78] Roozbeh Dargazany and Mikhail Itskov. Constitutive modeling of the mullins effect and
cyclic stress softening in filled elastomers. Physical Review E, 88(1):012602, 2013.

[79] A.E. Ehret, M. Itskov, and H. Schmid. Numerical integration on the sphere and its effect
on the material symmetry of constitutive equations-Â a comparative study. International
Journal for Numerical Methods in Engineering, 81:189, 2010.

[80] Erwan Verron. Questioning numerical integration methods for microsphere (and mi-

croplane) constitutive equations. Mechanics of Materials, 89:216–228, 2015.

[81] Vahid Morovati and Roozbeh Dargazany. Micro-mechanical modeling of the stress soften-
ing in double-network hydrogels. International Journal of Solids and Structures, 164:1–11,
2019.

[82] Vahid Morovati and Roozbeh Dargazany. Improved approximations of non-gaussian proba-
bility, force, and energy of a single polymer chain. Physical Review E, 99(5):052502, 2019.

[83] Vahid Morovati, Hamid Mohammadi, and Roozbeh Dargazany. A generalized approach to
improve approximation of inverse langevin function. In ASME 2018 International Mechani-
cal Engineering Congress and Exposition, pages V009T12A029–V009T12A029. American
Society of Mechanical Engineers, 2018.

[84] Vahid Morovati, Hamid Mohammadi, and Roozbeh Dargazany. A generalized approach
to generate optimized approximations of the inverse langevin function. Mathematics and
Mechanics of Solids, 24(7):2047–2059, 2019.

[85] R. Dargazany and M. Itskov. A network evolution model for the anisotropic Mullins effect in
carbon black filled rubbers. International Journal of Solids and Structures, 46:2967, 2009.

[86] R. Dargazany and M. Itskov. Constitutive modeling of the mullins effect and cyclic stress

softening in filled elastomers. Physical Review E, 88(1):012602, 2013.

[87] Vahid Morovati and Roozbeh Dargazany. Net v1. 0: A framework to simulate permanent
damage in elastomers under quasi-static deformations. SoftwareX, 10:100229, 2019.

[88] J.S. Bergström and M.C. Boyce. Mechanical behavior of particle filled elastomers. Rubber

Chem. Technol., 72:633, 1999.

[89] R. Dargazany, V. N. Khiêm, U. Navrath, and M. Itskov. Network evolution model of
anisotropic stress softening in filled rubber-like materials; parameter identification and fi-
nite element implementation. Journal of Mechanics of Materials and Structures, 7(8):861,
2013.

147

[90] P-Y Le Gac, P Davies, and D Choqueuse. Evaluation of long term behaviour of polymers for
offshore oil and gas applications. Oil & Gas Science and Technology Revue d’IFP Energies
nouvelles, 70(2):279–289, 2015.

[91] A Stevenson. On the durability of rubber/metal bonds in seawater. International Journal of

Adhesion and adhesives, 5(2):81–91, 1985.

[92] Guilherme Machado, Grégory Chagnon, and Denis Favier. Analysis of the isotropic models
of the mullins effect based on filled silicone rubber experimental results. Mechanics of
Materials, 42(9):841–851, 2010.

[93] Yunlong Li, Quan Wang, and Shijie Wang. A review on enhancement of mechanical and
tribological properties of polymer composites reinforced by carbon nanotubes and graphene
sheet: molecular dynamics simulations. Composites Part B: Engineering, 160:348–361,
2019.

[94] Alborz Izadi, Branton Toback, Mayank Sinha, Andrew Millar, Michael Bigelow, Sara
Roccabianca, and Rebecca Anthony. Mechanical and optical properties of stretch-
able silicon nanocrystal/polydimethylsiloxane nanocomposites. physica status solidi (a),
217(17):2000015, 2020.

[95] Alborz Izadi, Mayank Sinha, Cameron Papson, Sara Roccabianca, and Rebecca Anthony.
Mechanical behavior of sinc layers on pdms: effects of layer thickness, pdms modulus, and
sinc surface functionality. RSC Adv., 10:39087–39091, 2020.

[96] Baohui Shi, Man Zhang, Shengkai Liu, Baozhong Sun, and Bohong Gu. Multi-scale ageing
mechanisms of 3d four directional and five directional braided composites’ impact frac-
ture behaviors under thermo-oxidative environment. International Journal of Mechanical
Sciences, 155:50–65, 2019.

[97] Man Wang, Xiaowei Xu, Jin Ji, Yang Yang, Jianfeng Shen, and Mingxin Ye. The hy-
grothermal aging process and mechanism of the novolac epoxy resin. Composites Part B:
Engineering, 107:1–8, 2016.

[98] Lin Sang, Chuo Wang, Yingying Wang, and Wenbin Hou. Effects of hydrothermal aging
on moisture absorption and property prediction of short carbon fiber reinforced polyamide
6 composites. Composites Part B: Engineering, 153:306–314, 2018.

[99] Robert Bernstein, Dora K Derzon, and Kenneth T Gillen. Nylon 6.6 accelerated aging stud-
ies: thermal–oxidative degradation and its interaction with hydrolysis. Polymer Degradation
and Stability, 88(3):480–488, 2005.

[100] Bruno Musil, Michael Johlitz, and Alexander Lion. On the ageing behaviour of nbr: chemo-
mechanical experiments, modelling and simulation of tension set. Continuum Mechanics
and Thermodynamics, pages 1–17, 2018.

[101] Adelina David, Jing Huang, Emmanuel Richaud, and Pierre Yves LE Gac. Impact of ther-
mal oxidation on mechanical behaviour of polydicylopentadiene: Case of non-diffusion
limited oxidation. Polymer Degradation and Stability, page 109294, 2020.

148

[102] Jing Huang, Pierre Yves LE Gac, and Emmanuel Richaud. Thermal oxidation of poly
(dicyclopentadiene)-effect of phenolic and hindered amine stabilizers. Polymer Degradation
and Stability, page 109267, 2020.

[103] Rami Bouaziz, Komla Dela Ahose, Stéphane Lejeunes, Dominique Eyheramendy, and
Franck Sosson. Characterization and modeling of filled rubber submitted to thermal ag-
ing. International Journal of Solids and Structures, 169:122–140, 2019.

[104] Ricardo Baumhardt-Neto and Marco-A De Paoli. Mechanical degradation of polypropylene:

effect of uv irradiation. Polymer degradation and stability, 40(1):59–64, 1993.

[105] J-F Larché, P-O Bussière, and J-L Gardette. Photo-oxidation of acrylic-urethane thermoset
networks. relating materials properties to changes of chemical structure. Polymer degrada-
tion and stability, 96(8):1438–1444, 2011.

[106] Pierre-Olivier Bussiere, Agnes Rivaton, Sandrine Thérias, and Jean-Luc Gardette. Multi-
scale investigation of the poly (n-vinylcarbazole) photoageing mechanism. The Journal of
Physical Chemistry B, 116(2):802–812, 2012.

[107] Mathew C Celina. Review of polymer oxidation and its relationship with materials per-
formance and lifetime prediction. Polymer Degradation and Stability, 98(12):2419–2429,
2013.

[108] Mat Celina, Kenneth T Gillen, and RA Assink. Accelerated aging and lifetime prediction:
review of non-arrhenius behaviour due to two competing processes. Polymer Degradation
and stability, 90(3):395–404, 2005.

[109] L Steinke, U Veltin, M Flamm, A Lion, and M Celina. Numerical analysis of the heteroge-
neous ageing of rubber products. Constitutive Models for Rubber, 7:155–160, 2011.

[110] Kenneth T Gillen, Robert Bernstein, and Mathew Celina. Non-arrhenius behavior for ox-
idative degradation of chlorosulfonated polyethylene materials. Polymer degradation and
stability, 87(2):335–346, 2005.

[111] M Celina, J Wise, DK Ottesen, KT Gillen, and RL Clough. Oxidation profiles of thermally

aged nitrile rubber. Polymer Degradation and Stability, 60(2-3):493–504, 1998.

[112] Kenneth T Gillen, Mathew Celina, Roger L Clough, and Jonathan Wise. Extrapolation of
accelerated aging data-arrhenius or erroneous? Trends in polymer science, 8(5):250–257,
1997.

[113] John A Shaw, Alan S Jones, and Alan S Wineman. Chemorheological response of elas-
tomers at elevated temperatures: experiments and simulations. Journal of the Mechanics
and Physics of Solids, 53(12):2758–2793, 2005.

[114] N Rabanizada, F Lupberger, M Johlitz, and A Lion. Experimental investigation of the dy-
namic mechanical behaviour of chemically aged elastomers. Archive of Applied Mechanics,
85(8):1011–1023, 2015.

149

[115] RJ Pazur, I Lopez-Carreon, CG Porter, A Herzig, and M Johlitz. Thermal ageing of
peroxide-cured nbr, part i: Effect on bulk properties. In Constitutive Models for Rubber
XI: Proceedings of the 11th European Conference on Constitutive Models for Rubber (EC-
CMR 2019), June 25-27, 2019, Nantes, France, pages 581–586. CRC Press, 2019.

[116] AC Vieira, JC Vieira, JM Ferra, FD Magalhães, RM Guedes, and AT Marques. Mechanical
study of pla–pcl fibers during in vitro degradation. Journal of the mechanical behavior of
biomedical materials, 4(3):451–460, 2011.

[117] Vahid Morovati and Roozbeh Dargazany. Micro-mechanical modeling of the stress soften-
ing in double-network hydrogels. International Journal of Solids and Structures, 2019.

[118] Aref Ghaderi, Vahid Morovati, and Roozbeh Dargazany. A physics-informed assembly of
feed-forward neural network engines to predict inelasticity in cross-linked polymers. arXiv
preprint arXiv:2007.03067, 2020.

[119] Aref Ghaderi, Vahid Morovati, and Roozbeh Dargazany. A bayesian surrogate con-
arXiv preprint

stitutive model to estimate failure probability of rubber-like materials.
arXiv:2010.13241, 2020.

[120] T Ha-Anh and Toan Vu-Khanh. Prediction of mechanical properties of polychloroprene

during thermo-oxidative aging. Polymer testing, 24(6):775–780, 2005.

[121] A Lion and M Johlitz. On the representation of chemical ageing of rubber in continuum
mechanics. International Journal of Solids and Structures, 49(10):1227–1240, 2012.

[122] Michael Johlitz. On the representation of ageing phenomena. The Journal of Adhesion,

88(7):620–648, 2012.

[123] Benedikt Dippel, Michael Johlitz, and Alexander Lion. Ageing of polymer bonds: a cou-
pled chemomechanical modelling approach. Continuum Mechanics and Thermodynamics,
26(3):247–257, 2014.

[124] Shabnam Konica and Trisha Sain. A thermodynamically consistent chemo-mechanically
coupled large deformation model for polymer oxidation. Journal of the Mechanics and
Physics of Solids, page 103858, 2019.

[125] TRISHA SAIN and SHABNAM KONICA. A coupled thermo-chemo-mechanical model
for high temperature oxidations in polymers and polymer composites. In Proceedings of the
American Society for Composites—Thirty-third Technical Conference, 2018.

[126] Yiyi Zhang, Yi Li, Hanbo Zheng, Mengzhao Zhu, Jiefeng Liu, Tao Yang, Chaohai Zhang,
and Yang Li. Microscopic reaction mechanism of the production of methanol during the
thermal aging of cellulosic insulating paper. Cellulose, 27(5):2455–2467, 2020.

[127] Khashayar Mohammadi, Ashkan Ali Madadi, Zahra Bajalan, and Hossein Nejat Pishkenari.
Analysis of mechanical and thermal properties of carbon and silicon nanomaterials using a
coarse-grained molecular dynamics method. International Journal of Mechanical Sciences,
page 106112, 2020.

150

[128] C Schlomka, J Ihlemann, and C Naumann. Simulation of oxidative aging processes in
elastomer components using a dynamic network model. In Constitutive Models for Rubber
X, pages 77–82. CRC Press, 2017.

[129] C Naumann and J Ihlemann. A dynamic network model to simulate chemical aging pro-

cesses in elastomers. Constitutive Models for Rubber IX, 9:39, 2015.

[130] Darcy Beurle, Markus André, Udo Nackenhorst, and Rodrigue Desmorat. Micromechani-
cally motivated model for oxidation ageing of elastomers. In Virtual Design and Validation,
pages 271–290. Springer, 2020.

[131] Hamid Mohammadi and Roozbeh Dargazany. Micro-mechanical model for thermo-
In ASME 2018 International Mechanical Engineering
oxidative aging of elastomers.
Congress and Exposition, pages V009T12A028–V009T12A028. American Society of Me-
chanical Engineers, 2018.

[132] Amir Bahrololumi and Roozbeh Dargazany. A micro-mechanical approach to model hy-
drolytic aging in adhesives. In American Society of Mechanical Engineers (ASME) 2019
International Mechanical Engineering Congress and Exposition, 2019.

[133] Hamid Mohammadi and Roozbeh Dargazany. A micro-mechanical approach to model ther-
mal induced aging in elastomers. International Journal of Plasticity, 118:1–16, 2019.

[134] Mariacristina Gagliardi, Pietro Lenarda, and Marco Paggi. A reaction-diffusion formulation
to simulate eva polymer degradation in environmental and accelerated ageing conditions.
Solar Energy Materials and Solar Cells, 164:93–106, 2017.

[135] Pellegrino Musto, Giuseppe Ragosta, Mario Abbate, and Gennaro Scarinzi. Photo-oxidation
of high performance epoxy networks: correlation between the molecular mechanisms of
degradation and the viscoelastic and mechanical response. Macromolecules, 41(15):5729–
5743, 2008.

[136] Yoichi Kodera and Benjamin J McCoy. Distribution kinetics of radical mechanisms: re-

versible polymer decomposition. AIChE journal, 43(12):3205–3214, 1997.

[137] Gottfried W Ehrenstein and Sonja Pongratz. Beständigkeit von Kunststoffen. Hanser

München, 2007.

[138] Jacques Verdu. Oxydative ageing of polymers. John Wiley & Sons, 2012.

[139] Hesam Khajehsaeid. Mullins thresholds in context of the network alteration theories. Inter-

national Journal of Mechanical Sciences, 123:43–53, 2017.

[140] Yoshihiro Tomita, Keisuke Azuma, and Masato Naito. Computational evaluation of strain-
rate-dependent deformation behavior of rubber and carbon-black-filled rubber under mono-
tonic and cyclic straining. International Journal of Mechanical Sciences, 50(5):856–868,
2008.

151

[141] Yoshihiro Tomita, Wei Lu, Masato Naito, and Yasuhiro Furutani. Numerical evaluation
of micro-to macroscopic mechanical behavior of carbon-black-filled rubber. International
Journal of Mechanical Sciences, 48(2):108–116, 2006.

[142] Vahid Morovati and Roozbeh Dargazany. Improved non-gaussian statistical theory of poly-

mer physicsfor short chains. Physical Review E, 2019.

[143] G. Marckmann, E. Verron, L. Gornet, G. Chagnon, P. Charrier, and P. Fort. A theory of
network alteration for the Mullins effect. Journal of the Mechanics and Physics of Solids.,
50:2011, 2002.

[144] AR Khoei, A Rezaei Sameti, and Y Nikravesh Kazerooni. A continuum-atomistic multi-
scale technique for nonlinear behavior of nano-materials. International Journal of Mechan-
ical Sciences, 148:191–208, 2018.

[145] Jason Mulderrig, Bin Li, and Nikolaos Bouklas. Affine and non-affine microsphere models
for chain scission in polydisperse elastomer networks. Mechanics of Materials, 160:103857,
2021.

[146] Shabnam Konica and Trisha Sain. A reaction-driven evolving network theory coupled
with phase-field fracture to model polymer oxidative aging. Journal of the Mechanics and
Physics of Solids, 150:104347, 2021.

[147] Maha Zaghdoudi, Anja Kömmling, Matthias Jaunich, and Dietmar Wolff. Erroneous or
arrhenius: A degradation rate-based model for epdm during homogeneous ageing. Polymers,
12(9):2152, 2020.

[148] Shabnam Konica and Trisha Sain. A homogenized large deformation constitutive model for
high temperature oxidation in fiber-reinforced polymer composites. Mechanics of Materials,
160:103994, 2021.

[149] Amina Dinari, Fahmi Zaïri, Makram Chaabane, Jewan Ismail, and Tarek Benameur.
Thermo-oxidative stress relaxation in carbon-filled sbr. Plastics, Rubber and Composites,
pages 1–16, 2021.

[150] Makram Chaabane, Ning Ding, and Fahmi Zaïri. An approach to assess the thermal aging
effects on the coupling between inelasticity and network alteration in filled rubbers. Inter-
national Journal of Non-Linear Mechanics, page 103783, 2021.

[151] Zisheng Liao, Mokarram Hossain, Xiaohu Yao, Markus Mehnert, and Paul Steinmann. On
thermo-viscoelastic experimental characterization and numerical modelling of vhb polymer.
International Journal of Non-Linear Mechanics, 118:103263, 2020.

[152] Qiang Guo and Fahmi Zaïri. A physically-based thermo-mechanical model for stretch-
induced crystallizable rubbers: Crystallization thermodynamics and chain-network crystal-
lization anisotropy. International Journal of Plasticity, 131:102724, 2020.

[153] Qiang Guo, Fahmi Zaïri, and Xinglin Guo. Thermodynamics and mechanics of stretch-

induced crystallization in rubbers. Physical Review E, 97(5):052501, 2018.

152

[154] Zeang Zhao, Dong Wu, Ming Lei, Qiang Zhang, Panding Wang, and Hongshuai Lei. Me-
chanical behaviors and the equivalent network model of self-similar multinetwork elas-
tomers. International Journal of Solids and Structures, page 111135, 2021.

[155] Mahrez Saadedine, Fahmi Zaïri, Nourdine Ouali, Abderrahman Tamoud, and Amar Mes-
bah. A micromechanics-based model for visco-super-elastic hydrogel-based nanocompos-
ites. International Journal of Plasticity, page 103042, 2021.

[156] Mariusz Pyrz and Fahmi Zairi. Identification of viscoplastic parameters of phenomenolog-
ical constitutive equations for polymers by deterministic and evolutionary approach. Mod-
elling and Simulation in Materials Science and Engineering, 15(2):85, 2007.

[157] Qiang Guo and Fahmi Zaïri. A micromechanics-based model for deformation-induced dam-
age and failure in elastomeric media. International Journal of Plasticity, page 102976, 2021.

[158] M Ben Hassine, Moussa Naït-Abdelaziz, Fahmi Zaïri, Xavier Colin, C Tourcher, and Gre-
gory Marque. Time to failure prediction in rubber components subjected to thermal ageing:
A combined approach based upon the intrinsic defect concept and the fracture mechanics.
Mechanics of Materials, 79:15–24, 2014.

[159] Yongqiang Du, Jian Zheng, and Guibo Yu. Influence of thermally-accelerated aging on the
dynamic mechanical properties of htpb coating and crosslinking density-modified model for
the payne effect. Polymers, 12(2):403, 2020.

[160] Fang Xie, Liwu Liu, Xiaobo Gong, Longnan Huang, Jinsong Leng, and Yanju Liu. Effects
of accelerated aging on thermal, mechanical and shape memory properties of cyanate-based
shape memory polymer: I vacuum ultraviolet radiation. Polymer Degradation and Stability,
138:91–97, 2017.

[161] Alok Ranjan Verma and B Subba Reddy. Accelerated aging studies of silicon-rubber based

polymeric insulators used for hv transmission lines. Polymer Testing, 62:124–131, 2017.

[162] Ana P Cysne Barbosa, Ana Paula P Fulco, Erick SS Guerra, Franscisco K Arakaki, Márcio
Tosatto, Maria Carolina B Costa, and José Daniel D Melo. Accelerated aging effects on
carbon fiber/epoxy composites. Composites Part B: Engineering, 110:298–306, 2017.

[163] Maha Zaghdoudi, Anja Kömmling, Matthias Jaunich, and Dietmar Wolff. Erroneous or
arrhenius: A degradation rate-based model for epdm during homogeneous ageing. Polymers,
12(9):2152, 2020.

[164] A Simar, M Gigliotti, JC Grandidier, and I Ammar-Khodja. Evidence of thermo-oxidation
phenomena occurring during hygrothermal aging of thermosetting resins for rtm composite
applications. Composites Part A: Applied Science and Manufacturing, 66:175–182, 2014.

[165] Adam Quintana and Mathew C Celina. Overview of dlo modeling and approaches to pre-
dict heterogeneous oxidative polymer degradation. Polymer Degradation and Stability,
149:173–191, 2018.

153

[166] Erik Linde, Fritjof Nilsson, Matija Barrett, Mikael S Hedenqvist, and Mathew C Celina.
Time-and feedback-dependent dlo phenomena in oxidative polymer aging. Polymer Degra-
dation and Stability, 189:109598, 2021.

[167] Pei Ni and Ronald W Thring. Synthesis of polyurethanes from solvolysis lignin using a
polymerization catalyst: mechanical and thermal properties. International Journal of Poly-
meric Materials, 52(8):685–707, 2003.

[168] P.J. Flory and J.R. Rehner. Statistical mechanics of cross-linked polymer networks i. rub-

berlike elasticity. Journal of Chemical Physics, 11:512, 1943.

[169] Morgan Deroiné, Antoine Le Duigou, Yves-Marie Corre, Pierre-Yves Le Gac, Peter Davies,
Guy César, and Stéphane Bruzaud. Accelerated ageing of polylactide in aqueous environ-
ments: Comparative study between distilled water and seawater. Polymer degradation and
stability, 108:319–329, 2014.

[170] Morgan Deroiné, Antoine Le Duigou, Yves-Marie Corre, Pierre-Yves Le Gac, Peter Davies,
Guy César, and Stéphane Bruzaud. Seawater accelerated ageing of poly (3-hydroxybutyrate-
co-3-hydroxyvalerate). Polymer degradation and stability, 105:237–247, 2014.

[171] Hanen Mahjoubi, Fahmi Zaïri, and Zoubeir Tourki. A micro-macro constitutive model for
strain-induced molecular ordering in biopolymers: application to polylactide over a wide
range of temperatures. International Journal of Plasticity, 123:38–55, 2019.

[172] Reza Rastak and Christian Linder.

A non-affine micro-macro approach to strain-
crystallizing rubber-like materials. Journal of the Mechanics and Physics of Solids, 111:67–
99, 2018.

[173] Benjamin C Marchi and Ellen M Arruda. Generalized error-minimizing, rational inverse
langevin approximations. Mathematics and Mechanics of Solids, 0(0):1081286517754131,
2018.

[174] Vu Ngoc Khiêm and Mikhail Itskov. Analytical network-averaging of the tube model::
Rubber elasticity. Journal of the Mechanics and Physics of Solids, 95:254–269, 2016.

[175] Vahid Morovati, Amir Bahrololoumi, and Roozbeh Dargazany. Fatigue-induced stress-
In-
softening in cross-linked multi-network elastomers: Effect of damage accumulation.
ternational Journal of Plasticity, 142:102993, 2021.

[176] Michael R Buche and Meredith N Silberstein. Chain breaking in the statistical mechanical
constitutive theory of polymer networks. Journal of the Mechanics and Physics of Solids,
156:104593, 2021.

[177] Aref Ghaderi, Vahid Morovati, and Roozbeh Dargazany. A bayesian surrogate constitutive
model to estimate failure probability of elastomers. Mechanics of Materials, 162:104044,
2021.

154

[178] Ieuan Collins, Mokarram Hossain, Wulf Dettmer, and Ian Masters. Flexible membrane
structures for wave energy harvesting: A review of the developments, materials and compu-
tational modelling approaches. Renewable and Sustainable Energy Reviews, 151:111478,
2021.

[179] Chennakesava Kadapa, Zhanfeng Li, Mokarram Hossain, and Jiong Wang. On the advan-
tages of mixed formulation and higher-order elements for computational morphoelasticity.
Journal of the Mechanics and Physics of Solids, 148:104289, 2021.

[180] Bruno Musil, Michael Johlitz, and Alexander Lion. On the ageing behaviour of nbr: chemo-
mechanical experiments, modelling and simulation of tension set. Continuum Mechanics
and Thermodynamics, 32(2):369–385, 2020.

[181] Moussa Nait Abdelaziz, Georges Ayoub, Xavier Colin, Mouna Benhassine, and
M Mouwakeh. New developments in fracture of rubbers: Predictive tools and influence
of thermal aging. International Journal of Solids and Structures, 165:127–136, 2019.

[182] Mariaenrica Frigione and Alvaro Rodríguez-Prieto. Can accelerated aging procedures pre-
dict the long term behavior of polymers exposed to different environments? Polymers,
13(16):2688, 2021.

[183] Makram Chaabane, Ning Ding, and Fahmi Zaïri. An approach to assess the thermal aging
effects on the coupling between inelasticity and network alteration in filled rubbers. Inter-
national Journal of Non-Linear Mechanics, 136:103783, 2021.

[184] Narueporn Payungwong, Surakit Tuampoemsab, Porntip Rojruthai, and Jitladda Sakdapi-
panich. The role of model fatty acid and protein on thermal aging and ozone resistance of
peroxide vulcanized natural rubber. Journal of Rubber Research, pages 1–11, 2021.

[185] Réda Kadri, M Nait Abdelaziz, Bruno Fayolle, M Ben Hassine, and Jean-François Witz. A
unified mechanical based approach to fracture properties estimates of rubbers subjected to
aging. International Journal of Solids and Structures, 234:111305, 2022.

[186] Rami Bouaziz, Laurianne Truffault, Rouslan Borisov, Cristian Ovalle, Lucien Laiarinan-
drasana, Guillaume Miquelard-Garnier, and Bruno Fayolle. Elastic properties of polychloro-
prene rubbers in tension and compression during ageing. Polymers, 12(10):2354, 2020.

[187] Maryam Shakiba and Aimane Najmeddine.

Physics-based constitutive equation for
thermo-chemical aging in elastomers based on crosslink density evolution. arXiv preprint
arXiv:2104.09001, 2021.

[188] Camille Loyer, Gilles Regnier, V Duval, and Emmanuel Richaud. Multiscale study of poly
(butylene terephthalate) hydrolysis. Polymer Degradation and Stability, 192:109690, 2021.

[189] Aref Ghaderi, Vahid Morovati, Yang Chen, and Roozbeh Dargazany. A physics-informed
Interna-

multi-agents model to predict thermo-oxidative/hydrolytic aging of elastomers.
tional Journal of Mechanical Sciences, page 107236, 2022.

155

[190] Aimane Najmeddine, Zhen Xu, Guoliang Liu, Alan R Esker, and Maryam Shakiba. Physics-
based constitutive modeling of photo-oxidative aging in semi-crystalline polymers based on
chemical characterization techniques. arXiv preprint arXiv:2108.07143, 2021.

[191] Hamid Mohammadi, Vahid Morovati, Abd-Elrahman Korayem, Emad Poshtan, and
Roozbeh Dargazany. Constitutive modeling of elastomers during photo-and thermo-
oxidative aging. Polymer Degradation and Stability, 191:109663, 2021.

[192] Naima Rezig, Tassadit Bellahcene, Méziane Aberkane, and Moussa Nait Abdelaziz.
Thermo-oxidative ageing of a sbr rubber: effects on mechanical and chemical properties.
Journal of Polymer Research, 27(11):1–13, 2020.

[193] Benedikt Dippel, Michael Johlitz, and Alexander Lion. Ageing of polymer bonds: a cou-
pled chemomechanical modelling approach. Continuum Mechanics and Thermodynamics,
26(3):247–257, 2014.

[194] Hamza Lamnii, M Nait Abdelaziz, Georges Ayoub, Xavier Colin, and U Maschke. Experi-
mental investigation and modeling attempt on the effects of ultraviolet aging on the fatigue
behavior of an ldpe semi-crystalline polymer. International Journal of Fatigue, 142:105952,
2021.

[195] Wenbo Luo, Youjian Huang, Boyuan Yin, Xia Jiang, and Xiaoling Hu. Fatigue life assess-
ment of filled rubber by hysteresis induced self-heating temperature. Polymers, 12(4):846,
2020.

[196] Mahdi Ghazavi, Seyed Farhad Abdollahi, and M Emin Kutay. Implementation of nchrp 9-
44a fatigue endurance limit prediction model in mechanistic-empirical asphalt pavement
analysis web application. Transportation Research Record, page 03611981221076114,
2022.

[197] Zhongtong Wang, Jingda Tang, Ruobing Bai, Wenlei Zhang, Tongda Lian, Tongqing Lu,
and Tiejun Wang. A phenomenological model for shakedown of tough hydrogels under
cyclic loads. Journal of Applied Mechanics, 85(9), 2018.

[198] A Mousa, US Ishiaku, and ZA Mohd Ishak. Thermo-oxidative aging and fatigue behav-
ior of dynamically vulcanized pvc/enr thermoplastic elastomers. International Journal of
Polymeric Materials, 51(11):967–980, 2002.

[199] Weifeng Chen, Fei Xue, Dunji Yu, Weiwei Yu, Xu Chen, et al. Effect of thermal aging on
the low cycle fatigue behavior of z3cn20. 09m cast duplex stainless steel. Materials Science
and Engineering: A, 646:263–271, 2015.

[200] Jian Zhang, Valerian Hirschberg, and Denis Rodrigue. Mechanical fatigue of biodegradable
polymers: A study on polylactic acid (pla), polybutylene succinate (pbs) and polybutylene
adipate terephthalate (pbat). International Journal of Fatigue, 159:106798, 2022.

[201] PYDI Hariprasadarao, H Ravisankar, and V CHITTARANJAN Das. Effect of thermal aging
on mechanical properties and fatigue behavior of pp/hdpe thermoplastic polymer nanoclay
composites. Int. J. Mech. Prod. Eng. Res. Dev, 8:349–358, 2018.

156

[202] Amélie Malpot, Fabienne Touchard, and Sébastien Bergamo. Influence of moisture on the
fatigue behaviour of a woven thermoplastic composite used for automotive application. Ma-
terials & Design, 98:12–19, 2016.

[203] Hsoung-Wei Chou, Jong-Shin Huang, and Shih-Ting Lin. Effects of thermal aging on
Journal of applied polymer science,

fatigue of carbon black–reinforced epdm rubber.
103(2):1244–1251, 2007.

[204] Berkan Tan and Erol Feyzullaho˘glu. The investigation of effects of environmental condi-
tions on fatigue strength and spring stiffness of rubber materials. Journal of Elastomers &
Plastics, 53(6):632–652, 2021.

[205] Abbas Tcharkhtchi, Sedigheh Farzaneh, Sofiane Abdallah-Elhirtsi, Bardia Esmaeillou, Fa-
bien Nony, and A Baron. Thermal aging effect on mechanical properties of polyurethane.
International Journal of Polymer Analysis and Characterization, 19(7):571–584, 2014.

[206] Byungwoo Moon, Keekeun Kim, Keontae Park, Soo Park, and Chang-Sung Seok. Fatigue
life prediction of tire sidewall using modified arrhenius equation. Mechanics of Materials,
147:103405, 2020.

[207] Christoph Neuhaus, Alexander Lion, Michael Johlitz, Paul Heuler, Matthias Barkhoff, and
Frank Duisen. Fatigue behaviour of an elastomer under consideration of ageing effects.
International Journal of Fatigue, 104:72–80, 2017.

[208] Alexandre Jouan and Andrei Constantinescu. A modified dissipated energy fatigue criterion
to consider the thermo-oxidative ageing of electrically conductive silicone adhesive joints.
International Journal of Fatigue, 116:68–79, 2018.

[209] Xiaolei Wang, Kun Yang, and Ping Zhang. Evaluation of the aging coefficient and the ag-
ing lifetime of carbon black-filled styrene-isoprene-butadiene rubber after thermal-oxidative
aging. Composites Science and Technology, page 109258, 2022.

[210] Jieying Zhi, Qinglin Wang, Mengjie Zhang, Zhenze Zhou, Anna Liu, and Yuxi Jia. Coupled
analysis on hyper-viscoelastic mechanical behavior and macromolecular network alteration
of rubber during thermo-oxidative aging process. Polymer, 171:15–24, 2019.

[211] Maha Zaghdoudi, Anja Kömmling, Matthias Jaunich, and Dietmar Wolff. Scission, cross-
linking, and physical relaxation during thermal degradation of elastomers. Polymers,
11(8):1280, 2019.

[212] Tesfaldet Mengistu and Richard J Pazur.

The thermal oxidation of hydrogenated
acrylonitrile-co-butadiene rubber from ambient to 150° c. Polymer Degradation and Stabil-
ity, 188:109574, 2021.

[213] Zeng Zhang, Hui Li, Haiyue Zhou, Liying Zhao, Shikun Li, Hongzhen Wang, Zenglin
Wang, and Zaifeng Li. Modification and improvement of aging resistance for hnbr/graphite
composites. Journal of Macromolecular Science, Part A, pages 1–10, 2022.

157

[214] Jannik Voges and Daniel Juhre. A phase field model for thermo-oxidative aging in cracked

polymers. PAMM, 19(1):e201900384, 2019.

[215] Mathew Celina, Erik Linde, Douglas Brunson, Adam Quintana, and Nicholas Giron.
Overview of accelerated aging and polymer degradation kinetics for combined radiation-
thermal environments. Polymer Degradation and Stability, 166:353–378, 2019.

[216] G. Ayoub, F. Zaïri, M. Naït-Abdelaziz, and J.M. Gloaguen. Modeling the low-cycle fatigue
behavior of visco-hyperelastic elastomeric materials using a new network alteration theory:
Application to styrene-butadiene rubber. Journal of the Mechanics and Physics of Solids,
59(2):473–495, 2011.

[217] Weitao Lou, Chaoyang Xie, and Xuefei Guan. Thermal-aging constitutive model for a
silicone rubber foam under compression. Polymer Degradation and Stability, page 109873,
2022.

[218] AF Blanchard and D Parkinson. Breakage of carbon-rubber networks by applied stress.

Rubber Chemistry and Technology, 25(4):808–842, 1952.

[219] Rubber Houwink. Slipping of molecules during the deformation of reinforced rubber. Rub-

ber Chemistry and Technology, 29(3):888–893, 1956.

[220] David E Hanson, Marilyn Hawley, Robert Houlton, Kiran Chitanvis, Philip Rae, E Bruce
Orler, and Debra A Wrobleski. Stress softening experiments in silica-filled polydimethyl-
siloxane provide insight into a mechanism for the mullins effect. Polymer, 46(24):10989–
10995, 2005.

[221] G Kraus, CW Childers, and KW Rollmann. Stress softening in carbon black-reinforced
Journal of Applied Polymer Science,

vulcanizates. strain rate and temperature effects.
10(2):229–244, 1966.

[222] M.F. Beatty. An average-stretch full-network model for rubber elasticity. Journal of Elas-

ticity, 70:65, 2003.

[223] Georges Ayoub, Fahmi Zaïri, Moussa Naït-Abdelaziz, Jean-Michel Gloaguen, and G Kridli.
A visco-hyperelastic damage model for cyclic stress-softening, hysteresis and permanent set
in rubber using the network alteration theory. International Journal of Plasticity, 54:19–33,
2014.

[224] G. Chagnon, E. Verron, and L. Marckmann, G.and Gornet. Developement of new constitu-
tive equations for the mullins effect in rubber using the network alteration theory. Interna-
tional Journal of Solids and Structures, 43:6817, 2006.

[225] VN Khiêm, J-B Le Cam, Sylvain Charlès, and M Itskov. Thermodynamics of strain-induced
crystallization in filled natural rubber under uni-and biaxial loadings, part ii: Physically-
based constitutive theory. Journal of the Mechanics and Physics of Solids, 159:104712,
2022.

158

[226] Ramin Akbari, Vahid Morovati, and Roozbeh Dargazany. Reverse physically motivated
frameworks for investigation of strain energy function in rubber-like elasticity. International
Journal of Mechanical Sciences, 221:107110, 2022.

[227] M. Klüppel. The role of disorder in filler reinforcement of elastomers on various length

scales. Advances in Polymer Science, 164:1, 2003.

[228] Jaemin Kim, Erik Mailand, Ida Ang, Mahmut Selman Sakar, and Nikolaos Bouklas. A
model for 3d deformation and reconstruction of contractile microtissues. Soft Matter,
17(45):10198–10209, 2021.

[229] Markus Mehnert, Mokarram Hossain, and Paul Steinmann. A complete thermo-electro-
viscoelastic characterization of dielectric elastomers, part ii: Continuum modeling ap-
proach. Journal of the Mechanics and Physics of Solids, 157:104625, 2021.

[230] Chennakesava Kadapa and Mokarram Hossain. A unified numerical approach for soft to
hard magneto-viscoelastically coupled polymers. Mechanics of Materials, page 104207,
2022.

[231] Y Nikravesh, A Rezaei Sameti, and AR Khoei. An atomistic–continuum multiscale analysis
for heterogeneous nanomaterials and its application in nanoporous gold foams. Applied
Mathematical Modelling, 2022.

159